00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include "emfft.h"
00036 #include "cmp.h"
00037 #include "aligner.h"
00038 #include "averager.h"
00039 #include "emdata.h"
00040 #include "processor.h"
00041 #include "util.h"
00042 #include "symmetry.h"
00043 #include <gsl/gsl_multimin.h>
00044 #include "plugins/aligner_template.h"
00045
00046 #ifdef EMAN2_USING_CUDA
00047 #include "cuda/cuda_processor.h"
00048 #include "cuda/cuda_cmp.h"
00049 #endif
00050
00051 #ifdef SPARX_USING_CUDA
00052 #include <sparx/cuda/cuda_ccf.h>
00053 #endif
00054
00055 #define EMAN2_ALIGNER_DEBUG 0
00056
00057 using namespace EMAN;
00058
00059 const string TranslationalAligner::NAME = "translational";
00060 const string RotationalAligner::NAME = "rotational";
00061 const string RotationalAlignerIterative::NAME = "rotational_iterative";
00062 const string RotatePrecenterAligner::NAME = "rotate_precenter";
00063 const string RotateTranslateAligner::NAME = "rotate_translate";
00064 const string RotateTranslateScaleAligner::NAME = "rotate_translate_scale";
00065 const string RotateTranslateAlignerIterative::NAME = "rotate_translate_iterative";
00066 const string RotateTranslateScaleAlignerIterative::NAME = "rotate_trans_scale_iter";
00067 const string RotateTranslateAlignerPawel::NAME = "rotate_translate_resample";
00068 const string RotateTranslateBestAligner::NAME = "rotate_translate_best";
00069 const string RotateFlipAligner::NAME = "rotate_flip";
00070 const string RotateFlipAlignerIterative::NAME = "rotate_flip_iterative";
00071 const string RotateTranslateFlipAligner::NAME = "rotate_translate_flip";
00072 const string RotateTranslateFlipScaleAligner::NAME = "rotate_trans_flip_scale";
00073 const string RotateTranslateFlipAlignerIterative::NAME = "rotate_translate_flip_iterative";
00074 const string RotateTranslateFlipScaleAlignerIterative::NAME = "rotate_trans_flip_scale_iter";
00075 const string RotateTranslateFlipAlignerPawel::NAME = "rotate_translate_flip_resample";
00076 const string RTFExhaustiveAligner::NAME = "rtf_exhaustive";
00077 const string RTFSlowExhaustiveAligner::NAME = "rtf_slow_exhaustive";
00078 const string RefineAligner::NAME = "refine";
00079 const string RefineAlignerCG::NAME = "refinecg";
00080 const string SymAlignProcessorQuat::NAME = "symalignquat";
00081 const string SymAlignProcessor::NAME = "symalign";
00082 const string Refine3DAlignerGrid::NAME = "refine_3d_grid";
00083 const string Refine3DAlignerQuaternion::NAME = "refine_3d";
00084 const string RT3DGridAligner::NAME = "rotate_translate_3d_grid";
00085 const string RT3DSphereAligner::NAME = "rotate_translate_3d";
00086 const string RT3DSymmetryAligner::NAME = "rotate_symmetry_3d";
00087 const string FRM2DAligner::NAME = "frm2d";
00088 const string ScaleAligner::NAME = "scale";
00089
00090
00091 template <> Factory < Aligner >::Factory()
00092 {
00093 force_add<TranslationalAligner>();
00094 force_add<RotationalAligner>();
00095 force_add<RotationalAlignerIterative>();
00096 force_add<RotatePrecenterAligner>();
00097 force_add<RotateTranslateAligner>();
00098 force_add<RotateTranslateScaleAligner>();
00099 force_add<RotateTranslateAlignerIterative>();
00100 force_add<RotateTranslateScaleAlignerIterative>();
00101 force_add<RotateTranslateAlignerPawel>();
00102 force_add<RotateFlipAligner>();
00103 force_add<RotateFlipAlignerIterative>();
00104 force_add<RotateTranslateFlipAligner>();
00105 force_add<RotateTranslateFlipScaleAligner>();
00106 force_add<RotateTranslateFlipAlignerIterative>();
00107 force_add<RotateTranslateFlipScaleAlignerIterative>();
00108 force_add<RotateTranslateFlipAlignerPawel>();
00109 force_add<RTFExhaustiveAligner>();
00110 force_add<RTFSlowExhaustiveAligner>();
00111 force_add<SymAlignProcessor>();
00112 force_add<RefineAligner>();
00113 force_add<RefineAlignerCG>();
00114 force_add<SymAlignProcessorQuat>();
00115 force_add<Refine3DAlignerGrid>();
00116 force_add<Refine3DAlignerQuaternion>();
00117 force_add<RT3DGridAligner>();
00118 force_add<RT3DSphereAligner>();
00119 force_add<RT3DSymmetryAligner>();
00120 force_add<FRM2DAligner>();
00121 force_add<ScaleAligner>();
00122
00123 }
00124
00125 vector<Dict> Aligner::xform_align_nbest(EMData *, EMData *, const unsigned int, const string &, const Dict&) const
00126 {
00127 vector<Dict> solns;
00128 return solns;
00129 }
00130
00131 EMData* ScaleAlignerABS::align_using_base(EMData * this_img, EMData * to,
00132 const string & cmp_name, const Dict& cmp_params) const
00133 {
00134
00135 float min = params.set_default("min",0.95f);
00136 float max = params.set_default("max",1.05f);
00137 float step = params.set_default("step",0.01f);
00138
00139
00140 Transform t = Transform();
00141 t.set_scale(max);
00142
00143
00144 float* oridata = this_img->get_data();
00145
00146
00147 Processor* proc = Factory <Processor>::get("xform", Dict());
00148 TransformProcessor* xform = dynamic_cast<TransformProcessor*>(proc);
00149
00150
00151 EMData* result = 0;
00152
00153 float bestscore = 1.0e37;
00154
00155 for(float i = max; i > min; i-=step){
00156
00157
00158 float* des_data = xform->transform(this_img,t);
00159 this_img->set_data(des_data);
00160 this_img->update();
00161
00162
00163 EMData* aligned = this_img->align(basealigner, to, basealigner_params, cmp_name, cmp_params);
00164 float score = aligned->cmp(cmp_name, to, cmp_params);
00165 if(score < bestscore){
00166 bestscore = score;
00167
00168 if(result != 0) delete result;
00169 result = aligned;
00170 result->set_attr("scalefactor",i);
00171 }else{
00172 delete aligned;
00173 }
00174
00175 delete des_data;
00176
00177 t.set_scale(i);
00178
00179
00180 this_img->set_data(oridata);
00181 }
00182
00183 if (!result) throw UnexpectedBehaviorException("Alignment score is infinity! Something is seriously wrong with the data!");
00184 if (proc != 0) delete proc;
00185
00186 return result;
00187
00188 };
00189
00190 EMData* ScaleAligner::align(EMData * this_img, EMData *to,
00191 const string& cmp_name, const Dict& cmp_params) const
00192 {
00193
00194
00195 float min = params.set_default("min",0.95f);
00196 float max = params.set_default("max",1.05f);
00197 float step = params.set_default("step",0.01f);
00198
00199 Transform t = Transform();
00200 t.set_scale(max);
00201 float* oridata = this_img->get_data();
00202
00203
00204 Processor* proc = Factory <Processor>::get("xform", Dict());
00205 TransformProcessor* xform = dynamic_cast<TransformProcessor*>(proc);
00206
00207
00208 float bestscale = 1.0;
00209 float bestscore = 1.0e37;
00210
00211 for(float i = max; i > min; i-=step){
00212
00213 float* des_data = xform->transform(this_img,t);
00214 this_img->set_data(des_data);
00215 this_img->update();
00216
00217
00218 float score = this_img->cmp(cmp_name, to, cmp_params);
00219 if(score < bestscore){
00220 bestscore = score;
00221 bestscale = i;
00222 }
00223
00224 delete des_data;
00225
00226 t.set_scale(i);
00227
00228
00229 this_img->set_data(oridata);
00230 }
00231
00232
00233
00234
00235 t.set_scale(bestscale);
00236 EMData* result = this_img->process("xform",Dict("transform",&t));
00237 result->set_attr("scalefactor",bestscale);
00238 if (proc != 0) delete proc;
00239
00240 return result;
00241
00242 }
00243
00244
00245
00246
00247
00248 EMData *TranslationalAligner::align(EMData * this_img, EMData *to,
00249 const string&, const Dict&) const
00250 {
00251 if (!this_img) {
00252 return 0;
00253 }
00254
00255 if (to && !EMUtil::is_same_size(this_img, to))
00256 throw ImageDimensionException("Images must be the same size to perform translational alignment");
00257
00258 EMData *cf = 0;
00259 int nx = this_img->get_xsize();
00260 int ny = this_img->get_ysize();
00261 int nz = this_img->get_zsize();
00262
00263 int masked = params.set_default("masked",0);
00264 int useflcf = params.set_default("useflcf",0);
00265 bool use_cpu = true;
00266
00267 #ifdef EMAN2_USING_CUDA
00268 if(EMData::usecuda == 1) {
00269
00270
00271
00272
00273
00274
00275
00276 }
00277 #endif // EMAN2_USING_CUDA
00278
00279 if (use_cpu) {
00280 if (useflcf) cf = this_img->calc_flcf(to);
00281 else cf = this_img->calc_ccf(to);
00282 }
00283
00284
00285 if (masked) {
00286 EMData *msk=this_img->process("threshold.notzero");
00287 EMData *sqr=to->process("math.squared");
00288 EMData *cfn=msk->calc_ccf(sqr);
00289 cfn->process_inplace("math.sqrt");
00290 float *d1=cf->get_data();
00291 float *d2=cfn->get_data();
00292 for (size_t i=0; i<(size_t)nx*ny*nz; ++i) {
00293 if (d2[i]!=0) d1[i]/=d2[i];
00294 }
00295 cf->update();
00296 delete msk;
00297 delete sqr;
00298 delete cfn;
00299 }
00300
00301 int maxshiftx = params.set_default("maxshift",-1);
00302 int maxshifty = params["maxshift"];
00303 int maxshiftz = params["maxshift"];
00304 int nozero = params["nozero"];
00305
00306 if (maxshiftx <= 0) {
00307 maxshiftx = nx / 4;
00308 maxshifty = ny / 4;
00309 maxshiftz = nz / 4;
00310 }
00311
00312 if (maxshiftx > nx / 2 - 1) maxshiftx = nx / 2 - 1;
00313 if (maxshifty > ny / 2 - 1) maxshifty = ny / 2 - 1;
00314 if (maxshiftz > nz / 2 - 1) maxshiftz = nz / 2 - 1;
00315
00316 if (nx == 1) maxshiftx = 0;
00317 if (ny == 1) maxshifty = 0;
00318 if (nz == 1) maxshiftz = 0;
00319
00320
00321 if (nozero) {
00322 cf->zero_corner_circulant(1);
00323 }
00324
00325 IntPoint peak;
00326 #ifdef EMAN2_USING_CUDA
00327 if (!use_cpu) {
00328 cout << "USe CUDA TA 2" << endl;
00329 if (nozero) throw UnexpectedBehaviorException("Nozero is not yet supported in CUDA");
00330 CudaPeakInfo* data = calc_max_location_wrap_cuda(cf->getcudarwdata(), cf->get_xsize(), cf->get_ysize(), cf->get_zsize(), maxshiftx, maxshifty, maxshiftz);
00331 peak = IntPoint(data->px,data->py,data->pz);
00332 free(data);
00333 }
00334 #endif // EMAN2_USING_CUDA
00335
00336 if (use_cpu) {
00337 peak = cf->calc_max_location_wrap(maxshiftx, maxshifty, maxshiftz);
00338 }
00339
00340 Vec3f cur_trans = Vec3f ( (float)-peak[0], (float)-peak[1], (float)-peak[2]);
00341
00342
00343 if (!to) {
00344 cur_trans /= 2.0f;
00345 int intonly = params.set_default("intonly",false);
00346 if (intonly) {
00347 cur_trans[0] = floor(cur_trans[0] + 0.5f);
00348 cur_trans[1] = floor(cur_trans[1] + 0.5f);
00349 cur_trans[2] = floor(cur_trans[2] + 0.5f);
00350 }
00351 }
00352
00353 if( cf ){
00354 delete cf;
00355 cf = 0;
00356 }
00357
00358 Dict params("trans",static_cast< vector<int> >(cur_trans));
00359 if (use_cpu){
00360 cf=this_img->process("xform.translate.int",params);
00361 }
00362 Transform t;
00363 t.set_trans(cur_trans);
00364
00365 #ifdef EMAN2_USING_CUDA
00366 if (!use_cpu) {
00367 cout << "USe CUDA TA 3" << endl;
00368
00369 cf = this_img->process("xform",Dict("transform",&t));
00370 }
00371 #endif // EMAN2_USING_CUDA
00372
00373 if ( nz != 1 ) {
00374
00375
00376 cf->set_attr("xform.align3d",&t);
00377 } else if ( ny != 1 ) {
00378
00379 cur_trans[2] = 0;
00380 t.set_trans(cur_trans);
00381 cf->set_attr("xform.align2d",&t);
00382 }
00383 return cf;
00384 }
00385
00386 EMData * RotationalAligner::align_180_ambiguous(EMData * this_img, EMData * to, int rfp_mode,int zscore) {
00387
00388
00389 EMData* this_img_rfp, * to_rfp;
00390 if (rfp_mode == 0) {
00391 this_img_rfp = this_img->make_rotational_footprint_e1();
00392 to_rfp = to->make_rotational_footprint_e1();
00393 } else if (rfp_mode == 1) {
00394 this_img_rfp = this_img->make_rotational_footprint();
00395 to_rfp = to->make_rotational_footprint();
00396 } else if (rfp_mode == 2) {
00397 this_img_rfp = this_img->make_rotational_footprint_cmc();
00398 to_rfp = to->make_rotational_footprint_cmc();
00399 } else {
00400 throw InvalidParameterException("rfp_mode must be 0,1 or 2");
00401 }
00402 int this_img_rfp_nx = this_img_rfp->get_xsize();
00403
00404
00405 EMData *cf = this_img_rfp->calc_ccfx(to_rfp, 0, this_img->get_ysize(),false,false,zscore);
00406
00407
00408
00409
00410
00411
00412
00413
00414 delete this_img_rfp; this_img_rfp = 0;
00415 delete to_rfp; to_rfp = 0;
00416
00417
00418 float *data = cf->get_data();
00419
00420 float peak = 0;
00421 int peak_index = 0;
00422 Util::find_max(data, this_img_rfp_nx, &peak, &peak_index);
00423
00424 if( cf ) {
00425 delete cf;
00426 cf = 0;
00427 }
00428 float rot_angle = (float) (peak_index * 180.0f / this_img_rfp_nx);
00429
00430
00431 Transform tmp(Dict("type","2d","alpha",rot_angle));
00432 cf=this_img->process("xform",Dict("transform",(Transform*)&tmp));
00433
00434
00435
00436 cf->set_attr("xform.align2d",&tmp);
00437 return cf;
00438 }
00439
00440 EMData *RotationalAligner::align(EMData * this_img, EMData *to,
00441 const string& cmp_name, const Dict& cmp_params) const
00442 {
00443 if (!to) throw InvalidParameterException("Can not rotational align - the image to align to is NULL");
00444
00445 #ifdef EMAN2_USING_CUDA
00446 if(EMData::usecuda == 1) {
00447
00448
00449 }
00450 #endif
00451
00452
00453 int rfp_mode = params.set_default("rfp_mode",2);
00454 int zscore = params.set_default("zscore",0);
00455 int ambig180 = params.set_default("ambig180",0);
00456
00457 EMData* rot_aligned = RotationalAligner::align_180_ambiguous(this_img,to,rfp_mode,zscore);
00458 Transform * tmp = rot_aligned->get_attr("xform.align2d");
00459 Dict rot = tmp->get_rotation("2d");
00460 float rotate_angle_solution = rot["alpha"];
00461 delete tmp;
00462
00463
00464 if (ambig180) {
00465 return rot_aligned;
00466 }
00467
00468 EMData *rot_align_180 = rot_aligned->process("math.rotate.180");
00469
00470
00471 float rot_cmp = rot_aligned->cmp(cmp_name, to, cmp_params);
00472 float rot_180_cmp = rot_align_180->cmp(cmp_name, to, cmp_params);
00473
00474
00475 float score = 0.0;
00476 EMData* result = NULL;
00477 if (rot_cmp < rot_180_cmp){
00478 result = rot_aligned;
00479 score = rot_cmp;
00480 delete rot_align_180; rot_align_180 = 0;
00481 } else {
00482 result = rot_align_180;
00483 score = rot_180_cmp;
00484 delete rot_aligned; rot_aligned = 0;
00485 rotate_angle_solution = rotate_angle_solution-180.0f;
00486 }
00487
00488
00489
00490 Transform tmp2(Dict("type","2d","alpha",rotate_angle_solution));
00491 result->set_attr("xform.align2d",&tmp2);
00492 return result;
00493 }
00494
00495
00496 EMData *RotatePrecenterAligner::align(EMData * this_img, EMData *to,
00497 const string&, const Dict&) const
00498 {
00499 if (!to) {
00500 return 0;
00501 }
00502
00503 int ny = this_img->get_ysize();
00504 int size = Util::calc_best_fft_size((int) (M_PI * ny * 1.5));
00505 EMData *e1 = this_img->unwrap(4, ny * 7 / 16, size, 0, 0, 1);
00506 EMData *e2 = to->unwrap(4, ny * 7 / 16, size, 0, 0, 1);
00507 EMData *cf = e1->calc_ccfx(e2, 0, ny);
00508
00509 float *data = cf->get_data();
00510
00511 float peak = 0;
00512 int peak_index = 0;
00513 Util::find_max(data, size, &peak, &peak_index);
00514 float a = (float) ((1.0f - 1.0f * peak_index / size) * 180. * 2);
00515
00516 Transform rot;
00517 rot.set_rotation(Dict("type","2d","alpha",(float)(a*180./M_PI)));
00518 EMData* rslt = this_img->process("xform",Dict("transform",&rot));
00519 rslt->set_attr("xform.align2d",&rot);
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530 if( e1 )
00531 {
00532 delete e1;
00533 e1 = 0;
00534 }
00535
00536 if( e2 )
00537 {
00538 delete e2;
00539 e2 = 0;
00540 }
00541
00542 if (cf) {
00543 delete cf;
00544 cf = 0;
00545 }
00546 return rslt;
00547 }
00548
00549 EMData *RotationalAlignerIterative::align(EMData * this_img, EMData *to,
00550 const string &, const Dict&) const
00551 {
00552 int r1 = params.set_default("r1",-1);
00553 int r2 = params.set_default("r2",-1);
00554
00555 EMData * to_polar = to->unwrap(r1,r2,-1,0,0,true);
00556 EMData * this_img_polar = this_img->unwrap(r1,r2,-1,0,0,true);
00557 int this_img_polar_nx = this_img_polar->get_xsize();
00558
00559 EMData * cf = this_img_polar->calc_ccfx(to_polar, 0, this_img->get_ysize());
00560
00561
00562 delete to_polar; to_polar = 0;
00563 delete this_img_polar; this_img_polar = 0;
00564
00565 float * data = cf->get_data();
00566 float peak = 0;
00567 int peak_index = 0;
00568 Util::find_max(data, this_img_polar_nx, &peak, &peak_index);
00569
00570 delete cf; cf = 0;
00571 float rot_angle = (float) (peak_index * 360.0f / this_img_polar_nx);
00572
00573
00574
00575 Transform tmp(Dict("type","2d","alpha",rot_angle));
00576 EMData * rotimg=this_img->process("xform",Dict("transform",(Transform*)&tmp));
00577 rotimg->set_attr("xform.align2d",&tmp);
00578
00579 return rotimg;
00580
00581 }
00582
00583 EMData *RotateTranslateAlignerIterative::align(EMData * this_img, EMData *to,
00584 const string & cmp_name, const Dict& cmp_params) const
00585 {
00586 int maxiter = params.set_default("maxiter", 3);
00587
00588 Dict trans_params;
00589 trans_params["intonly"] = 0;
00590 trans_params["maxshift"] = params.set_default("maxshift", -1);
00591 trans_params["useflcf"] = params.set_default("useflcf",0);
00592 trans_params["nozero"] = params.set_default("nozero",false);
00593
00594 Dict rot_params;
00595 rot_params["r1"] = params.set_default("r1", -1);
00596 rot_params["r2"] = params.set_default("r2", -1);
00597
00598 Transform t;
00599 EMData * moving_img = this_img;
00600 for(int it = 0; it < maxiter; it++){
00601
00602
00603 EMData * trans_align = moving_img->align("translational", to, trans_params, cmp_name, cmp_params);
00604 Transform * tt = trans_align->get_attr("xform.align2d");
00605 t = *tt*t;
00606 delete tt;
00607
00608
00609 EMData * rottrans_align = trans_align->align("rotational_iterative", to, rot_params, cmp_name, cmp_params);
00610 Transform * rt = rottrans_align->get_attr("xform.align2d");
00611 t = *rt*t;
00612 delete trans_align; trans_align = 0;
00613 delete rottrans_align; rottrans_align = 0;
00614 delete rt;
00615
00616
00617 if(it > 0){delete moving_img;}
00618
00619 moving_img = this_img->process("xform",Dict("transform",&t));
00620 }
00621
00622
00623 moving_img->set_attr("xform.align2d", &t);
00624
00625 return moving_img;
00626 }
00627
00628 EMData *RotateTranslateScaleAlignerIterative::align(EMData * this_img, EMData *to,
00629 const string & cmp_name, const Dict& cmp_params) const
00630 {
00631
00632
00633 basealigner_params["maxshift"] = params.set_default("maxshift", -1);
00634 basealigner_params["r1"] = params.set_default("r1",-1);
00635 basealigner_params["r2"] = params.set_default("r2",-1);
00636 basealigner_params["maxiter"] = params.set_default("maxiter",3);
00637 basealigner_params["nozero"] = params.set_default("nozero",false);
00638 basealigner_params["useflcf"] = params.set_default("useflcf",0);
00639
00640
00641 return align_using_base(this_img, to, cmp_name, cmp_params);
00642
00643 }
00644
00645 EMData *RotateTranslateAlignerPawel::align(EMData * this_img, EMData *to,
00646 const string & cmp_name, const Dict& cmp_params) const
00647 {
00648 if (cmp_name != "dot" && cmp_name != "ccc") throw InvalidParameterException("Resample aligner only works for dot and ccc");
00649
00650 int maxtx = params.set_default("tx", 0);
00651 int maxty = params.set_default("ty", 0);
00652 int r1 = params.set_default("r1",-1);
00653 int r2 = params.set_default("r2",-1);
00654
00655 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");
00656 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");
00657
00658
00659 float best_peak = -1.0e37;
00660 int best_peak_index = 0;
00661 int best_tx = 0;
00662 int best_ty = 0;
00663 int polarxsize = 0;
00664
00665 for(int x = -maxtx; x <= maxtx; x++){
00666 for(int y = -maxty; y <= maxty; y++){
00667
00668 EMData * to_polar = to->unwrap(r1,r2,-1,0,0,true);
00669 EMData * this_img_polar = this_img->unwrap(r1,r2,-1,x,y,true);
00670 EMData * cf = this_img_polar->calc_ccfx(to_polar, 0, this_img_polar->get_ysize());
00671
00672 polarxsize = this_img_polar->get_xsize();
00673
00674
00675 delete to_polar; to_polar = 0;
00676 delete this_img_polar; this_img_polar = 0;
00677
00678 float *data = cf->get_data();
00679 float peak = 0;
00680 int peak_index = 0;
00681 Util::find_max(data, polarxsize, &peak, &peak_index);
00682 delete cf; cf = 0;
00683
00684 if(peak > best_peak) {
00685 best_peak = peak;
00686 best_peak_index = peak_index;
00687 best_tx = x;
00688 best_ty = y;
00689 }
00690 }
00691 }
00692
00693 float rot_angle = (float) (best_peak_index * 360.0f / polarxsize);
00694
00695
00696 Transform tmptt(Dict("type","2d","alpha",0,"tx",-best_tx,"ty",-best_ty));
00697 Transform tmprot(Dict("type","2d","alpha",rot_angle,"tx",0,"ty",0));
00698 Transform total = tmprot*tmptt;
00699 EMData* rotimg=this_img->process("xform",Dict("transform",(Transform*)&total));
00700 rotimg->set_attr("xform.align2d",&total);
00701
00702 return rotimg;
00703
00704 }
00705
00706 EMData *RotateTranslateAligner::align(EMData * this_img, EMData *to,
00707 const string & cmp_name, const Dict& cmp_params) const
00708 {
00709
00710 #ifdef EMAN2_USING_CUDA
00711 if(EMData::usecuda == 1) {
00712
00713
00714 }
00715 #endif
00716
00717
00718 int zscore = params.set_default("zscore",0);
00719 int rfp_mode = params.set_default("rfp_mode",2);
00720 EMData *rot_align = RotationalAligner::align_180_ambiguous(this_img,to,rfp_mode,zscore);
00721 Transform * tmp = rot_align->get_attr("xform.align2d");
00722 Dict rot = tmp->get_rotation("2d");
00723 float rotate_angle_solution = rot["alpha"];
00724 delete tmp;
00725
00726 EMData *rot_align_180 = rot_align->process("math.rotate.180");
00727
00728 Dict trans_params;
00729 trans_params["intonly"] = 0;
00730 trans_params["maxshift"] = params.set_default("maxshift", -1);
00731 trans_params["useflcf"]=params.set_default("useflcf",0);
00732
00733
00734 trans_params["nozero"] = params.set_default("nozero",false);
00735 EMData* rot_trans = rot_align->align("translational", to, trans_params, cmp_name, cmp_params);
00736 if( rot_align ) {
00737 delete rot_align;
00738 rot_align = 0;
00739 }
00740
00741
00742 EMData* rot_180_trans = rot_align_180->align("translational", to, trans_params, cmp_name, cmp_params);
00743 if( rot_align_180 ) {
00744 delete rot_align_180;
00745 rot_align_180 = 0;
00746 }
00747
00748
00749 float cmp1 = rot_trans->cmp(cmp_name, to, cmp_params);
00750 float cmp2 = rot_180_trans->cmp(cmp_name, to, cmp_params);
00751
00752 EMData *result = 0;
00753 if (cmp1 < cmp2) {
00754 if( rot_180_trans ) {
00755 delete rot_180_trans;
00756 rot_180_trans = 0;
00757 }
00758 result = rot_trans;
00759 }
00760 else {
00761 if( rot_trans ) {
00762 delete rot_trans;
00763 rot_trans = 0;
00764 }
00765 result = rot_180_trans;
00766 rotate_angle_solution -= 180.f;
00767 }
00768
00769 Transform* t = result->get_attr("xform.align2d");
00770 t->set_rotation(Dict("type","2d","alpha",rotate_angle_solution));
00771 result->set_attr("xform.align2d",t);
00772 delete t;
00773
00774 return result;
00775 }
00776
00777
00778 EMData *RotateTranslateScaleAligner::align(EMData * this_img, EMData *to,
00779 const string & cmp_name, const Dict& cmp_params) const
00780 {
00781
00782
00783 basealigner_params["maxshift"] = params.set_default("maxshift", -1);
00784 basealigner_params["rfp_mode"] = params.set_default("rfp_mode",2);
00785 basealigner_params["useflcf"] = params.set_default("useflcf",0);
00786 basealigner_params["zscore"] = params.set_default("zscore",0);
00787
00788
00789 return align_using_base(this_img, to, cmp_name, cmp_params);
00790
00791 }
00792
00793 EMData* RotateTranslateFlipAligner::align(EMData * this_img, EMData *to,
00794 const string & cmp_name, const Dict& cmp_params) const
00795 {
00796
00797 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));
00798 EMData *rot_trans_align = this_img->align("rotate_translate",to,rt_params,cmp_name, cmp_params);
00799
00800
00801 EMData *flipped = params.set_default("flip", (EMData *) 0);
00802 bool delete_flag = false;
00803 if (flipped == 0) {
00804 flipped = to->process("xform.flip", Dict("axis", "x"));
00805 delete_flag = true;
00806 }
00807
00808 EMData * rot_trans_align_flip = this_img->align("rotate_translate", flipped, rt_params, cmp_name, cmp_params);
00809 Transform * t = rot_trans_align_flip->get_attr("xform.align2d");
00810 t->set_mirror(true);
00811 rot_trans_align_flip->set_attr("xform.align2d",t);
00812 delete t;
00813
00814
00815 float cmp1 = rot_trans_align->cmp(cmp_name, to, cmp_params);
00816 float cmp2 = rot_trans_align_flip->cmp(cmp_name, flipped, cmp_params);
00817
00818 if (delete_flag){
00819 if(flipped) {
00820 delete flipped;
00821 flipped = 0;
00822 }
00823 }
00824
00825 EMData *result = 0;
00826 if (cmp1 < cmp2 ) {
00827
00828 if( rot_trans_align_flip ) {
00829 delete rot_trans_align_flip;
00830 rot_trans_align_flip = 0;
00831 }
00832 result = rot_trans_align;
00833 }
00834 else {
00835 if( rot_trans_align ) {
00836 delete rot_trans_align;
00837 rot_trans_align = 0;
00838 }
00839 result = rot_trans_align_flip;
00840 result->process_inplace("xform.flip",Dict("axis","x"));
00841 }
00842
00843 return result;
00844 }
00845
00846 EMData *RotateTranslateFlipScaleAligner::align(EMData * this_img, EMData *to,
00847 const string & cmp_name, const Dict& cmp_params) const
00848 {
00849
00850
00851 basealigner_params["flip"] = params.set_default("flip", (EMData *) 0);
00852 basealigner_params["maxshift"] = params.set_default("maxshift", -1);
00853 basealigner_params["rfp_mode"] = params.set_default("rfp_mode",2);
00854 basealigner_params["useflcf"] = params.set_default("useflcf",0);
00855 basealigner_params["zscore"] = params.set_default("zscore",0);
00856
00857
00858 return align_using_base(this_img, to, cmp_name, cmp_params);
00859
00860 }
00861
00862 EMData* RotateTranslateFlipAlignerIterative::align(EMData * this_img, EMData *to,
00863 const string & cmp_name, const Dict& cmp_params) const
00864 {
00865
00866 Dict rt_params("maxshift", params["maxshift"],"r1",params.set_default("r1",-1),"r2",params.set_default("r2",-1));
00867 EMData *rot_trans_align = this_img->align("rotate_translate_iterative",to,rt_params,cmp_name, cmp_params);
00868
00869
00870 EMData *flipped = params.set_default("flip", (EMData *) 0);
00871 bool delete_flag = false;
00872 if (flipped == 0) {
00873 flipped = to->process("xform.flip", Dict("axis", "x"));
00874 delete_flag = true;
00875 }
00876
00877 EMData * rot_trans_align_flip = this_img->align("rotate_translate_iterative", flipped, rt_params, cmp_name, cmp_params);
00878 Transform* t = rot_trans_align_flip->get_attr("xform.align2d");
00879 t->set_mirror(true);
00880 rot_trans_align_flip->set_attr("xform.align2d",t);
00881 delete t;
00882
00883
00884 float cmp1 = rot_trans_align->cmp(cmp_name, to, cmp_params);
00885 float cmp2 = rot_trans_align_flip->cmp(cmp_name, flipped, cmp_params);
00886
00887 if (delete_flag){
00888 if(flipped) {
00889 delete flipped;
00890 flipped = 0;
00891 }
00892 }
00893
00894 EMData *result = 0;
00895 if (cmp1 < cmp2 ) {
00896
00897 if( rot_trans_align_flip ) {
00898 delete rot_trans_align_flip;
00899 rot_trans_align_flip = 0;
00900 }
00901 result = rot_trans_align;
00902 }
00903 else {
00904 if( rot_trans_align ) {
00905 delete rot_trans_align;
00906 rot_trans_align = 0;
00907 }
00908 result = rot_trans_align_flip;
00909 result->process_inplace("xform.flip",Dict("axis","x"));
00910 }
00911
00912 return result;
00913 }
00914
00915 EMData *RotateTranslateFlipScaleAlignerIterative::align(EMData * this_img, EMData *to,
00916 const string & cmp_name, const Dict& cmp_params) const
00917 {
00918
00919
00920 basealigner_params["flip"] = params.set_default("flip", (EMData *) 0);
00921 basealigner_params["maxshift"] = params.set_default("maxshift", -1);
00922 basealigner_params["r1"] = params.set_default("r1",-1);
00923 basealigner_params["r2"] = params.set_default("r2",-1);
00924 basealigner_params["maxiter"] = params.set_default("maxiter",3);
00925
00926
00927 return align_using_base(this_img, to, cmp_name, cmp_params);
00928
00929 }
00930
00931 EMData *RotateTranslateFlipAlignerPawel::align(EMData * this_img, EMData *to,
00932 const string & cmp_name, const Dict& cmp_params) const
00933 {
00934 if (cmp_name != "dot" && cmp_name != "ccc") throw InvalidParameterException("Resample aligner only works for dot and ccc");
00935
00936 int maxtx = params.set_default("tx", 0);
00937 int maxty = params.set_default("ty", 0);
00938 int r1 = params.set_default("r1",-1);
00939 int r2 = params.set_default("r2",-1);
00940
00941 if(this_img->get_xsize()/2 - 1 - r2 - maxtx <= 0 || (r2 == -1 && maxtx > 0)){
00942 cout << "\nRunTimeError: nx/2 - 1 - r2 - tx must be greater than or = 0\n" << endl;
00943 throw InvalidParameterException("nx/2 - 1 - r2 - tx must be greater than or = 0");
00944 }
00945 if(this_img->get_ysize()/2 - 1 - r2 - maxty <= 0 || (r2 == -1 && maxty > 0)){
00946 cout << "\nRunTimeError:ny/2 - 1 - r2 - ty must be greater than or = 0\n" << endl;
00947 throw InvalidParameterException("ny/2 - 1 - r2 - ty must be greater than or = 0");
00948 }
00949
00950
00951 float best_peak = -1.0e37;
00952 int best_peak_index = 0;
00953 int best_tx = 0;
00954 int best_ty = 0;
00955 int polarxsize = 0;
00956 bool flip = false;
00957
00958 for(int x = -maxtx; x <= maxtx; x++){
00959 for(int y = -maxty; y <= maxty; y++){
00960
00961 EMData * to_polar = to->unwrap(r1,r2,-1,0,0,true);
00962 EMData * this_img_polar = this_img->unwrap(r1,r2,-1,x,y,true);
00963 EMData * cfflip = this_img_polar->calc_ccfx(to_polar, 0, this_img_polar->get_ysize(), false, true);
00964 EMData * cf = this_img_polar->calc_ccfx(to_polar, 0, this_img_polar->get_ysize());
00965
00966 polarxsize = this_img_polar->get_xsize();
00967
00968
00969 delete to_polar; to_polar = 0;
00970 delete this_img_polar; this_img_polar = 0;
00971
00972 float *data = cf->get_data();
00973 float peak = 0;
00974 int peak_index = 0;
00975 Util::find_max(data, polarxsize, &peak, &peak_index);
00976 delete cf; cf = 0;
00977
00978 if(peak > best_peak) {
00979 best_peak = peak;
00980 best_peak_index = peak_index;
00981 best_tx = x;
00982 best_ty = y;
00983 flip = false;
00984 }
00985
00986 data = cfflip->get_data();
00987 Util::find_max(data, polarxsize, &peak, &peak_index);
00988 delete cfflip; cfflip = 0;
00989
00990 if(peak > best_peak) {
00991 best_peak = peak;
00992 best_peak_index = peak_index;
00993 best_tx = x;
00994 best_ty = y;
00995 flip = true;
00996 }
00997 }
00998 }
00999
01000 float rot_angle = (float) (best_peak_index * 360.0f / polarxsize);
01001
01002
01003 Transform tmptt(Dict("type","2d","alpha",0,"tx",-best_tx,"ty",-best_ty));
01004 Transform tmprot(Dict("type","2d","alpha",rot_angle,"tx",0,"ty",0));
01005 Transform total = tmprot*tmptt;
01006 EMData * rotimg=this_img->process("xform",Dict("transform",(Transform*)&total));
01007 rotimg->set_attr("xform.align2d",&total);
01008 if(flip == true) {
01009 rotimg->process_inplace("xform.flip",Dict("axis", "x"));
01010 }
01011
01012 return rotimg;
01013
01014 }
01015
01016 EMData *RotateFlipAligner::align(EMData * this_img, EMData *to,
01017 const string& cmp_name, const Dict& cmp_params) const
01018 {
01019 Dict rot_params("rfp_mode",params.set_default("rfp_mode",2));
01020 EMData *r1 = this_img->align("rotational", to, rot_params,cmp_name, cmp_params);
01021
01022
01023 EMData* flipped =to->process("xform.flip", Dict("axis", "x"));
01024 EMData *r2 = this_img->align("rotational", flipped,rot_params, cmp_name, cmp_params);
01025 Transform* t = r2->get_attr("xform.align2d");
01026 t->set_mirror(true);
01027 r2->set_attr("xform.align2d",t);
01028 delete t;
01029
01030 float cmp1 = r1->cmp(cmp_name, to, cmp_params);
01031 float cmp2 = r2->cmp(cmp_name, flipped, cmp_params);
01032
01033 delete flipped; flipped = 0;
01034
01035 EMData *result = 0;
01036
01037 if (cmp1 < cmp2) {
01038 if( r2 )
01039 {
01040 delete r2;
01041 r2 = 0;
01042 }
01043 result = r1;
01044 }
01045 else {
01046 if( r1 )
01047 {
01048 delete r1;
01049 r1 = 0;
01050 }
01051 result = r2;
01052 result->process_inplace("xform.flip",Dict("axis","x"));
01053 }
01054
01055 return result;
01056 }
01057
01058 EMData *RotateFlipAlignerIterative::align(EMData * this_img, EMData *to,
01059 const string& cmp_name, const Dict& cmp_params) const
01060 {
01061 Dict rot_params("r1",params.set_default("r1",-1),"r2",params.set_default("r2",-1));
01062 EMData * r1 = this_img->align("rotational_iterative", to, rot_params,cmp_name, cmp_params);
01063
01064 EMData * flipped =to->process("xform.flip", Dict("axis", "x"));
01065 EMData * r2 = this_img->align("rotational_iterative", flipped,rot_params, cmp_name, cmp_params);
01066 Transform* t = r2->get_attr("xform.align2d");
01067 t->set_mirror(true);
01068 r2->set_attr("xform.align2d",t);
01069 delete t;
01070
01071 float cmp1 = r1->cmp(cmp_name, to, cmp_params);
01072 float cmp2 = r2->cmp(cmp_name, flipped, cmp_params);
01073
01074 delete flipped; flipped = 0;
01075
01076 EMData *result = 0;
01077
01078 if (cmp1 < cmp2) {
01079 if( r2 )
01080 {
01081 delete r2;
01082 r2 = 0;
01083 }
01084 result = r1;
01085 }
01086 else {
01087 if( r1 )
01088 {
01089 delete r1;
01090 r1 = 0;
01091 }
01092 result = r2;
01093 result->process_inplace("xform.flip",Dict("axis","x"));
01094 }
01095
01096 return result;
01097 }
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113 EMData *RTFExhaustiveAligner::align(EMData * this_img, EMData *to,
01114 const string & cmp_name, const Dict& cmp_params) const
01115 {
01116 EMData *flip = params.set_default("flip", (EMData *) 0);
01117 int maxshift = params.set_default("maxshift", this_img->get_xsize()/8);
01118 if (maxshift < 2) throw InvalidParameterException("maxshift must be greater than or equal to 2");
01119
01120 int ny = this_img->get_ysize();
01121 int xst = (int) floor(2 * M_PI * ny);
01122 xst = Util::calc_best_fft_size(xst);
01123
01124 Dict d("n",2);
01125 EMData *to_shrunk_unwrapped = to->process("math.medianshrink",d);
01126
01127 int to_copy_r2 = to_shrunk_unwrapped->get_ysize() / 2 - 2 - maxshift / 2;
01128 EMData *tmp = to_shrunk_unwrapped->unwrap(4, to_copy_r2, xst / 2, 0, 0, true);
01129 if( to_shrunk_unwrapped )
01130 {
01131 delete to_shrunk_unwrapped;
01132 to_shrunk_unwrapped = 0;
01133 }
01134 to_shrunk_unwrapped = tmp;
01135
01136 EMData *to_shrunk_unwrapped_copy = to_shrunk_unwrapped->copy();
01137 EMData* to_unwrapped = to->unwrap(4, to->get_ysize() / 2 - 2 - maxshift, xst, 0, 0, true);
01138 EMData *to_unwrapped_copy = to_unwrapped->copy();
01139
01140 bool delete_flipped = true;
01141 EMData *flipped = 0;
01142 if (flip) {
01143 delete_flipped = false;
01144 flipped = flip;
01145 }
01146 else {
01147 flipped = to->process("xform.flip", Dict("axis", "x"));
01148 }
01149 EMData *to_shrunk_flipped_unwrapped = flipped->process("math.medianshrink",d);
01150 tmp = to_shrunk_flipped_unwrapped->unwrap(4, to_copy_r2, xst / 2, 0, 0, true);
01151 if( to_shrunk_flipped_unwrapped )
01152 {
01153 delete to_shrunk_flipped_unwrapped;
01154 to_shrunk_flipped_unwrapped = 0;
01155 }
01156 to_shrunk_flipped_unwrapped = tmp;
01157 EMData *to_shrunk_flipped_unwrapped_copy = to_shrunk_flipped_unwrapped->copy();
01158 EMData* to_flip_unwrapped = flipped->unwrap(4, to->get_ysize() / 2 - 2 - maxshift, xst, 0, 0, true);
01159 EMData* to_flip_unwrapped_copy = to_flip_unwrapped->copy();
01160
01161 if (delete_flipped && flipped != 0) {
01162 delete flipped;
01163 flipped = 0;
01164 }
01165
01166 EMData *this_shrunk_2 = this_img->process("math.medianshrink",d);
01167
01168 float bestval = FLT_MAX;
01169 float bestang = 0;
01170 int bestflip = 0;
01171 float bestdx = 0;
01172 float bestdy = 0;
01173
01174 int half_maxshift = maxshift / 2;
01175
01176 int ur2 = this_shrunk_2->get_ysize() / 2 - 2 - half_maxshift;
01177 for (int dy = -half_maxshift; dy <= half_maxshift; dy += 1) {
01178 for (int dx = -half_maxshift; dx <= half_maxshift; dx += 1) {
01179 #ifdef _WIN32
01180 if (_hypot(dx, dy) <= half_maxshift) {
01181 #else
01182 if (hypot(dx, dy) <= half_maxshift) {
01183 #endif
01184 EMData *uw = this_shrunk_2->unwrap(4, ur2, xst / 2, dx, dy, true);
01185 EMData *uwc = uw->copy();
01186 EMData *a = uw->calc_ccfx(to_shrunk_unwrapped);
01187
01188 uwc->rotate_x(a->calc_max_index());
01189 float cm = uwc->cmp(cmp_name, to_shrunk_unwrapped_copy, cmp_params);
01190 if (cm < bestval) {
01191 bestval = cm;
01192 bestang = (float) (2.0 * M_PI * a->calc_max_index() / a->get_xsize());
01193 bestdx = (float)dx;
01194 bestdy = (float)dy;
01195 bestflip = 0;
01196 }
01197
01198
01199 if( a )
01200 {
01201 delete a;
01202 a = 0;
01203 }
01204 if( uw )
01205 {
01206 delete uw;
01207 uw = 0;
01208 }
01209 if( uwc )
01210 {
01211 delete uwc;
01212 uwc = 0;
01213 }
01214 uw = this_shrunk_2->unwrap(4, ur2, xst / 2, dx, dy, true);
01215 uwc = uw->copy();
01216 a = uw->calc_ccfx(to_shrunk_flipped_unwrapped);
01217
01218 uwc->rotate_x(a->calc_max_index());
01219 cm = uwc->cmp(cmp_name, to_shrunk_flipped_unwrapped_copy, cmp_params);
01220 if (cm < bestval) {
01221 bestval = cm;
01222 bestang = (float) (2.0 * M_PI * a->calc_max_index() / a->get_xsize());
01223 bestdx = (float)dx;
01224 bestdy = (float)dy;
01225 bestflip = 1;
01226 }
01227
01228 if( a )
01229 {
01230 delete a;
01231 a = 0;
01232 }
01233
01234 if( uw )
01235 {
01236 delete uw;
01237 uw = 0;
01238 }
01239 if( uwc )
01240 {
01241 delete uwc;
01242 uwc = 0;
01243 }
01244 }
01245 }
01246 }
01247 if( this_shrunk_2 )
01248 {
01249 delete this_shrunk_2;
01250 this_shrunk_2 = 0;
01251 }
01252 if( to_shrunk_unwrapped )
01253 {
01254 delete to_shrunk_unwrapped;
01255 to_shrunk_unwrapped = 0;
01256 }
01257 if( to_shrunk_unwrapped_copy )
01258 {
01259 delete to_shrunk_unwrapped_copy;
01260 to_shrunk_unwrapped_copy = 0;
01261 }
01262 if( to_shrunk_flipped_unwrapped )
01263 {
01264 delete to_shrunk_flipped_unwrapped;
01265 to_shrunk_flipped_unwrapped = 0;
01266 }
01267 if( to_shrunk_flipped_unwrapped_copy )
01268 {
01269 delete to_shrunk_flipped_unwrapped_copy;
01270 to_shrunk_flipped_unwrapped_copy = 0;
01271 }
01272 bestdx *= 2;
01273 bestdy *= 2;
01274 bestval = FLT_MAX;
01275
01276 float bestdx2 = bestdx;
01277 float bestdy2 = bestdy;
01278
01279
01280
01281 for (float dy = bestdy2 - 3; dy <= bestdy2 + 3; dy += 1.0 ) {
01282 for (float dx = bestdx2 - 3; dx <= bestdx2 + 3; dx += 1.0 ) {
01283
01284 #ifdef _WIN32
01285 if (_hypot(dx, dy) <= maxshift) {
01286 #else
01287 if (hypot(dx, dy) <= maxshift) {
01288 #endif
01289 EMData *uw = this_img->unwrap(4, this_img->get_ysize() / 2 - 2 - maxshift, xst, (int)dx, (int)dy, true);
01290 EMData *uwc = uw->copy();
01291 EMData *a = uw->calc_ccfx(to_unwrapped);
01292
01293 uwc->rotate_x(a->calc_max_index());
01294 float cm = uwc->cmp(cmp_name, to_unwrapped_copy, cmp_params);
01295
01296 if (cm < bestval) {
01297 bestval = cm;
01298 bestang = (float)(2.0 * M_PI * a->calc_max_index() / a->get_xsize());
01299 bestdx = dx;
01300 bestdy = dy;
01301 bestflip = 0;
01302 }
01303
01304 if( a )
01305 {
01306 delete a;
01307 a = 0;
01308 }
01309 if( uw )
01310 {
01311 delete uw;
01312 uw = 0;
01313 }
01314 if( uwc )
01315 {
01316 delete uwc;
01317 uwc = 0;
01318 }
01319 uw = this_img->unwrap(4, this_img->get_ysize() / 2 - 2 - maxshift, xst, (int)dx, (int)dy, true);
01320 uwc = uw->copy();
01321 a = uw->calc_ccfx(to_flip_unwrapped);
01322
01323 uwc->rotate_x(a->calc_max_index());
01324 cm = uwc->cmp(cmp_name, to_flip_unwrapped_copy, cmp_params);
01325
01326 if (cm < bestval) {
01327 bestval = cm;
01328 bestang = (float)(2.0 * M_PI * a->calc_max_index() / a->get_xsize());
01329 bestdx = dx;
01330 bestdy = dy;
01331 bestflip = 1;
01332 }
01333
01334 if( a )
01335 {
01336 delete a;
01337 a = 0;
01338 }
01339 if( uw )
01340 {
01341 delete uw;
01342 uw = 0;
01343 }
01344 if( uwc )
01345 {
01346 delete uwc;
01347 uwc = 0;
01348 }
01349 }
01350 }
01351 }
01352 if( to_unwrapped ) {delete to_unwrapped;to_unwrapped = 0;}
01353 if( to_shrunk_unwrapped ) { delete to_shrunk_unwrapped; to_shrunk_unwrapped = 0;}
01354 if (to_unwrapped_copy) { delete to_unwrapped_copy; to_unwrapped_copy = 0; }
01355 if (to_flip_unwrapped) { delete to_flip_unwrapped; to_flip_unwrapped = 0; }
01356 if (to_flip_unwrapped_copy) { delete to_flip_unwrapped_copy; to_flip_unwrapped_copy = 0;}
01357
01358 bestang *= (float)EMConsts::rad2deg;
01359 Transform t(Dict("type","2d","alpha",(float)bestang));
01360 t.set_pre_trans(Vec2f(-bestdx,-bestdy));
01361 if (bestflip) {
01362 t.set_mirror(true);
01363 }
01364
01365 EMData* ret = this_img->process("xform",Dict("transform",&t));
01366 ret->set_attr("xform.align2d",&t);
01367
01368 return ret;
01369 }
01370
01371
01372 EMData *RTFSlowExhaustiveAligner::align(EMData * this_img, EMData *to,
01373 const string & cmp_name, const Dict& cmp_params) const
01374 {
01375
01376 EMData *flip = params.set_default("flip", (EMData *) 0);
01377 int maxshift = params.set_default("maxshift", -1);
01378
01379 EMData *flipped = 0;
01380
01381 bool delete_flipped = true;
01382 if (flip) {
01383 delete_flipped = false;
01384 flipped = flip;
01385 }
01386 else {
01387 flipped = to->process("xform.flip", Dict("axis", "x"));
01388 }
01389
01390 int nx = this_img->get_xsize();
01391
01392 if (maxshift < 0) {
01393 maxshift = nx / 10;
01394 }
01395
01396 float angle_step = params.set_default("angstep", 0.0f);
01397 if ( angle_step == 0 ) angle_step = atan2(2.0f, (float)nx);
01398 else {
01399 angle_step *= (float)EMConsts::deg2rad;
01400 }
01401 float trans_step = params.set_default("transtep",1.0f);
01402
01403 if (trans_step <= 0) throw InvalidParameterException("transstep must be greater than 0");
01404 if (angle_step <= 0) throw InvalidParameterException("angstep must be greater than 0");
01405
01406
01407 Dict shrinkfactor("n",2);
01408 EMData *this_img_shrink = this_img->process("math.medianshrink",shrinkfactor);
01409 EMData *to_shrunk = to->process("math.medianshrink",shrinkfactor);
01410 EMData *flipped_shrunk = flipped->process("math.medianshrink",shrinkfactor);
01411
01412 int bestflip = 0;
01413 float bestdx = 0;
01414 float bestdy = 0;
01415
01416 float bestang = 0;
01417 float bestval = FLT_MAX;
01418
01419 int half_maxshift = maxshift / 2;
01420
01421
01422 for (int dy = -half_maxshift; dy <= half_maxshift; ++dy) {
01423 for (int dx = -half_maxshift; dx <= half_maxshift; ++dx) {
01424 if (hypot(dx, dy) <= maxshift) {
01425 for (float ang = -angle_step * 2.0f; ang <= (float)2 * M_PI; ang += angle_step * 4.0f) {
01426 EMData v(*this_img_shrink);
01427 Transform t(Dict("type","2d","alpha",static_cast<float>(ang*EMConsts::rad2deg)));
01428 t.set_trans((float)dx,(float)dy);
01429 v.transform(t);
01430
01431
01432 float lc = v.cmp(cmp_name, to_shrunk, cmp_params);
01433
01434 if (lc < bestval) {
01435 bestval = lc;
01436 bestang = ang;
01437 bestdx = (float)dx;
01438 bestdy = (float)dy;
01439 bestflip = 0;
01440 }
01441
01442 lc = v.cmp(cmp_name,flipped_shrunk , cmp_params);
01443 if (lc < bestval) {
01444 bestval = lc;
01445 bestang = ang;
01446 bestdx = (float)dx;
01447 bestdy = (float)dy;
01448 bestflip = 1;
01449 }
01450 }
01451 }
01452 }
01453 }
01454
01455 if( to_shrunk )
01456 {
01457 delete to_shrunk;
01458 to_shrunk = 0;
01459 }
01460 if( flipped_shrunk )
01461 {
01462 delete flipped_shrunk;
01463 flipped_shrunk = 0;
01464 }
01465 if( this_img_shrink )
01466 {
01467 delete this_img_shrink;
01468 this_img_shrink = 0;
01469 }
01470
01471 bestdx *= 2;
01472 bestdy *= 2;
01473 bestval = FLT_MAX;
01474
01475 float bestdx2 = bestdx;
01476 float bestdy2 = bestdy;
01477 float bestang2 = bestang;
01478
01479 for (float dy = bestdy2 - 3; dy <= bestdy2 + 3; dy += trans_step) {
01480 for (float dx = bestdx2 - 3; dx <= bestdx2 + 3; dx += trans_step) {
01481 if (hypot(dx, dy) <= maxshift) {
01482 for (float ang = bestang2 - angle_step * 6.0f; ang <= bestang2 + angle_step * 6.0f; ang += angle_step) {
01483 EMData v(*this_img);
01484 Transform t(Dict("type","2d","alpha",static_cast<float>(ang*EMConsts::rad2deg)));
01485 t.set_trans(dx,dy);
01486 v.transform(t);
01487
01488
01489 float lc = v.cmp(cmp_name, to, cmp_params);
01490
01491 if (lc < bestval) {
01492 bestval = lc;
01493 bestang = ang;
01494 bestdx = dx;
01495 bestdy = dy;
01496 bestflip = 0;
01497 }
01498
01499 lc = v.cmp(cmp_name, flipped, cmp_params);
01500
01501 if (lc < bestval) {
01502 bestval = lc;
01503 bestang = ang;
01504 bestdx = dx;
01505 bestdy = dy;
01506 bestflip = 1;
01507 }
01508 }
01509 }
01510 }
01511 }
01512
01513 if (delete_flipped) { delete flipped; flipped = 0; }
01514
01515 bestang *= (float)EMConsts::rad2deg;
01516 Transform t(Dict("type","2d","alpha",(float)bestang));
01517 t.set_trans(bestdx,bestdy);
01518
01519 if (bestflip) {
01520 t.set_mirror(true);
01521 }
01522
01523 EMData* rslt = this_img->process("xform",Dict("transform",&t));
01524 rslt->set_attr("xform.align2d",&t);
01525
01526 return rslt;
01527 }
01528
01529 EMData* SymAlignProcessor::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
01530 {
01531
01532
01533 float dphi = params.set_default("dphi",10.f);
01534 float lphi = params.set_default("lphi",0.0f);
01535 float uphi = params.set_default("uphi",359.9f);
01536
01537 Dict d;
01538 d["inc_mirror"] = true;
01539 d["delta"] = params.set_default("delta",10.f);
01540
01541
01542 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
01543 vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
01544
01545
01546 vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
01547
01548 float bestquality = 0.0f;
01549 EMData* bestimage = 0;
01550 for(vector<Transform>::const_iterator trans_it = transforms.begin(); trans_it != transforms.end(); trans_it++) {
01551 Dict tparams = trans_it->get_params("eman");
01552 Transform t(tparams);
01553 for( float phi = lphi; phi < uphi; phi += dphi ) {
01554 tparams["phi"] = phi;
01555 t.set_rotation(tparams);
01556
01557
01558 Averager* imgavg = Factory<Averager>::get((string)params.set_default("avger","mean"));
01559
01560 for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
01561 Transform sympos = (*it)*t;
01562 EMData* transformed = this_img->process("xform",Dict("transform",&sympos));
01563 imgavg->add_image(transformed);
01564 delete transformed;
01565 }
01566
01567 EMData* symptcl=imgavg->finish();
01568 delete imgavg;
01569
01570 float quality = symptcl->get_attr("sigma");
01571 cout << quality << " " << phi << endl;
01572 if(quality > bestquality) {
01573 bestquality = quality;
01574 bestimage = symptcl;
01575 } else {
01576 delete symptcl;
01577 }
01578 }
01579 }
01580 if(sym != 0) delete sym;
01581
01582 return bestimage;
01583 }
01584
01585 static double refalifn(const gsl_vector * v, void *params)
01586 {
01587 Dict *dict = (Dict *) params;
01588
01589 double x = gsl_vector_get(v, 0);
01590 double y = gsl_vector_get(v, 1);
01591 double a = gsl_vector_get(v, 2);
01592
01593 EMData *this_img = (*dict)["this"];
01594 EMData *with = (*dict)["with"];
01595 bool mirror = (*dict)["mirror"];
01596
01597 Transform t(Dict("type","2d","alpha",static_cast<float>(a)));
01598 t.set_trans((float)x,(float)y);
01599 t.set_mirror(mirror);
01600 if (v->size>3) {
01601 float sca=(float)gsl_vector_get(v, 3);
01602 if (sca<.7 || sca>1.3) return 1.0e20;
01603 t.set_scale((float)gsl_vector_get(v, 3));
01604 }
01605 EMData *tmp = this_img->process("xform",Dict("transform",&t));
01606 if (dict->has_key("mask")) tmp->mult(*(EMData *)((*dict)["mask"]));
01607
01608
01609 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
01610 double result = c->cmp(tmp,with);
01611
01612 if (tmp != 0) delete tmp;
01613
01614 return result;
01615 }
01616
01617 static void refalidf(const gsl_vector * v, void *params,gsl_vector * df) {
01618
01619
01620
01621 static double lstep[4] = { 0.05, 0.05, 0.1, 0.01 };
01622
01623 gsl_vector *vc = gsl_vector_alloc(v->size);
01624 gsl_vector_memcpy(vc,v);
01625
01626 double f = refalifn(v,params);
01627 for (unsigned int i=0; i<v->size; i++) {
01628 double *vp = gsl_vector_ptr(vc,i);
01629 *vp+=lstep[i];
01630 double f2 = refalifn(vc,params);
01631 *vp-=lstep[i];
01632
01633 gsl_vector_set(df,i,(f2-f)/lstep[i]);
01634 }
01635
01636 gsl_vector_free(vc);
01637 return;
01638 }
01639
01640 static void refalifdf(const gsl_vector * v, void *params, double * f, gsl_vector * df) {
01641
01642
01643
01644 static double lstep[4] = { 0.05, 0.05, 0.1, 0.01 };
01645
01646 gsl_vector *vc = gsl_vector_alloc(v->size);
01647 gsl_vector_memcpy(vc,v);
01648
01649 *f = refalifn(v,params);
01650 for (unsigned int i=0; i<v->size; i++) {
01651 double *vp = gsl_vector_ptr(vc,i);
01652 *vp+=lstep[i];
01653 double f2 = refalifn(vc,params);
01654 *vp-=lstep[i];
01655
01656 gsl_vector_set(df,i,(f2-*f)/lstep[i]);
01657 }
01658
01659 gsl_vector_free(vc);
01660 return;
01661
01662 }
01663
01664 static double refalifnfast(const gsl_vector * v, void *params)
01665 {
01666 Dict *dict = (Dict *) params;
01667 EMData *this_img = (*dict)["this"];
01668 EMData *img_to = (*dict)["with"];
01669 bool mirror = (*dict)["mirror"];
01670
01671 double x = gsl_vector_get(v, 0);
01672 double y = gsl_vector_get(v, 1);
01673 double a = gsl_vector_get(v, 2);
01674
01675 double r = this_img->dot_rotate_translate(img_to, (float)x, (float)y, (float)a, mirror);
01676 int nsec = this_img->get_xsize() * this_img->get_ysize();
01677 double result = 1.0 - r / nsec;
01678
01679
01680 return result;
01681 }
01682
01683 EMData *RefineAligner::align(EMData * this_img, EMData *to,
01684 const string & cmp_name, const Dict& cmp_params) const
01685 {
01686
01687 if (!to) {
01688 return 0;
01689 }
01690
01691 EMData *result;
01692 int mode = params.set_default("mode", 0);
01693 float saz = 0.0;
01694 float sdx = 0.0;
01695 float sdy = 0.0;
01696 float sscale = 1.0;
01697 bool mirror = false;
01698 Transform* t;
01699 if (params.has_key("xform.align2d") ) {
01700 t = params["xform.align2d"];
01701 Dict params = t->get_params("2d");
01702 saz = params["alpha"];
01703 sdx = params["tx"];
01704 sdy = params["ty"];
01705 mirror = params["mirror"];
01706 sscale = params["scale"];
01707 } else {
01708 t = new Transform();
01709 }
01710
01711
01712 if ((float)(this_img->get_attr("sigma"))==0.0 || (float)(to->get_attr("sigma"))==0.0) {
01713 result = this_img->process("xform",Dict("transform",t));
01714 result->set_attr("xform.align2d",t);
01715 delete t;
01716 return result;
01717 }
01718
01719 float stepx = params.set_default("stepx",1.0f);
01720 float stepy = params.set_default("stepy",1.0f);
01721
01722 float stepaz = params.set_default("stepaz",5.0f);
01723 float stepscale = params.set_default("stepscale",0.0f);
01724
01725 int np = 3;
01726 if (stepscale!=0.0) np++;
01727 Dict gsl_params;
01728 gsl_params["this"] = this_img;
01729 gsl_params["with"] = to;
01730 gsl_params["snr"] = params["snr"];
01731 gsl_params["mirror"] = mirror;
01732 if (params.has_key("mask")) gsl_params["mask"]=params["mask"];
01733
01734 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
01735 gsl_vector *ss = gsl_vector_alloc(np);
01736
01737
01738 gsl_vector_set(ss, 0, stepx);
01739 gsl_vector_set(ss, 1, stepy);
01740 gsl_vector_set(ss, 2, stepaz);
01741 if (stepscale!=0.0) gsl_vector_set(ss,3,stepscale);
01742
01743 gsl_vector *x = gsl_vector_alloc(np);
01744 gsl_vector_set(x, 0, sdx);
01745 gsl_vector_set(x, 1, sdy);
01746 gsl_vector_set(x, 2, saz);
01747 if (stepscale!=0.0) gsl_vector_set(x,3,1.0);
01748
01749 Cmp *c = 0;
01750
01751 gsl_multimin_function minex_func;
01752 if (mode == 2) {
01753 minex_func.f = &refalifnfast;
01754 }
01755 else {
01756 c = Factory < Cmp >::get(cmp_name, cmp_params);
01757 gsl_params["cmp"] = (void *) c;
01758 minex_func.f = &refalifn;
01759 }
01760
01761 minex_func.n = np;
01762 minex_func.params = (void *) &gsl_params;
01763
01764 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
01765 gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
01766
01767 int rval = GSL_CONTINUE;
01768 int status = GSL_SUCCESS;
01769 int iter = 1;
01770
01771 float precision = params.set_default("precision",0.04f);
01772 int maxiter = params.set_default("maxiter",28);
01773
01774
01775
01776
01777 while (rval == GSL_CONTINUE && iter < maxiter) {
01778 iter++;
01779 status = gsl_multimin_fminimizer_iterate(s);
01780 if (status) {
01781 break;
01782 }
01783 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
01784 }
01785
01786 int maxshift = params.set_default("maxshift",-1);
01787
01788 if (maxshift <= 0) {
01789 maxshift = this_img->get_xsize() / 4;
01790 }
01791 float fmaxshift = static_cast<float>(maxshift);
01792 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)) )
01793 {
01794
01795 Transform tsoln(Dict("type","2d","alpha",(float)gsl_vector_get(s->x, 2)));
01796 tsoln.set_mirror(mirror);
01797 tsoln.set_trans((float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1));
01798 if (stepscale!=0.0) tsoln.set_scale((float)gsl_vector_get(s->x, 3));
01799 result = this_img->process("xform",Dict("transform",&tsoln));
01800 result->set_attr("xform.align2d",&tsoln);
01801 } else {
01802
01803 result = this_img->process("xform",Dict("transform",t));
01804 result->set_attr("xform.align2d",t);
01805 }
01806
01807 delete t;
01808 t = 0;
01809
01810 gsl_vector_free(x);
01811 gsl_vector_free(ss);
01812 gsl_multimin_fminimizer_free(s);
01813
01814 if (c != 0) delete c;
01815 return result;
01816 }
01817
01818 EMData *RefineAlignerCG::align(EMData * this_img, EMData *to,
01819 const string & cmp_name, const Dict& cmp_params) const
01820 {
01821
01822 if (!to) {
01823 return 0;
01824 }
01825
01826 EMData *result;
01827 int mode = params.set_default("mode", 0);
01828 float saz = 0.0;
01829 float sdx = 0.0;
01830 float sdy = 0.0;
01831 float sscale = 1.0;
01832 bool mirror = false;
01833 Transform* t;
01834 if (params.has_key("xform.align2d") ) {
01835 t = params["xform.align2d"];
01836 Dict params = t->get_params("2d");
01837 saz = params["alpha"];
01838 sdx = params["tx"];
01839 sdy = params["ty"];
01840 mirror = params["mirror"];
01841 sscale = params["scale"];
01842 } else {
01843 t = new Transform();
01844 }
01845
01846
01847 if ((float)(this_img->get_attr("sigma"))==0.0 || (float)(to->get_attr("sigma"))==0.0) {
01848 result = this_img->process("xform",Dict("transform",t));
01849 result->set_attr("xform.align2d",t);
01850 delete t;
01851 return result;
01852 }
01853
01854 float step = params.set_default("step",0.1f);
01855 float stepscale = params.set_default("stepscale",0.0f);
01856
01857 int np = 3;
01858 if (stepscale!=0.0) np++;
01859 Dict gsl_params;
01860 gsl_params["this"] = this_img;
01861 gsl_params["with"] = to;
01862 gsl_params["snr"] = params["snr"];
01863 gsl_params["mirror"] = mirror;
01864 if (params.has_key("mask")) gsl_params["mask"]=params["mask"];
01865
01866 const gsl_multimin_fdfminimizer_type *T = gsl_multimin_fdfminimizer_vector_bfgs;
01867
01868 gsl_vector *x = gsl_vector_alloc(np);
01869 gsl_vector_set(x, 0, sdx);
01870 gsl_vector_set(x, 1, sdy);
01871 gsl_vector_set(x, 2, saz);
01872 if (stepscale!=0.0) gsl_vector_set(x,3,1.0);
01873
01874 Cmp *c = 0;
01875
01876 gsl_multimin_function_fdf minex_func;
01877 if (mode == 2) {
01878 minex_func.f = &refalifnfast;
01879 }
01880 else {
01881 c = Factory < Cmp >::get(cmp_name, cmp_params);
01882 gsl_params["cmp"] = (void *) c;
01883 minex_func.f = &refalifn;
01884 }
01885
01886 minex_func.df = &refalidf;
01887 minex_func.fdf = &refalifdf;
01888 minex_func.n = np;
01889 minex_func.params = (void *) &gsl_params;
01890
01891 gsl_multimin_fdfminimizer *s = gsl_multimin_fdfminimizer_alloc(T, np);
01892 gsl_multimin_fdfminimizer_set(s, &minex_func, x, step, 0.001f);
01893
01894 int rval = GSL_CONTINUE;
01895 int status = GSL_SUCCESS;
01896 int iter = 1;
01897
01898 float precision = params.set_default("precision",0.02f);
01899 int maxiter = params.set_default("maxiter",12);
01900 int verbose = params.set_default("verbose",0);
01901
01902
01903
01904
01905 while (rval == GSL_CONTINUE && iter < maxiter) {
01906 iter++;
01907 status = gsl_multimin_fdfminimizer_iterate(s);
01908 if (status) {
01909 break;
01910 }
01911 rval = gsl_multimin_test_gradient (s->gradient, precision);
01912
01913 }
01914
01915 int maxshift = params.set_default("maxshift",-1);
01916
01917 if (maxshift <= 0) {
01918 maxshift = this_img->get_xsize() / 4;
01919 }
01920 float fmaxshift = static_cast<float>(maxshift);
01921 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)) )
01922 {
01923 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));
01924 Transform tsoln(Dict("type","2d","alpha",(float)gsl_vector_get(s->x, 2)));
01925 tsoln.set_mirror(mirror);
01926 tsoln.set_trans((float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1));
01927 if (stepscale!=0.0) tsoln.set_scale((float)gsl_vector_get(s->x, 3));
01928 result = this_img->process("xform",Dict("transform",&tsoln));
01929 result->set_attr("xform.align2d",&tsoln);
01930 } else {
01931 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));
01932 result = this_img->process("xform",Dict("transform",t));
01933 result->set_attr("xform.align2d",t);
01934 }
01935
01936 delete t;
01937 t = 0;
01938
01939 gsl_vector_free(x);
01940 gsl_multimin_fdfminimizer_free(s);
01941
01942 if (c != 0) delete c;
01943 return result;
01944 }
01945
01946 static 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)
01947 {
01948 Vec3f normal(n0,n1,n2);
01949 normal.normalize();
01950
01951 float omega = spincoeff*sqrt(n0*n0 + n1*n1 + n2*n2);
01952 Dict d;
01953 d["type"] = "spin";
01954 d["omega"] = omega;
01955 d["n1"] = normal[0];
01956 d["n2"] = normal[1];
01957 d["n3"] = normal[2];
01958
01959
01960 Transform q(d);
01961 q.set_trans((float)x,(float)y,(float)z);
01962
01963 q = q*(*t);
01964
01965 return q;
01966 }
01967
01968 static double symquat(const gsl_vector * v, void *params)
01969 {
01970 Dict *dict = (Dict *) params;
01971
01972 double n0 = gsl_vector_get(v, 0);
01973 double n1 = gsl_vector_get(v, 1);
01974 double n2 = gsl_vector_get(v, 2);
01975 double x = gsl_vector_get(v, 3);
01976 double y = gsl_vector_get(v, 4);
01977 double z = gsl_vector_get(v, 5);
01978
01979 EMData* volume = (*dict)["volume"];
01980 float spincoeff = (*dict)["spincoeff"];
01981 Transform* t = (*dict)["transform"];
01982
01983 Transform soln = refalin3d_perturbquat(t,spincoeff,(float)n0,(float)n1,(float)n2,(float)x,(float)y,(float)z);
01984
01985 EMData *tmp = volume->process("xform",Dict("transform",&soln));
01986 EMData *symtmp = tmp->process("xform.applysym",Dict("sym",(*dict)["sym"]));
01987 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
01988 double result = c->cmp(symtmp,tmp);
01989 delete tmp;
01990 delete symtmp;
01991
01992
01993 return result;
01994 }
01995
01996 static double refalifn3dquat(const gsl_vector * v, void *params)
01997 {
01998 Dict *dict = (Dict *) params;
01999
02000 double n0 = gsl_vector_get(v, 0);
02001 double n1 = gsl_vector_get(v, 1);
02002 double n2 = gsl_vector_get(v, 2);
02003 double x = gsl_vector_get(v, 3);
02004 double y = gsl_vector_get(v, 4);
02005 double z = gsl_vector_get(v, 5);
02006
02007 EMData *this_img = (*dict)["this"];
02008 EMData *with = (*dict)["with"];
02009
02010 Transform* t = (*dict)["transform"];
02011 float spincoeff = (*dict)["spincoeff"];
02012
02013 Transform soln = refalin3d_perturbquat(t,spincoeff,(float)n0,(float)n1,(float)n2,(float)x,(float)y,(float)z);
02014
02015 EMData *tmp = this_img->process("xform",Dict("transform",&soln));
02016 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
02017 double result = c->cmp(tmp,with);
02018 if ( tmp != 0 ) delete tmp;
02019
02020
02021 return result;
02022 }
02023
02024 EMData* SymAlignProcessorQuat::align(EMData * volume, EMData *to, const string & cmp_name, const Dict& cmp_params) const
02025 {
02026
02027 Transform* t;
02028 if (params.has_key("xform.align3d") ) {
02029 t = params["xform.align3d"];
02030 }else {
02031 t = new Transform();
02032 }
02033
02034 float sdi = 0.0;
02035 float sdj = 0.0;
02036 float sdk = 0.0;
02037 float sdx = 0.0;
02038 float sdy = 0.0;
02039 float sdz = 0.0;
02040
02041 float spincoeff = params.set_default("spin_coeff",10.0f);
02042
02043 int np = 6;
02044 Dict gsl_params;
02045 gsl_params["volume"] = volume;
02046 gsl_params["transform"] = t;
02047 gsl_params["sym"] = params.set_default("sym","c1");
02048 gsl_params["spincoeff"] = spincoeff;
02049
02050 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
02051 gsl_vector *ss = gsl_vector_alloc(np);
02052
02053 float stepi = params.set_default("stepn0",1.0f);
02054 float stepj = params.set_default("stepn1",1.0f);
02055 float stepk = params.set_default("stepn2",1.0f);
02056 float stepx = params.set_default("stepx",1.0f);
02057 float stepy = params.set_default("stepy",1.0f);
02058 float stepz = params.set_default("stepz",1.0f);
02059
02060 gsl_vector_set(ss, 0, stepi);
02061 gsl_vector_set(ss, 1, stepj);
02062 gsl_vector_set(ss, 2, stepk);
02063 gsl_vector_set(ss, 3, stepx);
02064 gsl_vector_set(ss, 4, stepy);
02065 gsl_vector_set(ss, 5, stepz);
02066
02067 gsl_vector *x = gsl_vector_alloc(np);
02068 gsl_vector_set(x, 0, sdi);
02069 gsl_vector_set(x, 1, sdj);
02070 gsl_vector_set(x, 2, sdk);
02071 gsl_vector_set(x, 3, sdx);
02072 gsl_vector_set(x, 4, sdy);
02073 gsl_vector_set(x, 5, sdz);
02074
02075 gsl_multimin_function minex_func;
02076 Cmp *c = Factory < Cmp >::get(cmp_name, cmp_params);
02077 gsl_params["cmp"] = (void *) c;
02078 minex_func.f = &symquat;
02079 minex_func.n = np;
02080 minex_func.params = (void *) &gsl_params;
02081 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
02082 gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
02083
02084 int rval = GSL_CONTINUE;
02085 int status = GSL_SUCCESS;
02086 int iter = 1;
02087
02088 float precision = params.set_default("precision",0.01f);
02089 int maxiter = params.set_default("maxiter",100);
02090 while (rval == GSL_CONTINUE && iter < maxiter) {
02091 iter++;
02092 status = gsl_multimin_fminimizer_iterate(s);
02093 if (status) {
02094 break;
02095 }
02096 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
02097 }
02098
02099 int maxshift = params.set_default("maxshift",-1);
02100
02101 if (maxshift <= 0) {
02102 maxshift = volume->get_xsize() / 4;
02103 }
02104 float fmaxshift = static_cast<float>(maxshift);
02105
02106 EMData *result;
02107 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))
02108 {
02109 float n0 = (float)gsl_vector_get(s->x, 0);
02110 float n1 = (float)gsl_vector_get(s->x, 1);
02111 float n2 = (float)gsl_vector_get(s->x, 2);
02112 float x = (float)gsl_vector_get(s->x, 3);
02113 float y = (float)gsl_vector_get(s->x, 4);
02114 float z = (float)gsl_vector_get(s->x, 5);
02115
02116 Transform tsoln = refalin3d_perturbquat(t,spincoeff,n0,n1,n2,x,y,z);
02117
02118 result = volume->process("xform",Dict("transform",&tsoln));
02119 result->set_attr("xform.align3d",&tsoln);
02120 EMData *tmpsym = result->process("xform.applysym",Dict("sym",gsl_params["sym"]));
02121 result->set_attr("score", result->cmp(cmp_name,tmpsym,cmp_params));
02122 delete tmpsym;
02123 } else {
02124 result = volume->process("xform",Dict("transform",t));
02125 result->set_attr("xform.align3d",t);
02126 result->set_attr("score",0.0);
02127 }
02128
02129 gsl_vector_free(x);
02130 gsl_vector_free(ss);
02131 gsl_multimin_fminimizer_free(s);
02132
02133 if (c != 0) delete c;
02134 delete t;
02135
02136 return result;
02137 }
02138
02139 EMData* Refine3DAlignerQuaternion::align(EMData * this_img, EMData *to,
02140 const string & cmp_name, const Dict& cmp_params) const
02141 {
02142
02143 if (!to || !this_img) throw NullPointerException("Input image is null");
02144
02145 if (to->get_ndim() != 3 || this_img->get_ndim() != 3) throw ImageDimensionException("The Refine3D aligner only works for 3D images");
02146
02147 #ifdef EMAN2_USING_CUDA
02148 if(EMData::usecuda == 1) {
02149 if(!this_img->getcudarwdata()) this_img->copy_to_cuda();
02150 if(!to->getcudarwdata()) to->copy_to_cuda();
02151 }
02152 #endif
02153
02154 float sdi = 0.0;
02155 float sdj = 0.0;
02156 float sdk = 0.0;
02157 float sdx = 0.0;
02158 float sdy = 0.0;
02159 float sdz = 0.0;
02160 bool mirror = false;
02161
02162 Transform* t;
02163 if (params.has_key("xform.align3d") ) {
02164
02165
02166
02167
02168 t = params["xform.align3d"];
02169 }else {
02170 t = new Transform();
02171 }
02172
02173 float spincoeff = params.set_default("spin_coeff",10.0f);
02174
02175 int np = 6;
02176 Dict gsl_params;
02177 gsl_params["this"] = this_img;
02178 gsl_params["with"] = to;
02179 gsl_params["snr"] = params["snr"];
02180 gsl_params["mirror"] = mirror;
02181 gsl_params["transform"] = t;
02182 gsl_params["spincoeff"] = spincoeff;
02183 Dict altered_cmp_params(cmp_params);
02184
02185 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
02186 gsl_vector *ss = gsl_vector_alloc(np);
02187
02188 float stepi = params.set_default("stepn0",1.0f);
02189 float stepj = params.set_default("stepn1",1.0f);
02190 float stepk = params.set_default("stepn2",1.0f);
02191 float stepx = params.set_default("stepx",1.0f);
02192 float stepy = params.set_default("stepy",1.0f);
02193 float stepz = params.set_default("stepz",1.0f);
02194
02195
02196 gsl_vector_set(ss, 0, stepi);
02197 gsl_vector_set(ss, 1, stepj);
02198 gsl_vector_set(ss, 2, stepk);
02199 gsl_vector_set(ss, 3, stepx);
02200 gsl_vector_set(ss, 4, stepy);
02201 gsl_vector_set(ss, 5, stepz);
02202
02203 gsl_vector *x = gsl_vector_alloc(np);
02204 gsl_vector_set(x, 0, sdi);
02205 gsl_vector_set(x, 1, sdj);
02206 gsl_vector_set(x, 2, sdk);
02207 gsl_vector_set(x, 3, sdx);
02208 gsl_vector_set(x, 4, sdy);
02209 gsl_vector_set(x, 5, sdz);
02210
02211 gsl_multimin_function minex_func;
02212 Cmp *c = Factory < Cmp >::get(cmp_name, altered_cmp_params);
02213
02214 gsl_params["cmp"] = (void *) c;
02215 minex_func.f = &refalifn3dquat;
02216
02217 minex_func.n = np;
02218 minex_func.params = (void *) &gsl_params;
02219
02220 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
02221 gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
02222
02223 int rval = GSL_CONTINUE;
02224 int status = GSL_SUCCESS;
02225 int iter = 1;
02226
02227 float precision = params.set_default("precision",0.01f);
02228 int maxiter = params.set_default("maxiter",100);
02229 while (rval == GSL_CONTINUE && iter < maxiter) {
02230 iter++;
02231 status = gsl_multimin_fminimizer_iterate(s);
02232 if (status) {
02233 break;
02234 }
02235 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
02236 }
02237
02238 int maxshift = params.set_default("maxshift",-1);
02239
02240 if (maxshift <= 0) {
02241 maxshift = this_img->get_xsize() / 4;
02242 }
02243 float fmaxshift = static_cast<float>(maxshift);
02244
02245 EMData *result;
02246 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))
02247 {
02248 float n0 = (float)gsl_vector_get(s->x, 0);
02249 float n1 = (float)gsl_vector_get(s->x, 1);
02250 float n2 = (float)gsl_vector_get(s->x, 2);
02251 float x = (float)gsl_vector_get(s->x, 3);
02252 float y = (float)gsl_vector_get(s->x, 4);
02253 float z = (float)gsl_vector_get(s->x, 5);
02254
02255 Transform tsoln = refalin3d_perturbquat(t,spincoeff,n0,n1,n2,x,y,z);
02256
02257 result = this_img->process("xform",Dict("transform",&tsoln));
02258 result->set_attr("xform.align3d",&tsoln);
02259 result->set_attr("score", result->cmp(cmp_name,to,cmp_params));
02260
02261
02262 } else {
02263 result = this_img->process("xform",Dict("transform",t));
02264 result->set_attr("xform.align3d",t);
02265 result->set_attr("score",0.0);
02266 }
02267
02268
02269 delete t;
02270 gsl_vector_free(x);
02271 gsl_vector_free(ss);
02272 gsl_multimin_fminimizer_free(s);
02273
02274 if (c != 0) delete c;
02275
02276 return result;
02277 }
02278
02279 EMData* Refine3DAlignerGrid::align(EMData * this_img, EMData *to,
02280 const string & cmp_name, const Dict& cmp_params) const
02281 {
02282 if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
02283 throw ImageDimensionException("This aligner only works for 3D images");
02284 }
02285
02286 Transform* t;
02287 if (params.has_key("xform.align3d") ) {
02288
02289
02290
02291
02292 t = params["xform.align3d"];
02293 }else {
02294 t = new Transform();
02295 }
02296
02297 int searchx = 0;
02298 int searchy = 0;
02299 int searchz = 0;
02300 bool dotrans = params.set_default("dotrans",1);
02301 if (params.has_key("search")) {
02302 vector<string> check;
02303 check.push_back("searchx");
02304 check.push_back("searchy");
02305 check.push_back("searchz");
02306 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
02307 if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
02308 }
02309 int search = params["search"];
02310 searchx = search;
02311 searchy = search;
02312 searchz = search;
02313 } else {
02314 searchx = params.set_default("searchx",3);
02315 searchy = params.set_default("searchy",3);
02316 searchz = params.set_default("searchz",3);
02317 }
02318
02319 float delta = params.set_default("delta",1.0f);
02320 float range = params.set_default("range",10.0f);
02321 bool verbose = params.set_default("verbose",false);
02322
02323 bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
02324 EMData * tofft = 0;
02325 if(dotrans || tomography){
02326 tofft = to->do_fft();
02327 }
02328
02329 #ifdef EMAN2_USING_CUDA
02330 if(EMData::usecuda == 1) {
02331 if(!this_img->getcudarodata()) this_img->copy_to_cudaro();
02332 if(!to->getcudarwdata()) to->copy_to_cuda();
02333 if(to->getcudarwdata()){if(tofft) tofft->copy_to_cuda();}
02334 }
02335 #endif
02336
02337 Dict d;
02338 d["type"] = "eman";
02339 Dict best;
02340
02341 best["score"] = 1.0e37;
02342 bool use_cpu = true;
02343 Transform tran = Transform();
02344 Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
02345
02346 for ( float alt = 0; alt < range; alt += delta) {
02347
02348 float saz = 360;
02349 if(alt != 0) saz = delta/sin(alt*M_PI/180.0f);
02350 for ( float az = 0; az < 360; az += saz ){
02351 if (verbose) {
02352 cout << "Trying angle alt " << alt << " az " << az << endl;
02353 }
02354
02355 for( float phi = -range-az; phi < range-az; phi += delta ) {
02356 d["alt"] = alt;
02357 d["phi"] = phi;
02358 d["az"] = az;
02359 Transform tr(d);
02360 tr = tr*(*t);
02361
02362
02363 EMData* transformed;
02364 transformed = this_img->process("xform",Dict("transform",&tr));
02365
02366
02367 float score = 0.0f;
02368 if(dotrans || tomography){
02369 EMData* ccf = transformed->calc_ccf(tofft);
02370 #ifdef EMAN2_USING_CUDA
02371 if(EMData::usecuda == 1){
02372 use_cpu = false;
02373 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
02374 tran.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
02375
02376
02377 tr = tran*tr;
02378 if (tomography) {
02379 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
02380 score = -(data->peak - stats.x)/sqrt(stats.y);
02381 } else {
02382 score = -data->peak;
02383 }
02384 delete data;
02385 }
02386 #endif
02387 if(use_cpu){
02388 if(tomography) ccf->process_inplace("normalize");
02389
02390
02391
02392 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
02393 tran.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
02394 score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
02395 tr = tran*tr;
02396
02397 }
02398 delete ccf; ccf =0;
02399 delete transformed; transformed = 0;
02400 }
02401
02402 if(!tomography){
02403 if(!transformed) transformed = this_img->process("xform",Dict("transform",&tr));
02404 score = c->cmp(to,transformed);
02405 delete transformed; transformed = 0;
02406 }
02407
02408 if(score < float(best["score"])) {
02409
02410 best["score"] = score;
02411 best["xform.align3d"] = &tr;
02412 }
02413 }
02414 }
02415 }
02416
02417 if(tofft) {delete tofft; tofft = 0;}
02418 if (c != 0) delete c;
02419
02420
02421 EMData* best_match = this_img->process("xform",Dict("transform", best["xform.align3d"]));
02422 best_match->set_attr("xform.align3d", best["xform.align3d"]);
02423 best_match->set_attr("score", float(best["score"]));
02424
02425 return best_match;
02426
02427 }
02428
02429 EMData* RT3DGridAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
02430 {
02431
02432 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
02433
02434 Dict t;
02435 Transform* tr = (Transform*) alis[0]["xform.align3d"];
02436 t["transform"] = tr;
02437 EMData* soln = this_img->process("xform",t);
02438 soln->set_attr("xform.align3d",tr);
02439 delete tr; tr = 0;
02440
02441 return soln;
02442
02443 }
02444
02445 vector<Dict> RT3DGridAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const {
02446
02447 if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
02448 throw ImageDimensionException("This aligner only works for 3D images");
02449 }
02450
02451 int searchx = 0;
02452 int searchy = 0;
02453 int searchz = 0;
02454
02455 bool dotrans = params.set_default("dotrans",1);
02456 if (params.has_key("search")) {
02457 vector<string> check;
02458 check.push_back("searchx");
02459 check.push_back("searchy");
02460 check.push_back("searchz");
02461 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
02462 if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
02463 }
02464 int search = params["search"];
02465 searchx = search;
02466 searchy = search;
02467 searchz = search;
02468 } else {
02469 searchx = params.set_default("searchx",3);
02470 searchy = params.set_default("searchy",3);
02471 searchz = params.set_default("searchz",3);
02472 }
02473
02474 Transform* initxform;
02475 if (params.has_key("initxform") ) {
02476
02477
02478
02479
02480 initxform = params["initxform"];
02481 }else {
02482 initxform = new Transform();
02483 }
02484
02485 float lalt = params.set_default("alt0",0.0f);
02486 float laz = params.set_default("az0",0.0f);
02487 float lphi = params.set_default("phi0",0.0f);
02488 float ualt = params.set_default("alt1",180.0f);
02489 float uphi = params.set_default("phi1",360.0f);
02490 float uaz = params.set_default("az1",360.0f);
02491 float dalt = params.set_default("dalt",10.f);
02492 float daz = params.set_default("daz",10.f);
02493 float dphi = params.set_default("dphi",10.f);
02494 bool verbose = params.set_default("verbose",false);
02495
02496
02497 Dict altered_cmp_params(cmp_params);
02498 if (cmp_name == "ccc.tomo") {
02499 altered_cmp_params.set_default("searchx", searchx);
02500 altered_cmp_params.set_default("searchy", searchy);
02501 altered_cmp_params.set_default("searchz", searchz);
02502 altered_cmp_params.set_default("norm", true);
02503 }
02504
02505 vector<Dict> solns;
02506 if (nsoln == 0) return solns;
02507 for (unsigned int i = 0; i < nsoln; ++i ) {
02508 Dict d;
02509 d["score"] = 1.e24;
02510 Transform t;
02511 d["xform.align3d"] = &t;
02512 solns.push_back(d);
02513 }
02514
02515 bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
02516 EMData * tofft = 0;
02517 if(dotrans || tomography){
02518 tofft = to->do_fft();
02519 }
02520
02521 #ifdef EMAN2_USING_CUDA
02522 if(EMData::usecuda == 1) {
02523 if(!this_img->getcudarodata()) this_img->copy_to_cudaro();
02524 if(!to->getcudarwdata()) to->copy_to_cuda();
02525 if(to->getcudarwdata()){if(tofft) tofft->copy_to_cuda();}
02526 }
02527 #endif
02528
02529 Dict d;
02530 d["type"] = "eman";
02531 Transform trans = Transform();
02532 Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
02533 bool use_cpu = true;
02534 for ( float alt = lalt; alt <= ualt; alt += dalt) {
02535
02536
02537 for ( float az = laz; az < uaz; az += daz ){
02538 if (verbose) {
02539 cout << "Trying angle alt " << alt << " az " << az << endl;
02540 }
02541 for( float phi = lphi; phi < uphi; phi += dphi ) {
02542 d["alt"] = alt;
02543 d["phi"] = phi;
02544 d["az"] = az;
02545 Transform t(d);
02546 t = t*(*initxform);
02547 EMData* transformed = this_img->process("xform",Dict("transform",&t));
02548
02549
02550 float best_score = 0.0f;
02551 if(dotrans || tomography){
02552 EMData* ccf = transformed->calc_ccf(tofft);
02553 #ifdef EMAN2_USING_CUDA
02554 if(EMData::usecuda == 1){
02555 use_cpu = false;;
02556 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
02557 trans.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
02558 t = trans*t;
02559 if (tomography) {
02560 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
02561 best_score = -(data->peak - stats.x)/sqrt(stats.y);
02562 } else {
02563 best_score = -data->peak;
02564 }
02565 delete data;
02566 }
02567 #endif
02568 if(use_cpu){
02569 if(tomography) ccf->process_inplace("normalize");
02570 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
02571 trans.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
02572 t = trans*t;
02573 best_score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
02574 }
02575 delete ccf; ccf =0;
02576 delete transformed; transformed = 0;
02577 }
02578
02579 if(!tomography){
02580 if(!transformed) transformed = this_img->process("xform",Dict("transform",&t));
02581 best_score = c->cmp(to,transformed);
02582 delete transformed; transformed = 0;
02583 }
02584
02585 unsigned int j = 0;
02586 for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
02587 if ( (float)(*it)["score"] > best_score ) {
02588 vector<Dict>::reverse_iterator rit = solns.rbegin();
02589 copy(rit+1,solns.rend()-j,rit);
02590 Dict& d = (*it);
02591 d["score"] = best_score;
02592 d["xform.align3d"] = &t;
02593 break;
02594 }
02595 }
02596 }
02597 }
02598 }
02599
02600 if(tofft) {delete tofft; tofft = 0;}
02601 if (c != 0) delete c;
02602
02603 return solns;
02604
02605 }
02606
02607 EMData* RT3DSphereAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
02608 {
02609
02610 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
02611
02612 Dict t;
02613 Transform* tr = (Transform*) alis[0]["xform.align3d"];
02614 t["transform"] = tr;
02615 EMData* soln = this_img->process("xform",t);
02616 soln->set_attr("xform.align3d",tr);
02617 delete tr; tr = 0;
02618
02619 return soln;
02620
02621 }
02622
02623 vector<Dict> RT3DSphereAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const {
02624
02625 if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
02626 throw ImageDimensionException("This aligner only works for 3D images");
02627 }
02628
02629 int searchx = 0;
02630 int searchy = 0;
02631 int searchz = 0;
02632
02633 bool dotrans = params.set_default("dotrans",1);
02634 if (params.has_key("search")) {
02635 vector<string> check;
02636 check.push_back("searchx");
02637 check.push_back("searchy");
02638 check.push_back("searchz");
02639 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
02640 if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
02641 }
02642 int search = params["search"];
02643 searchx = search;
02644 searchy = search;
02645 searchz = search;
02646 } else {
02647 searchx = params.set_default("searchx",3);
02648 searchy = params.set_default("searchy",3);
02649 searchz = params.set_default("searchz",3);
02650 }
02651
02652 Transform* initxform;
02653 if (params.has_key("initxform") ) {
02654
02655
02656
02657
02658 initxform = params["initxform"];
02659 }else {
02660 initxform = new Transform();
02661 }
02662
02663 float lphi = params.set_default("phi0",0.0f);
02664 float uphi = params.set_default("phi1",360.0f);
02665 float dphi = params.set_default("dphi",10.f);
02666 float threshold = params.set_default("threshold",0.f);
02667 if (threshold < 0.0f) throw InvalidParameterException("The threshold parameter must be greater than or equal to zero");
02668 bool verbose = params.set_default("verbose",false);
02669
02670
02671 Dict altered_cmp_params(cmp_params);
02672 if (cmp_name == "ccc.tomo") {
02673 altered_cmp_params.set_default("searchx", searchx);
02674 altered_cmp_params.set_default("searchy", searchy);
02675 altered_cmp_params.set_default("searchz", searchz);
02676 altered_cmp_params.set_default("norm", true);
02677 }
02678
02679 vector<Dict> solns;
02680 if (nsoln == 0) return solns;
02681 for (unsigned int i = 0; i < nsoln; ++i ) {
02682 Dict d;
02683 d["score"] = 1.e24;
02684 Transform t;
02685 d["xform.align3d"] = &t;
02686 solns.push_back(d);
02687 }
02688
02689 Dict d;
02690 d["inc_mirror"] = true;
02691 if ( params.has_key("delta") && params.has_key("n") ) {
02692 throw InvalidParameterException("The delta and n parameters are mutually exclusive in the RT3DSphereAligner aligner");
02693 } else if ( params.has_key("n") ) {
02694 d["n"] = params["n"];
02695 } else {
02696
02697 d["delta"] = params.set_default("delta",10.f);
02698 }
02699
02700 if ((string)params.set_default("orientgen","eman")=="eman") d["perturb"]=0;
02701 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
02702 vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
02703
02704 bool tomography = (cmp_name == "ccc.tomo") ? 1 : 0;
02705
02706
02707 EMData * this_imgfft = 0;
02708 if(dotrans || tomography){
02709 this_imgfft = this_img->do_fft();
02710 }
02711
02712 #ifdef EMAN2_USING_CUDA
02713 if(EMData::usecuda == 1) {
02714 cout << "Using CUDA for 3D alignment" << endl;
02715 if(!to->getcudarodata()) to->copy_to_cudaro();
02716 if(!this_img->getcudarwdata()) this_img->copy_to_cuda();
02717 if(this_imgfft) this_imgfft->copy_to_cuda();
02718 }
02719 #endif
02720
02721 Transform trans = Transform();
02722 Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
02723
02724 bool use_cpu = true;
02725 for(vector<Transform>::const_iterator trans_it = transforms.begin(); trans_it != transforms.end(); trans_it++) {
02726 Dict params = trans_it->get_params("eman");
02727
02728 if (verbose) {
02729 float alt = params["alt"];
02730 float az = params["az"];
02731 cout << "Trying angle alt: " << alt << " az: " << az << endl;
02732 }
02733
02734 for( float phi = lphi; phi < uphi; phi += dphi ) {
02735 params["phi"] = phi;
02736 Transform t(params);
02737 t = t*(*initxform);
02738
02739 EMData* transformed;
02740 transformed = to->process("xform",Dict("transform",&t));
02741
02742
02743 float best_score = 0.0f;
02744
02745 if(dotrans || tomography){
02746 EMData* ccf = transformed->calc_ccf(this_imgfft);
02747 #ifdef EMAN2_USING_CUDA
02748 if(EMData::usecuda == 1){
02749
02750 use_cpu = false;
02751 CudaPeakInfo* data = calc_max_location_wrap_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize(), searchx, searchy, searchz);
02752 trans.set_trans((float)-data->px, (float)-data->py, (float)-data->pz);
02753 t = trans*t;
02754 if (tomography) {
02755 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
02756 best_score = -(data->peak - stats.x)/sqrt(stats.y);
02757 } else {
02758 best_score = -data->peak;
02759 }
02760 delete data;
02761 }
02762 #endif
02763 if(use_cpu){
02764
02765 if(tomography) ccf->process_inplace("normalize");
02766 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
02767 trans.set_trans((float)-point[0], (float)-point[1], (float)-point[2]);
02768 t = trans*t;
02769 best_score = -ccf->get_value_at_wrap(point[0], point[1], point[2]);
02770 }
02771 delete ccf; ccf =0;
02772 delete transformed; transformed = 0;
02773 }
02774
02775 if(!tomography){
02776 if(!transformed) transformed = to->process("xform",Dict("transform",&t));
02777 best_score = c->cmp(this_img,transformed);
02778 delete transformed; transformed = 0;
02779 }
02780
02781 unsigned int j = 0;
02782
02783 for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
02784 if ( (float)(*it)["score"] > best_score ) {
02785 vector<Dict>::reverse_iterator rit = solns.rbegin();
02786 copy(rit+1,solns.rend()-j,rit);
02787 Dict& d = (*it);
02788 d["score"] = best_score;
02789 t.invert();
02790 d["xform.align3d"] = &t;
02791 break;
02792 }
02793 }
02794
02795 }
02796 }
02797
02798 if(this_imgfft) {delete this_imgfft; this_imgfft = 0;}
02799 if(sym!=0) delete sym;
02800 if (c != 0) delete c;
02801
02802 return solns;
02803
02804 }
02805
02806
02807 EMData* RT3DSymmetryAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
02808 {
02809
02810 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
02811
02812 Transform* tr = (Transform*) alis[0]["xform.align3d"];
02813 EMData* soln = this_img->process("xform",Dict("transform",tr));
02814 soln->set_attr("xform.align3d",tr);
02815 delete tr; tr = 0;
02816
02817 return soln;
02818
02819 }
02820
02821 vector<Dict> RT3DSymmetryAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const
02822 {
02823
02824 bool verbose = params.set_default("verbose",false);
02825 Transform* ixform;
02826 if (params.has_key("transform") ) {
02827 ixform = params["transform"];
02828 }else{
02829 ixform = new Transform();
02830 }
02831
02832
02833 vector<Dict> solns;
02834 if (nsoln == 0) return solns;
02835 for (unsigned int i = 0; i < nsoln; ++i ) {
02836 Dict d;
02837 d["score"] = 1.e24;
02838 Transform t;
02839 d["xform.align3d"] = &t;
02840 solns.push_back(d);
02841 }
02842
02843 #ifdef EMAN2_USING_CUDA
02844 if(EMData::usecuda == 1) {
02845 cout << "Using CUDA for 3D sym alignment" << endl;
02846 if(!this_img->getcudarwdata()) this_img->copy_to_cudaro();
02847 if(!to->getcudarwdata()) to->copy_to_cuda();
02848 }
02849 #endif
02850
02851
02852 vector<Transform> syms = Symmetry3D::get_symmetries((string)params.set_default("sym","icos"));
02853 Cmp* c = Factory <Cmp>::get(cmp_name, cmp_params);
02854
02855 float score = 0.0f;
02856 for ( vector<Transform>::const_iterator symit = syms.begin(); symit != syms.end(); ++symit ) {
02857
02858 Transform sympos = (*symit)*(*ixform);
02859 EMData* transformed = this_img->process("xform",Dict("transform", &sympos));
02860 score = c->cmp(transformed,to);
02861 delete transformed; transformed = 0;
02862
02863 if (verbose) {
02864 Dict rots = sympos.get_rotation("eman");
02865 cout <<"Score is: " << score << " az " << float(rots["az"]) << " alt " << float(rots["alt"]) << " phi " << float(rots["phi"]) << endl;
02866 }
02867
02868 unsigned int j = 0;
02869 for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
02870 if ( (float)(*it)["score"] > score ) {
02871 vector<Dict>::reverse_iterator rit = solns.rbegin();
02872 copy(rit+1,solns.rend()-j,rit);
02873 Dict& d = (*it);
02874 d["score"] = score;
02875 d["xform.align3d"] = &sympos;
02876 break;
02877 }
02878 }
02879 }
02880
02881 if (c != 0) delete c;
02882
02883 return solns;
02884 }
02885
02886 namespace {
02887 float 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)
02888 {
02889 int size=rsize;
02890 float dx,dy;
02891 int bw=size/2;
02892 int MAXR=this_img->get_ysize()/2;
02893
02894 unsigned long tsize=2*size;
02895 unsigned long ind1=0, ind2=0, ind3=0, ind4=0, ind41=0;
02896 unsigned long index0=0;
02897 int i=0, j=0, m=0, n=0, r=0;
02898 int loop_rho=0, rho_best=0;
02899
02900 float* gnr2 = new float[size*2];
02901 float* maxcor = new float[size+1];
02902
02903 int p_max=p_max_input;
02904 float* result = new float[5*(p_max+1)];
02905 float* cr=new float[size*(bw+1)];
02906 float* ci=new float[size*(bw+1)];
02907 EMData *data_in=new EMData;
02908 data_in->set_complex(true);
02909 data_in->set_fftodd(false);
02910 data_in->set_ri(true);
02911 data_in->set_size(size+2,size,1);
02912 float *in=data_in->get_data();
02913
02914 float *self_sampl_fft = selfpcsfft->get_data();
02915
02916 float maxcor_sofar=0.0f;
02917 int p=0;
02918
02919 for(p=0; p<=p_max; ++p){
02920 ind1=p*size*bw;
02921 for (i=0;i<size;++i)
02922 for (j=0;j<bw+1;++j){
02923 cr[i*(bw+1)+j]=0.0;
02924 ci[i*(bw+1)+j]=0.0;
02925 }
02926 for(n=0;n<bw;++n){
02927 ind2=(ind1+n);
02928 index0=n*(bw+1);
02929 for(r=0;r<=MAXR;++r) {
02930 ind3=(ind2+r*bw)*size;
02931 for(m=0;m<size;m++){
02932 ind4=(ind3+m)*2;
02933 ind41=ind4+1;
02934 gnr2[2*m]=frm2dhhat[ind4];
02935 gnr2[2*m+1]=frm2dhhat[ind41];
02936 }
02937 for(m=0;m<bw;++m){
02938 float tempr=self_sampl_fft[2*m+r*(size+2)]*r;
02939 float tempi=self_sampl_fft[2*m+1+r*(size+2)]*r;
02940 float gnr2_r=gnr2[2*m];
02941 float gnr2_i=gnr2[2*m+1];
02942 cr[n*(bw+1)+m]+=gnr2_r*tempr+gnr2_i*tempi;
02943 ci[n*(bw+1)+m]+=gnr2_i*tempr-gnr2_r*tempi;
02944 if(n!=0){
02945 if(m!= 0){
02946 int ssize=tsize-2*m;
02947 int ssize1=ssize+1;
02948 float gnr2_r=gnr2[ssize];
02949 float gnr2_i=gnr2[ssize1];
02950 cr[(size-n)*(bw+1)+m]+=gnr2_r*tempr-gnr2_i*tempi;
02951 ci[(size-n)*(bw+1)+m]-=gnr2_i*tempr+gnr2_r*tempi;
02952 }
02953 else{
02954 cr[(size-n)*(bw+1)+m]+=*(gnr2)*tempr-*(gnr2+1)*tempi;
02955 ci[(size-n)*(bw+1)+m]-=*(gnr2+1)*tempr+*(gnr2)*tempi;
02956 }
02957 }
02958 }
02959 }
02960 }
02961 for (int cii=0; cii<size*(bw+1);++cii){
02962 in[2*cii]=cr[cii];
02963 in[2*cii+1]=ci[cii];
02964
02965 }
02966
02967 EMData *data_out;
02968 data_out=data_in->do_ift();
02969 float *c=data_out->get_data();
02970 float tempr=0.0f, corre_fcs=999.0f;
02971
02972 int n_best=0, m_best=0;
02973 float temp=-100.0f;
02974 for(n=0;n<size;++n){
02975 for(m=0;m<size;++m){
02976 temp=c[n*size+m];
02977 if(temp>tempr) {
02978 tempr=temp;
02979 n_best=n;
02980 m_best=m;
02981 }
02982 }
02983 }
02984 delete data_out;
02985
02986 float corre,Phi2,Phi,Tx,Ty,Vx, Vy;
02987
02988
02989
02990
02991
02992 Phi2=n_best*M_PI/bw;
02993 Phi=m_best*M_PI/bw;
02994 Vx=p*cos(Phi);
02995 Vy=-p*sin(Phi);
02996 Tx=Vx+(floor(com_this_x+0.5f)-floor(com_with_x+0.5f));
02997 Ty=Vy+(floor(com_this_y+0.5f)-floor(com_with_y+0.5f));
02998
02999 dx=-Tx;
03000 dy=-Ty;
03001
03002 EMData *this_tmp=this_img->copy();
03003 this_tmp->rotate(-(Phi2-Phi)*180.0f,0.0f,0.0f);
03004 this_tmp->translate(dx,dy,0.0);
03005
03006 corre=this_tmp->cmp(cmp_name,to,cmp_params);
03007
03008 delete this_tmp;
03009 if(corre<=corre_fcs) {
03010 corre_fcs=corre;
03011 result[0+5*p] = float(p);
03012 result[1+5*p] = corre;
03013 result[2+5*p] = (Phi2-Phi)*180.0f;
03014 result[3+5*p] = Tx;
03015 result[4+5*p] = Ty;
03016 }
03017 maxcor[p]=corre_fcs;
03018 if(corre_fcs<maxcor_sofar) {
03019 maxcor_sofar=corre_fcs;
03020 rho_best=p;
03021 }
03022 if(p>=4){
03023 if(maxcor[p] < maxcor[p-1] && maxcor[p-1] < maxcor[p-2]&& maxcor[p-2] < maxcor[p-3] && maxcor[p-3] < maxcor[p-4]){
03024 loop_rho=1;
03025 break;
03026 }
03027 }
03028 }
03029
03030 if(loop_rho == 1)
03031 p=p+1;
03032 int rb5=5*rho_best;
03033 float fsc = result[1+rb5];
03034 float ang_keep = result[2+rb5];
03035 float Tx = result[3+rb5];
03036 float Ty = result[4+rb5];
03037 delete[] gnr2;
03038 delete[] maxcor;
03039 delete[] result;
03040 delete cr;
03041 cr=0;
03042 delete ci;
03043 ci=0;
03044 delete data_in;
03045 dx = -Tx;
03046 dy = -Ty;
03047 this_img->rotate(-ang_keep,0,0);
03048 this_img->translate(dx,dy,0.0);
03049
03050
03051 Transform tsoln(Dict("type","2d","alpha",ang_keep));
03052 tsoln.set_trans(dx,dy);
03053 this_img->set_attr("xform.align2d",&tsoln);
03054 #ifdef DEBUG
03055 float fsc_best=this_img->cmp(cmp_name,to,cmp_params);
03056 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);
03057 #endif
03058 return fsc;
03059 }
03060 }
03061
03062
03063 EMData *FRM2DAligner::align(EMData * this_img, EMData * to,
03064 const string & cmp_name, const Dict& cmp_params) const
03065 {
03066 if (!this_img) {
03067 return 0;
03068 }
03069 if (to && !EMUtil::is_same_size(this_img, to))
03070 throw ImageDimensionException("Images must be the same size to perform translational alignment");
03071
03072 int nx=this_img->get_xsize();
03073 int ny=this_img->get_ysize();
03074 int size =(int)floor(M_PI*ny/4.0);
03075 size =Util::calc_best_fft_size(size);
03076 int MAXR=ny/2;
03077
03078 EMData *this_temp=this_img->copy();
03079 FloatPoint com_test,com_test1;
03080 com_test=this_temp->calc_center_of_mass();
03081 float com_this_x=com_test[0];
03082 float com_this_y=com_test[1];
03083 delete this_temp;
03084
03085
03086 EMData *that_temp=to->copy();
03087 com_test1=that_temp->calc_center_of_mass();
03088 float com_with_x=com_test1[0];
03089 float com_with_y=com_test1[1];
03090 delete that_temp;
03091
03092 EMData *avg_frm=to->copy();
03093 float dx,dy;
03094
03095
03096
03097 EMData *withpcs=avg_frm->unwrap_largerR(0,MAXR,size,float(MAXR));
03098
03099 EMData *withpcsfft=withpcs->oneDfftPolar(size, float(MAXR), float(MAXR));
03100
03101 float *sampl_fft=withpcsfft->get_data();
03102 delete avg_frm;
03103 delete withpcs;
03104
03105 int bw=size/2;
03106 unsigned long ind1=0, ind2=0, ind3=0, ind4=0, ind41=0;
03107 float pi2=2.0*M_PI, r2;
03108
03109 EMData *data_in=new EMData;
03110 data_in->set_complex(true);
03111 data_in->set_ri(1);
03112 data_in->set_size(2*size,1,1);
03113 float * comp_in=data_in->get_data();
03114
03115 int p_max=3;
03116 float *frm2dhhat=0;
03117
03118 if( (frm2dhhat=(float *)malloc((p_max+1)*(size+2)*bw*size*2* sizeof(float)))==NULL){
03119 cout <<"Error in allocating memory 13. \n";
03120 exit(1);
03121 }
03122
03123 float *sb=0, *cb=0;
03124 if((sb=new float[size])==NULL||(cb=new float[size])==NULL) {
03125 cout <<"can't allocate more memory, terminating. \n";
03126 exit(1);
03127 }
03128 for(int i=0;i<size;++i) {
03129 float beta=i*M_PI/bw;
03130 sb[i]=sin(beta);
03131 cb[i]=cos(beta);
03132 }
03133
03134 for(int p=0; p<=p_max; ++p){
03135 ind1=p*size*bw;
03136 float pp2=(float)(p*p);
03137 for(int n=0;n<bw;++n){
03138 ind2=ind1+n;
03139 for(int r=0;r<=MAXR;++r) {
03140 ind3=(ind2+r*bw)*size;
03141 float rr2=(float)(r*r);
03142 float rp2=(float)(r*p);
03143 for(int i=0;i<size;++i){
03144 r2=std::sqrt((float)(rr2+pp2-2.0*rp2*cb[i]));
03145 int r1=(int)floor(r2+0.5f);
03146 if(r1>MAXR){
03147 comp_in[2*i]=0.0f;
03148 comp_in[2*i+1]=0.0f;
03149 }
03150 else{
03151 float gn_r=sampl_fft[2*n+r1*(size+2)];
03152 float gn_i=sampl_fft[2*n+1+r1*(size+2)];
03153 float sb2, cb2, cn, sn;
03154 if(n!=0){
03155 if(r2 != 0.0){
03156 sb2=r*sb[i]/r2;
03157 cb2=(r*cb[i]-p)/r2;
03158 }
03159 else{
03160 sb2=0.0;
03161 cb2=1.0;
03162 }
03163 if(sb2>1.0) sb2= 1.0f;
03164 if(sb2<-1.0)sb2=-1.0f;
03165 if(cb2>1.0) cb2= 1.0f;
03166 if(cb2<-1.0)cb2=-1.0f;
03167 float beta2=atan2(sb2,cb2);
03168 if(beta2<0.0) beta2+=pi2;
03169 float nb2=n*beta2;
03170 cn=cos(nb2);
03171 sn=sin(nb2);
03172 }
03173 else{
03174 cn=1.0f; sn=0.0f;
03175 }
03176 comp_in[2*i]=cn*gn_r-sn*gn_i;
03177 comp_in[2*i+1]=-(cn*gn_i+sn*gn_r);
03178 }
03179 }
03180 EMData *data_out;
03181 data_out=data_in->do_fft();
03182 float * comp_out=data_out->get_data();
03183 for(int m=0;m<size;m++){
03184 ind4=(ind3+m)*2;
03185 ind41=ind4+1;
03186 frm2dhhat[ind4]=comp_out[2*m];
03187 frm2dhhat[ind41]=comp_out[2*m+1];
03188 }
03189 delete data_out;
03190 }
03191 }
03192 }
03193
03194 delete[] sb;
03195 delete[] cb;
03196 delete data_in;
03197 delete withpcsfft;
03198
03199 float dot_frm0=0.0f, dot_frm1=0.0f;
03200 EMData *da_nFlip=0, *da_yFlip=0, *dr_frm=0;
03201
03202 for (int iFlip=0;iFlip<=1;++iFlip){
03203 if (iFlip==0){dr_frm=this_img->copy(); da_nFlip=this_img->copy();}
03204 else {dr_frm=this_img->copy(); da_yFlip=this_img->copy();}
03205 if(iFlip==1) {com_this_x=nx-com_this_x; }
03206
03207 dx=-(com_this_x-nx/2);
03208 dy=-(com_this_y-ny/2);
03209 dr_frm->translate(dx,dy,0.0);
03210 EMData *selfpcs = dr_frm->unwrap_largerR(0,MAXR,size, (float)MAXR);
03211
03212 EMData *selfpcsfft = selfpcs->oneDfftPolar(size, (float)MAXR, (float)MAXR);
03213 delete selfpcs;
03214 delete dr_frm;
03215 if(iFlip==0)
03216 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);
03217 else
03218 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);
03219 delete selfpcsfft;
03220 }
03221
03222 delete[] frm2dhhat;
03223 if(dot_frm0 <=dot_frm1) {
03224 #ifdef DEBUG
03225 printf("best_corre=%f, no flip\n",dot_frm0);
03226 #endif
03227 delete da_yFlip;
03228 return da_nFlip;
03229 }
03230 else {
03231 #ifdef DEBUG
03232 printf("best_corre=%f, flipped\n",dot_frm1);
03233 #endif
03234 delete da_nFlip;
03235 return da_yFlip;
03236 }
03237 }
03238
03239 #ifdef SPARX_USING_CUDA
03240
03241 CUDA_Aligner::CUDA_Aligner(int id) {
03242 image_stack = NULL;
03243 image_stack_filtered = NULL;
03244 ccf = NULL;
03245 if (id != -1) cudasetup(id);
03246 }
03247
03248 void CUDA_Aligner::finish() {
03249 if (image_stack) free(image_stack);
03250 if (image_stack_filtered) free(image_stack_filtered);
03251 if (ccf) free(ccf);
03252 image_stack = NULL;
03253 image_stack_filtered = NULL;
03254 ccf = NULL;
03255 }
03256
03257 void CUDA_Aligner::setup(int nima, int nx, int ny, int ring_length, int nring, int ou, float step, int kx, int ky, bool ctf) {
03258
03259 NIMA = nima;
03260 NX = nx;
03261 NY = ny;
03262 RING_LENGTH = ring_length;
03263 NRING = nring;
03264 STEP = step;
03265 KX = kx;
03266 KY = ky;
03267 OU = ou;
03268 CTF = ctf;
03269
03270 image_stack = (float *)malloc(NIMA*NX*NY*sizeof(float));
03271 if (CTF == 1) image_stack_filtered = (float *)malloc(NIMA*NX*NY*sizeof(float));
03272 ccf = (float *)malloc(2*(2*KX+1)*(2*KY+1)*NIMA*(RING_LENGTH+2)*sizeof(float));
03273 }
03274
03275 void CUDA_Aligner::insert_image(EMData *image, int num) {
03276
03277 int base_address = num*NX*NY;
03278
03279 for (int y=0; y<NY; y++)
03280 for (int x=0; x<NX; x++)
03281 image_stack[base_address+y*NX+x] = (*image)(x, y);
03282 }
03283
03284 void CUDA_Aligner::filter_stack(vector<float> ctf_params) {
03285
03286 float *params;
03287
03288 params = (float *)malloc(NIMA*6*sizeof(float));
03289
03290 for (int i=0; i<NIMA*6; i++) params[i] = ctf_params[i];
03291
03292 filter_image(image_stack, image_stack_filtered, NIMA, NX, NY, params);
03293
03294 free(params);
03295 }
03296
03297 void CUDA_Aligner::sum_oe(vector<float> ctf_params, vector<float> ali_params, EMData *ave1, EMData *ave2) {
03298
03299 float *ctf_p, *ali_p, *av1, *av2;
03300
03301 ctf_p = (float *)malloc(NIMA*6*sizeof(float));
03302 ali_p = (float *)malloc(NIMA*4*sizeof(float));
03303
03304 if (CTF == 1) {
03305 for (int i=0; i<NIMA*6; i++) ctf_p[i] = ctf_params[i];
03306 }
03307 for (int i=0; i<NIMA*4; i++) ali_p[i] = ali_params[i];
03308
03309 av1 = ave1->get_data();
03310 av2 = ave2->get_data();
03311
03312 rot_filt_sum(image_stack, NIMA, NX, NY, CTF, ctf_p, ali_p, av1, av2);
03313
03314 free(ctf_p);
03315 free(ali_p);
03316 }
03317
03318 vector<float> CUDA_Aligner::alignment_2d(EMData *ref_image_em, vector<float> sx_list, vector<float> sy_list, int silent) {
03319
03320 float *ref_image, max_ccf;
03321 int base_address, ccf_offset;
03322 float ts, tm;
03323 float ang, sx = 0, sy = 0, mirror, co, so, sxs, sys;
03324 float *sx2, *sy2;
03325 vector<float> align_result;
03326
03327 sx2 = (float *)malloc(NIMA*sizeof(float));
03328 sy2 = (float *)malloc(NIMA*sizeof(float));
03329
03330 ref_image = ref_image_em->get_data();
03331
03332 for (int i=0; i<NIMA; i++) {
03333 sx2[i] = sx_list[i];
03334 sy2[i] = sy_list[i];
03335 }
03336
03337 if (CTF == 1) {
03338 calculate_ccf(image_stack_filtered, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, silent);
03339 } else {
03340 calculate_ccf(image_stack, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, silent);
03341 }
03342
03343 ccf_offset = NIMA*(RING_LENGTH+2)*(2*KX+1)*(2*KY+1);
03344
03345 for (int im=0; im<NIMA; im++) {
03346 max_ccf = -1.0e22;
03347 for (int kx=-KX; kx<=KX; kx++) {
03348 for (int ky=-KY; ky<=KY; ky++) {
03349 base_address = (((ky+KY)*(2*KX+1)+(kx+KX))*NIMA+im)*(RING_LENGTH+2);
03350 for (int l=0; l<RING_LENGTH; l++) {
03351 ts = ccf[base_address+l];
03352 tm = ccf[base_address+l+ccf_offset];
03353 if (ts > max_ccf) {
03354 ang = float(l)/RING_LENGTH*360.0;
03355 sx = -kx*STEP;
03356 sy = -ky*STEP;
03357 mirror = 0;
03358 max_ccf = ts;
03359 }
03360 if (tm > max_ccf) {
03361 ang = float(l)/RING_LENGTH*360.0;
03362 sx = -kx*STEP;
03363 sy = -ky*STEP;
03364 mirror = 1;
03365 max_ccf = tm;
03366 }
03367 }
03368 }
03369 }
03370 co = cos(ang*M_PI/180.0);
03371 so = -sin(ang*M_PI/180.0);
03372 sxs = sx*co - sy*so;
03373 sys = sx*so + sy*co;
03374
03375 align_result.push_back(ang);
03376 align_result.push_back(sxs);
03377 align_result.push_back(sys);
03378 align_result.push_back(mirror);
03379 }
03380
03381 free(sx2);
03382 free(sy2);
03383
03384 return align_result;
03385 }
03386
03387
03388 vector<float> CUDA_Aligner::ali2d_single_iter(EMData *ref_image_em, vector<float> ali_params, float csx, float csy, int silent, float delta) {
03389
03390 float *ref_image, max_ccf;
03391 int base_address, ccf_offset;
03392 float ts, tm;
03393 float ang = 0.0, sx = 0.0, sy = 0.0, co, so, sxs, sys;
03394 int mirror;
03395 float *sx2, *sy2;
03396 vector<float> align_result;
03397
03398 sx2 = (float *)malloc(NIMA*sizeof(float));
03399 sy2 = (float *)malloc(NIMA*sizeof(float));
03400
03401 ref_image = ref_image_em->get_data();
03402
03403 for (int i=0; i<NIMA; i++) {
03404 ang = ali_params[i*4]/180.0*M_PI;
03405 sx = (ali_params[i*4+3] < 0.5)?(ali_params[i*4+1]-csx):(ali_params[i*4+1]+csx);
03406 sy = ali_params[i*4+2]-csy;
03407 co = cos(ang);
03408 so = sin(ang);
03409 sx2[i] = -(sx*co-sy*so);
03410 sy2[i] = -(sx*so+sy*co);
03411 }
03412
03413 if (CTF == 1) {
03414 calculate_ccf(image_stack_filtered, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, silent);
03415 } else {
03416 calculate_ccf(image_stack, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, silent);
03417 }
03418
03419 ccf_offset = NIMA*(RING_LENGTH+2)*(2*KX+1)*(2*KY+1);
03420
03421 float sx_sum = 0.0f;
03422 float sy_sum = 0.0f;
03423
03424 int dl;
03425 dl = static_cast<int>(delta/360.0*RING_LENGTH);
03426 if (dl<1) { dl = 1; }
03427
03428 for (int im=0; im<NIMA; im++) {
03429 max_ccf = -1.0e22;
03430 for (int kx=-KX; kx<=KX; kx++) {
03431 for (int ky=-KY; ky<=KY; ky++) {
03432 base_address = (((ky+KY)*(2*KX+1)+(kx+KX))*NIMA+im)*(RING_LENGTH+2);
03433 for (int l=0; l<RING_LENGTH; l+=dl) {
03434 ts = ccf[base_address+l];
03435 tm = ccf[base_address+l+ccf_offset];
03436 if (ts > max_ccf) {
03437 ang = float(l)/RING_LENGTH*360.0;
03438 sx = -kx*STEP;
03439 sy = -ky*STEP;
03440 mirror = 0;
03441 max_ccf = ts;
03442 }
03443 if (tm > max_ccf) {
03444 ang = float(l)/RING_LENGTH*360.0;
03445 sx = -kx*STEP;
03446 sy = -ky*STEP;
03447 mirror = 1;
03448 max_ccf = tm;
03449 }
03450 }
03451 }
03452 }
03453 co = cos(ang*M_PI/180.0);
03454 so = -sin(ang*M_PI/180.0);
03455
03456 sxs = (sx-sx2[im])*co-(sy-sy2[im])*so;
03457 sys = (sx-sx2[im])*so+(sy-sy2[im])*co;
03458
03459
03460
03461 align_result.push_back(ang);
03462 align_result.push_back(sxs);
03463 align_result.push_back(sys);
03464 align_result.push_back(mirror);
03465
03466 if (mirror == 0) { sx_sum += sxs; } else { sx_sum -= sxs; }
03467 sy_sum += sys;
03468 }
03469
03470 align_result.push_back(sx_sum);
03471 align_result.push_back(sy_sum);
03472
03473 free(sx2);
03474 free(sy2);
03475
03476 return align_result;
03477 }
03478
03479
03480 CUDA_multiref_aligner::CUDA_multiref_aligner(int id) {
03481 image_stack = NULL;
03482 ref_image_stack = NULL;
03483 ref_image_stack_filtered = NULL;
03484 ccf = NULL;
03485 ctf_params = NULL;
03486 ali_params = NULL;
03487 cudasetup(id);
03488 }
03489
03490
03491 void CUDA_multiref_aligner::finish() {
03492 if (image_stack) free(image_stack);
03493 if (ref_image_stack) free(ref_image_stack);
03494 if (ref_image_stack_filtered) free(ref_image_stack_filtered);
03495 if (ccf) free(ccf);
03496 if (ctf_params) free(ctf_params);
03497 if (ali_params) free(ali_params);
03498 image_stack = NULL;
03499 ref_image_stack = NULL;
03500 ref_image_stack_filtered = NULL;
03501 ccf = NULL;
03502 ctf_params = NULL;
03503 ali_params = NULL;
03504 }
03505
03506 void CUDA_multiref_aligner::setup(int nima, int nref, int nx, int ny, int ring_length, int nring, int ou, float step, int kx, int ky, bool ctf) {
03507
03508 NIMA = nima;
03509 NREF = nref;
03510 NX = nx;
03511 NY = ny;
03512 RING_LENGTH = ring_length;
03513 NRING = nring;
03514 STEP = step;
03515 KX = kx;
03516 KY = ky;
03517 OU = ou;
03518 CTF = ctf;
03519
03520
03521 MAX_IMAGE_BATCH = 10;
03522
03523 image_stack = (float *)malloc(NIMA*NX*NY*sizeof(float));
03524 ref_image_stack = (float *)malloc(NREF*NX*NY*sizeof(float));
03525 if (CTF == 1) ref_image_stack_filtered = (float *)malloc(NREF*NX*NY*sizeof(float));
03526 ccf = (float *)malloc(2*(2*KX+1)*(2*KY+1)*NREF*(RING_LENGTH+2)*MAX_IMAGE_BATCH*sizeof(float));
03527 }
03528
03529 void CUDA_multiref_aligner::setup_params(vector<float> all_ali_params, vector<float> all_ctf_params) {
03530
03531 ali_params = (float *)malloc(NIMA*4*sizeof(float));
03532 for (int i=0; i<NIMA*4; i++) ali_params[i] = all_ali_params[i];
03533 if (CTF == 1) {
03534 ctf_params = (float *)malloc(NIMA*6*sizeof(float));
03535 for (int i=0; i<NIMA*6; i++) ctf_params[i] = all_ctf_params[i];
03536 }
03537 }
03538
03539 void CUDA_multiref_aligner::insert_image(EMData *image, int num) {
03540
03541 int base_address = num*NX*NY;
03542
03543 for (int y=0; y<NY; y++)
03544 for (int x=0; x<NX; x++)
03545 image_stack[base_address+y*NX+x] = (*image)(x, y);
03546 }
03547
03548 void CUDA_multiref_aligner::insert_ref_image(EMData *image, int num) {
03549
03550 int base_address = num*NX*NY;
03551
03552 for (int y=0; y<NY; y++)
03553 for (int x=0; x<NX; x++)
03554 ref_image_stack[base_address+y*NX+x] = (*image)(x, y);
03555 }
03556
03557 vector<float> CUDA_multiref_aligner::multiref_ali2d(int silent) {
03558
03559 float *ctf_params_ref = (float *)malloc(NREF*6*sizeof(float));
03560 float *sx2 = (float *)malloc(NIMA*sizeof(float));
03561 float *sy2 = (float *)malloc(NIMA*sizeof(float));
03562 vector<float> align_results;
03563 int ccf_offset = NREF*(RING_LENGTH+2)*(2*KX+1)*(2*KY+1);
03564
03565 vector<int> batch_size;
03566 vector<int> batch_begin;
03567
03568 if (CTF == 1) {
03569 float previous_defocus = ctf_params[0];
03570 int current_size = 1;
03571 for (int i=1; i<NIMA; i++) {
03572 if (ctf_params[i*6] != previous_defocus || current_size >= MAX_IMAGE_BATCH) {
03573 batch_size.push_back(current_size);
03574 current_size = 1;
03575 previous_defocus = ctf_params[i*6];
03576 } else current_size++;
03577 }
03578 batch_size.push_back(current_size);
03579 } else {
03580 batch_size.resize(NIMA/MAX_IMAGE_BATCH, MAX_IMAGE_BATCH);
03581 if (NIMA%MAX_IMAGE_BATCH != 0) batch_size.push_back(NIMA%MAX_IMAGE_BATCH);
03582 }
03583 int num_batch = batch_size.size();
03584 batch_begin.resize(num_batch, 0);
03585 for (int i=1; i<num_batch; i++) batch_begin[i] = batch_size[i-1]+batch_begin[i-1];
03586 assert(batch_begin[num_batch-1]+batch_size[num_batch-1] == NIMA-1);
03587
03588 for (int i=0; i<NIMA; i++) {
03589 float ang = ali_params[i*4]/180.0*M_PI;
03590 float sx = ali_params[i*4+1];
03591 float sy = ali_params[i*4+2];
03592 float co = cos(ang);
03593 float so = sin(ang);
03594 sx2[i] = -(sx*co-sy*so);
03595 sy2[i] = -(sx*so+sy*co);
03596 }
03597
03598 for (int i=0; i<num_batch; i++) {
03599 if (CTF == 1) {
03600 for (int p=0; p<NREF; p++)
03601 for (int q=0; q<6; q++)
03602 ctf_params_ref[p*6+q] = ctf_params[batch_begin[i]*6+q];
03603 filter_image(ref_image_stack, ref_image_stack_filtered, NREF, NX, NY, ctf_params_ref);
03604 calculate_multiref_ccf(image_stack+batch_begin[i]*NX*NY, ref_image_stack_filtered, ccf, batch_size[i], NREF, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY,
03605 sx2+batch_begin[i], sy2+batch_begin[i], silent);
03606 } else {
03607 calculate_multiref_ccf(image_stack+batch_begin[i]*NX*NY, ref_image_stack, ccf, batch_size[i], NREF, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY,
03608 sx2+batch_begin[i], sy2+batch_begin[i], silent);
03609 }
03610
03611 for (int j=0; j<batch_size[i]; j++) {
03612 for (int im=0; im<NREF; im++) {
03613 float max_ccf = -1.0e22;
03614 float ang = 0.0, sx = 0.0, sy = 0.0;
03615 int mirror = 0;
03616 for (int kx=-KX; kx<=KX; kx++) {
03617 for (int ky=-KY; ky<=KY; ky++) {
03618 int base_address = (((ky+KY)*(2*KX+1)+(kx+KX))*NREF+im)*(RING_LENGTH+2)+ccf_offset*2*j;
03619 for (int l=0; l<RING_LENGTH; l++) {
03620 float ts = ccf[base_address+l];
03621 float tm = ccf[base_address+l+ccf_offset];
03622 if (ts > max_ccf) {
03623 ang = 360.0-float(l)/RING_LENGTH*360.0;
03624 sx = -kx*STEP;
03625 sy = -ky*STEP;
03626 mirror = 0;
03627 max_ccf = ts;
03628 }
03629 if (tm > max_ccf) {
03630 ang = float(l)/RING_LENGTH*360.0;
03631 sx = -kx*STEP;
03632 sy = -ky*STEP;
03633 mirror = 1;
03634 max_ccf = tm;
03635 }
03636 }
03637 }
03638 }
03639 float co = cos(ang*M_PI/180.0);
03640 float so = -sin(ang*M_PI/180.0);
03641
03642 int img_num = batch_begin[i]+j;
03643 float sxs = (sx-sx2[img_num])*co-(sy-sy2[img_num])*so;
03644 float sys = (sx-sx2[img_num])*so+(sy-sy2[img_num])*co;
03645
03646 align_results.push_back(max_ccf);
03647 align_results.push_back(ang);
03648 align_results.push_back(sxs);
03649 align_results.push_back(sys);
03650 align_results.push_back(mirror);
03651 }
03652 }
03653 }
03654
03655 free(ctf_params_ref);
03656 free(sx2);
03657 free(sy2);
03658
03659 return align_results;
03660 }
03661
03662 #endif
03663
03664
03665 void EMAN::dump_aligners()
03666 {
03667 dump_factory < Aligner > ();
03668 }
03669
03670 map<string, vector<string> > EMAN::dump_aligners_list()
03671 {
03672 return dump_factory_list < Aligner > ();
03673 }