40 #define LEAST_SQUARES 2 64 double greyScale = p[0];
65 double greyShift = p[1];
81 if (greyScale!=1 || greyShift!=0)
119 double rot0, rotF, tilt0,
tiltF, psi0, psiF;
129 bool usePowell, onlyShift, useFRM, copyGeo, copyGray, store,
wrap;
139 addUsageLine(
"Align two volumes varying orientation, position and scale");
140 addParamsLine(
" --i1 <volume1> : the first volume to align");
141 addParamsLine(
" --i2 <volume2> : the second one");
142 addParamsLine(
" [--rot <rot0=0> <rotF=0> <step_rot=1>] : in degrees");
143 addParamsLine(
" [--tilt <tilt0=0> <tiltF=0> <step_tilt=1>] : in degrees");
144 addParamsLine(
" [--psi <psi0=0> <psiF=0> <step_psi=1>] : in degrees");
145 addParamsLine(
" [--scale <sc0=1> <scF=1> <step_sc=1>] : size scale margin");
146 addParamsLine(
" [--grey_scale <sc0=1> <scF=1> <step_sc=1>] : grey scale margin");
147 addParamsLine(
" requires --least_squares;");
148 addParamsLine(
" [--grey_shift <sh0=0> <shF=0> <step_sh=1>] : grey shift margin");
149 addParamsLine(
" requires --least_squares;");
150 addParamsLine(
" [-z <z0=0> <zF=0> <step_z=1>] : Z position in pixels");
151 addParamsLine(
" [-y <y0=0> <yF=0> <step_y=1>] : Y position in pixels");
152 addParamsLine(
" [-x <x0=0> <xF=0> <step_x=1>] : X position in pixels");
153 addParamsLine(
" [--show_fit] : Show fitness values");
154 addParamsLine(
" [--apply <file=\"\">] : Apply best movement to --i2 and store results in this file");
155 addParamsLine(
" [--covariance] : Covariance fitness criterion");
156 addParamsLine(
" [--least_squares] : LS fitness criterion");
157 addParamsLine(
" [--local] : Use local optimizer instead of exhaustive search");
158 addParamsLine(
" [--frm <maxFreq=0.25> <maxShift=10> <tilt0=-90> <tiltF=90>] : Use Fast Rotational Matching, if tilt0 and tiltF are left as -90 and 90 they do not effect the results of the alignment");
159 addParamsLine(
" : Maximum frequency is in digital frequencies (<0.5)");
160 addParamsLine(
" : Maximum shift is in pixels");
161 addParamsLine(
" :+ See Y. Chen, et al. Fast and accurate reference-free alignment of subtomograms. JSB, 182: 235-245 (2013)");
162 addParamsLine(
" [--onlyShift] : Only shift");
163 addParamsLine(
" [--dontScale] : Do not look for scale changes");
164 addParamsLine(
" [--copyGeo <file=\"\">] : copy transformation matrix in a txt file. ('A' matrix elements)");
165 addParamsLine(
" [--copyGray <file=\"\">] : copy gray scale and shift txt file.)");
166 addParamsLine(
" [--store <file=\"\">] : copy angles and shifts to a txt file.");
167 addParamsLine(
" [--dontWrap] : Do not wrap input2 when aligning to input1");
168 addParamsLine(
" == Mask Options == ");
170 addExampleLine(
"Typically you first look for a rough approximation of the alignment using exhaustive search. For instance, for a global rotational alignment use",
false);
171 addExampleLine(
"xmipp_volume_align --i1 volume1.vol --i2 volume2.vol --rot 0 360 15 --tilt 0 180 15 --psi 0 360 15");
172 addExampleLine(
"Then, assume the best alignment is obtained for rot=45, tilt=60, psi=90",
false);
173 addExampleLine(
"Now you perform a local search to refine the estimation and apply",
false);
174 addExampleLine(
"xmipp_volume_align --i1 volume1.vol --i2 volume2.vol --rot 45 --tilt 60 --psi 90 --local --apply volume2aligned.vol");
175 addSeeAlsoLine(
"xmipp_volumeset_align");
181 fn1 = getParam(
"--i1");
182 fn2 = getParam(
"--i2");
184 rot0 = getDoubleParam(
"--rot",0);
185 rotF = getDoubleParam(
"--rot",1);
186 step_rot = getDoubleParam(
"--rot",2);
188 tilt0 = getDoubleParam(
"--tilt",0);
189 tiltF = getDoubleParam(
"--tilt",1);
190 step_tilt = getDoubleParam(
"--tilt",2);
192 psi0 = getDoubleParam(
"--psi",0);
193 psiF = getDoubleParam(
"--psi",1);
194 step_psi = getDoubleParam(
"--psi",2);
196 scale0 = getDoubleParam(
"--scale",0);
197 scaleF = getDoubleParam(
"--scale",1);
198 step_scale = getDoubleParam(
"--scale",2);
200 grey_scale0 = getDoubleParam(
"--grey_scale",0);
201 grey_scaleF = getDoubleParam(
"--grey_scale",1);
202 step_grey = getDoubleParam(
"--grey_scale",2);
204 grey_shift0 = getDoubleParam(
"--grey_shift",0);
205 grey_shiftF = getDoubleParam(
"--grey_shift",1);
206 step_grey_shift = getDoubleParam(
"--grey_shift",2);
208 z0 = getDoubleParam(
"-z",0);
209 zF = getDoubleParam(
"-z",1);
210 step_z = getDoubleParam(
"-z",2);
212 y0 = getDoubleParam(
"-y",0);
213 yF = getDoubleParam(
"-y",1);
214 step_y = getDoubleParam(
"-y",2);
216 x0 = getDoubleParam(
"-x",0);
217 xF = getDoubleParam(
"-x",1);
218 step_x = getDoubleParam(
"-x",2);
220 mask_enabled = checkParam(
"--mask");
222 mask.
read(argc, argv);
224 usePowell = checkParam(
"--local");
225 useFRM = checkParam(
"--frm");
228 maxFreq=getDoubleParam(
"--frm",0);
229 maxShift=getIntParam(
"--frm",1);
230 starting_tilt=getIntParam(
"--frm",2);
231 ending_tilt=getIntParam(
"--frm",3);
233 onlyShift = checkParam(
"--onlyShift");
234 wrap = !checkParam(
"--dontWrap");
246 if (step_grey_shift == 0)
255 tell = checkParam(
"--show_fit");
256 apply = checkParam(
"--apply");
257 fnOut = getParam(
"--apply");
258 copyGeo = checkParam(
"--copyGeo");
259 copyGray = checkParam(
"--copyGray");
260 fnGeo = getParam(
"--copyGeo");
261 fnGray = getParam(
"--copyGray");
262 store = checkParam(
"--store");
263 fnStore = getParam(
"--store");
264 dontScale = checkParam(
"--dontScale");
266 if (checkParam(
"--covariance"))
270 else if (checkParam(
"--least_squares"))
286 params.
V1().setXmippOrigin();
288 params.
V2().setXmippOrigin();
292 double best_fit = 1e38;
306 if (!usePowell && !useFRM)
312 if (grey_scale0 != grey_scaleF)
313 times *=
FLOOR(1 + (grey_scaleF - grey_scale0) / step_grey);
314 if (grey_shift0 != grey_shiftF)
315 times *=
FLOOR(1 + (grey_shiftF - grey_shift0) / step_grey_shift);
317 times *=
FLOOR(1 + (rotF - rot0) / step_rot);
319 times *=
FLOOR(1 + (tiltF - tilt0) / step_tilt);
321 times *=
FLOOR(1 + (psiF - psi0) / step_psi);
322 if (scale0 != scaleF)
323 times *=
FLOOR(1 + (scaleF - scale0) / step_scale);
325 times *=
FLOOR(1 + (zF - z0) / step_z);
327 times *=
FLOOR(1 + (yF - y0) / step_y);
329 times *=
FLOOR(1 + (xF - x0) / step_x);
333 std::cout <<
"#grey_factor rot tilt psi scale z y x fitness\n";
337 int step_time =
CEIL((
double)times / 60.0);
340 for (
double grey_scale = grey_scale0; grey_scale <= grey_scaleF ; grey_scale += step_grey)
341 for (
double grey_shift = grey_shift0; grey_shift <= grey_shiftF ; grey_shift += step_grey_shift)
342 for (
double rot = rot0; rot <= rotF ; rot += step_rot)
343 for (
double tilt = tilt0; tilt <= tiltF ; tilt += step_tilt)
344 for (
double psi = psi0;
psi <= psiF ;
psi += step_psi)
345 for (
double scale = scale0; scale <= scaleF ; scale += step_scale)
346 for (
ZZ(r) = z0;
ZZ(r) <=
zF ;
ZZ(r) += step_z)
347 for (
YY(r) =
y0;
YY(r) <=
yF ;
YY(r) += step_y)
348 for (
XX(r) =
x0;
XX(r) <=
xF ;
XX(r) += step_x)
351 trial(0) = grey_scale;
352 trial(1) = grey_shift;
363 if (fit < best_fit || first)
369 std::cout <<
"Best so far\n";
373 std::cout << trial <<
" " << fit << std::endl;
375 if (++itime % step_time == 0)
413 double rot,tilt,
psi,
x,
y,
z,score;
415 if(starting_tilt!=-90 || ending_tilt!=90){
416 std::cout<<
"you are compensating for the missing wedge, the first volume should be rotated with 90 degrees about the y-axis"<<std::endl;
417 PyObject * pSTMMclass =
Python::getClassRef(
"sh_alignment.tompy.filter",
"SingleTiltWedge");
418 PyObject * arglist = Py_BuildValue(
"(ii)", starting_tilt , ending_tilt);
419 PyObject * SingleTiltWedgeMask = PyObject_CallObject(pSTMMclass, arglist);
421 alignVolumesFRM(pFunc, params.
V2(), params.
V1(), SingleTiltWedgeMask, rot,tilt,
psi,
x,
y,
z,score,A,maxShift,maxFreq,params.
mask_ptr);
422 std::cout<<
"If you intend to apply transform using xmipp_transform_geometry, use --inverse flag (if it was not present before), or remove it (if it was present before)"<<std::endl;
423 Py_DECREF(SingleTiltWedgeMask);
425 Py_DECREF(pSTMMclass);
428 alignVolumesFRM(pFunc, params.
V1(), params.
V2(), Py_None, rot,tilt,
psi,
x,
y,
z,score,A,maxShift,maxFreq,params.
mask_ptr);
445 std::cout <<
"The best correlation is for\n" 446 <<
"Scale : " << best_align(5) << std::endl
447 <<
"Translation (X,Y,Z) : " << best_align(8)
448 <<
" " << best_align(7) <<
" " << best_align(6)
450 <<
"Rotation (rot,tilt,psi): " 451 << best_align(2) <<
" " << best_align(3) <<
" " 452 << best_align(4) << std::endl
453 <<
"Best grey scale : " << best_align(0) << std::endl
454 <<
"Best grey shift : " << best_align(1) << std::endl
455 <<
"Fitness value : " << best_fit << std::endl;
457 XX(r) = best_align(8);
458 YY(r) = best_align(7);
459 ZZ(r) = best_align(6);
478 std::cout <<
"xmipp_transform_geometry will require the following values" 479 <<
"\n Angles: " << best_align(2) <<
" " 480 << best_align(3) <<
" " << best_align(4)
481 <<
"\n Shifts: " << A(0,3) <<
" " << A(1,3) <<
" " << A(2,3)
485 std::ofstream outputGeo (fnGeo.c_str());
486 outputGeo << A(0,0) <<
"\n" << A(0,1) <<
"\n" << A(0,2) <<
"\n" << A(0,3) <<
"\n" 487 << A(1,0) <<
"\n" << A(1,1) <<
"\n" << A(1,2) <<
"\n" << A(1,3) <<
"\n" 488 << A(2,0) <<
"\n" << A(2,1) <<
"\n" << A(2,2) <<
"\n" << A(2,3) <<
"\n" 489 << A(3,0) <<
"\n" << A(3,1) <<
"\n" << A(3,2) <<
"\n" << A(3,3) <<
"\n" 495 std::ofstream outputGray (fnGray.c_str());
496 outputGray << best_align(0) <<
"\n" 497 << best_align(1) <<
"\n" 503 std::ofstream outputStore (fnStore.c_str());
504 outputStore << best_align(2) <<
", " << best_align(3) <<
", " << best_align(4) <<
", " << A(0,3) <<
", " 505 << A(1,3) <<
", " << A(2,3) <<
", " << best_fit << std::endl;
511 params.
V2()=params.
Vaux();
void init_progress_bar(long total)
double rms(const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=nullptr, MultidimArray< double > *Contributions=nullptr)
double wrapperFitness(double *p, void *params)
void Euler_angles2matrix(T alpha, T beta, T gamma, Matrix2D< T > &A, bool homogeneous)
Matrix1D< double > vectorR3(double x, double y, double z)
void write(const FileName &name="", size_t select_img=ALL_IMAGES, bool isStack=false, int mode=WRITE_OVERWRITE, CastWriteMode castMode=CW_CAST, int _swapWrite=0)
double correlationIndex(const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=NULL, MultidimArray< double > *Contributions=NULL)
void applyTransformation(const MultidimArray< double > &V2, MultidimArray< double > &Vaux, double *p, bool wrap)
void initPythonAndNumpy()
static void defineParams(XmippProgram *program, int allowed_data_types=ALL_KINDS, const char *prefix=nullptr, const char *comment=nullptr, bool moreOptions=false)
void read(int argc, const char **argv)
PyObject * getFunctionRef(const std::string &moduleName, const std::string &funcName)
void progress_bar(long rlen)
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
#define DIRECT_MULTIDIM_ELEM(v, n)
const MultidimArray< int > * mask_ptr
PyObject * getClassRef(const std::string &moduleName, const std::string &className)
void generate_mask(bool apply_geo=false)
#define MATRIX1D_ARRAY(v)
void alignVolumesFRM(PyObject *pFunc, const MultidimArray< double > &Iref, MultidimArray< double > &I, PyObject *Imask, double &rot, double &tilt, double &psi, double &x, double &y, double &z, double &score, Matrix2D< double > &A, int maxshift, double maxFreq, const MultidimArray< int > *mask)
double psi(const double x)
int read(const FileName &name, DataMode datamode=DATA, size_t select_img=ALL_IMAGES, bool mapData=false, int mode=WRITE_READONLY)
const MultidimArray< int > & get_binary_mask() const
double fitness(double *p)