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