Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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

Generated on Tue Jul 12 13:49:24 2011 for EMAN2 by  doxygen 1.3.9.1