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

Generated on Mon Mar 7 18:15:24 2011 for EMAN2 by  doxygen 1.3.9.1