Splines

Collaboration diagram for Splines:

Functions
double	spatial_Bspline03 (const Matrix1D< double > &r)
double	Bspline03LUT (double x)
double	spatial_Bspline03LUT (const Matrix1D< double > &r)
double	sum_spatial_Bspline03_Grid (const Grid &grid)
double	spatial_Bspline03_proj (const Matrix1D< double > &r, const Matrix1D< double > &u)
void	spatial_Bspline032voxels (const GridVolume &vol_splines, MultidimArray< double > *vol_voxels, int Zdim=0, int Ydim=0, int Xdim=0)

Variables
const int	BSPLINE03_SUBSAMPLING = 2000

Detailed Description

Function Documentation

Bspline03LUT()

double Bspline03LUT

(

double

x

)

Value of a Bspline of order 3 in a Look-Up Table.

Definition at line 30 of file splines.cpp.

{
    static bool firstCall = true;
    static MultidimArray<double> table(BSPLINE03_SUBSAMPLING);
    static const double deltax = 2.0 / BSPLINE03_SUBSAMPLING;
    static const double ideltax = 1.0 / deltax;
    if (firstCall)
    {
        FOR_ALL_ELEMENTS_IN_ARRAY1D(table)
        table(i) = Bspline03(i * deltax);
        firstCall = false;
    }
    auto i = (size_t)round(fabs(x) * ideltax);
    if (i >= XSIZE(table)) return 0;
    else return table(i);
}

spatial_Bspline03()

double spatial_Bspline03 ( const Matrix1D< double > &  r )

inline

Spline value. This function returns the value of a Bspline of order 3 at a given point

Definition at line 44 of file splines.h.

{
    if (-2 <= XX(r) && XX(r) < 2 &&
        -2 <= YY(r) && YY(r) < 2 &&
        -2 <= ZZ(r) && ZZ(r) < 2)
        return Bspline03(XX(r))*Bspline03(YY(r))*Bspline03(ZZ(r));
    else return 0.0;
}

spatial_Bspline032voxels()

void spatial_Bspline032voxels	(	const GridVolume &	vol_splines,
		MultidimArray< double > *	vol_voxels,
		int	Zdim = 0,
		int	Ydim = 0,
		int	Xdim = 0
	)

Splines —> Voxels. The voxel size is defined in the same grid as the spline volume.

However, you might give a size, usually set to 0, i.e., no external size. If no size is provided a size is produced such that all spline centers fit into the output volume.

Definition at line 386 of file splines.cpp.

{
    // Resize and set starting corner .......................................
    if (Zdim == 0 || Ydim == 0 || Xdim == 0)
    {
        Matrix1D<double> size = vol_splines.grid(0).highest -
                                vol_splines.grid(0).lowest;
        Matrix1D<double> corner = vol_splines.grid(0).lowest;
        (*vol_voxels).initZeros(CEIL(ZZ(size)), CEIL(YY(size)), CEIL(XX(size)));
        STARTINGX(*vol_voxels) = FLOOR(XX(corner));
        STARTINGY(*vol_voxels) = FLOOR(YY(corner));
        STARTINGZ(*vol_voxels) = FLOOR(ZZ(corner));
    }
    else
    {
        (*vol_voxels).initZeros(Zdim, Ydim, Xdim);
        (*vol_voxels).setXmippOrigin();
    }

    // Convert each subvolume ...............................................
    for (size_t i = 0; i < vol_splines.VolumesNo(); i++)
    {
        spatial_Bspline032voxels_SimpleGrid(vol_splines(i)(), vol_splines.grid(i),
                                            vol_voxels);
#ifdef DEBUG
        std::cout << "Spline grid no " << i << " stats: ";
        vol_splines(i)().printStats();
        std::cout << std::endl;
        std::cout << "So far vol stats: ";
        (*vol_voxels).printStats();
        std::cout << std::endl;
        VolumeXmipp save;
        save() = *vol_voxels;
        save.write((std::string)"PPPvoxels" + integerToString(i));
#endif
    }

    // Now normalise the resulting volume ..................................
    double inorm = 1.0 / sum_spatial_Bspline03_Grid(vol_splines.grid());
    FOR_ALL_ELEMENTS_IN_ARRAY3D(*vol_voxels)
    A3D_ELEM(*vol_voxels, k, i, j) *= inorm;

    // Set voxels outside interest region to minimum value .................
    double R = vol_splines.grid(0).get_interest_radius();
    if (R != -1)
    {
        double R2 = (R - 6) * (R - 6);

        // Compute minimum value within sphere
        double min_val = A3D_ELEM(*vol_voxels, 0, 0, 0);
        FOR_ALL_ELEMENTS_IN_ARRAY3D(*vol_voxels)
        if (j*j + i*i + k*k <= R2 - 4)
            min_val = XMIPP_MIN(min_val, A3D_ELEM(*vol_voxels, k, i, j));

        // Substitute minimum value
        R2 = (R - 2) * (R - 2);
        FOR_ALL_ELEMENTS_IN_ARRAY3D(*vol_voxels)
        if (j*j + i*i + k*k >= R2)
            A3D_ELEM(*vol_voxels, k, i, j) = min_val;
    }
}

spatial_Bspline03_proj()

double spatial_Bspline03_proj	(	const Matrix1D< double > &	r,
		const Matrix1D< double > &	u
	)

Spline projection. This function returns the value of the spline line integral through a straight line which passes through the point r with direction u

Definition at line 103 of file splines.cpp.

{
    // Avoids divisions by zero and allows orthogonal rays computation
    static Matrix1D<double> ur(3);
    if (XX(u) == 0) XX(ur) = XMIPP_EQUAL_ACCURACY;
    else XX(ur) = XX(u);
    if (YY(u) == 0) YY(ur) = XMIPP_EQUAL_ACCURACY;
    else YY(ur) = YY(u);
    if (ZZ(u) == 0) ZZ(ur) = XMIPP_EQUAL_ACCURACY;
    else ZZ(ur) = ZZ(u);

    // Some precalculated variables
    double x_sign = SGN(XX(ur));
    double y_sign = SGN(YY(ur));
    double z_sign = SGN(ZZ(ur));

    // Compute the minimum and maximum alpha for the ray
    double alpha_xmin = (-2 - XX(r)) / XX(ur);
    double alpha_xmax = (2 - XX(r)) / XX(ur);
    double alpha_ymin = (-2 - YY(r)) / YY(ur);
    double alpha_ymax = (2 - YY(r)) / YY(ur);
    double alpha_zmin = (-2 - ZZ(r)) / ZZ(ur);
    double alpha_zmax = (2 - ZZ(r)) / ZZ(ur);

    double alpha_min = XMIPP_MAX(XMIPP_MIN(alpha_xmin, alpha_xmax),
                           XMIPP_MIN(alpha_ymin, alpha_ymax));
    alpha_min = XMIPP_MAX(alpha_min, XMIPP_MIN(alpha_zmin, alpha_zmax));
    double alpha_max = XMIPP_MIN(XMIPP_MAX(alpha_xmin, alpha_xmax),
                           XMIPP_MAX(alpha_ymin, alpha_ymax));
    alpha_max = XMIPP_MIN(alpha_max, XMIPP_MAX(alpha_zmin, alpha_zmax));
    if (alpha_max - alpha_min < XMIPP_EQUAL_ACCURACY) return 0.0;

#ifdef DEBUG
    std::cout << "Pixel:  " << r.transpose() << std::endl
    << "Dir:    " << ur.transpose() << std::endl
    << "Alpha x:" << alpha_xmin << " " << alpha_xmax << std::endl
    << "        " << (r + alpha_xmin*ur).transpose() << std::endl
    << "        " << (r + alpha_xmax*ur).transpose() << std::endl
    << "Alpha y:" << alpha_ymin << " " << alpha_ymax << std::endl
    << "        " << (r + alpha_ymin*ur).transpose() << std::endl
    << "        " << (r + alpha_ymax*ur).transpose() << std::endl
    << "Alpha z:" << alpha_zmin << " " << alpha_zmax << std::endl
    << "        " << (r + alpha_zmin*ur).transpose() << std::endl
    << "        " << (r + alpha_zmax*ur).transpose() << std::endl
    << "alpha  :" << alpha_min  << " " << alpha_max  << std::endl
    << std::endl;
#endif

    // Compute the first point in the volume intersecting the ray
    static Matrix1D<double> v(3);
    V3_BY_CT(v, ur, alpha_min);
    V3_PLUS_V3(v, r, v);

#ifdef DEBUG
    std::cout << "First entry point: " << v.transpose() << std::endl;
    std::cout << "   Alpha_min: " << alpha_min << std::endl;
#endif

    // Follow the ray
    double alpha = alpha_min;
    double ray_sum = 0;
    do
    {
        double alpha_x = (XX(v) + x_sign - XX(r)) / XX(ur);
        double alpha_y = (YY(v) + y_sign - YY(r)) / YY(ur);
        double alpha_z = (ZZ(v) + z_sign - ZZ(r)) / ZZ(ur);

        // Which dimension will ray move next step into?, it isn't necessary to be only
        // one.
        double diffx = ABS(alpha - alpha_x);
        double diffy = ABS(alpha - alpha_y);
        double diffz = ABS(alpha - alpha_z);
        double diff_alpha = XMIPP_MIN(XMIPP_MIN(diffx, diffy), diffz);
        ray_sum += spatial_Bspline03_integral(r, ur, alpha, alpha + diff_alpha);

        // Update alpha and the next entry point
        if (ABS(diff_alpha - diffx) <= XMIPP_EQUAL_ACCURACY) alpha = alpha_x;
        if (ABS(diff_alpha - diffy) <= XMIPP_EQUAL_ACCURACY) alpha = alpha_y;
        if (ABS(diff_alpha - diffz) <= XMIPP_EQUAL_ACCURACY) alpha = alpha_z;
        XX(v) += diff_alpha * XX(ur);
        YY(v) += diff_alpha * YY(ur);
        ZZ(v) += diff_alpha * ZZ(ur);

#ifdef DEBUG
        std::cout << "Alpha x,y,z: " << alpha_x << " " << alpha_y
        << " " << alpha_z << " ---> " << alpha << std::endl;

        std::cout << "    Next entry point: " << v.transpose() << std::endl
        << "    diff_alpha: " << diff_alpha << std::endl
        << "    ray_sum: " << ray_sum << std::endl
        << "    Alfa tot: " << alpha << "alpha_max: " << alpha_max <<
        std::endl;

#endif
    }
    while ((alpha_max - alpha) > XMIPP_EQUAL_ACCURACY);
    return ray_sum;
}

spatial_Bspline03LUT()

double spatial_Bspline03LUT ( const Matrix1D< double > &  r )

inline

Spline value with Look-Up Table. This function returns the value of a Bspline of order 3 at a given point

Definition at line 60 of file splines.h.

{
    if (-2 <= XX(r) && XX(r) < 2 &&
        -2 <= YY(r) && YY(r) < 2 &&
        -2 <= ZZ(r) && ZZ(r) < 2)
        return Bspline03LUT(XX(r))*Bspline03LUT(YY(r))*Bspline03LUT(ZZ(r));
    else return 0.0;
}

sum_spatial_Bspline03_Grid()

double sum_spatial_Bspline03_Grid

(

const Grid &

grid

)

Sum of a single spline over a grid. As a normalisation factor, the sum of the splines values over a given grid is needed. This function puts a spline at coordinate (0,0,0) and sums all the splines values at points of the grid which are inside the spline. It doesn't matter if the grid is compound or not.

Definition at line 72 of file splines.cpp.

{
    double sum = 0;
    for (size_t i = 0; i < grid.GridsNo(); i++)
        sum += sum_spatial_Bspline03_SimpleGrid(grid(i));
    return sum;
}

Variable Documentation

BSPLINE03_SUBSAMPLING

const int BSPLINE03_SUBSAMPLING = 2000

Definition at line 53 of file splines.h.