tomo_detect_missing_wedge.cpp File Reference

#include "tomo_detect_missing_wedge.h"
#include "core/geometry.h"
#include "core/matrix2d.h"
#include "core/xmipp_fftw.h"
#include "data/numerical_tools.h"

Include dependency graph for tomo_detect_missing_wedge.cpp:

Go to the source code of this file.

Classes
class	WedgeSolver

Functions
double	evaluatePlane (double rot, double tilt, const MultidimArray< double > *V, const MultidimArray< double > *Vmag, double maxFreq, double planeWidth, int direction, MultidimArray< double > *Vdraw=nullptr, bool setPos=false, double rotPos=0, double tiltPos=0)
double	wrapperFitnessDetectMissingWedge (double *p, void *extraArgs)
void	lookForPlane (const MultidimArray< double > *V, const MultidimArray< double > *Vmag, double maxFreq, double planeWidth, int direction, double &rot, double &tilt, bool setPos=false, double rotPos=0, double tiltPos=0)
void	drawWedge (double rotPos, double tiltPos, double rotNeg, double tiltNeg, const MultidimArray< double > *V, const MultidimArray< double > *Vmag, MultidimArray< double > *Vdraw)

Function Documentation

drawWedge()

void drawWedge	(	double	rotPos,
		double	tiltPos,
		double	rotNeg,
		double	tiltNeg,
		const MultidimArray< double > *	V,
		const MultidimArray< double > *	Vmag,
		MultidimArray< double > *	Vdraw
	)

Definition at line 245 of file tomo_detect_missing_wedge.cpp.

{
    Matrix2D<double> Epos, Eneg;
    Euler_angles2matrix(rotPos,tiltPos,0.,Epos);
    Euler_angles2matrix(rotNeg,tiltNeg,0.,Eneg);

    Matrix1D<double> freq(3), freqPos, freqNeg;
    Matrix1D<int> idx(3);
    Vdraw->initZeros(*Vmag);
    FOR_ALL_ELEMENTS_IN_ARRAY3D(*Vdraw)
    {
        // Frequency in the coordinate system of the volume
        VECTOR_R3(idx,j,i,k);
        FFT_idx2digfreq(*V,idx,freq);

        // Frequency in the coordinate system of the plane
        freqPos=Epos*freq;
        freqNeg=Eneg*freq;
        if (ZZ(freqPos)<0 || ZZ(freqNeg)>0)
            (*Vdraw)(k,i,j)+=1;
    }
}

evaluatePlane()

double evaluatePlane	(	double	rot,
		double	tilt,
		const MultidimArray< double > *	V,
		const MultidimArray< double > *	Vmag,
		double	maxFreq,
		double	planeWidth,
		int	direction,
		MultidimArray< double > *	Vdraw = nullptr,
		bool	setPos = false,
		double	rotPos = 0,
		double	tiltPos = 0
	)

Definition at line 33 of file tomo_detect_missing_wedge.cpp.

{
    if (rot<0 || rot>360 || tilt<-90 || tilt>90)
        return 0;

    Matrix2D<double> E, Einv;
    Euler_angles2matrix(rot,tilt,0.,E);
    Einv=E.transpose();

    if (setPos)
    {
        Matrix2D<double> Epos;
        Euler_angles2matrix(rotPos,tiltPos,0.,Epos);
        double angle=acos(E(2,0)*Epos(2,0)+E(2,1)*Epos(2,1)+E(2,2)*Epos(2,2));
        angle=RAD2DEG(angle);
        if (fabs(angle)<20 || fabs(180-angle)<20)
            return 0;
    }

    size_t N=XMIPP_MAX(XSIZE(*Vmag),YSIZE(*Vmag)/2);
    N=XMIPP_MAX(N,ZSIZE(*Vmag)/2);
    double df=0.5/N;
    Matrix1D<double> freq(3), freqp(3);
    Matrix1D<int> idx(3);
    double sumNeg=0, sumPos=0;
    int Nneg=0, Npos=0;
    double maxFreq2=maxFreq*maxFreq;
    auto iPlaneWidth=(int)ceil(planeWidth);
    for (int ix=0; ix<=N; ix++)
    {
        XX(freq)=ix*df;
        double fx2=XX(freq)*XX(freq);
        if (fx2>maxFreq2)
            continue;
        for (int iy=-(int)N; iy<=N; iy++)
        {
            YY(freq)=iy*df;
            double fx2fy2=fx2+YY(freq)*YY(freq);
            if (fx2fy2>maxFreq2)
                continue;
            for (int iz=-iPlaneWidth; iz<=iPlaneWidth; iz++)
            {
                if (iz==0 || ix==0 || iy==0)
                    continue;

                // Frequency in the coordinate system of the plane
                ZZ(freq)=iz*df;

                // Frequency in the coordinate system of the volume
                SPEED_UP_temps012;
                M3x3_BY_V3x1(freqp,Einv,freq);
                bool inverted=false;
                if (XX(freqp)<0)
                {
                    XX(freqp)=-XX(freqp);
                    YY(freqp)=-YY(freqp);
                    ZZ(freqp)=-ZZ(freqp);
                    inverted=true;
                }

                // Get the corresponding index
                DIGFREQ2FFT_IDX(ZZ(freqp), ZSIZE(*V), ZZ(idx));
                DIGFREQ2FFT_IDX(YY(freqp), YSIZE(*V), YY(idx));
                DIGFREQ2FFT_IDX(XX(freqp), XSIZE(*V), XX(idx));
                if (XX(idx) < STARTINGX(*Vmag) || XX(idx) > FINISHINGX(*Vmag) ||
                    YY(idx) < STARTINGY(*Vmag) || YY(idx) > FINISHINGY(*Vmag) ||
                    ZZ(idx) < STARTINGZ(*Vmag) || ZZ(idx) > FINISHINGZ(*Vmag))
                    continue;

                // Make the corresponding sums
                bool negativeSum;
                if (direction==1)
                    negativeSum=iz<0;
                else
                    negativeSum=iz>0;
                double val=A3D_ELEM(*Vmag,ZZ(idx),YY(idx),XX(idx));
                if ((negativeSum && !inverted) || (!negativeSum && inverted)) // XOR
                {
                    sumNeg+=val;
                    Nneg++;
                    if (Vdraw!=nullptr)
                        (*Vdraw)(idx)=2*direction*val;
                }
                else
                {
                    sumPos+=val;
                    Npos++;
                    if (Vdraw!=nullptr)
                        (*Vdraw)(idx)=1.0/2.0*direction*val;
                }
            }
        }
    }
    if (fabs(Nneg-Npos)/(0.5*(Nneg+Npos))>0.5)
        // If there is a difference of more than 50%
        return 1e38;
    if (Nneg!=0)
        sumNeg/=Nneg;
    else
        return 1e38;
    if (Npos!=0)
        sumPos/=Npos;
    else
        return 1e38;

    return -(sumPos-sumNeg);
}

lookForPlane()

void lookForPlane	(	const MultidimArray< double > *	V,
		const MultidimArray< double > *	Vmag,
		double	maxFreq,
		double	planeWidth,
		int	direction,
		double &	rot,
		double &	tilt,
		bool	setPos = false,
		double	rotPos = 0,
		double	tiltPos = 0
	)

Definition at line 202 of file tomo_detect_missing_wedge.cpp.

{
    // Optimize with DE
    int length=2;
    int Npop=50;

    double min_allowed[] = {0,-90};
    double max_allowed[] = {360,90};

    auto * solver = new WedgeSolver(length,length*Npop,V,Vmag,maxFreq,planeWidth,direction);
    solver->Setup( min_allowed, max_allowed, stBest2Bin, 0.5, 0.8);

    if (setPos)
    {
        solver->setPos=true;
        solver->rotPos=rotPos;
        solver->tiltPos=tiltPos;
    }

    solver->Solve(50);

    double * bestSolution = solver->Solution();

    // Refine with Powell
    Matrix1D<double> x(2), steps(2);
    double fitness;
    int iter;
    steps.initConstant(1);
    x(0)=bestSolution[0];
    x(1)=bestSolution[1];
    powellOptimizer(x,1,2,&wrapperFitnessDetectMissingWedge,solver,0.01,fitness,iter,steps,true);
    rot=x(0);
    tilt=x(1);

    delete solver;

    solver = nullptr;
}

wrapperFitnessDetectMissingWedge()

double wrapperFitnessDetectMissingWedge	(	double *	p,
		void *	extraArgs
	)

Definition at line 195 of file tomo_detect_missing_wedge.cpp.

{
    bool dummy;
    auto *wegde_solver=(WedgeSolver *) extraArgs;
    return wegde_solver->EnergyFunction(p+1,dummy);
}