angular_continuous_assign2.cpp File Reference

#include "angular_continuous_assign2.h"
#include "core/transformations.h"
#include "core/xmipp_image_extension.h"
#include "core/xmipp_image_generic.h"
#include "data/mask.h"

Include dependency graph for angular_continuous_assign2.cpp:

Go to the source code of this file.

Functions
double	tranformImage (ProgAngularContinuousAssign2 *prm, double rot, double tilt, double psi, double a, double b, const Matrix2D< double > &A, double deltaDefocusU, double deltaDefocusV, double deltaDefocusAngle, int degree)
double	continuous2cost (double *x, void *_prm)

Function Documentation

continuous2cost()

double continuous2cost	(	double *	x,
		void *	_prm
	)

Definition at line 348 of file angular_continuous_assign2.cpp.

{
    double a=x[1];
    double b=x[2];
    double deltax=x[3];
    double deltay=x[4];
    double scalex=x[5];
    double scaley=x[6];
    double scaleAngle=x[7];
    double deltaRot=x[8];
    double deltaTilt=x[9];
    double deltaPsi=x[10];
    double deltaDefocusU=x[11];
    double deltaDefocusV=x[12];
    double deltaDefocusAngle=x[13];
    auto *prm=(ProgAngularContinuousAssign2 *)_prm;
    if (prm->maxShift>0 && deltax*deltax+deltay*deltay>prm->maxShift*prm->maxShift)
        return 1e38;
    if (fabs(scalex)>prm->maxScale || fabs(scaley)>prm->maxScale)
        return 1e38;
    if (fabs(deltaRot)>prm->maxAngularChange || fabs(deltaTilt)>prm->maxAngularChange || fabs(deltaPsi)>prm->maxAngularChange)
        return 1e38;
    if (fabs(a-prm->old_grayA)>prm->maxA)
        return 1e38;
    if (fabs(b)>prm->maxB*prm->Istddev)
        return 1e38;
    if (fabs(deltaDefocusU)>prm->maxDefocusChange || fabs(deltaDefocusV)>prm->maxDefocusChange)
        return 1e38;
//  MAT_ELEM(prm->A,0,0)=1+scalex;
//  MAT_ELEM(prm->A,1,1)=1+scaley;
// In Matlab
//  syms sx sy t
//  R=[cos(t) -sin(t); sin(t) cos(t)]
//  S=[1+sx 0; 0 1+sy]
//  simple(transpose(R)*S*R)
//  [ sx - sx*sin(t)^2 + sy*sin(t)^2 + 1,            -sin(2*t)*(sx/2 - sy/2)]
//  [            -sin(2*t)*(sx/2 - sy/2), sy + sx*sin(t)^2 - sy*sin(t)^2 + 1]
    double sin2_t=sin(scaleAngle)*sin(scaleAngle);
    double sin_2t=sin(2*scaleAngle);
    MAT_ELEM(prm->A,0,0)=1+scalex+(scaley-scalex)*sin2_t;
    MAT_ELEM(prm->A,0,1)=0.5*(scaley-scalex)*sin_2t;
    MAT_ELEM(prm->A,1,0)=MAT_ELEM(prm->A,0,1);
    MAT_ELEM(prm->A,1,1)=1+scaley-(scaley-scalex)*sin2_t;
    MAT_ELEM(prm->A,0,2)=prm->old_shiftX+deltax;
    MAT_ELEM(prm->A,1,2)=prm->old_shiftY+deltay;
    return tranformImage(prm,prm->old_rot+deltaRot, prm->old_tilt+deltaTilt, prm->old_psi+deltaPsi,
            a, b, prm->A, deltaDefocusU, deltaDefocusV, deltaDefocusAngle, xmipp_transformation::LINEAR);
}

tranformImage()

double tranformImage	(	ProgAngularContinuousAssign2 *	prm,
		double	rot,
		double	tilt,
		double	psi,
		double	a,
		double	b,
		const Matrix2D< double > &	A,
		double	deltaDefocusU,
		double	deltaDefocusV,
		double	deltaDefocusAngle,
		int	degree
	)

Definition at line 251 of file angular_continuous_assign2.cpp.

{
    if (prm->hasCTF)
    {
        double defocusU=prm->old_defocusU+deltaDefocusU;
        double defocusV;
        if (prm->sameDefocus){
            defocusV=defocusU;
        }
        else
            defocusV=prm->old_defocusV+deltaDefocusV;
        double angle=prm->old_defocusAngle+deltaDefocusAngle;
        if (defocusU!=prm->currentDefocusU || defocusV!=prm->currentDefocusV || angle!=prm->currentAngle)
            prm->updateCTFImage(defocusU,defocusV,angle);
    }
    projectVolume(*(prm->projector), prm->P, (int)XSIZE(prm->I()), (int)XSIZE(prm->I()),  rot, tilt, psi, (const MultidimArray<double> *)prm->ctfImage);
    if (prm->old_flip)
    {
        MAT_ELEM(A,0,0)*=-1;
        MAT_ELEM(A,0,1)*=-1;
        MAT_ELEM(A,0,2)*=-1;
    }

    applyGeometry(degree,prm->Ifilteredp(),prm->Ifiltered(),A,xmipp_transformation::IS_NOT_INV,xmipp_transformation::DONT_WRAP,0.);
    const MultidimArray<double> &mP=prm->P();
    const MultidimArray<int> &mMask2D=prm->mask2D;
    MultidimArray<double> &mIfilteredp=prm->Ifilteredp();
    MultidimArray<double> &mE=prm->E();
    mE.initZeros();
    if (prm->contCost==CONTCOST_L1)
    {
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mMask2D)
        {
            if (DIRECT_MULTIDIM_ELEM(mMask2D,n))
            {
                DIRECT_MULTIDIM_ELEM(mIfilteredp,n)=DIRECT_MULTIDIM_ELEM(mIfilteredp,n);
                double val=(a*DIRECT_MULTIDIM_ELEM(mP,n)+b)-DIRECT_MULTIDIM_ELEM(mIfilteredp,n);
                DIRECT_MULTIDIM_ELEM(mE,n)=val;
            }
            else
                DIRECT_MULTIDIM_ELEM(mIfilteredp,n)=0;
        }
    }
    else
    {
        FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mMask2D)
        {
            if (DIRECT_MULTIDIM_ELEM(mMask2D,n))
            {
                double val=DIRECT_MULTIDIM_ELEM(mP,n)-DIRECT_MULTIDIM_ELEM(mIfilteredp,n);
                DIRECT_MULTIDIM_ELEM(mE,n)=val;
            }
            else
                DIRECT_MULTIDIM_ELEM(mIfilteredp,n)=0;
        }
    }

    double cost=0.0;
    if (prm->contCost==CONTCOST_L1)
    {
        FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(mE)
            cost+=fabs(DIRECT_A2D_ELEM(mE,i,j));
        cost*=prm->iMask2Dsum;
    }
    else
        cost=-correlationIndex(mIfilteredp,mP,&mMask2D);

    //covarianceMatrix(mE, prm->C);
    //double div=computeCovarianceMatrixDivergence(prm->C0,prm->C)/MAT_XSIZE(prm->C);
#ifdef DEBUG
    std::cout << "A=" << A << std::endl;
    Image<double> save;
    save()=a*prm->P()+b;
    save.write("PPPtheo.xmp");
    save()=prm->Ifilteredp();
    save.write("PPPfilteredp.xmp");
    save()=prm->Ifiltered();
    save.write("PPPfiltered.xmp");
    save()=prm->E();
    save.write("PPPe.xmp");
    //save()=prm->C;
    //save.write("PPPc.xmp");
    //save()=prm->C0;
    //save.write("PPPc0.xmp");
    //std::cout << "Cost=" << cost << " Div=" << div << " avg=" << avg << std::endl;
    std::cout << "Cost=" << cost << std::endl;
    std::cout << "Press any key" << std::endl;
    char c; std::cin >> c;
#endif
    //return div+prm->penalization*fabs(avg);
    //return cost+prm->penalization*fabs(avg);
    //return div;
    return cost;
}