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