Parameters Projection for Tomography

Collaboration diagram for Parameters Projection for Tomography:

Classes
class	ParametersProjectionTomography
struct	project_thread_params

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))

Functions
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=nullptr, double ray_length=-1.0, int thread_id=-1, int num_threads=1)
void	projectVolume (MultidimArray< double > &V, Projection &P, int Ydim, int Xdim, double rot, double tilt, double psi, const Matrix1D< double > *roffset=nullptr)
void	projectVolumeOffCentered (MultidimArray< double > &V, Projection &P, int Ydim, int Xdim)
void	singleWBP (MultidimArray< double > &V, Projection &P)
void	count_eqs_in_projection (GridVolumeT< int > &GVNeq, const Basis &basis, Projection &read_proj)
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=ARTK)
template<class T >
void *	project_SimpleGridThread (void *params)
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)

Variables
barrier_t	project_barrier
pthread_mutex_t	project_mutex
project_thread_params *	project_threads
constexpr int	FORWARD = 1
constexpr int	BACKWARD = 0
constexpr int	ARTK = 1
constexpr int	CAVK = 2
constexpr int	COUNT_EQ = 3
constexpr int	CAV = 4
constexpr int	CAVARTK = 5
const int	ART_PIXEL_SUBSAMPLING = 2

Detailed Description

Macro Definition Documentation

x0

#define x0   STARTINGX(IMGMATRIX(*proj))

Definition at line 255 of file projection.h.

xDim

#define xDim   XSIZE(IMGMATRIX(*proj))

Definition at line 259 of file projection.h.

xF

#define xF   FINISHINGX(IMGMATRIX(*proj))

Definition at line 256 of file projection.h.

y0

#define y0   STARTINGY(IMGMATRIX(*proj))

Definition at line 257 of file projection.h.

yDim

#define yDim   YSIZE(IMGMATRIX(*proj))

Definition at line 260 of file projection.h.

yF

#define yF   FINISHINGY(IMGMATRIX(*proj))

Definition at line 258 of file projection.h.

Function Documentation

count_eqs_in_projection()

void count_eqs_in_projection	(	GridVolumeT< int > &	GVNeq,
		const Basis &	basis,
		Projection &	read_proj
	)

Count equations in volume. For Component AVeraing (CAV), the number of equations in which each basis is involved is needed.

Definition at line 1206 of file projection.cpp.

{
    for (size_t i = 0; i < GVNeq.VolumesNo(); i++)
        project_SimpleGrid(&(GVNeq(i)), &(GVNeq.grid(i)), &basis,
                           &read_proj, &read_proj, FORWARD, COUNT_EQ, NULL, NULL);
}

project_Crystal_Volume()

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 = ARTK
	)

Project a crystal basis volume. This function projects a crystal deformed basis volume, ie, in the documentation volume g. However the angles given must be those for volume f, the undeformed one. You must supply the deformed lattice vectors, and the matrix to pass from the deformed to the undeformed vectors (D and Dinv). a=D*ai;

Valid eq_modes are ARTK, CAVARTK and CAV.

Definition at line 1157 of file projection.cpp.

{
    // If projecting forward initialise projections
    if (FORW)
    {
        proj.reset(Ydim, Xdim);
        proj.setAngles(rot, tilt, psi);
        norm_proj().resize(proj());
    }

    // Project each subvolume
    for (size_t i = 0; i < vol.VolumesNo(); i++)
    {
        project_Crystal_SimpleGrid(vol(i), vol.grid(i), basis,
                                   proj, norm_proj, shift, aint, bint, D, Dinv, mask, FORW, eq_mode);

#ifdef DEBUG

        Image<double> save;
        save = norm_proj;
        if (FORW)
            save.write((std::string)"PPPnorm_FORW" + (char)(48 + i));
        else
            save.write((std::string)"PPPnorm_BACK" + (char)(48 + i));
#endif

    }
}

project_GridVolume()

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
	)

Projection of a Grid Volume.

Project a grid volume with a basis. The Grid volume is projected onto a projection plane defined by (rot, tilt, psi) (1st, 2nd and 3rd Euler angles). The projection is previously is resized to Ydim x Xdim and initialized to 0. The projection itself, from now on, will keep the Euler angles.

FORWARD process: Each volume of the grid is projected on to the projection plane. The output is the projection itself and a normalising image, the normalising image is the projection of the same grid supposing that all basis are of value 1. This normalising image is used by the ART process

BACKWARD process: During the backward process the normalising projection contains the correction image to apply to the volume (in the ART sense). The output is the volume itself, the projection image is useless in this case, and the normalising projection is not modified at all.

As for the mode, valid modes are ARTK, CAVK, COUNT_EQ, CAVARTK.

M is the matrix corresponding to the projection process.

Definition at line 1863 of file projection.cpp.

{

    // If projecting forward initialize projections
    if (FORW)
    {
        proj.reset(Ydim, Xdim);
        proj.setAngles(rot, tilt, psi);
        norm_proj().initZeros(proj());
    }

    if( threads > 1 )
    {
        for( int c = 0 ; c < threads ; c++ )
        {
            project_threads[c].global_proj = &proj;
            project_threads[c].global_norm_proj = &norm_proj;
            project_threads[c].FORW = FORW;
            project_threads[c].eq_mode = eq_mode;
            project_threads[c].basis = &basis;
            project_threads[c].M = M;
            project_threads[c].rot = rot;
            project_threads[c].tilt = tilt;
            project_threads[c].psi = psi;
            project_threads[c].ray_length = ray_length;
        }
    }


#ifdef DEBUG_LITTLE
    if (FORW)
    {
        std::cout << "Number of volumes: " << vol.VolumesNo() << std::endl
        << "YdimxXdim: " << Ydim << "x" << Xdim << std::endl;
        for (int i = 0; i < vol.VolumesNo(); i++)
            std::cout << "Volume " << i << std::endl << vol.grid(i) << std::endl;
    }

#endif

    // Project each subvolume
    for (size_t i = 0; i < vol.VolumesNo(); i++)
    {
        Image<int> *VNeq;
        if (GVNeq != NULL)
            VNeq = &((*GVNeq)(i));
        else
            VNeq = NULL;

        if( threads > 1 )
        {
            for( int c = 0 ; c < threads ; c++ )
            {
                project_threads[c].vol = &(vol(i));
                project_threads[c].grid = &(vol.grid(i));
                project_threads[c].VNeq = VNeq;
                project_threads[c].mask = mask;
            }

            barrier_wait( &project_barrier );

            // Here is being processed the volume by the threads

            barrier_wait( &project_barrier );
        }
        else
        {
            // create no thread to do the job
            project_SimpleGrid(&(vol(i)), &(vol.grid(i)), &basis,
                               &proj, &norm_proj, FORW, eq_mode,
                               VNeq, M, mask, ray_length);
        }

#ifdef DEBUG
        Image<double> save;
        save = norm_proj;
        if (FORW)
            save.write((std::string)"PPPnorm_FORW" + (char)(48 + i));
        else
            save.write((std::string)"PPPnorm_BACK" + (char)(48 + i));
#endif

    }
}

project_SimpleGrid()

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
	)

Projection of a Simple Grid. Valid eq_modes are ARTK, CAVARTK and CAV.

Definition at line 1314 of file projection.cpp.

{
    Matrix1D<double> zero(3);                // Origin (0,0,0)
    Matrix1D<double> prjPix(3);       // Position of the pixel within the projection
    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> prjDir(3);       // Projection direction

    Matrix1D<double> actprj(3);       // Coord: Actual position inside
    // the projection plane
    Matrix1D<double> beginZ(3);       // Coord: Plane coordinates of the
    // projection of the 3D point
    // (z0,YY(lowest),XX(lowest))
    Matrix1D<double> univ_beginY(3);  // Coord: coordinates of the
    // grid point
    // (z0,y0,XX(lowest))
    Matrix1D<double> univ_beginZ(3);  // Coord: coordinates of the
    // grid point
    // (z0,YY(lowest),XX(lowest))
    Matrix1D<double> beginY(3);       // Coord: Plane coordinates of the
    // projection of the 3D point
    // (z0,y0,XX(lowest))
    double XX_footprint_size=0.;                // The footprint is supposed
    double YY_footprint_size=0.;                // to be defined between
    double ZZ_footprint_size=0.;
    // (-vmax,+vmax) in the Y axis,
    // and (-umax,+umax) in the X axis
    // This footprint size is the
    // R2 vector (umax,vmax).

    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
    int           foot_V1=0;
    int           foot_U1=0;      // Img Coord: coordinate (in
    // an image fashion, not in an
    // oversampled image fashion)
    // inside the blobprint of the
    // corner1
    int        foot_V;
    int        foot_U;            // Img Coord: coordinate
    int        foot_W = 0;
    // corresponding to the blobprint
    // point which matches with this
    // pixel position
    int           Vsampling=0;
    int           Usampling=0;  // Sampling rate in Y and X
    // directions respectively
    // inside the blobprint
    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
    bool   isVolPSF = false;    // Blob footprint is VolumePSF

    // Prepare system matrix for printing ...................................
    if (M != NULL)
        M->initZeros(YSIZE((*proj)())*XSIZE((*proj)()), grid->get_number_of_samples());

    // 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.
    // actprj is reused with a different purpose
    VECTOR_R3(actprj, 1, 0, 0);
    grid->Gdir_project_to_plane(actprj, proj->euler, prjX);
    VECTOR_R3(actprj, 0, 1, 0);
    grid->Gdir_project_to_plane(actprj, proj->euler, prjY);
    VECTOR_R3(actprj, 0, 0, 1);
    grid->Gdir_project_to_plane(actprj, proj->euler, prjZ);

    // This time the origin of the grid is in the universal coord system
    Uproject_to_plane(grid->origin, proj->euler, prjOrigin);

    // Get the projection direction .........................................
    proj->euler.getRow(2, prjDir);

    // Footprint size .......................................................
    // The following vectors are integer valued vectors, but they are
    // stored as real ones to make easier operations with other vectors
    if (basis->type == Basis::blobs)
    {
        XX_footprint_size = basis->blobprint.umax();
        YY_footprint_size = basis->blobprint.vmax();
        Usampling         = basis->blobprint.Ustep();
        Vsampling         = basis->blobprint.Vstep();

        // Set the limit for grid points out of PSF
        if (basis->VolPSF != NULL)
        {
            isVolPSF = true;
            ZZ_footprint_size = basis->blobprint.wmax();
        }
    }
    else if (basis->type == Basis::voxels || basis->type == Basis::splines)
    {
        YY_footprint_size = XX_footprint_size = CEIL(basis->maxLength());
        Usampling = Vsampling = 0;
    }
    XX_footprint_size += XMIPP_EQUAL_ACCURACY;
    YY_footprint_size += XMIPP_EQUAL_ACCURACY;

    // 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

    // This type conversion gives more speed

    auto ZZ_lowest = (int) ZZ(grid->lowest);

    if( thread_id != -1 )
        ZZ_lowest += thread_id;

    auto YY_lowest = (int) YY(grid->lowest);
    auto XX_lowest = (int) XX(grid->lowest);
    auto ZZ_highest = (int) ZZ(grid->highest);
    auto YY_highest = (int) YY(grid->highest);
    auto XX_highest = (int) XX(grid->highest);

    beginZ = (double)XX_lowest * prjX + (double)YY_lowest * prjY + (double)ZZ_lowest * prjZ + prjOrigin;

    // Check if in VSSNR
    bool VSSNR_mode = (ray_length == basis->maxLength());

#ifdef DEBUG_LITTLE

    int condition;
    condition = thread_id==1;
    if (condition || numthreads==1)
    {
        std::cout << "Equation mode " << eq_mode << std::endl;
        std::cout << "Footprint size " << YY_footprint_size << "x"
        << XX_footprint_size << std::endl;
        std::cout << "rot=" << proj->rot() << " tilt=" << proj->tilt()
        << " psi=" << proj->psi() << std::endl;
        std::cout << "Euler matrix " << proj->euler;
        std::cout << "Projection direction " << prjDir << std::endl;
        std::cout << *grid;
        std::cout << "image limits (" << x0 << "," << y0 << ") (" << xF << ","
        << yF << ")\n";
        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 << "beginZ(coord)  " << beginZ.transpose()    << std::endl;
        std::cout << "lowest         " << XX_lowest  << " " << YY_lowest
        << " " << XX_lowest  << std::endl;
        std::cout << "highest        " << XX_highest << " " << YY_highest
        << " " << XX_highest << std::endl;
        std::cout << "Stats of input basis volume ";
        (*vol)().printStats();
        std::cout << std::endl;
        std::cout.flush();
    }
#endif

    // Compute the grid lattice vectors in space ............................
    Matrix2D<double> grid_basis(3, 3);
    grid_basis = grid->basis * grid->relative_size;
    Matrix1D<double> gridX(3);  // Direction of the grid lattice vectors
    Matrix1D<double> gridY(3);  // in universal coordinates
    Matrix1D<double> gridZ(3);
    Matrix1D<double> univ_position(3);

    grid_basis.getCol(0, gridX);
    grid_basis.getCol(1, gridY);
    grid_basis.getCol(2, gridZ);

    univ_beginZ = (double)XX_lowest * gridX + (double)YY_lowest * gridY + (double)ZZ_lowest * gridZ + grid->origin;

    int number_of_basis = 0;

    for (k = ZZ_lowest; k <= ZZ_highest; k += numthreads)
    {
        // Corner of the row defined by Y
        beginY = beginZ;
        univ_beginY = univ_beginZ;
        for (i = YY_lowest; i <= YY_highest; i++)
        {
            // First point in the row
            actprj = beginY;
            univ_position = univ_beginY;
            for (j = XX_lowest; j <= XX_highest; j++)
            {
                // Ray length interesting
                bool ray_length_interesting = true;
                double zCenterDist;
                double z = 0; // z = 0 standard value for non 3D blobprints

                if (ray_length != -1 || isVolPSF)
                {
                    zCenterDist = point_plane_distance_3D(univ_position, zero,
                                                          proj->direction);
                    if (ray_length != -1)
                        ray_length_interesting = (ABS(zCenterDist) <= ray_length);

                    // Points out of 3DPSF
                    if (isVolPSF && ray_length_interesting)
                    {
                        ray_length_interesting = (ABS(zCenterDist) <= ZZ_footprint_size);
                        // There is still missing the shift of the volume from focal plane
                        z = zCenterDist; // + shiftZ
                    }
                }

                if (grid->is_interesting(univ_position) &&
                    ray_length_interesting)
                {
                    // Be careful that you cannot skip any basis, although its
                    // value be 0, because it is useful for norm_proj
#ifdef DEBUG
                    condition = thread_id == 1;
                    if (condition)
                    {
                        printf("\nProjecting grid coord (%d,%d,%d) ", j, i, k);
                        std::cout << "Vol there = " << VOLVOXEL(*vol, k, i, j) << std::endl;
                        printf(" 3D universal position (%f,%f,%f) \n",
                               XX(univ_position), YY(univ_position), ZZ(univ_position));
                        std::cout << " Center of the basis proj (2D) " << XX(actprj) << "," << YY(actprj) << std::endl;
                        Matrix1D<double> aux;
                        Uproject_to_plane(univ_position, proj->euler, aux);
                        std::cout << " Center of the basis proj (more accurate) " << aux.transpose() << std::endl;
                    }
#endif

                    // Search for integer corners for this basis
                    XX_corner1 = CEIL(XMIPP_MAX(STARTINGX(IMGMATRIX(*proj)), XX(actprj) - XX_footprint_size));
                    YY_corner1 = CEIL(XMIPP_MAX(STARTINGY(IMGMATRIX(*proj)), YY(actprj) - YY_footprint_size));
                    XX_corner2 = FLOOR(XMIPP_MIN(FINISHINGX(IMGMATRIX(*proj)), XX(actprj) + XX_footprint_size));
                    YY_corner2 = FLOOR(XMIPP_MIN(FINISHINGY(IMGMATRIX(*proj)), YY(actprj) + YY_footprint_size));

#ifdef DEBUG

                    if (condition)
                    {
                        std::cout << "Clipped and rounded Corner 1 " << XX_corner1
                        << " " << YY_corner1 << " " << std::endl;
                        std::cout << "Clipped and rounded Corner 2 " << XX_corner2
                        << " " << YY_corner2 << " " << std::endl;
                    }
#endif

                    // Check if the basis falls outside the projection plane
                    // COSS: I have changed here
                    if (XX_corner1 <= XX_corner2 && YY_corner1 <= YY_corner2)
                    {
                        // Compute the index of the basis for corner1
                        if (basis->type == Basis::blobs)
                        {
                            OVER2IMG(basis->blobprint, (double)YY_corner1 - YY(actprj),
                                     (double)XX_corner1 - XX(actprj), foot_V1, foot_U1);
                            if (isVolPSF != false)
                                OVER2IMG_Z(basis->blobprint, z, foot_W);
                        }

                        if (!FORW)
                            vol_corr = 0;

                        // Effectively project this basis
                        // N_eq=(YY_corner2-YY_corner1+1)*(XX_corner2-XX_corner1+1);
                        N_eq = 0;
                        foot_V = foot_V1;
                        for (int y = YY_corner1; y <= YY_corner2; y++)
                        {
                            foot_U = foot_U1;
                            for (int x = XX_corner1; x <= XX_corner2; x++)
                            {
                                if (!((mask != NULL) && A2D_ELEM(*mask,y,x)<0.5))
                                {
#ifdef DEBUG
                                    if (condition)
                                    {
                                        std::cout << "Position in projection (" << x << ","
                                        << y << ") ";
                                        double y, x;
                                        if (basis->type == Basis::blobs)
                                        {
                                            std::cout << "in footprint ("
                                            << foot_U << "," << foot_V << ")";
                                            IMG2OVER(basis->blobprint, foot_V, foot_U, y, x);
                                            std::cout << " (d= " << sqrt(y*y + x*x) << ") ";
                                            fflush(stdout);
                                        }
                                    }
#endif
                                    double a;
                                    double a2;
                                    // Check if volumetric interpolation (i.e., SSNR)
                                    if (VSSNR_mode)
                                    {
                                        // This is the VSSNR case
                                        // Get the pixel position in the universal coordinate
                                        // system
                                        SPEED_UP_temps012;
                                        VECTOR_R3(prjPix, x, y, z);
                                        M3x3_BY_V3x1(prjPix, proj->eulert, prjPix);
                                        V3_MINUS_V3(prjPix, prjPix, univ_position);
                                        a = basis->valueAt(prjPix);
                                        a2 = a * a;
                                    }
                                    else
                                    {
                                        // This is normal reconstruction from projections
                                        if (basis->type == Basis::blobs)
                                        {
                                            // Projection of a blob
                                            a = VOLVOXEL(basis->blobprint, foot_W, foot_V, foot_U);
                                            a2 = VOLVOXEL(basis->blobprint2, foot_W, foot_V, foot_U);

                                        }
                                        else
                                        {
                                            // Projection of other bases
                                            // If the basis is big enough, then
                                            // it is not necessary to integrate at several
                                            // places. Big enough is being greater than
                                            // 1.41 which is the maximum separation
                                            // between two pixels
                                            if (XX_footprint_size > 1.41)
                                            {
                                                // Get the pixel in universal coordinates
                                                SPEED_UP_temps012;
                                                VECTOR_R3(prjPix, x, y, 0);
                                                // Express the point in a local coordinate system
                                                M3x3_BY_V3x1(prjPix, proj->eulert, prjPix);
#ifdef DEBUG

                                                if (condition)
                                                    std::cout << " in volume coord ("
                                                    << prjPix.transpose() << ")";
#endif

                                                V3_MINUS_V3(prjPix, prjPix, univ_position);
#ifdef DEBUG

                                                if (condition)
                                                    std::cout << " in voxel coord ("
                                                    << prjPix.transpose() << ")";
#endif

                                                a = basis->projectionAt(prjDir, prjPix);
                                                a2 = a * a;
                                            }
                                            else
                                            {
                                                // If the basis is too small (of the
                                                // range of the voxel), then it is
                                                // necessary to sample in a few places
                                                const double p0 = 1.0 / (2 * ART_PIXEL_SUBSAMPLING) - 0.5;
                                                const double pStep = 1.0 / ART_PIXEL_SUBSAMPLING;
                                                const double pAvg = 1.0 / (ART_PIXEL_SUBSAMPLING * ART_PIXEL_SUBSAMPLING);
                                                int ii;
                                                int jj;
                                                double px;
                                                double py;
                                                a = 0;
#ifdef DEBUG

                                                if (condition)
                                                    std::cout << std::endl;
#endif

                                                for (ii = 0, px = p0; ii < ART_PIXEL_SUBSAMPLING; ii++, px += pStep)
                                                    for (jj = 0, py = p0; jj < ART_PIXEL_SUBSAMPLING; jj++, py += pStep)
                                                    {
#ifdef DEBUG
                                                        if (condition)
                                                            std::cout << "    subsampling (" << ii << ","
                                                            << jj << ") ";
#endif

                                                        SPEED_UP_temps012;
                                                        // Get the pixel in universal coordinates
                                                        VECTOR_R3(prjPix, x + px, y + py, 0);
                                                        // Express the point in a local coordinate system
                                                        M3x3_BY_V3x1(prjPix, proj->eulert, prjPix);
#ifdef DEBUG

                                                        if (condition)
                                                            std::cout << " in volume coord ("
                                                            << prjPix.transpose() << ")";
#endif

                                                        V3_MINUS_V3(prjPix, prjPix, univ_position);
#ifdef DEBUG

                                                        if (condition)
                                                            std::cout << " in voxel coord ("
                                                            << prjPix.transpose() << ")";
#endif

                                                        a += basis->projectionAt(prjDir, prjPix);
#ifdef DEBUG

                                                        if (condition)
                                                            std::cout << " partial a="
                                                            << basis->projectionAt(prjDir, prjPix)
                                                            << std::endl;
#endif

                                                    }
                                                a *= pAvg;
                                                a2 = a * a;
#ifdef DEBUG

                                                if (condition)
                                                    std::cout << "   Finally ";
#endif

                                            }
                                        }
                                    }
#ifdef DEBUG
                                    if (condition)
                                        std::cout << "=" << a << " , " << a2;
#endif

                                    if (FORW)
                                    {
                                        switch (eq_mode)
                                        {
                                        case CAVARTK:
                                        case ARTK:
                                            IMGPIXEL(*proj, y, x) += VOLVOXEL(*vol, k, i, j) * a;
                                            IMGPIXEL(*norm_proj, y, x) += a2;
                                            if (M != NULL)
                                            {
                                                int py;
                                                int px;
                                                (*proj)().toPhysical(y, x, py, px);
                                                int number_of_pixel = py * XSIZE((*proj)()) + px;
                                                dMij(*M, number_of_pixel, number_of_basis) = a;
                                            }
                                            break;
                                        case CAVK:
                                            IMGPIXEL(*proj, y, x) += VOLVOXEL(*vol, k, i, j) * a;
                                            IMGPIXEL(*norm_proj, y, x) += a2 * N_eq;
                                            break;
                                        case COUNT_EQ:
                                            VOLVOXEL(*vol, k, i, j)++;
                                            break;
                                        case CAV:
                                            IMGPIXEL(*proj, y, x) += VOLVOXEL(*vol, k, i, j) * a;
                                            IMGPIXEL(*norm_proj, y, x) += a2 *
                                                                          VOLVOXEL(*VNeq, k, i, j);
                                            break;
                                        }

#ifdef DEBUG
                                        if (condition)
                                        {
                                            std::cout << " proj= " << IMGPIXEL(*proj, y, x)
                                            << " norm_proj=" << IMGPIXEL(*norm_proj, y, x) << std::endl;
                                            std::cout.flush();
                                        }
#endif

                                    }
                                    else
                                    {
                                        vol_corr += IMGPIXEL(*norm_proj, y, x) * a;
                                        if (a != 0)
                                            N_eq++;
#ifdef DEBUG

                                        if (condition)
                                        {
                                            std::cout << " corr_img= " << IMGPIXEL(*norm_proj, y, x)
                                            << " correction=" << vol_corr << std::endl;
                                            std::cout.flush();
                                        }
#endif

                                    }
                                }
                                // Prepare for next operation
                                foot_U += Usampling;
                            }
                            foot_V += Vsampling;
                        } // Project this basis

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

#ifdef DEBUG

                            if (condition)
                            {
                                printf("\nFinal value at (%d,%d,%d) ", j, i, k);
                                std::cout << " = " << VOLVOXEL(*vol, k, i, j) << std::endl;
                                std::cout.flush();
                            }
#endif

                        }
                    } // If not collapsed
                    number_of_basis++;
                } // If interesting

                // Prepare for next iteration
                V2_PLUS_V2(actprj, actprj, prjX);
                V3_PLUS_V3(univ_position, univ_position, gridX);
            }
            V2_PLUS_V2(beginY, beginY, prjY);
            V3_PLUS_V3(univ_beginY, univ_beginY, gridY);
        }
        V2_PLUS_V2(beginZ, beginZ, prjZ * numthreads );
        V3_PLUS_V3(univ_beginZ, univ_beginZ, gridZ * numthreads);
    }
}

project_SimpleGridThread()

template<class T >

void* project_SimpleGridThread

(

void *

params

)

Threaded projection for simple grids

Definition at line 1215 of file projection.cpp.

{
    auto * thread_data = (project_thread_params *)params;

    Image<T> * vol;
    const SimpleGrid * grid;

    int thread_id = thread_data->thread_id;
    int threads_count = thread_data->threads_count;

    const Basis * basis;

    Projection * proj;
    Projection * norm_proj;

    auto * forw_proj = new Projection();
    auto * forw_norm_proj = new Projection();

    Projection * global_proj;
    Projection * global_norm_proj;

    int FORW;
    int eq_mode;
    const Image<int> *VNeq=NULL;
    Matrix2D<double> *M=NULL;
    const MultidimArray<int> *mask=NULL;
    double ray_length;

    double rot;
    double tilt;
    double psi;

    do
    {
        barrier_wait( &project_barrier );

        if( thread_data->destroy == true )
            break;

        vol = thread_data->vol;
        grid = thread_data->grid;
        basis = thread_data->basis;
        FORW = thread_data->FORW;
        eq_mode = thread_data->eq_mode;
        VNeq = thread_data->VNeq;
        M = thread_data->M;
        mask = thread_data->mask;
        ray_length = thread_data->ray_length;
        global_proj = thread_data->global_proj;
        global_norm_proj = thread_data->global_norm_proj;

        rot = thread_data->rot;
        tilt = thread_data->tilt;
        psi = thread_data->psi;

        if( FORW )
        {
            proj = forw_proj;
            norm_proj = forw_norm_proj;

            proj->reset(YSIZE( (*global_proj)() ), XSIZE( (*global_proj)() ));
            proj->setAngles(rot, tilt, psi);
            (*norm_proj)().initZeros((*proj)());
        }
        else
        {
            proj = global_proj;
            norm_proj = global_norm_proj;
        }

        project_SimpleGrid( vol, grid,
                            basis,
                            proj,
                            norm_proj, FORW, eq_mode,
                            VNeq, M, mask,
                            ray_length,
                            thread_id,
                            threads_count
                          );

        if( FORW )
        {

            pthread_mutex_lock( &project_mutex );

            (*global_proj)() = (*global_proj)() + (*proj)();
            (*global_norm_proj)() = (*global_norm_proj)() + (*norm_proj)();

            pthread_mutex_unlock( &project_mutex );
        }

        barrier_wait( &project_barrier );
    }
    while(1);

    return (void *)NULL;
}

projectVolume()

void projectVolume	(	MultidimArray< double > &	V,
		Projection &	P,
		int	Ydim,
		int	Xdim,
		double	rot,
		double	tilt,
		double	psi,
		const Matrix1D< double > *	roffset = nullptr
	)

From voxel volumes. The voxel volume is projected onto a projection plane defined by (rot, tilt, psi) (1st, 2nd and 3rd Euler angles) . The projection is previously is resized to Ydim x Xdim and initialized to 0. The projection itself, from now on, will keep the Euler angles.

The offset is a 3D vector specifying the offset that must be applied when going from the projection space to the universal space

rproj=E*r+roffset => r=E^t (rproj-roffset)

Set it to NULL if you don't want to use it

Definition at line 337 of file projection.cpp.

{
    SPEED_UP_temps012;

    // Initialise projection
    P.reset(Ydim, Xdim);
    P.setAngles(rot, tilt, psi);

    // Compute the distance for this line crossing one voxel
    int x_0 = STARTINGX(V);
    int x_F = FINISHINGX(V);
    int y_0 = STARTINGY(V);
    int y_F = FINISHINGY(V);
    int z_0 = STARTINGZ(V);
    int z_F = FINISHINGZ(V);

    // Distances in X and Y between the center of the projection pixel begin
    // computed and each computed ray
    double step = 1.0 / 3.0;

    // Avoids divisions by zero and allows orthogonal rays computation
    if (XX(P.direction) == 0)
        XX(P.direction) = XMIPP_EQUAL_ACCURACY;
    if (YY(P.direction) == 0)
        YY(P.direction) = XMIPP_EQUAL_ACCURACY;
    if (ZZ(P.direction) == 0)
        ZZ(P.direction) = XMIPP_EQUAL_ACCURACY;

    // Some precalculated variables
    int x_sign = SGN(XX(P.direction));
    int y_sign = SGN(YY(P.direction));
    int z_sign = SGN(ZZ(P.direction));
    double half_x_sign = 0.5 * x_sign;
    double half_y_sign = 0.5 * y_sign;
    double half_z_sign = 0.5 * z_sign;
    double iXXP_direction=1.0/XX(P.direction);
    double iYYP_direction=1.0/YY(P.direction);
    double iZZP_direction=1.0/ZZ(P.direction);

    MultidimArray<double> &mP = P();
    Matrix1D<double>  v(3);
    Matrix1D<double> r_p(3); // r_p are the coordinates of the
    // pixel being projected in the
    // coordinate system attached to the
    // projection
    Matrix1D<double> p1(3);  // coordinates of the pixel in the
    // universal space
    Matrix1D<double> p1_shifted(3); // shifted half a pixel
    FOR_ALL_ELEMENTS_IN_ARRAY2D(mP)
    {
        double ray_sum = 0.0;    // Line integral value

        // Computes 4 different rays for each pixel.
        for (int rays_per_pixel = 0; rays_per_pixel < 4; rays_per_pixel++)
        {
            // universal coordinate system
            switch (rays_per_pixel)
            {
            case 0:
                VECTOR_R3(r_p, j - step, i - step, 0);
                break;
            case 1:
                VECTOR_R3(r_p, j - step, i + step, 0);
                break;
            case 2:
                VECTOR_R3(r_p, j + step, i - step, 0);
                break;
            case 3:
                VECTOR_R3(r_p, j + step, i + step, 0);
                break;
            }

            // Express r_p in the universal coordinate system
            if (roffset!=nullptr)
                r_p-=*roffset;
            M3x3_BY_V3x1(p1, P.eulert, r_p);
            XX(p1_shifted)=XX(p1)-half_x_sign;
            YY(p1_shifted)=YY(p1)-half_y_sign;
            ZZ(p1_shifted)=ZZ(p1)-half_z_sign;

            // Compute the minimum and maximum alpha for the ray
            // intersecting the given volume
            double alpha_xmin = (x_0 - 0.5 - XX(p1))* iXXP_direction;
            double alpha_xmax = (x_F + 0.5 - XX(p1))* iXXP_direction;
            double alpha_ymin = (y_0 - 0.5 - YY(p1))* iYYP_direction;
            double alpha_ymax = (y_F + 0.5 - YY(p1))* iYYP_direction;
            double alpha_zmin = (z_0 - 0.5 - ZZ(p1))* iZZP_direction;
            double alpha_zmax = (z_F + 0.5 - ZZ(p1))* iZZP_direction;

            double auxMin;
            double auxMax;
            if (alpha_xmin<alpha_xmax)
            {
                auxMin=alpha_xmin;
                auxMax=alpha_xmax;
            }
            else
            {
                auxMin=alpha_xmax;
                auxMax=alpha_xmin;
            }
            double alpha_min=auxMin;
            double alpha_max=auxMax;
            if (alpha_ymin<alpha_ymax)
            {
                auxMin=alpha_ymin;
                auxMax=alpha_ymax;
            }
            else
            {
                auxMin=alpha_ymax;
                auxMax=alpha_ymin;
            }
            alpha_min=fmax(auxMin,alpha_min);
            alpha_max=fmin(auxMax,alpha_max);
            if (alpha_zmin<alpha_zmax)
            {
                auxMin=alpha_zmin;
                auxMax=alpha_zmax;
            }
            else
            {
                auxMin=alpha_zmax;
                auxMax=alpha_zmin;
            }
            alpha_min=fmax(auxMin,alpha_min);
            alpha_max=fmin(auxMax,alpha_max);
            if (alpha_max - alpha_min < XMIPP_EQUAL_ACCURACY)
                continue;

#ifdef DEBUG

            std::cout << "Pixel:  " << r_p.transpose() << std::endl
            << "Univ:   " << p1.transpose() << std::endl
            << "Dir:    " << P.direction.transpose() << std::endl
            << "Alpha x:" << alpha_xmin << " " << alpha_xmax << std::endl
            << "   " << (p1 + alpha_xmin*P.direction).transpose() << std::endl
            << "   " << (p1 + alpha_xmax*P.direction).transpose() << std::endl
            << "Alpha y:" << alpha_ymin << " " << alpha_ymax << std::endl
            << "   " << (p1 + alpha_ymin*P.direction).transpose() << std::endl
            << "   " << (p1 + alpha_ymax*P.direction).transpose() << std::endl
            << "Alpha z:" << alpha_zmin << " " << alpha_zmax << std::endl
            << "   " << (p1 + alpha_zmin*P.direction).transpose() << std::endl
            << "   " << (p1 + alpha_zmax*P.direction).transpose() << std::endl
            << "alpha  :" << alpha_min  << " " << alpha_max  << std::endl
            << std::endl;
#endif

            // Compute the first point in the volume intersecting the ray
            double zz_idxd;
            double yy_idxd;
            double xx_idxd;
            int    zz_idx;
            int    yy_idx;
            int    xx_idx;
            V3_BY_CT(v, P.direction, alpha_min);
            V3_PLUS_V3(v, p1, v);

            // Compute the index of the first voxel
            xx_idx = ROUND(XX(v));
            yy_idx = ROUND(YY(v));
            zz_idx = ROUND(ZZ(v));

            xx_idxd = xx_idx = CLIP(xx_idx, x_0, x_F);
            yy_idxd = yy_idx = CLIP(yy_idx, y_0, y_F);
            zz_idxd = zz_idx = CLIP(zz_idx, z_0, z_F);

#ifdef DEBUG

            std::cout << "First voxel: " << v.transpose() << std::endl;
            std::cout << "   First index: " << idx.transpose() << std::endl;
            std::cout << "   Alpha_min: " << alpha_min << std::endl;
#endif

            // Follow the ray
            double alpha = alpha_min;
            do
            {
#ifdef DEBUG
                std::cout << " \n\nCurrent Value: " << V(zz_idx, yy_idx, xx_idx) << std::endl;
#endif

                double alpha_x = (xx_idxd - XX(p1_shifted))* iXXP_direction;
                double alpha_y = (yy_idxd - YY(p1_shifted))* iYYP_direction;
                double alpha_z = (zz_idxd - ZZ(p1_shifted))* iZZP_direction;

                // Which dimension will ray move next step into?, it isn't necessary to be only
                // one.
                double diffx = fabs(alpha-alpha_x);
                double diffy = fabs(alpha-alpha_y);
                double diffz = fabs(alpha-alpha_z);
                int diff_source=0;
                double diff_alpha=diffx;
                if (diffy<diff_alpha)
                {
                    diff_source=1;
                    diff_alpha=diffy;
                }
                if (diffz<diff_alpha)
                {
                    diff_source=2;
                    diff_alpha=diffz;
                }
                ray_sum += diff_alpha * A3D_ELEM(V, zz_idx, yy_idx, xx_idx);

                switch (diff_source)
                {
                case 0:
                    alpha = alpha_x;
                    xx_idx += x_sign;
                    xx_idxd = xx_idx;
                    break;
                case 1:
                    alpha = alpha_y;
                    yy_idx += y_sign;
                    yy_idxd = yy_idx;
                    break;
                default:
                    alpha = alpha_z;
                    zz_idx += z_sign;
                    zz_idxd = zz_idx;
                }

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

                XX(v) += diff_alpha * XX(P.direction);
                YY(v) += diff_alpha * YY(P.direction);
                ZZ(v) += diff_alpha * ZZ(P.direction);

                std::cout << "    Next entry point: " << v.transpose() << std::endl
                << "    Index: " << idx.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);
        } // for

        A2D_ELEM(IMGMATRIX(P), i, j) = ray_sum * 0.25;
#ifdef DEBUG

        std::cout << "Assigning P(" << i << "," << j << ")=" << ray_sum << std::endl;
#endif

    }
}

projectVolumeOffCentered()

void projectVolumeOffCentered	(	MultidimArray< double > &	V,
		Projection &	P,
		int	Ydim,
		int	Xdim
	)

From voxel volumes, off-centered tilt axis. This routine projects a volume that is rotating (angle) degrees around the axis defined by the two angles (axisRot,axisTilt) and that passes through the point raxis. The projection can be further inplane rotated and shifted through the parameters (inplaneRot) and (rinplane).

All vectors involved must be 3D.

The projection model is rproj=H Rinplane Raxis r + Rinplane (I-Raxis) raxis + rinplane

Where Raxis is the 3D rotation matrix given by the axis and the angle.

Definition at line 594 of file projection.cpp.

{
    Matrix1D<double> roffset(3);
    P.getShifts(XX(roffset), YY(roffset), ZZ(roffset));

    projectVolume(V, P, Ydim, Xdim, P.rot(), P.tilt(), P.psi(), &roffset);
}

singleWBP()

void singleWBP	(	MultidimArray< double > &	V,
		Projection &	P
	)

Single Weighted Back Projection. Projects a single particle into a voxels volume by updating its components this way: Voxel(i,j,k) = Voxel(i,j,k) + Pixel( x,y) * Distance.

Where:

• Voxel(i,j,k) is the voxel the ray is crossing.
• Pixel( y,z ) is the value of the pixel where the ray departs from.
• Distance is the distance the ray goes through the Voxel.

Definition at line 609 of file projection.cpp.

{
    SPEED_UP_temps012;

    // Compute the distance for this line crossing one voxel
    int x_0 = STARTINGX(V);
    int x_F = FINISHINGX(V);
    int y_0 = STARTINGY(V);
    int y_F = FINISHINGY(V);
    int z_0 = STARTINGZ(V);
    int z_F = FINISHINGZ(V);

    // Avoids divisions by zero and allows orthogonal rays computation
    if (XX(P.direction) == 0)
        XX(P.direction) = XMIPP_EQUAL_ACCURACY;
    if (YY(P.direction) == 0)
        YY(P.direction) = XMIPP_EQUAL_ACCURACY;
    if (ZZ(P.direction) == 0)
        ZZ(P.direction) = XMIPP_EQUAL_ACCURACY;

    // Some precalculated variables
    int x_sign = SGN(XX(P.direction));
    int y_sign = SGN(YY(P.direction));
    int z_sign = SGN(ZZ(P.direction));
    double half_x_sign = 0.5 * x_sign;
    double half_y_sign = 0.5 * y_sign;
    double half_z_sign = 0.5 * z_sign;

    MultidimArray<double> &mP = P();
    Matrix1D<double> r_p(3); // r_p are the coordinates of the
    // pixel being projected in the
    // coordinate system attached to the
    // projection
    Matrix1D<double> p1(3);  // coordinates of the pixel in the
    Matrix1D<double>  v(3);
    Matrix1D<int>    idx(3);
    FOR_ALL_ELEMENTS_IN_ARRAY2D(mP)
    {
        // Computes 4 different rays for each pixel.
        VECTOR_R3(r_p, j, i, 0);

        // Express r_p in the universal coordinate system
        M3x3_BY_V3x1(p1, P.eulert, r_p);

        // Compute the minimum and maximum alpha for the ray
        // intersecting the given volume
        double alpha_xmin = (x_0 - 0.5 - XX(p1)) / XX(P.direction);
        double alpha_xmax = (x_F + 0.5 - XX(p1)) / XX(P.direction);
        double alpha_ymin = (y_0 - 0.5 - YY(p1)) / YY(P.direction);
        double alpha_ymax = (y_F + 0.5 - YY(p1)) / YY(P.direction);
        double alpha_zmin = (z_0 - 0.5 - ZZ(p1)) / ZZ(P.direction);
        double alpha_zmax = (z_F + 0.5 - ZZ(p1)) / ZZ(P.direction);

        double alpha_min = fmax(fmin(alpha_xmin, alpha_xmax),
                                fmin(alpha_ymin, alpha_ymax));
        alpha_min = fmax(alpha_min, fmin(alpha_zmin, alpha_zmax));
        double alpha_max = fmin(fmax(alpha_xmin, alpha_xmax),
                                fmax(alpha_ymin, alpha_ymax));
        alpha_max = fmin(alpha_max, fmax(alpha_zmin, alpha_zmax));
        if (alpha_max - alpha_min < XMIPP_EQUAL_ACCURACY)
            continue;

        // Compute the first point in the volume intersecting the ray
        V3_BY_CT(v, P.direction, alpha_min);
        V3_PLUS_V3(v, p1, v);

        // Compute the index of the first voxel
        XX(idx) = CLIP(ROUND(XX(v)), x_0, x_F);
        YY(idx) = CLIP(ROUND(YY(v)), y_0, y_F);
        ZZ(idx) = CLIP(ROUND(ZZ(v)), z_0, z_F);


#ifdef DEBUG

        std::cout << "First voxel: " << v.transpose() << std::endl;
        std::cout << "   First index: " << idx.transpose() << std::endl;
        std::cout << "   Alpha_min: " << alpha_min << std::endl;
#endif

        // Follow the ray
        double alpha = alpha_min;
        do
        {
#ifdef DEBUG
            std::cout << " \n\nCurrent Value: " << V(ZZ(idx), YY(idx), XX(idx)) << std::endl;
#endif

            double alpha_x = (XX(idx) + half_x_sign - XX(p1)) / XX(P.direction);
            double alpha_y = (YY(idx) + half_y_sign - YY(p1)) / YY(P.direction);
            double alpha_z = (ZZ(idx) + half_z_sign - ZZ(p1)) / ZZ(P.direction);

            // 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 = fmin(fmin(diffx, diffy), diffz);

            A3D_ELEM(V, ZZ(idx), YY(idx), XX(idx)) += diff_alpha * A2D_ELEM(P(), i, j);

            if (ABS(diff_alpha - diffx) <= XMIPP_EQUAL_ACCURACY)
            {
                alpha = alpha_x;
                XX(idx) += x_sign;
            }
            if (ABS(diff_alpha - diffy) <= XMIPP_EQUAL_ACCURACY)
            {
                alpha = alpha_y;
                YY(idx) += y_sign;
            }
            if (ABS(diff_alpha - diffz) <= XMIPP_EQUAL_ACCURACY)
            {
                alpha = alpha_z;
                ZZ(idx) += z_sign;
            }

        }
        while ((alpha_max - alpha) > XMIPP_EQUAL_ACCURACY);
#ifdef DEBUG

        std::cout << "Assigning P(" << i << "," << j << ")=" << ray_sum << std::endl;
#endif

    }
}

Variable Documentation

ART_PIXEL_SUBSAMPLING

const int ART_PIXEL_SUBSAMPLING = 2

Definition at line 296 of file projection.h.

ARTK

constexpr int ARTK = 1

Definition at line 177 of file projection.h.

BACKWARD

constexpr int BACKWARD = 0

Definition at line 175 of file projection.h.

CAV

constexpr int CAV = 4

Definition at line 180 of file projection.h.

CAVARTK

constexpr int CAVARTK = 5

Definition at line 181 of file projection.h.

CAVK

constexpr int CAVK = 2

Definition at line 178 of file projection.h.

COUNT_EQ

constexpr int COUNT_EQ = 3

Definition at line 179 of file projection.h.

FORWARD

constexpr int FORWARD = 1

Definition at line 174 of file projection.h.

project_barrier

barrier_t project_barrier

Definition at line 46 of file projection.cpp.

project_mutex

pthread_mutex_t project_mutex

Definition at line 47 of file projection.cpp.

project_threads

project_thread_params* project_threads

Definition at line 48 of file projection.cpp.