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

Generated on Thu Nov 17 12:43:44 2011 for EMAN2 by  doxygen 1.3.9.1