projection.cpp File Reference

#include "projection.h"
#include "core/geometry.h"
#include "core/metadata_vec.h"
#include "core/matrix2d.h"
#include "core/xmipp_image.h"
#include "core/xmipp_macros.h"
#include "core/xmipp_program.h"
#include "core/transformations.h"
#include "data/fourier_projection.h"
#include "data/basis.h"

Include dependency graph for projection.cpp:

Go to the source code of this file.

Macros
#define	x0   STARTINGX(IMGMATRIX(proj))
#define	xF   FINISHINGX(IMGMATRIX(proj))
#define	y0   STARTINGY(IMGMATRIX(proj))
#define	yF   FINISHINGY(IMGMATRIX(proj))
#define	xDim   XSIZE(IMGMATRIX(proj))
#define	yDim   YSIZE(IMGMATRIX(proj))
#define	wrap_as_Crystal(x, y, xw, yw)
#define	xc   XX(rc)
#define	yc   YY(rc)
#define	x   XX(r)
#define	y   YY(r)

Functions
void	projectVolume (MultidimArray< double > &V, Projection &P, int Ydim, int Xdim, double rot, double tilt, double psi, const Matrix1D< double > *roffset)
void	projectVolumeOffCentered (MultidimArray< double > &V, Projection &P, int Ydim, int Xdim)
void	singleWBP (MultidimArray< double > &V, Projection &P)
void	project_Crystal_SimpleGrid (Image< double > &vol, const SimpleGrid &grid, const Basis &basis, Projection &proj, Projection &norm_proj, const Matrix1D< double > &shift, const Matrix1D< double > &aint, const Matrix1D< double > &bint, const Matrix2D< double > &D, const Matrix2D< double > &Dinv, const MultidimArray< int > &mask, int FORW, int eq_mode)
void	project_Crystal_Volume (GridVolume &vol, const Basis &basis, Projection &proj, Projection &norm_proj, int Ydim, int Xdim, double rot, double tilt, double psi, const Matrix1D< double > &shift, const Matrix1D< double > &aint, const Matrix1D< double > &bint, const Matrix2D< double > &D, const Matrix2D< double > &Dinv, const MultidimArray< int > &mask, int FORW, int eq_mode)
void	count_eqs_in_projection (GridVolumeT< int > &GVNeq, const Basis &basis, Projection &read_proj)
template<class T >
void *	project_SimpleGridThread (void *params)
template<class T >
void	project_SimpleGrid (Image< T > *vol, const SimpleGrid *grid, const Basis *basis, Projection *proj, Projection *norm_proj, int FORW, int eq_mode, const Image< int > *VNeq, Matrix2D< double > *M, const MultidimArray< int > *mask, double ray_length, int thread_id, int numthreads)
template<class T >
void	project_GridVolume (GridVolumeT< T > &vol, const Basis &basis, Projection &proj, Projection &norm_proj, int Ydim, int Xdim, double rot, double tilt, double psi, int FORW, int eq_mode, GridVolumeT< int > *GVNeq, Matrix2D< double > *M, const MultidimArray< int > *mask, double ray_length, int threads)
template void *	project_SimpleGridThread< double > (void *)
template void	project_GridVolume< double > (GridVolumeT< double > &, Basis const &, Projection &, Projection &, int, int, double, double, double, int, int, GridVolumeT< int > *, Matrix2D< double > *, MultidimArray< int > const *, double, int)

Variables
barrier_t	project_barrier
pthread_mutex_t	project_mutex = PTHREAD_MUTEX_INITIALIZER
project_thread_params *	project_threads

Macro Definition Documentation

wrap_as_Crystal

#define wrap_as_Crystal	(		x,
			y,
			xw,
			yw
	)

Value:

xw=(int) intWRAP(x,x0,xF); \
    yw=(int) intWRAP(y,y0,yF);

Definition at line 769 of file projection.cpp.

x

double * x   XX(r)

Definition at line 2473 of file numerical_recipes.cpp.

x0

#define x0   STARTINGX(IMGMATRIX(proj))

Definition at line 37 of file projection.cpp.

xc

#define xc   XX(rc)

xDim

#define xDim   XSIZE(IMGMATRIX(proj))

Definition at line 41 of file projection.cpp.

xF

#define xF   FINISHINGX(IMGMATRIX(proj))

Definition at line 38 of file projection.cpp.

y

static void y   YY(r)

Definition at line 8508 of file numerical_recipes.cpp.

y0

#define y0   STARTINGY(IMGMATRIX(proj))

Definition at line 39 of file projection.cpp.

yc

#define yc   YY(rc)

yDim

#define yDim   YSIZE(IMGMATRIX(proj))

Definition at line 42 of file projection.cpp.

yF

#define yF   FINISHINGY(IMGMATRIX(proj))

Definition at line 40 of file projection.cpp.

Function Documentation

project_Crystal_SimpleGrid()

void project_Crystal_SimpleGrid	(	Image< double > &	vol,
		const SimpleGrid &	grid,
		const Basis &	basis,
		Projection &	proj,
		Projection &	norm_proj,
		const Matrix1D< double > &	shift,
		const Matrix1D< double > &	aint,
		const Matrix1D< double > &	bint,
		const Matrix2D< double > &	D,
		const Matrix2D< double > &	Dinv,
		const MultidimArray< int > &	mask,
		int	FORW,
		int	eq_mode
	)

Definition at line 773 of file projection.cpp.

{
    Matrix1D<double> prjX(3);                // Coordinate: Projection of the
    Matrix1D<double> prjY(3);                // 3 grid vectors
    Matrix1D<double> prjZ(3);
    Matrix1D<double> prjOrigin(3);           // Coordinate: Where in the 2D
    // projection plane the origin of
    // the grid projects
    Matrix1D<double> prjaint(3);             // Coordinate: Projection of the
    Matrix1D<double> prjbint(3);             // 2 deformed lattice vectors
    Matrix1D<double> prjDir;                 // Direction of projection
    // in the deformed space

    Matrix1D<double> actprj(3);              // Coord: Actual position inside
    // the projection plane
    Matrix1D<double> defactprj(3);           // Coord: Actual position inside
    // the deformed projection plane
    Matrix1D<double> beginZ(3);              // Coord: Plane coordinates of the
    // projection of the 3D point
    // (z0,YY(lowest),XX(lowest))
    Matrix1D<double> beginY(3);              // Coord: Plane coordinates of the
    // projection of the 3D point
    // (z0,y0,XX(lowest))
    Matrix1D<double> footprint_size(2);      // The footprint is supposed
    // to be defined between
    // (-vmax,+vmax) in the Y axis,
    // and (-umax,+umax) in the X axis
    // This footprint size is the
    // R2 vector (umax,vmax).
    Matrix1D<double> deffootprint_size(2);   // The footprint size
    // in the deformed space
    int XX_corner2;
    int XX_corner1;                          // Coord: Corners of the
    int YY_corner2;
    int YY_corner1;                          // footprint when it is projected
    // onto the projection plane
    Matrix1D<double> rc(2);
    Matrix1D<double> r(2);                   // Position vector which will
    // move from corner1 to corner2.
    // In rc the wrapping is not
    // considered, while it is in r
    int           foot_V;
    int           foot_U;                    // Img Coord: coordinate
    // corresponding to the blobprint
    // point which matches with this
    // pixel position
    double        vol_corr=0.;               // Correction for a volume element
    int           N_eq;                      // Number of equations in which
    // a blob is involved
    int           i;
    int           j;
    int           k;                         // volume element indexes
    SPEED_UP_temps01;                        // Auxiliary variables for
    // fast multiplications

    // Check that it is a blob volume .......................................
    if (basis.type != Basis::blobs)
        REPORT_ERROR(ERR_VALUE_INCORRECT, "project_Crystal_SimpleGrid: Cannot project other than "
                     "blob volumes");

    // Compute the deformed direction of projection .........................
    Matrix2D<double> Eulerg;
    Eulerg = proj.euler * D;

    // Compute deformation in the projection space ..........................
    // The following two vectors are defined in the deformed volume space
    VECTOR_R3(actprj, XX(aint), YY(aint), 0);
    grid.Gdir_project_to_plane(actprj, Eulerg, prjaint);
    VECTOR_R3(actprj, XX(bint), YY(bint), 0);
    grid.Gdir_project_to_plane(actprj, Eulerg, prjbint);
    //#define DEBUG_LITTLE
#ifdef DEBUG_LITTLE

    double rot, tilt, psi;
    Euler_matrix2angles(proj.euler, rot, tilt, psi);
    std::cout << "rot= "  << rot << " tilt= " << tilt
    << " psi= " << psi << std::endl;
    std::cout << "D\n" << D << std::endl;
    std::cout << "Eulerf\n" << proj.euler << std::endl;
    std::cout << "Eulerg\n" << Eulerg << std::endl;
    std::cout << "aint    " << aint.transpose() << std::endl;
    std::cout << "bint    " << bint.transpose() << std::endl;
    std::cout << "prjaint " << prjaint.transpose() << std::endl;
    std::cout << "prjbint " << prjbint.transpose() << std::endl;
    std::cout.flush();
#endif
    // Project grid axis ....................................................
    // These vectors ((1,0,0),(0,1,0),...) are referred to the grid
    // coord. system and the result is a 2D vector in the image plane
    // The unit size within the image plane is 1, ie, the same as for
    // the universal grid.
    // Be careful that these grid vectors are defined in the deformed
    // volume space, and the projection are defined in the deformed
    // projections,
    VECTOR_R3(actprj, 1, 0, 0);
    grid.Gdir_project_to_plane(actprj, Eulerg, prjX);
    VECTOR_R3(actprj, 0, 1, 0);
    grid.Gdir_project_to_plane(actprj, Eulerg, prjY);
    VECTOR_R3(actprj, 0, 0, 1);
    grid.Gdir_project_to_plane(actprj, Eulerg, prjZ);

    // This time the origin of the grid is in the universal coord system
    // but in the deformed space
    Uproject_to_plane(grid.origin, Eulerg, prjOrigin);

    // This is a correction used by the crystallographic symmetries
    prjOrigin += XX(shift) * prjaint + YY(shift) * prjbint;

#ifdef DEBUG_LITTLE

    std::cout << "prjX      " << prjX.transpose() << std::endl;
    std::cout << "prjY      " << prjY.transpose() << std::endl;
    std::cout << "prjZ      " << prjZ.transpose() << std::endl;
    std::cout << "prjOrigin " << prjOrigin.transpose() << std::endl;
    std::cout.flush();
#endif

    // Now I will impose that prja becomes (Xdim,0) and prjb, (0,Ydim)
    // A is a matrix such that
    // A*prjaint=(Xdim,0)'
    // A*prjbint=(0,Ydim)'
    Matrix2D<double> A(2, 2);
    Matrix2D<double> Ainv(2, 2);
    A(0, 0) = YY(prjbint) * xDim;
    A(0, 1) = -XX(prjbint) * xDim;
    A(1, 0) = -YY(prjaint) * yDim;
    A(1, 1) = XX(prjaint) * yDim;
    double nor = 1 / (XX(prjaint) * YY(prjbint) - XX(prjbint) * YY(prjaint));
    M2x2_BY_CT(A, A, nor);
    M2x2_INV(Ainv, A);

#ifdef DEBUG_LITTLE

    std::cout << "A\n" << A << std::endl;
    std::cout << "Ainv\n" << Ainv << std::endl;
    std::cout << "Check that A*prjaint=(Xdim,0)    "
    << (A*vectorR2(XX(prjaint), YY(prjaint))).transpose() << std::endl;
    std::cout << "Check that A*prjbint=(0,Ydim)    "
    << (A*vectorR2(XX(prjbint), YY(prjbint))).transpose() << std::endl;
    std::cout << "Check that Ainv*(Xdim,0)=prjaint "
    << (Ainv*vectorR2(xDim, 0)).transpose() << std::endl;
    std::cout << "Check that Ainv*(0,Ydim)=prjbint "
    << (Ainv*vectorR2(0, yDim)).transpose() << std::endl;
    std::cout.flush();
#endif

    // Footprint size .......................................................
    // The following vectors are integer valued vectors, but they are
    // stored as real ones to make easier operations with other vectors.
    // Look out that this footprint size is in the deformed projection,
    // that is why it is not a square footprint but a parallelogram
    // and we have to look for the smaller rectangle which surrounds it
    XX(footprint_size) = basis.blobprint.umax();
    YY(footprint_size) = basis.blobprint.vmax();
    N_eq = (2 * basis.blobprint.umax() + 1) * (2 * basis.blobprint.vmax() + 1);
    Matrix1D<double> c1(3);
    Matrix1D<double> c2(3);
    XX(c1)           = XX(footprint_size);
    YY(c1)           = YY(footprint_size);
    XX(c2)           = -XX(c1);
    YY(c2)           = YY(c1);
    M2x2_BY_V2x1(c1, A, c1);
    M2x2_BY_V2x1(c2, A, c2);
    XX(deffootprint_size) = fmax(fabs(XX(c1)), fabs(XX(c2)));
    YY(deffootprint_size) = fmax(fabs(YY(c1)), fabs(YY(c2)));

#ifdef DEBUG_LITTLE

    std::cout << "Footprint_size " << footprint_size.transpose() << std::endl;
    std::cout << "c1: " << c1.transpose() << std::endl;
    std::cout << "c2: " << c2.transpose() << std::endl;
    std::cout << "Deformed Footprint_size " << deffootprint_size.transpose()
    << std::endl;
    std::cout.flush();
#endif

    // This type conversion gives more speed
    auto ZZ_lowest = (int) ZZ(grid.lowest);
    int YY_lowest = XMIPP_MAX((int) YY(grid.lowest), STARTINGY(mask));
    int XX_lowest = XMIPP_MAX((int) XX(grid.lowest), STARTINGX(mask));
    auto ZZ_highest = (int) ZZ(grid.highest);
    int YY_highest = XMIPP_MIN((int) YY(grid.highest), FINISHINGY(mask));
    int XX_highest = XMIPP_MIN((int) XX(grid.highest), FINISHINGX(mask));

    // Project the whole grid ...............................................
    // Corner of the plane defined by Z. These coordinates try to be within
    // the valid indexes of the projection (defined between STARTING and
    // FINISHING values, but it could be that they may lie outside.
    // These coordinates need not to be integer, in general, they are
    // real vectors.
    // The vectors returned by the projection routines are R3 but we
    // are only interested in their first 2 components, ie, in the
    // in-plane components
    beginZ = (double)XX_lowest * prjX + (double)YY_lowest * prjY + (double)ZZ_lowest * prjZ + prjOrigin;

#ifdef DEBUG_LITTLE

    std::cout << "BeginZ     " << beginZ.transpose()             << std::endl;
    std::cout << "Mask       ";
    mask.printShape();
    std::cout << std::endl;
    std::cout << "Vol        ";
    vol().printShape();
    std::cout << std::endl;
    std::cout << "Proj       ";
    proj().printShape();
    std::cout << std::endl;
    std::cout << "Norm Proj  ";
    norm_proj().printShape();
    std::cout << std::endl;
    //std::cout << "Footprint  ";
    //footprint().printShape();
    //std::cout << std::endl;
    //std::cout << "Footprint2 ";
    //footprint2().printShape();
    //std::cout << std::endl;
#endif

    Matrix1D<double> grid_index(3);
    for (k = ZZ_lowest; k <= ZZ_highest; k++)
    {
        // Corner of the row defined by Y
        beginY = beginZ;
        for (i = YY_lowest; i <= YY_highest; i++)
        {
            // First point in the row
            actprj = beginY;
            for (j = XX_lowest; j <= XX_highest; j++)
            {
                VECTOR_R3(grid_index, j, i, k);
#ifdef DEBUG

                std::cout << "Visiting " << i << " " << j << std::endl;
#endif

                // Be careful that you cannot skip any blob, although its
                // value be 0, because it is useful for norm_proj
                // unless it doesn't belong to the reconstruction mask
                if (A2D_ELEM(mask, i, j) && grid.is_interesting(grid_index))
                {
                    // Look for the position in the deformed projection
                    M2x2_BY_V2x1(defactprj, A, actprj);

                    // Search for integer corners for this blob
                    XX_corner1 = CEIL(XX(defactprj) - XX(deffootprint_size));
                    YY_corner1 = CEIL(YY(defactprj) - YY(deffootprint_size));
                    XX_corner2 = FLOOR(XX(defactprj) + XX(deffootprint_size));
                    YY_corner2 = FLOOR(YY(defactprj) + YY(deffootprint_size));

#ifdef DEBUG

                    std::cout << "  k= " << k << " i= " << i << " j= " << j << std::endl;
                    std::cout << "  Actual position: " << actprj.transpose() << std::endl;
                    std::cout << "  Deformed position: " << defactprj.transpose() << std::endl;
                    std::cout << "  Blob corners: (" << XX_corner1 << "," << YY_corner1
                    << ") (" << XX_corner2 << "," << YY_corner2 << ")\n";
                    std::cout.flush();
#endif

                    // Now there is no way that both corners collapse into a line,
                    // since the projection points are wrapped

                    if (!FORW)
                        vol_corr = 0;

                    // Effectively project this blob
                    // (xc,yc) is the position of the considered pixel
                    // in the crystal undeformed projection
                    // When wrapping and deformation are considered then x and
                    // y are used
#define xc XX(rc)
#define yc YY(rc)
#define x  XX(r)
#define y  YY(r)

                    for (y = YY_corner1; y <= YY_corner2; y++)
                        for (x = XX_corner1; x <= XX_corner2; x++)
                        {
                            // Compute position in undeformed space and
                            // corresponding sample in the blob footprint
                            M2x2_BY_V2x1(rc, Ainv, r);
                            OVER2IMG(basis.blobprint, yc - YY(actprj), xc - XX(actprj),
                                     foot_V, foot_U);

#ifdef DEBUG

                            std::cout << "    Studying: " << r.transpose()
                            << "--> " << rc.transpose()
                            << "dist (" << xc - XX(actprj)
                            << "," << yc - YY(actprj) << ") --> "
                            << foot_U << " " << foot_V << std::endl;
                            std::cout.flush();
#endif

                            if (IMGMATRIX(basis.blobprint).outside(foot_V, foot_U))
                                continue;

                            // Wrap positions
                            int yw;
                            int xw;
                            wrap_as_Crystal(x, y, xw, yw);
#ifdef DEBUG

                            std::cout << "      After wrapping " << xw << " " << yw << std::endl;
                            std::cout << "      Value added = " << VOLVOXEL(vol, k, i, j) *
                            IMGPIXEL(basis.blobprint, foot_V, foot_U) << " Blob value = "
                            << IMGPIXEL(basis.blobprint, foot_V, foot_U)
                            << " Blob^2 " << IMGPIXEL(basis.blobprint2, foot_V, foot_U)
                            << std::endl;
                            std::cout.flush();
#endif

                            if (FORW)
                            {
                                IMGPIXEL(proj, yw, xw) += VOLVOXEL(vol, k, i, j) *
                                                          IMGPIXEL(basis.blobprint, foot_V, foot_U);
                                switch (eq_mode)
                                {
                                case CAVARTK:
                                case ARTK:
                                    IMGPIXEL(norm_proj, yw, xw) +=
                                        IMGPIXEL(basis.blobprint2, foot_V, foot_U);
                                    break;
                                }
                            }
                            else
                            {
                                vol_corr += IMGPIXEL(norm_proj, yw, xw) *
                                            IMGPIXEL(basis.blobprint, foot_V, foot_U);
                            }

                        }

#ifdef DEBUG_INTERMIDIATE
                    proj.write("inter.xmp");
                    norm_proj.write("inter_norm.xmp");
                    std::cout << "Press any key\n";
                    char c;
                    std::cin >> c;
#endif

                    if (!FORW)
                        switch (eq_mode)
                        {
                        case ARTK:
                            VOLVOXEL(vol, k, i, j) += vol_corr;
                            break;
                        case CAVARTK:
                            VOLVOXEL(vol, k, i, j) += vol_corr / N_eq;
                            break;
                        }
                }

                // Prepare for next iteration
                V2_PLUS_V2(actprj, actprj, prjX);
            }
            V2_PLUS_V2(beginY, beginY, prjY);
        }
        V2_PLUS_V2(beginZ, beginZ, prjZ);
    }
    //#define DEBUG_AT_THE_END
#ifdef DEBUG_AT_THE_END
    proj.write("inter.xmp");
    norm_proj.write("inter_norm.xmp");
    std::cout << "Press any key\n";
    char c;
    std::cin >> c;
#endif

#undef xc
#undef yc
#undef x
#undef y
#undef wrap_as_Crystal
}

project_GridVolume< double >()

template void project_GridVolume< double >	(	GridVolumeT< double > &	,
		Basis const &	,
		Projection &	,
		Projection &	,
		int	,
		int	,
		double	,
		double	,
		double	,
		int	,
		int	,
		GridVolumeT< int > *	,
		Matrix2D< double > *	,
		MultidimArray< int > const *	,
		double	,
		int
	)

project_SimpleGridThread< double >()

template void* project_SimpleGridThread< double >

(

void *

)