project_real_shears.cpp File Reference

#include "project_real_shears.h"
#include "core/bilib/getputd.h"
#include "core/bilib/changebasis.h"
#include "core/bilib/configs.h"
#include "core/bilib/kernel.h"
#include "core/transformations.h"
#include "data/fourier_projection.h"

Include dependency graph for project_real_shears.cpp:

Go to the source code of this file.

Macros
#define	DBL_EPSILON   2.2204460492503130808472633361816000000000e-16

Functions
void	convertAngles (double &phi, double &theta, double &psi)
	Transforms angles from (Ry, Rz, Ry) to (Rx, Ry, Rz) system. Returns possible error. More...
void	projectionRealShears2 (MultidimArray< double > &CoefVolume, double absscale, Matrix2D< double > &Binv, Matrix1D< double > &BinvCscaled, int *arr, MultidimArray< double > &projection)
void	projectionRealShears1 (RealShearsInfo &Data, double phi, double theta, double psi, double shiftX, double shiftY, MultidimArray< double > &projection)
void	projectVolume (RealShearsInfo &Data, Projection &P, int Ydim, int Xdim, double rot, double tilt, double psi, double shiftX, double shiftY)
	Make projection. More...

Macro Definition Documentation

DBL_EPSILON

#define DBL_EPSILON   2.2204460492503130808472633361816000000000e-16

Definition at line 40 of file project_real_shears.cpp.

Function Documentation

convertAngles()

void convertAngles	(	double &	phi,
		double &	theta,
		double &	psi
	)

Transforms angles from (Ry, Rz, Ry) to (Rx, Ry, Rz) system. Returns possible error.

Definition at line 127 of file project_real_shears.cpp.

{
    Matrix2D<double> E;
    Euler_angles2matrix(phi,theta,psi,E,false);
    double A00=MAT_ELEM(E,0,0);
    double A02=MAT_ELEM(E,0,2);
    double A10=MAT_ELEM(E,1,0);
    double A12=MAT_ELEM(E,1,2);
    double A20=MAT_ELEM(E,2,0);
    double A21=MAT_ELEM(E,2,1);
    double A22=MAT_ELEM(E,2,2);

    double abs_cosay = sqrt(A22*A22+A21*A21);

    //Results angles
    double ax, ay, az;

    if(abs_cosay > 0.0)
    {
        double sign_cosay;

        ax = atan2(-A21, A22);
        az = atan2(-A10, A00);

        if(fabs(cos(ax)) == 0.0)
            sign_cosay = SGN(-A21/sin(ax));
        else
        {
            if (cos(ax) > 0.0)
                sign_cosay = SGN(A22);
            else
                sign_cosay = -SGN(A22);
        }

        ay  = atan2(A20, sign_cosay * abs_cosay);
    }
    else
    {
        //Let's consider the matrix as a rotation around Z
        if (SGN(A20) > 0.0)
        {
            ax = 0.0;
            ay  = PI/2.0;
            az = atan2(A12, -A02);
        }
        else
        {
            ax = 0.0;
            ay  = -PI/2.0;
            az = atan2(-A12, A02);
        }
    }

    phi   = ax;
    theta = ay;
    psi   = az;
}

projectionRealShears1()

void projectionRealShears1	(	RealShearsInfo &	Data,
		double	phi,
		double	theta,
		double	psi,
		double	shiftX,
		double	shiftY,
		MultidimArray< double > &	projection
	)

Definition at line 267 of file project_real_shears.cpp.

{
    Matrix2D<double> R;
    R.initIdentity(4);
    MAT_ELEM(R,0,3)=shiftX;
    MAT_ELEM(R,1,3)=shiftY;
    Matrix2D<double> B=Data.Ac*R;
    rotation3DMatrix(-RAD2DEG(phi),'X',R,true);
    B=B*R;
    rotation3DMatrix(RAD2DEG(theta),'Y',R,true);
    B=B*R;
    rotation3DMatrix(-RAD2DEG(psi),'Z',R,true);
    B=B*R;
    B=B*Data.Acinv;
    Matrix2D<double> Binv;
    B.inv(Binv);
    Matrix1D<double> BinvC;
    Binv.getCol(2,BinvC);

    double  scale_x, scale_y, scale_z;
    if (XX(BinvC) != 0.0)
        scale_x = 1.0 / XX(BinvC);
    else
        scale_x = 1e38;

    if (YY(BinvC) != 0.0)
        scale_y = 1.0 / YY(BinvC);
    else
        scale_y = 1e38;

    if (ZZ(BinvC) != 0.0)
        scale_z = 1.0 / ZZ(BinvC);
    else
        scale_z = 1e38;
    double m_x = fabs(scale_x);
    double m_y = fabs(scale_y);
    double m_z = fabs(scale_z);

    int     arr[3];
    double minm = m_x;
    double scale = scale_x;
    arr[0] = 0;
    arr[1] = 1;
    arr[2] = 2;
    MultidimArray<double> *Coef_xyz = &Data.Coef_x;
    if (m_y < minm)
    {
        minm = m_y;
        scale = scale_y;
        arr[0] = 1;
        arr[1] = 0;
        arr[2] = 2;
        Coef_xyz = &Data.Coef_y;
    }
    if (m_z < minm)
    {
        minm = m_z;
        scale = scale_z;
        arr[0] = 2;
        arr[1] = 0;
        arr[2] = 1;
        Coef_xyz = &Data.Coef_z;
    }
    Matrix1D<double> BinvCscaled=BinvC;
    BinvCscaled*=scale;

    projectionRealShears2(*Coef_xyz, minm, Binv, BinvCscaled, arr, projection);
}

projectionRealShears2()

void projectionRealShears2	(	MultidimArray< double > &	CoefVolume,
		double	absscale,
		Matrix2D< double > &	Binv,
		Matrix1D< double > &	BinvCscaled,
		int *	arr,
		MultidimArray< double > &	projection
	)

Definition at line 187 of file project_real_shears.cpp.

{
    double  K[4], Arg[4], Ki[4];
    int Xdim=XSIZE(projection);
    double *ptrProjection=MULTIDIM_ARRAY(projection);
    size_t CC1 = Xdim * Xdim;

    Ki[0]=dMn(Binv,3);
    Ki[1]=dMn(Binv,7);
    Ki[2]=dMn(Binv,11);
    Ki[3]=dMn(Binv,15);
    for (int i = 0; i < Xdim; i++)
    {
        K[0]=Ki[0];
        K[1]=Ki[1];
        K[2]=Ki[2];
        K[3]=Ki[3];
        for (int n = 0; n < Xdim; n++)
        {
            double Proj = 0.0;
            for (int ksi = 0; ksi < Xdim; ksi++)
            {
                size_t CC2 = CC1 * ksi;
                double sc = (double) ksi - K[arr[0]];
                Arg[0]=VEC_ELEM(BinvCscaled,0)*sc+K[0];
                Arg[1]=VEC_ELEM(BinvCscaled,1)*sc+K[1];
                Arg[2]=VEC_ELEM(BinvCscaled,2)*sc+K[2];
                Arg[3]=VEC_ELEM(BinvCscaled,3)*sc+K[3];
                double g = Arg[arr[1]];
                double h = Arg[arr[2]];

                auto l1 = (int)ceil(g-2.0);
                int l2 = l1 + 3;
                auto m1 = (int)ceil(h-2.0);
                int m2 = m1 + 3;
                double columns = 0.0;
                double aux;
                for (int m = m1; m <= m2; m++)
                {
                    if (m < Xdim && m > -1)
                    {
                        size_t CC3 = CC2 + Xdim * m;
                        double rows = 0.0;
                        for (int l = l1; l <= l2; l++)
                        {
                            if (l < Xdim && l > -1)
                            {
                                double gminusl = g - (double) l;
                                BSPLINE03(aux,gminusl);
                                double Coeff = DIRECT_A1D_ELEM(CoefVolume,CC3 + l);
                                rows += Coeff * aux;
                            }
                        }
                        double hminusm = h - (double) m;
                        BSPLINE03(aux,hminusm);
                        columns +=  rows * aux;
                    }
                }
                Proj += columns;
            }
            *ptrProjection++ = Proj*absscale;
            K[0]+=dMn(Binv,0);
            K[1]+=dMn(Binv,4);
            K[2]+=dMn(Binv,8);
            K[3]+=dMn(Binv,12);
        }
        Ki[0]+=dMn(Binv,1);
        Ki[1]+=dMn(Binv,5);
        Ki[2]+=dMn(Binv,9);
        Ki[3]+=dMn(Binv,13);
    }
}

projectVolume()

void projectVolume	(	RealShearsInfo &	Data,
		Projection &	P,
		int	Ydim,
		int	Xdim,
		double	rot,
		double	tilt,
		double	psi,
		double	shiftX,
		double	shiftY
	)

Make projection.

Definition at line 339 of file project_real_shears.cpp.

{
    P.reset(Data.Xdim,Data.Xdim);
    double Phi=-psi;
    double Theta=-tilt;
    double Psi=-rot;
    convertAngles(Phi, Theta, Psi);
    projectionRealShears1(Data, Phi, Theta, Psi, shiftX, shiftY, P());
    if (Ydim!=Data.Xdim || Xdim!=Data.Xdim)
        P().selfWindow(FIRST_XMIPP_INDEX(Ydim),FIRST_XMIPP_INDEX(Xdim),
                       LAST_XMIPP_INDEX(Ydim),LAST_XMIPP_INDEX(Xdim));
}