## Automatically adapted for numpy.oldnumeric Jul 23, 2007 by 

import numpy.oldnumeric as Numeric
import numpy

class Grid3D:
    """Class to represent 3D discret volumetric data.
The data has to be a Numeric array od shape 3
"""
    
    def __init__(self, data, origin, stepSize, header, crystal=None):

        assert isinstance(data, Numeric.ArrayType)
        assert len(data.shape)==3
        assert len(origin)==3
        assert len(stepSize)==3

        # header dict
        self.header = header
        
        # actual values at the grid points
        self.data = data

        # number of grid points along 3 dimensions
        self.dimensions = data.shape

        # 3D coordinates of lower left front corner (i.e. (0,0,0) grid point) 
        self.origin = origin

        # grid step size
        self.stepSize = stepSize

        # crystal object, used to cnvert between fractioanl and cartesian space
        # will be an instance mglutil.math.crystal.Crystal
        self.crystal = crystal


    def gridPoints(self, cutLow, cutHigh):
        """
        get 3D coordinates of selected the grid points and corresponding values

        coords, values, indi, indj, indk <- gridPoints(cutLow, cutHigh)
        
        cutLow: low cutoff value for selecting grid points
        cutHigh: high cutoff value for selecting points
    
        returns:
            coords: list of 3D coordinates of the grid with values in
                    specified range
            values: values at these grid points
            indi:   list of i indices
            indj:   list of j indices
            indk:   list of k indices
        """
        # array of boolean for values larger
        gr = numpy.greater(self.data, cutLow)
        # array of boolean for values smaller
        ls = numpy.less(self.data, cutHigh)
        # array of boolean for large than cutLow and smaller than cutHigh
        inrange = gr*ls
        # list of indices
        indi, indj, indk = numpy.nonzero( inrange )

        coords = []
        values = []
        data = self.data
        ox, oy, oz = self.getOriginReal()
        step = self.getStepSizeReal()
        for i, j, k in zip( indi, indj, indk):
            values.append( data[i][j][k] )
            coords.append( (ox+i*step[0], oy+j*step[1], oz+k*step[2] ) )
        
        return coords, values, indi, indj, indk


    def centerPoint(self):
        """compute and return a center point of this grid"""
        orig = self.origin
        dims = self.dimensions
        size = self.stepSize
        pt = (orig[0] + (dims[0]*size[0])*0.5,
              orig[1] + (dims[1]*size[1])*0.5,
              orig[2] + (dims[2]*size[2])*0.5)
        if self.crystal:
            pt = grid3D.crystal.toCartesian(pt)
        return pt
        
    def stats(self):
        """returns the mimn, max, mean and standard deviation of a list of values"""
        return self.data.min(), self.data.max(), self.data.mean(), self.data.std()
    

    def transpose(self, copy=False):
        """
        transpose the grid i.e. switch between fortran and C
        """
        art = numpy.ascontiguousarray(numpy.transpose(self.data),
                                      self.data.dtype.char)
        if copy is False:
            self.data = art
            self.dimensions = art.shape
            self.stepSize = self.stepSize[::-1]
            newgrid = self
        else:
            newgrid = self.__class__( art, self.origin, self.stepSize[::-1],
                                      self.header)
        return newgrid


    def normalize(self):
        """if data is given the values in data are modified, else the Grid3D
data attribute is used as data
"""
        data = self.data

        mymin, mymax, mymean, mystdev = self.stats()
        for i in xrange(data.shape[0]):
            for j in xrange(data.shape[1]):
                for k in xrange(data.shape[2]):
                    data[i][j][k] = (data[i][j][k]-mymean)/mystdev
        self.header['amean'] = mymean
        self.header['arms'] = mystdev
        self.header['amax'] = mymax
        self.header['amin'] = mymin


    def getOriginReal(self):
        """return the origin of the grid in real space"""
        origin = self.origin
        if self.crystal:
            origin = self.crystal.toCartesian(self.origin)
        return origin


    def getStepSizeReal(self):
        """return the stepsize of the grid in real space"""
        stepSize = self.stepSize
        from math import sqrt
        if self.crystal:
            vx = self.crystal.toCartesian( (stepSize[0], 0, 0) )
            vy = self.crystal.toCartesian( (0, stepSize[1], 0) )
            vz = self.crystal.toCartesian( (0, 0, stepSize[2]) )
            stepSize = (sqrt(vx[0]*vx[0]+vx[1]*vx[1]+vx[2]*vx[2]),
                        sqrt(vy[0]*vy[0]+vy[1]*vy[1]+vy[2]*vy[2]),
                        sqrt(vz[0]*vz[0]+vz[1]*vz[1]+vz[2]*vz[2]))
        return stepSize


    def overlapWithBox(self, grid3D, Box):
        """Compute the indices of self.data representing the corners
of a box in the grid.
"""
        inds1 = [0,0,0]
        inds2 = [0,0,0]

        stepSize = grid3D.getStepSizeReal()

        origin = grid3D.getOriginReal()

        bside = (Box.xside, Box.yside, Box.zside) 

        # compute the box indices in related to the grid's origin
        # find the center of the box, go to each edge, then set the edges
        # in relation to the grid3D origin
        for i in (0,1,2):
            inds1[i] = ( Box.center[i] - bside[i]*0.5 - origin[i] ) / stepSize[i]
            inds2[i] = ( Box.center[i] + bside[i]*0.5 - origin[i] ) / stepSize[i]
            
        return inds1,inds2
        
    def overlapWithGrid(self, grid):
        """Compute the indices into self.data representing the overlap between
this grid and the one provided as an argument.  Both grids have to have the
same stepSize.
"""
        assert self.stepSize[0]==grid.stepSize[0]
        assert self.stepSize[1]==grid.stepSize[1]
        assert self.stepSize[2]==grid.stepSize[2]

        bc1 = [0,0,0]
        bc2 = [0,0,0]

        stepSize = self.getStepSizeReal()
        origin = self.getOriginReal()
        origin1 = grid.getOriginReal()

        shape1 = grid.data.shape
        for i in (0,1,2):
            bc1[i] = ( origin1[i] - origin[i] ) / stepSize[i]
            bc2[i] = bc1[i] + shape1[i] - 1

        shape = self.data.shape
        for i in (0,1,2):
            bc1[i] = int(round(bc1[i]))
            if bc1[i] < 0: bc1[i]=0
            if bc1[i] > shape[i]: bc1[i]=shape[i]
            bc2[i] = int(round(bc2[i]))
            if bc2[i] < 0: bc2[i]=0
            if bc2[i] > shape[i]: bc2[i]=shape[i]
            
        return bc1, bc2


    def get2DOrthoSlice(self, axis, sliceNum):
        """Return a 2D array corresponding to an axis aligned slice orthogonal
to axis at index along this axis.  We also return the real vertices
corresponding the corner of the quad for this slice.
axis cand be 'x', 'y' or 'z'
sliceNum varies from 0 to max(extent(axis))
"""

        if axis=='x': axisInd = 0
        elif axis=='y': axisInd = 1
        elif axis=='z': axisInd = 2
        else:
            raise ValueError ("axis can only be 'x', 'y', or 'z', got %s"%\
                              str(axis))

        dims = self.dimensions
        dx, dy, dz = (dims[0]-1,dims[1]-1,dims[2]-1)
        if sliceNum < 0:
            raise ValueError ("sliceNum has to be positive")
          
        if sliceNum > dims[axisInd]:
            raise ValueError ("sliceNum for axis %s too large, max %d got %d"%\
                              (axis, dims[axisInd], sliceNum))

        sx, sy, sz = self.stepSize
        ox, oy, oz = self.origin
        ex, ey, ez = (ox+sx*dx, oy+sy*dy, oz+sz*dz)
        if axis=='x':
            slice = Numeric.array(self.data[sliceNum,:,:])
            x = ox + (sx*0.5)+(sliceNum*sx)
            vertices = [ [x,oy,oz], [x,ey,oz], [x,ey,ez], [x,oy,ez] ]

        elif axis=='y':
            slice = Numeric.array(self.data[:,sliceNum,:])
            y = oy + (sy*0.5)+(sliceNum*sy)
            #vertices = [ [ox,y,oz], [ox,y,ez], [ex,y,ez], [ex,y,oz] ]
            vertices = [ [ox,y,oz], [ex,y,oz], [ex,y,ez], [ox,y,ez] ]

        elif axis=='z':
            slice = Numeric.array(self.data[:,:,sliceNum])
            z = oz + (sz*0.5)+(sliceNum*sz)
            vertices = [ [ox,oy,z], [ex,oy,z], [ex,ey,z], [ox,ey,z] ]

        if self.crystal:
            nv = []
            for v in vertices:
                nv.append( tuple(self.crystal.toCartesian(v)) )
            vertices = nv

        return slice, vertices


class Grid3DD(Grid3D):
    def __init__(self, data, origin, stepSize, header, crystal=None):
        assert data.dtype.char in [Numeric.Float64, Numeric.Float]
        Grid3D.__init__(self, data, origin, stepSize, header, crystal)
                        

class Grid3DF(Grid3D):
    def __init__(self, data, origin, stepSize, header, crystal=None):
        assert data.dtype.char in [Numeric.Float32, Numeric.Float16,
                                   Numeric.Float8, Numeric.Float0]
        Grid3D.__init__(self, data, origin, stepSize, header, crystal)
                        

class Grid3DI(Grid3D):
    def __init__(self, data, origin, stepSize, header, crystal=None):
        assert data.dtype.char in [Numeric.Int32, Numeric.Int]
        Grid3D.__init__(self, data, origin, stepSize, header, crystal)


class Grid3DSI(Grid3D):
    def __init__(self, data, origin, stepSize, header, crystal=None):
        assert data.dtype.char == Numeric.Int16
        Grid3D.__init__(self, data, origin, stepSize, header, crystal)


class Grid3DUI(Grid3D):
    def __init__(self, data, origin, stepSize, header, crystal=None):
        assert data.dtype.char in [Numeric.UInt32, Numeric.UInt]
        Grid3D.__init__(self, data, origin, stepSize, header, crystal)


class Grid3DUSI(Grid3D):
    def __init__(self, data, origin, stepSize, header, crystal=None):
        assert data.dtype.char == Numeric.UInt16
        Grid3D.__init__(self, data, origin, stepSize, header, crystal)


class Grid3DUC(Grid3D):
    def __init__(self, data, origin, stepSize, header, crystal=None):
        assert data.dtype.char in [Numeric.Int0, Numeric.UInt8]
        Grid3D.__init__(self, data, origin, stepSize, header, crystal)
        

ArrayTypeToGrid = {
    Numeric.Float64:Grid3DD, Numeric.Float:Grid3DD,
    Numeric.Float32:Grid3DF, Numeric.Float16:Grid3DF,
    Numeric.Float8:Grid3DF, Numeric.Float0:Grid3DF,
    Numeric.Int32:Grid3DI, Numeric.Int:Grid3DI,
    Numeric.Int16:Grid3DSI,
    Numeric.UInt8:Grid3DUC,
    }

GridTypeToArray = {
    'Grid3DD'  :Numeric.Float,
    'Grid3DF'  :Numeric.Float16,
    'Grid3DF'  :Numeric.Float0,
    'Grid3DI'  :Numeric.Int,
    'Grid3DSI' :Numeric.Int16,
    'Grid3DUC' :Numeric.UInt8,
}

