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