"""Regular grids together with values at the vertices."""

import numpy as np
cimport numpy as np
from scipy.ndimage import map_coordinates
from numpy.fft import fftshift, ifftshift, ifft2
from PIL import Image

cimport cython

from cython_gsl cimport *
gsl_rng_env_setup()
cdef:
    const gsl_rng_type * T
    gsl_rng * r
T = gsl_rng_default
r = gsl_rng_alloc (T)


from cryoem.data.grid import Grids
from cryoem.data.grid cimport Grid2D, Grid3D
from cryoem.data import mrc
from cryoem.data.cnv import Cnvarray
from cryoem.model.transformation cimport Transformation3D


cdef class GridDataReal:
    """Regular grid with real values at vertices."""

    def __cinit__(self, grid, values):
        """Initialise with Grid object and values on the grid.

        Args:
            grid: Grid object
            values: (N,) array, values at grid vertices.
        """

        self._grid = grid
        self._values = values

    @property
    def grid(self):
        return self._grid

    @property
    def values(self):
        return np.array(self._values)

    @values.setter
    def values(self, values_new):
        self._values = values_new

    @property
    def dimension(self):
        return self.grid.dimension

    @property
    def size(self):
        return len(self._values)

    @classmethod
    def create_from_function(cls, grid, function):
        """Sample function on grid, use values to initialise.

        Args:
            grid: Grid object
            function: callable object

        Returns:
            GridData object with function evaluated at grid points.
        """

        values = function(grid.points)

        return cls(grid, values.flatten())

    @classmethod
    def create_constant(cls, grid, value=0.0):
        """Initialise grid with constant value.

        Args:
            grid: Grid object.
            value: Constant value at all grid points, default to 0.0.

        Returns:
            GridData object with constant values.
        """

        N = np.prod(grid.shape)
        values = np.ones((N,), dtype='double') * value

        return cls(grid, values)


    def project_centered(self, transformation, shape, spacing):
        """Transform and project to lower dimensional grid.

        Args:
            transformation: Transformation object
            shape: The number of grid points in each direction
                as an integer.
            spacing: The grid spacing as a double.

        Returns:
            GridDataReal object of one lower dimension with projected data.
            The projected grid is centered.

        The current grid is assumed to be centered and should
        have dimension 2 or 3.

        ## Note that the transformation is applied to the reference frame,
        ## not to the object that is being projected.
        The transformation is applied to the object being projected,
        not the reference frame.
        """

        # Create rotated 2D or 3D grid on which to interpolate.
        # This could be different from the existing grid,
        # could be finer for instance.
        grid = Grids.create_centered(self.grid.dimension, shape, spacing)

        # Transform target grid points.
        # Take the inverse to ensure that the actual transformation
        # is applied to the object, not the reference frame.
        transformed = transformation.inverse_transform(grid.points)

        # First scale and shift them to be in the
        # funny coordinate system needed for 'map_coordinates'.
        transformed *= 1.0 / self.grid.spacing
        ## transformed += (self.grid.shape + 1.0) / 2.0
        transformed += (self.grid.shape - 1.0) / 2.0

        # Now interpolate from the existing values to this new grid.
        values = map_coordinates(self.reshaped(), transformed.T, order=3)

        # Sum over the z-coordinate.
        v = values.reshape(grid.shape)
        values_projected = v.sum(axis=-1) * grid.spacing[0]

        # Create final grid.
        grid_projected = Grids.create_centered(self.grid.dimension - 1,
                                               shape, spacing)

        return GridDataReal(grid_projected,
                            values_projected.reshape(-1, 1).flatten())



    def reshaped(self):
        """Returns values reshaped according to grid."""

        # Copy memoryview contents to new numpy array.
        numpy_array = np.zeros(len(self._values))
        numpy_array[:] = self._values[:]

        return numpy_array.reshape(self.grid.shape)

    def save(self, parent_group, name='griddata'):
        """Save to group in HDF5 database.

        Assume that the group parent_group[name] doesn't exist.
        ## Later: Delete group if it exists.

        Args:
            parent_group: Group within which to create the griddata group.
            name: Name of the new griddata group.
        """

        ## if name in parent_group:
        ##     print 'Removing group %s from parent group %s' % (
        ##         name, parent_group.name)
        ##     del parent_group[name]
        group = parent_group.create_group(name)
        self.grid.save(group, 'grid')
        dataset = group.create_dataset('values',
                                       (<object>self._values).shape,
                                       np.double)
        dataset[...] = self._values

    @classmethod
    def load(cls, parent_group, name='griddata'):
        """Loads GridData from group in HDF5 database."""

        group = parent_group[name]
        grid = Grids.load(group, 'grid')
        values = group['values'][:]

        return cls(grid, values)

    def save_mrc(self, mrc_filename):
        """Save griddata as mrc file."""

        mrc.write(self.reshaped(), mrc_filename,
                  spacing=self.grid.spacing,
                  origin=self.grid.origin)

    @classmethod
    def load_mrc(cls, mrc_filename):
        """Load griddata from mrc file."""

        data, spacing, origin = mrc.read(mrc_filename)
        grid = Grids.create(3, data.shape, np.array(spacing),
                            np.array(origin))

        return cls(grid, data.astype('double').flatten())


    def to_counts(self, total_counts):
        """Convert to count data by scaling and rounding."""

        total = np.sum(self.values)
        values_int = np.round((
            self.values *
            (np.ones(len(self.values)) * total_counts / total)))

        return GridDataInt(self.grid, values_int.astype('int'))

    def to_real(self):
        """To support common interface with GridDataInt."""

        return self

    def to_spatial_domain(self):
        """Convert values and grid to spatial domain.

        Assume that DC is in the center.

        Assume that in both domains the signal is real and symmetric.
        """

        Xs = fftshift(self.reshaped()) # DC to the corner.
        xs = ifft2(Xs) # Inverse FFT.
        x = ifftshift(xs) # Make it symmetric around center.
        x = np.real(x) # Assume that original data was symmetric.

        return GridDataReal(self.grid.to_spatial_grid(), x.flatten())

    cpdef GridDataReal add(self, GridDataReal other):
        """Adds values of 'other' to this.

        Does not create any new memory.

        Created to act as accumulator when adding griddatas.

        Args:
            other: GridDataReal object with same grid shape.

        Return:
            This same object, with modified values.
        """

        cdef:
            int i, n
            double[:] other_values

        n = len(self._values)
        other_values = other._values
        for i in range(n):
            self._values[i] += other_values[i]

        return self

    cpdef GridDataReal add_squared(self, GridDataReal other):
        """Adds squared values of 'other' to this.

        Does not create any new memory.

        Created to act as accumulator when adding griddatas.

        Args:
            other: GridDataReal object with same grid shape.

        Return:
            This same object, with modified values.
        """

        cdef:
            int i, n
            double[:] other_values

        n = len(self._values)
        other_values = other._values
        for i in range(n):
            self._values[i] += (other_values[i] * other_values[i])

        return self

    @classmethod
    def sum_all(cls, griddatas):
        """Return sum of griddatas.

        Sums are computed in-place.

        Args:
            griddatas: List of GridDataReal objects, with same grids.

        Return:
            GridDataReal object.
        """

        g = griddatas[0]
        accumulator = cls(g.grid, np.zeros(g.values.shape))

        return reduce(lambda x, y: x.add(y), griddatas, accumulator)

    @classmethod
    def sum_squares(cls, griddatas):
        """Return sum of squares of griddatas.

        Sums are computed in-place.

        Args:
            griddatas: List of GridDataReal objects, with same grids.

        Return:
            GridDataReal object.
        """

        g = griddatas[0]
        accumulator = cls(g.grid, np.zeros(g.values.shape))

        return reduce(lambda x, y: x.add_squared(y), griddatas, accumulator)

    @classmethod
    def average(cls, griddatas, weights=None):
        """Return average of griddatas.

        Averages are computed in-place.

        Args:
            griddatas: List of GridDataReal objects, with same grids.
            weights: Double array, same length as num objects.

        Return:
            GridDataReal object.
        """

        if weights is not None:
            w = weights / weights.sum()
            g = griddatas[0]
            accumulator = cls(g.grid, np.zeros(g.values.shape))
            s = reduce(lambda x, y: x.add(y[0].scale(float(y[1]))),
                       zip(griddatas, w), accumulator)
            return s
        else:
            return cls.sum_all(griddatas).scale(1.0 / len(griddatas))

    @classmethod
    def std(cls, griddatas, average=None):
        """Return voxel-wise standard deviation of griddatas.

        """

        if not average:
            average = GridDataReal.average(griddatas)

        minus_average = GridDataReal(
            average.grid, -np.array(average.values))

        # Subtract average from each griddata.
        centered = [gd.add_slow(minus_average) for gd in griddatas]

        # Compute mean sum of squares of centered griddatas.
        variance = cls.sum_squares(centered).scale(1.0 / len(centered))

        # Compute square roots.
        sqrt_gd = GridDataReal(
            variance.grid, np.sqrt(np.array(variance.values)))

        return sqrt_gd

    def add_slow(self, other):
        """Return the sum of this and other as new GridDataReal object.

        Args:
            other: GridDataReal object with same grid shape.

        Return:
            New GridDataReal object, element-wise sum.
        """

        return GridDataReal(
            self.grid, np.array(self.values) + np.array(other.values))

    def scale(self, alpha):
        """Return GridDataReal with scaled values."""

        return GridDataReal(self.grid, np.array(self.values) * alpha)

    def center_mass(self):
        """Set the grid origin so that the center of mass is at (0,0,0)."""

        grid = self.grid
        p = grid.points # (N, 3)
        v = self.values
        c = np.dot(v, p) / np.sum(v) # (3, ), center of mass
        pixel_center = c / np.array(grid.spacing)
        print 'Center of mass (pixels) before: ', pixel_center

        print np.array(grid.origin)

        grid.set_origin(np.array(grid.origin) - c)

        print np.array(grid.origin)

        p = grid.points # (N, 3)
        v = self.values
        c = np.dot(v, p) / np.sum(v) # (3, ), center of mass
        pixel_center = c / np.array(grid.spacing)
        print 'Center of mass (pixels) after: ', pixel_center

    def clip(self, threshold=0.0):
        """Set values below threshold to 0."""

        v = np.array(self.values)
        v[v < threshold] = 0.0

        return GridDataReal(self.grid, v)

    def negate(self):
        """Take negative of values."""

        return GridDataReal(self.grid, -np.array(self.values))

    def resize_centered(self, new_shape):
        """Return GridDataReal with resized grid.

        The grid still covers the same area.

        Args:
            new_shape: (2, ) integer array with dimensions of new images.
        """

        # Create new (smaller) grid.
        g = Grids.resize_centered(self.grid, new_shape)
        # Resize values.
        I = Image.fromarray(self.reshaped())
        I = I.resize(new_shape, Image.ANTIALIAS)
        gd = GridDataReal(g, np.asarray(I).flatten().astype('double'))

        return gd

    cpdef add_noise(self, double std):
        """Add Gaussian noise to values."""

        cdef:
            int i, n

        n = len(self._values)
        for i in range(n):
            self._values[i] += gsl_ran_gaussian(r, 1) * std

    def add_noise_snr(self, snr):
        """Add Gaussian noise with given signal-to-noise ratio."""

        ## signal_variance = self.variance()
        ## noise_variance = signal_variance / snr
        ## noise_std = np.sqrt(noise_variance)
        noise_std = self.noise_std_from_snr(snr)
        # print 'noise std = %1.3f' % noise_std
        self.add_noise(noise_std)

    def noise_std_from_snr(self, snr):
        """Computes the std required for the noise for a given snr.

        Assumes that the current image is clean.
        The snr (signal-to-noise ratio) is defined as the signal variance
        divided by the noise variance.
        """

        signal_variance = self.variance()
        noise_variance = signal_variance / snr

        return np.sqrt(noise_variance)

    def copy(self):
        """Return deep copy."""

        return GridDataReal(self.grid.copy(), np.array(self.values).copy())

    def variance(self):
        """Return variance of values.

        Note that this can be artificially reduced by using a
        larger background (assuming no noise).
        """

        v = np.array(self.values)
        v -= v.mean()
        return np.sum(v * v) / len(v)

    def normalise(self, subtract_mean=True):
        """Return GridDataReal object with mean 0 and norm 1.

        Original object is not modified.
        """

        v = np.array(self.values)
        if subtract_mean:
            v -= v.mean()
        v /= np.sqrt(np.sum(v * v))

        return GridDataReal(self.grid.copy(), v)

    def convolve(self, psf):
        """Convolve with psf. In-place.

        Args:
            psf: GridDataReal object.
        """

        y = self.reshaped()

        F = Cnvarray(psf.reshaped(), y.shape, 'same')

        self._values = F.cnv(y).flatten()

    ## def thick(self, depth):
    ##     """Return 3D griddata from 2D by creating copies

    @classmethod
    def from_stack(cls, stack, spacing):
        """Return list of GridDataReal objects from 3D stack of images.

        Args:
            stack: (s, s, n), for instance (64, 64, 1000)
            spacing: grid spacing, pixel size, in Angstrom
        """

        print stack.shape[:2]
        print spacing
        grid = Grids.create_centered(2, np.array(stack.shape[:2]), spacing)
        gds = []
        for i in range(stack.shape[2]):
            x = stack[:,:,i].astype('double').flatten()
            gds.append(cls(grid, x))

        return gds

    @classmethod
    def from_stack4D(cls, stack, spacing):
        """Return list of GridDataReal objects from 4D stack of images.

        The motivation for the extra dimension is to store defocus groups.

        Args:
            stack: (s, s, m, n), for instance (64, 64, 24, 7)
            spacing: grid spacing, pixel size, in Angstrom
        """

        print stack.shape[:2]
        print spacing
        grid = Grids.create_centered(2, np.array(stack.shape[:2]), spacing)
        gds = []
        for i in range(stack.shape[2]):
            for j in range(stack.shape[3]):
                x = stack[:,:,i,j].astype('double').flatten()
                gds.append(cls(grid, x))

        return gds

    # def is_non_negative(self):
    #     """Return false if the image has negative values."""
        
    #     self.values
        
        



cdef class GridDataInt:
    """Regular grid with integer values at vertices."""

    def __cinit__(self, Grid grid, values):
        """Initialise with Grid object and values on the grid.

        Args:
            grid: Grid object
            values: (N,) int array, values at grid vertices.
        """

        self._grid = grid
        self._values = values
        self._values_sum = np.sum(values)
        self._cumsum = np.cumsum(values)

    @property
    def grid(self):
        return self._grid

    @property
    def values(self):
        return self._values

    @values.setter
    def values(self, values_new):
        self._values = values_new
        # Also update the sum and cumsum?

    @property
    def values_sum(self):
        return self._values_sum

    @property
    def cumsum(self):
        return self._cumsum

    @property
    def dimension(self):
        return self._grid.dimension

    @property
    def size(self):
        return len(self._values)

    @classmethod
    def create_from_function(cls, grid, function):
        """Sample function on grid, use values to initialise.

        Args:
            grid: Grid object
            function: callable object

        Returns:
            GridData object with function evaluated at grid points.
        """

        values = function(grid.points)

        return cls(grid, values)

    @classmethod
    def create_constant(cls, grid, value=0):
        """Initialise grid with constant value.

        Args:
            grid: Grid object.
            value: Constant value at all grid points, default to 0.0.

        Returns:
            GridData object with constant values.
        """

        values = np.ones((grid.size,), dtype='int') * value

        return cls(grid, values)

    @classmethod
    def create_binned(cls, grid, points):
        """Return GridDataInt representing histogram of points on grid.
        
        Args:
            grid: Grid object
            points: (N, D), where D is the dimension of the grid.

        Currently only implemented for 1D.
        """

        if grid.dimension > 1:
            raise ValueError(
                'Can only create binned data for 1D points.')

        dx = grid.spacing[0]
        N = grid.size
        a = grid.origin - dx/2.0
        b = a + dx*N
        values, edges = np.histogram(points, bins=N, range=(a, b))

        return cls(grid, values)

    def reshaped(self):
        """Returns values reshaped according to grid.

        Creates a copy of the data.
        """

        x = np.zeros(self.size, dtype='int')
        x[...] = self._values

        return x.reshape(self._grid.shape)

    def save(self, parent_group, name):
        """Save to group in HDF5 database.

        Assume that the group parent_group[name] doesn't exist.

        Args:
            parent_group: Group within which to create the griddata group.
            name: Name of the new griddata group.
        """

        group = parent_group.create_group(name)
        self._grid.save(group, 'grid')
        dataset = group.create_dataset('values',
                                       (<object>self._values).shape,
                                       np.int)
        dataset[...] = self._values

    @classmethod
    def load(cls, parent_group, name):
        """Loads GridData from group in HDF5 database."""

        group = parent_group[name]
        grid = Grids.load(group, 'grid')
        values = group['values'][:]

        return cls(grid, values)

    def to_real(self):
        """Convert to GridDataReal object."""

        return GridDataReal(self._grid,
                            np.array(self._values).astype('float'))

    def subtract(self, other):
        """Subtract count values of other GridDataInt object from this one.
        """

        return GridDataInt(self._grid,
                           np.array(self.values) - np.array(other.values))

    def scale(self, alpha):
        """Return GridDataReal with scaled values."""

        return GridDataReal(self.grid, np.array(self.values) * alpha)

    def shift_int(self, dx, dy):
        """Return GridDataInt with shifted values.

        First convert dx and dy from Angstrom to pixels.
        Then round to nearest integers, and then shift with those
        values.
        """

        a = np.array(self.reshaped())
        dxi = np.int(np.round(dx / self.grid.spacing[0]))
        dyi = np.int(np.round(dy / self.grid.spacing[1]))

        values = np.roll(np.roll(a, dxi, axis=0), dyi, axis=1).flatten()

        shifted = GridDataInt(self.grid, values)

        return shifted

    def to_histogram(self):
        """Return points for plotting 1D histogram.

        Return:
            (4*N, 2) array, where N is the number of grid points.
        """
        
        N = self.size

        ps = np.zeros([4*N, 2], dtype='float')
        gridpoints = self.grid.points

        dx = self.grid.spacing[0]
        y = 0.0
        for i in range(N):
            x = gridpoints[i] - 0.5*dx
            v = self.values[i]

            ps[(4*i):(4*(i+1)), :] = np.array(
                [[x, y], [x, y+v], [x+dx, y+v], [x+dx, y]])

        return ps

    def clip(self, threshold=0):
        """Set values below threshold to 0.

        Returns a copy with clipped values. Original not changed.
        """

        v = np.array(self.values)
        v[v < threshold] = 0

        return GridDataInt(self.grid, v)


cdef class GridDatasNested:
    """Nested family of GridDataInt objects parametrised by beta.

    Generated from a single GridDataReal object, and can supply
    all the GridDataInt objects that could be constructed from it
    up to a certain maximum count level.
    """

    def __cinit__(self, GridDataReal griddata_real, long C):
        """
        Args:
            griddata_real: Original image, scale arbitrary.
            C: Desired number of counts. Actual max counts will be
              slightly lower.
        """

        self._griddata_real = griddata_real
        self._C = C

        self.initialise()

    def initialise(self):
        """Analyse the given real griddata to simplify later conversion."""
        cdef:
            int point_index

        min_value = np.min(self._griddata_real.values)
        if min_value < 0.0:
            raise ValueError('Image should be non-negative.')

        total_real = np.sum(self._griddata_real._values)
        if (total_real <= 0):
            raise ValueError('Total image value is %f, ' % (total_real) +
                             'should be positive.')

        min_threshold = total_real / self._C
        # 'min_threshold' is the intensity discretisation level to
        # use to obtain roughly 'self._C' counts.
        # Higher levels would yield fewer counts.

        levels = np.empty(self._C, dtype='f')
        point_indices = np.empty(self._C, dtype='int')

        num_points = len(self._griddata_real._values)
        count_index = 0
        for point_index in range(num_points):

            value = self._griddata_real._values[point_index]
            p_levels = int(value / min_threshold)
            # 'p_levels' is the number of beta values at which
            # the number of counts contributed by this point will change.

            for i in range(p_levels):
                levels[count_index] = value / (i+1)
                point_indices[count_index] = point_index
                count_index += 1
        # Now levels is a list of threshold values where there is a change
        # in the number of counts.
        # By sorting the point_indices according to their levels,
        # we obtain the order in which the points should be added
        # to obtain all griddata objects for lower beta values.

        n = count_index
        # Actual maximum counts (usually slightly lower
        # than the specified number of counts 'self._C')
        sorted_point_indices = point_indices[np.argsort(levels[:n])[::-1]]
        # Length 'n'.
        self._max_counts = n

        self._indices = sorted_point_indices

        self._default = self.from_num_counts(self._max_counts)

    @property
    def grid(self):
        return self._griddata_real._grid

    cpdef GridDataInt from_num_counts(self, int num_counts):
        """Returns a GridDataInt object with the specified number of counts.

        Args:
            num_counts: Between 0 and

        This would correspond to beta = num_counts / total_counts.
        """

        cdef:
            int i
            long[::1] values

        values = np.zeros(len(self._griddata_real._values), dtype='int')
        for i in range(num_counts):
            values[self._indices[i]] += 1
        return GridDataInt(self._griddata_real.grid, values)

    cpdef int beta_to_counts(self, double beta):
        """Return the corresponding number of counts for this beta."""

        return int(beta * self._max_counts)

    cpdef long[::1] beta_array_to_counts(self, double[::1] betas):
        """Return the number of counts for each beta."""

        return np.array([self.beta_to_counts(beta) for beta in betas])

    cpdef GridDataInt from_beta(self, double beta):
        """Returns a GridDataInt object at the given inverse temperature.

        Args:
            beta: Inverse temperature, should be between 0 and 1.
        """

        if beta==1.0:
            return self._default
        else:
            num_counts = self.beta_to_counts(beta)
            ## int(beta * self._max_counts)
            return self.from_num_counts(num_counts)


cdef class GridDataList:

    def __cinit__(self, *args):
        """
        Args:
            data: List of GridDataInt or GridDataReal objects.
        """
        self._data = args[0]

    cpdef get(self, index, beta):

        return self._data[index]

    cpdef get_num_projections(self):
        return len(self._data)

    cpdef as_list(self):
        """Returns a list of GridDataInt or GridDataReal objects."""

        return [self.get(i, 1.0)
                for i in range(self.get_num_projections())]

    def to_nested_list(self, C):
        """Converts to GridDataNestedList.

        Args:
            C: Desired number of counts.

        If the objects in the list are GridDataInt, then they're
        first converted to GridDataReal, and then to GridDataNested.
        The latter step is done in the constructor for GridDataNestedList.
        """

        return GridDataNestedList([d.to_real() for d in self._data], C)

    @classmethod
    def load_mrc(cls, mrc_filename, spacing=1.0):
        """Load list of GridDataReal objects from mrc file.

        Args:
            spacing: Pixel size in Angstrom, example: 4.
        """

        data, _spacing, _origin = mrc.read(mrc_filename)

        if np.isnan(np.min(data)):
            raise ValueError('nan values in mrc file')
        data = data.astype('double')
        if np.isnan(np.min(data)):
            raise ValueError('nan values in mrc file')

        nx, ny, n = data.shape
        grid = Grids.create_centered(2, np.array([nx, ny]), spacing)

        gds = [GridDataReal(grid, data[:, :, i].flatten())
               for i in range(n)]

        return cls(gds)

    def to_stack(self):
        """Convert to 3D stack."""

        return np.dstack([d.reshaped() for d in self._data])

    def save_mrc(self, mrc_filename):
        """Save list of GridDataReal/Int objects as stacked mrc file."""

        stack = self.to_stack()
        mrc.write(stack, mrc_filename)

    def clip(self, threshold=0.0):
        """Set values below threshold to 0."""

        self._data = [d.clip(threshold) for d in self._data]


cdef class GridDataNestedList(GridDataList):

    def __cinit__(self, *args):
        """
        Args:
            data: List of GridDataReal objects.
            C: Desired number of counts.
        """
        data_original = args[0]
        C = args[1]
        self._data = [GridDatasNested(d, C) for d in data_original]

    cpdef get(self, index, beta):
        """Return GridDataInt object."""

        return self._data[index].from_beta(beta)

    def to_nested_list(self, C):
        """Return this object, to share interface with GridDataList."""

        return self

    def get_nested(self, index):
        """Return GridDatasNested object."""

        return self._data[index]

    def save_mrc(self, mrc_filename):
        """Save list of GridDataReal/Int objects as stacked mrc file."""

        stack = np.dstack([d.from_beta(1.0).reshaped() for d in self._data])

        mrc.write(stack, mrc_filename)


cdef class GridDataProjector:
    """For projecting 3D grids to 2D grids multiple times, and fast."""

    def __cinit__(self, GridDataReal reference, int shape, double spacing,
                  int order):
        """
        Args:
            reference: 3D object.
            shape: The number of grid points in each direction as an integer.
            spacing: The grid spacing as a double.
            order: 3D spline interpolation order. Between 1 and 3.

        N: Number of 3D grid points.
        grid3D: Transformed 3D grid onto which to interpolate.
        transformed: Transformed coordinates of all 3D grid points.
        transformed_transpose: Just the transpose of 'transformed'.
        evaluated3D: 1D array with interpolated values on rotated 3D grid.
        grid2D: Grid on which to project.
        projection: 2D projection.
        """

        cdef:
            int j
            Grid3D g

        self._reference = reference
        self._shape = shape
        self._spacing = spacing
        self._order = order
        self._N = shape * shape * shape
        self._grid3D = Grids.create_centered(3, shape, spacing)
        self._transformed = np.zeros((self._N, 3), dtype='double')
        self._transformed_transpose = np.zeros((3, self._N), dtype='double')
        self._evaluated3D = np.zeros((self._N, ), dtype='double')
        self._grid2D = Grids.create_centered(2, shape, spacing)
        self._projection = GridDataReal.create_constant(self._grid2D)


        g = self._reference._grid
        self._alpha = np.empty((3, ), dtype='double')
        self._beta = np.empty((3, ), dtype='double')
        for j in range(3):
            self._alpha[j] = 1.0 / g._spacing[j]
            self._beta[j] = (g._shape[j] - 1.0) / 2.0

    cpdef GridDataReal project(self, Transformation3D transformation):
        """Apply given transformation, and project to 2D.

        Args:
            transformation: Rotation and translation applied to 3D object.

        Returns:
            GridDataReal object with 2D projection. The same pointer
            will be returned everytime, just the values will be changed.
            So make a deep copy in between to keep them in memory.
            This way should be faster when writing to disk.
        """

        cdef:
            int i, j, k, index_2D, index_3D
            double alpha_j, beta_j
            double[:] alpha, beta
            double[:, :] transformed, transformed_transpose
            Transformation3D inverse_transformation

        alpha = self._alpha
        beta = self._beta
        transformed = self._transformed
        transformed_transpose = self._transformed_transpose

        # Transform target grid points.
        # Take the inverse to ensure that the actual transformation
        # is applied to the object, not the reference frame.
        inverse_transformation = transformation.invert()
        inverse_transformation.transform(
            self._grid3D.points, self._transformed)

        # First scale and shift them to be in the
        # funny coordinate system needed for 'map_coordinates'.
        for i in range(self._N):
            for j in range(3):
                transformed_transpose[j, i] = (
                    transformed[i, j] * alpha[j] + beta[j])

        # Now interpolate from the existing values to this new grid.
        map_coordinates(self._reference.reshaped(),
                        self._transformed_transpose,
                        np.asarray(self._evaluated3D),
                        order=self._order)

        # Sum over the z-coordinate.
        index_3D = 0
        index_2D = 0
        for i in range(self._shape):
            for j in range(self._shape):
                self._projection._values[index_2D] = 0.0
                for k in range(self._shape):
                    self._projection._values[index_2D] += (
                        self._evaluated3D[index_3D])
                    index_3D += 1
                self._projection._values[index_2D] *= self._grid3D._spacing[2]
                index_2D += 1

        return self._projection
