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