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