reconstruct_art_pseudo.cpp File Reference

#include <fstream>
#include "reconstruct_art_pseudo.h"
#include "core/histogram.h"
#include "core/geometry.h"
#include "core/xmipp_image.h"
#include "core/xmipp_funcs.h"

Include dependency graph for reconstruct_art_pseudo.cpp:

Go to the source code of this file.

Functions
void	project_Pseudo (const std::vector< Matrix1D< double > > &atomPosition, std::vector< double > &atomWeight, double sigma, MultidimArray< double > &proj, MultidimArray< double > &norm_proj, Matrix2D< double > &Euler, double shiftX, double shiftY, const std::vector< double > &lambda, const std::vector< Matrix2D< double > > &NMA, int direction, const Matrix1D< double > &gaussianProjectionTable, const Matrix1D< double > &gaussianProjectionTable2)

Variables
constexpr int	FORWARD = 1
constexpr int	BACKWARD = -1

Function Documentation

project_Pseudo()

void project_Pseudo	(	const std::vector< Matrix1D< double > > &	atomPosition,
		std::vector< double > &	atomWeight,
		double	sigma,
		MultidimArray< double > &	proj,
		MultidimArray< double > &	norm_proj,
		Matrix2D< double > &	Euler,
		double	shiftX,
		double	shiftY,
		const std::vector< double > &	lambda,
		const std::vector< Matrix2D< double > > &	NMA,
		int	direction,
		const Matrix1D< double > &	gaussianProjectionTable,
		const Matrix1D< double > &	gaussianProjectionTable2
	)

Projection of a pseudoatom volume

Definition at line 38 of file reconstruct_art_pseudo.cpp.

{
    // Project all pseudo atoms ............................................
    int nmax=atomPosition.size();
    Matrix1D<double> actprj(3);
    double sigma4=4*sigma;
    Matrix1D<double> actualAtomPosition;
    int lambdaSize=lambda.size();
    for (int n=0; n<nmax; n++)
    {
        actualAtomPosition=atomPosition[n];
        double weight=atomWeight[n];
        for (int mode=0; mode<lambdaSize; mode++)
        {
            const Matrix2D<double> &NMAmode=NMA[mode];
            double lambdam=lambda[mode];
            FOR_ALL_ELEMENTS_IN_MATRIX1D(actualAtomPosition)
            VEC_ELEM(actualAtomPosition,i)+=
                lambdam*MAT_ELEM(NMAmode,n,i);
        }

        Uproject_to_plane(actualAtomPosition, Euler, actprj);
        XX(actprj)+=shiftX;
        YY(actprj)+=shiftY;
#ifdef DEBUG

        bool condition = true;
        if (condition)
        {
            std::cout << "Projecting point " << atomPositions[n].transpose() << std::endl;
            std::cout << "Vol there = " << atomWeight[n] << std::endl;
            std::cout << " Center of the basis proj (2D) " << XX(actprj) << "," << YY(actprj) << std::endl;
        }
#endif

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

#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
        if (XX_corner1 <= XX_corner2 && YY_corner1 <= YY_corner2)
        {
            double vol_corr=0;

            // Effectively project this basis
            for (int y = YY_corner1; y <= YY_corner2; y++)
            {
                double y_diff2=y-YY(actprj);
                y_diff2=y_diff2*y_diff2;
                for (int x = XX_corner1; x <= XX_corner2; x++)
                {
                    double x_diff2=x-XX(actprj);
                    x_diff2=x_diff2*x_diff2;
                    double r=sqrt(x_diff2+y_diff2);
                    double didx=r*1000;
                    int idx=ROUND(didx);
                    double a  = VEC_ELEM(gaussianProjectionTable,idx);
                    double a2 = VEC_ELEM(gaussianProjectionTable2,idx);
#ifdef DEBUG

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

                    if (direction==FORWARD)
                    {
                        A2D_ELEM(proj, y, x) += weight * a;
                        A2D_ELEM(norm_proj, y, x) += a2;
#ifdef DEBUG

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

                    }
                    else
                    {
                        vol_corr += A2D_ELEM(norm_proj, y, x) * a;
#ifdef DEBUG

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

                    }
                }
            }

            if (direction==BACKWARD)
            {
                atomWeight[n] += vol_corr;
#ifdef DEBUG

                if (condition)
                {
                    std::cout << "\nFinal value at ( " << n << ") = "
                    << atomWeight[n] << std::endl;
                }
#endif

            }
        } // If not collapsed
    }
}

Variable Documentation

BACKWARD

constexpr int BACKWARD = -1

Definition at line 35 of file reconstruct_art_pseudo.cpp.

FORWARD

constexpr int FORWARD = 1

Definition at line 34 of file reconstruct_art_pseudo.cpp.