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