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
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 <sparx/cuda/cuda_ccf.h>
00047 #endif
00048
00049 #define EMAN2_ALIGNER_DEBUG 0
00050
00051 using namespace EMAN;
00052
00053 const string TranslationalAligner::NAME = "translational";
00054 const string RotationalAligner::NAME = "rotational";
00055 const string RotatePrecenterAligner::NAME = "rotate_precenter";
00056 const string RotateTranslateAligner::NAME = "rotate_translate";
00057 const string RotateTranslateBestAligner::NAME = "rotate_translate_best";
00058 const string RotateFlipAligner::NAME = "rotate_flip";
00059 const string RotateTranslateFlipAligner::NAME = "rotate_translate_flip";
00060 const string RTFExhaustiveAligner::NAME = "rtf_exhaustive";
00061 const string RTFSlowExhaustiveAligner::NAME = "rtf_slow_exhaustive";
00062 const string RefineAligner::NAME = "refine";
00063 const string Refine3DAligner::NAME = "refine.3d";
00064 const string RT3DGridAligner::NAME = "rt.3d.grid";
00065 const string RT3DSphereAligner::NAME = "rt.3d.sphere";
00066
00067 template <> Factory < Aligner >::Factory()
00068 {
00069 force_add<TranslationalAligner>();
00070 force_add<RotationalAligner>();
00071 force_add<RotatePrecenterAligner>();
00072 force_add<RotateTranslateAligner>();
00073 force_add<RotateFlipAligner>();
00074 force_add<RotateTranslateFlipAligner>();
00075 force_add<RTFExhaustiveAligner>();
00076 force_add<RTFSlowExhaustiveAligner>();
00077 force_add<RefineAligner>();
00078 force_add<Refine3DAligner>();
00079 force_add<RT3DGridAligner>();
00080 force_add<RT3DSphereAligner>();
00081
00082 }
00083
00084
00085
00086
00087
00088 EMData *TranslationalAligner::align(EMData * this_img, EMData *to,
00089 const string&, const Dict&) const
00090 {
00091 if (!this_img) {
00092 return 0;
00093 }
00094
00095 if (to && !EMUtil::is_same_size(this_img, to))
00096 throw ImageDimensionException("Images must be the same size to perform translational alignment");
00097
00098 EMData *cf = 0;
00099 int nx = this_img->get_xsize();
00100 int ny = this_img->get_ysize();
00101 int nz = this_img->get_zsize();
00102
00103 int masked = params.set_default("masked",0);
00104 int useflcf = params.set_default("useflcf",0);
00105 bool use_cpu = true;
00106 #ifdef EMAN2_USING_CUDA
00107 if (this_img->gpu_operation_preferred() ) {
00108
00109 use_cpu = false;
00110 cf = this_img->calc_ccf_cuda(to,false,false);
00111 }
00112 #endif // EMAN2_USING_CUDA
00113 if (use_cpu) {
00114 if (useflcf) cf = this_img->calc_flcf(to);
00115 else cf = this_img->calc_ccf(to);
00116 }
00117
00118
00119 if (masked) {
00120 EMData *msk=this_img->process("threshold.notzero");
00121 EMData *sqr=to->process("math.squared");
00122 EMData *cfn=msk->calc_ccf(sqr);
00123 cfn->process_inplace("math.sqrt");
00124 float *d1=cf->get_data();
00125 float *d2=cfn->get_data();
00126 for (int i=0; i<nx*ny*nz; i++) {
00127 if (d2[i]!=0) d1[i]/=d2[i];
00128 }
00129 cf->update();
00130 delete msk;
00131 delete sqr;
00132 delete cfn;
00133 }
00134
00135
00136
00137 int maxshiftx = params.set_default("maxshift",-1);
00138 int maxshifty = params["maxshift"];
00139 int maxshiftz = params["maxshift"];
00140 int nozero = params["nozero"];
00141
00142 if (maxshiftx <= 0) {
00143 maxshiftx = nx / 4;
00144 maxshifty = ny / 4;
00145 maxshiftz = nz / 4;
00146 }
00147
00148 if (maxshiftx > nx / 2 - 1) maxshiftx = nx / 2 - 1;
00149 if (maxshifty > ny / 2 - 1) maxshifty = ny / 2 - 1;
00150 if (maxshiftz > nz / 2 - 1) maxshiftz = nz / 2 - 1;
00151
00152 if (nx == 1) maxshiftx = 0;
00153 if (ny == 1) maxshifty = 0;
00154 if (nz == 1) maxshiftz = 0;
00155
00156
00157 if (nozero) {
00158 cf->zero_corner_circulant(1);
00159 }
00160
00161 IntPoint peak;
00162 #ifdef EMAN2_USING_CUDA
00163 if (!use_cpu) {
00164 EMDataForCuda tmp = cf->get_data_struct_for_cuda();
00165 int* p = calc_max_location_wrap_cuda(&tmp,maxshiftx, maxshifty, maxshiftz);
00166 peak = IntPoint(p[0],p[1],p[2]);
00167 free(p);
00168 }
00169 #endif // EMAN2_USING_CUDA
00170 if (use_cpu) {
00171 peak = cf->calc_max_location_wrap(maxshiftx, maxshifty, maxshiftz);
00172 }
00173 Vec3f cur_trans = Vec3f ( (float)-peak[0], (float)-peak[1], (float)-peak[2]);
00174
00175 if (!to) {
00176 cur_trans /= 2.0f;
00177 int intonly = params.set_default("intonly",false);
00178 if (intonly) {
00179 cur_trans[0] = floor(cur_trans[0] + 0.5f);
00180 cur_trans[1] = floor(cur_trans[1] + 0.5f);
00181 cur_trans[2] = floor(cur_trans[2] + 0.5f);
00182 }
00183 }
00184
00185 if( cf ){
00186 delete cf;
00187 cf = 0;
00188 }
00189 Dict params("trans",static_cast< vector<int> >(cur_trans));
00190 cf=this_img->process("math.translate.int",params);
00191 Transform t;
00192 t.set_trans(cur_trans);
00193 if ( nz != 1 ) {
00194
00195
00196 cf->set_attr("xform.align3d",&t);
00197 } else if ( ny != 1 ) {
00198
00199 cur_trans[2] = 0;
00200 t.set_trans(cur_trans);
00201 cf->set_attr("xform.align2d",&t);
00202 }
00203
00204 return cf;
00205 }
00206
00207 EMData * RotationalAligner::align_180_ambiguous(EMData * this_img, EMData * to, int rfp_mode) {
00208
00209
00210 EMData* this_img_rfp, * to_rfp;
00211 if (rfp_mode == 0) {
00212 this_img_rfp = this_img->make_rotational_footprint_e1();
00213 to_rfp = to->make_rotational_footprint_e1();
00214 } else if (rfp_mode == 1) {
00215 this_img_rfp = this_img->make_rotational_footprint();
00216 to_rfp = to->make_rotational_footprint();
00217 } else if (rfp_mode == 2) {
00218 this_img_rfp = this_img->make_rotational_footprint_cmc();
00219 to_rfp = to->make_rotational_footprint_cmc();
00220 } else {
00221 throw InvalidParameterException("rfp_mode must be 0,1 or 2");
00222 }
00223 int this_img_rfp_nx = this_img_rfp->get_xsize();
00224
00225
00226 EMData *cf = this_img_rfp->calc_ccfx(to_rfp, 0, this_img->get_ysize());
00227
00228
00229 delete this_img_rfp; this_img_rfp = 0;
00230 delete to_rfp; to_rfp = 0;
00231
00232
00233 float *data = cf->get_data();
00234 float peak = 0;
00235 int peak_index = 0;
00236 Util::find_max(data, this_img_rfp_nx, &peak, &peak_index);
00237
00238 if( cf ) {
00239 delete cf;
00240 cf = 0;
00241 }
00242 float rot_angle = (float) (peak_index * 180.0f / this_img_rfp_nx);
00243
00244
00245 Transform tmp(Dict("type","2d","alpha",rot_angle));
00246 cf=this_img->process("xform",Dict("transform",(Transform*)&tmp));
00247
00248
00249
00250 cf->set_attr("xform.align2d",&tmp);
00251 return cf;
00252 }
00253
00254 EMData *RotationalAligner::align(EMData * this_img, EMData *to,
00255 const string& cmp_name, const Dict& cmp_params) const
00256 {
00257 if (!to) throw InvalidParameterException("Can not rotational align - the image to align to is NULL");
00258
00259
00260 int rfp_mode = params.set_default("rfp_mode",0);
00261 EMData* rot_aligned = RotationalAligner::align_180_ambiguous(this_img,to,rfp_mode);
00262 Transform * tmp = rot_aligned->get_attr("xform.align2d");
00263 Dict rot = tmp->get_rotation("2d");
00264 float rotate_angle_solution = rot["alpha"];
00265 delete tmp;
00266
00267 EMData *rot_align_180 = rot_aligned->process("math.rotate.180");
00268
00269
00270 float rot_cmp = rot_aligned->cmp(cmp_name, to, cmp_params);
00271 float rot_180_cmp = rot_align_180->cmp(cmp_name, to, cmp_params);
00272
00273 float score = 0.0;
00274 EMData* result = NULL;
00275 if (rot_cmp < rot_180_cmp){
00276 result = rot_aligned;
00277 score = rot_cmp;
00278 delete rot_align_180; rot_align_180 = 0;
00279 } else {
00280 result = rot_align_180;
00281 score = rot_180_cmp;
00282 delete rot_aligned; rot_aligned = 0;
00283 rotate_angle_solution = rotate_angle_solution-180.0f;
00284 }
00285
00286
00287
00288 Transform tmp2(Dict("type","2d","alpha",rotate_angle_solution));
00289 result->set_attr("xform.align2d",&tmp2);
00290 return result;
00291 }
00292
00293
00294 EMData *RotatePrecenterAligner::align(EMData * this_img, EMData *to,
00295 const string&, const Dict&) const
00296 {
00297 if (!to) {
00298 return 0;
00299 }
00300
00301 int ny = this_img->get_ysize();
00302 int size = Util::calc_best_fft_size((int) (M_PI * ny * 1.5));
00303 EMData *e1 = this_img->unwrap(4, ny * 7 / 16, size, 0, 0, 1);
00304 EMData *e2 = to->unwrap(4, ny * 7 / 16, size, 0, 0, 1);
00305 EMData *cf = e1->calc_ccfx(e2, 0, ny);
00306
00307 float *data = cf->get_data();
00308
00309 float peak = 0;
00310 int peak_index = 0;
00311 Util::find_max(data, size, &peak, &peak_index);
00312 float a = (float) ((1.0f - 1.0f * peak_index / size) * 180. * 2);
00313
00314 Transform rot;
00315 rot.set_rotation(Dict("type","2d","alpha",(float)(a*180./M_PI)));
00316 EMData* rslt = this_img->process("xform",Dict("transform",&rot));
00317 rslt->set_attr("xform.align2d",&rot);
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328 if( e1 )
00329 {
00330 delete e1;
00331 e1 = 0;
00332 }
00333
00334 if( e2 )
00335 {
00336 delete e2;
00337 e2 = 0;
00338 }
00339
00340 if (cf) {
00341 delete cf;
00342 cf = 0;
00343 }
00344
00345 return rslt;
00346 }
00347
00348 EMData *RotateTranslateAligner::align(EMData * this_img, EMData *to,
00349 const string & cmp_name, const Dict& cmp_params) const
00350 {
00351
00352 int rfp_mode = params.set_default("rfp_mode",0);
00353 EMData *rot_align = RotationalAligner::align_180_ambiguous(this_img,to,rfp_mode);
00354 Transform * tmp = rot_align->get_attr("xform.align2d");
00355 Dict rot = tmp->get_rotation("2d");
00356 float rotate_angle_solution = rot["alpha"];
00357 delete tmp;
00358
00359 EMData *rot_align_180 = rot_align->copy();
00360 rot_align_180->process_inplace("math.rotate.180");
00361
00362 Dict trans_params;
00363 trans_params["intonly"] = 0;
00364 trans_params["maxshift"] = params.set_default("maxshift", -1);
00365 trans_params["useflcf"]=params.set_default("useflcf",0);
00366
00367
00368 trans_params["nozero"] = params.set_default("nozero",false);
00369 EMData* rot_trans = rot_align->align("translational", to, trans_params, cmp_name, cmp_params);
00370 if( rot_align ) {
00371 delete rot_align;
00372 rot_align = 0;
00373 }
00374
00375
00376 EMData* rot_180_trans = rot_align_180->align("translational", to, trans_params, cmp_name, cmp_params);
00377 if( rot_align_180 ) {
00378 delete rot_align_180;
00379 rot_align_180 = 0;
00380 }
00381
00382
00383 float cmp1 = rot_trans->cmp(cmp_name, to, cmp_params);
00384 float cmp2 = rot_180_trans->cmp(cmp_name, to, cmp_params);
00385
00386 EMData *result = 0;
00387 if (cmp1 < cmp2) {
00388 if( rot_180_trans ) {
00389 delete rot_180_trans;
00390 rot_180_trans = 0;
00391 }
00392 result = rot_trans;
00393 }
00394 else {
00395 if( rot_trans ) {
00396 delete rot_trans;
00397 rot_trans = 0;
00398 }
00399 result = rot_180_trans;
00400 rotate_angle_solution -= 180.f;
00401 }
00402
00403 Transform* t = result->get_attr("xform.align2d");
00404 t->set_rotation(Dict("type","2d","alpha",rotate_angle_solution));
00405 result->set_attr("xform.align2d",t);
00406 delete t;
00407
00408 return result;
00409 }
00410
00411
00412
00413
00414 EMData* RotateTranslateFlipAligner::align(EMData * this_img, EMData *to,
00415 const string & cmp_name, const Dict& cmp_params) const
00416 {
00417
00418 Dict rt_params("maxshift", params["maxshift"], "rfp_mode", params.set_default("rfp_mode",0),"useflcf",params.set_default("useflcf",0));
00419 EMData *rot_trans_align = this_img->align("rotate_translate",to,rt_params,cmp_name, cmp_params);
00420
00421
00422 EMData *flipped = params.set_default("flip", (EMData *) 0);
00423 bool delete_flag = false;
00424 if (flipped == 0) {
00425 flipped = to->process("xform.flip", Dict("axis", "x"));
00426 delete_flag = true;
00427 }
00428
00429 EMData * rot_trans_align_flip = this_img->align("rotate_translate", flipped, rt_params, cmp_name, cmp_params);
00430 Transform* t = rot_trans_align_flip->get_attr("xform.align2d");
00431 t->set_mirror(true);
00432 rot_trans_align_flip->set_attr("xform.align2d",t);
00433 delete t;
00434
00435
00436 float cmp1 = rot_trans_align->cmp(cmp_name, to, cmp_params);
00437 float cmp2 = rot_trans_align_flip->cmp(cmp_name, flipped, cmp_params);
00438
00439 if (delete_flag){
00440 if(flipped) {
00441 delete flipped;
00442 flipped = 0;
00443 }
00444 }
00445
00446 EMData *result = 0;
00447 if (cmp1 < cmp2 ) {
00448
00449 if( rot_trans_align_flip ) {
00450 delete rot_trans_align_flip;
00451 rot_trans_align_flip = 0;
00452 }
00453 result = rot_trans_align;
00454 }
00455 else {
00456 if( rot_trans_align ) {
00457 delete rot_trans_align;
00458 rot_trans_align = 0;
00459 }
00460 result = rot_trans_align_flip;
00461 result->process_inplace("xform.flip",Dict("axis","x"));
00462 }
00463
00464 return result;
00465 }
00466
00467
00468
00469
00470 EMData *RotateFlipAligner::align(EMData * this_img, EMData *to,
00471 const string& cmp_name, const Dict& cmp_params) const
00472 {
00473 Dict rot_params("rfp_mode",params.set_default("rfp_mode",0));
00474 EMData *r1 = this_img->align("rotational", to, rot_params,cmp_name, cmp_params);
00475
00476
00477 EMData* flipped =to->process("xform.flip", Dict("axis", "x"));
00478 EMData *r2 = this_img->align("rotational", flipped,rot_params, cmp_name, cmp_params);
00479 Transform* t = r2->get_attr("xform.align2d");
00480 t->set_mirror(true);
00481 r2->set_attr("xform.align2d",t);
00482 delete t;
00483
00484 float cmp1 = r1->cmp(cmp_name, to, cmp_params);
00485 float cmp2 = r2->cmp(cmp_name, flipped, cmp_params);
00486
00487 delete flipped; flipped = 0;
00488
00489 EMData *result = 0;
00490
00491 if (cmp1 < cmp2) {
00492 if( r2 )
00493 {
00494 delete r2;
00495 r2 = 0;
00496 }
00497 result = r1;
00498 }
00499 else {
00500 if( r1 )
00501 {
00502 delete r1;
00503 r1 = 0;
00504 }
00505 result = r2;
00506 result->process_inplace("xform.flip",Dict("axis","x"));
00507 }
00508
00509 return result;
00510 }
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527 EMData *RTFExhaustiveAligner::align(EMData * this_img, EMData *to,
00528 const string & cmp_name, const Dict& cmp_params) const
00529 {
00530 EMData *flip = params.set_default("flip", (EMData *) 0);
00531 int maxshift = params.set_default("maxshift", this_img->get_xsize()/8);
00532 if (maxshift < 2) throw InvalidParameterException("maxshift must be greater than or equal to 2");
00533
00534 int ny = this_img->get_ysize();
00535 int xst = (int) floor(2 * M_PI * ny);
00536 xst = Util::calc_best_fft_size(xst);
00537
00538 Dict d("n",2);
00539 EMData *to_shrunk_unwrapped = to->process("math.medianshrink",d);
00540
00541 int to_copy_r2 = to_shrunk_unwrapped->get_ysize() / 2 - 2 - maxshift / 2;
00542 EMData *tmp = to_shrunk_unwrapped->unwrap(4, to_copy_r2, xst / 2, 0, 0, true);
00543 if( to_shrunk_unwrapped )
00544 {
00545 delete to_shrunk_unwrapped;
00546 to_shrunk_unwrapped = 0;
00547 }
00548 to_shrunk_unwrapped = tmp;
00549
00550 EMData *to_shrunk_unwrapped_copy = to_shrunk_unwrapped->copy();
00551 EMData* to_unwrapped = to->unwrap(4, to->get_ysize() / 2 - 2 - maxshift, xst, 0, 0, true);
00552 EMData *to_unwrapped_copy = to_unwrapped->copy();
00553
00554 bool delete_flipped = true;
00555 EMData *flipped = 0;
00556 if (flip) {
00557 delete_flipped = false;
00558 flipped = flip;
00559 }
00560 else {
00561 flipped = to->process("xform.flip", Dict("axis", "x"));
00562 }
00563 EMData *to_shrunk_flipped_unwrapped = flipped->process("math.medianshrink",d);
00564 tmp = to_shrunk_flipped_unwrapped->unwrap(4, to_copy_r2, xst / 2, 0, 0, true);
00565 if( to_shrunk_flipped_unwrapped )
00566 {
00567 delete to_shrunk_flipped_unwrapped;
00568 to_shrunk_flipped_unwrapped = 0;
00569 }
00570 to_shrunk_flipped_unwrapped = tmp;
00571 EMData *to_shrunk_flipped_unwrapped_copy = to_shrunk_flipped_unwrapped->copy();
00572 EMData* to_flip_unwrapped = flipped->unwrap(4, to->get_ysize() / 2 - 2 - maxshift, xst, 0, 0, true);
00573 EMData* to_flip_unwrapped_copy = to_flip_unwrapped->copy();
00574
00575 if (delete_flipped && flipped != 0) {
00576 delete flipped;
00577 flipped = 0;
00578 }
00579
00580 EMData *this_shrunk_2 = this_img->process("math.medianshrink",d);
00581
00582 float bestval = FLT_MAX;
00583 float bestang = 0;
00584 int bestflip = 0;
00585 float bestdx = 0;
00586 float bestdy = 0;
00587
00588 int half_maxshift = maxshift / 2;
00589
00590 int ur2 = this_shrunk_2->get_ysize() / 2 - 2 - half_maxshift;
00591 for (int dy = -half_maxshift; dy <= half_maxshift; dy += 1) {
00592 for (int dx = -half_maxshift; dx <= half_maxshift; dx += 1) {
00593 #ifdef _WIN32
00594 if (_hypot(dx, dy) <= half_maxshift) {
00595 #else
00596 if (hypot(dx, dy) <= half_maxshift) {
00597 #endif
00598 EMData *uw = this_shrunk_2->unwrap(4, ur2, xst / 2, dx, dy, true);
00599 EMData *uwc = uw->copy();
00600 EMData *a = uw->calc_ccfx(to_shrunk_unwrapped);
00601
00602 uwc->rotate_x(a->calc_max_index());
00603 float cm = uwc->cmp(cmp_name, to_shrunk_unwrapped_copy, cmp_params);
00604 if (cm < bestval) {
00605 bestval = cm;
00606 bestang = (float) (2.0 * M_PI * a->calc_max_index() / a->get_xsize());
00607 bestdx = (float)dx;
00608 bestdy = (float)dy;
00609 bestflip = 0;
00610 }
00611
00612
00613 if( a )
00614 {
00615 delete a;
00616 a = 0;
00617 }
00618 if( uw )
00619 {
00620 delete uw;
00621 uw = 0;
00622 }
00623 if( uwc )
00624 {
00625 delete uwc;
00626 uwc = 0;
00627 }
00628 uw = this_shrunk_2->unwrap(4, ur2, xst / 2, dx, dy, true);
00629 uwc = uw->copy();
00630 a = uw->calc_ccfx(to_shrunk_flipped_unwrapped);
00631
00632 uwc->rotate_x(a->calc_max_index());
00633 cm = uwc->cmp(cmp_name, to_shrunk_flipped_unwrapped_copy, cmp_params);
00634 if (cm < bestval) {
00635 bestval = cm;
00636 bestang = (float) (2.0 * M_PI * a->calc_max_index() / a->get_xsize());
00637 bestdx = (float)dx;
00638 bestdy = (float)dy;
00639 bestflip = 1;
00640 }
00641
00642 if( a )
00643 {
00644 delete a;
00645 a = 0;
00646 }
00647
00648 if( uw )
00649 {
00650 delete uw;
00651 uw = 0;
00652 }
00653 if( uwc )
00654 {
00655 delete uwc;
00656 uwc = 0;
00657 }
00658 }
00659 }
00660 }
00661 if( this_shrunk_2 )
00662 {
00663 delete this_shrunk_2;
00664 this_shrunk_2 = 0;
00665 }
00666 if( to_shrunk_unwrapped )
00667 {
00668 delete to_shrunk_unwrapped;
00669 to_shrunk_unwrapped = 0;
00670 }
00671 if( to_shrunk_unwrapped_copy )
00672 {
00673 delete to_shrunk_unwrapped_copy;
00674 to_shrunk_unwrapped_copy = 0;
00675 }
00676 if( to_shrunk_flipped_unwrapped )
00677 {
00678 delete to_shrunk_flipped_unwrapped;
00679 to_shrunk_flipped_unwrapped = 0;
00680 }
00681 if( to_shrunk_flipped_unwrapped_copy )
00682 {
00683 delete to_shrunk_flipped_unwrapped_copy;
00684 to_shrunk_flipped_unwrapped_copy = 0;
00685 }
00686 bestdx *= 2;
00687 bestdy *= 2;
00688 bestval = FLT_MAX;
00689
00690 float bestdx2 = bestdx;
00691 float bestdy2 = bestdy;
00692
00693
00694
00695 for (float dy = bestdy2 - 3; dy <= bestdy2 + 3; dy += 1.0 ) {
00696 for (float dx = bestdx2 - 3; dx <= bestdx2 + 3; dx += 1.0 ) {
00697
00698 #ifdef _WIN32
00699 if (_hypot(dx, dy) <= maxshift) {
00700 #else
00701 if (hypot(dx, dy) <= maxshift) {
00702 #endif
00703 EMData *uw = this_img->unwrap(4, this_img->get_ysize() / 2 - 2 - maxshift, xst, (int)dx, (int)dy, true);
00704 EMData *uwc = uw->copy();
00705 EMData *a = uw->calc_ccfx(to_unwrapped);
00706
00707 uwc->rotate_x(a->calc_max_index());
00708 float cm = uwc->cmp(cmp_name, to_unwrapped_copy, cmp_params);
00709
00710 if (cm < bestval) {
00711 bestval = cm;
00712 bestang = (float)(2.0 * M_PI * a->calc_max_index() / a->get_xsize());
00713 bestdx = dx;
00714 bestdy = dy;
00715 bestflip = 0;
00716 }
00717
00718 if( a )
00719 {
00720 delete a;
00721 a = 0;
00722 }
00723 if( uw )
00724 {
00725 delete uw;
00726 uw = 0;
00727 }
00728 if( uwc )
00729 {
00730 delete uwc;
00731 uwc = 0;
00732 }
00733 uw = this_img->unwrap(4, this_img->get_ysize() / 2 - 2 - maxshift, xst, (int)dx, (int)dy, true);
00734 uwc = uw->copy();
00735 a = uw->calc_ccfx(to_flip_unwrapped);
00736
00737 uwc->rotate_x(a->calc_max_index());
00738 cm = uwc->cmp(cmp_name, to_flip_unwrapped_copy, cmp_params);
00739
00740 if (cm < bestval) {
00741 bestval = cm;
00742 bestang = (float)(2.0 * M_PI * a->calc_max_index() / a->get_xsize());
00743 bestdx = dx;
00744 bestdy = dy;
00745 bestflip = 1;
00746 }
00747
00748 if( a )
00749 {
00750 delete a;
00751 a = 0;
00752 }
00753 if( uw )
00754 {
00755 delete uw;
00756 uw = 0;
00757 }
00758 if( uwc )
00759 {
00760 delete uwc;
00761 uwc = 0;
00762 }
00763 }
00764 }
00765 }
00766 if( to_unwrapped ) {delete to_unwrapped;to_unwrapped = 0;}
00767 if( to_shrunk_unwrapped ) { delete to_shrunk_unwrapped; to_shrunk_unwrapped = 0;}
00768 if (to_unwrapped_copy) { delete to_unwrapped_copy; to_unwrapped_copy = 0; }
00769 if (to_flip_unwrapped) { delete to_flip_unwrapped; to_flip_unwrapped = 0; }
00770 if (to_flip_unwrapped_copy) { delete to_flip_unwrapped_copy; to_flip_unwrapped_copy = 0;}
00771
00772 bestang *= (float)EMConsts::rad2deg;
00773 Transform t(Dict("type","2d","alpha",(float)bestang));
00774 t.set_pre_trans(Vec2f(-bestdx,-bestdy));
00775 if (bestflip) {
00776 t.set_mirror(true);
00777 }
00778
00779 EMData* ret = this_img->process("xform",Dict("transform",&t));
00780 ret->set_attr("xform.align2d",&t);
00781
00782 return ret;
00783 }
00784
00785
00786 EMData *RTFSlowExhaustiveAligner::align(EMData * this_img, EMData *to,
00787 const string & cmp_name, const Dict& cmp_params) const
00788 {
00789
00790 EMData *flip = params.set_default("flip", (EMData *) 0);
00791 int maxshift = params.set_default("maxshift", -1);
00792
00793 EMData *flipped = 0;
00794
00795 bool delete_flipped = true;
00796 if (flip) {
00797 delete_flipped = false;
00798 flipped = flip;
00799 }
00800 else {
00801 flipped = to->process("xform.flip", Dict("axis", "x"));
00802 }
00803
00804 int nx = this_img->get_xsize();
00805
00806 if (maxshift < 0) {
00807 maxshift = nx / 10;
00808 }
00809
00810 float angle_step = params.set_default("angstep", 0.0f);
00811 if ( angle_step == 0 ) angle_step = atan2(2.0f, (float)nx);
00812 else {
00813 angle_step *= (float)EMConsts::deg2rad;
00814 }
00815 float trans_step = params.set_default("transtep",1.0f);
00816
00817 if (trans_step <= 0) throw InvalidParameterException("transstep must be greater than 0");
00818 if (angle_step <= 0) throw InvalidParameterException("angstep must be greater than 0");
00819
00820
00821 Dict shrinkfactor("n",2);
00822 EMData *this_img_shrink = this_img->process("math.medianshrink",shrinkfactor);
00823 EMData *to_shrunk = to->process("math.medianshrink",shrinkfactor);
00824 EMData *flipped_shrunk = flipped->process("math.medianshrink",shrinkfactor);
00825
00826 int bestflip = 0;
00827 float bestdx = 0;
00828 float bestdy = 0;
00829
00830 float bestang = 0;
00831 float bestval = FLT_MAX;
00832
00833 int half_maxshift = maxshift / 2;
00834
00835
00836 for (int dy = -half_maxshift; dy <= half_maxshift; ++dy) {
00837 for (int dx = -half_maxshift; dx <= half_maxshift; ++dx) {
00838 if (hypot(dx, dy) <= maxshift) {
00839 for (float ang = -angle_step * 2.0f; ang <= (float)2 * M_PI; ang += angle_step * 4.0f) {
00840 EMData v(*this_img_shrink);
00841 Transform t(Dict("type","2d","alpha",static_cast<float>(ang*EMConsts::rad2deg)));
00842 t.set_trans((float)dx,(float)dy);
00843 v.transform(t);
00844
00845
00846 float lc = v.cmp(cmp_name, to_shrunk, cmp_params);
00847
00848 if (lc < bestval) {
00849 bestval = lc;
00850 bestang = ang;
00851 bestdx = (float)dx;
00852 bestdy = (float)dy;
00853 bestflip = 0;
00854 }
00855
00856 lc = v.cmp(cmp_name,flipped_shrunk , cmp_params);
00857 if (lc < bestval) {
00858 bestval = lc;
00859 bestang = ang;
00860 bestdx = (float)dx;
00861 bestdy = (float)dy;
00862 bestflip = 1;
00863 }
00864 }
00865 }
00866 }
00867 }
00868
00869 if( to_shrunk )
00870 {
00871 delete to_shrunk;
00872 to_shrunk = 0;
00873 }
00874 if( flipped_shrunk )
00875 {
00876 delete flipped_shrunk;
00877 flipped_shrunk = 0;
00878 }
00879 if( this_img_shrink )
00880 {
00881 delete this_img_shrink;
00882 this_img_shrink = 0;
00883 }
00884
00885 bestdx *= 2;
00886 bestdy *= 2;
00887 bestval = FLT_MAX;
00888
00889 float bestdx2 = bestdx;
00890 float bestdy2 = bestdy;
00891 float bestang2 = bestang;
00892
00893 for (float dy = bestdy2 - 3; dy <= bestdy2 + 3; dy += trans_step) {
00894 for (float dx = bestdx2 - 3; dx <= bestdx2 + 3; dx += trans_step) {
00895 if (hypot(dx, dy) <= maxshift) {
00896 for (float ang = bestang2 - angle_step * 6.0f; ang <= bestang2 + angle_step * 6.0f; ang += angle_step) {
00897 EMData v(*this_img);
00898 Transform t(Dict("type","2d","alpha",static_cast<float>(ang*EMConsts::rad2deg)));
00899 t.set_trans(dx,dy);
00900 v.transform(t);
00901
00902
00903 float lc = v.cmp(cmp_name, to, cmp_params);
00904
00905 if (lc < bestval) {
00906 bestval = lc;
00907 bestang = ang;
00908 bestdx = dx;
00909 bestdy = dy;
00910 bestflip = 0;
00911 }
00912
00913 lc = v.cmp(cmp_name, flipped, cmp_params);
00914
00915 if (lc < bestval) {
00916 bestval = lc;
00917 bestang = ang;
00918 bestdx = dx;
00919 bestdy = dy;
00920 bestflip = 1;
00921 }
00922 }
00923 }
00924 }
00925 }
00926
00927 if (delete_flipped) { delete flipped; flipped = 0; }
00928
00929 bestang *= (float)EMConsts::rad2deg;
00930 Transform t(Dict("type","2d","alpha",(float)bestang));
00931 t.set_trans(bestdx,bestdy);
00932
00933 if (bestflip) {
00934 t.set_mirror(true);
00935 }
00936
00937 EMData* rslt = this_img->process("xform",Dict("transform",&t));
00938 rslt->set_attr("xform.align2d",&t);
00939
00940 return rslt;
00941 }
00942
00943
00944
00945 static double refalifn(const gsl_vector * v, void *params)
00946 {
00947 Dict *dict = (Dict *) params;
00948
00949 double x = gsl_vector_get(v, 0);
00950 double y = gsl_vector_get(v, 1);
00951 double a = gsl_vector_get(v, 2);
00952
00953 EMData *this_img = (*dict)["this"];
00954 EMData *with = (*dict)["with"];
00955 bool mirror = (*dict)["mirror"];
00956
00957
00958
00959
00960
00961
00962 Transform t(Dict("type","2d","alpha",static_cast<float>(a)));
00963
00964
00965
00966 t.set_trans((float)x,(float)y);
00967 t.set_mirror(mirror);
00968 EMData *tmp = this_img->process("xform",Dict("transform",&t));
00969
00970 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
00971 double result = c->cmp(tmp,with);
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00987
00988
00989
00990
00991
00992
00993 if ( tmp != 0 ) delete tmp;
00994
00995 return result;
00996 }
00997
00998 static double refalifnfast(const gsl_vector * v, void *params)
00999 {
01000 Dict *dict = (Dict *) params;
01001 EMData *this_img = (*dict)["this"];
01002 EMData *img_to = (*dict)["with"];
01003 bool mirror = (*dict)["mirror"];
01004
01005 double x = gsl_vector_get(v, 0);
01006 double y = gsl_vector_get(v, 1);
01007 double a = gsl_vector_get(v, 2);
01008
01009 double r = this_img->dot_rotate_translate(img_to, (float)x, (float)y, (float)a, mirror);
01010 int nsec = this_img->get_xsize() * this_img->get_ysize();
01011 double result = 1.0 - r / nsec;
01012
01013
01014 return result;
01015 }
01016
01017
01018 EMData *RefineAligner::align(EMData * this_img, EMData *to,
01019 const string & cmp_name, const Dict& cmp_params) const
01020 {
01021
01022 if (!to) {
01023 return 0;
01024 }
01025
01026 int mode = params.set_default("mode", 0);
01027 float saz = 0.0;
01028 float sdx = 0.0;
01029 float sdy = 0.0;
01030 bool mirror = false;
01031 Transform* t;
01032 if (params.has_key("xform.align2d") ) {
01033 t = params["xform.align2d"];
01034 Dict params = t->get_params("2d");
01035 saz = params["alpha"];
01036 sdx = params["tx"];
01037 sdy = params["ty"];
01038 mirror = params["mirror"];
01039
01040 } else {
01041 t = new Transform();
01042 }
01043
01044 int np = 3;
01045 Dict gsl_params;
01046 gsl_params["this"] = this_img;
01047 gsl_params["with"] = to;
01048 gsl_params["snr"] = params["snr"];
01049 gsl_params["mirror"] = mirror;
01050
01051
01052
01053 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
01054 gsl_vector *ss = gsl_vector_alloc(np);
01055
01056 float stepx = params.set_default("stepx",1.0f);
01057 float stepy = params.set_default("stepy",1.0f);
01058
01059 float stepaz = params.set_default("stepaz",5.0f);
01060
01061 gsl_vector_set(ss, 0, stepx);
01062 gsl_vector_set(ss, 1, stepy);
01063 gsl_vector_set(ss, 2, stepaz);
01064
01065 gsl_vector *x = gsl_vector_alloc(np);
01066 gsl_vector_set(x, 0, sdx);
01067 gsl_vector_set(x, 1, sdy);
01068 gsl_vector_set(x, 2, saz);
01069
01070 Cmp *c = 0;
01071
01072 gsl_multimin_function minex_func;
01073 if (mode == 2) {
01074 minex_func.f = &refalifnfast;
01075 }
01076 else {
01077 c = Factory < Cmp >::get(cmp_name, cmp_params);
01078 gsl_params["cmp"] = (void *) c;
01079 minex_func.f = &refalifn;
01080 }
01081
01082 minex_func.n = np;
01083 minex_func.params = (void *) &gsl_params;
01084
01085 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
01086 gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
01087
01088 int rval = GSL_CONTINUE;
01089 int status = GSL_SUCCESS;
01090 int iter = 1;
01091
01092 float precision = params.set_default("precision",0.04f);
01093 int maxiter = params.set_default("maxiter",28);
01094
01095
01096
01097
01098 while (rval == GSL_CONTINUE && iter < maxiter) {
01099 iter++;
01100 status = gsl_multimin_fminimizer_iterate(s);
01101 if (status) {
01102 break;
01103 }
01104 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
01105 }
01106
01107 int maxshift = params.set_default("maxshift",-1);
01108
01109 if (maxshift <= 0) {
01110 maxshift = this_img->get_xsize() / 4;
01111 }
01112 float fmaxshift = static_cast<float>(maxshift);
01113 EMData *result;
01114 if ( fmaxshift >= fabs((float)gsl_vector_get(s->x, 0)) && fmaxshift >= fabs((float)gsl_vector_get(s->x, 1)) )
01115 {
01116
01117 Transform tsoln(Dict("type","2d","alpha",(float)gsl_vector_get(s->x, 2)));
01118 tsoln.set_mirror(mirror);
01119 tsoln.set_trans((float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1));
01120 result = this_img->process("xform",Dict("transform",&tsoln));
01121 result->set_attr("xform.align2d",&tsoln);
01122 } else {
01123
01124 result = this_img->process("xform",Dict("transform",t));
01125 result->set_attr("xform.align2d",t);
01126 }
01127
01128 delete t;
01129 t = 0;
01130
01131 gsl_vector_free(x);
01132 gsl_vector_free(ss);
01133 gsl_multimin_fminimizer_free(s);
01134
01135 if ( c != 0 ) delete c;
01136 return result;
01137 }
01138
01139 static Transform refalin3d_perturb(const Transform*const t, const float& delta, const float& arc, const float& phi, const float& x, const float& y, const float& z)
01140 {
01141 Dict orig_params = t->get_params("eman");
01142 float orig_phi = orig_params["phi"];
01143 float orig_x = orig_params["tx"];
01144 float orig_y = orig_params["ty"];
01145 float orig_z = orig_params["tz"];
01146 orig_params["phi"] = 0;
01147 orig_params["tx"] = 0;
01148 orig_params["ty"] = 0;
01149 orig_params["tz"] = 0;
01150 Transform t_no_phi(orig_params);
01151
01152 Vec3f zz(0,0,1);
01153
01154 Vec3f vv = t_no_phi.transpose()*zz;
01155 Vec3f normal = Vec3f(-vv[2],0,-vv[0]);
01156
01157 normal.normalize();
01158
01159 Dict d;
01160 d["type"] = "spin";
01161 d["Omega"] = arc;
01162 d["n1"] = vv[0];
01163 d["n2"] = vv[1];
01164 d["n3"] = vv[2];
01165
01166 Transform q(d);
01167
01168 Vec3f rot_vec = q*normal;
01169 rot_vec.normalize();
01170
01171 Dict e;
01172 e["type"] = "spin";
01173 e["Omega"] = delta;
01174 e["n1"] = rot_vec[0];
01175 e["n2"] = rot_vec[1];
01176 e["n3"] = rot_vec[2];
01177
01178 Transform perturb(e);
01179
01180 Dict g;
01181 g["type"] = "eman";
01182 g["alt"] = 0;
01183 g["az"] = 0;
01184 g["phi"] = 0*phi+orig_phi;
01185
01186 Transform phi_rot(g);
01187 Transform soln = t_no_phi*perturb*phi_rot;
01188 soln.set_trans(x+orig_x,y+orig_y,z+orig_z);
01189
01190 Dict params = soln.get_params("eman");
01191 return soln;
01192 }
01193
01194 static double refalifn3d(const gsl_vector * v, void *params)
01195 {
01196 Dict *dict = (Dict *) params;
01197 double x = gsl_vector_get(v, 0);
01198 double y = gsl_vector_get(v, 1);
01199 double z = gsl_vector_get(v, 2);
01200 double arc = gsl_vector_get(v, 3);
01201 double delta = gsl_vector_get(v, 4);
01202 double phi = gsl_vector_get(v, 5);
01203 EMData *this_img = (*dict)["this"];
01204 EMData *with = (*dict)["with"];
01205
01206
01207 Transform* t = (*dict)["transform"];
01208
01209 Transform soln = refalin3d_perturb(t,(float)delta,(float)arc,(float)phi,(float)x,(float)y,(float)z);
01210
01211 EMData *tmp = this_img->process("xform",Dict("transform",&soln));
01212 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
01213 double result = c->cmp(tmp,with);
01214 if ( tmp != 0 ) delete tmp;
01215 delete t; t = 0;
01216
01217 return result;
01218 }
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250 EMData* Refine3DAligner::align(EMData * this_img, EMData *to,
01251 const string & cmp_name, const Dict& cmp_params) const
01252 {
01253
01254 if (!to || !this_img) throw NullPointerException("Input image is null");
01255
01256 if (to->get_ndim() != 3 || this_img->get_ndim() != 3) throw ImageDimensionException("The Refine3D aligner only works for 3D images");
01257
01258 float saz = 0.0;
01259 float sphi = 0.0;
01260 float salt = 0.0;
01261 float sdx = 0.0;
01262 float sdy = 0.0;
01263 float sdz = 0.0;
01264 bool mirror = false;
01265 Transform* t;
01266 if (params.has_key("xform.align3d") ) {
01267
01268
01269
01270
01271 t = params["xform.align3d"];
01272 }else {
01273 t = new Transform();
01274 }
01275
01276 int np = 6;
01277 Dict gsl_params;
01278 gsl_params["this"] = this_img;
01279 gsl_params["with"] = to;
01280 gsl_params["snr"] = params["snr"];
01281 gsl_params["mirror"] = mirror;
01282 gsl_params["transform"] = t;
01283
01284 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
01285 gsl_vector *ss = gsl_vector_alloc(np);
01286
01287 float stepx = params.set_default("stepx",1.0f);
01288 float stepy = params.set_default("stepy",1.0f);
01289 float stepz = params.set_default("stepz",1.0f);
01290
01291 float half_circle_step = 180.0f;
01292 float stepphi = params.set_default("stephi",5.0f);
01293 float stepdelta = params.set_default("stepdelta",5.0f);
01294
01295 gsl_vector_set(ss, 0, stepx);
01296 gsl_vector_set(ss, 1, stepy);
01297 gsl_vector_set(ss, 2, stepz);
01298 gsl_vector_set(ss, 3, half_circle_step);
01299 gsl_vector_set(ss, 4, stepdelta);
01300 gsl_vector_set(ss, 5, stepphi);
01301
01302 gsl_vector *x = gsl_vector_alloc(np);
01303 gsl_vector_set(x, 0, sdx);
01304 gsl_vector_set(x, 1, sdy);
01305 gsl_vector_set(x, 2, sdz);
01306 gsl_vector_set(x, 3, saz);
01307 gsl_vector_set(x, 4, salt);
01308 gsl_vector_set(x, 5, sphi);
01309
01310 gsl_multimin_function minex_func;
01311 Cmp *c = Factory < Cmp >::get(cmp_name, cmp_params);
01312 gsl_params["cmp"] = (void *) c;
01313 minex_func.f = &refalifn3d;
01314
01315 minex_func.n = np;
01316 minex_func.params = (void *) &gsl_params;
01317
01318 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
01319 gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
01320
01321 int rval = GSL_CONTINUE;
01322 int status = GSL_SUCCESS;
01323 int iter = 1;
01324
01325 float precision = params.set_default("precision",0.04f);
01326 int maxiter = params.set_default("maxiter",60);
01327 while (rval == GSL_CONTINUE && iter < maxiter) {
01328 iter++;
01329 status = gsl_multimin_fminimizer_iterate(s);
01330 if (status) {
01331 break;
01332 }
01333 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
01334 }
01335
01336 int maxshift = params.set_default("maxshift",-1);
01337
01338 if (maxshift <= 0) {
01339 maxshift = this_img->get_xsize() / 4;
01340 }
01341 float fmaxshift = static_cast<float>(maxshift);
01342 EMData *result;
01343 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))
01344 {
01345 float x = (float)gsl_vector_get(s->x, 0);
01346 float y = (float)gsl_vector_get(s->x, 1);
01347 float z = (float)gsl_vector_get(s->x, 2);
01348
01349 float arc = (float)gsl_vector_get(s->x, 3);
01350 float delta = (float)gsl_vector_get(s->x, 4);
01351 float phi = (float)gsl_vector_get(s->x, 5);
01352
01353 Transform tsoln = refalin3d_perturb(t,delta,arc,phi,x,y,z);
01354
01355 result = this_img->process("xform",Dict("transform",&tsoln));
01356 result->set_attr("xform.align3d",&tsoln);
01357
01358 } else {
01359 result = this_img->process("xform",Dict("transform",t));
01360 result->set_attr("xform.align3d",t);
01361 }
01362
01363 delete t;
01364 t = 0;
01365
01366 gsl_vector_free(x);
01367 gsl_vector_free(ss);
01368 gsl_multimin_fminimizer_free(s);
01369
01370 if ( c != 0 ) delete c;
01371 return result;
01372 }
01373
01374 EMData* RT3DGridAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
01375 {
01376
01377 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
01378
01379 Dict t;
01380 Transform* tr = (Transform*) alis[0]["xform.align3d"];
01381 t["transform"] = tr;
01382 EMData* soln = this_img->process("xform",t);
01383 soln->set_attr("xform.align3d",tr);
01384 delete tr; tr = 0;
01385
01386 return soln;
01387
01388 }
01389
01390 vector<Dict> RT3DGridAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const {
01391
01392 if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
01393 throw ImageDimensionException("This aligner only works for 3D images");
01394 }
01395
01396 int searchx = 0;
01397 int searchy = 0;
01398 int searchz = 0;
01399
01400 if (params.has_key("search")) {
01401 vector<string> check;
01402 check.push_back("searchx");
01403 check.push_back("searchy");
01404 check.push_back("searchz");
01405 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
01406 if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
01407 }
01408 int search = params["search"];
01409 searchx = search;
01410 searchy = search;
01411 searchz = search;
01412 } else {
01413 searchx = params.set_default("searchx",3);
01414 searchy = params.set_default("searchy",3);
01415 searchz = params.set_default("searchz",3);
01416 }
01417
01418 float ralt = params.set_default("ralt",180.f);
01419 float rphi = params.set_default("rphi",180.f);
01420 float raz = params.set_default("raz",180.f);
01421 float dalt = params.set_default("dalt",10.f);
01422 float daz = params.set_default("daz",10.f);
01423 float dphi = params.set_default("dphi",10.f);
01424 float threshold = params.set_default("threshold",0.f);
01425 if (threshold < 0.0f) throw InvalidParameterException("The threshold parameter must be greater than or equal to zero");
01426 bool verbose = params.set_default("verbose",false);
01427
01428 vector<Dict> solns;
01429 if (nsoln == 0) return solns;
01430 for (unsigned int i = 0; i < nsoln; ++i ) {
01431 Dict d;
01432 d["score"] = 1.e24;
01433 Transform t;
01434 d["xform.align3d"] = &t;
01435 solns.push_back(d);
01436 }
01437 Dict d;
01438 d["type"] = "eman";
01439 for ( float alt = 0.0f; alt <= ralt; alt += dalt) {
01440
01441
01442 if (verbose) {
01443 cout << "Trying angle " << alt << endl;
01444 }
01445
01446 float begin_az = -raz;
01447 float end_az = raz;
01448 if (alt == 0.0f) {
01449 begin_az = 0.0f;
01450 end_az = 0.0f;
01451 }
01452
01453 for ( float az = begin_az; az <= end_az; az += daz ){
01454 for( float phi = -rphi-az; phi <= rphi-az; phi += dphi ) {
01455 d["alt"] = alt;
01456 d["phi"] = phi;
01457 d["az"] = az;
01458 Transform t(d);
01459 EMData* transformed = this_img->process("xform",Dict("transform",&t));
01460 EMData* ccf = transformed->calc_ccf(to);
01461
01462 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
01463 Dict altered_cmp_params(cmp_params);
01464 if (cmp_name == "dot.tomo") {
01465 altered_cmp_params["ccf"] = ccf;
01466 altered_cmp_params["tx"] = point[0];
01467 altered_cmp_params["ty"] = point[1];
01468 altered_cmp_params["tz"] = point[2];
01469 }
01470
01471 float best_score = transformed->cmp(cmp_name,to,altered_cmp_params);
01472 delete transformed; transformed =0;
01473 delete ccf; ccf = 0;
01474
01475 unsigned int j = 0;
01476 for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
01477 if ( (float)(*it)["score"] > best_score ) {
01478 vector<Dict>::reverse_iterator rit = solns.rbegin();
01479 copy(rit+1,solns.rend()-j,rit);
01480 Dict& d = (*it);
01481 d["score"] = best_score;
01482 d["xform.align3d"] = &t;
01483 break;
01484 }
01485 }
01486 }
01487 }
01488 }
01489
01490 return solns;
01491
01492 }
01493
01494 EMData* RT3DSphereAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
01495 {
01496
01497 vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
01498
01499 Dict t;
01500 Transform* tr = (Transform*) alis[0]["xform.align3d"];
01501 t["transform"] = tr;
01502 EMData* soln = this_img->process("xform",t);
01503 soln->set_attr("xform.align3d",tr);
01504 delete tr; tr = 0;
01505
01506 return soln;
01507
01508 }
01509
01510 vector<Dict> RT3DSphereAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const {
01511
01512 if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
01513 throw ImageDimensionException("This aligner only works for 3D images");
01514 }
01515
01516 int searchx = 0;
01517 int searchy = 0;
01518 int searchz = 0;
01519
01520 if (params.has_key("search")) {
01521 vector<string> check;
01522 check.push_back("searchx");
01523 check.push_back("searchy");
01524 check.push_back("searchz");
01525 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
01526 if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
01527 }
01528 int search = params["search"];
01529 searchx = search;
01530 searchy = search;
01531 searchz = search;
01532 } else {
01533 searchx = params.set_default("searchx",3);
01534 searchy = params.set_default("searchy",3);
01535 searchz = params.set_default("searchz",3);
01536 }
01537
01538 float rphi = params.set_default("rphi",180.f);
01539 float dphi = params.set_default("dphi",10.f);
01540 float threshold = params.set_default("threshold",0.f);
01541 if (threshold < 0.0f) throw InvalidParameterException("The threshold parameter must be greater than or equal to zero");
01542 bool verbose = params.set_default("verbose",false);
01543
01544 vector<Dict> solns;
01545 if (nsoln == 0) return solns;
01546 for (unsigned int i = 0; i < nsoln; ++i ) {
01547 Dict d;
01548 d["score"] = 1.e24;
01549 Transform t;
01550 d["xform.align3d"] = &t;
01551 solns.push_back(d);
01552 }
01553
01554 Dict d;
01555 d["inc_mirror"] = true;
01556 if ( params.has_key("delta") && params.has_key("n") ) {
01557 throw InvalidParameterException("The delta and n parameters are mutually exclusive in the RT3DSphereAligner aligner");
01558 } else if ( params.has_key("n") ) {
01559 d["n"] = params["n"];
01560 } else {
01561
01562 d["delta"] = params.set_default("delta",10.f);
01563 }
01564
01565 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
01566 vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
01567
01568 float verbose_alt = -1.0f;;
01569 for(vector<Transform>::const_iterator trans_it = transforms.begin(); trans_it != transforms.end(); trans_it++) {
01570 Dict params = trans_it->get_params("eman");
01571 float az = params["az"];
01572 if (verbose) {
01573 float alt = params["alt"];
01574 if ( alt != verbose_alt ) {
01575 verbose_alt = alt;
01576 cout << "Trying angle " << alt << endl;
01577 }
01578 }
01579 for( float phi = -rphi-az; phi <= rphi-az; phi += dphi ) {
01580 params["phi"] = phi;
01581 Transform t(params);
01582 EMData* transformed = this_img->process("xform",Dict("transform",&t));
01583 EMData* ccf = transformed->calc_ccf(to);
01584 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
01585 Dict altered_cmp_params(cmp_params);
01586 if (cmp_name == "dot.tomo") {
01587 altered_cmp_params["ccf"] = ccf;
01588 altered_cmp_params["tx"] = point[0];
01589 altered_cmp_params["ty"] = point[1];
01590 altered_cmp_params["tz"] = point[2];
01591 }
01592
01593 float best_score = transformed->cmp(cmp_name,to,altered_cmp_params);
01594 delete transformed; transformed =0;
01595 delete ccf; ccf =0;
01596
01597 unsigned int j = 0;
01598 for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
01599 if ( (float)(*it)["score"] > best_score ) {
01600 vector<Dict>::reverse_iterator rit = solns.rbegin();
01601 copy(rit+1,solns.rend()-j,rit);
01602 Dict& d = (*it);
01603 d["score"] = best_score;
01604 d["xform.align3d"] = &t;
01605 break;
01606 }
01607 }
01608
01609 }
01610 }
01611 delete sym; sym = 0;
01612 return solns;
01613
01614 }
01615
01616 CUDA_Aligner::CUDA_Aligner() {
01617 image_stack = NULL;
01618 image_stack_filtered = NULL;
01619 ccf = NULL;
01620 }
01621
01622 #ifdef EMAN2_USING_CUDA
01623 void CUDA_Aligner::finish() {
01624 if (image_stack) delete image_stack;
01625 if (image_stack_filtered) delete image_stack_filtered;
01626 if (ccf) delete ccf;
01627 }
01628
01629 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) {
01630
01631 NIMA = nima;
01632 NX = nx;
01633 NY = ny;
01634 RING_LENGTH = ring_length;
01635 NRING = nring;
01636 STEP = step;
01637 KX = kx;
01638 KY = ky;
01639 OU = ou;
01640 CTF = ctf;
01641
01642 image_stack = (float *)malloc(NIMA*NX*NY*sizeof(float));
01643 if (CTF == 1) image_stack_filtered = (float *)malloc(NIMA*NX*NY*sizeof(float));
01644 ccf = (float *)malloc(2*(2*KX+1)*(2*KY+1)*NIMA*(RING_LENGTH+2)*sizeof(float));
01645 }
01646
01647 void CUDA_Aligner::insert_image(EMData *image, int num) {
01648
01649 int base_address = num*NX*NY;
01650
01651 for (int y=0; y<NY; y++)
01652 for (int x=0; x<NX; x++)
01653 image_stack[base_address+y*NX+x] = (*image)(x, y);
01654 }
01655
01656 void CUDA_Aligner::filter_stack(vector<float> ctf_params, int id) {
01657
01658 float *params;
01659
01660 params = (float *)malloc(NIMA*6*sizeof(float));
01661
01662 for (int i=0; i<NIMA*6; i++) params[i] = ctf_params[i];
01663
01664 filter_image(image_stack, image_stack_filtered, NIMA, NX, NY, params, id);
01665
01666 free(params);
01667 }
01668
01669 void CUDA_Aligner::sum_oe(vector<float> ctf_params, vector<float> ali_params, EMData *ave1, EMData *ave2, int id) {
01670
01671 float *ctf_p, *ali_p, *av1, *av2;
01672
01673 ctf_p = (float *)malloc(NIMA*6*sizeof(float));
01674 ali_p = (float *)malloc(NIMA*4*sizeof(float));
01675
01676 if (CTF == 1) {
01677 for (int i=0; i<NIMA*6; i++) ctf_p[i] = ctf_params[i];
01678 }
01679 for (int i=0; i<NIMA*4; i++) ali_p[i] = ali_params[i];
01680
01681 av1 = ave1->get_data();
01682 av2 = ave2->get_data();
01683
01684 rot_filt_sum(image_stack, NIMA, NX, NY, CTF, ctf_p, ali_p, av1, av2, id);
01685
01686 free(ctf_p);
01687 free(ali_p);
01688 }
01689
01690 vector<float> CUDA_Aligner::alignment_2d(EMData *ref_image_em, vector<float> sx_list, vector<float> sy_list, int id, int silent) {
01691
01692 float *ref_image, max_ccf;
01693 int base_address, ccf_offset;
01694 float ts, tm;
01695 float ang, sx, sy, mirror, co, so, sxs, sys;
01696 float *sx2, *sy2;
01697 vector<float> align_result;
01698
01699 sx2 = (float *)malloc(NIMA*sizeof(float));
01700 sy2 = (float *)malloc(NIMA*sizeof(float));
01701
01702 ref_image = ref_image_em->get_data();
01703
01704 for (int i=0; i<NIMA; i++) {
01705 sx2[i] = sx_list[i];
01706 sy2[i] = sy_list[i];
01707 }
01708
01709 if (CTF == 1) {
01710 calculate_ccf(image_stack_filtered, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, id, silent);
01711 } else {
01712 calculate_ccf(image_stack, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, id, silent);
01713 }
01714
01715 ccf_offset = NIMA*(RING_LENGTH+2)*(2*KX+1)*(2*KY+1);
01716
01717 for (int im=0; im<NIMA; im++) {
01718 max_ccf = -1.0e22;
01719 for (int kx=-KX; kx<=KX; kx++) {
01720 for (int ky=-KY; ky<=KY; ky++) {
01721 base_address = (((ky+KY)*(2*KX+1)+(kx+KX))*NIMA+im)*(RING_LENGTH+2);
01722 for (int l=0; l<RING_LENGTH; l++) {
01723 ts = ccf[base_address+l];
01724 tm = ccf[base_address+l+ccf_offset];
01725 if (ts > max_ccf) {
01726 ang = float(l)/RING_LENGTH*360.0;
01727 sx = -kx*STEP;
01728 sy = -ky*STEP;
01729 mirror = 0;
01730 max_ccf = ts;
01731 }
01732 if (tm > max_ccf) {
01733 ang = float(l)/RING_LENGTH*360.0;
01734 sx = -kx*STEP;
01735 sy = -ky*STEP;
01736 mirror = 1;
01737 max_ccf = tm;
01738 }
01739 }
01740 }
01741 }
01742 co = cos(ang*M_PI/180.0);
01743 so = -sin(ang*M_PI/180.0);
01744 sxs = sx*co - sy*so;
01745 sys = sx*so + sy*co;
01746
01747 align_result.push_back(ang);
01748 align_result.push_back(sxs);
01749 align_result.push_back(sys);
01750 align_result.push_back(mirror);
01751 }
01752
01753 free(sx2);
01754 free(sy2);
01755
01756 return align_result;
01757 }
01758
01759
01760 vector<float> CUDA_Aligner::ali2d_single_iter(EMData *ref_image_em, vector<float> ali_params, float csx, float csy, int id, int silent) {
01761
01762 float *ref_image, max_ccf;
01763 int base_address, ccf_offset;
01764 float ts, tm;
01765 float ang, sx, sy, mirror, co, so, sxs, sys;
01766 float *sx2, *sy2;
01767 vector<float> align_result;
01768
01769 sx2 = (float *)malloc(NIMA*sizeof(float));
01770 sy2 = (float *)malloc(NIMA*sizeof(float));
01771
01772 ref_image = ref_image_em->get_data();
01773
01774 for (int i=0; i<NIMA; i++) {
01775 ang = ali_params[i*4]/180.0*M_PI;
01776 sx = (ali_params[i*4+3] < 0.5)?(ali_params[i*4+1]-csx):(ali_params[i*4+1]+csx);
01777 sy = ali_params[i*4+2]-csy;
01778 co = cos(ang);
01779 so = sin(ang);
01780 sx2[i] = -(sx*co-sy*so);
01781 sy2[i] = -(sx*so+sy*co);
01782 }
01783
01784 if (CTF == 1) {
01785 calculate_ccf(image_stack_filtered, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, id, silent);
01786 } else {
01787 calculate_ccf(image_stack, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, id, silent);
01788 }
01789
01790 ccf_offset = NIMA*(RING_LENGTH+2)*(2*KX+1)*(2*KY+1);
01791
01792 float sx_sum = 0.0f;
01793 float sy_sum = 0.0f;
01794
01795 for (int im=0; im<NIMA; im++) {
01796 max_ccf = -1.0e22;
01797 for (int kx=-KX; kx<=KX; kx++) {
01798 for (int ky=-KY; ky<=KY; ky++) {
01799 base_address = (((ky+KY)*(2*KX+1)+(kx+KX))*NIMA+im)*(RING_LENGTH+2);
01800 for (int l=0; l<RING_LENGTH; l++) {
01801 ts = ccf[base_address+l];
01802 tm = ccf[base_address+l+ccf_offset];
01803 if (ts > max_ccf) {
01804 ang = float(l)/RING_LENGTH*360.0;
01805 sx = -kx*STEP;
01806 sy = -ky*STEP;
01807 mirror = 0;
01808 max_ccf = ts;
01809 }
01810 if (tm > max_ccf) {
01811 ang = float(l)/RING_LENGTH*360.0;
01812 sx = -kx*STEP;
01813 sy = -ky*STEP;
01814 mirror = 1;
01815 max_ccf = tm;
01816 }
01817 }
01818 }
01819 }
01820 co = cos(ang*M_PI/180.0);
01821 so = -sin(ang*M_PI/180.0);
01822
01823 sxs = (sx-sx2[im])*co-(sy-sy2[im])*so;
01824 sys = (sx-sx2[im])*so+(sy-sy2[im])*co;
01825
01826 align_result.push_back(ang);
01827 align_result.push_back(sxs);
01828 align_result.push_back(sys);
01829 align_result.push_back(mirror);
01830
01831 if (mirror == 0) { sx_sum += sxs; } else { sx_sum -= sxs; }
01832 sy_sum += sys;
01833 }
01834
01835 align_result.push_back(sx_sum);
01836 align_result.push_back(sy_sum);
01837
01838 free(sx2);
01839 free(sy2);
01840
01841 return align_result;
01842 }
01843
01844
01845 #endif
01846
01847 void EMAN::dump_aligners()
01848 {
01849 dump_factory < Aligner > ();
01850 }
01851
01852 map<string, vector<string> > EMAN::dump_aligners_list()
01853 {
01854 return dump_factory_list < Aligner > ();
01855 }