project_xr (project for tilt series)

Collaboration diagram for project_xr (project for tilt series):

Classes
struct	XrayThreadArgument
	Data struct to be passed to threads. More...
class	XrayProjPhantom
class	ProgXrayProject
struct	CIGTArgument
struct	PXSGTArgument

Functions
void	XrayRotateAndProjectVolumeOffCentered (XrayProjPhantom &side, XRayPSF &psf, Projection &P, Projection &standardP, int Ydim, int Xdim)
void	projectXrayVolume (MultidimArray< double > &muVol, MultidimArray< double > &IgeoVol, XRayPSF &psf, Projection &P, MultidimArray< double > *projNorm=nullptr, ThreadManager *ThrMgr=nullptr)
void	threadXrayProject (ThreadArgument &thArg)
	Thread Job to generate an X-ray microscope projection. More...
void	calculateIgeo (MultidimArray< double > &muVol, double sampling, MultidimArray< double > &IgeoVol, MultidimArray< double > &cumMu, MultidimArray< double > &IgeoZb, int nThreads=1, ThreadManager *ThrMgr=nullptr)
	Calculate the volume of the information of Igeometric at each plane of the phantom. More...
void	calculateIgeoThread (ThreadArgument &thArg)
void	projectXrayGridVolume (MultidimArray< double > &muVol, XRayPSF &psf, MultidimArray< double > &IgeoVol, Projection &proj, MultidimArray< double > *projNorm, int FORW, ThreadManager *thMgr)
	Project as in an X-ray microscope using a grids and blobs. More...
void	projectXraySimpleGrid (MultidimArray< double > *vol, const XRayPSF &psf, MultidimArray< double > *IgeoVol, Projection *proj, MultidimArray< double > &projNorm, int FORW, int threadId=-1, int nThreads=1)
void	projectXraySimpleGridThread (ThreadArgument &thArg)

Variables
Mutex	mutex
Barrier *	barrier
ThreadManager *	thMgr

Detailed Description

Function Documentation

calculateIgeo()

void calculateIgeo	(	MultidimArray< double > &	muVol,
		double	sampling,
		MultidimArray< double > &	IgeoVol,
		MultidimArray< double > &	cumMu,
		MultidimArray< double > &	IgeoZb,
		int	nThreads = 1,
		ThreadManager *	ThrMgr = nullptr
	)

Calculate the volume of the information of Igeometric at each plane of the phantom.

Definition at line 616 of file project_xray.cpp.

{
    // Checking sizes
    if ( (XSIZE(muVol) != XSIZE(IgeoZb)) || (YSIZE(muVol) != YSIZE(IgeoZb)) )
        REPORT_ERROR(ERR_MULTIDIM_SIZE, "CalculateIgeo: IgeoZb size does not match muVol size.");

    IgeoVol.initZeros(muVol);
    cumMu.initZeros(YSIZE(muVol), XSIZE(muVol));

    ThreadManager * myThMgr;

    if (ThrMgr == nullptr)
        myThMgr = new ThreadManager(nThreads);
    else
        myThMgr = ThrMgr;

    CIGTArgument iGeoArgs;
    iGeoArgs.samplingZ = sampling;
    iGeoArgs.muVol = &muVol;
    iGeoArgs.IgeoVol = &IgeoVol;
    iGeoArgs.cumMu = &cumMu;
    iGeoArgs.IgeoZb = & IgeoZb;
    iGeoArgs.td = new ThreadTaskDistributor(YSIZE(muVol),YSIZE(muVol)/nThreads/2);

    myThMgr->run(calculateIgeoThread,&iGeoArgs);

    delete iGeoArgs.td;

    if (ThrMgr == nullptr)
        delete myThMgr;
}

calculateIgeoThread()

void calculateIgeoThread

(

ThreadArgument &

thArg

)

Definition at line 650 of file project_xray.cpp.

{
    auto *dataThread = (CIGTArgument*) thArg.data;

    double sampling = dataThread->samplingZ;
    MultidimArray<double> &muVol = *(dataThread->muVol);
    MultidimArray<double> &IgeoVol = *(dataThread->IgeoVol);
    MultidimArray<double> &cumMu = *(dataThread->cumMu);
    MultidimArray<double> &IgeoZb = *(dataThread->IgeoZb);
    ParallelTaskDistributor * td = dataThread->td;

    // Resizing
    //    IgeoVol.resize(muVol, false);

    size_t first, last;
    while (td->getTasks(first, last))
    {
        //        first--;
        //        last--;
        for (size_t i = first; i <= last; ++i)
        {
            //            for (int j=0; j<XSIZE(muVol); ++j)
            //                dAkij(IgeoVol,0,i,j) = dAij(IgeoZb,i,j)*exp(-dAkij(muVol,0,i,j)*sampling);
            //
            //            for (int k = 1; k < ZSIZE(muVol); ++k)
            //                for (int j=0; j<XSIZE(muVol); ++j)
            //                    dAkij(IgeoVol,k,i,j) = dAkij(IgeoVol,k-1,i,j)*exp(-dAkij(muVol,k,i,j)*sampling);

            for (size_t k = 0; k < ZSIZE(muVol); ++k)
                for (size_t j=0; j<XSIZE(muVol); ++j)
                {
                    dAij(cumMu,i,j) += dAkij(muVol,k,i,j)*sampling;
                    dAkij(IgeoVol,k,i,j) = dAij(IgeoZb,i,j)*exp(-dAij(cumMu,i,j));
                }
        }
    }
}

projectXrayGridVolume()

void projectXrayGridVolume	(	MultidimArray< double > &	muVol,
		XRayPSF &	psf,
		MultidimArray< double > &	IgeoVol,
		Projection &	proj,
		MultidimArray< double > *	projNorm,
		int	FORW,
		ThreadManager *	thMgr
	)

Project as in an X-ray microscope using a grids and blobs.

Parameters

proj	Vol with the Igeometrical distribution along specimen volume

Definition at line 688 of file project_xray.cpp.

{
    XrayThreadArgument projThrData;
    Barrier barrier(thMgr->threads);

    projThrData.psf = &psf;
    projThrData.muVol = &muVol;
    projThrData.IgeoVol = &IgeoVol;
    projThrData.projOut = &proj;
    projThrData.projNorm = projNorm;
    projThrData.forw = FORW;
    projThrData.phantomSlabIdx = nullptr;
    projThrData.psfSlicesIdx = nullptr;

    //    //Create the job handler to distribute thread jobs
    //    size_t blockSize, numberOfJobs = psfSlicesIdx.size() ;
    //    blockSize = 1;
    //    ThreadTaskDistributor td(numberOfJobs, blockSize);
    //    projThrData.td = &td;
    //    projThrData.barrier = &barrier;

    //the really really final project routine, I swear by Snoopy.
    //    project_xr(psf,side.rotPhantomVol,P, idxSlice);

    thMgr->run(projectXraySimpleGridThread,(void*) &projThrData);

}

projectXraySimpleGrid()

void projectXraySimpleGrid	(	MultidimArray< double > *	vol,
		const XRayPSF &	psf,
		MultidimArray< double > *	IgeoVol,
		Projection *	proj,
		MultidimArray< double > &	projNorm,
		int	FORW,
		int	threadId = -1,
		int	nThreads = 1
	)

ARTK Equation mode

Definition at line 747 of file project_xray.cpp.

{

    MultidimArray<double> psfVol;
    psfVol.alias(psf.psfVol());


    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;                // The footprint is supposed
    double YY_footprint_size;                // 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).

    int XX_corner2, XX_corner1;              // Coord: Corners of the
    int YY_corner2, YY_corner1;              // footprint when it is projected
    // onto the projection plane
    // 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, j, k;                   // volume element indexes

    // 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);
    //    Uproject_to_plane(actprj, proj->euler, prjX);
    //    VECTOR_R3(actprj, 0, 1, 0);
    //    Uproject_to_plane(actprj, proj->euler, prjY);
    //    VECTOR_R3(actprj, 0, 0, 1);
    //    Uproject_to_plane(actprj, proj->euler, prjZ);
    proj->euler.getCol(0, prjX);
    proj->euler.getCol(1, prjY);
    proj->euler.getCol(2, prjZ);

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

    // 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
    XX_footprint_size = FINISHINGX(psfVol) + XMIPP_EQUAL_ACCURACY;
    YY_footprint_size = FINISHINGY(psfVol) + 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

    int ZZ_lowest = STARTINGZ(*vol);
    if( threadId != -1 )
        ZZ_lowest += threadId;
    int YY_lowest = STARTINGY(*vol);
    int XX_lowest = STARTINGX(*vol);
    int ZZ_highest = FINISHINGZ(*vol);
    int YY_highest = FINISHINGY(*vol);
    int XX_highest = FINISHINGX(*vol);

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


#ifdef DEBUG_LITTLE

    int condition;
    condition = threadId==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

    for (k = ZZ_lowest; k <= ZZ_highest; k += numthreads)
    {
        // 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++)
            {
                // Ray length interesting
                bool ray_length_interesting = true;

                //                // Points out of 3DPSF
                //                ray_length_interesting = (ABS(zCenterDist) <= ZZ_footprint_size);
                //                // There is still missing the shift of the volume from focal plane

                if (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 = threadId == 1;
                    if (condition)
                    {
                        printf("\nProjecting grid coord (%d,%d,%d) ", j, i, k);
                        std::cout << "Vol there = " << A3D_ELEM(*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 voxel falls outside the projection plane
                    // COSS: I have changed here
                    if (XX_corner1 <= XX_corner2 && YY_corner1 <= YY_corner2
                        && INSIDEZ(*IgeoVol,ZZ(actprj)))
                    {

                        double IgeoValue = A3D_ELEM(*IgeoVol, (int)ZZ(actprj), (int)YY(actprj), (int)XX(actprj));

                        if (!FORW)
                            vol_corr = 0;

                        N_eq = 0;
                        for (int y = YY_corner1; y <= YY_corner2; y++)
                        {
                            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, a2;

                                    a = A3D_ELEM(psfVol, (int)ZZ(actprj), y, x) * IgeoValue;
                                    a2 = a * a;

#ifdef DEBUG

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

                                    if (FORW)
                                    {
                                        IMGPIXEL(*proj, y, x) += A3D_ELEM(*vol, k, i, j) * a;
                                        A2D_ELEM(projNorm , y, x) += a2;

#ifdef DEBUG

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

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

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

                                    }
                                }// close possible mask constraint
                            }
                        } // Project this basis

                        if (!FORW)
                        {
                            A3D_ELEM(*vol, k, i, j) += vol_corr;

#ifdef DEBUG

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

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

                // Prepare for next iteration
                V3_PLUS_V3(actprj, actprj, prjX);
            }
            V3_PLUS_V3(beginY, beginY, prjY);
        }
        V3_PLUS_V3(beginZ, beginZ, prjZ * numthreads );
    }
}

projectXraySimpleGridThread()

void projectXraySimpleGridThread

(

ThreadArgument &

thArg

)

Definition at line 726 of file project_xray.cpp.

{

    int threadId = thArg.thread_id;

    auto *dataThread = (XrayThreadArgument*) thArg.data;
    const XRayPSF &psf = *(dataThread->psf);
    MultidimArray<double> *muVol =  (dataThread->muVol);
    MultidimArray<double> *IgeoVol =  (dataThread->IgeoVol);
    Projection *proj = (dataThread->projOut);
    MultidimArray<double> &projNorm = *(dataThread->projNorm);
    int forw = dataThread->forw;

    projectXraySimpleGrid(muVol, psf, IgeoVol, proj, projNorm , forw,
                          threadId, thArg.threads);




}

projectXrayVolume()

void projectXrayVolume	(	MultidimArray< double > &	muVol,
		MultidimArray< double > &	IgeoVol,
		XRayPSF &	psf,
		Projection &	P,
		MultidimArray< double > *	projNorm = nullptr,
		ThreadManager *	ThrMgr = nullptr
	)

Definition at line 334 of file project_xray.cpp.

{

    /* Now we have to split the phantom rotated volume into slabs, taking into
     * account the possible different resolution and Z size.
     *
     * PhantomSlabIdx are the indexes of rotVol defining the slabs according to
     * the splitting done to PSFVol.
     *
     * psfSlicesIdx are the indexes of the z planes of rotVol whose psf are
     * used for the slabs.
     */
    std::vector<int> phantomSlabIdx, psfSlicesIdx;

    // Search for the PSFslab of the beginning of the volume
    auto firstSlab = (int)(STARTINGZ(muVol)*psf.dzo/psf.dzoPSF);

    if (!XMIPP_EQUAL_ZERO(psf.slabThr))
    {
        for (size_t kk = 0; kk < psf.slabIndex.size(); ++kk)
        {
            if (firstSlab < psf.slabIndex[kk])
            {
                firstSlab = kk-1;
                break;
            }
        }

        // Searching the equivalent index in rotvol for the indexes for the Slabs of PSFVol
        phantomSlabIdx.push_back(STARTINGZ(muVol));

        for (size_t kk = firstSlab+1; kk < psf.slabIndex.size(); ++kk)
        {
            auto tempK = (int)(psf.slabIndex[kk] * psf.dzoPSF / psf.dzo);

            if (tempK <= FINISHINGZ(muVol))
            {
                phantomSlabIdx.push_back(tempK);
                int tempKK = psf.slabIndex[kk-1];
                auto psfMeanSlice = (int)((tempK + tempKK)*0.5 * psf.dzoPSF / psf.dzo);
                psfSlicesIdx.push_back(psfMeanSlice);
            }
            else
            {
                phantomSlabIdx.push_back(FINISHINGZ(muVol));
                int psfMeanSlice = (tempK + phantomSlabIdx[kk-1])/2;
                psfSlicesIdx.push_back(psfMeanSlice);
                continue;
            }
        }
    }
    else
    {
        for (int k = STARTINGZ(muVol); k <= FINISHINGZ(muVol); ++k)
        {
            phantomSlabIdx.push_back(k);
            psfSlicesIdx.push_back(k);
        }
        // To define the last range
        phantomSlabIdx.push_back(FINISHINGZ(muVol)+1);
    }

    XrayThreadArgument projThrData;
    Barrier barrier(thMgr->threads);

    projThrData.psf = &psf;
    projThrData.muVol = &muVol;
    projThrData.IgeoVol = &IgeoVol;
    projThrData.projOut = &P;
    projThrData.projNorm = projNorm;
    projThrData.phantomSlabIdx = & phantomSlabIdx;
    projThrData.psfSlicesIdx = & psfSlicesIdx;

    //Create the job handler to distribute thread jobs
    size_t blockSize, numberOfJobs = psfSlicesIdx.size() ;
    blockSize = 1;
    ThreadTaskDistributor td(numberOfJobs, blockSize);
    projThrData.td = &td;
    projThrData.barrier = &barrier;

    //the really really final project routine, I swear by Snoopy.
    //    project_xr(psf,side.rotPhantomVol,P, idxSlice);

    thMgr->run(threadXrayProject,(void*) &projThrData);
    //    P.write("projection_after_threads.spi");
}

threadXrayProject()

void threadXrayProject

(

ThreadArgument &

thArg

)

Thread Job to generate an X-ray microscope projection.

Thread Job to generate an X-ray microscope projection.

Calculate projNorm

Definition at line 440 of file project_xray.cpp.

{

    int thread_id = thArg.thread_id;

    auto *dataThread = (XrayThreadArgument*) thArg.data;
    const XRayPSF &psf = *(dataThread->psf);
    MultidimArray<double> &muVol =  *(dataThread->muVol);
    MultidimArray<double> &IgeoVol =  *(dataThread->IgeoVol);
    Image<double> &imOutGlobal = *(dataThread->projOut);
    MultidimArray<double> * projNorm = dataThread->projNorm;
    std::vector<int> &phantomSlabIdx = *(dataThread->phantomSlabIdx);
    std::vector<int> &psfSlicesIdx = *(dataThread->psfSlicesIdx);
    ParallelTaskDistributor * td = dataThread->td;
    Barrier * barrier = dataThread->barrier;

    size_t first = 0, last = 0;

    if (thread_id == 0)
    {
        muVol.setXmippOrigin();
        MULTIDIM_ARRAY(imOutGlobal).resizeNoCopy(psf.Niy, psf.Nix);
        MULTIDIM_ARRAY(imOutGlobal).initZeros();
        MULTIDIM_ARRAY(imOutGlobal).setXmippOrigin();
    }

    barrier->wait();

    MultidimArray<double> imOut(psf.Niy, psf.Nix), projNormTemp(psf.Niy, psf.Nix);
    imOut.setXmippOrigin();
    projNormTemp.setXmippOrigin();

    MultidimArray<double> imTemp(psf.Noy, psf.Nox),intExp(psf.Noy, psf.Nox),imTempSc(imOut),*imTempP=nullptr;
    intExp.setXmippOrigin();
    imTemp.setXmippOrigin();
    imTempSc.setXmippOrigin();

    // Parallel thread jobs
    while (td->getTasks(first, last))
    {
        //#define DEBUG
#ifdef DEBUG
        Image<double> _Im;
#endif

        int kini = phantomSlabIdx[first];
        int kend = phantomSlabIdx[last + 1] - 1 ;

        //        double deltaZSlab =  psf.dzo*(kend-kini+1);

        intExp.initZeros();

        FOR_ALL_ELEMENTS_IN_ARRAY2D(intExp)
        {
            for (int k = kini; k <= kend; k++)
            {
                double tempValue = A3D_ELEM(muVol,k,i,j);
                A2D_ELEM(intExp,i,j) += tempValue;
            }
            A2D_ELEM(imTemp,i,j) = A3D_ELEM(IgeoVol,kini,i,j)*(1 - exp( -A2D_ELEM(intExp,i,j) * psf.dzo));
        }
#ifdef DEBUG
        _Im().alias(intExp);
        _Im.write("projectXR-intExp.spi");

        _Im().alias(imTemp);
        _Im.write("projectXR-imTemp.spi");
#endif

        switch (psf.AdjustType)
        {
        case PSFXR_INT:
            imTempP = &imTempSc;
            scaleToSize(xmipp_transformation::LINEAR,*imTempP,imTemp,psf.Nix,psf.Niy);
            break;

        case PSFXR_STD:
            imTempP = &imTemp;
            break;

        case PSFXR_ZPAD:
            //    (*imTempSc).resize(imTemp);
            imTempP = &imTempSc;
            imTemp.window(*imTempP,-ROUND(psf.Niy/2)+1,-ROUND(psf.Nix/2)+1,ROUND(psf.Niy/2)-1,ROUND(psf.Nix/2)-1);
            break;
        }

#ifdef DEBUG
        _Im().alias(*imTempP);
        _Im.write("projectXR-imTempEsc_before.spi");
#endif

        psf.applyOTF(*imTempP, imOutGlobal.Zoff()- psfSlicesIdx[first]*psf.dzo);

#ifdef DEBUG

        _Im.write("projectXR-imTempEsc_after.spi");
#endif

        // Adding to the thread temporal imOut
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(*imTempP)
        dAij(imOut,i,j) += dAij(*imTempP,i,j);

#ifdef DEBUG

        _Im().alias(imOut);
        _Im.write("projectXR-imout.spi");
#endif


        if (projNorm != nullptr)
        {

            FOR_ALL_ELEMENTS_IN_ARRAY2D(imTemp)
            A2D_ELEM(imTemp,i,j) = A3D_ELEM(IgeoVol,kini,i,j) * psf.dzo;

            //            IgeoVol.getSlice(kini, imTemp);
            psf.applyOTF(imTemp, imOutGlobal.Zoff()- psfSlicesIdx[first]*psf.dzo);

            FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(imTemp)
            dAij(projNormTemp,i,j) += dAij(imTemp,i,j)*dAij(imTemp,i,j);
        }


    }

    //Lock to update the total addition and multiply for the pixel width (to avoid multiply several times)
    mutex.lock();
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(imOut)
    dAij(MULTIDIM_ARRAY(imOutGlobal),i,j) += dAij(imOut,i,j);

    if (projNorm != nullptr)
    {
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(projNormTemp)
        dAij(*projNorm,i,j) += dAij(projNormTemp,i,j);
    }
    mutex.unlock();
#ifdef DEBUG

    imOutGlobal.write("projectXR-imoutGlobal.spi");
#endif

    barrier->wait();

    if (thread_id==0)
    {
        // The output is the addition of the  accumulated differences, the real projection is background illumination minus this
        //        MULTIDIM_ARRAY(imOutGlobal) = 1.0- MULTIDIM_ARRAY(imOutGlobal);

        switch (psf.AdjustType)
        {
        case PSFXR_STD:
            // Do nothing;
            break;
        case PSFXR_INT:
            selfScaleToSize(xmipp_transformation::LINEAR,imOutGlobal(), psf.Nox, psf.Noy);
            break;
        case PSFXR_ZPAD:
            MULTIDIM_ARRAY(imOutGlobal).selfWindow(-ROUND(psf.Noy/2)+1,-ROUND(psf.Nox/2)+1,ROUND(psf.Noy/2)-1,ROUND(psf.Nox/2)-1);
            break;
        }

        MULTIDIM_ARRAY(imOutGlobal).setXmippOrigin();

#ifdef DEBUG

        imOutGlobal.write("projectXR-imoutGlobal_negative.spi");
#endif

    }
    //    std::cerr << "th" << thread_id << ": Thread Finished" <<std::endl;

#undef DEBUG
}

XrayRotateAndProjectVolumeOffCentered()

void XrayRotateAndProjectVolumeOffCentered	(	XrayProjPhantom &	side,
		XRayPSF &	psf,
		Projection &	P,
		Projection &	standardP,
		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.

Off-centered not implemented. Rotations are around volume center

Calculation of IgeoVol to optimize process

Definition at line 210 of file project_xray.cpp.

{

    int iniXdim, iniYdim, iniZdim, newXdim, newYdim, newZdim, rotXdim, rotZdim;

    iniXdim = XSIZE(MULTIDIM_ARRAY(phantom.iniVol));
    iniYdim = YSIZE(MULTIDIM_ARRAY(phantom.iniVol));
    iniZdim = ZSIZE(MULTIDIM_ARRAY(phantom.iniVol));

    // Projection offset in pixels
    Matrix1D<double> offsetNV(3);
    P.getShifts(XX(offsetNV), YY(offsetNV), ZZ(offsetNV));

    //    xOffsetN = P.Xoff()/psf.dxo;
    //    yOffsetN = P.Yoff()/psf.dxo;


    /* If Xdim is greater than Zdim, then as we rotate the volume the rotated Zdim must be
     * great enough to cover all the volume.
     */
    /* By now, we are going to keep it unused */
    //    if (true)
    //    {
    //        double radi, theta0, theta, tilt;
    //
    //        radi = sqrt(iniZdim*iniZdim + iniXdim*iniXdim);
    //        theta0 = atan2(iniZdim, iniXdim);
    //        tilt = ((P.tilt() < 90)? P.tilt(): 180 - P.tilt())* PI / 180;
    //
    //        rotZdim = XMIPP_MAX(ROUND(radi * sin(ABS(tilt) + ABS(theta0))), iniZdim);
    //        rotXdim = XMIPP_MAX(ROUND(radi * cos(ABS(tilt) - ABS(theta0))), iniXdim);
    //        rotYdim = iniYdim;
    //
    //        //    rotZdim = iniZdim;
    //        //    rotXdim = iniXdim;
    //
    //        newXdim = rotXdim + 2*ABS(XX(offsetNV));
    //        newYdim = iniYdim + 2*ABS(YY(offsetNV));
    //        newZdim = rotZdim;
    //    }
    //    else
    rotXdim = iniXdim;
    rotZdim = iniZdim;

    newXdim = (int)(rotXdim + 2*fabs(XX(offsetNV)));
    newYdim = (int)(iniYdim + 2*fabs(YY(offsetNV)));
    newZdim = rotZdim;

    // We set the dimensions only to obtain the values of starting X,Y,Z
    //    phantom.rotVol.setDimensions(rotXdim, rotYdim, rotZdim, 1);
    phantom.rotVol.setDimensions(newXdim, newYdim, newZdim, 1);
    phantom.rotVol.setXmippOrigin();

    if (psf.verbose > 1)
    {
        if (newXdim != iniXdim || newYdim != iniYdim)
        {
            std::cout << std::endl;
            std::cout << "--------------------------------" << std::endl;
            std::cout << "XrayProject::Volume_offCentered:" << std::endl;
            std::cout << "--------------------------------" << std::endl;
            std::cout << "(X,Y,Z) shifts = (" << XX(offsetNV) << "," << YY(offsetNV) << ","
            << ZZ(offsetNV) << ") um" << std::endl;
            std::cout << "Image resize (Nx,Ny): (" << iniXdim << "," << iniYdim << ") --> ("
            << newXdim << "," << newYdim << ") " << std::endl;
        }
#ifdef DEBUG

        std::cout <<"yoffsetN "<< YY(offsetNV) <<std::endl;
        std::cout <<"xoffsetN "<< XX(offsetNV) <<std::endl;
        std::cout <<"yinit    " << yinit  <<std::endl;
        std::cout <<"yend     "    << yend  <<std::endl;
        std::cout <<"xinit    "   << xinit  <<std::endl;
        std::cout <<"xend     "    << xend  <<std::endl;
        std::cout <<"zinit    "   << zinit  <<std::endl;
        std::cout <<"zend     "    << zend  <<std::endl;
#endif

    }



    //    phantom.rotVol.setMmap(true);
    phantom.rotVol.resizeNoCopy(newZdim, newYdim, newXdim);
    phantom.rotVol.setXmippOrigin();

    // Rotate volume ....................................................

    Matrix2D<double> R, T;

    ZZ(offsetNV) = 0; // We are not interested in apply this Zshift
    translation3DMatrix(offsetNV, T);
    Euler_angles2matrix(P.rot(), P.tilt(), P.psi(), R, true);

    double outside = 0; //phantom.iniVol.getPixel(0,0,0,0);
    MULTIDIM_ARRAY(phantom.iniVol).setXmippOrigin();

    applyGeometry(xmipp_transformation::LINEAR, phantom.rotVol, MULTIDIM_ARRAY(phantom.iniVol), T*R, xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP, outside);

    psf.adjustParam(phantom.rotVol);


    MultidimArray<double> IgeoVol, IgeoZb;
    IgeoZb.resize(1, 1, YSIZE(phantom.rotVol), XSIZE(phantom.rotVol),false);
    IgeoZb.initConstant(1.);

    calculateIgeo(phantom.rotVol, psf.dzo, IgeoVol, MULTIDIM_ARRAY(standardP), IgeoZb, psf.nThr, thMgr);



    projectXrayVolume(phantom.rotVol, IgeoVol, psf, P, nullptr, thMgr);

    int outXDim = XMIPP_MIN(Xdim,iniXdim);
    int outYDim = XMIPP_MIN(Ydim,iniYdim);

    P().selfWindow(-ROUND(outYDim/2),
                   -ROUND(outXDim/2),
                   -ROUND(outYDim/2) + outYDim -1,
                   -ROUND(outXDim/2) + outXDim -1);

}//XrayRotateAndProjectVolumeOffCentered

Variable Documentation

barrier

Barrier* barrier

Definition at line 77 of file project_xray.cpp.

mutex

Mutex mutex

Definition at line 76 of file project_xray.cpp.

thMgr

ThreadManager* thMgr

Definition at line 78 of file project_xray.cpp.