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

aligner.cpp

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

Generated on Fri Aug 10 16:35:22 2012 for EMAN2 by  doxygen 1.3.9.1