EMAN2
aligner.cpp
Go to the documentation of this file.
1/*
2 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
3 * Copyright (c) 2000-2006 Baylor College of Medicine
4 *
5 * This software is issued under a joint BSD/GNU license. You may use the
6 * source code in this file under either license. However, note that the
7 * complete EMAN2 and SPARX software packages have some GPL dependencies,
8 * so you are responsible for compliance with the licenses of these packages
9 * if you opt to use BSD licensing. The warranty disclaimer below holds
10 * in either instance.
11 *
12 * This complete copyright notice must be included in any revised version of the
13 * source code. Additional authorship citations may be added, but existing
14 * author citations must be preserved.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 *
30 * */
31#include <random>
32
33#include "emfft.h"
34#include "cmp.h"
35#include "aligner.h"
36#include "averager.h"
37#include "emdata.h"
38#include "processor.h"
39#include "util.h"
40#include "symmetry.h"
41#include <gsl/gsl_multimin.h>
43
44#ifdef EMAN2_USING_CUDA
45 #include "cuda/cuda_processor.h"
46 #include "cuda/cuda_cmp.h"
47#endif
48
49#ifdef SPARX_USING_CUDA
50 #include <sparx/cuda/cuda_ccf.h>
51#endif
52
53#define EMAN2_ALIGNER_DEBUG 0
54
55using namespace EMAN;
56
57const string TranslationalAligner::NAME = "translational";
58const string RotationalAligner::NAME = "rotational";
59const string RotationalAlignerBispec::NAME = "rotational_bispec";
60const string RotationalAlignerIterative::NAME = "rotational_iterative";
61const string RotatePrecenterAligner::NAME = "rotate_precenter";
62const string RotateTranslateAligner::NAME = "rotate_translate";
63const string RotateTranslateAlignerBispec::NAME = "rotate_translate_bispec";
64const string RotateTranslateScaleAligner::NAME = "rotate_translate_scale";
65const string RotateTranslateAlignerIterative::NAME = "rotate_translate_iterative";
66const string RotateTranslateScaleAlignerIterative::NAME = "rotate_trans_scale_iter";
67const string RotateTranslateAlignerPawel::NAME = "rotate_translate_resample";
68const string RotateTranslateBestAligner::NAME = "rotate_translate_best";
69const string RT2DTreeAligner::NAME = "rotate_translate_tree";
70const string RotateFlipAligner::NAME = "rotate_flip";
71const string RotateFlipAlignerIterative::NAME = "rotate_flip_iterative";
72const string RotateTranslateFlipAligner::NAME = "rotate_translate_flip";
73const string RotateTranslateFlipScaleAligner::NAME = "rotate_trans_flip_scale";
74const string RotateTranslateFlipAlignerIterative::NAME = "rotate_translate_flip_iterative";
75const string RotateTranslateFlipScaleAlignerIterative::NAME = "rotate_trans_flip_scale_iter";
76const string RotateTranslateFlipAlignerPawel::NAME = "rotate_translate_flip_resample";
77const string RTFExhaustiveAligner::NAME = "rtf_exhaustive";
78const string RTFSlowExhaustiveAligner::NAME = "rtf_slow_exhaustive";
79const string RefineAligner::NAME = "refine";
80const string RefineAlignerCG::NAME = "refinecg";
81const string SymAlignProcessorQuat::NAME = "symalignquat";
82const string SymAlignProcessor::NAME = "symalign";
83const string Refine3DAlignerGrid::NAME = "refine_3d_grid";
84const string Refine3DAlignerQuaternion::NAME = "refine_3d";
85const string RT3DGridAligner::NAME = "rotate_translate_3d_grid";
86const string RT3DSphereAligner::NAME = "rotate_translate_3d";
87const string RT3DTreeAligner::NAME = "rotate_translate_3d_tree";
88const string RT3DLocalTreeAligner::NAME = "rotate_translate_3d_local_tree";
89const string RT2Dto3DTreeAligner::NAME = "rotate_translate_2d_to_3d_tree";
90const string RT3DSymmetryAligner::NAME = "rotate_symmetry_3d";
91const string FRM2DAligner::NAME = "frm2d";
92const string ScaleAligner::NAME = "scale";
93
94
96{
97 force_add<TranslationalAligner>();
98 force_add<RotationalAligner>();
99 force_add<RotationalAlignerBispec>();
100 force_add<RotationalAlignerIterative>();
101 force_add<RotatePrecenterAligner>();
102 force_add<RotateTranslateAligner>();
103 force_add<RotateTranslateAlignerBispec>();
104 force_add<RotateTranslateScaleAligner>();
105 force_add<RotateTranslateAlignerIterative>();
106 force_add<RotateTranslateScaleAlignerIterative>();
107 force_add<RotateTranslateAlignerPawel>();
108 force_add<RT2DTreeAligner>();
109 force_add<RotateFlipAligner>();
110 force_add<RotateFlipAlignerIterative>();
111 force_add<RotateTranslateFlipAligner>();
112 force_add<RotateTranslateFlipScaleAligner>();
113 force_add<RotateTranslateFlipAlignerIterative>();
114 force_add<RotateTranslateFlipScaleAlignerIterative>();
115 force_add<RotateTranslateFlipAlignerPawel>();
116 force_add<RTFExhaustiveAligner>();
117 force_add<RTFSlowExhaustiveAligner>();
118 force_add<SymAlignProcessor>();
119 force_add<RefineAligner>();
120 force_add<RefineAlignerCG>();
121 force_add<SymAlignProcessorQuat>();
122 force_add<Refine3DAlignerGrid>();
123 force_add<Refine3DAlignerQuaternion>();
124 force_add<RT3DGridAligner>();
125 force_add<RT3DSphereAligner>();
126 force_add<RT2Dto3DTreeAligner>();
127 force_add<RT3DTreeAligner>();
128 force_add<RT3DLocalTreeAligner>();
129 force_add<RT3DSymmetryAligner>();
130 force_add<FRM2DAligner>();
131 force_add<ScaleAligner>();
132// force_add<XYZAligner>();
133}
134
135vector<Dict> Aligner::xform_align_nbest(EMData *, EMData *, const unsigned int, const string &, const Dict&) const
136{
137 vector<Dict> solns;
138 return solns;
139}
140
142 const string & cmp_name, const Dict& cmp_params) const
143{
144 //get the scale range
145 float min = params.set_default("min",0.95f);
146 float max = params.set_default("max",1.05f);
147 float step = params.set_default("step",0.01f);
148
149 // crate the starting transform
150 Transform t = Transform();
151 t.set_scale(max);
152
153 //save orignal data
154 float* oridata = this_img->get_data();
155
156 //get the transform processor and cast to correct factory product
157 Processor* proc = Factory <Processor>::get("xform", Dict());
158 TransformProcessor* xform = dynamic_cast<TransformProcessor*>(proc);
159
160 // Using the following method we only create one EMdata object. If I just used the processor, then I would create many EMdata objects
161 EMData* result = 0;
162// float bestscore = numeric_limits<float>::infinity();
163 float bestscore = 1.0e37;
164
165 for(float i = max; i > min; i-=step){
166
167 //scale the image
168 float* des_data = xform->transform(this_img,t);
169 this_img->set_data(des_data);
170 this_img->update();
171
172 //check compairsion
173 EMData* aligned = this_img->align(basealigner, to, basealigner_params, cmp_name, cmp_params);
174 float score = aligned->cmp(cmp_name, to, cmp_params);
175 if(score < bestscore){
176 bestscore = score;
177 //If we just reassign w/o cleaing up we will get memory leaks!
178 if(result != 0) delete result;
179 result = aligned;
180 Transform *tt=(Transform*)(result->get_attr("xform.align2d"));
181 tt->set_scale(i);
182 result->set_attr("xform.align2d",tt);
183// result->set_attr("scalefactor",i); // I don't know why this was done instead of updating the transform, but it breaks a bunch of things
184 }else{
185 delete aligned;
186 }
187 //clean up scaled image data
188 delete des_data;
189
190 t.set_scale(i);
191
192 //reset original data
193 this_img->set_data(oridata);
194 }
195
196 if (!result) throw UnexpectedBehaviorException("Alignment score is infinity! Something is seriously wrong with the data!");
197 if (proc != 0) delete proc;
198
199 return result;
200
201};
202
203EMData* ScaleAligner::align(EMData * this_img, EMData *to,
204 const string& cmp_name, const Dict& cmp_params) const
205{
206
207 //get the scale range
208 float min = params.set_default("min",0.95f);
209 float max = params.set_default("max",1.05f);
210 float step = params.set_default("step",0.01f);
211
212 Transform t = Transform();
213 t.set_scale(max);
214 float* oridata = this_img->get_data();
215
216 //get the transform processor and cast to correct factory product
217 Processor* proc = Factory <Processor>::get("xform", Dict());
218 TransformProcessor* xform = dynamic_cast<TransformProcessor*>(proc);
219
220 // Using the following method we only create one EMdata object. If I just used the processor, then I would create many EMdata objects
221 float bestscale = 1.0;
222 float bestscore = 1.0e37;
223
224 for(float i = max; i > min; i-=step){
225
226 float* des_data = xform->transform(this_img,t);
227 this_img->set_data(des_data);
228 this_img->update();
229
230 //check compairsion
231 float score = this_img->cmp(cmp_name, to, cmp_params);
232 if(score < bestscore){
233 bestscore = score;
234 bestscale = i;
235 }
236 //clean up scaled image data
237 delete des_data;
238
239 t.set_scale(i);
240
241 //reset original data
242 this_img->set_data(oridata);
243 }
244
245
246
247 //Return scaled image
248 t.set_scale(bestscale);
249 EMData* result = this_img->process("xform",Dict("transform",&t));
250 result->set_attr("scalefactor",bestscale);
251 if (proc != 0) delete proc;
252
253 return result;
254
255}
256
257// Note, the translational aligner assumes that the correlation image
258// generated by the calc_ccf function is centered on the bottom left corner
259// That is, if you did at calc_cff using identical images, the
260// peak would be at 0,0
262 const string&, const Dict&) const
263{
264 if (!this_img) {
265 return 0;
266 }
267
268 if (to && !EMUtil::is_same_size(this_img, to))
269 throw ImageDimensionException("Images must be the same size to perform translational alignment");
270
271 EMData *cf = 0;
272 int nx = this_img->get_xsize();
273 int ny = this_img->get_ysize();
274 int nz = this_img->get_zsize();
275
276 int masked = params.set_default("masked",0);
277 int useflcf = params.set_default("useflcf",0);
278 bool use_cpu = true;
279
280#ifdef EMAN2_USING_CUDA
281 if(EMData::usecuda == 1) {
282 //if(!this_img->getcudarwdata()) this_img->copy_to_cuda();
283 //if(to && !to->getcudarwdata()) to->copy_to_cuda();
284 //if (masked) throw UnexpectedBehaviorException("Masked is not yet supported in CUDA");
285 //if (useflcf) throw UnexpectedBehaviorException("Useflcf is not yet supported in CUDA");
286 //cout << "Translate on GPU" << endl;
287 //use_cpu = false;
288 //cf = this_img->calc_ccf(to);
289 }
290#endif // EMAN2_USING_CUDA
291
292 if (use_cpu) {
293 if (useflcf) cf = this_img->calc_flcf(to);
294 else cf = this_img->calc_ccf(to);
295 }
296 //return cf;
297 // This is too expensive, esp for CUDA(we we can fix later
298 if (masked) {
299 EMData *msk=this_img->process("threshold.notzero");
300 EMData *sqr=to->process("math.squared");
301 EMData *cfn=msk->calc_ccf(sqr);
302 cfn->process_inplace("math.sqrt");
303 float *d1=cf->get_data();
304 float *d2=cfn->get_data();
305 for (size_t i=0; i<(size_t)nx*ny*nz; ++i) {
306 if (d2[i]!=0) d1[i]/=d2[i];
307 }
308 cf->update();
309 delete msk;
310 delete sqr;
311 delete cfn;
312 }
313
314 int maxshiftx = params.set_default("maxshift",-1);
315 int maxshifty = params["maxshift"];
316 int maxshiftz = params["maxshift"];
317 int nozero = params["nozero"];
318
319 if (maxshiftx <= 0) {
320 maxshiftx = nx / 4;
321 maxshifty = ny / 4;
322 maxshiftz = nz / 4;
323 }
324
325 if (maxshiftx > nx / 2 - 1) maxshiftx = nx / 2 - 1;
326 if (maxshifty > ny / 2 - 1) maxshifty = ny / 2 - 1;
327 if (maxshiftz > nz / 2 - 1) maxshiftz = nz / 2 - 1;
328
329 if (nx == 1) maxshiftx = 0; // This is justhere for completeness really... plus it saves errors
330 if (ny == 1) maxshifty = 0;
331 if (nz == 1) maxshiftz = 0;
332
333 // If nozero the portion of the image in the center (and its 8-connected neighborhood) is zeroed
334 if (nozero) {
336 }
337
338 IntPoint peak;
339#ifdef EMAN2_USING_CUDA
340 if (!use_cpu) {
341 cout << "USe CUDA TA 2" << endl;
342 if (nozero) throw UnexpectedBehaviorException("Nozero is not yet supported in CUDA");
343 CudaPeakInfo* data = calc_max_location_wrap_cuda(cf->getcudarwdata(), cf->get_xsize(), cf->get_ysize(), cf->get_zsize(), maxshiftx, maxshifty, maxshiftz);
344 peak = IntPoint(data->px,data->py,data->pz);
345 free(data);
346 }
347#endif // EMAN2_USING_CUDA
348
349 float maxvalue=0;
350 if (use_cpu) {
351 peak = cf->calc_max_location_wrap(maxshiftx, maxshifty, maxshiftz, &maxvalue);
352 }
353 //cout << -peak[0] << " " << -peak[1] << " " << -peak[2] << endl;
354 Vec3f cur_trans = Vec3f ( (float)-peak[0], (float)-peak[1], (float)-peak[2]);
355 //cout << peak[0] << " " << peak[1] << endl;
356
357 if (!to) {
358 cur_trans /= 2.0f; // If aligning theimage to itself then only go half way -
359 int intonly = params.set_default("intonly",false);
360 if (intonly) {
361 cur_trans[0] = floor(cur_trans[0] + 0.5f);
362 cur_trans[1] = floor(cur_trans[1] + 0.5f);
363 cur_trans[2] = floor(cur_trans[2] + 0.5f);
364 }
365 }
366
367 if( cf ){
368 delete cf;
369 cf = 0;
370 }
371
372 Dict params("trans",static_cast< vector<int> >(cur_trans));
373 if (use_cpu){
374 cf=this_img->process("xform.translate.int",params);
375 }
376 Transform t;
377 t.set_trans(cur_trans);
378
379#ifdef EMAN2_USING_CUDA
380 if (!use_cpu) {
381 cout << "USe CUDA TA 3" << endl;
382 //this will work just fine....
383 cf = this_img->process("xform",Dict("transform",&t));
384 }
385#endif // EMAN2_USING_CUDA
386
387 if ( nz != 1 ) {
388// Transform* t = get_set_align_attr("xform.align3d",cf,this_img);
389// t->set_trans(cur_trans);
390 cf->set_attr("xform.align3d",&t);
391 cf->set_attr("score_align",-maxvalue);
392 } else if ( ny != 1 ) {
393 //Transform* t = get_set_align_attr("xform.align2d",cf,this_img);
394 cur_trans[2] = 0; // just make sure of it
395 t.set_trans(cur_trans);
396 cf->set_attr("xform.align2d",&t);
397 cf->set_attr("score_align",-maxvalue);
398 }
399 return cf;
400}
401
402EMData * RotationalAlignerBispec::align(EMData * this_img, EMData *to, const string& cmp_name, const Dict& cmp_params) const {
403 // Make translationally invariant rotational footprints
404 EMData* this_img_bispec, * to_bispec;
405 int rfp = params.set_default("rfpn",4);
406 int size = params.set_default("size",32);
407 int harmonic = params.set_default("harmonic",0);
408 if (harmonic) {
409 this_img_bispec=this_img->process("math.harmonicpow",Dict("rfp",1));
410 to_bispec=to->process("math.harmonicpow",Dict("rfp",1));
411 } else {
412 this_img_bispec=this_img->process("math.bispectrum.slice",Dict("rfp",rfp,"size",size));
413 to_bispec=to->process("math.bispectrum.slice",Dict("rfp",rfp,"size",size));
414 }
415 int this_img_rfp_nx = this_img_bispec->get_xsize();
416
417 // Do row-wise correlation, returning a sum.
418 EMData *cf = this_img_bispec->calc_ccfx(to_bispec, 0, this_img->get_ysize(),false,false,false);
419
420 // Delete them, they're no longer needed
421 delete this_img_bispec;
422 delete to_bispec;
423
424 // Now solve the rotational alignment by finding the max in the column sum
425 float *data = cf->get_data();
426
427 float peak = 0;
428 int peak_index = 0;
429 Util::find_max(data, this_img_rfp_nx, &peak, &peak_index);
430
431 delete cf;
432
433 float rot_angle = (float) (peak_index * 360.0f / this_img_rfp_nx);
434
435 // Return the result
436 Transform tmp(Dict("type","2d","alpha",rot_angle));
437 cf=this_img->process("xform",Dict("transform",(Transform*)&tmp));
438// Transform* t = get_set_align_attr("xform.align2d",cf,this_img);
439// Dict d("type","2d","alpha",rot_angle);
440// t->set_rotation(d);
441 cf->set_attr("xform.align2d",&tmp);
442 return cf;
443}
444
445
446EMData * RotationalAligner::align_180_ambiguous(EMData * this_img, EMData * to, int rfp_mode,int zscore) {
447
448 // Make translationally invariant rotational footprints
449 EMData* this_img_rfp, * to_rfp;
450 if (rfp_mode == 0) {
451 this_img_rfp = this_img->make_rotational_footprint_e1();
452 to_rfp = to->make_rotational_footprint_e1();
453 } else if (rfp_mode == 1) {
454 this_img_rfp = this_img->make_rotational_footprint();
455 to_rfp = to->make_rotational_footprint();
456 } else if (rfp_mode == 2) {
457 this_img_rfp = this_img->make_rotational_footprint_cmc();
458 to_rfp = to->make_rotational_footprint_cmc();
459 } else {
460 throw InvalidParameterException("rfp_mode must be 0,1 or 2");
461 }
462 int this_img_rfp_nx = this_img_rfp->get_xsize();
463
464 // Do row-wise correlation, returning a sum.
465 EMData *cf = this_img_rfp->calc_ccfx(to_rfp, 0, this_img->get_ysize(),false,false,zscore);
466// cf->process_inplace("normalize");
467// cf->write_image("ralisum.hdf",-1);
468//
469// EMData *cf2 = this_img_rfp->calc_ccfx(to_rfp, 0, this_img->get_ysize(),true);
470// cf2->write_image("ralistack.hdf",-1);
471// delete cf2;
472
473 // Delete them, they're no longer needed
474 delete this_img_rfp; this_img_rfp = 0;
475 delete to_rfp; to_rfp = 0;
476
477 // Now solve the rotational alignment by finding the max in the column sum
478 float *data = cf->get_data();
479
480 float peak = 0;
481 int peak_index = 0;
482 Util::find_max(data, this_img_rfp_nx, &peak, &peak_index);
483
484 if( cf ) {
485 delete cf;
486 cf = 0;
487 }
488 float rot_angle = (float) (peak_index * 180.0f / this_img_rfp_nx);
489
490 // Return the result
491 Transform tmp(Dict("type","2d","alpha",rot_angle));
492 cf=this_img->process("xform",Dict("transform",(Transform*)&tmp));
493// Transform* t = get_set_align_attr("xform.align2d",cf,this_img);
494// Dict d("type","2d","alpha",rot_angle);
495// t->set_rotation(d);
496 cf->set_attr("xform.align2d",&tmp);
497 return cf;
498}
499
501 const string& cmp_name, const Dict& cmp_params) const
502{
503 if (!to) throw InvalidParameterException("Can not rotational align - the image to align to is NULL");
504
505#ifdef EMAN2_USING_CUDA
506 if(EMData::usecuda == 1) {
507 //if(!this_img->getcudarwdata()) this_img->copy_to_cuda();
508 //if(!to->getcudarwdata()) to->copy_to_cuda();
509 }
510#endif
511
512 // Perform 180 ambiguous alignment
513 int rfp_mode = params.set_default("rfp_mode",2);
514 int zscore = params.set_default("zscore",0);
515 int ambig180 = params.set_default("ambig180",0);
516
517 EMData* rot_aligned = RotationalAligner::align_180_ambiguous(this_img,to,rfp_mode,zscore);
518 Transform * tmp = rot_aligned->get_attr("xform.align2d");
519 Dict rot = tmp->get_rotation("2d");
520 float rotate_angle_solution = rot["alpha"];
521 delete tmp;
522
523 // Don't resolve the 180 degree ambiguity here
524 if (ambig180) {
525 return rot_aligned;
526 }
527
528 EMData *rot_align_180 = rot_aligned->process("math.rotate.180");
529
530 // Generate the comparison metrics for both rotational candidates
531 float rot_cmp = rot_aligned->cmp(cmp_name, to, cmp_params);
532 float rot_180_cmp = rot_align_180->cmp(cmp_name, to, cmp_params);
533
534 // Decide on the result
535 float score = 0.0;
536 EMData* result = NULL;
537 if (rot_cmp < rot_180_cmp){
538 result = rot_aligned;
539 score = rot_cmp;
540 delete rot_align_180; rot_align_180 = 0;
541 } else {
542 result = rot_align_180;
543 score = rot_180_cmp;
544 delete rot_aligned; rot_aligned = 0;
545 rotate_angle_solution = rotate_angle_solution-180.0f;
546 }
547
548// Transform* t = get_align_attr("xform.align2d",result);
549// t->set_rotation(Dict("type","2d","alpha",rotate_angle_solution));
550 Transform tmp2(Dict("type","2d","alpha",rotate_angle_solution));
551 result->set_attr("xform.align2d",&tmp2);
552 return result;
553}
554
555
557 const string&, const Dict&) const
558{
559 if (!to) {
560 return 0;
561 }
562
563 int ny = this_img->get_ysize();
564 int size = Util::calc_best_fft_size((int) (M_PI * ny * 1.5));
565 EMData *e1 = this_img->unwrap(4, ny * 7 / 16, size, 0, 0, 1);
566 EMData *e2 = to->unwrap(4, ny * 7 / 16, size, 0, 0, 1);
567 EMData *cf = e1->calc_ccfx(e2, 0, ny);
568
569 float *data = cf->get_data();
570
571 float peak = 0;
572 int peak_index = 0;
573 Util::find_max(data, size, &peak, &peak_index);
574 float a = (float) ((1.0f - 1.0f * peak_index / size) * 180. * 2);
575
576 Transform rot;
577 rot.set_rotation(Dict("type","2d","alpha",(float)(a*180./M_PI)));
578 EMData* rslt = this_img->process("xform",Dict("transform",&rot));
579 rslt->set_attr("xform.align2d",&rot);
580//
581// Transform* t = get_set_align_attr("xform.align2d",rslt,this_img);
582// t->set_rotation(Dict("type","2d","alpha",-a));
583//
584// EMData* result this_img->transform(Dict("type","2d","alpha",(float)(a*180./M_PI)));
585//
586// cf->set_attr("xform.align2d",t);
587// delete t;
588// cf->update();
589
590 if( e1 )
591 {
592 delete e1;
593 e1 = 0;
594 }
595
596 if( e2 )
597 {
598 delete e2;
599 e2 = 0;
600 }
601
602 if (cf) {
603 delete cf;
604 cf = 0;
605 }
606 return rslt;
607}
608
610 const string &, const Dict&) const
611{
612 int r1 = params.set_default("r1",-1);
613 int r2 = params.set_default("r2",-1);
614 //to start lest try the original size image. If needed, we can pad it....
615 EMData * to_polar = to->unwrap(r1,r2,-1,0,0,true);
616 EMData * this_img_polar = this_img->unwrap(r1,r2,-1,0,0,true);
617 int this_img_polar_nx = this_img_polar->get_xsize();
618
619 EMData * cf = this_img_polar->calc_ccfx(to_polar, 0, this_img->get_ysize());
620
621 //take out the garbage
622 delete to_polar; to_polar = 0;
623 delete this_img_polar; this_img_polar = 0;
624
625 float * data = cf->get_data();
626 float peak = 0;
627 int peak_index = 0;
628 Util::find_max(data, this_img_polar_nx, &peak, &peak_index);
629
630 delete cf; cf = 0;
631 float rot_angle = (float) (peak_index * 360.0f / this_img_polar_nx);
632
633 //return the result
634 //cout << rot_angle << endl;
635 Transform tmp(Dict("type","2d","alpha",rot_angle));
636 EMData * rotimg=this_img->process("xform",Dict("transform",(Transform*)&tmp));
637 rotimg->set_attr("xform.align2d",&tmp);
638
639 return rotimg;
640
641}
642
644 const string & cmp_name, const Dict& cmp_params) const
645{
646 int maxiter = params.set_default("maxiter", 3);
647
648 Dict trans_params;
649 trans_params["intonly"] = 0;
650 trans_params["maxshift"] = params.set_default("maxshift", -1);
651 trans_params["useflcf"] = params.set_default("useflcf",0);
652 trans_params["nozero"] = params.set_default("nozero",false);
653
654 Dict rot_params;
655 rot_params["r1"] = params.set_default("r1", -1);
656 rot_params["r2"] = params.set_default("r2", -1);
657
658 Transform t;
659 EMData * moving_img = this_img;
660 for(int it = 0; it < maxiter; it++){
661
662 // First do a translational alignment
663 EMData * trans_align = moving_img->align("translational", to, trans_params, cmp_name, cmp_params);
664 Transform * tt = trans_align->get_attr("xform.align2d");
665 t = *tt*t;
666 delete tt;
667
668 //now do rotation
669 EMData * rottrans_align = trans_align->align("rotational_iterative", to, rot_params, cmp_name, cmp_params);
670 Transform * rt = rottrans_align->get_attr("xform.align2d");
671 t = *rt*t;
672 delete trans_align; trans_align = 0;
673 delete rottrans_align; rottrans_align = 0;
674 delete rt;
675
676 //this minimizes interpolation errors (all images that are futher processed will be interpolated at most twice)
677 if(it > 0){delete moving_img;}
678
679 moving_img = this_img->process("xform",Dict("transform",&t)); //iterate
680 }
681
682 //write the total transformation; Avoids interpolation erros
683 moving_img->set_attr("xform.align2d", &t);
684
685 return moving_img;
686}
687
689 const string & cmp_name, const Dict& cmp_params) const
690{
691
692 //Basically copy params into rotate_translate
693 basealigner_params["maxshift"] = params.set_default("maxshift", -1);
694 basealigner_params["r1"] = params.set_default("r1",-1);
695 basealigner_params["r2"] = params.set_default("r2",-1);
696 basealigner_params["maxiter"] = params.set_default("maxiter",3);
697 basealigner_params["nozero"] = params.set_default("nozero",false);
698 basealigner_params["useflcf"] = params.set_default("useflcf",0);
699
700 //return the correct results
701 return align_using_base(this_img, to, cmp_name, cmp_params);
702
703}
704
706 const string & cmp_name, const Dict& cmp_params) const
707{
708 if (cmp_name != "dot" && cmp_name != "ccc") throw InvalidParameterException("Resample aligner only works for dot and ccc");
709
710 int maxtx = params.set_default("tx", 0);
711 int maxty = params.set_default("ty", 0);
712 int r1 = params.set_default("r1",-1);
713 int r2 = params.set_default("r2",-1);
714
715 if(this_img->get_xsize()/2 - 1 - r2 - maxtx <= 0 || (r2 == -1 && maxtx > 0)) throw InvalidParameterException("nx/2 - 1 - r2 - tx must be greater than or = 0");
716 if(this_img->get_ysize()/2 - 1 - r2 - maxty <= 0 || (r2 == -1 && maxty > 0)) throw InvalidParameterException("ny/2 - 1 - r2 - ty must be greater than or = 0");
717
718// float best_peak = -numeric_limits<float>::infinity();
719 float best_peak = -1.0e37;
720 int best_peak_index = 0;
721 int best_tx = 0;
722 int best_ty = 0;
723 int polarxsize = 0;
724
725 for(int x = -maxtx; x <= maxtx; x++){
726 for(int y = -maxty; y <= maxty; y++){
727
728 EMData * to_polar = to->unwrap(r1,r2,-1,0,0,true);
729 EMData * this_img_polar = this_img->unwrap(r1,r2,-1,x,y,true);
730 EMData * cf = this_img_polar->calc_ccfx(to_polar, 0, this_img_polar->get_ysize());
731
732 polarxsize = this_img_polar->get_xsize();
733
734 //take out the garbage
735 delete to_polar; to_polar = 0;
736 delete this_img_polar; this_img_polar = 0;
737
738 float *data = cf->get_data();
739 float peak = 0;
740 int peak_index = 0;
741 Util::find_max(data, polarxsize, &peak, &peak_index);
742 delete cf; cf = 0;
743
744 if(peak > best_peak) {
745 best_peak = peak;
746 best_peak_index = peak_index;
747 best_tx = x;
748 best_ty = y;
749 }
750 }
751 }
752
753 float rot_angle = (float) (best_peak_index * 360.0f / polarxsize);
754
755 //return the result
756 Transform tmptt(Dict("type","2d","alpha",0,"tx",-best_tx,"ty",-best_ty));
757 Transform tmprot(Dict("type","2d","alpha",rot_angle,"tx",0,"ty",0));
758 Transform total = tmprot*tmptt;
759 EMData* rotimg=this_img->process("xform",Dict("transform",(Transform*)&total));
760 rotimg->set_attr("xform.align2d",&total);
761
762 return rotimg;
763
764}
765
767 const string & cmp_name, const Dict& cmp_params) const
768{
769
770#ifdef EMAN2_USING_CUDA
771 if(EMData::usecuda == 1) {
772 //if(!this_img->getcudarwdata()) this_img->copy_to_cuda();
773 //if(!to->getcudarwdata()) to->copy_to_cuda();
774 }
775#endif
776
777 // Get the 180 degree ambiguously rotationally aligned and its 180 degree rotation counterpart
778 int zscore = params.set_default("zscore",0);
779 int rfp_mode = params.set_default("rfp_mode",2);
780 EMData *rot_align = RotationalAligner::align_180_ambiguous(this_img,to,rfp_mode,zscore);
781 Transform * tmp = rot_align->get_attr("xform.align2d");
782 Dict rot = tmp->get_rotation("2d");
783 float rotate_angle_solution = rot["alpha"];
784 delete tmp;
785
786 EMData *rot_align_180 = rot_align->process("math.rotate.180");
787
788 Dict trans_params;
789 trans_params["intonly"] = 0;
790 trans_params["maxshift"] = params.set_default("maxshift", -1);
791 trans_params["useflcf"]=params.set_default("useflcf",0);
792
793 // Do the first case translational alignment
794 trans_params["nozero"] = params.set_default("nozero",false);
795 EMData* rot_trans = rot_align->align("translational", to, trans_params, cmp_name, cmp_params);
796 if( rot_align ) { // Clean up
797 delete rot_align;
798 rot_align = 0;
799 }
800
801 // Do the second case translational alignment
802 EMData* rot_180_trans = rot_align_180->align("translational", to, trans_params, cmp_name, cmp_params);
803 if( rot_align_180 ) { // Clean up
804 delete rot_align_180;
805 rot_align_180 = 0;
806 }
807
808 // Finally decide on the result
809 float cmp1 = rot_trans->cmp(cmp_name, to, cmp_params);
810 float cmp2 = rot_180_trans->cmp(cmp_name, to, cmp_params);
811
812 EMData *result = 0;
813 if (cmp1 < cmp2) { // All comparators are defined so default return is "smaller is better"
814 if( rot_180_trans ) {
815 delete rot_180_trans;
816 rot_180_trans = 0;
817 }
818 result = rot_trans;
819 }
820 else {
821 if( rot_trans ) {
822 delete rot_trans;
823 rot_trans = 0;
824 }
825 result = rot_180_trans;
826 rotate_angle_solution -= 180.f;
827 }
828
829 Transform* t = result->get_attr("xform.align2d");
830 t->set_rotation(Dict("type","2d","alpha",rotate_angle_solution));
831 result->set_attr("xform.align2d",t);
832 delete t;
833
834 return result;
835}
836
837EMData *RotateTranslateAlignerBispec::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
838{
839
840 // Get the 180 degree ambiguously rotationally aligned and its 180 degree rotation counterpart
841 int zscore = params.set_default("zscore",0);
842 int rfp = params.set_default("rfpn",4); // rfpn and size were determined empirically. rfpn<4 shows obvious alignment shifts
843 int size = params.set_default("size",16); // size=32 and size=16 seem generally to give equivalent results
844 int harmonic = params.set_default("harmonic",0);
845 EMData *rot_align = this_img->align("rotational_bispec", to,Dict("rfpn",rfp,"size",size,"harmonic",harmonic));
846 Transform * tmp = rot_align->get_attr("xform.align2d");
847 Dict rot = tmp->get_rotation("2d");
848 float rotate_angle_solution = rot["alpha"];
849 delete tmp;
850
851 Dict trans_params;
852 trans_params["intonly"] = false;
853 trans_params["maxshift"] = params.set_default("maxshift", -1);
854 trans_params["useflcf"] = params.set_default("useflcf",0);
855 trans_params["nozero"] = params.set_default("nozero",false);
856 EMData* rot_trans = rot_align->align("translational", to, trans_params, cmp_name, cmp_params);
857 delete rot_align;
858
859 Transform* t = rot_trans->get_attr("xform.align2d");
860 t->set_rotation(Dict("type","2d","alpha",rotate_angle_solution));
861 delete rot_trans;
862
863 EMData *result=this_img->process("xform",Dict("transform",t)); // final object only has one interpolation
864 result->set_attr("xform.align2d",t);
865 delete t;
866
867 return result;
868}
869
870
872 const string & cmp_name, const Dict& cmp_params) const
873{
874
875 //Basically copy params into rotate_translate
876 basealigner_params["maxshift"] = params.set_default("maxshift", -1);
877 basealigner_params["rfp_mode"] = params.set_default("rfp_mode",2);
878 basealigner_params["useflcf"] = params.set_default("useflcf",0);
879
880 //return the correct results
881 return align_using_base(this_img, to, cmp_name, cmp_params);
882
883}
884
885EMData* RotateTranslateFlipAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
886{
887 EMData *flipped = params.set_default("flip", (EMData *) 0);
888 bool delete_flag = false;
889 if (flipped == 0) {
890 flipped = to->process("xform.flip", Dict("axis", "x"));
891 delete_flag = true;
892 }
893
894 EMData *rot_trans_align_flip=0;
895 EMData *rot_trans_align=0;
896 int usebispec=params.set_default("usebispec",0);
897 int useharmonic=params.set_default("useharmonic",0);
898
899
900 if (usebispec || useharmonic) {
901 // Get the non flipped rotational, tranlsationally aligned image
902 Dict rt_params("maxshift", params["maxshift"], "useflcf",params.set_default("useflcf",0),"harmonic",useharmonic);
903 rot_trans_align = this_img->align("rotate_translate_bispec",to,rt_params,cmp_name, cmp_params);
904
905 // Do the same alignment, but using the flipped version of the image
906 rot_trans_align_flip = this_img->align("rotate_translate_bispec", flipped, rt_params, cmp_name, cmp_params);
907 Transform * t = rot_trans_align_flip->get_attr("xform.align2d");
908 t->set_mirror(true);
909 rot_trans_align_flip->set_attr("xform.align2d",t);
910 delete t;
911 }
912 else {
913 // Get the non flipped rotational, tranlsationally aligned image
914 Dict rt_params("maxshift", params["maxshift"], "rfp_mode", params.set_default("rfp_mode",2),"useflcf",params.set_default("useflcf",0),"zscore",params.set_default("zscore",0));
915 rot_trans_align = this_img->align("rotate_translate",to,rt_params,cmp_name, cmp_params);
916
917 // Do the same alignment, but using the flipped version of the image
918 rot_trans_align_flip = this_img->align("rotate_translate", flipped, rt_params, cmp_name, cmp_params);
919 Transform * t = rot_trans_align_flip->get_attr("xform.align2d");
920 t->set_mirror(true);
921 rot_trans_align_flip->set_attr("xform.align2d",t);
922 delete t;
923 }
924
925 // Now finally decide on what is the best answer
926 float cmp1 = rot_trans_align->cmp(cmp_name, to, cmp_params);
927 float cmp2 = rot_trans_align_flip->cmp(cmp_name, flipped, cmp_params);
928
929 if (delete_flag){
930 if(flipped) {
931 delete flipped;
932 flipped = 0;
933 }
934 }
935
936 EMData *result = 0;
937 if (cmp1 < cmp2 ) {
938
939 if( rot_trans_align_flip ) {
940 delete rot_trans_align_flip;
941 rot_trans_align_flip = 0;
942 }
943 result = rot_trans_align;
944 }
945 else {
946 if( rot_trans_align ) {
947 delete rot_trans_align;
948 rot_trans_align = 0;
949 }
950 result = rot_trans_align_flip;
951 result->process_inplace("xform.flip",Dict("axis","x"));
952 }
953
954 return result;
955}
956
958 const string & cmp_name, const Dict& cmp_params) const
959{
960
961 //Basically copy params into rotate_translate
962 basealigner_params["flip"] = params.set_default("flip", (EMData *) 0);
963 basealigner_params["maxshift"] = params.set_default("maxshift", -1);
964 basealigner_params["rfp_mode"] = params.set_default("rfp_mode",2);
965 basealigner_params["useflcf"] = params.set_default("useflcf",0);
966 basealigner_params["zscore"] = params.set_default("zscore",0);
967
968 //return the correct results
969 return align_using_base(this_img, to, cmp_name, cmp_params);
970
971}
972
974 const string & cmp_name, const Dict& cmp_params) const
975{
976 // Get the non flipped rotational, tranlsationally aligned image
977 Dict rt_params("maxshift", params["maxshift"],"r1",params.set_default("r1",-1),"r2",params.set_default("r2",-1));
978 EMData *rot_trans_align = this_img->align("rotate_translate_iterative",to,rt_params,cmp_name, cmp_params);
979
980 // Do the same alignment, but using the flipped version of the image
981 EMData *flipped = params.set_default("flip", (EMData *) 0);
982 bool delete_flag = false;
983 if (flipped == 0) {
984 flipped = to->process("xform.flip", Dict("axis", "x"));
985 delete_flag = true;
986 }
987
988 EMData * rot_trans_align_flip = this_img->align("rotate_translate_iterative", flipped, rt_params, cmp_name, cmp_params);
989 Transform* t = rot_trans_align_flip->get_attr("xform.align2d");
990 t->set_mirror(true);
991 rot_trans_align_flip->set_attr("xform.align2d",t);
992 delete t;
993
994 // Now finally decide on what is the best answer
995 float cmp1 = rot_trans_align->cmp(cmp_name, to, cmp_params);
996 float cmp2 = rot_trans_align_flip->cmp(cmp_name, flipped, cmp_params);
997
998 if (delete_flag){
999 if(flipped) {
1000 delete flipped;
1001 flipped = 0;
1002 }
1003 }
1004
1005 EMData *result = 0;
1006 if (cmp1 < cmp2 ) {
1007
1008 if( rot_trans_align_flip ) {
1009 delete rot_trans_align_flip;
1010 rot_trans_align_flip = 0;
1011 }
1012 result = rot_trans_align;
1013 }
1014 else {
1015 if( rot_trans_align ) {
1016 delete rot_trans_align;
1017 rot_trans_align = 0;
1018 }
1019 result = rot_trans_align_flip;
1020 result->process_inplace("xform.flip",Dict("axis","x"));
1021 }
1022
1023 return result;
1024}
1025
1027 const string & cmp_name, const Dict& cmp_params) const
1028{
1029
1030 //Basically copy params into rotate_translate
1031 basealigner_params["flip"] = params.set_default("flip", (EMData *) 0);
1032 basealigner_params["maxshift"] = params.set_default("maxshift", -1);
1033 basealigner_params["r1"] = params.set_default("r1",-1);
1034 basealigner_params["r2"] = params.set_default("r2",-1);
1035 basealigner_params["maxiter"] = params.set_default("maxiter",3);
1036
1037 //return the correct results
1038 return align_using_base(this_img, to, cmp_name, cmp_params);
1039
1040}
1041
1043 const string & cmp_name, const Dict& cmp_params) const
1044{
1045 if (cmp_name != "dot" && cmp_name != "ccc") throw InvalidParameterException("Resample aligner only works for dot and ccc");
1046
1047 int maxtx = params.set_default("tx", 0);
1048 int maxty = params.set_default("ty", 0);
1049 int r1 = params.set_default("r1",-1);
1050 int r2 = params.set_default("r2",-1);
1051
1052 if(this_img->get_xsize()/2 - 1 - r2 - maxtx <= 0 || (r2 == -1 && maxtx > 0)){
1053 cout << "\nRunTimeError: nx/2 - 1 - r2 - tx must be greater than or = 0\n" << endl; // For some reason the expection message is not being print, stupid C++
1054 throw InvalidParameterException("nx/2 - 1 - r2 - tx must be greater than or = 0");
1055 }
1056 if(this_img->get_ysize()/2 - 1 - r2 - maxty <= 0 || (r2 == -1 && maxty > 0)){
1057 cout << "\nRunTimeError:ny/2 - 1 - r2 - ty must be greater than or = 0\n" << endl; // For some reason the expection message is not being print, stupid C++
1058 throw InvalidParameterException("ny/2 - 1 - r2 - ty must be greater than or = 0");
1059 }
1060
1061// float best_peak = -numeric_limits<float>::infinity();
1062 float best_peak = -1.0e37;
1063 int best_peak_index = 0;
1064 int best_tx = 0;
1065 int best_ty = 0;
1066 int polarxsize = 0;
1067 bool flip = false;
1068
1069 for(int x = -maxtx; x <= maxtx; x++){
1070 for(int y = -maxty; y <= maxty; y++){
1071
1072 EMData * to_polar = to->unwrap(r1,r2,-1,0,0,true);
1073 EMData * this_img_polar = this_img->unwrap(r1,r2,-1,x,y,true);
1074 EMData * cfflip = this_img_polar->calc_ccfx(to_polar, 0, this_img_polar->get_ysize(), false, true);
1075 EMData * cf = this_img_polar->calc_ccfx(to_polar, 0, this_img_polar->get_ysize());
1076
1077 polarxsize = this_img_polar->get_xsize();
1078
1079 //take out the garbage
1080 delete to_polar; to_polar = 0;
1081 delete this_img_polar; this_img_polar = 0;
1082
1083 float *data = cf->get_data();
1084 float peak = 0;
1085 int peak_index = 0;
1086 Util::find_max(data, polarxsize, &peak, &peak_index);
1087 delete cf; cf = 0;
1088
1089 if(peak > best_peak) {
1090 best_peak = peak;
1091 best_peak_index = peak_index;
1092 best_tx = x;
1093 best_ty = y;
1094 flip = false;
1095 }
1096
1097 data = cfflip->get_data();
1098 Util::find_max(data, polarxsize, &peak, &peak_index);
1099 delete cfflip; cfflip = 0;
1100
1101 if(peak > best_peak) {
1102 best_peak = peak;
1103 best_peak_index = peak_index;
1104 best_tx = x;
1105 best_ty = y;
1106 flip = true;
1107 }
1108 }
1109 }
1110
1111 float rot_angle = (float) (best_peak_index * 360.0f / polarxsize);
1112
1113 //return the result
1114 Transform tmptt(Dict("type","2d","alpha",0,"tx",-best_tx,"ty",-best_ty));
1115 Transform tmprot(Dict("type","2d","alpha",rot_angle,"tx",0,"ty",0));
1116 Transform total = tmprot*tmptt;
1117 EMData * rotimg=this_img->process("xform",Dict("transform",(Transform*)&total));
1118 rotimg->set_attr("xform.align2d",&total);
1119 if(flip == true) {
1120 rotimg->process_inplace("xform.flip",Dict("axis", "x"));
1121 }
1122
1123 return rotimg;
1124
1125}
1126
1128 const string& cmp_name, const Dict& cmp_params) const
1129{
1130 Dict rot_params("rfp_mode",params.set_default("rfp_mode",2));
1131 EMData *r1 = this_img->align("rotational", to, rot_params,cmp_name, cmp_params);
1132
1133
1134 EMData* flipped =to->process("xform.flip", Dict("axis", "x"));
1135 EMData *r2 = this_img->align("rotational", flipped,rot_params, cmp_name, cmp_params);
1136 Transform* t = r2->get_attr("xform.align2d");
1137 t->set_mirror(true);
1138 r2->set_attr("xform.align2d",t);
1139 delete t;
1140
1141 float cmp1 = r1->cmp(cmp_name, to, cmp_params);
1142 float cmp2 = r2->cmp(cmp_name, flipped, cmp_params);
1143
1144 delete flipped; flipped = 0;
1145
1146 EMData *result = 0;
1147
1148 if (cmp1 < cmp2) {
1149 if( r2 )
1150 {
1151 delete r2;
1152 r2 = 0;
1153 }
1154 result = r1;
1155 }
1156 else {
1157 if( r1 )
1158 {
1159 delete r1;
1160 r1 = 0;
1161 }
1162 result = r2;
1163 result->process_inplace("xform.flip",Dict("axis","x"));
1164 }
1165
1166 return result;
1167}
1168
1170 const string& cmp_name, const Dict& cmp_params) const
1171{
1172 Dict rot_params("r1",params.set_default("r1",-1),"r2",params.set_default("r2",-1));
1173 EMData * r1 = this_img->align("rotational_iterative", to, rot_params,cmp_name, cmp_params);
1174
1175 EMData * flipped =to->process("xform.flip", Dict("axis", "x"));
1176 EMData * r2 = this_img->align("rotational_iterative", flipped,rot_params, cmp_name, cmp_params);
1177 Transform* t = r2->get_attr("xform.align2d");
1178 t->set_mirror(true);
1179 r2->set_attr("xform.align2d",t);
1180 delete t;
1181
1182 float cmp1 = r1->cmp(cmp_name, to, cmp_params);
1183 float cmp2 = r2->cmp(cmp_name, flipped, cmp_params);
1184
1185 delete flipped; flipped = 0;
1186
1187 EMData *result = 0;
1188
1189 if (cmp1 < cmp2) {
1190 if( r2 )
1191 {
1192 delete r2;
1193 r2 = 0;
1194 }
1195 result = r1;
1196 }
1197 else {
1198 if( r1 )
1199 {
1200 delete r1;
1201 r1 = 0;
1202 }
1203 result = r2;
1204 result->process_inplace("xform.flip",Dict("axis","x"));
1205 }
1206
1207 return result;
1208}
1209
1210// David Woolford says FIXME
1211// You will note the excessive amount of EMData copying that's going in this function
1212// This is because functions that are operating on the EMData objects are changing them
1213// and if you do not use copies the whole algorithm breaks. I did not have time to go
1214// through and rectify this situation.
1215// David Woolford says - this problem is related to the fact that many functions that
1216// take EMData pointers as arguments do not take them as constant pointers to constant
1217// objects, instead they are treated as raw (completely changeable) pointers. This means
1218// it's hard to track down which functions are changing the EMData objects, because they
1219// all do (in name). If this behavior is unavoidable then ignore this comment, however if possible it would
1220// be good to make things const as much as possible. For example in alignment, technically
1221// the argument EMData objects (raw pointers) should not be altered... should they?
1222//
1223// But const can be very annoying sometimes...
1225 const string & cmp_name, const Dict& cmp_params) const
1226{
1227 EMData *flip = params.set_default("flip", (EMData *) 0);
1228 int maxshift = params.set_default("maxshift", this_img->get_xsize()/8);
1229 if (maxshift < 2) throw InvalidParameterException("maxshift must be greater than or equal to 2");
1230
1231 int ny = this_img->get_ysize();
1232 int xst = (int) floor(2 * M_PI * ny);
1233 xst = Util::calc_best_fft_size(xst);
1234
1235 Dict d("n",2);
1236 EMData *to_shrunk_unwrapped = to->process("math.medianshrink",d);
1237
1238 int to_copy_r2 = to_shrunk_unwrapped->get_ysize() / 2 - 2 - maxshift / 2;
1239 EMData *tmp = to_shrunk_unwrapped->unwrap(4, to_copy_r2, xst / 2, 0, 0, true);
1240 if( to_shrunk_unwrapped )
1241 {
1242 delete to_shrunk_unwrapped;
1243 to_shrunk_unwrapped = 0;
1244 }
1245 to_shrunk_unwrapped = tmp;
1246
1247 EMData *to_shrunk_unwrapped_copy = to_shrunk_unwrapped->copy();
1248 EMData* to_unwrapped = to->unwrap(4, to->get_ysize() / 2 - 2 - maxshift, xst, 0, 0, true);
1249 EMData *to_unwrapped_copy = to_unwrapped->copy();
1250
1251 bool delete_flipped = true;
1252 EMData *flipped = 0;
1253 if (flip) {
1254 delete_flipped = false;
1255 flipped = flip;
1256 }
1257 else {
1258 flipped = to->process("xform.flip", Dict("axis", "x"));
1259 }
1260 EMData *to_shrunk_flipped_unwrapped = flipped->process("math.medianshrink",d);
1261 tmp = to_shrunk_flipped_unwrapped->unwrap(4, to_copy_r2, xst / 2, 0, 0, true);
1262 if( to_shrunk_flipped_unwrapped )
1263 {
1264 delete to_shrunk_flipped_unwrapped;
1265 to_shrunk_flipped_unwrapped = 0;
1266 }
1267 to_shrunk_flipped_unwrapped = tmp;
1268 EMData *to_shrunk_flipped_unwrapped_copy = to_shrunk_flipped_unwrapped->copy();
1269 EMData* to_flip_unwrapped = flipped->unwrap(4, to->get_ysize() / 2 - 2 - maxshift, xst, 0, 0, true);
1270 EMData* to_flip_unwrapped_copy = to_flip_unwrapped->copy();
1271
1272 if (delete_flipped && flipped != 0) {
1273 delete flipped;
1274 flipped = 0;
1275 }
1276
1277 EMData *this_shrunk_2 = this_img->process("math.medianshrink",d);
1278
1279 float bestval = FLT_MAX;
1280 float bestang = 0;
1281 int bestflip = 0;
1282 float bestdx = 0;
1283 float bestdy = 0;
1284
1285 int half_maxshift = maxshift / 2;
1286
1287 int ur2 = this_shrunk_2->get_ysize() / 2 - 2 - half_maxshift;
1288 for (int dy = -half_maxshift; dy <= half_maxshift; dy += 1) {
1289 for (int dx = -half_maxshift; dx <= half_maxshift; dx += 1) {
1290 if (hypot(dx, dy) <= half_maxshift) {
1291 EMData *uw = this_shrunk_2->unwrap(4, ur2, xst / 2, dx, dy, true);
1292 EMData *uwc = uw->copy();
1293 EMData *a = uw->calc_ccfx(to_shrunk_unwrapped);
1294
1295 uwc->rotate_x(a->calc_max_index());
1296 float cm = uwc->cmp(cmp_name, to_shrunk_unwrapped_copy, cmp_params);
1297 if (cm < bestval) {
1298 bestval = cm;
1299 bestang = (float) (2.0 * M_PI * a->calc_max_index() / a->get_xsize());
1300 bestdx = (float)dx;
1301 bestdy = (float)dy;
1302 bestflip = 0;
1303 }
1304
1305
1306 if( a )
1307 {
1308 delete a;
1309 a = 0;
1310 }
1311 if( uw )
1312 {
1313 delete uw;
1314 uw = 0;
1315 }
1316 if( uwc )
1317 {
1318 delete uwc;
1319 uwc = 0;
1320 }
1321 uw = this_shrunk_2->unwrap(4, ur2, xst / 2, dx, dy, true);
1322 uwc = uw->copy();
1323 a = uw->calc_ccfx(to_shrunk_flipped_unwrapped);
1324
1325 uwc->rotate_x(a->calc_max_index());
1326 cm = uwc->cmp(cmp_name, to_shrunk_flipped_unwrapped_copy, cmp_params);
1327 if (cm < bestval) {
1328 bestval = cm;
1329 bestang = (float) (2.0 * M_PI * a->calc_max_index() / a->get_xsize());
1330 bestdx = (float)dx;
1331 bestdy = (float)dy;
1332 bestflip = 1;
1333 }
1334
1335 if( a )
1336 {
1337 delete a;
1338 a = 0;
1339 }
1340
1341 if( uw )
1342 {
1343 delete uw;
1344 uw = 0;
1345 }
1346 if( uwc )
1347 {
1348 delete uwc;
1349 uwc = 0;
1350 }
1351 }
1352 }
1353 }
1354 if( this_shrunk_2 )
1355 {
1356 delete this_shrunk_2;
1357 this_shrunk_2 = 0;
1358 }
1359 if( to_shrunk_unwrapped )
1360 {
1361 delete to_shrunk_unwrapped;
1362 to_shrunk_unwrapped = 0;
1363 }
1364 if( to_shrunk_unwrapped_copy )
1365 {
1366 delete to_shrunk_unwrapped_copy;
1367 to_shrunk_unwrapped_copy = 0;
1368 }
1369 if( to_shrunk_flipped_unwrapped )
1370 {
1371 delete to_shrunk_flipped_unwrapped;
1372 to_shrunk_flipped_unwrapped = 0;
1373 }
1374 if( to_shrunk_flipped_unwrapped_copy )
1375 {
1376 delete to_shrunk_flipped_unwrapped_copy;
1377 to_shrunk_flipped_unwrapped_copy = 0;
1378 }
1379 bestdx *= 2;
1380 bestdy *= 2;
1381 bestval = FLT_MAX;
1382
1383 float bestdx2 = bestdx;
1384 float bestdy2 = bestdy;
1385 // Note I tried steps less than 1.0 (sub pixel precision) and it actually appeared detrimental
1386 // So my advice is to stick with dx += 1.0 etc unless you really are looking to fine tune this
1387 // algorithm
1388 for (float dy = bestdy2 - 3; dy <= bestdy2 + 3; dy += 1.0 ) {
1389 for (float dx = bestdx2 - 3; dx <= bestdx2 + 3; dx += 1.0 ) {
1390
1391 if (hypot(dx, dy) <= maxshift) {
1392 EMData *uw = this_img->unwrap(4, this_img->get_ysize() / 2 - 2 - maxshift, xst, (int)dx, (int)dy, true);
1393 EMData *uwc = uw->copy();
1394 EMData *a = uw->calc_ccfx(to_unwrapped);
1395
1396 uwc->rotate_x(a->calc_max_index());
1397 float cm = uwc->cmp(cmp_name, to_unwrapped_copy, cmp_params);
1398
1399 if (cm < bestval) {
1400 bestval = cm;
1401 bestang = (float)(2.0 * M_PI * a->calc_max_index() / a->get_xsize());
1402 bestdx = dx;
1403 bestdy = dy;
1404 bestflip = 0;
1405 }
1406
1407 if( a )
1408 {
1409 delete a;
1410 a = 0;
1411 }
1412 if( uw )
1413 {
1414 delete uw;
1415 uw = 0;
1416 }
1417 if( uwc )
1418 {
1419 delete uwc;
1420 uwc = 0;
1421 }
1422 uw = this_img->unwrap(4, this_img->get_ysize() / 2 - 2 - maxshift, xst, (int)dx, (int)dy, true);
1423 uwc = uw->copy();
1424 a = uw->calc_ccfx(to_flip_unwrapped);
1425
1426 uwc->rotate_x(a->calc_max_index());
1427 cm = uwc->cmp(cmp_name, to_flip_unwrapped_copy, cmp_params);
1428
1429 if (cm < bestval) {
1430 bestval = cm;
1431 bestang = (float)(2.0 * M_PI * a->calc_max_index() / a->get_xsize());
1432 bestdx = dx;
1433 bestdy = dy;
1434 bestflip = 1;
1435 }
1436
1437 if( a )
1438 {
1439 delete a;
1440 a = 0;
1441 }
1442 if( uw )
1443 {
1444 delete uw;
1445 uw = 0;
1446 }
1447 if( uwc )
1448 {
1449 delete uwc;
1450 uwc = 0;
1451 }
1452 }
1453 }
1454 }
1455 if( to_unwrapped ) {delete to_unwrapped;to_unwrapped = 0;}
1456 if( to_shrunk_unwrapped ) { delete to_shrunk_unwrapped; to_shrunk_unwrapped = 0;}
1457 if (to_unwrapped_copy) { delete to_unwrapped_copy; to_unwrapped_copy = 0; }
1458 if (to_flip_unwrapped) { delete to_flip_unwrapped; to_flip_unwrapped = 0; }
1459 if (to_flip_unwrapped_copy) { delete to_flip_unwrapped_copy; to_flip_unwrapped_copy = 0;}
1460
1461 bestang *= (float)EMConsts::rad2deg;
1462 Transform t(Dict("type","2d","alpha",(float)bestang));
1463 t.set_pre_trans(Vec2f(-bestdx,-bestdy));
1464 if (bestflip) {
1465 t.set_mirror(true);
1466 }
1467
1468 EMData* ret = this_img->process("xform",Dict("transform",&t));
1469 ret->set_attr("xform.align2d",&t);
1470
1471 return ret;
1472}
1473
1474
1476 const string & cmp_name, const Dict& cmp_params) const
1477{
1478
1479 EMData *flip = params.set_default("flip", (EMData *) 0);
1480 int maxshift = params.set_default("maxshift", -1);
1481
1482 EMData *flipped = 0;
1483
1484 bool delete_flipped = true;
1485 if (flip) {
1486 delete_flipped = false;
1487 flipped = flip;
1488 }
1489 else {
1490 flipped = to->process("xform.flip", Dict("axis", "x"));
1491 }
1492
1493 int nx = this_img->get_xsize();
1494
1495 if (maxshift < 0) {
1496 maxshift = nx / 10;
1497 }
1498
1499 float angle_step = params.set_default("angstep", 0.0f);
1500 if ( angle_step == 0 ) angle_step = atan2(2.0f, (float)nx);
1501 else {
1502 angle_step *= (float)EMConsts::deg2rad; //convert to radians
1503 }
1504 float trans_step = params.set_default("transtep",1.0f);
1505
1506 if (trans_step <= 0) throw InvalidParameterException("transstep must be greater than 0");
1507 if (angle_step <= 0) throw InvalidParameterException("angstep must be greater than 0");
1508
1509
1510 Dict shrinkfactor("n",2);
1511 EMData *this_img_shrink = this_img->process("math.medianshrink",shrinkfactor);
1512 EMData *to_shrunk = to->process("math.medianshrink",shrinkfactor);
1513 EMData *flipped_shrunk = flipped->process("math.medianshrink",shrinkfactor);
1514
1515 int bestflip = 0;
1516 float bestdx = 0;
1517 float bestdy = 0;
1518
1519 float bestang = 0;
1520 float bestval = FLT_MAX;
1521
1522 int half_maxshift = maxshift / 2;
1523
1524
1525 for (int dy = -half_maxshift; dy <= half_maxshift; ++dy) {
1526 for (int dx = -half_maxshift; dx <= half_maxshift; ++dx) {
1527 if (hypot(dx, dy) <= maxshift) {
1528 for (float ang = -angle_step * 2.0f; ang <= (float)2 * M_PI; ang += angle_step * 4.0f) {
1529 EMData v(*this_img_shrink);
1530 Transform t(Dict("type","2d","alpha",static_cast<float>(ang*EMConsts::rad2deg)));
1531 t.set_trans((float)dx,(float)dy);
1532 v.transform(t);
1533// v.rotate_translate(ang*EMConsts::rad2deg, 0.0f, 0.0f, (float)dx, (float)dy, 0.0f);
1534
1535 float lc = v.cmp(cmp_name, to_shrunk, cmp_params);
1536
1537 if (lc < bestval) {
1538 bestval = lc;
1539 bestang = ang;
1540 bestdx = (float)dx;
1541 bestdy = (float)dy;
1542 bestflip = 0;
1543 }
1544
1545 lc = v.cmp(cmp_name,flipped_shrunk , cmp_params);
1546 if (lc < bestval) {
1547 bestval = lc;
1548 bestang = ang;
1549 bestdx = (float)dx;
1550 bestdy = (float)dy;
1551 bestflip = 1;
1552 }
1553 }
1554 }
1555 }
1556 }
1557
1558 if( to_shrunk )
1559 {
1560 delete to_shrunk;
1561 to_shrunk = 0;
1562 }
1563 if( flipped_shrunk )
1564 {
1565 delete flipped_shrunk;
1566 flipped_shrunk = 0;
1567 }
1568 if( this_img_shrink )
1569 {
1570 delete this_img_shrink;
1571 this_img_shrink = 0;
1572 }
1573
1574 bestdx *= 2;
1575 bestdy *= 2;
1576 bestval = FLT_MAX;
1577
1578 float bestdx2 = bestdx;
1579 float bestdy2 = bestdy;
1580 float bestang2 = bestang;
1581
1582 for (float dy = bestdy2 - 3; dy <= bestdy2 + 3; dy += trans_step) {
1583 for (float dx = bestdx2 - 3; dx <= bestdx2 + 3; dx += trans_step) {
1584 if (hypot(dx, dy) <= maxshift) {
1585 for (float ang = bestang2 - angle_step * 6.0f; ang <= bestang2 + angle_step * 6.0f; ang += angle_step) {
1586 EMData v(*this_img);
1587 Transform t(Dict("type","2d","alpha",static_cast<float>(ang*EMConsts::rad2deg)));
1588 t.set_trans(dx,dy);
1589 v.transform(t);
1590// v.rotate_translate(ang*EMConsts::rad2deg, 0.0f, 0.0f, (float)dx, (float)dy, 0.0f);
1591
1592 float lc = v.cmp(cmp_name, to, cmp_params);
1593
1594 if (lc < bestval) {
1595 bestval = lc;
1596 bestang = ang;
1597 bestdx = dx;
1598 bestdy = dy;
1599 bestflip = 0;
1600 }
1601
1602 lc = v.cmp(cmp_name, flipped, cmp_params);
1603
1604 if (lc < bestval) {
1605 bestval = lc;
1606 bestang = ang;
1607 bestdx = dx;
1608 bestdy = dy;
1609 bestflip = 1;
1610 }
1611 }
1612 }
1613 }
1614 }
1615
1616 if (delete_flipped) { delete flipped; flipped = 0; }
1617
1618 bestang *= (float)EMConsts::rad2deg;
1619 Transform t(Dict("type","2d","alpha",(float)bestang));
1620 t.set_trans(bestdx,bestdy);
1621
1622 if (bestflip) {
1623 t.set_mirror(true);
1624 }
1625
1626 EMData* rslt = this_img->process("xform",Dict("transform",&t));
1627 rslt->set_attr("xform.align2d",&t);
1628
1629 return rslt;
1630}
1631
1632EMData* SymAlignProcessor::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
1633{
1634
1635 // Set parms
1636 float dphi = params.set_default("dphi",10.f);
1637 float lphi = params.set_default("lphi",0.0f);
1638 float uphi = params.set_default("uphi",359.9f);
1639
1640 Dict d;
1641 d["inc_mirror"] = true;
1642 d["delta"] = params.set_default("delta",10.f);
1643
1644 //Genrate points on a sphere in an asymmetric unit
1645 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
1646 vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
1647
1648 //Genrate symmetry related orritenations
1649 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
1650
1651 float bestquality = 0.0f;
1652 EMData* bestimage = 0;
1653 for(vector<Transform>::const_iterator trans_it = transforms.begin(); trans_it != transforms.end(); trans_it++) {
1654 Dict tparams = trans_it->get_params("eman");
1655 Transform t(tparams);
1656 for( float phi = lphi; phi < uphi; phi += dphi ) {
1657 tparams["phi"] = phi;
1658 t.set_rotation(tparams);
1659
1660 //Get the averagaer
1661 Averager* imgavg = Factory<Averager>::get((string)params.set_default("avger","mean"));
1662 //Now make the averages
1663 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
1664 Transform sympos = (*it)*t;
1665 EMData* transformed = this_img->process("xform",Dict("transform",&sympos));
1666 imgavg->add_image(transformed);
1667 delete transformed;
1668 }
1669
1670 EMData* symptcl=imgavg->finish();
1671 delete imgavg;
1672 //See which average is the best
1673 float quality = symptcl->get_attr("sigma");
1674 cout << quality << " " << phi << endl;
1675 if(quality > bestquality) {
1676 bestquality = quality;
1677 bestimage = symptcl;
1678 } else {
1679 delete symptcl;
1680 }
1681 }
1682 }
1683 if(sym != 0) delete sym;
1684
1685 return bestimage;
1686}
1687
1688static double refalifn(const gsl_vector * v, void *params)
1689{
1690 Dict *dict = (Dict *) params;
1691
1692 double x = gsl_vector_get(v, 0);
1693 double y = gsl_vector_get(v, 1);
1694 double a = gsl_vector_get(v, 2);
1695
1696 EMData *this_img = (*dict)["this"];
1697 EMData *with = (*dict)["with"];
1698 bool mirror = (*dict)["mirror"];
1699
1700 Transform t(Dict("type","2d","alpha",static_cast<float>(a)));
1701 t.set_trans((float)x,(float)y);
1702 t.set_mirror(mirror);
1703 if (v->size>3) {
1704 float sca=(float)gsl_vector_get(v, 3);
1705 if (sca<.7 || sca>1.3) return 1.0e20;
1706 t.set_scale((float)gsl_vector_get(v, 3));
1707 }
1708 EMData *tmp = this_img->process("xform",Dict("transform",&t));
1709 if (dict->has_key("mask")) tmp->mult(*(EMData *)((*dict)["mask"]));
1710
1711// printf("GSL %f %f %f %d %f\n",x,y,a,mirror,(float)gsl_vector_get(v, 3));
1712 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
1713 double result = c->cmp(tmp,with);
1714
1715 if (tmp != 0) delete tmp;
1716
1717 return result;
1718}
1719
1720static void refalidf(const gsl_vector * v, void *params,gsl_vector * df) {
1721 // we do this using a simple local difference estimate due to the expense of the calculation.
1722 // The step has to be large enough for the similarity metric
1723 // To provide an accurate change in value.
1724 static double lstep[4] = { 0.05, 0.05, 0.1, 0.01 };
1725
1726 gsl_vector *vc = gsl_vector_alloc(v->size);
1727 gsl_vector_memcpy(vc,v);
1728
1729 double f = refalifn(v,params);
1730 for (unsigned int i=0; i<v->size; i++) {
1731 // double *vp = gsl_vector_ptr(vc,i);
1732 double vp = gsl_vector_get(vc,i);
1733 // *vp+=lstep[i];
1734 gsl_vector_set(vc,i,vp+lstep[i]);
1735 double f2 = refalifn(vc,params);
1736 // *vp-=lstep[i];
1737 gsl_vector_set(vc,i,vp);
1738
1739 gsl_vector_set(df,i,(f2-f)/lstep[i]);
1740 }
1741
1742 gsl_vector_free(vc);
1743 return;
1744}
1745
1746static void refalifdf(const gsl_vector * v, void *params, double * f, gsl_vector * df) {
1747 // we do this using a simple local difference estimate due to the expense of the calculation.
1748 // The step has to be large enough for the similarity metric
1749 // To provide an accurate change in value.
1750 static double lstep[4] = { 0.05, 0.05, 0.1, 0.01 };
1751
1752 gsl_vector *vc = gsl_vector_alloc(v->size);
1753 gsl_vector_memcpy(vc,v);
1754
1755 *f = refalifn(v,params);
1756 for (unsigned int i=0; i<v->size; i++) {
1757 // double *vp = gsl_vector_ptr(vc,i);
1758 double vp = gsl_vector_get(vc,i);
1759 // *vp+=lstep[i];
1760 gsl_vector_set(vc,i,vp+lstep[i]);
1761 double f2 = refalifn(vc,params);
1762 // *vp-=lstep[i];
1763 gsl_vector_set(vc,i,vp);
1764
1765 gsl_vector_set(df,i,(f2-*f)/lstep[i]);
1766 }
1767
1768 gsl_vector_free(vc);
1769 return;
1770}
1771
1772static double refalifnfast(const gsl_vector * v, void *params)
1773{
1774 Dict *dict = (Dict *) params;
1775 EMData *this_img = (*dict)["this"];
1776 EMData *img_to = (*dict)["with"];
1777 bool mirror = (*dict)["mirror"];
1778
1779 double x = gsl_vector_get(v, 0);
1780 double y = gsl_vector_get(v, 1);
1781 double a = gsl_vector_get(v, 2);
1782
1783 double r = this_img->dot_rotate_translate(img_to, (float)x, (float)y, (float)a, mirror);
1784 int nsec = this_img->get_xsize() * this_img->get_ysize();
1785 double result = 1.0 - r / nsec;
1786
1787// cout << result << " x " << x << " y " << y << " az " << a << endl;
1788 return result;
1789}
1790
1791EMData *RefineAligner::align(EMData * this_img, EMData *to,
1792 const string & cmp_name, const Dict& cmp_params) const
1793{
1794
1795 if (!to) {
1796 return 0;
1797 }
1798
1799 EMData *result;
1800 int mode = params.set_default("mode", 0);
1801 float saz = 0.0;
1802 float sdx = 0.0;
1803 float sdy = 0.0;
1804 float sscale = 1.0;
1805 bool mirror = false;
1806 Transform* t;
1807 if (params.has_key("xform.align2d") ) {
1808 t = params["xform.align2d"];
1809 Dict params = t->get_params("2d");
1810 saz = params["alpha"];
1811 sdx = params["tx"];
1812 sdy = params["ty"];
1813 mirror = params["mirror"];
1814 sscale = params["scale"];
1815 } else {
1816 t = new Transform(); // is the identity
1817 }
1818
1819 // We do this to prevent the GSL routine from crashing on an invalid alignment
1820 if ((float)(this_img->get_attr("sigma"))==0.0 || (float)(to->get_attr("sigma"))==0.0) {
1821 result = this_img->process("xform",Dict("transform",t));
1822 result->set_attr("xform.align2d",t);
1823 delete t;
1824 return result;
1825 }
1826
1827 float stepx = params.set_default("stepx",1.0f);
1828 float stepy = params.set_default("stepy",1.0f);
1829 // Default step is 5 degree - note in EMAN1 it was 0.1 radians
1830 float stepaz = params.set_default("stepaz",5.0f);
1831 float stepscale = params.set_default("stepscale",0.0f);
1832
1833 int np = 3;
1834 if (stepscale!=0.0) np++;
1835 Dict gsl_params;
1836 gsl_params["this"] = this_img;
1837 gsl_params["with"] = to;
1838 gsl_params["snr"] = params["snr"];
1839 gsl_params["mirror"] = mirror;
1840 if (params.has_key("mask")) gsl_params["mask"]=params["mask"];
1841
1842 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
1843 gsl_vector *ss = gsl_vector_alloc(np);
1844
1845
1846 gsl_vector_set(ss, 0, stepx);
1847 gsl_vector_set(ss, 1, stepy);
1848 gsl_vector_set(ss, 2, stepaz);
1849 if (stepscale!=0.0) gsl_vector_set(ss,3,stepscale);
1850
1851 gsl_vector *x = gsl_vector_alloc(np);
1852 gsl_vector_set(x, 0, sdx);
1853 gsl_vector_set(x, 1, sdy);
1854 gsl_vector_set(x, 2, saz);
1855 if (stepscale!=0.0) gsl_vector_set(x,3,1.0);
1856
1857 Cmp *c = 0;
1858
1859 gsl_multimin_function minex_func;
1860 if (mode == 2) {
1861 minex_func.f = &refalifnfast;
1862 }
1863 else {
1864 c = Factory < Cmp >::get(cmp_name, cmp_params);
1865 gsl_params["cmp"] = (void *) c;
1866 minex_func.f = &refalifn;
1867 }
1868
1869 minex_func.n = np;
1870 minex_func.params = (void *) &gsl_params;
1871
1872 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
1873 gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
1874
1875 int rval = GSL_CONTINUE;
1876 int status = GSL_SUCCESS;
1877 int iter = 1;
1878
1879 float precision = params.set_default("precision",0.04f);
1880 int maxiter = params.set_default("maxiter",28);
1881
1882// printf("Refine sx=%1.2f sy=%1.2f sa=%1.2f prec=%1.4f maxit=%d\n",stepx,stepy,stepaz,precision,maxiter);
1883// printf("%1.2f %1.2f %1.1f ->",(float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1),(float)gsl_vector_get(s->x, 2));
1884
1885 while (rval == GSL_CONTINUE && iter < maxiter) {
1886 iter++;
1887 status = gsl_multimin_fminimizer_iterate(s);
1888 if (status) {
1889 break;
1890 }
1891 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
1892 }
1893
1894 int maxshift = params.set_default("maxshift",-1);
1895
1896 if (maxshift <= 0) {
1897 maxshift = this_img->get_xsize() / 4;
1898 }
1899 float fmaxshift = static_cast<float>(maxshift);
1900 if ( fmaxshift >= fabs((float)gsl_vector_get(s->x, 0)) && fmaxshift >= fabs((float)gsl_vector_get(s->x, 1)) && (stepscale==0 || (((float)gsl_vector_get(s->x, 3))<1.3 && ((float)gsl_vector_get(s->x, 3))<0.7)) )
1901 {
1902// printf(" Refine good %1.2f %1.2f %1.1f\n",(float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1),(float)gsl_vector_get(s->x, 2));
1903 Transform tsoln(Dict("type","2d","alpha",(float)gsl_vector_get(s->x, 2)));
1904 tsoln.set_mirror(mirror);
1905 tsoln.set_trans((float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1));
1906 if (stepscale!=0.0) tsoln.set_scale((float)gsl_vector_get(s->x, 3));
1907 result = this_img->process("xform",Dict("transform",&tsoln));
1908 result->set_attr("xform.align2d",&tsoln);
1909 } else { // The refine aligner failed - this shift went beyond the max shift
1910// printf(" Refine Failed %1.2f %1.2f %1.1f\n",(float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1),(float)gsl_vector_get(s->x, 2));
1911 result = this_img->process("xform",Dict("transform",t));
1912 result->set_attr("xform.align2d",t);
1913 }
1914
1915 delete t;
1916 t = 0;
1917
1918 gsl_vector_free(x);
1919 gsl_vector_free(ss);
1920 gsl_multimin_fminimizer_free(s);
1921
1922 if (c != 0) delete c;
1923 return result;
1924}
1925
1927 const string & cmp_name, const Dict& cmp_params) const
1928{
1929
1930 if (!to) {
1931 return 0;
1932 }
1933
1934 EMData *result;
1935 int mode = params.set_default("mode", 0);
1936 float saz = 0.0;
1937 float sdx = 0.0;
1938 float sdy = 0.0;
1939 float sscale = 1.0;
1940 bool mirror = false;
1941 Transform* t;
1942 if (params.has_key("xform.align2d") ) {
1943 t = params["xform.align2d"];
1944 Dict params = t->get_params("2d");
1945 saz = params["alpha"];
1946 sdx = params["tx"];
1947 sdy = params["ty"];
1948 mirror = params["mirror"];
1949 sscale = params["scale"];
1950 } else {
1951 t = new Transform(); // is the identity
1952 }
1953
1954 // We do this to prevent the GSL routine from crashing on an invalid alignment
1955 if ((float)(this_img->get_attr("sigma"))==0.0 || (float)(to->get_attr("sigma"))==0.0) {
1956 result = this_img->process("xform",Dict("transform",t));
1957 result->set_attr("xform.align2d",t);
1958 delete t;
1959 return result;
1960 }
1961
1962 float step = params.set_default("step",0.1f);
1963 float stepscale = params.set_default("stepscale",0.0f);
1964
1965 int np = 3;
1966 if (stepscale!=0.0) np++;
1967 Dict gsl_params;
1968 gsl_params["this"] = this_img;
1969 gsl_params["with"] = to;
1970 gsl_params["snr"] = params["snr"];
1971 gsl_params["mirror"] = mirror;
1972 if (params.has_key("mask")) gsl_params["mask"]=params["mask"];
1973
1974 const gsl_multimin_fdfminimizer_type *T = gsl_multimin_fdfminimizer_vector_bfgs;
1975
1976 gsl_vector *x = gsl_vector_alloc(np);
1977 gsl_vector_set(x, 0, sdx);
1978 gsl_vector_set(x, 1, sdy);
1979 gsl_vector_set(x, 2, saz);
1980 if (stepscale!=0.0) gsl_vector_set(x,3,1.0);
1981
1982 Cmp *c = 0;
1983
1984 gsl_multimin_function_fdf minex_func;
1985 if (mode == 2) {
1986 minex_func.f = &refalifnfast;
1987 }
1988 else {
1989 c = Factory < Cmp >::get(cmp_name, cmp_params);
1990 gsl_params["cmp"] = (void *) c;
1991 minex_func.f = &refalifn;
1992 }
1993
1994 minex_func.df = &refalidf;
1995 minex_func.fdf = &refalifdf;
1996 minex_func.n = np;
1997 minex_func.params = (void *) &gsl_params;
1998
1999 gsl_multimin_fdfminimizer *s = gsl_multimin_fdfminimizer_alloc(T, np);
2000 gsl_multimin_fdfminimizer_set(s, &minex_func, x, step, 0.001f);
2001
2002 int rval = GSL_CONTINUE;
2003 int status = GSL_SUCCESS;
2004 int iter = 1;
2005
2006 float precision = params.set_default("precision",0.02f);
2007 int maxiter = params.set_default("maxiter",12);
2008 int verbose = params.set_default("verbose",0);
2009
2010// printf("Refine sx=%1.2f sy=%1.2f sa=%1.2f prec=%1.4f maxit=%d\n",stepx,stepy,stepaz,precision,maxiter);
2011// printf("%1.2f %1.2f %1.1f ->",(float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1),(float)gsl_vector_get(s->x, 2));
2012
2013 while (rval == GSL_CONTINUE && iter < maxiter) {
2014 iter++;
2015 status = gsl_multimin_fdfminimizer_iterate(s);
2016 if (status) {
2017 break;
2018 }
2019 rval = gsl_multimin_test_gradient (s->gradient, precision);
2020// if (verbose>2) printf("GSL %d. %1.3f %1.3f %1.3f %1.3f\n",iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1),gsl_vector_get(s->x,2),s->gradient[0]);
2021 }
2022
2023 int maxshift = params.set_default("maxshift",-1);
2024
2025 if (maxshift <= 0) {
2026 maxshift = this_img->get_xsize() / 4;
2027 }
2028 float fmaxshift = static_cast<float>(maxshift);
2029 if ( fmaxshift >= fabs((float)gsl_vector_get(s->x, 0)) && fmaxshift >= fabs((float)gsl_vector_get(s->x, 1)) && (stepscale==0 || (((float)gsl_vector_get(s->x, 3))<1.3 && ((float)gsl_vector_get(s->x, 3))<0.7)) )
2030 {
2031 if (verbose>0) printf(" Refine good (%d) %1.2f %1.2f %1.1f\n",iter,(float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1),(float)gsl_vector_get(s->x, 2));
2032 Transform tsoln(Dict("type","2d","alpha",(float)gsl_vector_get(s->x, 2)));
2033 tsoln.set_mirror(mirror);
2034 tsoln.set_trans((float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1));
2035 if (stepscale!=0.0) tsoln.set_scale((float)gsl_vector_get(s->x, 3));
2036 result = this_img->process("xform",Dict("transform",&tsoln));
2037 result->set_attr("xform.align2d",&tsoln);
2038 } else { // The refine aligner failed - this shift went beyond the max shift
2039 if (verbose>1) printf(" Refine Failed %1.2f %1.2f %1.1f\n",(float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1),(float)gsl_vector_get(s->x, 2));
2040 result = this_img->process("xform",Dict("transform",t));
2041 result->set_attr("xform.align2d",t);
2042 }
2043
2044 delete t;
2045 t = 0;
2046
2047 gsl_vector_free(x);
2048 gsl_multimin_fdfminimizer_free(s);
2049
2050 if (c != 0) delete c;
2051 return result;
2052}
2053
2054static Transform refalin3d_perturbquat(const Transform*const t, const float& spincoeff, const float& n0, const float& n1, const float& n2, const float& x, const float& y, const float& z)
2055{
2056 Vec3f normal(n0,n1,n2);
2057 normal.normalize();
2058
2059 float omega = spincoeff*sqrt(n0*n0 + n1*n1 + n2*n2); // Here we compute the spin by the rotation axis vector length
2060 Dict d;
2061 d["type"] = "spin";
2062 d["omega"] = omega;
2063 d["n1"] = normal[0];
2064 d["n2"] = normal[1];
2065 d["n3"] = normal[2];
2066 //cout << omega << " " << normal[0] << " " << normal[1] << " " << normal[2] << " " << n0 << " " << n1 << " " << n2 << endl;
2067
2068 Transform q(d);
2069 q.set_trans((float)x,(float)y,(float)z);
2070
2071 q = q*(*t); //compose transforms
2072
2073 return q;
2074}
2075
2076static double symquat(const gsl_vector * v, void *params)
2077{
2078 Dict *dict = (Dict *) params;
2079
2080 double n0 = gsl_vector_get(v, 0);
2081 double n1 = gsl_vector_get(v, 1);
2082 double n2 = gsl_vector_get(v, 2);
2083 double x = gsl_vector_get(v, 3);
2084 double y = gsl_vector_get(v, 4);
2085 double z = gsl_vector_get(v, 5);
2086
2087 EMData* volume = (*dict)["volume"];
2088 float spincoeff = (*dict)["spincoeff"];
2089 Transform* t = (*dict)["transform"];
2090
2091 Transform soln = refalin3d_perturbquat(t,spincoeff,(float)n0,(float)n1,(float)n2,(float)x,(float)y,(float)z);
2092
2093 EMData *tmp = volume->process("xform",Dict("transform",&soln));
2094 EMData *symtmp = tmp->process("xform.applysym",Dict("sym",(*dict)["sym"]));
2095 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
2096 double result = c->cmp(symtmp,tmp);
2097 delete tmp;
2098 delete symtmp;
2099
2100 //cout << result << endl;
2101 return result;
2102}
2103
2104static double refalifn3dquat(const gsl_vector * v, void *params)
2105{
2106 Dict *dict = (Dict *) params;
2107
2108 double n0 = gsl_vector_get(v, 0);
2109 double n1 = gsl_vector_get(v, 1);
2110 double n2 = gsl_vector_get(v, 2);
2111 double x = gsl_vector_get(v, 3);
2112 double y = gsl_vector_get(v, 4);
2113 double z = gsl_vector_get(v, 5);
2114
2115 EMData *this_img = (*dict)["this"];
2116 EMData *with = (*dict)["with"];
2117
2118 Transform* t = (*dict)["transform"];
2119 float spincoeff = (*dict)["spincoeff"];
2120
2121 Transform soln = refalin3d_perturbquat(t,spincoeff,(float)n0,(float)n1,(float)n2,(float)x,(float)y,(float)z);
2122
2123 EMData *tmp = this_img->process("xform",Dict("transform",&soln));
2124 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
2125 double result = c->cmp(tmp,with);
2126 if ( tmp != 0 ) delete tmp;
2127
2128 //cout << result << endl;
2129 return result;
2130}
2131
2132EMData* SymAlignProcessorQuat::align(EMData * volume, EMData *to, const string & cmp_name, const Dict& cmp_params) const
2133{
2134 //Get pretransform
2135 Transform* t;
2136 if (params.has_key("xform.align3d") ) {
2137 t = params["xform.align3d"];
2138 }else {
2139 t = new Transform(); // is the identity
2140 }
2141
2142 float sdi = 0.0;
2143 float sdj = 0.0;
2144 float sdk = 0.0;
2145 float sdx = 0.0;
2146 float sdy = 0.0;
2147 float sdz = 0.0;
2148
2149 float spincoeff = params.set_default("spin_coeff",10.0f); // spin coefficient, controls speed of convergence (sort of)
2150
2151 int np = 6; // the number of dimensions
2152 Dict gsl_params;
2153 gsl_params["volume"] = volume;
2154 gsl_params["transform"] = t;
2155 gsl_params["sym"] = params.set_default("sym","c1");
2156 gsl_params["spincoeff"] = spincoeff;
2157
2158 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
2159 gsl_vector *ss = gsl_vector_alloc(np);
2160
2161 float stepi = params.set_default("stepn0",1.0f); // doesn't really matter b/c the vecor part will be normalized anyway
2162 float stepj = params.set_default("stepn1",1.0f); // doesn't really matter b/c the vecor part will be normalized anyway
2163 float stepk = params.set_default("stepn2",1.0f); // doesn't really matter b/c the vecor part will be normalized anyway
2164 float stepx = params.set_default("stepx",1.0f);
2165 float stepy = params.set_default("stepy",1.0f);
2166 float stepz = params.set_default("stepz",1.0f);
2167
2168 gsl_vector_set(ss, 0, stepi);
2169 gsl_vector_set(ss, 1, stepj);
2170 gsl_vector_set(ss, 2, stepk);
2171 gsl_vector_set(ss, 3, stepx);
2172 gsl_vector_set(ss, 4, stepy);
2173 gsl_vector_set(ss, 5, stepz);
2174
2175 gsl_vector *x = gsl_vector_alloc(np);
2176 gsl_vector_set(x, 0, sdi);
2177 gsl_vector_set(x, 1, sdj);
2178 gsl_vector_set(x, 2, sdk);
2179 gsl_vector_set(x, 3, sdx);
2180 gsl_vector_set(x, 4, sdy);
2181 gsl_vector_set(x, 5, sdz);
2182
2183 gsl_multimin_function minex_func;
2184 Cmp *c = Factory < Cmp >::get(cmp_name, cmp_params);
2185 gsl_params["cmp"] = (void *) c;
2186 minex_func.f = &symquat;
2187 minex_func.n = np;
2188 minex_func.params = (void *) &gsl_params;
2189 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
2190 gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
2191
2192 int rval = GSL_CONTINUE;
2193 int status = GSL_SUCCESS;
2194 int iter = 1;
2195
2196 float precision = params.set_default("precision",0.01f);
2197 int maxiter = params.set_default("maxiter",100);
2198 while (rval == GSL_CONTINUE && iter < maxiter) {
2199 iter++;
2200 status = gsl_multimin_fminimizer_iterate(s);
2201 if (status) {
2202 break;
2203 }
2204 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
2205 }
2206
2207 int maxshift = params.set_default("maxshift",-1);
2208
2209 if (maxshift <= 0) {
2210 maxshift = volume->get_xsize() / 4;
2211 }
2212 float fmaxshift = static_cast<float>(maxshift);
2213
2214 EMData *result;
2215 if ( fmaxshift >= (float)gsl_vector_get(s->x, 0) && fmaxshift >= (float)gsl_vector_get(s->x, 1) && fmaxshift >= (float)gsl_vector_get(s->x, 2))
2216 {
2217 float n0 = (float)gsl_vector_get(s->x, 0);
2218 float n1 = (float)gsl_vector_get(s->x, 1);
2219 float n2 = (float)gsl_vector_get(s->x, 2);
2220 float x = (float)gsl_vector_get(s->x, 3);
2221 float y = (float)gsl_vector_get(s->x, 4);
2222 float z = (float)gsl_vector_get(s->x, 5);
2223
2224 Transform tsoln = refalin3d_perturbquat(t,spincoeff,n0,n1,n2,x,y,z);
2225
2226 result = volume->process("xform",Dict("transform",&tsoln));
2227 result->set_attr("xform.align3d",&tsoln);
2228 EMData *tmpsym = result->process("xform.applysym",Dict("sym",gsl_params["sym"]));
2229 result->set_attr("score", result->cmp(cmp_name,tmpsym,cmp_params));
2230 delete tmpsym;
2231 } else { // The refine aligner failed - this shift went beyond the max shift
2232 result = volume->process("xform",Dict("transform",t));
2233 result->set_attr("xform.align3d",t);
2234 result->set_attr("score",0.0);
2235 }
2236
2237 gsl_vector_free(x);
2238 gsl_vector_free(ss);
2239 gsl_multimin_fminimizer_free(s);
2240
2241 if (c != 0) delete c;
2242 delete t;
2243
2244 return result;
2245}
2246
2248 const string & cmp_name, const Dict& cmp_params) const
2249{
2250
2251 if (!to || !this_img) throw NullPointerException("Input image is null"); // not sure if this is necessary, it was there before I started
2252
2253 if (to->get_ndim() != 3 || this_img->get_ndim() != 3) throw ImageDimensionException("The Refine3D aligner only works for 3D images");
2254
2255#ifdef EMAN2_USING_CUDA
2256 if(EMData::usecuda == 1) {
2257 if(!this_img->getcudarwdata()) this_img->copy_to_cuda();
2258 if(!to->getcudarwdata()) to->copy_to_cuda();
2259 }
2260#endif
2261
2262 float sdi = 0.0;
2263 float sdj = 0.0;
2264 float sdk = 0.0;
2265 float sdx = 0.0;
2266 float sdy = 0.0;
2267 float sdz = 0.0;
2268 bool mirror = false;
2269
2270 Transform* t;
2271 if (params.has_key("xform.align3d") ) {
2272 // Unlike the 2d refine aligner, this class doesn't require the starting transform's
2273 // parameters to form the starting guess. Instead the Transform itself
2274 // is perturbed carefully (using quaternion rotation) to overcome problems that arise
2275 // when you use orthogonally-based Euler angles
2276 t = params["xform.align3d"];
2277 }else {
2278 t = new Transform(); // is the identity
2279 }
2280
2281 float spincoeff = params.set_default("spin_coeff",10.0f); // spin coefficient, controls speed of convergence (sort of)
2282
2283 int np = 6; // the number of dimensions
2284 Dict gsl_params;
2285 gsl_params["this"] = this_img;
2286 gsl_params["with"] = to;
2287 gsl_params["snr"] = params["snr"];
2288 gsl_params["mirror"] = mirror;
2289 gsl_params["transform"] = t;
2290 gsl_params["spincoeff"] = spincoeff;
2291 Dict altered_cmp_params(cmp_params);
2292
2293 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
2294 gsl_vector *ss = gsl_vector_alloc(np);
2295
2296 float stepi = params.set_default("stepn0",1.0f); // doesn't really matter b/c the vecor part will be normalized anyway
2297 float stepj = params.set_default("stepn1",1.0f); // doesn't really matter b/c the vecor part will be normalized anyway
2298 float stepk = params.set_default("stepn2",1.0f); // doesn't really matter b/c the vecor part will be normalized anyway
2299 float stepx = params.set_default("stepx",1.0f);
2300 float stepy = params.set_default("stepy",1.0f);
2301 float stepz = params.set_default("stepz",1.0f);
2302
2303 //gsl_vector_set(ss, 0, stepw);
2304 gsl_vector_set(ss, 0, stepi);
2305 gsl_vector_set(ss, 1, stepj);
2306 gsl_vector_set(ss, 2, stepk);
2307 gsl_vector_set(ss, 3, stepx);
2308 gsl_vector_set(ss, 4, stepy);
2309 gsl_vector_set(ss, 5, stepz);
2310
2311 gsl_vector *x = gsl_vector_alloc(np);
2312 gsl_vector_set(x, 0, sdi);
2313 gsl_vector_set(x, 1, sdj);
2314 gsl_vector_set(x, 2, sdk);
2315 gsl_vector_set(x, 3, sdx);
2316 gsl_vector_set(x, 4, sdy);
2317 gsl_vector_set(x, 5, sdz);
2318
2319 gsl_multimin_function minex_func;
2320 Cmp *c = Factory < Cmp >::get(cmp_name, altered_cmp_params);
2321
2322 gsl_params["cmp"] = (void *) c;
2323 minex_func.f = &refalifn3dquat;
2324
2325 minex_func.n = np;
2326 minex_func.params = (void *) &gsl_params;
2327
2328 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
2329 gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
2330
2331 int rval = GSL_CONTINUE;
2332 int status = GSL_SUCCESS;
2333 int iter = 1;
2334
2335 float precision = params.set_default("precision",0.01f);
2336 int maxiter = params.set_default("maxiter",100);
2337 while (rval == GSL_CONTINUE && iter < maxiter) {
2338 iter++;
2339 status = gsl_multimin_fminimizer_iterate(s);
2340 if (status) {
2341 break;
2342 }
2343 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
2344 }
2345
2346 int maxshift = params.set_default("maxshift",-1);
2347
2348 if (maxshift <= 0) {
2349 maxshift = this_img->get_xsize() / 4;
2350 }
2351 float fmaxshift = static_cast<float>(maxshift);
2352
2353 EMData *result;
2354 if ( fmaxshift >= (float)gsl_vector_get(s->x, 0) && fmaxshift >= (float)gsl_vector_get(s->x, 1) && fmaxshift >= (float)gsl_vector_get(s->x, 2))
2355 {
2356 float n0 = (float)gsl_vector_get(s->x, 0);
2357 float n1 = (float)gsl_vector_get(s->x, 1);
2358 float n2 = (float)gsl_vector_get(s->x, 2);
2359 float x = (float)gsl_vector_get(s->x, 3);
2360 float y = (float)gsl_vector_get(s->x, 4);
2361 float z = (float)gsl_vector_get(s->x, 5);
2362
2363 Transform tsoln = refalin3d_perturbquat(t,spincoeff,n0,n1,n2,x,y,z);
2364
2365 result = this_img->process("xform",Dict("transform",&tsoln));
2366 result->set_attr("xform.align3d",&tsoln);
2367 result->set_attr("score", result->cmp(cmp_name,to,cmp_params));
2368
2369 //coda goes here
2370 } else { // The refine aligner failed - this shift went beyond the max shift
2371 result = this_img->process("xform",Dict("transform",t));
2372 result->set_attr("xform.align3d",t);
2373 result->set_attr("score",0.0);
2374 }
2375
2376 //EMData *result = this_img->process("xform",Dict("transform",t));
2377 delete t;
2378 gsl_vector_free(x);
2379 gsl_vector_free(ss);
2380 gsl_multimin_fminimizer_free(s);
2381
2382 if (c != 0) delete c;
2383
2384 return result;
2385}
2386
2388 const string & cmp_name, const Dict& cmp_params) const
2389{
2390 if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
2391 throw ImageDimensionException("This aligner only works for 3D images");
2392 }
2393
2394 Transform* t;
2395 if (params.has_key("xform.align3d") ) {
2396 // Unlike the 2d refine aligner, this class doesn't require the starting transform's
2397 // parameters to form the starting guess. Instead the Transform itself
2398 // is perturbed carefully (using quaternion rotation) to overcome problems that arise
2399 // when you use orthogonally-based Euler angles
2400 t = params["xform.align3d"];
2401 }else {
2402 t = new Transform(); // is the identity
2403 }
2404
2405 int searchx = 0;
2406 int searchy = 0;
2407 int searchz = 0;
2408 bool dotrans = params.set_default("dotrans",1);
2409 if (params.has_key("search")) {
2410 vector<string> check;
2411 check.push_back("searchx");
2412 check.push_back("searchy");
2413 check.push_back("searchz");
2414 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
2415 if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
2416 }
2417 int search = params["search"];
2418 searchx = search;
2419 searchy = search;
2420 searchz = search;
2421 } else {
2422 searchx = params.set_default("searchx",3);
2423 searchy = params.set_default("searchy",3);
2424 searchz = params.set_default("searchz",3);
2425 }
2426
2427 float delta = params.set_default("delta",1.0f);
2428 float range = params.set_default("range",10.0f);
2429 bool verbose = params.set_default("verbose",false);
2430
2431 bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
2432 EMData * tofft = 0;
2433 if(dotrans || tomography){
2434 tofft = to->do_fft();
2435 }
2436
2437#ifdef EMAN2_USING_CUDA
2438 if(EMData::usecuda == 1) {
2439 if(!this_img->getcudarodata()) this_img->copy_to_cudaro(); // This is safer
2440 if(!to->getcudarwdata()) to->copy_to_cuda();
2441 if(to->getcudarwdata()){if(tofft) tofft->copy_to_cuda();}
2442 }
2443#endif
2444
2445 Dict d;
2446 d["type"] = "eman"; // d is used in the loop below
2447 Dict best;
2448// best["score"] = numeric_limits<float>::infinity();
2449 best["score"] = 1.0e37;
2450 bool use_cpu = true;
2451 Transform tran = Transform();
2452 Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
2453
2454 for ( float alt = 0; alt < range; alt += delta) {
2455 // now compute a sane az step size
2456 float saz = 360;
2457 if(alt != 0) saz = delta/sin(alt*M_PI/180.0f); // This gives consistent az step sizes(arc lengths)
2458 for ( float az = 0; az < 360; az += saz ){
2459 if (verbose) {
2460 cout << "Trying angle alt " << alt << " az " << az << endl;
2461 }
2462 // account for any changes in az
2463 for( float phi = -range-az; phi < range-az; phi += delta ) {
2464 d["alt"] = alt;
2465 d["phi"] = phi;
2466 d["az"] = az;
2467 Transform tr(d);
2468 tr = tr*(*t); // compose transforms, this moves to the pole (aprox)
2469
2470 //EMData* transformed = this_img->process("xform",Dict("transform",&tr));
2471 EMData* transformed;
2472 transformed = this_img->process("xform",Dict("transform",&tr));
2473
2474 //need to do things a bit diffrent if we want to compare two tomos
2475 float score = 0.0f;
2476 if(dotrans || tomography){
2477 EMData* ccf = transformed->calc_ccf(tofft);
2478#ifdef EMAN2_USING_CUDA
2479 if(EMData::usecuda == 1){
2480 use_cpu = false;
2481 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
2482 tran.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
2483 //CudaPeakInfoFloat* data = calc_max_location_wrap_intp_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
2484 //tran.set_trans(-data->xintp, -data->yintp, -data->zintp);
2485 tr = tran*tr; // to reflect the fact that we have done a rotation first and THEN a transformation
2486 if (tomography) {
2487 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
2488 score = -(data->peak - stats.x)/sqrt(stats.y); // Normalize, this is better than calling the norm processor since we only need to normalize one point
2489 } else {
2490 score = -data->peak;
2491 }
2492 delete data;
2493 }
2494#endif
2495 if(use_cpu){
2496 if(tomography) ccf->process_inplace("normalize");
2497 //vector<float> fpoint = ccf->calc_max_location_wrap_intp(searchx,searchy,searchz);
2498 //tran.set_trans(-fpoint[0], -fpoint[1], -fpoint[2]);
2499 //score = -fpoint[3];
2500 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
2501 tran.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
2502 score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
2503 tr = tran*tr;// to reflect the fact that we have done a rotation first and THEN a transformation
2504
2505 }
2506 delete ccf; ccf =0;
2507 delete transformed; transformed = 0;// this is to stop a mem leak
2508 }
2509
2510 if(!tomography){
2511 if(!transformed) transformed = this_img->process("xform",Dict("transform",&tr)); // we are returning a map
2512 score = c->cmp(to,transformed);
2513 delete transformed; transformed = 0;// this is to stop a mem leak
2514 }
2515
2516 if(score < float(best["score"])) {
2517// printf("%f\n",score);
2518 best["score"] = score;
2519 best["xform.align3d"] = &tr; // I wonder if this will cause a mem leak?
2520 }
2521 }
2522 }
2523 }
2524
2525 if(tofft) {delete tofft; tofft = 0;}
2526 if (c != 0) delete c;
2527
2528 //make aligned map;
2529 EMData* best_match = this_img->process("xform",Dict("transform", best["xform.align3d"])); // we are returning a map
2530 best_match->set_attr("xform.align3d", best["xform.align3d"]);
2531 best_match->set_attr("score", float(best["score"]));
2532
2533 return best_match;
2534
2535}
2536
2537EMData* RT3DGridAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
2538{
2539
2540 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
2541
2542 Dict t;
2543 Transform* tr = (Transform*) alis[0]["xform.align3d"];
2544 t["transform"] = tr;
2545 EMData* soln = this_img->process("xform",t);
2546 soln->set_attr("xform.align3d",tr);
2547 delete tr; tr = 0;
2548
2549 return soln;
2550
2551}
2552
2553vector<Dict> RT3DGridAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const {
2554
2555 if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
2556 throw ImageDimensionException("This aligner only works for 3D images");
2557 }
2558
2559 int searchx = 0;
2560 int searchy = 0;
2561 int searchz = 0;
2562
2563 bool dotrans = params.set_default("dotrans",1);
2564 if (params.has_key("search")) {
2565 vector<string> check;
2566 check.push_back("searchx");
2567 check.push_back("searchy");
2568 check.push_back("searchz");
2569 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
2570 if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
2571 }
2572 int search = params["search"];
2573 searchx = search;
2574 searchy = search;
2575 searchz = search;
2576 } else {
2577 searchx = params.set_default("searchx",3);
2578 searchy = params.set_default("searchy",3);
2579 searchz = params.set_default("searchz",3);
2580 }
2581
2582 Transform* initxform;
2583 if (params.has_key("initxform") ) {
2584 // Unlike the 2d refine aligner, this class doesn't require the starting transform's
2585 // parameters to form the starting guess. Instead the Transform itself
2586 // is perturbed carefully (using quaternion rotation) to overcome problems that arise
2587 // when you use orthogonally-based Euler angles
2588 initxform = params["initxform"];
2589 }else {
2590 initxform = new Transform(); // is the identity
2591 }
2592
2593 float lalt = params.set_default("alt0",0.0f);
2594 float laz = params.set_default("az0",0.0f);
2595 float lphi = params.set_default("phi0",0.0f);
2596 float ualt = params.set_default("alt1",180.0f); // I am using 179.9 rather than 180 to avoid resampling
2597 float uphi = params.set_default("phi1",360.0f); // I am using 359.9 rather than 180 to avoid resampling 0 = 360 (for perodic functions)
2598 float uaz = params.set_default("az1",360.0f); // I am using 359.9 rather than 180 to avoid resampling 0 = 360 (for perodic functions)
2599 float dalt = params.set_default("dalt",10.f);
2600 float daz = params.set_default("daz",10.f);
2601 float dphi = params.set_default("dphi",10.f);
2602 bool verbose = params.set_default("verbose",false);
2603
2604 //in case we arre aligning tomos
2605 Dict altered_cmp_params(cmp_params);
2606 if (cmp_name == "ccc.tomo") {
2607 altered_cmp_params.set_default("searchx", searchx);
2608 altered_cmp_params.set_default("searchy", searchy);
2609 altered_cmp_params.set_default("searchz", searchz);
2610 altered_cmp_params.set_default("norm", true);
2611 }
2612
2613 vector<Dict> solns;
2614 if (nsoln == 0) return solns; // What was the user thinking?
2615 for (unsigned int i = 0; i < nsoln; ++i ) {
2616 Dict d;
2617 d["score"] = 1.e24;
2618 Transform t; // identity by default
2619 d["xform.align3d"] = &t; // deep copy is going on here
2620 solns.push_back(d);
2621 }
2622
2623 bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
2624 EMData * tofft = 0;
2625 if(dotrans || tomography){
2626 tofft = to->do_fft();
2627 }
2628
2629#ifdef EMAN2_USING_CUDA
2630 if(EMData::usecuda == 1) {
2631 if(!this_img->getcudarodata()) this_img->copy_to_cudaro(); // safer call
2632 if(!to->getcudarwdata()) to->copy_to_cuda();
2633 if(to->getcudarwdata()){if(tofft) tofft->copy_to_cuda();}
2634 }
2635#endif
2636
2637 Dict d;
2638 d["type"] = "eman"; // d is used in the loop below
2639 Transform trans = Transform();
2640 Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
2641 bool use_cpu = true;
2642 for ( float alt = lalt; alt <= ualt; alt += dalt) {
2643 // An optimization for the range of az is made at the top of the sphere
2644 // If you think about it, this is just a coarse way of making this approach slightly more efficient
2645 for ( float az = laz; az < uaz; az += daz ){
2646 if (verbose) {
2647 cout << "Trying angle alt " << alt << " az " << az << endl;
2648 }
2649 for( float phi = lphi; phi < uphi; phi += dphi ) {
2650 d["alt"] = alt;
2651 d["phi"] = phi;
2652 d["az"] = az;
2653 Transform t(d);
2654 t = t*(*initxform);
2655 EMData* transformed = this_img->process("xform",Dict("transform",&t));
2656
2657 //need to do things a bit diffrent if we want to compare two tomos
2658 float best_score = 0.0f;
2659 if(dotrans || tomography){
2660 EMData* ccf = transformed->calc_ccf(tofft);
2661#ifdef EMAN2_USING_CUDA
2662 if(EMData::usecuda == 1){
2663 use_cpu = false;;
2664 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
2665 trans.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
2666 t = trans*t; //composite transfrom to reflect the fact that we have done a rotation first and THEN a transformation
2667 if (tomography) {
2668 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
2669 best_score = -(data->peak - stats.x)/sqrt(stats.y); // Normalize, this is better than calling the norm processor since we only need to normalize one point
2670 } else {
2671 best_score = -data->peak;
2672 }
2673 delete data;
2674 }
2675#endif
2676 if(use_cpu){
2677 if(tomography) ccf->process_inplace("normalize");
2678 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
2679 trans.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
2680 t = trans*t; //composite transfrom to reflect the fact that we have done a rotation first and THEN a transformation
2681 best_score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
2682 }
2683 delete ccf; ccf =0;
2684 delete transformed; transformed = 0;
2685 }
2686
2687 if(!tomography){
2688 if(!transformed) transformed = this_img->process("xform",Dict("transform",&t));
2689 best_score = c->cmp(to,transformed);
2690 delete transformed; transformed = 0;
2691 }
2692
2693 unsigned int j = 0;
2694 for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
2695 if ( (float)(*it)["score"] > best_score ) { // Note greater than - EMAN2 preferes minimums as a matter of policy
2696 vector<Dict>::reverse_iterator rit = solns.rbegin();
2697 copy(rit+1,solns.rend()-j,rit);
2698 Dict& d = (*it);
2699 d["score"] = best_score;
2700 d["xform.align3d"] = &t;
2701 break;
2702 }
2703 }
2704 }
2705 }
2706 }
2707
2708 if(tofft) {delete tofft; tofft = 0;}
2709 if (c != 0) delete c;
2710
2711 return solns;
2712
2713}
2714
2715EMData* RT2DTreeAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
2716{
2717
2718 vector<Dict> alis = xform_align_nbest(this_img,to,2,cmp_name,cmp_params);
2719
2720 Dict t;
2721 Transform* tr = (Transform*) alis[0]["xform.align2d"];
2722
2723 t["transform"] = tr;
2724 EMData* soln = this_img->process("xform",t);
2725 soln->set_attr("xform.align2d",tr);
2726
2727 return soln;
2728
2729}
2730
2731vector<Dict> RT2DTreeAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nrsoln, const string & cmp_name, const Dict& cmp_params) const {
2732 if (nrsoln == 0) throw InvalidParameterException("ERROR (RT2DTreeAligner): nsoln must be >0"); // What was the user thinking?
2733
2734 int nsoln = nrsoln*2;
2735 if (nrsoln<8) nsoln=8; // we start with at least n solutions, but then gradually decrease with increasing scale
2736
2737 // !!!!!! IMPORTANT NOTE - we are inverting the order of 'this' and 'to' here to match convention in other aligners, to compensate
2738 // the Transform is inverted before being returned
2739 EMData *base_this;
2740 EMData *base_to;
2741 if (this_img->is_complex()) base_this=this_img->copy();
2742 else {
2743 base_this=this_img->do_fft();
2744 base_this->process_inplace("xform.phaseorigin.tocorner"); // This was originally .tocorner, taken from RT3DTree, but I think that's probably wrong...
2745 }
2746
2747 if (to->is_complex()) base_to=to->copy();
2748 else {
2749 base_to=to->do_fft();
2750 base_to->process_inplace("xform.phaseorigin.tocorner");
2751 }
2752
2753 int verbose = params.set_default("verbose",0);
2754 int maxshift = params.set_default("maxshift",-1);
2755 int flip = params.set_default("flip",1);
2756 float maxres = params.set_default("maxres",-1.0f);
2757 if (maxres<0.1) maxres=0.1;
2758
2759 if (base_this->get_xsize()!=base_this->get_ysize()+2 || base_to->get_xsize()!=base_to->get_ysize()+2 ) throw InvalidCallException("ERROR (RT2DTreeAligner): requires cubic images with even numbered box sizes");
2760
2761 base_this->process_inplace("xform.fourierorigin.tocenter"); // easier to chop out Fourier subvolumes
2762 base_to->process_inplace("xform.fourierorigin.tocenter");
2763
2764 float apix=(float)this_img->get_attr("apix_x");
2765 int ny=this_img->get_ysize();
2766 if (maxshift<0) maxshift=ny*3/8;
2767 float maxs = ny*apix/maxres;
2768
2769// int downsample=floor(ny/20); // Minimum shrunken box size is 20^3
2770
2771 vector<float> s_score(nsoln,0.0f);
2772 vector<Transform> s_xform(nsoln);
2773// vector<float> s_step(nsoln,7.5);
2774 if (verbose>0) printf("%d solutions\n",nsoln);
2775
2776 // We start with 32^3, 64^3 ...
2777 for (int sexp=5; sexp<10; sexp++) {
2778 int ss=pow(2.0,sexp);
2779 float rescale=2.0;
2780 if (ss>ny) {
2781 rescale=(float)ny/pow(2.0,sexp-1);
2782 ss=ny;
2783 }
2784 if (verbose>0) printf("\nSize %d\n",ss);
2785
2786 //ss=good_size(ny/ds);
2787 // Clearly these regions will be messed up on x=nyquist edge, but we are zeroing that pixel during rotation anyway
2788 EMData *small_this=base_this->get_clip(Region(0,(ny-ss)/2,ss+2,ss));
2789 EMData *small_to= base_to-> get_clip(Region(0,(ny-ss)/2,ss+2,ss));
2790 small_this->process_inplace("xform.fourierorigin.tocorner"); // after clipping back to canonical form
2791 float cut2=0.33*ny;
2792 if (maxs<cut2) cut2=maxs;
2793 small_this->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",3));
2794 small_this->process_inplace("filter.lowpass.gauss",Dict("cutoff_pixels",cut2));
2795 small_to->process_inplace("xform.fourierorigin.tocorner");
2796 small_to->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",3));
2797 small_to->process_inplace("filter.lowpass.gauss",Dict("cutoff_pixels",cut2));
2798
2799 // This is a solid estimate for very complete searching
2800 float astep = 360.0/int(M_PI/atan(1.25/ss)); // Decent angular step in degrees
2801
2802 // This insures we make at least one real effort at each level
2803 for (int i=0; i<nsoln; i++) {
2804 s_score[i]=-1.0e24; // reset the scores since the different scales will not match
2805// if (fabs(s_step[i])<astep/4.0) s_step[i]*=2.0;
2806 }
2807
2808
2809 // This is for the first loop, we do a full search in a (normally) heavily downsampled space
2810 if (sexp==5) {
2811 if (verbose>1) printf("stage 1 - ang step %1.2f, filter %1.3f\n",astep,cut2/(apix*ny));
2812 vector<Transform> transforms;
2813 for (float a=astep/2.0; a<359.9; a+=astep) {
2814 transforms.push_back(Transform(Dict("type","2d","alpha",a,"mirror",0)));
2815 if (flip) transforms.push_back(Transform(Dict("type","2d","alpha",a,"mirror",1)));
2816 }
2817 if (verbose>0) printf("%lu orientations to test\n",transforms.size());
2818
2819 // We iterate over all orientations
2820 for (unsigned int it=0; it<transforms.size(); it++) {
2821 Transform t = transforms[it];
2822
2823 // somewhat strangely, rotations are actually much more expensive than FFTs, so we use a CCF for translation
2824 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
2825 EMData *ccf=small_to->calc_ccf(stt);
2826
2827 int lmaxs=maxshift/ss;
2828 IntPoint ml(0,0,0);
2829 float sim=ccf->get_attr("minimum");
2830 int mir=t.get_mirror()?-1:1;
2831 for (int y=-lmaxs; y<=lmaxs; y++) {
2832 for (int x=-lmaxs; x<=lmaxs; x++) {
2833 float v=ccf->get_value_at_wrap(x,y);
2834 if (v>sim) {
2835 sim=v;
2836 ml[0]=x;
2837 ml[1]=y;
2838 ml[2]=0;
2839 }
2840 }
2841 }
2842
2843// IntPoint ml=ccf->calc_max_location_wrap();
2844
2845// Dict aap=t.get_params("2d");
2846// aap["tx"]=(int)ml[0];
2847// aap["ty"]=(int)ml[1];
2848// aap["tz"]=(int)ml[2];
2849
2850 t.set_params(Dict("tx",(int)ml[0],"ty",(int)ml[1],"tz",(int)ml[2]));
2851 delete stt;
2852 delete ccf;
2853// stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1)); // we have to do 1 slow transform here now that we have the translation
2854//
2855// sim=stt->cmp("ccc",small_to);
2856
2857 // could have used a priority queue, but would have required more infrastructure
2858 int worst=0;
2859 // First we find the worst solution in the list of possible best solutions, or the first
2860 // solution which is currently "empty"
2861//printf("a\n");
2862 for (int i=0; i<nsoln; i++) {
2863 if (s_score[i]==-1.0e24) { worst=i; break; }
2864 if (s_score[i]<s_score[worst]) worst=i;
2865 }
2866
2867 // If the current solution is better than the 'worst' of the previous solutions, then we
2868 // displace it. Note that there is no sorting performed here
2869 if (sim>s_score[worst]) {
2870 s_score[worst]=sim;
2871 s_xform[worst]=t;
2872 //printf("%f\t%f\t%d\n",s_score[worst],s_coverage[worst],worst);
2873 }
2874 }
2875
2876 if (verbose>2) {
2877 for (int i=0; i<nsoln; i++) {
2878 Dict aap=s_xform[i].get_params("2d");
2879 printf("%d) %d\t%1.1f\t%1.1f\t%1.2f\t%1.1f\n",i,(int)aap["mirror"],(float)aap["tx"],(float)aap["ty"],(float)aap["alpha"],s_score[i]);
2880 }
2881 }
2882
2883 }
2884 // Once we have our initial list of best locations, we just refine each possibility individually
2885 else {
2886 // We generate a search pattern around each existing solution
2887 if (verbose>1) printf("stage 2 (%1.2f, %1.2f, %d) filter: %1.3f\n",astep,rescale,nsoln,cut2/(apix*ny));
2888 for (int i=0; i<nsoln; i++) {
2889
2890 Dict aap=s_xform[i].get_params("2d");
2891 aap["tx"]=(float)aap["tx"]*rescale; // compensate for scaling steps, usually 2, but not in the last step
2892 aap["ty"]=(float)aap["ty"]*rescale;
2893 float alpha=aap["alpha"];
2894 for (float a=alpha-astep; a<=alpha+astep; a+=astep) {
2895 Transform t=Transform(Dict("type","2d","alpha",a,"mirror",aap["mirror"]));
2896 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
2897 EMData *ccf=small_to->calc_ccf(stt);
2898// ccf->process_inplace("xform.phaseorigin.tocenter");
2899 int bx,by;
2900 float ba;
2901 int mir=t.get_mirror()?-1:1;
2902 int tx=aap["tx"];
2903 int ty=aap["ty"];
2904 for (int y=-4+ty; y<=4+ty; y++) {
2905 for (int x=-4+tx; x<=4+tx; x++) {
2906 float v=ccf->get_value_at_wrap(x,y);
2907 if (v>s_score[i]) {
2908 t.set_trans(x,y);
2909 s_xform[i]=t;
2910 s_score[i]=v;
2911 }
2912 }
2913 }
2914 delete stt;
2915 delete ccf;
2916 }
2917
2918 if (verbose>2) {
2919 aap=s_xform[i].get_params("2d");
2920 printf("%d) %d\t%1.1f\t%1.1f\t%1.2f\t%1.1f\n",i,(int)aap["mirror"],(float)aap["tx"],(float)aap["ty"],(float)aap["alpha"],s_score[i]);
2921 }
2922 }
2923 }
2924 // lazy earlier in defining s_ vectors, so lazy here too and inefficiently sorting
2925 // We are sorting inside the outermost loop so we can decrease the number of solutions
2926 // before we get to the finest precision
2927 for (unsigned int i=0; i<nsoln-1; i++) {
2928 for (unsigned int j=i+1; j<nsoln; j++) {
2929 if (s_score[i]<s_score[j]) {
2930 float t=s_score[i]; s_score[i]=s_score[j]; s_score[j]=t;
2931 Transform tt=s_xform[i]; s_xform[i]=s_xform[j]; s_xform[j]=tt;
2932 }
2933 }
2934 }
2935
2936 // At each level of sampling we (potentially) decrease the number of answers we check in detail
2937 // assuming we are gradually homing in on the best solution
2938 nsoln/=2;
2939 if (nsoln<nrsoln) nsoln=nrsoln;
2940
2941
2942 delete small_this;
2943 delete small_to;
2944 if (ss==ny) break;
2945 }
2946
2947 delete base_this;
2948 delete base_to;
2949
2950
2951 // initialize results
2952 vector<Dict> solns;
2953 for (unsigned int i = 0; i < nrsoln; ++i ) {
2954 Dict d;
2955 d["score"] = -s_score[i];
2956// s_xform[i].invert(); // this is because we inverted the order of the input images above to match convention
2957 if (s_xform[i].get_mirror()) {
2958 Vec3f t=s_xform[i].get_trans();
2959 t[0]*=-1.0;
2960 s_xform[i].set_trans(t);
2961 }
2962 d["xform.align2d"] = &s_xform[i];
2963 solns.push_back(d);
2964 }
2965 if (verbose>1) {
2966 Dict aap=s_xform[0].get_params("2d");
2967 printf("Final: %d\t%1.1f\t%1.1f\t%1.2f\n",(int)aap["mirror"],(float)aap["tx"],(float)aap["ty"],(float)aap["alpha"]);
2968 }
2969
2970 return solns;
2971}
2972
2973EMData* RT2Dto3DTreeAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
2974{
2975
2976 int nsoln=8;
2977 if (params.has_key("initxform"))
2978 {
2980 const vector< Transform > xfs=params["initxform"];
2981 nsoln=xfs.size();
2982 }
2983
2984 vector<Dict> alis = xform_align_nbest(this_img,to,nsoln,cmp_name,cmp_params);
2985 Transform *tr = (Transform*) alis[0]["xform.align3d"];
2986
2987// t["transform"] = tr;
2988 // seems strange to transform the 2d image so just return the original
2989 EMData* soln = this_img->copy(); //process("xform",t);
2990 soln->set_attr("xform.align3d",tr);
2991 soln->set_attr("xform.projection",tr);
2992 soln->set_attr("score",alis[0]["score"]);
2993
2994 return soln;
2995}
2996
2997// NOTE - if symmetry is applied, it is critical that "to" be the volume which is already aligned to the symmetry axes (ie - the reference)
2998vector<Dict> RT2Dto3DTreeAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nrsoln, const string & cmp_name, const Dict& cmp_params) const {
2999 if (this_img->get_zsize()!=1 || to->get_zsize()==1) throw InvalidParameterException("ERROR (RT2Dto3DTreeAligner): first image must be 2D and second 3D");
3000
3001 if (nrsoln == 0) throw InvalidParameterException("ERROR (RT2Dto3DTreeAligner): nsoln must be >0"); // What was the user thinking?
3002
3003
3004 int verbose = params.set_default("verbose",0);
3005 int nsoln = nrsoln*2;
3006 if (nrsoln<16) nsoln=32; // we start with at least 32 solutions, but then gradually decrease with increasing scale
3007 if (params.has_key("initxform")){// refine alignment
3008 const vector< Transform > xfs=params["initxform"];
3009 nsoln=xfs.size();
3010 }
3011
3012 vector<float> s_score(nsoln,1.0e24);
3013 vector<float> s_step(nsoln*3,7.5f);
3014 vector<Transform> s_xform(nsoln);
3015
3016 int curiter=-1;
3017 int sexp_start=4;
3018 int ny=this_img->get_ysize();
3019 int maxshift00=(int)params.set_default("maxshift",ny/4);
3020 float maxres = params.set_default("maxres",-1.0f);
3021 float apix=(float)this_img->get_attr("apix_x");
3022 int maxny=ny;
3023 int maxrescut=ny/2;
3024 if (maxres>0){
3025 maxrescut=int(ny*apix/maxres);
3026 maxny=maxrescut*3/2*2;
3027 if(maxny>ny) maxny=ny;
3028 }
3029 if (verbose>0)
3030 printf("\n\n*******\nmax resolution %1.2f, box size %d\n", maxres, maxny);
3031
3032
3033 float maxang=params.set_default("maxang",-1.0);
3034 Transform initxf;
3035
3036 if (params.has_key("initxform")){
3037 const vector< Transform > xfs=params["initxform"];
3038 initxf.set_params(xfs[0].get_params("eman"));
3039 for (unsigned int i=0; i<nsoln; i++){
3040 s_xform[i].set_params(xfs[i].get_params("eman"));
3041 }
3042 sexp_start=9;
3043
3044
3045 curiter=0;
3046 float s=2.5; // random initial direction to avoid bias, increased initial step 12/30/20
3047 if ((maxang>0) && (s>maxang/2)){
3048 s=maxang/2;
3049 }
3050 std::random_device rd;
3051 std::uniform_int_distribution<int> dist(0, 1);
3052 for (int i=0; i<nsoln*3; i++) {
3053 s_step[i]=s*(dist(rd))?1.0:-1.0;
3054 }
3055 }
3056
3057 if (verbose>0) printf("%d solutions\n",nsoln);
3058
3059 // !!!!!! IMPORTANT NOTE - we are inverting the order of this and to here to match convention in other aligners, and the transform that project the 3D image to the 2D image is returned.
3060 EMData *base_this;
3061 EMData *base_to;
3062 if (this_img->is_complex()) base_to=this_img->copy();
3063 else {
3064 base_to=this_img->do_fft();
3065 base_to->process_inplace("xform.phaseorigin.tocorner");
3066 }
3067
3068 if (to->is_complex()) base_this=to->copy();
3069 else {
3070 base_this=to->do_fft();
3071 base_this->process_inplace("xform.phaseorigin.tocorner");
3072 }
3073
3074
3075
3076 if (base_this->get_xsize()!=base_this->get_ysize()+2 || base_this->get_ysize()!=base_this->get_zsize()
3077 || base_to->get_xsize()!=base_to->get_ysize()+2 /*|| base_to->get_ysize()!=base_to->get_zsize()*/) throw InvalidCallException("ERROR (RT3DTreeAligner): requires cubic images with even numbered box sizes");
3078
3079 base_this->process_inplace("xform.fourierorigin.tocenter"); // easier to chop out Fourier subvolumes
3080 base_to->process_inplace("xform.fourierorigin.tocenter");
3081
3082 params["boxsize"]=ny;
3083
3084
3085 string axname[] = {"az","alt","phi"};
3086
3087 // We start with 32^3, 64^3 ...
3088
3089 for (int sexp=sexp_start; sexp<10; sexp++) {
3090 curiter++;
3091 int ss=pow(2.0,sexp);
3092 if (ss==16) ss=24; // 16 may be too small, but 32 takes too long...
3093 if (ss==32) ss=48; // 16 may be too small, but 32 takes too long...
3094 if (ss>maxny) ss=maxny;
3095
3096 int maxshift=maxshift00*ss/ny;
3097 if (maxshift00<0) maxshift=-1;
3098 if (verbose>0) printf("\nSize %d, maxshift %d\n",ss, maxshift);
3099
3100 EMData *small_this=base_this->get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss)); // This one is 3D
3101 EMData *small_to= base_to-> get_clip(Region(0,(ny-ss)/2,0,ss+2,ss,1)); // to is 2D
3102 small_this->process_inplace("xform.fourierorigin.tocorner"); // after clipping back to canonical form
3103 small_this->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
3104
3105 small_to->process_inplace("xform.fourierorigin.tocorner");
3106 small_to->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
3107
3108 if (ss<maxny){
3109 small_to->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.375f));
3110 small_this->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.375f));
3111 }
3112
3113 small_this->process_inplace("xform.fourierorigin.tocenter");
3114
3115 // This is a solid estimate for very complete searching, 2.5 is a bit arbitrary
3116 // make sure the altitude step hits 90 degrees, not absolutely necessary for this, but can't hurt
3117// float astep = 89.999/floor(pi/(2.5*2.0*atan(2.0/ss)));
3118 float astep = 89.999/floor(pi/(1.5*2.0*atan(2.0/ss)));
3119 if ((maxang>0) && (astep>maxang/2)){
3120 astep=maxang/2;
3121 }
3122
3123 // This is drawn from single particle analysis testing, which in that case insures that enough sampling to
3124 // reasonably fill Fourier space is achieved, but doesn't perfectly apply to SPT
3125// float astep = (float)(89.99/ceil(90.0*9.0/(8.0*sqrt((float)(4300.0/ss))))); // 8 is (3+speed) from SPA with speed=5
3126
3127 // This insures we make at least one real effort at each level
3128 for (int i=0; i<nsoln; i++) {
3129 s_score[i]=1.0e24; // reset the scores since the different scales will not match
3130 if (fabs(s_step[i*3+0])<astep/4.0) s_step[i*3+0]*=2.0;
3131 if (fabs(s_step[i*3+1])<astep/4.0) s_step[i*3+1]*=2.0;
3132 if (fabs(s_step[i*3+2])<astep/4.0) s_step[i*3+2]*=2.0;
3133 }
3134
3135 // This is for the first loop, we do a full search in a heavily downsampled space
3136 if (curiter==0) {
3137 // Genrate points on a sphere in an asymmetric unit
3138 if (verbose>1) printf("stage 1 - ang step %1.2f\n",astep);
3139 Dict d;
3140 d["inc_mirror"] = 1;
3141 d["delta"] = astep;
3142 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
3143 // We don't generate for phi, since this can produce a very large number of orientations
3144 vector<Transform> transforms = sym->gen_orientations("eman", d);
3145 if (verbose>0) printf("%d orientations to test (%lu)\n",(int)(transforms.size()*(360.0/astep)),transforms.size());
3146
3147 for (unsigned int it=0; it<transforms.size(); it++) {
3148 if (verbose>2) {
3149 printf(" %d/%lu \r",it,transforms.size());
3150 fflush(stdout);
3151 }
3152 for (float phi=0; phi<360.0; phi+=astep) {
3153 Transform t = transforms[it];
3154 Dict aap=t.get_params("eman");
3155 aap["phi"]=phi;
3156 aap["tx"]=0;
3157 aap["ty"]=0;
3158 aap["tz"]=0;
3159 t.set_params(aap);
3160// t.invert();
3161// aap=t.get_params("eman");
3162
3163 // somewhat strangely, rotations are actually much more expensive than FFTs, so we use a CCF for translation
3164 EMData *stt=small_this->project("gauss_fft", Dict("transform", EMObject(&t), "returnfft", 1));
3165
3166// EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
3167 EMData *ccf=small_to->calc_ccf(stt);
3168 IntPoint ml=ccf->calc_max_location_wrap(maxshift,maxshift,0);
3169
3170 aap["tx"]=(int)ml[0];
3171 aap["ty"]=(int)ml[1];
3172 aap["tz"]=0;//(int)ml[2];
3173 t.set_params(aap);
3174
3175 delete stt;
3176 delete ccf;
3177 stt=small_this->project("gauss_fft", Dict("transform", EMObject(&t), "returnfft", 1));
3178 // stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
3179 // we have to do 1 slow transform here now that we have the translation
3180
3181// float sim=stt->cmp("ccc.tomo.thresh",small_to,Dict("sigmaimg",sigmathis,"sigmawith",sigmato));
3182// float sim=stt->cmp("ccc.tomo.thresh",small_to);
3183 float sim=stt->cmp("frc",small_to);
3184// float sim=stt->cmp("fsc.tomo.auto",small_to,Dict("sigmaimg",sigmathisv,"sigmawith",sigmatov));
3185// float sim=stt->cmp("fsc.tomo.auto",small_to);
3186
3187 // We want to make sure our starting points are somewhat separated from each other, so we replace any angles too close to an existing angle
3188 // If we find an existing 'best' angle within range, then we either replace it or skip
3189 int worst=-1;
3190 float worstv=1.0e20;
3191 for (int i=0; i<nsoln; i++) {
3192 if (s_score[i]>1.0e20) continue; // hasn't been set yet
3193 Transform tdif=s_xform[i].inverse();
3194 tdif=tdif*t;
3195 float adif=tdif.get_rotation("spin")["omega"];
3196 if (adif<astep*4) {
3197 worst=i;
3198// printf("= %1.3f\n",adif);
3199 }
3200 }
3201
3202 // if we weren't close to an existing angle, then we find the lowest current score and use that
3203 if (worst==-1) {
3204 // First we find the worst solution in the list of possible best solutions, or the first
3205 // solution which is currently "empty"
3206 for (int i=0; i<nsoln; i++) {
3207 if (s_score[i]>1.0e20) { worst=i; break; }
3208 if (s_score[i]<worstv) {worst=i; worstv=s_score[i];}
3209// if (s_score[i]<s_score[worst]) worst=i;
3210 }
3211 }
3212
3213 // If the current solution is better than the 'worst' of the previous solutions, then we
3214 // displace it. Note that there is no sorting performed here
3215 if (sim<s_score[worst]) {
3216 s_score[worst]=sim;
3217// s_coverage[worst]=1;//stt->get_attr("fft_overlap");
3218 s_xform[worst]=t;
3219 }
3220 delete stt;
3221 }
3222 }
3223 if (verbose>2) printf("\n");
3224
3225
3226 }
3227 // Once we have our initial list of best locations, we just refine each possibility individually
3228 else {
3229 // We generate a search pattern around each existing solution
3230 if (verbose>1) printf("stage 2 (%1.2f)\n",astep);
3231 for (int i=0; i<nsoln; i++) {
3232
3233 if (verbose>2) {
3234 printf(" %d\t%d\r",i,nsoln);
3235 fflush(stdout);
3236 }
3237 // We work an axis at a time until we get where we want to be. Somewhat like a simplex
3238 int changed=1;
3239 Dict upd0;
3240 int r=testort(small_this,small_to,s_score,s_xform,i,upd0,initxf,maxshift, maxang);
3241 while (changed) {
3242 changed=0;
3243 for (int axis=0; axis<3; axis++) {
3244 if (fabs(s_step[i*3+axis])<astep/4.0) continue; // skip axes where we already have enough precision on this axis
3245 Dict upd;
3246 upd[axname[axis]]=s_step[i*3+axis];
3247 // when moving az, we move phi in the opposite direction by the same amount since the two are singular at alt=0
3248 // phi continues to move independently. I believe this should produce a more monotonic energy surface
3249 if (axis==0) upd[axname[2]]=-s_step[i*3+axis];
3250
3251 int r=testort(small_this,small_to,s_score,s_xform,i,upd,initxf,maxshift, maxang);
3252
3253 // If we fail, we reverse direction with a slightly smaller step and try that
3254 // Whether this fails or not, we move on to the next axis
3255 if (r) changed=1;
3256 else {
3257 s_step[i*3+axis]*=-0.75;
3258 upd[axname[axis]]=s_step[i*3+axis];
3259 r=testort(small_this,small_to,s_score,s_xform,i,upd,initxf,maxshift, maxang);
3260 if (r) changed=1;
3261 }
3262 if (verbose>4) printf("\nX %1.3f\t%1.3f\t%1.3f\t%d\t",s_step[i*3],s_step[i*3+1],s_step[i*3+2],changed);
3263 }
3264 if (verbose>3) {
3265 Dict aap=s_xform[i].get_params("eman");
3266 printf("\n%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t(%1.3f)",s_step[i*3],s_step[i*3+1],s_step[i*3+2],float(aap["az"]),float(aap["alt"]),float(aap["phi"]),s_score[i]);
3267 }
3268
3269 if (!changed) {
3270// for (int j=0; j<3; j++) s_step[i*3+j]*-0.75;
3271 changed=1;
3272 }
3273 if (fabs(s_step[i*3])<astep/4 && fabs(s_step[i*3+1])<astep/4 && fabs(s_step[i*3+2])<astep/4) changed=0;
3274 }
3275
3276 // Ouch, exhaustive (local) search
3277// for (int daz=-1; daz<=1; daz++) {
3278// for (int dalt=-1; dalt<=1; dalt++) {
3279// for (int dphi=-1; dphi<=1; dphi++) {
3280// Dict upd;
3281// upd["az"]=daz*astep;
3282// upd["alt"]=dalt*astep;
3283// upd["phi"]=dphi*astep;
3284// int r=testort(small_this,small_to,s_score,s_coverage,s_xform,i,upd);
3285// }
3286// }
3287// }
3288 }
3289 }
3290 // lazy earlier in defining s_ vectors, so lazy here too and inefficiently sorting
3291 // We are sorting inside the outermost loop so we can decrease the number of solutions
3292 // before we get to the finest precision
3293 for (unsigned int i=0; i<nsoln-1; i++) {
3294 for (unsigned int j=i+1; j<nsoln; j++) {
3295 if (s_score[i]>s_score[j]) {
3296 float t=s_score[i]; s_score[i]=s_score[j]; s_score[j]=t;
3297 Transform tt=s_xform[i]; s_xform[i]=s_xform[j]; s_xform[j]=tt;
3298 }
3299 }
3300 }
3301 for (unsigned int i=0; i<nsoln-1; i++) {
3302 Transform tt=s_xform[i];
3303 tt.set_trans(tt.get_trans()*(float)ny/(float)ss);
3304 s_xform[i]=tt;
3305 }
3306
3307 // At each level of sampling we (potentially) decrease the number of answers we check in detail
3308 // assuming we are gradually homing in on the best solution
3309 nsoln/=2;
3310 if (nsoln<nrsoln) nsoln=nrsoln;
3311
3312
3313 delete small_this;
3314 delete small_to;
3315 if (ss>=maxny && curiter>0) break;
3316 }
3317
3318 delete base_this;
3319 delete base_to;
3320
3321
3322 // initialize results
3323 vector<Dict> solns;
3324 for (unsigned int i = 0; i < nrsoln; ++i ) {
3325 Dict d;
3326 d["score"] = s_score[i];
3327// d["coverage"] = s_coverage[i];
3328// s_xform[i].invert(); // this is because we inverted the order of the input images above to match convention
3329 d["xform.align3d"] = &s_xform[i];
3330 d["xform.projection"] = &s_xform[i]; // we actually return the transform that project the 3D reference to 2D image.
3331 solns.push_back(d);
3332 }
3333
3334 return solns;
3335}
3336
3337// This is just to prevent redundancy. It takes the existing solution vectors as arguments, an a proposed update for
3338// vector i. It updates the vectors if the proposal makes an improvement, in which case it returns true
3339bool RT2Dto3DTreeAligner::testort(EMData *small_this,EMData *small_to,vector<float> &s_score, vector<Transform> &s_xform,int i,Dict &upd, Transform initxf, int maxshift, int maxang) const {
3340 Transform t;
3341 Dict aap=s_xform[i].get_params("eman");
3342
3343 int ny=small_this->get_ysize();
3344 if (params.has_key("initxform")){
3345 // when doing refinement, search around the given position
3346 const vector< Transform > xfs=params["initxform"];
3347 Dict xf0=xfs[i].get_params("eman");
3348 aap["tx"]=(float)xf0["tx"]*ny/(float)params["boxsize"];
3349 aap["ty"]=(float)xf0["ty"]*ny/(float)params["boxsize"];
3350
3351 }
3352 else{
3353 aap["tx"]=0;
3354 aap["ty"]=0;
3355 }
3356 aap["tz"]=0;
3357 for (Dict::const_iterator p=upd.begin(); p!=upd.end(); p++) {
3358 aap[p->first]=(float)aap[p->first]+(float)p->second;
3359 }
3360
3361 t.set_params(aap);
3362
3363
3364 if (maxang>0){
3365
3366 Transform tmp=initxf * t.inverse();
3367 float r=tmp.get_params("spin")["omega"];
3368 if(r>maxang)
3369 return false;
3370
3371 }
3372 // rotate in Fourier space then use a CCF to find translation
3373// EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
3374 EMData *stt=small_this->project("gauss_fft", Dict("transform", EMObject(&t), "returnfft", 1));
3375 EMData *ccf=small_to->calc_ccf(stt);
3376
3377
3378 IntPoint ml=ccf->calc_max_location_wrap(maxshift,maxshift,0);
3379 aap["tx"]=(int)aap["tx"]+(int)ml[0];
3380 aap["ty"]=(int)aap["ty"]+(int)ml[1];
3381
3382 t.set_params(aap);
3383 stt->translate(t.get_trans());
3384// EMData *st2=small_this->project("gauss_fft", Dict("transform", EMObject(&t), "returnfft", 1));
3385
3386 float sim;
3387 sim=stt->cmp("frc",small_to, Dict("pmin",(float)ny*0.125f, "pmax", (float)ny*0.33f));
3388
3389 delete ccf;
3390 delete stt;
3391
3392 // If the score is better than before, we update this particular best value
3393 if (sim<s_score[i]) {
3394 s_score[i]=sim;
3395// s_coverage[i]=1;//st2->get_attr("fft_overlap");
3396 s_xform[i]=t;
3397
3398 return true;
3399 }
3400 return false;
3401}
3402
3403EMData* RT3DTreeAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
3404{
3405
3406 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
3407
3408 Dict t;
3409 Transform* tr = (Transform*) alis[0]["xform.align3d"];
3410 t["transform"] = tr;
3411 EMData* soln = this_img->process("xform",t);
3412 soln->set_attr("xform.align3d",tr);
3413 soln->set_attr("score",alis[0]["score"]);
3414 soln->set_attr("coverage",alis[0]["coverage"]);
3415
3416 return soln;
3417
3418}
3419
3420// NOTE - if symmetry is applied, it is critical that "this" (as passed in to the algorithm) be the volume which is already aligned
3421// to the symmetry axes (ie - the reference). this is confusing as it is inverted internally in the algorithm, so the passed in
3422// "this" becomes "to" in the code to conform to the way the other aligners work.
3423vector<Dict> RT3DTreeAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nrsoln, const string & cmp_name, const Dict& cmp_params) const {
3424 if (nrsoln == 0) throw InvalidParameterException("ERROR (RT3DTreeAligner): nsoln must be >0"); // What was the user thinking?
3425
3426 int nsoln = nrsoln*2;
3427 if (nrsoln<32) nsoln=64; // we start with at least 32 solutions, but then gradually decrease with increasing scale
3428
3429 float sigmathis = params.set_default("sigmathis",0.01f);
3430 float sigmato = params.set_default("sigmato",0.01f);
3431 int verbose = params.set_default("verbose",0);
3432 float maxres = params.set_default("maxres",-1.0f);
3433 EMData *mask = params.set_default("mask",(EMData *)0);
3434
3435 // !!!!!! IMPORTANT NOTE - we are inverting the order of this and to here to match convention in other aligners, to compensate
3436 // the Transform is inverted before being returned
3437 EMData *base_this;
3438 EMData *base_to;
3439 if (mask) {
3440 if (this_img->is_complex()) {
3441 EMData *tmp = this_img->do_ift();
3442 tmp->process_inplace("xform.phaseorigin.tocenter");
3443 tmp->process_inplace("normalize.mask",Dict("mask",mask,"apply_mask",1));
3444 base_to=tmp->do_fft();
3445 base_to->process_inplace("xform.phaseorigin.tocorner");
3446 delete tmp;
3447 }
3448 else {
3449 EMData *tmp = this_img->process("normalize.mask",Dict("mask",mask,"apply_mask",1));
3450 base_to=tmp->do_fft();
3451 base_to->process_inplace("xform.phaseorigin.tocorner");
3452 delete tmp;
3453 }
3454
3455 if (to->is_complex()) {
3456 base_this=to->copy();
3457 }
3458 else {
3459 base_this=to->do_fft();
3460 base_this->process_inplace("xform.phaseorigin.tocorner");
3461 }
3462 }
3463 else {
3464 if (this_img->is_complex()) base_to=this_img->copy();
3465 else {
3466 base_to=this_img->do_fft();
3467 base_to->process_inplace("xform.phaseorigin.tocorner");
3468 }
3469
3470 if (to->is_complex()) base_this=to->copy();
3471 else {
3472 base_this=to->do_fft();
3473 base_this->process_inplace("xform.phaseorigin.tocorner");
3474 }
3475 }
3476
3477
3478 if (base_this->get_xsize()!=base_this->get_ysize()+2 || base_this->get_ysize()!=base_this->get_zsize()
3479 || base_to->get_xsize()!=base_to->get_ysize()+2 || base_to->get_ysize()!=base_to->get_zsize()) throw InvalidCallException("ERROR (RT3DTreeAligner): requires cubic images with even numbered box sizes");
3480
3481// base_this->process_inplace("mask.wedgefill", Dict("thresh_sigma", sigmathis));
3482// base_to->process_inplace("mask.wedgefill", Dict("thresh_sigma", sigmato));
3483
3484
3485 base_this->process_inplace("xform.fourierorigin.tocenter"); // easier to chop out Fourier subvolumes
3486 base_to->process_inplace("xform.fourierorigin.tocenter");
3487
3488 float apix=(float)this_img->get_attr("apix_x");
3489 int ny=this_img->get_ysize();
3490 params["boxsize"]=ny;
3491
3492 int maxny=ny;
3493 if (0)//maxres>0)
3494 maxny=4*int(ny*apix/maxres/2+1);
3495 if (verbose>0)
3496 printf("\n\n*******\nmax resolution %1.2f, box size %d\n", maxres, maxny);
3497
3498
3499 Transform initxf;
3500
3501// int downsample=floor(ny/20); // Minimum shrunken box size is 20^3
3502
3503 vector<float> s_score(nsoln,0.0f);
3504 vector<float> s_coverage(nsoln,0.0f);
3505 vector<float> s_step(nsoln*3,7.5f);
3506 vector<Transform> s_xform(nsoln);
3507 if (verbose>0) printf("%d solutions\n",nsoln);
3508
3509
3510 int curiter=-1;
3511 int sexp_start=4;
3512 if (params.has_key("initxform")){
3513 const vector< Transform > xfs=params["initxform"];
3514 nsoln=xfs.size();
3515 initxf.set_params(xfs[0].get_params("eman"));
3516 for (unsigned int i=0; i<nsoln; i++){
3517 s_xform[i].set_params(xfs[i].get_params("eman"));
3518 }
3519 sexp_start=5;
3520 curiter=0; // bypass the first round of global orientation search
3521// for (int i=0; i<nsoln*3; i++) {
3522// s_step[i]=.5;
3523// }
3524 }
3525 else if (params.has_key("maxang")) { // if we got maxshift or maxang without initxform, we need to set up identity xform
3526 nsoln=1; // Default transform should be identity, so just set number to 1
3527 curiter=0; // skip global search
3528 sexp_start=5;
3529 }
3530
3531 int maxshift00=(int)params.set_default("maxshift",ny/4);
3532 float maxang=params.set_default("maxang",-1.0);
3533
3534// float dstep[3] = {7.5,7.5,7.5}; // we take steps for each of the 3 angles, may be positive or negative
3535 string axname[] = {"az","alt","phi"};
3536 int lastss=24;
3537 // We start with 32^3, 64^3 ...
3538 for (int sexp=sexp_start; sexp<10; sexp++) {
3539 curiter++;
3540// for (int sexp=4; sexp<5; sexp++) {
3541 int ss=pow(2.0,sexp);
3542 if (ss>maxny) ss=maxny;
3543 if (ss<24) ss=24; // 16 may be too small, but 32 takes too long...
3544 else if (ss<48) ss=48; // 16 may be too small, but 32 takes too long...
3545 if (verbose>0) printf("\nSize %d\n",ss);
3546
3547 int maxshift=maxshift00*ss/ny;
3548 if (maxshift00<0) maxshift=-1;
3549
3550 //ss=good_size(ny/ds);
3551 EMData *small_this=base_this->get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
3552 EMData *small_to= base_to-> get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
3553 small_this->process_inplace("xform.fourierorigin.tocorner");
3554 small_this->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
3555 small_to->process_inplace("xform.fourierorigin.tocorner");
3556 small_to->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
3557
3558 if (ss<maxny){
3559 // skip lp filter at full sampling seems to help..
3560 small_this->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.33f));
3561 small_to->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.33f));
3562 }
3563
3564 // these are cached for speed in the comparator
3565 vector<float>sigmathisv=small_this->calc_radial_dist(ss/2,0,1,4);
3566 vector<float>sigmatov=small_to->calc_radial_dist(ss/2,0,1,4);
3567 for (int i=0; i<ss/2; i++) {
3568 sigmathisv[i]*=sigmathisv[i]*sigmathis;
3569 sigmatov[i]*=sigmatov[i]*sigmato;
3570 }
3571
3572 // This is a solid estimate for very complete searching, 2.5 is a bit arbitrary
3573 // make sure the altitude step hits 90 degrees, not absolutely necessary for this, but can't hurt
3574// float astep = 89.999/floor(pi/(2.5*2.0*atan(2.0/ss)));
3575 float astep = 89.999/floor(pi/(1.5*2.0*atan(2.0/ss)));
3576
3577 // This is drawn from single particle analysis testing, which in that case insures that enough sampling to
3578 // reasonably fill Fourier space is achieved, but doesn't perfectly apply to SPT
3579// float astep = (float)(89.99/ceil(90.0*9.0/(8.0*sqrt((float)(4300.0/ss))))); // 8 is (3+speed) from SPA with speed=5
3580
3581 // This insures we make at least one real effort at each level
3582 for (int i=0; i<s_score.size(); i++) {
3583 s_score[i]=1.0e24; // reset the scores since the different scales will not match
3584 if (fabs(s_step[i*3+0])<astep/4.0) s_step[i*3+0]*=2.0;
3585 if (fabs(s_step[i*3+1])<astep/4.0) s_step[i*3+1]*=2.0;
3586 if (fabs(s_step[i*3+2])<astep/4.0) s_step[i*3+2]*=2.0;
3587 }
3588
3589 // This is for the first loop, we do a full search in a heavily downsampled space
3590 if (curiter==0){ //s_coverage[0]==0.0f) {
3591 // Genrate points on a sphere in an asymmetric unit
3592 if (verbose>1) printf("stage 1 - ang step %1.2f\n",astep);
3593 Dict d;
3594 d["inc_mirror"] = true;
3595 d["delta"] = astep;
3596 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
3597 // We don't generate for phi, since this can produce a very large number of orientations
3598 vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
3599 if (verbose>0) printf("%d orientations to test (%lu)\n",(int)(transforms.size()*(360.0/astep)),transforms.size());
3600 if (transforms.size()<30) continue; // for very high symmetries we will go up to 32 instead of 24
3601
3602 // We iterate over all orientations in an asym triangle (alt & az) then deal with phi ourselves
3603// for (std::vector<Transform>::iterator t = transforms.begin(); t!=transforms.end(); ++t) { // iterator form was causing all sorts of problems
3604 for (unsigned int it=0; it<transforms.size(); it++) {
3605
3606 if (verbose>2) {
3607 printf(" %d/%lu \r",it,transforms.size());
3608 fflush(stdout);
3609 }
3610 for (float phi=0; phi<360.0; phi+=astep) {
3611 Transform t = transforms[it];
3612 Dict aap=t.get_params("eman");
3613 aap["phi"]=phi;
3614 aap["tx"]=0;
3615 aap["ty"]=0;
3616 aap["tz"]=0;
3617 t.set_params(aap);
3618 t.invert();
3619 aap=t.get_params("eman");
3620
3621 // somewhat strangely, rotations are actually much more expensive than FFTs, so we use a CCF for translation
3622 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
3623 EMData *ccf=small_to->calc_ccf(stt);
3624 IntPoint ml=ccf->calc_max_location_wrap();
3625
3626 aap["tx"]=(int)ml[0];
3627 aap["ty"]=(int)ml[1];
3628 aap["tz"]=(int)ml[2];
3629 t.set_params(aap);
3630 delete stt;
3631 delete ccf;
3632 stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1)); // we have to do 1 slow transform here now that we have the translation
3633
3634// float sim=stt->cmp("ccc.tomo.thresh",small_to,Dict("sigmaimg",sigmathis,"sigmawith",sigmato));
3635 float sim=stt->cmp("ccc.tomo.thresh",small_to);
3636
3637// float sim=stt->cmp("fsc.tomo.auto",small_to,Dict("sigmaimg",sigmathisv,"sigmawith",sigmatov));
3638// float sim=stt->cmp("fsc.tomo.auto",small_to);
3639
3640 // We want to make sure our starting points are somewhat separated from each other, so we replace any angles too close to an existing angle
3641 // If we find an existing 'best' angle within range, then we either replace it or skip
3642 int worst=-1;
3643 float worstv=1.0e20;
3644 for (int i=0; i<nsoln; i++) {
3645 if (s_coverage[i]==0.0) continue; // hasn't been set yet
3646 Transform tdif=s_xform[i].inverse();
3647 tdif=tdif*t;
3648 float adif=tdif.get_rotation("spin")["omega"];
3649 if (adif<astep*2.5) {
3650 worst=i;
3651// printf("= %1.3f\n",adif);
3652 }
3653 }
3654
3655 // if we weren't close to an existing angle, then we find the lowest current score and use that
3656 if (worst==-1) {
3657 // First we find the worst solution in the list of possible best solutions, or the first
3658 // solution which is currently "empty"
3659 for (int i=0; i<nsoln; i++) {
3660 if (s_coverage[i]==0.0) { worst=i; break; }
3661 if (s_score[i]<worstv) {worst=i; worstv=s_score[i];}
3662 }
3663 }
3664
3665 // If the current solution is better than the 'worst' of the previous solutions, then we
3666 // displace it. Note that there is no sorting performed here
3667 if (sim<s_score[worst]) {
3668 s_score[worst]=sim;
3669 s_coverage[worst]=stt->get_attr("fft_overlap");
3670 s_xform[worst]=t;
3671 //printf("%f\t%f\t%d\n",s_score[worst],s_coverage[worst],worst);
3672 }
3673 delete stt;
3674 }
3675 }
3676 if (verbose>2) printf("\n");
3677// for (int i=0; i<nsoln; i++) {
3678// Dict d=s_xform[i].get_params("eman");
3679// printf("%d, %f, %f, %f, %f\n",i,(float)d["alt"], (float)d["az"], (float)d["phi"], s_score[i]);
3680// }
3681
3682 }
3683 // Once we have our initial list of best locations, we just refine each possibility individually
3684 else {
3685 // We generate a search pattern around each existing solution
3686 float res=float(ny)/float(ss)*apix*2;
3687 if (verbose>1) printf("stage 2 - maxres %1.2f, astep %1.2f\n",res, astep);
3688 if (ss>=maxny && params.has_key("initxform")){
3689 // last iteration, add back the initial transform so we do not do worse than the last run
3690 s_xform[nsoln]=initxf;
3691 nsoln+=1;
3692 }
3693
3694 for (int i=0; i<nsoln; i++) {
3695
3696 if (verbose>2) {
3697 printf(" %d\t%d\r",i,nsoln);
3698 fflush(stdout);
3699 }
3700 // We work an axis at a time until we get where we want to be. Somewhat like a simplex
3701 int changed=1;
3702 Dict upd;
3703 testort(small_this,small_to,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift,mask);
3704 while (changed) {
3705 changed=0;
3706 for (int axis=0; axis<3; axis++) {
3707 if (fabs(s_step[i*3+axis])<astep/4.0) continue; // skip axes where we already have enough precision on this axis
3708 upd[axname[axis]]=s_step[i*3+axis];
3709 // when moving az, we move phi in the opposite direction by the same amount since the two are singular at alt=0
3710 // phi continues to move independently. I believe this should produce a more monotonic energy surface
3711 if (axis==0) upd[axname[2]]=-s_step[i*3+axis];
3712
3713 int r=testort(small_this,small_to,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift,mask);
3714
3715 // If we fail, we reverse direction with a slightly smaller step and try that
3716 // Whether this fails or not, we move on to the next axis
3717 if (r) changed=1;
3718 else {
3719 s_step[i*3+axis]*=-0.75;
3720 upd[axname[axis]]=s_step[i*3+axis];
3721 r=testort(small_this,small_to,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift,mask);
3722 if (r) changed=1;
3723 }
3724 if (verbose>4) printf("\nX %1.3f\t%1.3f\t%1.3f\t%d\t",s_step[i*3],s_step[i*3+1],s_step[i*3+2],changed);
3725 }
3726 if (verbose>3) {
3727 Dict aap=s_xform[i].get_params("eman");
3728 printf("\n%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t(%1.3f)",s_step[i*3],s_step[i*3+1],s_step[i*3+2],float(aap["az"]),float(aap["alt"]),float(aap["phi"]),s_score[i]);
3729 }
3730
3731 if (!changed) {
3732// for (int j=0; j<3; j++) s_step[i*3+j]*-0.75;
3733 changed=1;
3734 }
3735 if (fabs(s_step[i*3])<astep/4 && fabs(s_step[i*3+1])<astep/4 && fabs(s_step[i*3+2])<astep/4) changed=0;
3736 }
3737
3738 // Ouch, exhaustive (local) search
3739// for (int daz=-1; daz<=1; daz++) {
3740// for (int dalt=-1; dalt<=1; dalt++) {
3741// for (int dphi=-1; dphi<=1; dphi++) {
3742// Dict upd;
3743// upd["az"]=daz*astep;
3744// upd["alt"]=dalt*astep;
3745// upd["phi"]=dphi*astep;
3746// int r=testort(small_this,small_to,s_score,s_coverage,s_xform,i,upd);
3747// }
3748// }
3749// }
3750 }
3751 }
3752 // lazy earlier in defining s_ vectors, so lazy here too and inefficiently sorting
3753 // We are sorting inside the outermost loop so we can decrease the number of solutions
3754 // before we get to the finest precision
3755 for (unsigned int i=0; i<nsoln-1; i++) {
3756 for (unsigned int j=i+1; j<nsoln; j++) {
3757 if (s_score[i]>s_score[j]) {
3758 float t=s_score[i]; s_score[i]=s_score[j]; s_score[j]=t;
3759 t=s_coverage[i]; s_coverage[i]=s_coverage[j]; s_coverage[j]=t;
3760 Transform tt=s_xform[i]; s_xform[i]=s_xform[j]; s_xform[j]=tt;
3761 }
3762 }
3763 }
3764
3765 // At each level of sampling we (potentially) decrease the number of answers we check in detail
3766 // assuming we are gradually homing in on the best solution
3767 nsoln/=2;
3768 if (nsoln<nrsoln) nsoln=nrsoln;
3769
3770
3771 delete small_this;
3772 delete small_to;
3773
3774 lastss=ss;
3775 if (ss>=maxny && curiter>0) break;
3776 }
3777
3778 delete base_this;
3779 delete base_to;
3780
3781 // note the translations are wrong when the alignment stops before ss==ny
3782 if (lastss<ny){
3783 for (unsigned int i = 0; i < nrsoln; ++i ) {
3784 Transform tt=s_xform[i];
3785 Vec3f v=tt.get_trans();
3786 v=v*ny/lastss;
3787 tt.set_trans(v);
3788 s_xform[i]=tt;
3789
3790 }
3791 }
3792
3793 // initialize results
3794 vector<Dict> solns;
3795 for (unsigned int i = 0; i < nrsoln; ++i ) {
3796 Dict d;
3797 d["score"] = s_score[i];
3798 d["coverage"] = s_coverage[i];
3799 s_xform[i].invert(); // this is because we inverted the order of the input images above to match convention
3800 d["xform.align3d"] = &s_xform[i];
3801 solns.push_back(d);
3802 }
3803
3804 return solns;
3805}
3806
3807// This is just to prevent redundancy. It takes the existing solution vectors as arguments, an a proposed update for
3808// vector i. It updates the vectors if the proposal makes an improvement, in which case it returns true
3809bool RT3DTreeAligner::testort(EMData *small_this,EMData *small_to,vector<float> &sigmathisv,vector<float> &sigmatov,vector<float> &s_score, vector<float> &s_coverage,vector<Transform> &s_xform,int i,Dict &upd, Transform initxf, int maxshift, EMData* mask) const {
3810 Transform t;
3811 Dict aap=s_xform[i].get_params("eman");
3812 aap["tx"]=0;
3813 aap["ty"]=0;
3814 aap["tz"]=0;
3815 for (Dict::const_iterator p=upd.begin(); p!=upd.end(); p++) {
3816 aap[p->first]=(float)aap[p->first]+(float)p->second;
3817 }
3818
3819 t.set_params(aap);
3820 float maxang=params.set_default("maxang",-1.0);
3821 bool randphi=params.set_default("randphi",false);
3822 bool rand180=params.set_default("rand180",false);
3823
3824 if (maxang>0){
3825
3826 Transform tmp=initxf * t.inverse();
3827
3828 if (randphi){
3829 Dict td=t.get_params("eman");
3830 Dict ti=initxf.get_params("eman");
3831 td["phi"]=ti["phi"];
3832 Transform tmp1=Transform(td);
3833 tmp=initxf * tmp1.inverse();
3834 }
3835
3836 float r=tmp.get_params("spin")["omega"];
3837 if (rand180){
3838 r=r>90?(180-r):r;
3839 }
3840 if(r>maxang)
3841 return false;
3842
3843 }
3844
3845 int ny=small_this->get_ysize();
3846 if (params.has_key("initxform")){
3847 // when doing refinement, search around the given position
3848 t.set_trans(initxf.get_trans()*ny/(float)params["boxsize"]);
3849 aap=t.get_params("eman");
3850 }
3851
3852 // rotate in Fourier space then use a CCF to find translation
3853 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
3854 EMData *ccf=small_to->calc_ccf(stt);
3855// EMData *ccf=small_to->calc_flcf(stt);
3856 IntPoint ml=ccf->calc_max_location_wrap(maxshift, maxshift, maxshift);
3857
3858
3859 aap["tx"]=int(aap["tx"])+(int)ml[0];
3860 aap["ty"]=int(aap["ty"])+(int)ml[1];
3861 aap["tz"]=int(aap["tz"])+(int)ml[2];
3862 t.set_params(aap);
3863 EMData *st2=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1)); // we have to do 1 slow transform here now that we have the translation
3864
3865// float maxres = params.set_default("maxres",0.0f);
3866// float minres = params.set_default("minres",0.0f);
3867
3868 float sim=st2->cmp("fsc.tomo.auto",small_to,Dict("sigmaimg",sigmathisv,"sigmawith",sigmatov, "pmin", 4, "pmax",ny/3));
3869// float sim=st2->cmp("ccc.tomo.thresh",small_to,Dict("sigmaimg",sigmathis,"sigmawith",sigmato));
3870// printf("\nTESTORT %6.1f %6.1f %6.1f\t%4d %4d %4d\t%1.5g\t%1.5g %d (%d)",
3871// float(aap["az"]),float(aap["alt"]),float(aap["phi"]),int(aap["tx"]),int(aap["ty"]),int(aap["tz"]),sim,s_score[i],int(sim<s_score[i]),ccf->get_ysize());
3872
3873 delete ccf;
3874 // If the score is better than before, we update this particular best value
3875 if (sim<s_score[i]) {
3876 s_score[i]=sim;
3877 s_coverage[i]=st2->get_attr("fft_overlap");
3878 s_xform[i]=t;
3879 delete stt;
3880 delete st2;
3881 return true;
3882 }
3883 delete stt;
3884 delete st2;
3885 return false;
3886}
3887
3888EMData* RT3DLocalTreeAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
3889{
3890
3891 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
3892
3893 Dict t;
3894 Transform* tr = (Transform*) alis[0]["xform.align3d"];
3895 t["transform"] = tr;
3896 EMData* soln = this_img->process("xform",t);
3897 soln->set_attr("xform.align3d",tr);
3898 soln->set_attr("score",alis[0]["score"]);
3899 soln->set_attr("coverage",alis[0]["coverage"]);
3900
3901 return soln;
3902
3903}
3904
3905// NOTE - if symmetry is applied, it is critical that "this" (as passed in to the algorithm) be the volume which is already aligned
3906// to the symmetry axes (ie - the reference). this is confusing as it is inverted internally in the algorithm, so the passed in
3907// "this" becomes "to" in the code to conform to the way the other aligners work.
3908vector<Dict> RT3DLocalTreeAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nrsoln, const string & cmp_name, const Dict& cmp_params) const {
3909 if (nrsoln == 0) throw InvalidParameterException("ERROR (RT3DTreeAligner): nsoln must be >0"); // What was the user thinking?
3910
3911 int nsoln = nrsoln*2;
3912 if (nrsoln<32) nsoln=64; // we start with at least 32 solutions, but then gradually decrease with increasing scale
3913
3914 // !!!!!! IMPORTANT NOTE - we are inverting the order of this and to here to match convention in other aligners, to compensate
3915 // the Transform is inverted before being returned
3916 EMData *base_this;
3917 EMData *base_to;
3918 EMData *base_thissq;
3919 EMData *base_mask;
3920
3921 if (this_img->is_complex()) {
3922 base_to=this_img->copy();
3923 EMData *tmp = base_to->do_ift();
3924 tmp->process_inplace("threshold.notzero");
3925 base_mask=tmp->do_fft();
3926 delete tmp;
3927 }
3928 else {
3929 base_to=this_img->do_fft();
3930 base_to->process_inplace("xform.phaseorigin.tocorner");
3931 EMData *tmp = this_img->process("threshold.notzero");
3932 base_mask=tmp->do_fft();
3933 delete tmp;
3934 }
3935
3936 if (to->is_complex()) {
3937 base_this=to->copy();
3938 EMData *tmp = base_this->do_ift();
3939 tmp->process_inplace("math.squared");
3940 base_thissq = tmp->do_fft();
3941 delete tmp;
3942 }
3943 else {
3944 base_this=to->do_fft();
3945 base_this->process_inplace("xform.phaseorigin.tocorner");
3946 EMData *tmp = base_this->process("math.squared");
3947 base_thissq = tmp->do_fft();
3948 delete tmp;
3949 }
3950
3951 float sigmathis = params.set_default("sigmathis",0.01f);
3952 float sigmato = params.set_default("sigmato",0.01f);
3953 int verbose = params.set_default("verbose",0);
3954 float maxres = params.set_default("maxres",-1.0f);
3955
3956 if (base_this->get_xsize()!=base_this->get_ysize()+2 || base_this->get_ysize()!=base_this->get_zsize()
3957 || base_to->get_xsize()!=base_to->get_ysize()+2 || base_to->get_ysize()!=base_to->get_zsize()) throw InvalidCallException("ERROR (RT3DTreeAligner): requires cubic images with even numbered box sizes");
3958
3959// base_this->process_inplace("mask.wedgefill", Dict("thresh_sigma", sigmathis));
3960// base_to->process_inplace("mask.wedgefill", Dict("thresh_sigma", sigmato));
3961
3962
3963 base_this->process_inplace("xform.fourierorigin.tocenter"); // easier to chop out Fourier subvolumes
3964 base_to->process_inplace("xform.fourierorigin.tocenter");
3965 base_thissq->process_inplace("xform.fourierorigin.tocenter");
3966 base_mask->process_inplace("xform.fourierorigin.tocenter");
3967
3968 float apix=(float)this_img->get_attr("apix_x");
3969 int ny=this_img->get_ysize();
3970 params["boxsize"]=ny;
3971 int maxshift00=(int)params.set_default("maxshift",ny/4);
3972
3973 int maxny=ny;
3974 if (0)//maxres>0)
3975 maxny=4*int(ny*apix/maxres/2+1);
3976 if (verbose>0)
3977 printf("\n\n*******\nmax resolution %1.2f, box size %d\n", maxres, maxny);
3978
3979 float maxang=params.set_default("maxang",-1.0);
3980 Transform initxf;
3981
3982// int downsample=floor(ny/20); // Minimum shrunken box size is 20^3
3983
3984 vector<float> s_score(nsoln,0.0f);
3985 vector<float> s_coverage(nsoln,0.0f);
3986 vector<float> s_step(nsoln*3,7.5f);
3987 vector<Transform> s_xform(nsoln);
3988 if (verbose>0) printf("%d solutions\n",nsoln);
3989
3990
3991 int curiter=-1;
3992 int sexp_start=4;
3993 if (params.has_key("initxform")){
3994 const vector< Transform > xfs=params["initxform"];
3995 nsoln=xfs.size();
3996 initxf.set_params(xfs[0].get_params("eman"));
3997 for (unsigned int i=0; i<nsoln; i++){
3998 s_xform[i].set_params(xfs[i].get_params("eman"));
3999 }
4000 sexp_start=5;
4001 curiter=0;
4002// for (int i=0; i<nsoln*3; i++) {
4003// s_step[i]=.5;
4004// }
4005 }
4006
4007
4008// float dstep[3] = {7.5,7.5,7.5}; // we take steps for each of the 3 angles, may be positive or negative
4009 string axname[] = {"az","alt","phi"};
4010 int lastss=24;
4011 // We start with 32^3, 64^3 ...
4012 for (int sexp=sexp_start; sexp<10; sexp++) {
4013 curiter++;
4014// for (int sexp=4; sexp<5; sexp++) {
4015 int ss=pow(2.0,sexp);
4016 if (ss>maxny) ss=maxny;
4017 if (ss<24) ss=24; // 16 may be too small, but 32 takes too long...
4018 else if (ss<48) ss=48; // 16 may be too small, but 32 takes too long...
4019 if (verbose>0) printf("\nSize %d\n",ss);
4020
4021 int maxshift=maxshift00*ss/ny;
4022 if (maxshift00<0) maxshift=-1;
4023
4024 //ss=good_size(ny/ds);
4025 EMData *small_this=base_this->get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
4026 EMData *small_to= base_to-> get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
4027 EMData *small_thissq=base_thissq->get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
4028 EMData *small_mask= base_mask-> get_clip(Region(0,(ny-ss)/2,(ny-ss)/2,ss+2,ss,ss));
4029 small_this->process_inplace("xform.fourierorigin.tocorner");
4030 small_this->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
4031 small_to->process_inplace("xform.fourierorigin.tocorner");
4032 small_to->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
4033 small_mask->process_inplace("xform.fourierorigin.tocorner");
4034 small_mask->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
4035 small_thissq->process_inplace("xform.fourierorigin.tocorner");
4036 small_thissq->process_inplace("filter.highpass.gauss",Dict("cutoff_pixels",4));
4037
4038 if (ss<maxny){
4039 // skip lp filter at full sampling seems to help..
4040 small_this->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.33f));
4041 small_to->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.33f));
4042 }
4043
4044 // these are cached for speed in the comparator
4045 vector<float>sigmathisv=small_this->calc_radial_dist(ss/2,0,1,4);
4046 vector<float>sigmatov=small_to->calc_radial_dist(ss/2,0,1,4);
4047 for (int i=0; i<ss/2; i++) {
4048 sigmathisv[i]*=sigmathisv[i]*sigmathis;
4049 sigmatov[i]*=sigmatov[i]*sigmato;
4050 }
4051
4052 // This is a solid estimate for very complete searching, 2.5 is a bit arbitrary
4053 // make sure the altitude step hits 90 degrees, not absolutely necessary for this, but can't hurt
4054// float astep = 89.999/floor(pi/(2.5*2.0*atan(2.0/ss)));
4055 float astep = 89.999/floor(pi/(1.5*2.0*atan(2.0/ss)));
4056
4057 // This is drawn from single particle analysis testing, which in that case insures that enough sampling to
4058 // reasonably fill Fourier space is achieved, but doesn't perfectly apply to SPT
4059// float astep = (float)(89.99/ceil(90.0*9.0/(8.0*sqrt((float)(4300.0/ss))))); // 8 is (3+speed) from SPA with speed=5
4060
4061 // This insures we make at least one real effort at each level
4062 for (int i=0; i<s_score.size(); i++) {
4063 s_score[i]=1.0e24; // reset the scores since the different scales will not match
4064 if (fabs(s_step[i*3+0])<astep/4.0) s_step[i*3+0]*=2.0;
4065 if (fabs(s_step[i*3+1])<astep/4.0) s_step[i*3+1]*=2.0;
4066 if (fabs(s_step[i*3+2])<astep/4.0) s_step[i*3+2]*=2.0;
4067 }
4068
4069 // This is for the first loop, we do a full search in a heavily downsampled space
4070 if (curiter==0){ //s_coverage[0]==0.0f) {
4071 // Genrate points on a sphere in an asymmetric unit
4072 if (verbose>1) printf("stage 1 - ang step %1.2f\n",astep);
4073 Dict d;
4074 d["inc_mirror"] = true;
4075 d["delta"] = astep;
4076 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
4077 // We don't generate for phi, since this can produce a very large number of orientations
4078 vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
4079 if (verbose>0) printf("%d orientations to test (%lu)\n",(int)(transforms.size()*(360.0/astep)),transforms.size());
4080 if (transforms.size()<30) continue; // for very high symmetries we will go up to 32 instead of 24
4081
4082 // We iterate over all orientations in an asym triangle (alt & az) then deal with phi ourselves
4083// for (std::vector<Transform>::iterator t = transforms.begin(); t!=transforms.end(); ++t) { // iterator form was causing all sorts of problems
4084 for (unsigned int it=0; it<transforms.size(); it++) {
4085
4086 if (verbose>2) {
4087 printf(" %d/%lu \r",it,transforms.size());
4088 fflush(stdout);
4089 }
4090 for (float phi=0; phi<360.0; phi+=astep) {
4091 Transform t = transforms[it];
4092 Dict aap=t.get_params("eman");
4093 aap["phi"]=phi;
4094 aap["tx"]=0;
4095 aap["ty"]=0;
4096 aap["tz"]=0;
4097 t.set_params(aap);
4098 t.invert();
4099 aap=t.get_params("eman");
4100
4101 // somewhat strangely, rotations are actually much more expensive than FFTs, so we use a CCF for translation
4102 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
4103 EMData *sttsq=small_thissq->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
4104// EMData *ccf=small_to->calc_ccf(stt);
4105 EMData *ccf=small_to->calc_ccf_masked(stt,sttsq,small_mask);
4106 IntPoint ml=ccf->calc_max_location_wrap();
4107
4108 aap["tx"]=(int)ml[0];
4109 aap["ty"]=(int)ml[1];
4110 aap["tz"]=(int)ml[2];
4111 t.set_params(aap);
4112 delete stt;
4113 delete sttsq;
4114 delete ccf;
4115 stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1)); // we have to do 1 slow transform here now that we have the translation
4116
4117// float sim=stt->cmp("ccc.tomo.thresh",small_to,Dict("sigmaimg",sigmathis,"sigmawith",sigmato));
4118 float sim=stt->cmp("ccc.tomo.thresh",small_to);
4119
4120// float sim=stt->cmp("fsc.tomo.auto",small_to,Dict("sigmaimg",sigmathisv,"sigmawith",sigmatov));
4121// float sim=stt->cmp("fsc.tomo.auto",small_to);
4122
4123 // We want to make sure our starting points are somewhat separated from each other, so we replace any angles too close to an existing angle
4124 // If we find an existing 'best' angle within range, then we either replace it or skip
4125 int worst=-1;
4126 float worstv=1.0e20;
4127 for (int i=0; i<nsoln; i++) {
4128 if (s_coverage[i]==0.0) continue; // hasn't been set yet
4129 Transform tdif=s_xform[i].inverse();
4130 tdif=tdif*t;
4131 float adif=tdif.get_rotation("spin")["omega"];
4132 if (adif<astep*2.5) {
4133 worst=i;
4134// printf("= %1.3f\n",adif);
4135 }
4136 }
4137
4138 // if we weren't close to an existing angle, then we find the lowest current score and use that
4139 if (worst==-1) {
4140 // First we find the worst solution in the list of possible best solutions, or the first
4141 // solution which is currently "empty"
4142 for (int i=0; i<nsoln; i++) {
4143 if (s_coverage[i]==0.0) { worst=i; break; }
4144 if (s_score[i]<worstv) {worst=i; worstv=s_score[i];}
4145 }
4146 }
4147
4148 // If the current solution is better than the 'worst' of the previous solutions, then we
4149 // displace it. Note that there is no sorting performed here
4150 if (sim<s_score[worst]) {
4151 s_score[worst]=sim;
4152 s_coverage[worst]=stt->get_attr("fft_overlap");
4153 s_xform[worst]=t;
4154 //printf("%f\t%f\t%d\n",s_score[worst],s_coverage[worst],worst);
4155 }
4156 delete stt;
4157 }
4158 }
4159 if (verbose>2) printf("\n");
4160// for (int i=0; i<nsoln; i++) {
4161// Dict d=s_xform[i].get_params("eman");
4162// printf("%d, %f, %f, %f, %f\n",i,(float)d["alt"], (float)d["az"], (float)d["phi"], s_score[i]);
4163// }
4164
4165 }
4166 // Once we have our initial list of best locations, we just refine each possibility individually
4167 else {
4168 // We generate a search pattern around each existing solution
4169 float res=float(ny)/float(ss)*apix*2;
4170 if (verbose>1) printf("stage 2 - maxres %1.2f, astep %1.2f\n",res, astep);
4171 if (ss>=maxny && params.has_key("initxform")){
4172 // last iteration, add back the initial transform so we do not do worse than the last run
4173 s_xform[nsoln]=initxf;
4174 nsoln+=1;
4175 }
4176
4177 for (int i=0; i<nsoln; i++) {
4178
4179 if (verbose>2) {
4180 printf(" %d\t%d\r",i,nsoln);
4181 fflush(stdout);
4182 }
4183 // We work an axis at a time until we get where we want to be. Somewhat like a simplex
4184 int changed=1;
4185 Dict upd;
4186 testort(small_this,small_to,small_mask,small_thissq,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift);
4187 while (changed) {
4188 changed=0;
4189 for (int axis=0; axis<3; axis++) {
4190 if (fabs(s_step[i*3+axis])<astep/4.0) continue; // skip axes where we already have enough precision on this axis
4191 upd[axname[axis]]=s_step[i*3+axis];
4192 // when moving az, we move phi in the opposite direction by the same amount since the two are singular at alt=0
4193 // phi continues to move independently. I believe this should produce a more monotonic energy surface
4194 if (axis==0) upd[axname[2]]=-s_step[i*3+axis];
4195
4196 int r=testort(small_this,small_to,small_mask,small_thissq,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift);
4197
4198 // If we fail, we reverse direction with a slightly smaller step and try that
4199 // Whether this fails or not, we move on to the next axis
4200 if (r) changed=1;
4201 else {
4202 s_step[i*3+axis]*=-0.75;
4203 upd[axname[axis]]=s_step[i*3+axis];
4204 r=testort(small_this,small_to,small_mask,small_thissq,sigmathisv,sigmatov,s_score,s_coverage,s_xform,i,upd, initxf, maxshift);
4205 if (r) changed=1;
4206 }
4207 if (verbose>4) printf("\nX %1.3f\t%1.3f\t%1.3f\t%d\t",s_step[i*3],s_step[i*3+1],s_step[i*3+2],changed);
4208 }
4209 if (verbose>3) {
4210 Dict aap=s_xform[i].get_params("eman");
4211 printf("\n%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t%1.3f\t(%1.3f)",s_step[i*3],s_step[i*3+1],s_step[i*3+2],float(aap["az"]),float(aap["alt"]),float(aap["phi"]),s_score[i]);
4212 }
4213
4214 if (!changed) {
4215// for (int j=0; j<3; j++) s_step[i*3+j]*-0.75;
4216 changed=1;
4217 }
4218 if (fabs(s_step[i*3])<astep/4 && fabs(s_step[i*3+1])<astep/4 && fabs(s_step[i*3+2])<astep/4) changed=0;
4219 }
4220
4221 // Ouch, exhaustive (local) search
4222// for (int daz=-1; daz<=1; daz++) {
4223// for (int dalt=-1; dalt<=1; dalt++) {
4224// for (int dphi=-1; dphi<=1; dphi++) {
4225// Dict upd;
4226// upd["az"]=daz*astep;
4227// upd["alt"]=dalt*astep;
4228// upd["phi"]=dphi*astep;
4229// int r=testort(small_this,small_to,s_score,s_coverage,s_xform,i,upd);
4230// }
4231// }
4232// }
4233 }
4234 }
4235 // lazy earlier in defining s_ vectors, so lazy here too and inefficiently sorting
4236 // We are sorting inside the outermost loop so we can decrease the number of solutions
4237 // before we get to the finest precision
4238 for (unsigned int i=0; i<nsoln-1; i++) {
4239 for (unsigned int j=i+1; j<nsoln; j++) {
4240 if (s_score[i]>s_score[j]) {
4241 float t=s_score[i]; s_score[i]=s_score[j]; s_score[j]=t;
4242 t=s_coverage[i]; s_coverage[i]=s_coverage[j]; s_coverage[j]=t;
4243 Transform tt=s_xform[i]; s_xform[i]=s_xform[j]; s_xform[j]=tt;
4244 }
4245 }
4246 }
4247
4248 // At each level of sampling we (potentially) decrease the number of answers we check in detail
4249 // assuming we are gradually homing in on the best solution
4250 nsoln/=2;
4251 if (nsoln<nrsoln) nsoln=nrsoln;
4252
4253
4254 delete small_this;
4255 delete small_to;
4256 delete small_thissq;
4257 delete small_mask;
4258
4259 lastss=ss;
4260 if (ss>=maxny && curiter>0) break;
4261 }
4262
4263 delete base_this;
4264 delete base_to;
4265 delete base_thissq;
4266 delete base_mask;
4267
4268 // note the translations are wrong when the alignment stops before ss==ny
4269 if (lastss<ny){
4270 for (unsigned int i = 0; i < nrsoln; ++i ) {
4271 Transform tt=s_xform[i];
4272 Vec3f v=tt.get_trans();
4273 v=v*ny/lastss;
4274 tt.set_trans(v);
4275 s_xform[i]=tt;
4276
4277 }
4278 }
4279
4280 // initialize results
4281 vector<Dict> solns;
4282 for (unsigned int i = 0; i < nrsoln; ++i ) {
4283 Dict d;
4284 d["score"] = s_score[i];
4285 d["coverage"] = s_coverage[i];
4286 s_xform[i].invert(); // this is because we inverted the order of the input images above to match convention
4287 d["xform.align3d"] = &s_xform[i];
4288 solns.push_back(d);
4289 }
4290
4291 return solns;
4292}
4293
4294// This is just to prevent redundancy. It takes the existing solution vectors as arguments, an a proposed update for
4295// vector i. It updates the vectors if the proposal makes an improvement, in which case it returns true
4296bool RT3DLocalTreeAligner::testort(EMData *small_this,EMData *small_to,EMData *small_mask,EMData *small_thissq,vector<float> &sigmathisv,vector<float> &sigmatov,vector<float> &s_score, vector<float> &s_coverage,vector<Transform> &s_xform,int i,Dict &upd, Transform initxf, int maxshift) const {
4297 Transform t;
4298 Dict aap=s_xform[i].get_params("eman");
4299 aap["tx"]=0;
4300 aap["ty"]=0;
4301 aap["tz"]=0;
4302 for (Dict::const_iterator p=upd.begin(); p!=upd.end(); p++) {
4303 aap[p->first]=(float)aap[p->first]+(float)p->second;
4304 }
4305
4306 t.set_params(aap);
4307 float maxang=params.set_default("maxang",-1.0);
4308 bool randphi=params.set_default("randphi",false);
4309 bool rand180=params.set_default("rand180",false);
4310
4311 if (maxang>0){
4312
4313 Transform tmp=initxf * t.inverse();
4314
4315 if (randphi){
4316 Dict td=t.get_params("eman");
4317 Dict ti=initxf.get_params("eman");
4318 td["phi"]=ti["phi"];
4319 Transform tmp1=Transform(td);
4320 tmp=initxf * tmp1.inverse();
4321 }
4322
4323 float r=tmp.get_params("spin")["omega"];
4324 if (rand180){
4325 r=r>90?(180-r):r;
4326 }
4327 if(r>maxang)
4328 return false;
4329
4330 }
4331
4332 int ny=small_this->get_ysize();
4333 if (params.has_key("initxform")){
4334 // when doing refinement, search around the given position
4335 t.set_trans(initxf.get_trans()*ny/(float)params["boxsize"]);
4336 aap=t.get_params("eman");
4337 }
4338
4339 // rotate in Fourier space then use a CCF to find translation
4340 EMData *stt=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
4341 EMData *sttsq=small_thissq->process("xform",Dict("transform",EMObject(&t),"zerocorners",1));
4342// EMData *ccf=small_to->calc_ccf(stt);
4343 EMData *ccf=small_to->calc_ccf_masked(stt,sttsq,small_mask);
4344 IntPoint ml=ccf->calc_max_location_wrap(maxshift, maxshift, maxshift);
4345
4346
4347 aap["tx"]=int(aap["tx"])+(int)ml[0];
4348 aap["ty"]=int(aap["ty"])+(int)ml[1];
4349 aap["tz"]=int(aap["tz"])+(int)ml[2];
4350 t.set_params(aap);
4351 EMData *st2=small_this->process("xform",Dict("transform",EMObject(&t),"zerocorners",1)); // we have to do 1 slow transform here now that we have the translation
4352
4353// float maxres = params.set_default("maxres",0.0f);
4354// float minres = params.set_default("minres",0.0f);
4355
4356 float sim=st2->cmp("fsc.tomo.auto",small_to,Dict("sigmaimg",sigmathisv,"sigmawith",sigmatov, "pmin", 4, "pmax",ny/3));
4357// float sim=st2->cmp("ccc.tomo.thresh",small_to,Dict("sigmaimg",sigmathis,"sigmawith",sigmato));
4358// printf("\nTESTORT %6.1f %6.1f %6.1f\t%4d %4d %4d\t%1.5g\t%1.5g %d (%d)",
4359// float(aap["az"]),float(aap["alt"]),float(aap["phi"]),int(aap["tx"]),int(aap["ty"]),int(aap["tz"]),sim,s_score[i],int(sim<s_score[i]),ccf->get_ysize());
4360
4361 delete ccf;
4362 // If the score is better than before, we update this particular best value
4363 if (sim<s_score[i]) {
4364 s_score[i]=sim;
4365 s_coverage[i]=st2->get_attr("fft_overlap");
4366 s_xform[i]=t;
4367 delete stt;
4368 delete sttsq;
4369 delete st2;
4370 return true;
4371 }
4372 delete stt;
4373 delete sttsq;
4374 delete st2;
4375 return false;
4376}
4377
4378EMData* RT3DSphereAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
4379{
4380
4381 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
4382
4383 Dict t;
4384 Transform* tr = (Transform*) alis[0]["xform.align3d"];
4385 t["transform"] = tr;
4386 EMData* soln = this_img->process("xform",t);
4387 soln->set_attr("xform.align3d",tr);
4388 delete tr; tr = 0;
4389
4390 return soln;
4391
4392}
4393
4394vector<Dict> RT3DSphereAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const {
4395
4396 if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
4397 throw ImageDimensionException("This aligner only works for 3D images");
4398 }
4399
4400 int searchx = 0;
4401 int searchy = 0;
4402 int searchz = 0;
4403
4404 bool dotrans = params.set_default("dotrans",1);
4405 if (params.has_key("search")) {
4406 vector<string> check;
4407 check.push_back("searchx");
4408 check.push_back("searchy");
4409 check.push_back("searchz");
4410 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
4411 if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
4412 }
4413 int search = params["search"];
4414 searchx = search;
4415 searchy = search;
4416 searchz = search;
4417 } else {
4418 searchx = params.set_default("searchx",3);
4419 searchy = params.set_default("searchy",3);
4420 searchz = params.set_default("searchz",3);
4421 }
4422
4423 Transform* initxform;
4424 if (params.has_key("initxform") ) {
4425 // Unlike the 2d refine aligner, this class doesn't require the starting transform's
4426 // parameters to form the starting guess. Instead the Transform itself
4427 // is perturbed carefully (using quaternion rotation) to overcome problems that arise
4428 // when you use orthogonally-based Euler angles
4429 initxform = params["initxform"];
4430 }else {
4431 initxform = new Transform(); // is the identity
4432 }
4433
4434 float lphi = params.set_default("phi0",0.0f);
4435 float uphi = params.set_default("phi1",360.0f);
4436 float dphi = params.set_default("dphi",10.f);
4437 float threshold = params.set_default("threshold",0.f);
4438 if (threshold < 0.0f) throw InvalidParameterException("The threshold parameter must be greater than or equal to zero");
4439 bool verbose = params.set_default("verbose",false);
4440
4441 //in case we are aligning tomos
4442 Dict altered_cmp_params(cmp_params);
4443 if (cmp_name == "ccc.tomo") {
4444 altered_cmp_params.set_default("searchx", searchx);
4445 altered_cmp_params.set_default("searchy", searchy);
4446 altered_cmp_params.set_default("searchz", searchz);
4447 altered_cmp_params.set_default("norm", true);
4448 }
4449
4450 vector<Dict> solns;
4451 if (nsoln == 0) return solns; // What was the user thinking?
4452 for (unsigned int i = 0; i < nsoln; ++i ) {
4453 Dict d;
4454 d["score"] = 1.e24;
4455 Transform t; // identity by default
4456 d["xform.align3d"] = &t; // deep copy is going on here
4457 solns.push_back(d);
4458 }
4459
4460 Dict d;
4461 d["inc_mirror"] = true; // This should probably always be true for 3D case. If it ever changes we might have to make inc_mirror a parameter
4462 if ( params.has_key("delta") && params.has_key("n") ) {
4463 throw InvalidParameterException("The delta and n parameters are mutually exclusive in the RT3DSphereAligner aligner");
4464 } else if ( params.has_key("n") ) {
4465 d["n"] = params["n"];
4466 } else {
4467 // If they didn't specify either then grab the default delta - if they did supply delta we're still safe doing this
4468 d["delta"] = params.set_default("delta",10.f);
4469 }
4470
4471 if ((string)params.set_default("orientgen","eman")=="eman") d["perturb"]=0;
4472 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
4473 vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
4474
4475 bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
4476
4477 //precompute fixed FT, saves a LOT of time!!!
4478 EMData * this_imgfft = 0;
4479 if(dotrans || tomography){
4480 this_imgfft = this_img->do_fft();
4481 }
4482
4483#ifdef EMAN2_USING_CUDA
4484 if(EMData::usecuda == 1) {
4485 cout << "Using CUDA for 3D alignment" << endl;
4486 if(!to->getcudarodata()) to->copy_to_cudaro(); // Safer call
4487 if(!this_img->getcudarwdata()) this_img->copy_to_cuda();
4488 if(this_imgfft) this_imgfft->copy_to_cuda();
4489 }
4490#endif
4491
4492 Transform trans = Transform();
4493 Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
4494
4495 bool use_cpu = true;
4496 for(vector<Transform>::const_iterator trans_it = transforms.begin(); trans_it != transforms.end(); trans_it++) {
4497 Dict params = trans_it->get_params("eman");
4498
4499 if (verbose) {
4500 float alt = params["alt"];
4501 float az = params["az"];
4502 cout << "Trying angle alt: " << alt << " az: " << az << endl;
4503 }
4504
4505 for( float phi = lphi; phi < uphi; phi += dphi ) {
4506 params["phi"] = phi;
4507 Transform t(params);
4508 t = t*(*initxform);
4509
4510 EMData* transformed;
4511 transformed = to->process("xform",Dict("transform",&t));
4512
4513 //need to do things a bit diffrent if we want to compare two tomos
4514 float best_score = 0.0f;
4515 // Dotrans is effectievly ignored for tomography
4516 if(dotrans || tomography){
4517 EMData* ccf = transformed->calc_ccf(this_imgfft);
4518#ifdef EMAN2_USING_CUDA
4519 if(EMData::usecuda == 1){
4520 // I use the following code rather than ccc.tomo to avoid doing two CCCs
4521 use_cpu = false;
4522 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
4523 trans.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
4524 t = trans*t; //composite transform to reflect the fact that we have done a rotation first and THEN a transformation
4525 if (tomography) {
4526 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
4527 best_score = -(data->peak - stats.x)/sqrt(stats.y); // Normalize, this is better than calling the norm processor since we only need to normalize one point
4528 } else {
4529 best_score = -data->peak;
4530 }
4531 delete data;
4532 }
4533#endif
4534 if(use_cpu){
4535 // I use the following code rather than ccc.tomo to avoid doing two CCCs
4536 if(tomography) ccf->process_inplace("normalize");
4537 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
4538 trans.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
4539 t = trans*t; //composite transform to reflect the fact that we have done a rotation first and THEN a transformation
4540 best_score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
4541 }
4542 delete ccf; ccf =0;
4543 delete transformed; transformed = 0;// this is to stop a mem leak
4544 }
4545
4546 if(!tomography){
4547 if(!transformed) transformed = to->process("xform",Dict("transform",&t));
4548 best_score = c->cmp(this_img,transformed);
4549 delete transformed; transformed = 0;
4550 }
4551
4552 unsigned int j = 0;
4553 //cout << "alt " <<float(t.get_rotation("eman").get("alt")) << " az " << float(t.get_rotation("eman").get("az")) << " phi " << float(t.get_rotation("eman").get("phi")) << endl;
4554 for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
4555 if ( (float)(*it)["score"] > best_score ) { // Note greater than - EMAN2 preferes minimums as a matter of policy
4556 vector<Dict>::reverse_iterator rit = solns.rbegin();
4557 copy(rit+1,solns.rend()-j,rit);
4558 Dict& d = (*it);
4559 d["score"] = best_score;
4560 t.invert(); //We actually moved the ref onto the moving, so we need to invert to do the opposite(this is done b/c the ref is aligned to the sym axis, whereas the mvoing is not)
4561 d["xform.align3d"] = &t; // deep copy is going on here
4562 break;
4563 }
4564 }
4565
4566 }
4567 }
4568
4569 if(this_imgfft) {delete this_imgfft; this_imgfft = 0;}
4570 if(sym!=0) delete sym;
4571 if (c != 0) delete c;
4572
4573 return solns;
4574
4575}
4576
4577//Could refactor the code here......(But not really woth it)
4578EMData* RT3DSymmetryAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
4579{
4580
4581 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
4582
4583 Transform* tr = (Transform*) alis[0]["xform.align3d"];
4584 EMData* soln = this_img->process("xform",Dict("transform",tr));
4585 soln->set_attr("xform.align3d",tr);
4586 delete tr; tr = 0;
4587
4588 return soln;
4589
4590}
4591
4592vector<Dict> RT3DSymmetryAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const
4593{
4594
4595 bool verbose = params.set_default("verbose",false);
4596 Transform* ixform;
4597 if (params.has_key("transform") ) {
4598 ixform = params["transform"];
4599 }else{
4600 ixform = new Transform(); // is the identity
4601 }
4602
4603 //Initialize a soln dict
4604 vector<Dict> solns;
4605 if (nsoln == 0) return solns; // What was the user thinking?
4606 for (unsigned int i = 0; i < nsoln; ++i ) {
4607 Dict d;
4608 d["score"] = 1.e24;
4609 Transform t; // identity by default
4610 d["xform.align3d"] = &t; // deep copy is going on here
4611 solns.push_back(d);
4612 }
4613
4614 #ifdef EMAN2_USING_CUDA
4615 if(EMData::usecuda == 1) {
4616 cout << "Using CUDA for 3D sym alignment" << endl;
4617 if(!this_img->getcudarwdata()) this_img->copy_to_cudaro();
4618 if(!to->getcudarwdata()) to->copy_to_cuda();
4619 }
4620 #endif
4621
4622 //Generate symmetry related orientations
4623 vector<Transform> syms = Symmetry3D::get_symmetries((string)params.set_default("sym","icos"));
4624 Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
4625
4626 float score = 0.0f;
4627 for ( vector<Transform>::const_iterator symit = syms.begin(); symit != syms.end(); ++symit ) {
4628 //Here move to sym position and compute the score
4629 Transform sympos = (*symit)*(*ixform);
4630 EMData* transformed = this_img->process("xform",Dict("transform", &sympos));
4631 score = c->cmp(transformed,to);
4632 delete transformed; transformed = 0;
4633
4634 if (verbose) {
4635 Dict rots = sympos.get_rotation("eman");
4636 cout <<"Score is: " << score << " az " << float(rots["az"]) << " alt " << float(rots["alt"]) << " phi " << float(rots["phi"]) << endl;
4637 }
4638
4639 unsigned int j = 0;
4640 for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
4641 if ( (float)(*it)["score"] > score ) { // Note greater than - EMAN2 preferes minimums as a matter of policy
4642 vector<Dict>::reverse_iterator rit = solns.rbegin();
4643 copy(rit+1,solns.rend()-j,rit);
4644 Dict& d = (*it);
4645 d["score"] = score;
4646 d["xform.align3d"] = &sympos; // deep copy is going on here
4647 break;
4648 }
4649 }
4650 }
4651
4652 if (c != 0) delete c;
4653
4654 return solns;
4655}
4656
4657namespace {
4658float frm_2d_Align(EMData *this_img, EMData *to, float *frm2dhhat, EMData* selfpcsfft, int p_max_input,int rsize, float &com_this_x, float &com_this_y, float &com_with_x, float &com_with_y,const string & cmp_name, const Dict& cmp_params)
4659{
4660 int size=rsize;
4661 float dx,dy;
4662 int bw=size/2;
4663 int MAXR=this_img->get_ysize()/2;
4664 //int MAXR=size;
4665 unsigned long tsize=2*size;
4666 unsigned long ind1=0, ind2=0, ind3=0, ind4=0, ind41=0;
4667 unsigned long index0=0;
4668 int i=0, j=0, m=0, n=0, r=0;
4669 int loop_rho=0, rho_best=0;
4670
4671 float* gnr2 = new float[size*2];
4672 float* maxcor = new float[size+1]; // MAXR need change
4673
4674 int p_max=p_max_input;
4675 float* result = new float[5*(p_max+1)];
4676 float* cr=new float[size*(bw+1)];
4677 float* ci=new float[size*(bw+1)];
4678 EMData *data_in=new EMData;
4679 data_in->set_complex(true);
4680 data_in->set_fftodd(false);
4681 data_in->set_ri(true);
4682 data_in->set_size(size+2,size,1);
4683 float *in=data_in->get_data();
4684
4685 float *self_sampl_fft = selfpcsfft->get_data(); // ming f(r)
4686
4687 float maxcor_sofar=0.0f;
4688 int p=0;
4689
4690 for(p=0; p<=p_max; ++p){
4691 ind1=p*size*bw;
4692 for (i=0;i<size;++i)
4693 for (j=0;j<bw+1;++j){
4694 cr[i*(bw+1)+j]=0.0;
4695 ci[i*(bw+1)+j]=0.0;
4696 }
4697 for(n=0;n<bw;++n){ // loop for n
4698 ind2=(ind1+n);
4699 index0=n*(bw+1);
4700 for(r=0;r<=MAXR;++r) {
4701 ind3=(ind2+r*bw)*size;
4702 for(m=0;m<size;m++){ // take back hat{h(n,r,p)}(m)
4703 ind4=(ind3+m)*2;
4704 ind41=ind4+1;
4705 gnr2[2*m]=frm2dhhat[ind4];
4706 gnr2[2*m+1]=frm2dhhat[ind41];
4707 }
4708 for(m=0;m<bw;++m){
4709 float tempr=self_sampl_fft[2*m+r*(size+2)]*r;
4710 float tempi=self_sampl_fft[2*m+1+r*(size+2)]*r;
4711 float gnr2_r=gnr2[2*m];
4712 float gnr2_i=gnr2[2*m+1];
4713 cr[n*(bw+1)+m]+=gnr2_r*tempr+gnr2_i*tempi;
4714 ci[n*(bw+1)+m]+=gnr2_i*tempr-gnr2_r*tempi;
4715 if(n!=0){ // m,-n
4716 if(m!= 0){
4717 int ssize=tsize-2*m; // ssize = 2*size-2m
4718 int ssize1=ssize+1;
4719 float gnr2_r=gnr2[ssize];
4720 float gnr2_i=gnr2[ssize1];
4721 cr[(size-n)*(bw+1)+m]+=gnr2_r*tempr-gnr2_i*tempi;
4722 ci[(size-n)*(bw+1)+m]-=gnr2_i*tempr+gnr2_r*tempi;
4723 }
4724 else{
4725 cr[(size-n)*(bw+1)+m]+=*(gnr2)*tempr-*(gnr2+1)*tempi;
4726 ci[(size-n)*(bw+1)+m]-=*(gnr2+1)*tempr+*(gnr2)*tempi;
4727 }
4728 }
4729 }
4730 }
4731 }
4732 for (int cii=0; cii<size*(bw+1);++cii){
4733 in[2*cii]=cr[cii];
4734 in[2*cii+1]=ci[cii];
4735 //printf("cii=%d,in[2i+1]=%f\n",cii, cr[cii]);
4736 }
4737
4738 EMData *data_out;
4739 data_out=data_in->do_ift();
4740 float *c=data_out->get_data();
4741 float tempr=0.0f, corre_fcs=999.0f;
4742
4743 int n_best=0, m_best=0;
4744 float temp=-100.0f;
4745 for(n=0;n<size;++n){// move Tri_2D to Tri = c(phi,phi';rho)
4746 for(m=0;m<size;++m){
4747 temp=c[n*size+m];
4748 if(temp>tempr) {
4749 tempr=temp;
4750 n_best=n;
4751 m_best=m;
4752 }
4753 }
4754 }
4755 delete data_out;
4756
4757 float corre,Phi2,Phi,Tx,Ty,Vx, Vy;
4758
4759 //for (n_best=0;n_best<bw;n_best++)
4760 // for (m_best=0;m_best<2*bw;m_best++){
4761 //n_best=0;
4762 //m_best=70;
4763 Phi2=n_best*M_PI/bw; // ming this is reference image rotation angle
4764 Phi=m_best*M_PI/bw; // ming this is particle image rotation angle
4765 Vx=p*cos(Phi);//should use the angle of the centered one
4766 Vy=-p*sin(Phi);
4767 Tx=Vx+(floor(com_this_x+0.5f)-floor(com_with_x+0.5f));
4768 Ty=Vy+(floor(com_this_y+0.5f)-floor(com_with_y+0.5f));
4769
4770 dx=-Tx; // the Rota & Trans value (Tx,Ty, ang_keep) are for the projection image,
4771 dy=-Ty; // need to convert to raw image
4772
4773 EMData *this_tmp=this_img->copy();//ming change to to
4774 this_tmp->rotate(-(Phi2-Phi)*180.0f,0.0f,0.0f);
4775 this_tmp->translate(dx,dy,0.0);
4776
4777 corre=this_tmp->cmp(cmp_name,to,cmp_params);
4778 //printf("corre=%f\n",corre);
4779 delete this_tmp;
4780 if(corre<=corre_fcs) { //ming, cmp use smaller value stands for more similarity
4781 corre_fcs=corre;
4782 result[0+5*p] = float(p); // rho
4783 result[1+5*p] = corre; // correlation_fcs
4784 result[2+5*p] = (Phi2-Phi)*180.0f; // rotation angle
4785 result[3+5*p] = Tx; // Translation_x
4786 result[4+5*p] = Ty; // Translation_y
4787 }
4788 maxcor[p]=corre_fcs; // maximum correlation for current rho
4789 if(corre_fcs<maxcor_sofar) {
4790 maxcor_sofar=corre_fcs; // max correlation up to current rho
4791 rho_best=p; // the rho value with maxinum correlation value
4792 }
4793 if(p>=4){
4794 if(maxcor[p] < maxcor[p-1] && maxcor[p-1] < maxcor[p-2]&& maxcor[p-2] < maxcor[p-3] && maxcor[p-3] < maxcor[p-4]){
4795 loop_rho=1;
4796 break; //exit p loop
4797 }
4798 }
4799 } // end for p
4800 //}//test my method
4801 if(loop_rho == 1)
4802 p=p+1;
4803 int rb5=5*rho_best;
4804 float fsc = result[1+rb5];
4805 float ang_keep = result[2+rb5];
4806 float Tx = result[3+rb5];
4807 float Ty = result[4+rb5];
4808 delete[] gnr2;
4809 delete[] maxcor;
4810 delete[] result;
4811 delete[] cr;
4812 cr=0;
4813 delete[] ci;
4814 ci=0;
4815 delete data_in; // ming add
4816 dx = -Tx; // the Rota & Trans value (Tx,Ty, ang_keep) are for the projection image,
4817 dy = -Ty; // need to convert to raw image
4818 this_img->rotate(-ang_keep,0,0); // ming change this to this_img??
4819 this_img->translate(dx,dy,0.0); // ming change this to this_img
4820
4821
4822 Transform tsoln(Dict("type","2d","alpha",ang_keep));
4823 tsoln.set_trans(dx,dy);
4824 this_img->set_attr("xform.align2d",&tsoln);
4825#ifdef DEBUG
4826 float fsc_best=this_img->cmp(cmp_name,to,cmp_params);
4827 printf("rho_best=%d fsc=%f fsc_best=%f dx=%f dy=%f ang_keep=%f com_withx=%f com_selfx=%f com_selfy=%f\n",rho_best,fsc,fsc_best,dx,dy,ang_keep,com_with_x,com_this_x,com_this_y);
4828#endif
4829 return fsc; // return the fsc coefficients
4830} // FRM2D aligner sub_class
4831} // end namespace
4832
4833
4834EMData *FRM2DAligner::align(EMData * this_img, EMData * to,
4835 const string & cmp_name, const Dict& cmp_params) const
4836{
4837 if (!this_img) {
4838 return 0;
4839 }
4840 if (to && !EMUtil::is_same_size(this_img, to))
4841 throw ImageDimensionException("Images must be the same size to perform translational alignment");
4842
4843 int nx=this_img->get_xsize();
4844 int ny=this_img->get_ysize();
4845 int size =(int)floor(M_PI*ny/4.0);
4846 size =Util::calc_best_fft_size(size);//ming bestfftsize(size);
4847 int MAXR=ny/2;
4848 //int MAXR=size;
4849 EMData *this_temp=this_img->copy(); // ming change avg to to
4850 FloatPoint com_test,com_test1;
4851 com_test=this_temp->calc_center_of_mass();//ming add
4852 float com_this_x=com_test[0];
4853 float com_this_y=com_test[1];
4854 delete this_temp;
4855
4856
4857 EMData *that_temp=to->copy();
4858 com_test1=that_temp->calc_center_of_mass();
4859 float com_with_x=com_test1[0];
4860 float com_with_y=com_test1[1];
4861 delete that_temp;
4862
4863 EMData *avg_frm=to->copy();
4864 float dx,dy;
4865 //float dx=-(com_with_x-nx/2); //ming
4866 //float dy=-(com_with_y-ny/2); //ming
4867 //avg_frm->translate(dx,dy,0.0);
4868 EMData *withpcs=avg_frm->unwrap_largerR(0,MAXR,size,float(MAXR)); // ming, something wrong inside this subroutine
4869 //EMData *withpcs=avg_frm->unwrap(-1,-1,-1,0,0,1);
4870 EMData *withpcsfft=withpcs->oneDfftPolar(size, float(MAXR), float(MAXR));
4871
4872 float *sampl_fft=withpcsfft->get_data(); //
4873 delete avg_frm;
4874 delete withpcs;
4875
4876 int bw=size/2;
4877 unsigned long ind1=0, ind2=0, ind3=0, ind4=0, ind41=0;
4878 float pi2=2.0*M_PI, r2;
4879
4880 EMData *data_in=new EMData;
4881 data_in->set_complex(true);
4882 data_in->set_ri(1);
4883 data_in->set_size(2*size,1,1);
4884 float * comp_in=data_in->get_data();
4885
4886 int p_max=3;
4887 float *frm2dhhat=0;
4888
4889 if( (frm2dhhat=(float *)malloc((p_max+1)*(size+2)*bw*size*2* sizeof(float)))==NULL){
4890 cout <<"Error in allocating memory 13. \n";
4891 exit(1);
4892 }
4893 //printf("p_max=%d\n",p_max);
4894 float *sb=0, *cb=0; // sin(beta) and cos(beta) for get h_hat, 300>size
4895 if((sb=new float[size])==NULL||(cb=new float[size])==NULL) {
4896 cout <<"can't allocate more memory, terminating. \n";
4897 exit(1);
4898 }
4899 for(int i=0;i<size;++i) { // beta sampling, to calculate beta' and r'
4900 float beta=i*M_PI/bw;
4901 sb[i]=sin(beta);
4902 cb[i]=cos(beta);
4903 }
4904
4905 for(int p=0; p<=p_max; ++p){
4906 ind1=p*size*bw;
4907 float pp2=(float)(p*p);
4908 for(int n=0;n<bw;++n){ /* loop for n */
4909 ind2=ind1+n;
4910 for(int r=0;r<=MAXR;++r) {
4911 ind3=(ind2+r*bw)*size;
4912 float rr2=(float)(r*r);
4913 float rp2=(float)(r*p);
4914 for(int i=0;i<size;++i){ // beta sampling, to get beta' and r'
4915 r2=std::sqrt((float)(rr2+pp2-2.0*rp2*cb[i])); // r2->r'
4916 int r1=(int)floor(r2+0.5f); // for computing gn(r')
4917 if(r1>MAXR){
4918 comp_in[2*i]=0.0f;
4919 comp_in[2*i+1]=0.0f;
4920 }
4921 else{
4922 float gn_r=sampl_fft[2*n+r1*(size+2)]; // real part of gn(r')
4923 float gn_i=sampl_fft[2*n+1+r1*(size+2)]; // imaginary part of gn(r')
4924 float sb2, cb2, cn, sn;
4925 if(n!=0){
4926 if(r2 != 0.0){
4927 sb2=r*sb[i]/r2;
4928 cb2=(r*cb[i]-p)/r2;
4929 }
4930 else{
4931 sb2=0.0;
4932 cb2=1.0;
4933 }
4934 if(sb2>1.0) sb2= 1.0f;
4935 if(sb2<-1.0)sb2=-1.0f;
4936 if(cb2>1.0) cb2= 1.0f;
4937 if(cb2<-1.0)cb2=-1.0f;
4938 float beta2=atan2(sb2,cb2);
4939 if(beta2<0.0) beta2+=pi2;
4940 float nb2=n*beta2;
4941 cn=cos(nb2);
4942 sn=sin(nb2);
4943 }
4944 else{
4945 cn=1.0f; sn=0.0f;
4946 }
4947 comp_in[2*i]=cn*gn_r-sn*gn_i;
4948 comp_in[2*i+1]=-(cn*gn_i+sn*gn_r);
4949 }
4950 }
4951 EMData *data_out;
4952 data_out=data_in->do_fft();
4953 float * comp_out=data_out->get_data();
4954 for(int m=0;m<size;m++){ // store hat{h(n,r,p)}(m)
4955 ind4=(ind3+m)*2;
4956 ind41=ind4+1;
4957 frm2dhhat[ind4]=comp_out[2*m];
4958 frm2dhhat[ind41]=comp_out[2*m+1];
4959 }
4960 delete data_out;
4961 }
4962 }
4963 }
4964
4965 delete[] sb;
4966 delete[] cb;
4967 delete data_in;
4968 delete withpcsfft;
4969
4970 float dot_frm0=0.0f, dot_frm1=0.0f;
4971 EMData *da_nFlip=0, *da_yFlip=0, *dr_frm=0;
4972 //dr_frm=this_img->copy();
4973 for (int iFlip=0;iFlip<=1;++iFlip){
4974 if (iFlip==0){dr_frm=this_img->copy(); da_nFlip=this_img->copy();}
4975 else {dr_frm=this_img->copy(); da_yFlip=this_img->copy();}
4976 if(iFlip==1) {com_this_x=nx-com_this_x; } //ming // image mirror about Y axis, so y keeps the same
4977
4978 dx=-(com_this_x-nx/2); //ming
4979 dy=-(com_this_y-ny/2); //ming
4980 dr_frm->translate(dx,dy,0.0); // this
4981 EMData *selfpcs = dr_frm->unwrap_largerR(0,MAXR,size, (float)MAXR);
4982 //EMData *selfpcs=dr_frm->unwrap(-1,-1,-1,0,0,1);
4983 EMData *selfpcsfft = selfpcs->oneDfftPolar(size, (float)MAXR, (float)MAXR);
4984 delete selfpcs;
4985 delete dr_frm;
4986 if(iFlip==0)
4987 dot_frm0=frm_2d_Align(da_nFlip,to, frm2dhhat, selfpcsfft, p_max, size, com_this_x, com_this_y, com_with_x, com_with_y,cmp_name,cmp_params);
4988 else
4989 dot_frm1=frm_2d_Align(da_yFlip,to, frm2dhhat, selfpcsfft, p_max, size, com_this_x, com_this_y, com_with_x, com_with_y,cmp_name,cmp_params);
4990 delete selfpcsfft;
4991 }
4992
4993 delete[] frm2dhhat;
4994 if(dot_frm0 <=dot_frm1) {
4995#ifdef DEBUG
4996 printf("best_corre=%f, no flip\n",dot_frm0);
4997#endif
4998 delete da_yFlip;
4999 return da_nFlip;
5000 }
5001 else {
5002#ifdef DEBUG
5003 printf("best_corre=%f, flipped\n",dot_frm1);
5004#endif
5005 delete da_nFlip;
5006 return da_yFlip;
5007 }
5008}
5009
5010#ifdef SPARX_USING_CUDA
5011
5013 image_stack = NULL;
5014 image_stack_filtered = NULL;
5015 ccf = NULL;
5016 if (id != -1) cudasetup(id);
5017}
5018
5019void CUDA_Aligner::finish() {
5020 if (image_stack) free(image_stack);
5022 if (ccf) free(ccf);
5023 image_stack = NULL;
5024 image_stack_filtered = NULL;
5025 ccf = NULL;
5026}
5027
5028void CUDA_Aligner::setup(int nima, int nx, int ny, int ring_length, int nring, int ou, float step, int kx, int ky, bool ctf) {
5029
5030 NIMA = nima;
5031 NX = nx;
5032 NY = ny;
5033 RING_LENGTH = ring_length;
5034 NRING = nring;
5035 STEP = step;
5036 KX = kx;
5037 KY = ky;
5038 OU = ou;
5039 CTF = ctf;
5040
5041 image_stack = (float *)malloc(NIMA*NX*NY*sizeof(float));
5042 if (CTF == 1) image_stack_filtered = (float *)malloc(NIMA*NX*NY*sizeof(float));
5043 ccf = (float *)malloc(2*(2*KX+1)*(2*KY+1)*NIMA*(RING_LENGTH+2)*sizeof(float));
5044}
5045
5046void CUDA_Aligner::insert_image(EMData *image, int num) {
5047
5048 int base_address = num*NX*NY;
5049
5050 for (int y=0; y<NY; y++)
5051 for (int x=0; x<NX; x++)
5052 image_stack[base_address+y*NX+x] = (*image)(x, y);
5053}
5054
5055void CUDA_Aligner::filter_stack(vector<float> ctf_params) {
5056
5057 float *params;
5058
5059 params = (float *)malloc(NIMA*6*sizeof(float));
5060
5061 for (int i=0; i<NIMA*6; i++) params[i] = ctf_params[i];
5062
5063 filter_image(image_stack, image_stack_filtered, NIMA, NX, NY, params);
5064
5065 free(params);
5066}
5067
5068void CUDA_Aligner::sum_oe(vector<float> ctf_params, vector<float> ali_params, EMData *ave1, EMData *ave2) {
5069
5070 float *ctf_p, *ali_p, *av1, *av2;
5071
5072 ctf_p = (float *)malloc(NIMA*6*sizeof(float));
5073 ali_p = (float *)malloc(NIMA*4*sizeof(float));
5074
5075 if (CTF == 1) {
5076 for (int i=0; i<NIMA*6; i++) ctf_p[i] = ctf_params[i];
5077 }
5078 for (int i=0; i<NIMA*4; i++) ali_p[i] = ali_params[i];
5079
5080 av1 = ave1->get_data();
5081 av2 = ave2->get_data();
5082
5083 rot_filt_sum(image_stack, NIMA, NX, NY, CTF, ctf_p, ali_p, av1, av2);
5084
5085 free(ctf_p);
5086 free(ali_p);
5087}
5088
5089vector<float> CUDA_Aligner::alignment_2d(EMData *ref_image_em, vector<float> sx_list, vector<float> sy_list, int silent) {
5090
5091 float *ref_image, max_ccf;
5092 int base_address, ccf_offset;
5093 float ts, tm;
5094 float ang, sx = 0, sy = 0, mirror, co, so, sxs, sys;
5095 float *sx2, *sy2;
5096 vector<float> align_result;
5097
5098 sx2 = (float *)malloc(NIMA*sizeof(float));
5099 sy2 = (float *)malloc(NIMA*sizeof(float));
5100
5101 ref_image = ref_image_em->get_data();
5102
5103 for (int i=0; i<NIMA; i++) {
5104 sx2[i] = sx_list[i];
5105 sy2[i] = sy_list[i];
5106 }
5107
5108 if (CTF == 1) {
5109 calculate_ccf(image_stack_filtered, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, silent);
5110 } else {
5111 calculate_ccf(image_stack, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, silent);
5112 }
5113
5114 ccf_offset = NIMA*(RING_LENGTH+2)*(2*KX+1)*(2*KY+1);
5115
5116 for (int im=0; im<NIMA; im++) {
5117 max_ccf = -1.0e22;
5118 for (int kx=-KX; kx<=KX; kx++) {
5119 for (int ky=-KY; ky<=KY; ky++) {
5120 base_address = (((ky+KY)*(2*KX+1)+(kx+KX))*NIMA+im)*(RING_LENGTH+2);
5121 for (int l=0; l<RING_LENGTH; l++) {
5122 ts = ccf[base_address+l];
5123 tm = ccf[base_address+l+ccf_offset];
5124 if (ts > max_ccf) {
5125 ang = float(l)/RING_LENGTH*360.0;
5126 sx = -kx*STEP;
5127 sy = -ky*STEP;
5128 mirror = 0;
5129 max_ccf = ts;
5130 }
5131 if (tm > max_ccf) {
5132 ang = float(l)/RING_LENGTH*360.0;
5133 sx = -kx*STEP;
5134 sy = -ky*STEP;
5135 mirror = 1;
5136 max_ccf = tm;
5137 }
5138 }
5139 }
5140 }
5141 co = cos(ang*M_PI/180.0);
5142 so = -sin(ang*M_PI/180.0);
5143 sxs = sx*co - sy*so;
5144 sys = sx*so + sy*co;
5145
5146 align_result.push_back(ang);
5147 align_result.push_back(sxs);
5148 align_result.push_back(sys);
5149 align_result.push_back(mirror);
5150 }
5151
5152 free(sx2);
5153 free(sy2);
5154
5155 return align_result;
5156}
5157
5158
5159vector<float> CUDA_Aligner::ali2d_single_iter(EMData *ref_image_em, vector<float> ali_params, float csx, float csy, int silent, float delta) {
5160
5161 float *ref_image, max_ccf;
5162 int base_address, ccf_offset;
5163 float ts, tm;
5164 float ang = 0.0, sx = 0.0, sy = 0.0, co, so, sxs, sys;
5165 int mirror;
5166 float *sx2, *sy2;
5167 vector<float> align_result;
5168
5169 sx2 = (float *)malloc(NIMA*sizeof(float));
5170 sy2 = (float *)malloc(NIMA*sizeof(float));
5171
5172 ref_image = ref_image_em->get_data();
5173
5174 for (int i=0; i<NIMA; i++) {
5175 ang = ali_params[i*4]/180.0*M_PI;
5176 sx = (ali_params[i*4+3] < 0.5)?(ali_params[i*4+1]-csx):(ali_params[i*4+1]+csx);
5177 sy = ali_params[i*4+2]-csy;
5178 co = cos(ang);
5179 so = sin(ang);
5180 sx2[i] = -(sx*co-sy*so);
5181 sy2[i] = -(sx*so+sy*co);
5182 }
5183
5184 if (CTF == 1) {
5185 calculate_ccf(image_stack_filtered, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, silent);
5186 } else {
5187 calculate_ccf(image_stack, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, silent);
5188 }
5189
5190 ccf_offset = NIMA*(RING_LENGTH+2)*(2*KX+1)*(2*KY+1);
5191
5192 float sx_sum = 0.0f;
5193 float sy_sum = 0.0f;
5194
5195 int dl;
5196 dl = static_cast<int>(delta/360.0*RING_LENGTH);
5197 if (dl<1) { dl = 1; }
5198
5199 for (int im=0; im<NIMA; im++) {
5200 max_ccf = -1.0e22;
5201 for (int kx=-KX; kx<=KX; kx++) {
5202 for (int ky=-KY; ky<=KY; ky++) {
5203 base_address = (((ky+KY)*(2*KX+1)+(kx+KX))*NIMA+im)*(RING_LENGTH+2);
5204 for (int l=0; l<RING_LENGTH; l+=dl) {
5205 ts = ccf[base_address+l];
5206</