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

Generated on Thu May 3 10:06:23 2012 for EMAN2 by  doxygen 1.4.7