#include <image_assignment_tilt_pair.h>

Inheritance diagram for ProgassignmentTiltPair:

Collaboration diagram for ProgassignmentTiltPair:

Public Member Functions
void	defineParams ()

void	readParams ()

void	search_affine_transform (float u1x, float u1y, float u2x, float u2y, float u3x, float u3y, float t1x, float t1y, float t2x, float t2y, float t3x, float t3y, const Matrix1D< double > &ux, const Matrix1D< double > &uy, size_t Xdim, size_t Ydim, struct Delaunay_T &delaunay_tilt, int &bestInliers, Matrix2D< double > &A_coarse, Matrix1D< double > &T_coarse, bool contingency, int thrs)

void	findMaximumMinimum (const float u1, const float u2, const float u3, double &u_max, double &u_min)

bool	checkwindow (const float t1, const float t2, const float t3, const double u_max, const double u_min)

void	run ()

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	fntilt

FileName	fnuntilt

FileName	fndir

FileName	fnmic

int	mshift

double	particle_size

double	thr

double	tiltest

MetaDataVec	mdPartial

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

Definition at line 41 of file image_assignment_tilt_pair.h.

Member Function Documentation

◆ checkwindow()

bool ProgassignmentTiltPair::checkwindow	(	const float	t1,
		const float	t2,
		const float	t3,
		const double	u_max,
		const double	u_min
	)

Definition at line 209 of file image_assignment_tilt_pair.cpp.

 {
     double window_p, window_m;
     window_p = u_max + 1.2*mshift;
     window_m = u_min - 1.2*mshift;
 
 
     if (((t1<window_p) && (t2<window_p) && (t3<window_p)) && ((t1>window_m) && (t2>window_m) && (t3>window_m)))
         return true;
     else
         return false;
 
 }

◆ defineParams()

void ProgassignmentTiltPair::defineParams ( )

virtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 47 of file image_assignment_tilt_pair.cpp.

 {
     //usage
     addUsageLine("Validate a 3D reconstruction from its projections attending to directionality and spread of the angular assignments from a given significant value");
     //params
 
     addParamsLine("  [--untiltcoor <md_file=\"\">]    : Untilt coordinates");
     addParamsLine("  [--tiltcoor <md_file=\"\">]    : Tilt coordinates");
     addParamsLine("  [--tiltmicsize <img_file=\"\">]    : Tilt micrography");
     addParamsLine("  [--odir <outputDir=\".\">]   : Output directory");
     addParamsLine("  [--maxshift <s=1000>]   : Maximum shift");
     addParamsLine("  [--tiltangle <s=-1>]   : Tilt angle estimation, the method will look for the assignment in the interval of [tiltangle-15º, til_angle+15º]");
     addParamsLine("  [--particlesize <p=100>]   : Particle size");
     addParamsLine("  [--threshold <d=0.3>]      : If the distance between two points is lesser than threshold*particlesize, they will be the same point");
 }

◆ findMaximumMinimum()

void ProgassignmentTiltPair::findMaximumMinimum	(	const float	u1,
		const float	u2,
		const float	u3,
		double &	u_max,
		double &	u_min
	)

Definition at line 181 of file image_assignment_tilt_pair.cpp.

 {
 if (u1 > u2)
 {
     if (u1 > u3)
         u_max = u1;
     else
         u_max = u3;
 
     if (u2 > u3)
         u_min = u3;
     else
         u_min = u2;
 }
 else
 {
     if (u2 > u3)
         u_max = u2;
     else
         u_max = u3;
 
     if (u1 > u3)
         u_min = u3;
     else
         u_min = u1;
 }
 }

◆ readParams()

void ProgassignmentTiltPair::readParams ( )

virtual

Function in which each program will read parameters that it need. If some error occurs the usage will be printed out.

Reimplemented from XmippProgram.

Definition at line 35 of file image_assignment_tilt_pair.cpp.

 {
     fnuntilt = getParam("--untiltcoor");
     fntilt = getParam("--tiltcoor");
     fnmic = getParam("--tiltmicsize");
     fndir = getParam("--odir");
     mshift = getDoubleParam("--maxshift");
     tiltest = getDoubleParam("--tiltangle");
     particle_size = getDoubleParam("--particlesize");
     thr = getDoubleParam("--threshold");
 }

◆ run()

void ProgassignmentTiltPair::run ( )

virtual

This function will be start running the program. it also should be implemented by derived classes.

Reimplemented from XmippProgram.

Definition at line 225 of file image_assignment_tilt_pair.cpp.

 {
     std::cout << "Starting..." << std::endl;
 
     //LOAD METADATA and TRIANGULATIONS
     MetaDataVec md_untilt, md_tilt, mduntilt, mdtilt;
     size_t Ndim, Zdim, Ydim , Xdim, objId;
 
     getImageSize(fnmic, Xdim, Ydim, Zdim, Ndim);
 
     bool exist1, exist2;
     bool flag_assign = true;
 
     exist1 = fnuntilt.existsTrim();
     exist2 = fntilt.existsTrim();
 
     if ((exist1 && exist2) == true)
     {
     md_untilt.read(fnuntilt);
     md_tilt.read(fntilt);
 
     int x, y, len_u=0, len_t=0;
 
     //storing untilted points and creating Delaunay triangulation
     struct Delaunay_T delaunay_untilt;
     Matrix1D<double> ux(md_untilt.size()), uy(md_untilt.size()), tx(md_tilt.size()), ty(md_tilt.size());
     init_Delaunay( &delaunay_untilt, md_untilt.size());
     for (size_t objId : md_untilt.ids())
     {
         md_untilt.getValue(MDL_XCOOR, x, objId);
         md_untilt.getValue(MDL_YCOOR, y, objId);
 
         VEC_ELEM(ux,len_u) = x;
         VEC_ELEM(uy,len_u) = y;
 
         insert_Point( &delaunay_untilt, x, y);
         len_u = len_u +1;
     }
     create_Delaunay_Triangulation( &delaunay_untilt, 0);
     int tri_number_untilt = get_Number_Real_Faces(delaunay_untilt.dcel);
 
     //storing tilted points and creating Delaunay triangulation
     struct Delaunay_T delaunay_tilt;
     init_Delaunay( &delaunay_tilt, md_tilt.size());
     for (size_t objId : md_tilt.ids())
     {
         md_tilt.getValue(MDL_XCOOR, x, objId);
         md_tilt.getValue(MDL_YCOOR, y, objId);
 
         VEC_ELEM(tx,len_t) = x;
         VEC_ELEM(ty,len_t) = y;
 
         insert_Point( &delaunay_tilt, x, y);
         len_t = len_t +1;
     }
     int thrs;
     if (len_u<len_t)
     {
          thrs = len_u;
     }
     else
     {
         thrs = len_t;
     }
 
     if (len_u == 0 || len_t == 0)
     {
         flag_assign = false;
     }
 
     if (flag_assign == true)
     {
     create_Delaunay_Triangulation( &delaunay_tilt, 1);
     int tri_number_tilt = get_Number_Real_Faces(delaunay_tilt.dcel);
 
     //DETERMINING TRIANGLES AREAS
     struct Point_T p, q, r, u1, u2, u3, t1, t2, t3;
     MultidimArray<double> trig_untilt_area(tri_number_untilt+1), trig_tilt_area(tri_number_tilt+1), sortedArea;
     Matrix1D<double> trig_untilt_area_1D(tri_number_untilt+1), trig_tilt_area_1D(tri_number_tilt+1);
 
     for (int i=1; i<tri_number_untilt+1 ;i++)
     {
         get_Face_Points(&delaunay_untilt, i, &p, &q, &r); //i is the triangle number
         A1D_ELEM(trig_untilt_area,i-1) = triangle_area(p.x, p.y, q.x, q.y, r.x, r.y);
     }
 
     for (int i=1; i<tri_number_tilt+1 ;i++)
     {
         get_Face_Points(&delaunay_tilt, i, &p, &q, &r);
         A1D_ELEM(trig_tilt_area,i-1) = triangle_area(p.x, p.y, q.x, q.y, r.x, r.y);
     }
 
     trig_untilt_area.sort(sortedArea);
 
     double threshold_area = sortedArea( round((tri_number_untilt+1)*0.2) );
 
     //DEFINING PARAMETERS FOR SEARCH_AFFINE_TRANSFORM
     int bestInliers=0, bestInliers_con=0;
     Matrix2D<double> A_coarse;
     Matrix1D<double> T_coarse, def_T, t_test(2), u(2), dist_vec, dist_vec2;
 
     struct Point_T t_dist, t_closest;
 
     //TILT ANGLE ESTIMATION
     double tilt_max = (tiltest + 15)*PI/180, tilt_min = (tiltest - 15)*PI/180, cos_tilt_max, cos_tilt_min;
 
     if (tiltest == -1)
     {
         cos_tilt_max = 0;
         cos_tilt_min = 1;
     }
     else
     {
         cos_tilt_max= cos(tilt_max);
         cos_tilt_min= cos(tilt_min);
     }
 
 
     if (verbose==1)
     {
         std::cerr << "Exploring triangle matching" << std::endl;
         init_progress_bar(tri_number_untilt);
     }
 
 
     if (verbose==1)
         std::cerr << "Coarse Phase" << std::endl;
 
     for (int k=0; k<tri_number_untilt; k++)
     {
         if (trig_untilt_area(k) < threshold_area)
             continue;
 
         for (int j=0; j<tri_number_tilt; j++)
         {
             //std::cout << "Iteration k = " << k << "     Iteration j = " << j << std::endl;
             if (trig_untilt_area(k) < trig_tilt_area(j))
                 continue;
 
             if ( (trig_untilt_area(k)*cos_tilt_min < trig_tilt_area(j)) || (trig_untilt_area(k)*cos_tilt_max > trig_tilt_area(j)) )
                 continue;
 
             get_Face_Points(&delaunay_untilt, k, &u1, &u2, &u3);
             get_Face_Points(&delaunay_tilt, j, &t1, &t2, &t3);
 
 
             //New Idea
             double ux_max, uy_max, ux_min, uy_min;
             findMaximumMinimum(u1.x, u2.x, u3.x, ux_max, ux_min);
             findMaximumMinimum(u1.y, u2.y, u3.y, uy_max, uy_min);
 
             if (!checkwindow(t1.x, t2.x, t3.x, ux_max, ux_min))
                 continue;
             if (!checkwindow(t1.y, t2.y, t3.y, uy_max, uy_min))
                 continue;
 
             search_affine_transform(u1.x, u1.y, u2.x, u2.y, u3.x, u3.y, t1.x,
                     t1.y, t2.x, t2.y, t3.x, t3.y,
                     ux, uy, Xdim, Ydim, delaunay_tilt, bestInliers,
                     A_coarse, T_coarse, false, thrs);
         }
         if (verbose==1 && k%100==0)
             progress_bar(k);
     }
     std::cout << "Coarse Inliers = " << bestInliers << std::endl;
 
 
         if (verbose==1)
         std::cerr << "Refinement Phase" << std::endl;
 
         dist_vec.initConstant(len_u, -1);
         Matrix1D<double> t_clo(2*bestInliers), X;
         Matrix2D<double> invref, invW2;
         //PseudoInverseHelper pseudoInverter;
         invref.initZeros(2*bestInliers,6);
         invW2.initZeros(2*bestInliers,6);
         int count = 0;
 
         PseudoInverseHelper h;
         for (int k=0; k<len_u; k++)
         {
             VEC_ELEM(u,0) = VEC_ELEM(ux,k);
             VEC_ELEM(u,1) = VEC_ELEM(uy,k);
             t_test = A_coarse*u + T_coarse;
             t_dist.x = VEC_ELEM(t_test,0);
             t_dist.y = VEC_ELEM(t_test,1);
             double dist;
             if (!select_Closest_Point(&delaunay_tilt, &t_dist, &t_closest, &dist) )
             {
                 //std::cerr << "WARNING IN TRIANGULATION OR CLOSEST NEIGHBOUR (in Coarse Phase)" << std::endl;
                 //std::cout << "Maybe the warning involves the tilt coordinates " << "(" << t_dist.x << ", " << t_dist.y << ")" << std::endl;
                 continue;
             }
             if (dist<thr*particle_size)
             {
                 VEC_ELEM(t_clo,count) = t_closest.x;
                 VEC_ELEM(t_clo,count+bestInliers) = t_closest.y;
                 MAT_ELEM(invW2,count,0) = VEC_ELEM(u,0);
                 MAT_ELEM(invW2,count,1) = VEC_ELEM(u,1);
                 MAT_ELEM(invW2,count,4) = 1;
 
                 MAT_ELEM(invW2, count+bestInliers, 2) = VEC_ELEM(u,0);
                 MAT_ELEM(invW2, count+bestInliers, 3) = VEC_ELEM(u,1);
                 MAT_ELEM(invW2, count+bestInliers, 5) = 1;
 
                 count = count + 1;
             }
 
         }
 
 
         Matrix2D<double> def_A, A_con;
         Matrix1D<double> T_con;
         T_con.initZeros(2);
         A_con.initZeros(2,2);
         def_A.initZeros(2,2);
         def_T.initZeros(2);
         bool flag = false;
 
 
         if ((count >= bestInliers) && (count >= 0.2*thrs))
         {
             //std::cout << "Count > bestInliers" << std::endl;
             h.A=invW2;
 
             h.b=t_clo;
             solveLinearSystem(h,X);
             MAT_ELEM(def_A,0,0) = VEC_ELEM(X,0);
             MAT_ELEM(def_A,0,1) = VEC_ELEM(X,1);
             MAT_ELEM(def_A,1,0) = VEC_ELEM(X,2);
             MAT_ELEM(def_A,1,1) = VEC_ELEM(X,3);
             VEC_ELEM(def_T,0) = VEC_ELEM(X,4);
             VEC_ELEM(def_T,1) = VEC_ELEM(X,5);
             std::cout << "Best fitting is gotten using refinement phase" << std::endl;
             std::cout << "Refinement Inliers = " << count << std::endl;
             std::cout << "Coarse Inliers = " << bestInliers << std::endl;
             std::cout << "A" << def_A << std::endl;
             std::cout << "---------------------------" << std::endl;
             std::cout << "T" << def_T << std::endl;
             flag = true;
         }
         if (((count <= bestInliers) && (bestInliers >= 0.2*thrs)) && !flag)
         {
             std::cout << "Best fitting is gotten using coarse phase" << std::endl;
             std::cout << "Coarse Inliers = " << bestInliers << std::endl;
             std::cout << "Refinement Inliers = " << count << std::endl;
             def_A = A_coarse;
             def_T = T_coarse;
             std::cout << "A" << def_A << std::endl;
             std::cout << "---------------------------" << std::endl;
             std::cout << "T" << def_T << std::endl;
             flag = true;
         }
         if (((count <=bestInliers || count >=bestInliers) && ((count < 0.2*thrs) || (bestInliers < 0.2*thrs))) && !flag)
         {
             std::cerr << " Matching not found" << std::endl;
             std::cerr << " Starting contingency method" << std::endl;
             std::cerr << " This phase can take a long time" << std::endl;
             //Defining triangles in untilted space
             double window=4*particle_size;
             double window_t=window *(1-0.34);  //0.34 (approx 0.3) is cos(70), tilts higher than 70 degrees are not considered
             std::vector<Matrix1D<double> > allTrianglesU;
             Matrix1D<double> triangleu(7);
             for (int k = 0; k< len_u; k++)
             {
                 VEC_ELEM(triangleu,0)=VEC_ELEM(ux,k);
                 VEC_ELEM(triangleu,1)=VEC_ELEM(uy,k);
                 for (int j=k+1; j<len_u; j++)
                 {
                     if ( ( (VEC_ELEM(ux,j)>VEC_ELEM(ux,k)+window) || (VEC_ELEM(ux,j)<VEC_ELEM(ux,k)-window)) ||
                          ( (VEC_ELEM(uy,j)>VEC_ELEM(uy,k)+window) || (VEC_ELEM(uy,j)<VEC_ELEM(uy,k)-window)) )
                         continue;
                     VEC_ELEM(triangleu,2)=VEC_ELEM(ux,j);
                     VEC_ELEM(triangleu,3)=VEC_ELEM(uy,j);
                     for (int l=j+1; l<len_u; l++)
                     {
                         if ( ( (VEC_ELEM(ux,l)>VEC_ELEM(ux,k)+window) || (VEC_ELEM(ux,l)<VEC_ELEM(ux,k)-window)) ||
                              ( (VEC_ELEM(uy,l)>VEC_ELEM(uy,k)+window) || (VEC_ELEM(uy,l)<VEC_ELEM(uy,k)-window)) )
                             continue;
 
                         VEC_ELEM(triangleu,6) = triangle_area(VEC_ELEM(ux,k), VEC_ELEM(uy,k), VEC_ELEM(ux,j), VEC_ELEM(uy,j), VEC_ELEM(ux,l), VEC_ELEM(uy,l));
                         if (VEC_ELEM(triangleu,6) < threshold_area)
                             continue;
 
                         VEC_ELEM(triangleu,4) = VEC_ELEM(ux,l);
                         VEC_ELEM(triangleu,5) = VEC_ELEM(uy,l);
                         allTrianglesU.push_back(triangleu);
                     }
                 }
             }
             //Defining triangles in tilted space
             std::vector<Matrix1D<double> > allTrianglesT;
             Matrix1D<double> trianglet(7);
             for (int k = 0; k< len_t; k++)
             {
                 VEC_ELEM(trianglet,0)=VEC_ELEM(tx,k);
                 VEC_ELEM(trianglet,1)=VEC_ELEM(ty,k);
                 for (int j=k+1; j<len_t; j++)
                 {
                     if ( ( (VEC_ELEM(tx,j)>VEC_ELEM(tx,k)+window_t) || (VEC_ELEM(tx,j)<VEC_ELEM(tx,k)-window_t)) ||
                          ( (VEC_ELEM(ty,j)>VEC_ELEM(ty,k)+window_t) || (VEC_ELEM(ty,j)<VEC_ELEM(ty,k)-window_t)) )
                         continue;
                     VEC_ELEM(trianglet,2)=VEC_ELEM(tx,j);
                     VEC_ELEM(trianglet,3)=VEC_ELEM(ty,j);
                     for (int l=j+1; l<len_t; l++)
                     {
                         if ( ( (VEC_ELEM(tx,l)>VEC_ELEM(tx,k)+window_t) || (VEC_ELEM(tx,l)<VEC_ELEM(tx,k)-window_t)) ||
                              ( (VEC_ELEM(ty,l)>VEC_ELEM(ty,k)+window_t) || (VEC_ELEM(ty,l)<VEC_ELEM(ty,k)-window_t)) )
                             continue;
                         VEC_ELEM(trianglet,6) = triangle_area(VEC_ELEM(tx,k), VEC_ELEM(ty,k), VEC_ELEM(tx,j), VEC_ELEM(ty,j), VEC_ELEM(tx,l), VEC_ELEM(ty,l));
                         VEC_ELEM(trianglet,4) = VEC_ELEM(tx,l);
                         VEC_ELEM(trianglet,5) = VEC_ELEM(ty,l);
                         allTrianglesT.push_back(trianglet);
                     }
                 }
             }
 
             //Searching the affine application
 
             for (size_t ku=0; ku<allTrianglesU.size(); ku++)
             {
 
                 const Matrix1D<double> &triangleU=allTrianglesU[ku];
 
                 for (size_t kt=0; kt<allTrianglesT.size(); kt++)
                 {
                     const Matrix1D<double> &triangleT=allTrianglesT[kt];
 
                     if (VEC_ELEM(triangleU,6)<VEC_ELEM(triangleT,6))
                         continue;
 
                     search_affine_transform(VEC_ELEM(triangleU,0), VEC_ELEM(triangleU,1), VEC_ELEM(triangleU,2), VEC_ELEM(triangleU,3),
                                             VEC_ELEM(triangleU,4), VEC_ELEM(triangleU,5), VEC_ELEM(triangleT,0), VEC_ELEM(triangleT,1),
                                             VEC_ELEM(triangleT,2), VEC_ELEM(triangleT,3), VEC_ELEM(triangleT,4), VEC_ELEM(triangleT,5),
                                             ux, uy, Xdim, Ydim, delaunay_tilt, bestInliers_con,
                                             A_con, T_con, true, thrs);
                 }
             }
             std::cout << "Contingency Inliers = " << bestInliers_con << std::endl;
             if ((bestInliers_con > bestInliers) && (bestInliers_con>count))
             {
                 std::cout << "Best fitting is gotten using contingency phase" << std::endl;
                 std::cout << "Contingency Inliers = " << bestInliers_con << std::endl;
                 std::cout << "Refinement Inliers = " << count << std::endl;
                 def_A = A_con;
                 def_T = T_con;
                 std::cout << "A" << def_A << std::endl;
                 std::cout << "---------------------------" << std::endl;
                 std::cout << "T" << def_T << std::endl;
             }
 //          else
 //          {
 //              def_A = A_coarse;
 //              def_T = T_coarse;
 //          }
         }
 
 
         for (int k=0; k<len_u; k++)
         {
             VEC_ELEM(u,0) = VEC_ELEM(ux,k);
             VEC_ELEM(u,1) = VEC_ELEM(uy,k);
             t_test = def_A*u + def_T;
             t_dist.x = VEC_ELEM(t_test,0);
             t_dist.y = VEC_ELEM(t_test,1);
             double dist;
             if (!select_Closest_Point(&delaunay_tilt, &t_dist, &t_closest, &dist) )
             {
                 //std::cerr << "WARNING IN TRIANGULATION OR CLOSEST NEIGHBOUR (in writing phase)" << std::endl;
                 //std::cout << "Maybe the warning involves the tilt coordinates " << "(" << t_dist.x << ", " << t_dist.y << ")" << std::endl;
                 continue;
             }
 
             if (dist<thr*particle_size)
             {
                 objId = mduntilt.addObject();
                 mduntilt.setValue(MDL_XCOOR,(int) VEC_ELEM(u,0),objId);
                 mduntilt.setValue(MDL_YCOOR,(int) VEC_ELEM(u,1),objId);
 
                 objId = mdtilt.addObject();
                 mdtilt.setValue(MDL_XCOOR,(int) t_closest.x,objId);
                 mdtilt.setValue(MDL_YCOOR,(int) t_closest.y,objId);
             }
         }
         delete_Delaunay(&delaunay_untilt);
         delete_Delaunay(&delaunay_tilt);
     }
     }
     else
     {
         std::cout << "WARNING: input pos. file does not exist, therefore output files are empty" << std::endl;
     }
     if (flag_assign == false)
     {
         std::cout << "WARNING: input pos file is empty." << std::endl;
     }
 
     mduntilt.write((String)"particles@"+fndir+'/'+fnuntilt.getBaseName() + ".pos" );
     mdtilt.write((String)"particles@"+fndir+'/'+fntilt.getBaseName() + ".pos" );
 
 }

◆ search_affine_transform()

void ProgassignmentTiltPair::search_affine_transform	(	float	u1x,
		float	u1y,
		float	u2x,
		float	u2y,
		float	u3x,
		float	u3y,
		float	t1x,
		float	t1y,
		float	t2x,
		float	t2y,
		float	t3x,
		float	t3y,
		const Matrix1D< double > &	ux,
		const Matrix1D< double > &	uy,
		size_t	Xdim,
		size_t	Ydim,
		struct Delaunay_T &	delaunay_tilt,
		int &	bestInliers,
		Matrix2D< double > &	A_coarse,
		Matrix1D< double > &	T_coarse,
		bool	contingency,
		int	thrs
	)

Definition at line 63 of file image_assignment_tilt_pair.cpp.

 {
     double estimator, dist;
     struct Point_T t_dist, t_closest;
     Matrix1D<double> T;
     Matrix2D<double> invW;
     Matrix2D<double> A_matrix;
 
     A_matrix.initZeros(2,2);
 
     for (int i=0; i<6; i++)
     {
         if (i==0){
             def_affinity(u1x, u1y, u2x, u2y, u3x, u3y, t1x, t1y, t2x, t2y, t3x, t3y, A_matrix, T, invW);}
         if (i==1)
             def_affinity(u1x, u1y, u2x, u2y, u3x, u3y, t1x, t1y, t3x, t3y, t2x, t2y, A_matrix, T, invW);
         if (i==2)
             def_affinity(u1x, u1y, u2x, u2y, u3x, u3y, t2x, t2y, t1x, t1y, t3x, t3y, A_matrix, T, invW);
         if (i==3)
             def_affinity(u1x, u1y, u2x, u2y, u3x, u3y, t2x, t2y, t3x, t3y, t1x, t1y, A_matrix, T, invW);
         if (i==4)
             def_affinity(u1x, u1y, u2x, u2y, u3x, u3y, t3x, t3y, t1x, t1y, t2x, t2y, A_matrix, T, invW);
         if (i==5)
             def_affinity(u1x, u1y, u2x, u2y, u3x, u3y, t3x, t3y, t2x, t2y, t1x, t1y, A_matrix, T, invW);
 
         double trace_A = MAT_ELEM(A_matrix, 0, 0) + MAT_ELEM(A_matrix, 1, 1);
         double det_A = MAT_ELEM(A_matrix, 0, 0)*MAT_ELEM(A_matrix, 1, 1) - MAT_ELEM(A_matrix, 0, 1)*MAT_ELEM(A_matrix, 1, 0);
 
 
 //      if ( fabs(det_A - cos_tilt_max)>0.1)
 //          continue;
 
         double discriminant_A = 0.25*trace_A*trace_A - det_A;
 
         if (discriminant_A < DBL_EPSILON)
             continue;
 
         double sqrt_aux = sqrt(discriminant_A);
         double Eig_A1 = trace_A/2 + sqrt_aux;
         double Eig_A2 = trace_A/2 - sqrt_aux;
 
         if (Eig_A1<0.34 || Eig_A2<0.34 || fabs(Eig_A1-1)>0.05)
             continue;
 
         Matrix1D<double> u(2), dist_vec, dist_vec2, t_test(2);
         size_t len_u=VEC_XSIZE(ux);
 
         dist_vec.initConstant(len_u, -1);
         dist_vec2.initConstant(len_u, 10);
 
         for (size_t tt=0; tt<len_u; tt++)
         {
             VEC_ELEM(u,0) = VEC_ELEM(ux,tt);
             VEC_ELEM(u,1) = VEC_ELEM(uy,tt);
             t_test = A_matrix*u + T;
 
 
             if (VEC_ELEM(t_test,0)<0 || VEC_ELEM(t_test,0)>Xdim || VEC_ELEM(t_test,1)<0 || VEC_ELEM(t_test,1)>Ydim)
                 continue;
 
             t_dist.x = VEC_ELEM(t_test,0);
             t_dist.y = VEC_ELEM(t_test,1);
 
             if (!select_Closest_Point(&delaunay_tilt, &t_dist, &t_closest, &dist) )
             {
 
                 write_DCEL(delaunay_tilt.dcel, 0, "tiltdata_dcel.txt");
                 //std::cerr << "WARNING IN TRIANGULATION OR CLOSEST NEIGHBOUR" << std::endl;
                 //printf("Maybe the warning involves the tilt coordinates ( %f , %f ) \n", t_dist.x, t_dist.y);
                 continue;
             }
             VEC_ELEM(dist_vec,tt) = dist;
         }
 
         int counter = 0;
         size_t inliers2 = 0;
         double dist2 = 0;
 
         for (size_t kk=0; kk<len_u; kk++)
         {
             if (VEC_ELEM(dist_vec,kk)>-1)
             {
                 counter = counter + 1;
                 if (VEC_ELEM(dist_vec,kk) < thr*particle_size)
                 {
                     dist2 = VEC_ELEM(dist_vec,kk) + dist2;
                     inliers2 = inliers2 + 1;
                 }
             }
         }
 
         if (counter == 0)
             continue;
 
         if (inliers2 > bestInliers)
         {
             estimator = dist2/inliers2;
             bestInliers = inliers2;
             if (contingency)
                 if (bestInliers>0.3*thrs)
                     break;
 
             A_coarse = A_matrix;
             T_coarse = T;
             std::cout << bestInliers << std::endl;
         }
         else if ((inliers2 == bestInliers) && (dist2/inliers2 < estimator) )
         {
             estimator = dist2/inliers2;
             A_coarse = A_matrix;
             T_coarse = T;
         }
     }
 }