splines.cpp File Reference

#include "splines.h"
#include "integration.h"

Include dependency graph for splines.cpp:

Go to the source code of this file.

Macros
#define	x0   STARTINGX(*vol_voxels)
#define	xF   FINISHINGX(*vol_voxels)
#define	y0   STARTINGY(*vol_voxels)
#define	yF   FINISHINGY(*vol_voxels)
#define	z0   STARTINGZ(*vol_voxels)
#define	zF   FINISHINGZ(*vol_voxels)

Functions
double	Bspline03LUT (double x)
double	sum_spatial_Bspline03_SimpleGrid (const SimpleGrid &grid)
double	sum_spatial_Bspline03_Grid (const Grid &grid)
double	spatial_Bspline03_integralf (double alpha)
double	spatial_Bspline03_integral (const Matrix1D< double > &r, const Matrix1D< double > &u, double alpha0, double alphaF)
double	spatial_Bspline03_proj (const Matrix1D< double > &r, const Matrix1D< double > &u)
void	spatial_Bspline032voxels_SimpleGrid (const MultidimArray< double > &vol_splines, const SimpleGrid &grid, MultidimArray< double > *vol_voxels, const MultidimArray< double > *vol_mask=nullptr)
void	spatial_Bspline032voxels (const GridVolume &vol_splines, MultidimArray< double > *vol_voxels, int Zdim, int Ydim, int Xdim)

Macro Definition Documentation

x0

#define x0   STARTINGX(*vol_voxels)

xF

#define xF   FINISHINGX(*vol_voxels)

y0

#define y0   STARTINGY(*vol_voxels)

yF

#define yF   FINISHINGY(*vol_voxels)

z0

#define z0   STARTINGZ(*vol_voxels)

zF

#define zF   FINISHINGZ(*vol_voxels)

Function Documentation

spatial_Bspline032voxels_SimpleGrid()

void spatial_Bspline032voxels_SimpleGrid	(	const MultidimArray< double > &	vol_splines,
		const SimpleGrid &	grid,
		MultidimArray< double > *	vol_voxels,
		const MultidimArray< double > *	vol_mask = nullptr
	)

Definition at line 205 of file splines.cpp.

{
    Matrix1D<double> act_coord(3);           // Coord: Actual position inside
    // the voxel volume without deforming
    Matrix1D<double> beginZ(3);              // Coord: Voxel coordinates of the
    // blob at the 3D point
    // (z0,YY(lowest),XX(lowest))
    Matrix1D<double> beginY(3);              // Coord: Voxel coordinates of the
    // blob at the 3D point
    // (z0,y0,XX(lowest))
    Matrix1D<int> corner2(3), corner1(3);    // Coord: Corners of the blob in the voxel volume
    Matrix1D<double> gcurrent(3);            // Position in g of current point
    int           i, j, k;                   // Index within the blob volume
    bool          process;                   // True if this blob has to be
    // processed
    double spline_radius = 2 * sqrt(3.0);

    // Some aliases
#define x0 STARTINGX(*vol_voxels)
#define xF FINISHINGX(*vol_voxels)
#define y0 STARTINGY(*vol_voxels)
#define yF FINISHINGY(*vol_voxels)
#define z0 STARTINGZ(*vol_voxels)
#define zF FINISHINGZ(*vol_voxels)

#ifdef DEBUG
    bool condition = true;
    (*vol_voxels)().printShape();
    std::cout << std::endl;
    std::cout << "x0= " << x0 << " xF= " << xF << std::endl;
    std::cout << "y0= " << y0 << " yF= " << yF << std::endl;
    std::cout << "z0= " << z0 << " zF= " << zF << std::endl;
    std::cout << grid;
#endif

    // Convert the whole grid ...............................................
    // Corner of the plane defined by Z. These coordinates are in the
    // universal coord. system
    grid.grid2universe(grid.lowest, beginZ);

    Matrix1D<double> grid_index(3);
    for (k = (int) ZZ(grid.lowest); k <= (int) ZZ(grid.highest); k++)
    {
        // Corner of the row defined by Y
        beginY = beginZ;
        for (i = (int) YY(grid.lowest); i <= (int) YY(grid.highest); i++)
        {
            // First point in the row
            act_coord = beginY;
            for (j = (int) XX(grid.lowest); j <= (int) XX(grid.highest); j++)
            {
                VECTOR_R3(grid_index, j, i, k);
#ifdef DEBUG
                if (condition)
                {
                    printf("Dealing spline at (%d,%d,%d) = %f\n", j, i, k, A3D_ELEM(vol_splines, k, i, j));
                    std::cout << "Center of the blob      "
                    << act_coord.transpose() << std::endl;
                }
#endif

                // These two corners are also real valued
                process = true;
                if (XX(act_coord) >= xF) process = false;
                if (YY(act_coord) >= yF) process = false;
                if (ZZ(act_coord) >= zF) process = false;
                if (XX(act_coord) <= x0) process = false;
                if (YY(act_coord) <= y0) process = false;
                if (ZZ(act_coord) <= z0) process = false;
#ifdef DEBUG
                if (!process && condition) std::cout << "   It is outside output volume\n";
#endif
                if (!grid.is_interesting(act_coord))
                {
#ifdef DEBUG
                    if (process && condition) std::cout << "   It is not interesting\n";
#endif
                    process = false;
                }

                if (process)
                {
                    V3_PLUS_CT(corner1, act_coord, -spline_radius);
                    V3_PLUS_CT(corner2, act_coord, spline_radius);
#ifdef DEBUG
                    if (condition)
                    {
                        std::cout << "Corner 1 for this point " << corner1.transpose() << std::endl;
                        std::cout << "Corner 2 for this point " << corner2.transpose() << std::endl;
                    }
#endif

                    // Clip the corners to the volume borders
                    XX(corner1) = ROUND(CLIP(XX(corner1), x0, xF));
                    YY(corner1) = ROUND(CLIP(YY(corner1), y0, yF));
                    ZZ(corner1) = ROUND(CLIP(ZZ(corner1), z0, zF));
                    XX(corner2) = ROUND(CLIP(XX(corner2), x0, xF));
                    YY(corner2) = ROUND(CLIP(YY(corner2), y0, yF));
                    ZZ(corner2) = ROUND(CLIP(ZZ(corner2), z0, zF));
#ifdef DEBUG
                    if (condition)
                    {
                        std::cout << "Clipped and rounded Corner 1 " << corner1.transpose() << std::endl;
                        std::cout << "Clipped and rounded Corner 2 " << corner2.transpose() << std::endl;
                    }
#endif

                    // Effectively convert
                    for (int iz = ZZ(corner1); iz <= ZZ(corner2); iz++)
                    {
                        double intz=iz;
                        for (int iy = YY(corner1); iy <= YY(corner2); iy++)
                        {
                            double inty=iy;
                            for (int ix = XX(corner1); ix <= XX(corner2); ix++)
                            {
                                if (vol_mask != nullptr)
                                    if (!A3D_ELEM(*vol_mask, iz, iy, ix)) continue;
                                double intx=ix;

                                // Compute the spline value at this point
                                VECTOR_R3(gcurrent, static_cast<double>(intx),
                                static_cast<double>(inty), static_cast<double>(intz));
                                V3_MINUS_V3(gcurrent, act_coord, gcurrent);
                                double spline_value = spatial_Bspline03(gcurrent);
#ifdef DEBUG_MORE
                                if (condition)
                                {
                                    std::cout << "At (" << intx << ","
                                    << inty << "," << intz << ") value="
                                    << spline_value;
                                    std::cout.flush();
                                }
#endif

                                // Add at that position the corresponding spline value
                                A3D_ELEM(*vol_voxels, iz, iy, ix) +=
                                    A3D_ELEM(vol_splines, k, i, j) * spline_value;
#ifdef DEBUG_MORE
                                if (condition)
                                {
                                    std::cout << " adding " << A3D_ELEM(vol_splines, k, i, j)
                                    << " * " << value_spline << " = "
                                    << A3D_ELEM(vol_splines, k, i, j)*
                                    value_spline << std::endl;
                                    std::cout.flush();
                                }
#endif
                            }
                        }
                    }
                }

                // Prepare for next iteration
                XX(act_coord) = XX(act_coord) + grid.relative_size * grid.basis(0, 0);
                YY(act_coord) = YY(act_coord) + grid.relative_size * grid.basis(1, 0);
                ZZ(act_coord) = ZZ(act_coord) + grid.relative_size * grid.basis(2, 0);
            }
            XX(beginY) = XX(beginY) + grid.relative_size * grid.basis(0, 1);
            YY(beginY) = YY(beginY) + grid.relative_size * grid.basis(1, 1);
            ZZ(beginY) = ZZ(beginY) + grid.relative_size * grid.basis(2, 1);
        }
        XX(beginZ) = XX(beginZ) + grid.relative_size * grid.basis(0, 2);
        YY(beginZ) = YY(beginZ) + grid.relative_size * grid.basis(1, 2);
        ZZ(beginZ) = ZZ(beginZ) + grid.relative_size * grid.basis(2, 2);
    }
}

spatial_Bspline03_integral()

double spatial_Bspline03_integral	(	const Matrix1D< double > &	r,
		const Matrix1D< double > &	u,
		double	alpha0,
		double	alphaF
	)

Definition at line 92 of file splines.cpp.

{
    global_r = r;
    global_u = u;
    return integrateNewtonCotes(&spatial_Bspline03_integralf, alpha0, alphaF, 9);
}

spatial_Bspline03_integralf()

double spatial_Bspline03_integralf

(

double

alpha

)

Definition at line 84 of file splines.cpp.

{
    XX(global_aux) = XX(global_r) + alpha * XX(global_u);
    YY(global_aux) = YY(global_r) + alpha * YY(global_u);
    ZZ(global_aux) = ZZ(global_r) + alpha * ZZ(global_u);
    return spatial_Bspline03LUT(global_aux);
}

sum_spatial_Bspline03_SimpleGrid()

double sum_spatial_Bspline03_SimpleGrid

(

const SimpleGrid &

grid

)

Definition at line 48 of file splines.cpp.

{
    Matrix1D<double> gr(3), ur(3), corner1(3), corner2(3);
    int          i, j, k;
    double        sum = 0.0;

// Compute the limits of the spline in the grid coordinate system
    grid.universe2grid(vectorR3(-2.0, -2.0, -2.0), corner1);
    grid.universe2grid(vectorR3(2.0, 2.0, 2.0), corner2);

// Compute the sum in the points inside the grid
// The integer part of the vectors is taken for not picking points
// just in the border of the spline, which we know they are 0.
    for (i = (int)XX(corner1); i <= (int)XX(corner2); i++)
        for (j = (int)YY(corner1); j <= (int)YY(corner2); j++)
            for (k = (int)ZZ(corner1); k <= (int)ZZ(corner2); k++)
            {
                VECTOR_R3(gr, i, j, k);
                grid.grid2universe(gr, ur);
                sum += spatial_Bspline03(ur);
            }
    return sum;
}