#include <tomo_align_dual_tilt_series.h>

Inheritance diagram for ProgAlignDualTiltSeries:

Collaboration diagram for ProgAlignDualTiltSeries:

Public Member Functions
void	readParams ()
	Read parameters from command line. More...

void	show ()
	Show parameters. More...

void	defineParams ()
	Usage. More...

void	readDual ()
	Read dual series. More...

void	produceSideInfo ()
	Produce side info. More...

void	findParametersAt0degrees (bool rotateDual)
	Find parameters (shift+rotation) at 0 degrees. More...

double	distanceBetweenCommonLines (int refi, int dualj, const Matrix2D< double > &E, double X, double Y, double Z)
	Distance between a pair of common lines. More...

double	evaluateAlignment (const Matrix1D< double > &_alignment)
	Evaluate alignment. More...

double	optimizeAlignment ()
	Optimize current alignment. More...

void	alignDual ()
	Align dual. More...

void	run ()
	Run. More...

void	shiftProjectionInZ (MultidimArray< double > &I, int dualj, double Z) const
	Shift the projection in Z. More...

Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)

const char *	getParam (const char param, const char subparam, int arg=0)

int	getIntParam (const char *param, int arg=0)

int	getIntParam (const char param, const char subparam, int arg=0)

double	getDoubleParam (const char *param, int arg=0)

double	getDoubleParam (const char param, const char subparam, int arg=0)

float	getFloatParam (const char *param, int arg=0)

float	getFloatParam (const char param, const char subparam, int arg=0)

void	getListParam (const char *param, StringVector &list)

int	getCountParam (const char *param)

bool	checkParam (const char *param)

bool	existsParam (const char *param)

void	addParamsLine (const String &line)

void	addParamsLine (const char *line)

ParamDef *	getParamDef (const char *param) const

virtual void	quit (int exit_code=0) const

virtual int	tryRun ()

void	initProgress (size_t total, size_t stepBin=60)

void	setProgress (size_t value=0)

void	endProgress ()

void	processDefaultComment (const char param, const char left)

void	setDefaultComment (const char param, const char comment)

virtual void	initComments ()

void	setProgramName (const char *name)

void	addUsageLine (const char *line, bool verbatim=false)

void	clearUsage ()

void	addExampleLine (const char *example, bool verbatim=true)

void	addSeeAlsoLine (const char *seeAlso)

void	addKeywords (const char *keywords)

const char *	name () const

virtual void	usage (int verb=0) const

virtual void	usage (const String &param, int verb=2)

int	version () const

virtual void	show () const

virtual void	read (int argc, const char **argv, bool reportErrors=true)

virtual void	read (int argc, char **argv, bool reportErrors=true)

void	read (const String &argumentsLine)

	XmippProgram ()

	XmippProgram (int argc, const char **argv)

virtual	~XmippProgram ()

Public Attributes
FileName	fnRef
	Reference tilt series. More...

FileName	fnDual
	Dual tilt series. More...

FileName	fnOut
	Aligned tilt series. More...

double	scaleFactor
	Scale factor. More...

bool	debugging

MetaDataVec	SFRef

MetaDataVec	SFDual

MultidimArray< double >	tiltRef

MultidimArray< double >	tiltDual

std::vector< MultidimArray< double > >	imgRef

std::vector< MultidimArray< double > >	imgDual

bool	rotatedDual

FileName	fnRef0

FileName	fnDual0

Matrix1D< double >	alignment

Matrix1D< double >	normali

Matrix1D< double >	normalj

Matrix1D< double >	commonline

Matrix1D< double >	commonlinejE

Matrix1D< double >	commonlinei

Matrix1D< double >	commonlinej

MultidimArray< double >	profilei

MultidimArray< double >	profilej

MultidimArray< double >	I

Matrix2D< double >	E

Matrix2D< double >	Et

Matrix2D< double >	A

int	Xdim

int	Ydim

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Additional Inherited Members
Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Attributes inherited from XmippProgram
int	errorCode

ProgramDef *	progDef
	Program definition and arguments parser. More...

std::map< String, CommentList >	defaultComments

int	argc
	Original command line arguments. More...

const char **	argv

Detailed Description

Align dual tilt series parameters.

Definition at line 38 of file tomo_align_dual_tilt_series.h.

Member Function Documentation

◆ alignDual()

void ProgAlignDualTiltSeries::alignDual ( )

Align dual.

Definition at line 370 of file tomo_align_dual_tilt_series.cpp.

 {
 /* There might be a problem, the code seemingly working in TomoJ is
 
         for(int i=0;i<ts.length;i++){ // ts is an array of tilt series
             double[] ali=ts[i].getGlobalOrientation(); // I get the
                 alignment computed between the tilt series with common lines
             if(ali!=null){ // is there an alignment?
                 double cxtmp=ali[3]; //this is the shift of the center of
                                         volume original is (0,0,0)
                 double cytmp=ali[4];
                 //Now where will be this center of tilt axis in the volume?
                 DoubleMatrix2D et= TiltSeries.eulerAngles2Matrix(ali[0],
                     ali[1], ali[2]).viewDice(); // I compute the euler angles of the rotation
                                                     between the tilt series
                 cx[i]=cxtmp*et.getQuick(0,0)+cytmp*et.getQuick(0,1)+width/2;
                 // I compute the new position of the center of the volume + add width/2 for
                     correct indexation of arrays (in Xmipp should not be there)
  
                 cy[i]=cxtmp*et.getQuick(1,0)+cytmp*et.getQuick(1,1)+height/2;
 
 
             }else{ //there is no alignment so I define it as classical
                 cx[i]=width/2;
                 cy[i]=height/2;
             }
         }
 */
 
     Euler_angles2matrix(alignment(0),alignment(1),alignment(2),E);
     Matrix1D<double> shift3D=vectorR3(alignment(3),alignment(4),alignment(5));
     Matrix1D<double> shift2D=vectorR2(alignment(3),alignment(4));
     shift3D/=scaleFactor;
     shift2D/=scaleFactor;
     MetaDataVec SFout;
     int n=0;
     //std::cout << "Aligning dual" << std::endl;
     Image<double> Idual;
     Matrix2D<double> Edual;
     FileName fnImg, fn;
 
     for (size_t objId : SFDual.ids())
     {
         SFDual.getValue(MDL_IMAGE,fnImg,objId);
         Idual.read(fnImg);
         Idual().setXmippOrigin();
         if (rotatedDual)
             selfRotate(xmipp_transformation::BSPLINE3,Idual(),180,xmipp_transformation::DONT_WRAP);
         selfTranslate(xmipp_transformation::BSPLINE3,Idual(),shift2D,xmipp_transformation::DONT_WRAP);
         shiftProjectionInZ(Idual(), n, ZZ(shift3D));
         Euler_angles2matrix(0., tiltDual(n), 0., Edual);
         Edual=Edual*E;
         double rot, tilt, psi;
         Euler_matrix2angles(Edual, rot, tilt, psi);
         fn.compose(n+1,fnOut+".stk");
         Idual.write(fn);
         size_t id = SFout.addObject();
         SFout.setValue(MDL_IMAGE, fn, id);
         SFout.setValue(MDL_ANGLE_ROT, rot, id);
         SFout.setValue(MDL_ANGLE_TILT, tilt, id);
         SFout.setValue(MDL_ANGLE_PSI, psi, id);
         n++;
     }
     SFout.write(fnOut.removeDirectories()+".sel");
 }

◆ defineParams()

void ProgAlignDualTiltSeries::defineParams ( )

virtual

Usage.

Reimplemented from XmippProgram.

Definition at line 63 of file tomo_align_dual_tilt_series.cpp.

 {
     addUsageLine("Align two dual tilt series that have been previously internally aligned");
     addUsageLine("with xmipp_tomo_align_tilt_series");
     addSeeAlsoLine("tomo_align_tilt_series");
     addParamsLine("  --ref  <selfile>  : Reference tilt series");
     addParamsLine("                    : The selfile must contain the list of micrographs");
     addParamsLine("                    : and its tilt angles");
     addParamsLine("  --dual <selfile>  : Dual tilt series");
     addParamsLine(" [-o <rootname=\"\">]  : Rootname for the aligned tilt series");
     addParamsLine("                    : If not given, the dual tilt series+_aligned");
     addParamsLine(" [--scale <s=0.25>] : Scale for performing the common line comparisons");
 }

◆ distanceBetweenCommonLines()

double ProgAlignDualTiltSeries::distanceBetweenCommonLines	(	int	refi,
		int	dualj,
		const Matrix2D< double > &	E,
		double	X,
		double	Y,
		double	Z
	)

Distance between a pair of common lines.

Definition at line 218 of file tomo_align_dual_tilt_series.cpp.

 {
     SPEED_UP_temps012;
 
     // Compute the direction of the common line in the
     // universal coordinate system
     Euler_direction(0, tiltRef(refi), 0, normali);
     Euler_direction(0, tiltDual(dualj), 0, normalj);
     M3x3_BY_V3x1(normalj, Et, normalj);
 
     vectorProduct(normali,normalj,commonline);
     if (commonline.module()<XMIPP_EQUAL_ACCURACY)
         return -2;
     commonline.selfNormalize();
 
     // Compute the direction of the common line in each
     // of the projection coordinate systems
     Uproject_to_plane(commonline, 0, tiltRef(refi), 0, commonlinei);
     M3x3_BY_V3x1(commonlinejE, E, commonline);
     Uproject_to_plane(commonlinejE, 0, tiltDual(dualj), 0, commonlinej);
     commonlinei.selfNormalize();
     commonlinej.selfNormalize();
 
     // Compute the angle of the common line in i and j images
     double angi=-RAD2DEG(atan2(YY(commonlinei),XX(commonlinei)));
     double angj=-RAD2DEG(atan2(YY(commonlinej),XX(commonlinej)));
 
     I=imgRef[refi];
 #ifdef DEBUG
     debugging=true;
     Image<double> save;
     save()=I;
     save.write("PPPref.xmp");
 #endif
 
     selfRotate(xmipp_transformation::LINEAR,I, -angi, xmipp_transformation::DONT_WRAP);
     profilei.initZeros();
     FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
     A1D_ELEM(profilei,j)+=A2D_ELEM(I,i,j);
 #ifdef DEBUG
 
     save()=I;
     save.write("PPPrefaligned.xmp");
 #endif
 
     I=imgDual[dualj];
 #ifdef DEBUG
 
     save()=I;
     save.write("PPPdual.xmp");
 #endif
 
     // Modify Idual according to z
     shiftProjectionInZ(I, dualj, Z/scaleFactor);
 
     // Modify Idual according to x,y,angj
     // First translation, then rotation
     double cj=COSD(angj);
     double sj=SIND(angj);
     MAT_ELEM(A,0,0)=MAT_ELEM(A,1,1)=cj;
     MAT_ELEM(A,1,0)=sj;
     MAT_ELEM(A,0,1)=-sj;
     MAT_ELEM(A,0,2)=cj*X-sj*Y;
     MAT_ELEM(A,1,2)=sj*X+cj*Y;
     selfApplyGeometry(xmipp_transformation::LINEAR, I, A, xmipp_transformation::IS_NOT_INV, xmipp_transformation::DONT_WRAP);
     profilej.initZeros();
     FOR_ALL_ELEMENTS_IN_ARRAY2D(I)
     A1D_ELEM(profilej,j)+=A2D_ELEM(I,i,j);
 #ifdef DEBUG
 
     save()=I;
     save.write("PPPdualAligned.xmp");
 #endif
 
     if (debugging)
     {
         std::cout << "refi=" << refi << " dualj=" << dualj << std::endl
         << "   tiltRef=" << tiltRef(refi) << " tiltDual=" << tiltDual(dualj) << std::endl
         << "   normali=" << normali.transpose() << std::endl
         << "   normalj=" << normalj.transpose() << std::endl
         << "   commonline=" << commonline.transpose() << std::endl
         << "   commonlinei=" << commonlinei.transpose() << std::endl
         << "   commonlinejE=" << commonlinejE.transpose() << std::endl
         << "   commonlinej=" << commonlinej.transpose() << std::endl
         << "   angi=" << angi << " angj=" << angj << std::endl
         << "   A\n" << A << std::endl
         << "   corr=" << correlationIndex(profilei,profilej) << std::endl
         ;
 
         profilei.write("PPPi.txt");
         profilej.write("PPPj.txt");
 
         std::cout << "Press any key\n";
         char c;
         std::cin >> c;
     }
     return 1-correlationIndex(profilei,profilej);
 }

◆ evaluateAlignment()

double ProgAlignDualTiltSeries::evaluateAlignment ( const Matrix1D< double > & _alignment )

Evaluate alignment.

Definition at line 321 of file tomo_align_dual_tilt_series.cpp.

 {
     Euler_angles2matrix(_alignment(0),_alignment(1),_alignment(2),E);
     Et=E.transpose();
     double X=_alignment(3);
     double Y=_alignment(4);
     double Z=_alignment(5);
 
     int Nref=imgRef.size();
     int Ndual=imgDual.size();
     int Ncomparisons=0;
     double avgDistance=0;
     for (int nref=0; nref<Nref; nref++)
         for (int ndual=0; ndual<Ndual; ndual++)
         {
             double d=distanceBetweenCommonLines(nref,ndual,E,X,Y,Z);
             if (d>0)
             {
                 avgDistance+=d;
                 Ncomparisons++;
             }
         }
     if (Ncomparisons>0)
         avgDistance/=Ncomparisons;
     return avgDistance;
 }

◆ findParametersAt0degrees()

void ProgAlignDualTiltSeries::findParametersAt0degrees ( bool rotateDual )

Find parameters (shift+rotation) at 0 degrees.

Definition at line 158 of file tomo_align_dual_tilt_series.cpp.

 {
     Image<double> Iref, Idual;
     Iref.read(fnRef0);
     Idual.read(fnDual0);
     Iref().setXmippOrigin();
     Idual().setXmippOrigin();
     if (rotateDual)
         selfRotate(xmipp_transformation::BSPLINE3,Idual(),180);
     Matrix2D<double> M0,M0inv;
 #ifdef DEBUG
 
     Image<double> save;
     save()=Iref();
     save.write("PPPref0.xmp");
     save()=Idual();
     save.write("PPPdual0.xmp");
 #endif
 
     alignImages(Iref(), Idual(), M0);
     M0inv=M0.inv();
 #ifdef DEBUG
 
     save()=Idual();
     save.write("PPPdual0Aligned.xmp");
     std::cout << "M0\n" << M0 << std::endl;
     std::cout << "M0inv\n" << M0inv << std::endl;
     std::cout << "Press any key\n";
     char c;
     std::cin >> c;
 #endif
 
     alignment.resize(6);
     alignment(0)=-RAD2DEG(atan2(M0inv(1,0),M0inv(1,1)));  // rot
     alignment(1)=0;                                // tilt
     alignment(2)=0;                                // psi
     alignment(3)=-M0inv(0,2)*scaleFactor;          // X
     alignment(4)=-M0inv(1,2)*scaleFactor;          // Y
     alignment(5)=0;                                // Z
 
     /*
             alignment(0)=90;
             alignment(3)=-25*scaleFactor;          // X
             alignment(4)= 15*scaleFactor;          // Y
             alignment(5)=5;
     */
     if (rotateDual)
         std::cout << "First estimate (180) of (rot,tilt,psi,x,y,z)=\n"
         << alignment.transpose() << std::endl;
     else
         std::cout << "First estimate (0) of (rot,tilt,psi,x,y,z)=\n"
         << alignment.transpose() << std::endl;
     /*
             debugging=true;
     */
 }

◆ optimizeAlignment()

double ProgAlignDualTiltSeries::optimizeAlignment ( )

Optimize current alignment.

Definition at line 349 of file tomo_align_dual_tilt_series.cpp.

 {
     Matrix1D<double> steps;
     steps.resize(alignment);
     steps.initConstant(1);
     double fitness;
     int iter;
     powellOptimizer(alignment,1,6,&wrapperDualAligment,this,0.01,fitness,
                     iter,steps,true);
 #ifdef DEBUG
 
     debugging=true;
     evaluateAlignment(alignment);
     debugging=false;
 #endif
 
     return fitness;
 }

◆ produceSideInfo()

void ProgAlignDualTiltSeries::produceSideInfo ( )

Produce side info.

Definition at line 107 of file tomo_align_dual_tilt_series.cpp.

 {
     debugging=false;
     if (fnOut=="")
         fnOut=fnDual.withoutExtension()+"_aligned";
 
     std::cout << "Reading data...\n";
     SFRef.read(fnRef);
     SFDual.read(fnDual);
 
     // Read Reference series
     tiltRef.initZeros(SFRef.size());
     int i=0;
     Image<double> I;
     FileName fnImg;
     double minAbsTilt=1000;
     for (size_t objId: SFRef.ids())
     {
         SFRef.getValue(MDL_IMAGE,fnImg,objId);
         SFRef.getValue(MDL_ANGLE_TILT,tiltRef(i),objId);
         if (fabs(tiltRef(i))<minAbsTilt)
         {
             minAbsTilt=fabs(tiltRef(i));
             fnRef0=fnImg;
         }
         I.read(fnImg);
         selfScaleToSize(xmipp_transformation::BSPLINE3,I(),ROUND(YSIZE(I())*scaleFactor),
                         ROUND(XSIZE(I())*scaleFactor));
         I().setXmippOrigin();
         imgRef.push_back(I());
         i++;
     }
 
     // Read Dual series
     readDual();
 
     normali.resize(3);
     normalj.resize(3);
     commonline.resize(3);
     commonlinei.resize(3);
     commonlinej.resize(3);
     commonlinejE.resize(3);
     A.initIdentity(3);
     int minDim=XMIPP_MIN(Xdim,Ydim);
     profilei.resize(minDim);
     profilei.setXmippOrigin();
     profilej=profilei;
 }

◆ readDual()

void ProgAlignDualTiltSeries::readDual ( )

Read dual series.

Definition at line 78 of file tomo_align_dual_tilt_series.cpp.

 {
     imgDual.clear();
     tiltDual.initZeros(SFDual.size());
     int i=0;
     Image<double> I;
     FileName fnImg;
     double minAbsTilt=1000;
     for (size_t objId: SFDual.ids())
     {
         SFDual.getValue(MDL_IMAGE,fnImg, objId);
         SFDual.getValue(MDL_ANGLE_TILT,tiltDual(i), objId);
         if (fabs(tiltDual(i))<minAbsTilt)
         {
             minAbsTilt=fabs(tiltDual(i));
             fnDual0=fnImg;
         }
         I.read(fnImg);
         selfScaleToSize(xmipp_transformation::BSPLINE3,I(),ROUND(YSIZE(I())*scaleFactor),
                         ROUND(XSIZE(I())*scaleFactor));
         I().setXmippOrigin();
         Xdim=XSIZE(I());
         Ydim=YSIZE(I());
         imgDual.push_back(I());
         i++;
     }
 }

◆ readParams()

void ProgAlignDualTiltSeries::readParams ( )

virtual

Read parameters from command line.

Read parameters.

Reimplemented from XmippProgram.

Definition at line 43 of file tomo_align_dual_tilt_series.cpp.

 {
     fnRef=getParam("--ref");
     fnDual=getParam("--dual");
     fnOut=getParam("-o");
     scaleFactor=getDoubleParam("--scale");
 }

◆ run()

void ProgAlignDualTiltSeries::run ( )

virtual

Run.

Reimplemented from XmippProgram.

Definition at line 437 of file tomo_align_dual_tilt_series.cpp.

 {
     produceSideInfo();
 
     // Check if there is a rotation of 180 degrees
     findParametersAt0degrees(false);
     double objective0=evaluateAlignment(alignment);
     Matrix1D<double> alignmentBackup=alignment;
 
     int Ndual=imgDual.size();
     for (int ndual=0; ndual<Ndual; ndual++)
         selfRotate(xmipp_transformation::BSPLINE3,imgDual[ndual],180);
     findParametersAt0degrees(true);
     double objective180=evaluateAlignment(alignment);
     std::cout << "objective0=" << objective0 << " objective180=" << objective180 << std::endl;
     if (objective0<objective180)
     {
         readDual();
         rotatedDual=false;
         alignment=alignmentBackup;
     }
     else
         rotatedDual=true;
 
     optimizeAlignment();
     alignDual();
 }

◆ shiftProjectionInZ()

void ProgAlignDualTiltSeries::shiftProjectionInZ	(	MultidimArray< double > &	I,
		int	dualj,
		double	Z
	)		const

Shift the projection in Z.

Definition at line 465 of file tomo_align_dual_tilt_series.cpp.

 {
     if (Z==0)
         return;
     FourierTransformer transformer;
     MultidimArray< std::complex<double> > Ifft;
     Matrix1D<double> w(3);
     Matrix1D<int> idx(3);
     Matrix2D<double> Edual, Edualt;
     Euler_angles2matrix(0., tiltDual(dualj), 0., Edual);
     Edualt=Edual.transpose();
     transformer.FourierTransform(I,Ifft,false);
     FOR_ALL_ELEMENTS_IN_ARRAY2D(Ifft)
     {
         SPEED_UP_temps012;
         VECTOR_R2(idx,j,i);
         FFT_idx2digfreq(I,idx,w);
         ZZ(w)=0;
         M3x3_BY_V3x1(w,Edualt,w);
         double arg=ZZ(w)*Z;
         A2D_ELEM(Ifft,i,j)*=std::complex<double>(cos(arg),sin(arg));
     }
     transformer.inverseFourierTransform();
 }

◆ show()

void ProgAlignDualTiltSeries::show ( )

Show parameters.

Definition at line 52 of file tomo_align_dual_tilt_series.cpp.

 {
     std::cout
     << "Reference: " << fnRef       << std::endl
     << "Dual:      " << fnDual      << std::endl
     << "Output:    " << fnOut       << std::endl
     << "Scale:     " << scaleFactor << std::endl
     ;
 }