aligner.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "cmp.h"
00037 #include "aligner.h"
00038 #include "emdata.h"
00039 #include "processor.h"
00040 #include "util.h"
00041 #include "symmetry.h"
00042 #include <gsl/gsl_multimin.h>
00043 
00044 #ifdef EMAN2_USING_CUDA
00045         #include <sparx/cuda/cuda_ccf.h>
00046 #endif
00047 
00048 #define EMAN2_ALIGNER_DEBUG 0
00049 
00050 using namespace EMAN;
00051 
00052 template <> Factory < Aligner >::Factory()
00053 {
00054         force_add(&TranslationalAligner::NEW);
00055         force_add(&RotationalAligner::NEW);
00056         force_add(&RotatePrecenterAligner::NEW);
00057         force_add(&RotateTranslateAligner::NEW);
00058         force_add(&RotateFlipAligner::NEW);
00059         force_add(&RotateTranslateFlipAligner::NEW);
00060         force_add(&RTFExhaustiveAligner::NEW);
00061         force_add(&RTFSlowExhaustiveAligner::NEW);
00062         force_add(&RefineAligner::NEW);
00063         force_add(&Refine3DAligner::NEW);
00064         force_add(&RT3DGridAligner::NEW);
00065         force_add(&RT3DSphereAligner::NEW);
00066 }
00067 
00068 // Note, the translational aligner assumes that the correlation image
00069 // generated by the calc_ccf function is centered on the bottom left corner
00070 // That is, if you did at calc_cff using identical images, the
00071 // peak would be at 0,0
00072 EMData *TranslationalAligner::align(EMData * this_img, EMData *to,
00073                                         const string&, const Dict&) const
00074 {
00075         if (!this_img) {
00076                 return 0;
00077         }
00078 
00079         if (to && !EMUtil::is_same_size(this_img, to))
00080                 throw ImageDimensionException("Images must be the same size to perform translational alignment");
00081 
00082         EMData *cf = 0;
00083         int nx = this_img->get_xsize();
00084         int ny = this_img->get_ysize();
00085         int nz = this_img->get_zsize();
00086 
00087         int masked = params.set_default("masked",0);
00088         bool use_cpu = true;
00089 #ifdef EMAN2_USING_CUDA
00090         if (this_img->gpu_operation_preferred() ) {
00091 //              cout << "Translate on GPU" << endl;
00092                 use_cpu = false;
00093                 cf = this_img->calc_ccf_cuda(to,false,false);
00094         }
00095 #endif // EMAN2_USING_CUDA
00096         if (use_cpu) cf = this_img->calc_ccf(to);
00097 
00098         // This is too expensive
00099         if (masked) {
00100                 EMData *msk=this_img->process("threshold.notzero");
00101                 EMData *sqr=to->process("math.squared");
00102                 EMData *cfn=msk->calc_ccf(sqr);
00103                 cfn->process_inplace("math.sqrt");
00104                 float *d1=cf->get_data();
00105                 float *d2=cfn->get_data();
00106                 for (int i=0; i<nx*ny*nz; i++) {
00107                         if (d2[i]!=0) d1[i]/=d2[i];
00108                 }
00109                 cf->update();
00110                 delete msk;
00111                 delete sqr;
00112                 delete cfn;
00113         }
00114 
00115 //
00116 
00117         int maxshiftx = params.set_default("maxshift",-1);
00118         int maxshifty = params["maxshift"];
00119         int maxshiftz = params["maxshift"];
00120         int nozero = params["nozero"];
00121 
00122         if (maxshiftx <= 0) {
00123                 maxshiftx = nx / 4;
00124                 maxshifty = ny / 4;
00125                 maxshiftz = nz / 4;
00126         }
00127 
00128         if (maxshiftx > nx / 2 - 1) maxshiftx = nx / 2 - 1;
00129         if (maxshifty > ny / 2 - 1)     maxshifty = ny / 2 - 1;
00130         if (maxshiftz > nz / 2 - 1) maxshiftz = nz / 2 - 1;
00131 
00132         if (nx == 1) maxshiftx = 0; // This is justhere for completeness really... plus it saves errors
00133         if (ny == 1) maxshifty = 0;
00134         if (nz == 1) maxshiftz = 0;
00135 
00136         // If nozero the portion of the image in the center (and its 8-connected neighborhood) is zeroed
00137         if (nozero) {
00138                 cf->zero_corner_circulant(1);
00139         }
00140 
00141         IntPoint peak;
00142 #ifdef EMAN2_USING_CUDA
00143         if (!use_cpu) {
00144                 EMDataForCuda tmp = cf->get_data_struct_for_cuda();
00145                 int* p = calc_max_location_wrap_cuda(&tmp,maxshiftx, maxshifty, maxshiftz);
00146                 peak = IntPoint(p[0],p[1],p[2]);
00147                 free(p);
00148         }
00149 #endif // EMAN2_USING_CUDA
00150         if (use_cpu) {
00151                 peak = cf->calc_max_location_wrap(maxshiftx, maxshifty, maxshiftz);
00152         }
00153         Vec3f cur_trans = Vec3f ( (float)-peak[0], (float)-peak[1], (float)-peak[2]);
00154 
00155         if (!to) {
00156                 cur_trans /= 2.0f; // If aligning theimage to itself then only go half way -
00157                 int intonly = params.set_default("intonly",false);
00158                 if (intonly) {
00159                         cur_trans[0] = floor(cur_trans[0] + 0.5f);
00160                         cur_trans[1] = floor(cur_trans[1] + 0.5f);
00161                         cur_trans[2] = floor(cur_trans[2] + 0.5f);
00162                 }
00163         }
00164 
00165         if( cf ){
00166                 delete cf;
00167                 cf = 0;
00168         }
00169         Dict params("trans",static_cast< vector<int> >(cur_trans));
00170         cf=this_img->process("math.translate.int",params);
00171         Transform t;
00172         t.set_trans(cur_trans);
00173         if ( nz != 1 ) {
00174 //              Transform* t = get_set_align_attr("xform.align3d",cf,this_img);
00175 //              t->set_trans(cur_trans);
00176                 cf->set_attr("xform.align3d",&t);
00177         } else if ( ny != 1 ) {
00178                 //Transform* t = get_set_align_attr("xform.align2d",cf,this_img);
00179                 cur_trans[2] = 0; // just make sure of it
00180                 t.set_trans(cur_trans);
00181                 cf->set_attr("xform.align2d",&t);
00182         }
00183 
00184         return cf;
00185 }
00186 
00187 EMData * RotationalAligner::align_180_ambiguous(EMData * this_img, EMData * to, int rfp_mode) {
00188 
00189         // Make translationally invariant rotational footprints
00190         EMData* this_img_rfp, * to_rfp;
00191         if (rfp_mode == 0) {
00192                 this_img_rfp = this_img->make_rotational_footprint_e1();
00193                 to_rfp = to->make_rotational_footprint_e1();
00194         } else if (rfp_mode == 1) {
00195                 this_img_rfp = this_img->make_rotational_footprint();
00196                 to_rfp = to->make_rotational_footprint();
00197         } else if (rfp_mode == 2) {
00198                 this_img_rfp = this_img->make_rotational_footprint_cmc();
00199                 to_rfp = to->make_rotational_footprint_cmc();
00200         } else {
00201                 throw InvalidParameterException("rfp_mode must be 0,1 or 2");
00202         }
00203         int this_img_rfp_nx = this_img_rfp->get_xsize();
00204 
00205         // Do row-wise correlation, returning a sum.
00206         EMData *cf = this_img_rfp->calc_ccfx(to_rfp, 0, this_img->get_ysize());
00207 
00208         // Delete them, they're no longer needed
00209         delete this_img_rfp; this_img_rfp = 0;
00210         delete to_rfp; to_rfp = 0;
00211 
00212         // Now solve the rotational alignment by finding the max in the column sum
00213         float *data = cf->get_data();
00214         float peak = 0;
00215         int peak_index = 0;
00216         Util::find_max(data, this_img_rfp_nx, &peak, &peak_index);
00217 
00218         if( cf ) {
00219                 delete cf;
00220                 cf = 0;
00221         }
00222         float rot_angle = (float) (peak_index * 180.0f / this_img_rfp_nx);
00223 
00224         // Return the result
00225         Transform tmp(Dict("type","2d","alpha",rot_angle));
00226         cf=this_img->process("math.transform",Dict("transform",(Transform*)&tmp));
00227 //      Transform* t = get_set_align_attr("xform.align2d",cf,this_img);
00228 //      Dict d("type","2d","alpha",rot_angle);
00229 //      t->set_rotation(d);
00230         cf->set_attr("xform.align2d",&tmp);
00231         return cf;
00232 }
00233 
00234 EMData *RotationalAligner::align(EMData * this_img, EMData *to,
00235                         const string& cmp_name, const Dict& cmp_params) const
00236 {
00237         if (!to) throw InvalidParameterException("Can not rotational align - the image to align to is NULL");
00238 
00239         // Perform 180 ambiguous alignment
00240         int rfp_mode = params.set_default("rfp_mode",0);
00241         EMData* rot_aligned = RotationalAligner::align_180_ambiguous(this_img,to,rfp_mode);
00242         Transform * tmp = rot_aligned->get_attr("xform.align2d");
00243         Dict rot = tmp->get_rotation("2d");
00244         float rotate_angle_solution = rot["alpha"];
00245         delete tmp;
00246 
00247         EMData *rot_align_180 = rot_aligned->process("math.rotate.180");
00248 
00249         // Generate the comparison metrics for both rotational candidates
00250         float rot_cmp = rot_aligned->cmp(cmp_name, to, cmp_params);
00251         float rot_180_cmp = rot_align_180->cmp(cmp_name, to, cmp_params);
00252         // Decide on the result
00253         float score = 0.0;
00254         EMData* result = NULL;
00255         if (rot_cmp < rot_180_cmp){
00256                 result = rot_aligned;
00257                 score = rot_cmp;
00258                 delete rot_align_180; rot_align_180 = 0;
00259         } else {
00260                 result = rot_align_180;
00261                 score = rot_180_cmp;
00262                 delete rot_aligned; rot_aligned = 0;
00263                 rotate_angle_solution = rotate_angle_solution-180.0f;
00264         }
00265 
00266 //      Transform* t = get_align_attr("xform.align2d",result);
00267 //      t->set_rotation(Dict("type","2d","alpha",rotate_angle_solution));
00268         Transform tmp2(Dict("type","2d","alpha",rotate_angle_solution));
00269         result->set_attr("xform.align2d",&tmp2);
00270         return result;
00271 }
00272 
00273 
00274 EMData *RotatePrecenterAligner::align(EMData * this_img, EMData *to,
00275                         const string&, const Dict&) const
00276 {
00277         if (!to) {
00278                 return 0;
00279         }
00280 
00281         int ny = this_img->get_ysize();
00282         int size = Util::calc_best_fft_size((int) (M_PI * ny * 1.5));
00283         EMData *e1 = this_img->unwrap(4, ny * 7 / 16, size, 0, 0, 1);
00284         EMData *e2 = to->unwrap(4, ny * 7 / 16, size, 0, 0, 1);
00285         EMData *cf = e1->calc_ccfx(e2, 0, ny);
00286 
00287         float *data = cf->get_data();
00288 
00289         float peak = 0;
00290         int peak_index = 0;
00291         Util::find_max(data, size, &peak, &peak_index);
00292         float a = (float) ((1.0f - 1.0f * peak_index / size) * 180. * 2);
00293 
00294         Transform rot;
00295         rot.set_rotation(Dict("type","2d","alpha",(float)(a*180./M_PI)));
00296         EMData* rslt = this_img->process("math.transform",Dict("transform",&rot));
00297         rslt->set_attr("xform.align2d",&rot);
00298 //
00299 //      Transform* t = get_set_align_attr("xform.align2d",rslt,this_img);
00300 //      t->set_rotation(Dict("type","2d","alpha",-a));
00301 //
00302 //      EMData* result this_img->transform(Dict("type","2d","alpha",(float)(a*180./M_PI)));
00303 //
00304 //      cf->set_attr("xform.align2d",t);
00305 //      delete t;
00306 //      cf->update();
00307 
00308         if( e1 )
00309         {
00310                 delete e1;
00311                 e1 = 0;
00312         }
00313 
00314         if( e2 )
00315         {
00316                 delete e2;
00317                 e2 = 0;
00318         }
00319 
00320         if (cf) {
00321                 delete cf;
00322                 cf = 0;
00323         }
00324 
00325         return rslt;
00326 }
00327 
00328 EMData *RotateTranslateAligner::align(EMData * this_img, EMData *to,
00329                         const string & cmp_name, const Dict& cmp_params) const
00330 {
00331         // Get the 180 degree ambiguously rotationally aligned and its 180 degree rotation counterpart
00332         int rfp_mode = params.set_default("rfp_mode",0);
00333         EMData *rot_align  =  RotationalAligner::align_180_ambiguous(this_img,to,rfp_mode);
00334         Transform * tmp = rot_align->get_attr("xform.align2d");
00335         Dict rot = tmp->get_rotation("2d");
00336         float rotate_angle_solution = rot["alpha"];
00337         delete tmp;
00338 
00339         EMData *rot_align_180 = rot_align->copy();
00340         rot_align_180->process_inplace("math.rotate.180");
00341 
00342         Dict trans_params;
00343         trans_params["intonly"]  = 0;
00344         trans_params["maxshift"] = params.set_default("maxshift", -1);
00345 
00346         // Do the first case translational alignment
00347         trans_params["nozero"]   = params.set_default("nozero",false);
00348         EMData* rot_trans = rot_align->align("translational", to, trans_params, cmp_name, cmp_params);
00349         if( rot_align ) { // Clean up
00350                 delete rot_align;
00351                 rot_align = 0;
00352         }
00353 
00354         // Do the second case translational alignment
00355         EMData*  rot_180_trans = rot_align_180->align("translational", to, trans_params, cmp_name, cmp_params);
00356         if( rot_align_180 )     { // Clean up
00357                 delete rot_align_180;
00358                 rot_align_180 = 0;
00359         }
00360 
00361         // Finally decide on the result
00362         float cmp1 = rot_trans->cmp(cmp_name, to, cmp_params);
00363         float cmp2 = rot_180_trans->cmp(cmp_name, to, cmp_params);
00364 
00365         EMData *result = 0;
00366         if (cmp1 < cmp2) { // Assumes smaller is better - thus all comparitors should support "smaller is better"
00367                 if( rot_180_trans )     {
00368                         delete rot_180_trans;
00369                         rot_180_trans = 0;
00370                 }
00371                 result = rot_trans;
00372         }
00373         else {
00374                 if( rot_trans ) {
00375                         delete rot_trans;
00376                         rot_trans = 0;
00377                 }
00378                 result = rot_180_trans;
00379                 rotate_angle_solution -= 180.f;
00380         }
00381 
00382         Transform* t = result->get_attr("xform.align2d");
00383         t->set_rotation(Dict("type","2d","alpha",rotate_angle_solution));
00384         result->set_attr("xform.align2d",t);
00385         delete t;
00386 
00387         return result;
00388 }
00389 
00390 
00391 
00392 
00393 EMData* RotateTranslateFlipAligner::align(EMData * this_img, EMData *to,
00394                                                                                   const string & cmp_name, const Dict& cmp_params) const
00395 {
00396         // Get the non flipped rotational, tranlsationally aligned image
00397         Dict rt_params("maxshift", params["maxshift"], "rfp_mode", params.set_default("rfp_mode",0));
00398         EMData *rot_trans_align = this_img->align("rotate_translate",to,rt_params,cmp_name, cmp_params);
00399 
00400         // Do the same alignment, but using the flipped version of the image
00401         EMData *flipped = params.set_default("flip", (EMData *) 0);
00402         bool delete_flag = false;
00403         if (flipped == 0) {
00404                 flipped = to->process("xform.flip", Dict("axis", "x"));
00405                 delete_flag = true;
00406         }
00407 
00408         EMData * rot_trans_align_flip = this_img->align("rotate_translate", flipped, rt_params, cmp_name, cmp_params);
00409         Transform* t = rot_trans_align_flip->get_attr("xform.align2d");
00410         t->set_mirror(true);
00411         rot_trans_align_flip->set_attr("xform.align2d",t);
00412         delete t;
00413 
00414         // Now finally decide on what is the best answer
00415         float cmp1 = rot_trans_align->cmp(cmp_name, to, cmp_params);
00416         float cmp2 = rot_trans_align_flip->cmp(cmp_name, flipped, cmp_params);
00417 
00418         if (delete_flag){
00419                 if(flipped) {
00420                         delete flipped;
00421                         flipped = 0;
00422                 }
00423         }
00424 
00425         EMData *result = 0;
00426         if (cmp1 < cmp2 )  {
00427 
00428                 if( rot_trans_align_flip ) {
00429                         delete rot_trans_align_flip;
00430                         rot_trans_align_flip = 0;
00431                 }
00432                 result = rot_trans_align;
00433         }
00434         else {
00435                 if( rot_trans_align ) {
00436                         delete rot_trans_align;
00437                         rot_trans_align = 0;
00438                 }
00439                 result = rot_trans_align_flip;
00440                 result->process_inplace("xform.flip",Dict("axis","x"));
00441         }
00442 
00443         return result;
00444 }
00445 
00446 
00447 
00448 
00449 EMData *RotateFlipAligner::align(EMData * this_img, EMData *to,
00450                         const string& cmp_name, const Dict& cmp_params) const
00451 {
00452         Dict rot_params("rfp_mode",params.set_default("rfp_mode",0));
00453         EMData *r1 = this_img->align("rotational", to, rot_params,cmp_name, cmp_params);
00454 
00455 
00456         EMData* flipped =to->process("xform.flip", Dict("axis", "x"));
00457         EMData *r2 = this_img->align("rotational", flipped,rot_params, cmp_name, cmp_params);
00458         Transform* t = r2->get_attr("xform.align2d");
00459         t->set_mirror(true);
00460         r2->set_attr("xform.align2d",t);
00461         delete t;
00462 
00463         float cmp1 = r1->cmp(cmp_name, to, cmp_params);
00464         float cmp2 = r2->cmp(cmp_name, flipped, cmp_params);
00465 
00466         delete flipped; flipped = 0;
00467 
00468         EMData *result = 0;
00469 
00470         if (cmp1 < cmp2) {
00471                 if( r2 )
00472                 {
00473                         delete r2;
00474                         r2 = 0;
00475                 }
00476                 result = r1;
00477         }
00478         else {
00479                 if( r1 )
00480                 {
00481                         delete r1;
00482                         r1 = 0;
00483                 }
00484                 result = r2;
00485                 result->process_inplace("xform.flip",Dict("axis","x"));
00486         }
00487 
00488         return result;
00489 }
00490 
00491 
00492 // David Woolford says FIXME
00493 // You will note the excessive amount of EMData copying that's going in this function
00494 // This is because functions that are operating on the EMData objects are changing them
00495 // and if you do not use copies the whole algorithm breaks. I did not have time to go
00496 // through and rectify this situation.
00497 // David Woolford says - this problem is related to the fact that many functions that
00498 // take EMData pointers as arguments do not take them as constant pointers to constant
00499 // objects, instead they are treated as raw (completely changeable) pointers. This means
00500 // it's hard to track down which functions are changing the EMData objects, because they
00501 // all do (in name). If this behavior is unavoidable then ignore this comment, however if possible it would
00502 // be good to make things const as much as possible. For example in alignment, technically
00503 // the argument EMData objects (raw pointers) should not be altered... should they?
00504 //
00505 // But const can be very annoying sometimes...
00506 EMData *RTFExhaustiveAligner::align(EMData * this_img, EMData *to,
00507                         const string & cmp_name, const Dict& cmp_params) const
00508 {
00509         EMData *flip = params.set_default("flip", (EMData *) 0);
00510         int maxshift = params.set_default("maxshift", this_img->get_xsize()/8);
00511         if (maxshift < 2) throw InvalidParameterException("maxshift must be greater than or equal to 2");
00512 
00513         int ny = this_img->get_ysize();
00514         int xst = (int) floor(2 * M_PI * ny);
00515         xst = Util::calc_best_fft_size(xst);
00516 
00517         Dict d("n",2);
00518         EMData *to_shrunk_unwrapped = to->process("math.medianshrink",d);
00519 
00520         int to_copy_r2 = to_shrunk_unwrapped->get_ysize() / 2 - 2 - maxshift / 2;
00521         EMData *tmp = to_shrunk_unwrapped->unwrap(4, to_copy_r2, xst / 2, 0, 0, true);
00522         if( to_shrunk_unwrapped )
00523         {
00524                 delete to_shrunk_unwrapped;
00525                 to_shrunk_unwrapped = 0;
00526         }
00527         to_shrunk_unwrapped = tmp;
00528 
00529         EMData *to_shrunk_unwrapped_copy = to_shrunk_unwrapped->copy();
00530         EMData* to_unwrapped = to->unwrap(4, to->get_ysize() / 2 - 2 - maxshift, xst, 0, 0, true);
00531         EMData *to_unwrapped_copy = to_unwrapped->copy();
00532 
00533         bool delete_flipped = true;
00534         EMData *flipped = 0;
00535         if (flip) {
00536                 delete_flipped = false;
00537                 flipped = flip;
00538         }
00539         else {
00540                 flipped = to->process("xform.flip", Dict("axis", "x"));
00541         }
00542         EMData *to_shrunk_flipped_unwrapped = flipped->process("math.medianshrink",d);
00543         tmp = to_shrunk_flipped_unwrapped->unwrap(4, to_copy_r2, xst / 2, 0, 0, true);
00544         if( to_shrunk_flipped_unwrapped )
00545         {
00546                 delete to_shrunk_flipped_unwrapped;
00547                 to_shrunk_flipped_unwrapped = 0;
00548         }
00549         to_shrunk_flipped_unwrapped = tmp;
00550         EMData *to_shrunk_flipped_unwrapped_copy = to_shrunk_flipped_unwrapped->copy();
00551         EMData* to_flip_unwrapped = flipped->unwrap(4, to->get_ysize() / 2 - 2 - maxshift, xst, 0, 0, true);
00552         EMData* to_flip_unwrapped_copy = to_flip_unwrapped->copy();
00553 
00554         if (delete_flipped && flipped != 0) {
00555                 delete flipped;
00556                 flipped = 0;
00557         }
00558 
00559         EMData *this_shrunk_2 = this_img->process("math.medianshrink",d);
00560 
00561         float bestval = FLT_MAX;
00562         float bestang = 0;
00563         int bestflip = 0;
00564         float bestdx = 0;
00565         float bestdy = 0;
00566 
00567         int half_maxshift = maxshift / 2;
00568 
00569         int ur2 = this_shrunk_2->get_ysize() / 2 - 2 - half_maxshift;
00570         for (int dy = -half_maxshift; dy <= half_maxshift; dy += 1) {
00571                 for (int dx = -half_maxshift; dx <= half_maxshift; dx += 1) {
00572 #ifdef  _WIN32
00573                         if (_hypot(dx, dy) <= half_maxshift) {
00574 #else
00575                         if (hypot(dx, dy) <= half_maxshift) {
00576 #endif
00577                                 EMData *uw = this_shrunk_2->unwrap(4, ur2, xst / 2, dx, dy, true);
00578                                 EMData *uwc = uw->copy();
00579                                 EMData *a = uw->calc_ccfx(to_shrunk_unwrapped);
00580 
00581                                 uwc->rotate_x(a->calc_max_index());
00582                                 float cm = uwc->cmp(cmp_name, to_shrunk_unwrapped_copy, cmp_params);
00583                                 if (cm < bestval) {
00584                                         bestval = cm;
00585                                         bestang = (float) (2.0 * M_PI * a->calc_max_index() / a->get_xsize());
00586                                         bestdx = (float)dx;
00587                                         bestdy = (float)dy;
00588                                         bestflip = 0;
00589                                 }
00590 
00591 
00592                                 if( a )
00593                                 {
00594                                         delete a;
00595                                         a = 0;
00596                                 }
00597                                 if( uw )
00598                                 {
00599                                         delete uw;
00600                                         uw = 0;
00601                                 }
00602                                 if( uwc )
00603                                 {
00604                                         delete uwc;
00605                                         uwc = 0;
00606                                 }
00607                                 uw = this_shrunk_2->unwrap(4, ur2, xst / 2, dx, dy, true);
00608                                 uwc = uw->copy();
00609                                 a = uw->calc_ccfx(to_shrunk_flipped_unwrapped);
00610 
00611                                 uwc->rotate_x(a->calc_max_index());
00612                                 cm = uwc->cmp(cmp_name, to_shrunk_flipped_unwrapped_copy, cmp_params);
00613                                 if (cm < bestval) {
00614                                         bestval = cm;
00615                                         bestang = (float) (2.0 * M_PI * a->calc_max_index() / a->get_xsize());
00616                                         bestdx = (float)dx;
00617                                         bestdy = (float)dy;
00618                                         bestflip = 1;
00619                                 }
00620 
00621                                 if( a )
00622                                 {
00623                                         delete a;
00624                                         a = 0;
00625                                 }
00626 
00627                                 if( uw )
00628                                 {
00629                                         delete uw;
00630                                         uw = 0;
00631                                 }
00632                                 if( uwc )
00633                                 {
00634                                         delete uwc;
00635                                         uwc = 0;
00636                                 }
00637                         }
00638                 }
00639         }
00640         if( this_shrunk_2 )
00641         {
00642                 delete this_shrunk_2;
00643                 this_shrunk_2 = 0;
00644         }
00645         if( to_shrunk_unwrapped )
00646         {
00647                 delete to_shrunk_unwrapped;
00648                 to_shrunk_unwrapped = 0;
00649         }
00650         if( to_shrunk_unwrapped_copy )
00651         {
00652                 delete to_shrunk_unwrapped_copy;
00653                 to_shrunk_unwrapped_copy = 0;
00654         }
00655         if( to_shrunk_flipped_unwrapped )
00656         {
00657                 delete to_shrunk_flipped_unwrapped;
00658                 to_shrunk_flipped_unwrapped = 0;
00659         }
00660         if( to_shrunk_flipped_unwrapped_copy )
00661         {
00662                 delete to_shrunk_flipped_unwrapped_copy;
00663                 to_shrunk_flipped_unwrapped_copy = 0;
00664         }
00665         bestdx *= 2;
00666         bestdy *= 2;
00667         bestval = FLT_MAX;
00668 
00669         float bestdx2 = bestdx;
00670         float bestdy2 = bestdy;
00671         // Note I tried steps less than 1.0 (sub pixel precision) and it actually appeared detrimental
00672         // So my advice is to stick with dx += 1.0 etc unless you really are looking to fine tune this
00673         // algorithm
00674         for (float dy = bestdy2 - 3; dy <= bestdy2 + 3; dy += 1.0 ) {
00675                 for (float dx = bestdx2 - 3; dx <= bestdx2 + 3; dx += 1.0 ) {
00676 
00677 #ifdef  _WIN32
00678                         if (_hypot(dx, dy) <= maxshift) {
00679 #else
00680                         if (hypot(dx, dy) <= maxshift) {
00681 #endif
00682                                 EMData *uw = this_img->unwrap(4, this_img->get_ysize() / 2 - 2 - maxshift, xst, (int)dx, (int)dy, true);
00683                                 EMData *uwc = uw->copy();
00684                                 EMData *a = uw->calc_ccfx(to_unwrapped);
00685 
00686                                 uwc->rotate_x(a->calc_max_index());
00687                                 float cm = uwc->cmp(cmp_name, to_unwrapped_copy, cmp_params);
00688 
00689                                 if (cm < bestval) {
00690                                         bestval = cm;
00691                                         bestang = (float)(2.0 * M_PI * a->calc_max_index() / a->get_xsize());
00692                                         bestdx = dx;
00693                                         bestdy = dy;
00694                                         bestflip = 0;
00695                                 }
00696 
00697                                 if( a )
00698                                 {
00699                                         delete a;
00700                                         a = 0;
00701                                 }
00702                                 if( uw )
00703                                 {
00704                                         delete uw;
00705                                         uw = 0;
00706                                 }
00707                                 if( uwc )
00708                                 {
00709                                         delete uwc;
00710                                         uwc = 0;
00711                                 }
00712                                 uw = this_img->unwrap(4, this_img->get_ysize() / 2 - 2 - maxshift, xst, (int)dx, (int)dy, true);
00713                                 uwc = uw->copy();
00714                                 a = uw->calc_ccfx(to_flip_unwrapped);
00715 
00716                                 uwc->rotate_x(a->calc_max_index());
00717                                 cm = uwc->cmp(cmp_name, to_flip_unwrapped_copy, cmp_params);
00718 
00719                                 if (cm < bestval) {
00720                                         bestval = cm;
00721                                         bestang = (float)(2.0 * M_PI * a->calc_max_index() / a->get_xsize());
00722                                         bestdx = dx;
00723                                         bestdy = dy;
00724                                         bestflip = 1;
00725                                 }
00726 
00727                                 if( a )
00728                                 {
00729                                         delete a;
00730                                         a = 0;
00731                                 }
00732                                 if( uw )
00733                                 {
00734                                         delete uw;
00735                                         uw = 0;
00736                                 }
00737                                 if( uwc )
00738                                 {
00739                                         delete uwc;
00740                                         uwc = 0;
00741                                 }
00742                         }
00743                 }
00744         }
00745         if( to_unwrapped ) {delete to_unwrapped;to_unwrapped = 0;}
00746         if( to_shrunk_unwrapped ) {     delete to_shrunk_unwrapped;     to_shrunk_unwrapped = 0;}
00747         if (to_unwrapped_copy) { delete to_unwrapped_copy; to_unwrapped_copy = 0; }
00748         if (to_flip_unwrapped) { delete to_flip_unwrapped; to_flip_unwrapped = 0; }
00749         if (to_flip_unwrapped_copy) { delete to_flip_unwrapped_copy; to_flip_unwrapped_copy = 0;}
00750 
00751         bestang *= (float)EMConsts::rad2deg;
00752         Transform t(Dict("type","2d","alpha",(float)bestang));
00753         t.set_pre_trans(Vec2f(-bestdx,-bestdy));
00754         if (bestflip) {
00755                 t.set_mirror(true);
00756         }
00757 
00758         EMData* ret = this_img->process("math.transform",Dict("transform",&t));
00759         ret->set_attr("xform.align2d",&t);
00760 
00761         return ret;
00762 }
00763 
00764 
00765 EMData *RTFSlowExhaustiveAligner::align(EMData * this_img, EMData *to,
00766                         const string & cmp_name, const Dict& cmp_params) const
00767 {
00768 
00769         EMData *flip = params.set_default("flip", (EMData *) 0);
00770         int maxshift = params.set_default("maxshift", -1);
00771 
00772         EMData *flipped = 0;
00773 
00774         bool delete_flipped = true;
00775         if (flip) {
00776                 delete_flipped = false;
00777                 flipped = flip;
00778         }
00779         else {
00780                 flipped = to->process("xform.flip", Dict("axis", "x"));
00781         }
00782 
00783         int nx = this_img->get_xsize();
00784 
00785         if (maxshift < 0) {
00786                 maxshift = nx / 10;
00787         }
00788 
00789         float angle_step =  params.set_default("angstep", 0.0f);
00790         if ( angle_step == 0 ) angle_step = atan2(2.0f, (float)nx);
00791         else {
00792                 angle_step *= (float)EMConsts::deg2rad; //convert to radians
00793         }
00794         float trans_step =  params.set_default("transtep",1.0f);
00795 
00796         if (trans_step <= 0) throw InvalidParameterException("transstep must be greater than 0");
00797         if (angle_step <= 0) throw InvalidParameterException("angstep must be greater than 0");
00798 
00799 
00800         Dict shrinkfactor("n",2);
00801         EMData *this_img_shrink = this_img->process("math.medianshrink",shrinkfactor);
00802         EMData *to_shrunk = to->process("math.medianshrink",shrinkfactor);
00803         EMData *flipped_shrunk = flipped->process("math.medianshrink",shrinkfactor);
00804 
00805         int bestflip = 0;
00806         float bestdx = 0;
00807         float bestdy = 0;
00808 
00809         float bestang = 0;
00810         float bestval = FLT_MAX;
00811 
00812         int half_maxshift = maxshift / 2;
00813 
00814 
00815         for (int dy = -half_maxshift; dy <= half_maxshift; ++dy) {
00816                 for (int dx = -half_maxshift; dx <= half_maxshift; ++dx) {
00817                         if (hypot(dx, dy) <= maxshift) {
00818                                 for (float ang = -angle_step * 2.0f; ang <= (float)2 * M_PI; ang += angle_step * 4.0f) {
00819                                         EMData v(*this_img_shrink);
00820                                         Transform t(Dict("type","2d","alpha",static_cast<float>(ang*EMConsts::rad2deg)));
00821                                         t.set_trans((float)dx,(float)dy);
00822                                         v.transform(t);
00823 //                                      v.rotate_translate(ang*EMConsts::rad2deg, 0.0f, 0.0f, (float)dx, (float)dy, 0.0f);
00824 
00825                                         float lc = v.cmp(cmp_name, to_shrunk, cmp_params);
00826 
00827                                         if (lc < bestval) {
00828                                                 bestval = lc;
00829                                                 bestang = ang;
00830                                                 bestdx = (float)dx;
00831                                                 bestdy = (float)dy;
00832                                                 bestflip = 0;
00833                                         }
00834 
00835                                         lc = v.cmp(cmp_name,flipped_shrunk , cmp_params);
00836                                         if (lc < bestval) {
00837                                                 bestval = lc;
00838                                                 bestang = ang;
00839                                                 bestdx = (float)dx;
00840                                                 bestdy = (float)dy;
00841                                                 bestflip = 1;
00842                                         }
00843                                 }
00844                         }
00845                 }
00846         }
00847 
00848         if( to_shrunk )
00849         {
00850                 delete to_shrunk;
00851                 to_shrunk = 0;
00852         }
00853         if( flipped_shrunk )
00854         {
00855                 delete flipped_shrunk;
00856                 flipped_shrunk = 0;
00857         }
00858         if( this_img_shrink )
00859         {
00860                 delete this_img_shrink;
00861                 this_img_shrink = 0;
00862         }
00863 
00864         bestdx *= 2;
00865         bestdy *= 2;
00866         bestval = FLT_MAX;
00867 
00868         float bestdx2 = bestdx;
00869         float bestdy2 = bestdy;
00870         float bestang2 = bestang;
00871 
00872         for (float dy = bestdy2 - 3; dy <= bestdy2 + 3; dy += trans_step) {
00873                 for (float dx = bestdx2 - 3; dx <= bestdx2 + 3; dx += trans_step) {
00874                         if (hypot(dx, dy) <= maxshift) {
00875                                 for (float ang = bestang2 - angle_step * 6.0f; ang <= bestang2 + angle_step * 6.0f; ang += angle_step) {
00876                                         EMData v(*this_img);
00877                                         Transform t(Dict("type","2d","alpha",static_cast<float>(ang*EMConsts::rad2deg)));
00878                                         t.set_trans(dx,dy);
00879                                         v.transform(t);
00880 //                                      v.rotate_translate(ang*EMConsts::rad2deg, 0.0f, 0.0f, (float)dx, (float)dy, 0.0f);
00881 
00882                                         float lc = v.cmp(cmp_name, to, cmp_params);
00883 
00884                                         if (lc < bestval) {
00885                                                 bestval = lc;
00886                                                 bestang = ang;
00887                                                 bestdx = dx;
00888                                                 bestdy = dy;
00889                                                 bestflip = 0;
00890                                         }
00891 
00892                                         lc = v.cmp(cmp_name, flipped, cmp_params);
00893 
00894                                         if (lc < bestval) {
00895                                                 bestval = lc;
00896                                                 bestang = ang;
00897                                                 bestdx = dx;
00898                                                 bestdy = dy;
00899                                                 bestflip = 1;
00900                                         }
00901                                 }
00902                         }
00903                 }
00904         }
00905 
00906         if (delete_flipped) { delete flipped; flipped = 0; }
00907 
00908         bestang *= (float)EMConsts::rad2deg;
00909         Transform t(Dict("type","2d","alpha",(float)bestang));
00910         t.set_trans(bestdx,bestdy);
00911 
00912         if (bestflip) {
00913                 t.set_mirror(true);
00914         }
00915 
00916         EMData* rslt = this_img->process("math.transform",Dict("transform",&t));
00917         rslt->set_attr("xform.align2d",&t);
00918 
00919         return rslt;
00920 }
00921 
00922 
00923 
00924 static double refalifn(const gsl_vector * v, void *params)
00925 {
00926         Dict *dict = (Dict *) params;
00927 
00928         double x = gsl_vector_get(v, 0);
00929         double y = gsl_vector_get(v, 1);
00930         double a = gsl_vector_get(v, 2);
00931 
00932         EMData *this_img = (*dict)["this"];
00933         EMData *with = (*dict)["with"];
00934         bool mirror = (*dict)["mirror"];
00935 
00936 //      float mean = (float)this_img->get_attr("mean");
00937 //      if ( Util::goodf(&mean) ) {
00938 //              //cout << "tmps mean is nan even before rotation" << endl;
00939 //      }
00940 
00941         Transform t(Dict("type","2d","alpha",static_cast<float>(a)));
00942 //      Transform3D t3d(Transform3D::EMAN, (float)a, 0.0f, 0.0f);
00943 //      t3d.set_posttrans( (float) x, (float) y);
00944 //      tmp->rotate_translate(t3d);
00945         t.set_trans((float)x,(float)y);
00946         t.set_mirror(mirror);
00947         EMData *tmp = this_img->process("math.transform",Dict("transform",&t));
00948 
00949         Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
00950         double result = c->cmp(tmp,with);
00951 
00952         // DELETE AT SOME STAGE, USEFUL FOR PRERELEASE STUFF
00953         //      float test_result = (float)result;
00954 //      if ( Util::goodf(&test_result) ) {
00955 //              cout << "result " << result << " " << x << " " << y << " " << a << endl;
00956 //              cout << (float)this_img->get_attr("mean") << " " << (float)tmp->get_attr("mean") << " " << (float)with->get_attr("mean") << endl;
00957 //              tmp->write_image("tmp.hdf");
00958 //              with->write_image("with.hdf");
00959 //              this_img->write_image("this_img.hdf");
00960 //              EMData* t = this_img->copy();
00961 //              cout << (float)t->get_attr("mean") << endl;
00962 //              t->rotate_translate( t3d );
00963 //              cout << (float)t->get_attr("mean") << endl;
00964 //              cout << "exit" << endl;
00966 //              cout << (float)t->get_attr("mean") << endl;
00967 //              cout << "now exit" << endl;
00968 //              delete t;
00969 //      }
00970 
00971 
00972         if ( tmp != 0 ) delete tmp;
00973 
00974         return result;
00975 }
00976 
00977 static double refalifnfast(const gsl_vector * v, void *params)
00978 {
00979         Dict *dict = (Dict *) params;
00980         EMData *this_img = (*dict)["this"];
00981         EMData *img_to = (*dict)["with"];
00982         bool mirror = (*dict)["mirror"];
00983 
00984         double x = gsl_vector_get(v, 0);
00985         double y = gsl_vector_get(v, 1);
00986         double a = gsl_vector_get(v, 2);
00987 
00988         double r = this_img->dot_rotate_translate(img_to, (float)x, (float)y, (float)a, mirror);
00989         int nsec = this_img->get_xsize() * this_img->get_ysize();
00990         double result = 1.0 - r / nsec;
00991 
00992 //      cout << result << " x " << x << " y " << y << " az " << a <<  endl;
00993         return result;
00994 }
00995 
00996 
00997 EMData *RefineAligner::align(EMData * this_img, EMData *to,
00998         const string & cmp_name, const Dict& cmp_params) const
00999 {
01000 
01001         if (!to) {
01002                 return 0;
01003         }
01004 
01005         int mode = params.set_default("mode", 0);
01006         float saz = 0.0;
01007         float sdx = 0.0;
01008         float sdy = 0.0;
01009         bool mirror = false;
01010         Transform* t;
01011         if (params.has_key("xform.align2d") ) {
01012                 t = params["xform.align2d"];
01013                 Dict params = t->get_params("2d");
01014                 saz = params["alpha"];
01015                 sdx = params["tx"];
01016                 sdy = params["ty"];
01017                 mirror = params["mirror"];
01018 
01019         } else {
01020                 t = new Transform(); // is the identity
01021         }
01022 
01023         int np = 3;
01024         Dict gsl_params;
01025         gsl_params["this"] = this_img;
01026         gsl_params["with"] = to;
01027         gsl_params["snr"]  = params["snr"];
01028         gsl_params["mirror"] = mirror;
01029 
01030 
01031 
01032         const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
01033         gsl_vector *ss = gsl_vector_alloc(np);
01034 
01035         float stepx = params.set_default("stepx",1.0f);
01036         float stepy = params.set_default("stepy",1.0f);
01037         // Default step is 5 degree - note in EMAN1 it was 0.1 radians
01038         float stepaz = params.set_default("stepaz",5.0f);
01039 
01040         gsl_vector_set(ss, 0, stepx);
01041         gsl_vector_set(ss, 1, stepy);
01042         gsl_vector_set(ss, 2, stepaz);
01043 
01044         gsl_vector *x = gsl_vector_alloc(np);
01045         gsl_vector_set(x, 0, sdx);
01046         gsl_vector_set(x, 1, sdy);
01047         gsl_vector_set(x, 2, saz);
01048 
01049         Cmp *c = 0;
01050 
01051         gsl_multimin_function minex_func;
01052         if (mode == 2) {
01053                 minex_func.f = &refalifnfast;
01054         }
01055         else {
01056                 c = Factory < Cmp >::get(cmp_name, cmp_params);
01057                 gsl_params["cmp"] = (void *) c;
01058                 minex_func.f = &refalifn;
01059         }
01060 
01061         minex_func.n = np;
01062         minex_func.params = (void *) &gsl_params;
01063 
01064         gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
01065         gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
01066 
01067         int rval = GSL_CONTINUE;
01068         int status = GSL_SUCCESS;
01069         int iter = 1;
01070 
01071         float precision = params.set_default("precision",0.04f);
01072         int maxiter = params.set_default("maxiter",28);
01073 
01074 //      printf("Refine sx=%1.2f sy=%1.2f sa=%1.2f prec=%1.4f maxit=%d\n",stepx,stepy,stepaz,precision,maxiter);
01075 //      printf("%1.2f %1.2f %1.1f  ->",(float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1),(float)gsl_vector_get(s->x, 2));
01076 
01077         while (rval == GSL_CONTINUE && iter < maxiter) {
01078                 iter++;
01079                 status = gsl_multimin_fminimizer_iterate(s);
01080                 if (status) {
01081                         break;
01082                 }
01083                 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
01084         }
01085 
01086         int maxshift = params.set_default("maxshift",-1);
01087 
01088         if (maxshift <= 0) {
01089                 maxshift = this_img->get_xsize() / 4;
01090         }
01091         float fmaxshift = static_cast<float>(maxshift);
01092         EMData *result;
01093         if ( fmaxshift >= fabs((float)gsl_vector_get(s->x, 0)) && fmaxshift >= fabs((float)gsl_vector_get(s->x, 1))  )
01094         {
01095 //              printf(" Refine good %1.2f %1.2f %1.1f\n",(float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1),(float)gsl_vector_get(s->x, 2));
01096                 Transform  tsoln(Dict("type","2d","alpha",(float)gsl_vector_get(s->x, 2)));
01097                 tsoln.set_mirror(mirror);
01098                 tsoln.set_trans((float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1));
01099                 result = this_img->process("math.transform",Dict("transform",&tsoln));
01100                 result->set_attr("xform.align2d",&tsoln);
01101         } else { // The refine aligner failed - this shift went beyond the max shift
01102 //              printf(" Refine Failed %1.2f %1.2f %1.1f\n",(float)gsl_vector_get(s->x, 0),(float)gsl_vector_get(s->x, 1),(float)gsl_vector_get(s->x, 2));
01103                 result = this_img->process("math.transform",Dict("transform",t));
01104                 result->set_attr("xform.align2d",t);
01105         }
01106 
01107         delete t;
01108         t = 0;
01109 
01110         gsl_vector_free(x);
01111         gsl_vector_free(ss);
01112         gsl_multimin_fminimizer_free(s);
01113 
01114         if ( c != 0 ) delete c;
01115         return result;
01116 }
01117 
01118 static Transform refalin3d_perturb(const Transform*const t, const float& delta, const float& arc, const float& phi, const float& x, const float& y, const float& z)
01119 {
01120         Dict orig_params = t->get_params("eman");
01121         float orig_phi = orig_params["phi"];
01122         float orig_x = orig_params["tx"];
01123         float orig_y = orig_params["ty"];
01124         float orig_z = orig_params["tz"];
01125         orig_params["phi"] = 0;
01126         orig_params["tx"] = 0;
01127         orig_params["ty"] = 0;
01128         orig_params["tz"] = 0;
01129         Transform t_no_phi(orig_params);
01130 
01131         Vec3f zz(0,0,1);
01132 
01133         Vec3f vv  = t_no_phi.transpose()*zz;
01134         Vec3f normal = Vec3f(-vv[2],0,-vv[0]);
01135 
01136         normal.normalize();
01137 
01138         Dict d;
01139         d["type"] = "spin";
01140         d["Omega"] = arc;
01141         d["n1"] = vv[0];
01142         d["n2"] = vv[1];
01143         d["n3"] = vv[2];
01144 
01145         Transform q(d);
01146 
01147         Vec3f  rot_vec = q*normal;
01148         rot_vec.normalize();
01149 
01150         Dict e;
01151         e["type"] = "spin";
01152         e["Omega"] = delta;
01153         e["n1"] = rot_vec[0];
01154         e["n2"] = rot_vec[1];
01155         e["n3"] = rot_vec[2];
01156 
01157         Transform perturb(e);
01158 
01159         Dict g;
01160         g["type"] = "eman";
01161         g["alt"] = 0;
01162         g["az"] = 0;
01163         g["phi"] = 0*phi+orig_phi;
01164 
01165         Transform phi_rot(g);
01166         Transform soln = t_no_phi*perturb*phi_rot;
01167         soln.set_trans(x+orig_x,y+orig_y,z+orig_z);
01168 
01169         Dict params = soln.get_params("eman");
01170         return soln;
01171 }
01172 
01173 static double refalifn3d(const gsl_vector * v, void *params)
01174 {
01175         Dict *dict = (Dict *) params;
01176         double x = gsl_vector_get(v, 0);
01177         double y = gsl_vector_get(v, 1);
01178         double z = gsl_vector_get(v, 2);
01179         double arc = gsl_vector_get(v, 3);
01180         double delta = gsl_vector_get(v, 4);
01181         double phi = gsl_vector_get(v, 5);
01182         EMData *this_img = (*dict)["this"];
01183         EMData *with = (*dict)["with"];
01184 //      bool mirror = (*dict)["mirror"];
01185 
01186         Transform* t = (*dict)["transform"];
01187 
01188         Transform soln = refalin3d_perturb(t,(float)delta,(float)arc,(float)phi,(float)x,(float)y,(float)z);
01189 
01190         EMData *tmp = this_img->process("math.transform",Dict("transform",&soln));
01191         Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
01192         double result = c->cmp(tmp,with);
01193         if ( tmp != 0 ) delete tmp;
01194         delete t; t = 0;
01195 //      cout << result << " " << az << " " << alt << " " << phi << " " << x << " " << y << " " << z << endl;
01196         return result;
01197 }
01198 
01199 /*
01200 static double refalifn3d(const gsl_vector * v, void *params)
01201 {
01202         Dict *dict = (Dict *) params;
01203 
01204         double x = gsl_vector_get(v, 0);
01205         double y = gsl_vector_get(v, 1);
01206         double z = gsl_vector_get(v, 2);
01207         double az = gsl_vector_get(v, 3);
01208         double alt = gsl_vector_get(v, 4);
01209         double phi = gsl_vector_get(v, 5);
01210         EMData *this_img = (*dict)["this"];
01211         EMData *with = (*dict)["with"];
01212         bool mirror = (*dict)["mirror"];
01213 
01214         Dict d("type","eman","az",static_cast<float>(az));
01215         d["alt"] = static_cast<float>(alt);
01216         d["phi"] = static_cast<float>(phi);
01217         Transform t(d);
01218         t.set_trans((float)x,(float)y,(float)z);
01219         t.set_mirror(mirror);
01220         EMData *tmp = this_img->process("math.transform",Dict("transform",&t));
01221 
01222         Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
01223         double result = c->cmp(tmp,with);
01224         if ( tmp != 0 ) delete tmp;
01225 
01226         return result;
01227 }*/
01228 
01229 EMData* Refine3DAligner::align(EMData * this_img, EMData *to,
01230         const string & cmp_name, const Dict& cmp_params) const
01231 {
01232 
01233         if (!to || !this_img) throw NullPointerException("Input image is null"); // not sure if this is necessary, it was there before I started
01234 
01235         if (to->get_ndim() != 3 || this_img->get_ndim() != 3) throw ImageDimensionException("The Refine3D aligner only works for 3D images");
01236 
01237         float saz = 0.0;
01238         float sphi = 0.0;
01239         float salt = 0.0;
01240         float sdx = 0.0;
01241         float sdy = 0.0;
01242         float sdz = 0.0;
01243         bool mirror = false;
01244         Transform* t;
01245         if (params.has_key("xform.align3d") ) {
01246                 // Unlike the 2d refine aligner, this class doesn't require the starting transform's
01247                 // parameters to form the starting guess. Instead the Transform itself
01248                 // is perturbed carefully (using quaternion rotation) to overcome problems that arise
01249                 // when you use orthogonally-based Euler angles
01250                 t = params["xform.align3d"];
01251         }else {
01252                 t = new Transform(); // is the identity
01253         }
01254 
01255         int np = 6; // the number of dimensions
01256         Dict gsl_params;
01257         gsl_params["this"] = this_img;
01258         gsl_params["with"] = to;
01259         gsl_params["snr"]  = params["snr"];
01260         gsl_params["mirror"] = mirror;
01261         gsl_params["transform"] = t;
01262 
01263         const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
01264         gsl_vector *ss = gsl_vector_alloc(np);
01265 
01266         float stepx = params.set_default("stepx",1.0f);
01267         float stepy = params.set_default("stepy",1.0f);
01268         float stepz = params.set_default("stepz",1.0f);
01269         // Default step is 5 degree - note in EMAN1 it was 0.1 radians
01270         float half_circle_step = 180.0f; // This shouldn't be altered
01271         float stepphi = params.set_default("stephi",5.0f);
01272         float stepdelta = params.set_default("stepdelta",5.0f);
01273 
01274         gsl_vector_set(ss, 0, stepx);
01275         gsl_vector_set(ss, 1, stepy);
01276         gsl_vector_set(ss, 2, stepz);
01277         gsl_vector_set(ss, 3, half_circle_step);
01278         gsl_vector_set(ss, 4, stepdelta);
01279         gsl_vector_set(ss, 5, stepphi);
01280 
01281         gsl_vector *x = gsl_vector_alloc(np);
01282         gsl_vector_set(x, 0, sdx);
01283         gsl_vector_set(x, 1, sdy);
01284         gsl_vector_set(x, 2, sdz);
01285         gsl_vector_set(x, 3, saz);
01286         gsl_vector_set(x, 4, salt);
01287         gsl_vector_set(x, 5, sphi);
01288 
01289         gsl_multimin_function minex_func;
01290         Cmp *c = Factory < Cmp >::get(cmp_name, cmp_params);
01291         gsl_params["cmp"] = (void *) c;
01292         minex_func.f = &refalifn3d;
01293 
01294         minex_func.n = np;
01295         minex_func.params = (void *) &gsl_params;
01296 
01297         gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, np);
01298         gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
01299 
01300         int rval = GSL_CONTINUE;
01301         int status = GSL_SUCCESS;
01302         int iter = 1;
01303 
01304         float precision = params.set_default("precision",0.04f);
01305         int maxiter = params.set_default("maxiter",60);
01306         while (rval == GSL_CONTINUE && iter < maxiter) {
01307                 iter++;
01308                 status = gsl_multimin_fminimizer_iterate(s);
01309                 if (status) {
01310                         break;
01311                 }
01312                 rval = gsl_multimin_test_size(gsl_multimin_fminimizer_size(s), precision);
01313         }
01314 
01315         int maxshift = params.set_default("maxshift",-1);
01316 
01317         if (maxshift <= 0) {
01318                 maxshift = this_img->get_xsize() / 4;
01319         }
01320         float fmaxshift = static_cast<float>(maxshift);
01321         EMData *result;
01322         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))
01323         {
01324                 float x = (float)gsl_vector_get(s->x, 0);
01325                 float y = (float)gsl_vector_get(s->x, 1);
01326                 float z = (float)gsl_vector_get(s->x, 2);
01327 
01328                 float arc = (float)gsl_vector_get(s->x, 3);
01329                 float delta = (float)gsl_vector_get(s->x, 4);
01330                 float phi = (float)gsl_vector_get(s->x, 5);
01331 
01332                 Transform tsoln = refalin3d_perturb(t,delta,arc,phi,x,y,z);
01333 
01334                 result = this_img->process("math.transform",Dict("transform",&tsoln));
01335                 result->set_attr("xform.align3d",&tsoln);
01336 
01337         } else { // The refine aligner failed - this shift went beyond the max shift
01338                 result = this_img->process("math.transform",Dict("transform",t));
01339                 result->set_attr("xform.align3d",t);
01340         }
01341 
01342         delete t;
01343         t = 0;
01344 
01345         gsl_vector_free(x);
01346         gsl_vector_free(ss);
01347         gsl_multimin_fminimizer_free(s);
01348 
01349         if ( c != 0 ) delete c;
01350         return result;
01351 }
01352 
01353 EMData* RT3DGridAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
01354 {
01355 
01356         vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
01357 
01358         Dict t;
01359         Transform* tr = (Transform*) alis[0]["xform.align3d"];
01360         t["transform"] = tr;
01361         EMData* soln = this_img->process("math.transform",t);
01362         soln->set_attr("xform.align3d",tr);
01363         delete tr; tr = 0;
01364 
01365         return soln;
01366 
01367 }
01368 
01369 vector<Dict> RT3DGridAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const {
01370 
01371         if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
01372                 throw ImageDimensionException("This aligner only works for 3D images");
01373         }
01374 
01375         int searchx = 0;
01376         int searchy = 0;
01377         int searchz = 0;
01378 
01379         if (params.has_key("search")) {
01380                 vector<string> check;
01381                 check.push_back("searchx");
01382                 check.push_back("searchy");
01383                 check.push_back("searchz");
01384                 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
01385                         if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
01386                 }
01387                 int search  = params["search"];
01388                 searchx = search;
01389                 searchy = search;
01390                 searchz = search;
01391         } else {
01392                 searchx = params.set_default("searchx",3);
01393                 searchy = params.set_default("searchy",3);
01394                 searchz = params.set_default("searchz",3);
01395         }
01396 
01397         float ralt = params.set_default("ralt",180.f);
01398         float rphi = params.set_default("rphi",180.f);
01399         float raz = params.set_default("raz",180.f);
01400         float dalt = params.set_default("dalt",10.f);
01401         float daz = params.set_default("daz",10.f);
01402         float dphi = params.set_default("dphi",10.f);
01403         float threshold = params.set_default("threshold",0.f);
01404         if (threshold < 0.0f) throw InvalidParameterException("The threshold parameter must be greater than or equal to zero");
01405         bool verbose = params.set_default("verbose",false);
01406 
01407         vector<Dict> solns;
01408         if (nsoln == 0) return solns; // What was the user thinking?
01409         for (unsigned int i = 0; i < nsoln; ++i ) {
01410                 Dict d;
01411                 d["score"] = 1.e24;
01412                 Transform t; // identity by default
01413                 d["xform.align3d"] = &t; // deep copy is going on here
01414                 solns.push_back(d);
01415         }
01416         Dict d;
01417         d["type"] = "eman"; // d is used in the loop below
01418         for ( float alt = 0.0f; alt <= ralt; alt += dalt) {
01419                 // An optimization for the range of az is made at the top of the sphere
01420                 // If you think about it, this is just a coarse way of making this approach slightly more efficient
01421                 if (verbose) {
01422                         cout << "Trying angle " << alt << endl;
01423                 }
01424 
01425                 float begin_az = -raz;
01426                 float end_az = raz;
01427                 if (alt == 0.0f) {
01428                         begin_az = 0.0f;
01429                         end_az = 0.0f;
01430                 }
01431 
01432                 for ( float az = begin_az; az <= end_az; az += daz ){
01433                         for( float phi = -rphi-az; phi <= rphi-az; phi += dphi ) {
01434                                 d["alt"] = alt;
01435                                 d["phi"] = phi;
01436                                 d["az"] = az;
01437                                 Transform t(d);
01438                                 EMData* transformed = this_img->process("math.transform",Dict("transform",&t));
01439                                 EMData* ccf = transformed->calc_ccf(to);
01440 
01441                                 IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
01442                                 Dict altered_cmp_params(cmp_params);
01443                                 if (cmp_name == "dot.tomo") {
01444                                         altered_cmp_params["ccf"] = ccf;
01445                                         altered_cmp_params["tx"] = point[0];
01446                                         altered_cmp_params["ty"] = point[1];
01447                                         altered_cmp_params["tz"] = point[2];
01448                                 }
01449 
01450                                 float best_score = transformed->cmp(cmp_name,to,altered_cmp_params);
01451                                 delete transformed; transformed =0;
01452                                 delete ccf; ccf = 0;
01453 
01454                                 unsigned int j = 0;
01455                                 for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
01456                                         if ( (float)(*it)["score"] > best_score ) {  // Note greater than - EMAN2 preferes minimums as a matter of policy
01457                                                 vector<Dict>::reverse_iterator rit = solns.rbegin();
01458                                                 copy(rit+1,solns.rend()-j,rit);
01459                                                 Dict& d = (*it);
01460                                                 d["score"] = best_score;
01461                                                 d["xform.align3d"] = &t;
01462                                                 break;
01463                                         }
01464                                 }
01465                         }
01466                 }
01467         }
01468 
01469         return solns;
01470 
01471 }
01472 
01473 EMData* RT3DSphereAligner::align(EMData * this_img, EMData *to, const string & cmp_name, const Dict& cmp_params) const
01474 {
01475 
01476         vector<Dict> alis = xform_align_nbest(this_img,to,1,cmp_name,cmp_params);
01477 
01478         Dict t;
01479         Transform* tr = (Transform*) alis[0]["xform.align3d"];
01480         t["transform"] = tr;
01481         EMData* soln = this_img->process("math.transform",t);
01482         soln->set_attr("xform.align3d",tr);
01483         delete tr; tr = 0;
01484 
01485         return soln;
01486 
01487 }
01488 
01489 vector<Dict> RT3DSphereAligner::xform_align_nbest(EMData * this_img, EMData * to, const unsigned int nsoln, const string & cmp_name, const Dict& cmp_params) const {
01490 
01491         if ( this_img->get_ndim() != 3 || to->get_ndim() != 3 ) {
01492                 throw ImageDimensionException("This aligner only works for 3D images");
01493         }
01494 
01495         int searchx = 0;
01496         int searchy = 0;
01497         int searchz = 0;
01498 
01499         if (params.has_key("search")) {
01500                 vector<string> check;
01501                 check.push_back("searchx");
01502                 check.push_back("searchy");
01503                 check.push_back("searchz");
01504                 for(vector<string>::const_iterator cit = check.begin(); cit != check.end(); ++cit) {
01505                         if (params.has_key(*cit)) throw InvalidParameterException("The search parameter is mutually exclusive of the searchx, searchy, and searchz parameters");
01506                 }
01507                 int search  = params["search"];
01508                 searchx = search;
01509                 searchy = search;
01510                 searchz = search;
01511         } else {
01512                 searchx = params.set_default("searchx",3);
01513                 searchy = params.set_default("searchy",3);
01514                 searchz = params.set_default("searchz",3);
01515         }
01516 
01517         float rphi = params.set_default("rphi",180.f);
01518         float dphi = params.set_default("dphi",10.f);
01519         float threshold = params.set_default("threshold",0.f);
01520         if (threshold < 0.0f) throw InvalidParameterException("The threshold parameter must be greater than or equal to zero");
01521         bool verbose = params.set_default("verbose",false);
01522 
01523         vector<Dict> solns;
01524         if (nsoln == 0) return solns; // What was the user thinking?
01525         for (unsigned int i = 0; i < nsoln; ++i ) {
01526                 Dict d;
01527                 d["score"] = 1.e24;
01528                 Transform t; // identity by default
01529                 d["xform.align3d"] = &t; // deep copy is going on here
01530                 solns.push_back(d);
01531         }
01532 
01533         Dict d;
01534         d["inc_mirror"] = true; // This should probably always be true for 3D case. If it ever changes we might have to make inc_mirror a parameter
01535         if ( params.has_key("delta") && params.has_key("n") ) {
01536                 throw InvalidParameterException("The delta and n parameters are mutually exclusive in the RT3DSphereAligner aligner");
01537         } else if ( params.has_key("n") ) {
01538                 d["n"] = params["n"];
01539         } else {
01540                 // If they didn't specify either then grab the default delta - if they did supply delta we're still safe doing this
01541                 d["delta"] = params.set_default("delta",10.f);
01542         }
01543 
01544         Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
01545         vector<Transform> transforms = sym->gen_orientations((string)params.set_default("orientgen","eman"),d);
01546 
01547         float verbose_alt = -1.0f;;
01548         for(vector<Transform>::const_iterator trans_it = transforms.begin(); trans_it != transforms.end(); trans_it++) {
01549                 Dict params = trans_it->get_params("eman");
01550                 float az = params["az"];
01551                 if (verbose) {
01552                         float alt = params["alt"];
01553                         if ( alt != verbose_alt ) {
01554                                 verbose_alt = alt;
01555                                 cout << "Trying angle " << alt << endl;
01556                         }
01557                 }
01558                 for( float phi = -rphi-az; phi <= rphi-az; phi += dphi ) {
01559                         params["phi"] = phi;
01560                         Transform t(params);
01561                         EMData* transformed = this_img->process("math.transform",Dict("transform",&t));
01562                         EMData* ccf = transformed->calc_ccf(to);
01563                         IntPoint point = ccf->calc_max_location_wrap(searchx,searchy,searchz);
01564                         Dict altered_cmp_params(cmp_params);
01565                         if (cmp_name == "dot.tomo") {
01566                                 altered_cmp_params["ccf"] = ccf;
01567                                 altered_cmp_params["tx"] = point[0];
01568                                 altered_cmp_params["ty"] = point[1];
01569                                 altered_cmp_params["tz"] = point[2];
01570                         }
01571 
01572                         float best_score = transformed->cmp(cmp_name,to,altered_cmp_params);
01573                         delete transformed; transformed =0;
01574                         delete ccf; ccf =0;
01575 
01576                         unsigned int j = 0;
01577                         for ( vector<Dict>::iterator it = solns.begin(); it != solns.end(); ++it, ++j ) {
01578                                 if ( (float)(*it)["score"] > best_score ) { // Note greater than - EMAN2 preferes minimums as a matter of policy
01579                                         vector<Dict>::reverse_iterator rit = solns.rbegin();
01580                                         copy(rit+1,solns.rend()-j,rit);
01581                                         Dict& d = (*it);
01582                                         d["score"] = best_score;
01583                                         d["xform.align3d"] = &t; // deep copy is going on here
01584                                         break;
01585                                 }
01586                         }
01587 
01588                 }
01589         }
01590         delete sym; sym = 0;
01591         return solns;
01592 
01593 }
01594 
01595 CUDA_Aligner::CUDA_Aligner() {
01596         image_stack = NULL;
01597         image_stack_filtered = NULL;
01598         ccf = NULL;
01599 }
01600 
01601 #ifdef EMAN2_USING_CUDA
01602 void CUDA_Aligner::finish() {
01603         if (image_stack) delete image_stack;
01604         if (image_stack_filtered) delete image_stack_filtered;
01605         if (ccf) delete ccf;
01606 }
01607 
01608 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) {
01609 
01610         NIMA = nima;
01611         NX = nx;
01612         NY = ny;
01613         RING_LENGTH = ring_length;
01614         NRING = nring;
01615         STEP = step;
01616         KX = kx;
01617         KY = ky;
01618         OU = ou;
01619         CTF = ctf;
01620         
01621         image_stack = (float *)malloc(NIMA*NX*NY*sizeof(float));
01622         if (CTF == 1) image_stack_filtered = (float *)malloc(NIMA*NX*NY*sizeof(float));
01623         ccf = (float *)malloc(2*(2*KX+1)*(2*KY+1)*NIMA*(RING_LENGTH+2)*sizeof(float));
01624 }
01625 
01626 void CUDA_Aligner::insert_image(EMData *image, int num) {
01627 
01628         int base_address = num*NX*NY;
01629 
01630         for (int y=0; y<NY; y++)
01631                 for (int x=0; x<NX; x++)
01632                         image_stack[base_address+y*NX+x] = (*image)(x, y);
01633 }
01634 
01635 void CUDA_Aligner::filter_stack(vector<float> ctf_params, int id) {
01636         
01637         float *params;
01638         
01639         params = (float *)malloc(NIMA*6*sizeof(float)); 
01640         
01641         for (int i=0; i<NIMA*6; i++) params[i] = ctf_params[i];
01642 
01643         filter_image(image_stack, image_stack_filtered, NIMA, NX, NY, params, id);
01644 
01645         free(params);
01646 }
01647 
01648 void CUDA_Aligner::sum_oe(vector<float> ctf_params, vector<float> ali_params, EMData *ave1, EMData *ave2, int id) {
01649         
01650         float *ctf_p, *ali_p, *av1, *av2;
01651         
01652         ctf_p = (float *)malloc(NIMA*6*sizeof(float));
01653         ali_p = (float *)malloc(NIMA*4*sizeof(float));
01654         
01655         if (CTF == 1) {
01656                 for (int i=0; i<NIMA*6; i++)  ctf_p[i] = ctf_params[i];
01657         }
01658         for (int i=0; i<NIMA*4; i++)   ali_p[i] = ali_params[i];
01659         
01660         av1 = ave1->get_data();
01661         av2 = ave2->get_data();
01662         
01663         rot_filt_sum(image_stack, NIMA, NX, NY, CTF, ctf_p, ali_p, av1, av2, id);
01664         
01665         free(ctf_p);
01666         free(ali_p);
01667 }
01668 
01669 vector<float> CUDA_Aligner::alignment_2d(EMData *ref_image_em, vector<float> sx_list, vector<float> sy_list, int id, int silent) {
01670 
01671         float *ref_image, max_ccf;
01672         int base_address, ccf_offset;
01673         float ts, tm;
01674         float ang, sx, sy, mirror, co, so, sxs, sys;
01675         float *sx2, *sy2;
01676         vector<float> align_result;
01677 
01678         sx2 = (float *)malloc(NIMA*sizeof(float));
01679         sy2 = (float *)malloc(NIMA*sizeof(float));
01680 
01681         ref_image = ref_image_em->get_data();
01682         
01683         for (int i=0; i<NIMA; i++) {
01684                 sx2[i] = sx_list[i];
01685                 sy2[i] = sy_list[i];
01686         }
01687         
01688         if (CTF == 1) {
01689                 calculate_ccf(image_stack_filtered, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, id, silent);
01690         } else {
01691                 calculate_ccf(image_stack, ref_image, ccf, NIMA, NX, NY, RING_LENGTH, NRING, OU, STEP, KX, KY, sx2, sy2, id, silent);
01692         }
01693 
01694         ccf_offset = NIMA*(RING_LENGTH+2)*(2*KX+1)*(2*KY+1);
01695 
01696         for (int im=0; im<NIMA; im++) {
01697                 max_ccf = -1.0e22;
01698                 for (int kx=-KX; kx<=KX; kx++) {
01699                         for (int ky=-KY; ky<=KY; ky++) {
01700                                 base_address = (((ky+KY)*(2*KX+1)+(kx+KX))*NIMA+im)*(RING_LENGTH+2);
01701                                 for (int l=0; l<RING_LENGTH; l++) {
01702                                         ts = ccf[base_address+l];
01703                                         tm = ccf[base_address+l+ccf_offset];
01704                                         if (ts > max_ccf) {
01705                                                 ang = float(l)/RING_LENGTH*360.0;
01706                                                 sx = -kx*STEP;
01707                                                 sy = -ky*STEP;
01708                                                 mirror = 0;
01709                                                 max_ccf = ts;
01710                                         }
01711                                         if (tm > max_ccf) {
01712                                                 ang = float(l)/RING_LENGTH*360.0; 
01713                                                 sx = -kx*STEP;
01714                                                 sy = -ky*STEP;
01715                                                 mirror = 1;
01716                                                 max_ccf = tm;
01717                                         }
01718                                 }
01719                         }
01720                 }
01721                 co =  cos(ang*M_PI/180.0);
01722                 so = -sin(ang*M_PI/180.0);
01723                 sxs = sx*co - sy*so;
01724                 sys = sx*so + sy*co;
01725 
01726                 align_result.push_back(ang);
01727                 align_result.push_back(sxs);
01728                 align_result.push_back(sys);
01729                 align_result.push_back(mirror);
01730         }
01731         
01732         free(sx2);
01733         free(sy2);
01734         
01735         return align_result;
01736 }
01737 #endif
01738 
01739 void EMAN::dump_aligners()
01740 {
01741         dump_factory < Aligner > ();
01742 }
01743 
01744 map<string, vector<string> > EMAN::dump_aligners_list()
01745 {
01746         return dump_factory_list < Aligner > ();
01747 }

Generated on Sat Nov 21 02:19:14 2009 for EMAN2 by  doxygen 1.5.6