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