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

aligner.cpp File Reference

#include "emfft.h"
#include "cmp.h"
#include "aligner.h"
#include "averager.h"
#include "emdata.h"
#include "processor.h"
#include "util.h"
#include "symmetry.h"
#include <gsl/gsl_multimin.h>
#include "plugins/aligner_template.h"

Include dependency graph for aligner.cpp:

Include dependency graph

Go to the source code of this file.

Defines

#define EMAN2_ALIGNER_DEBUG   0
 
Id
aligner.cpp,v 1.272 2012/08/06 00:45:23 stevel Exp


Functions

double refalifn (const gsl_vector *v, void *params)
void refalidf (const gsl_vector *v, void *params, gsl_vector *df)
void refalifdf (const gsl_vector *v, void *params, double *f, gsl_vector *df)
double refalifnfast (const gsl_vector *v, void *params)
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)
double symquat (const gsl_vector *v, void *params)
double refalifn3dquat (const gsl_vector *v, void *params)
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)


Define Documentation

#define EMAN2_ALIGNER_DEBUG   0
 

Id
aligner.cpp,v 1.272 2012/08/06 00:45:23 stevel Exp

Definition at line 55 of file aligner.cpp.


Function Documentation

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
[static]
 

Definition at line 2871 of file aligner.cpp.

References EMAN::EMData::cmp(), EMAN::EMData::copy(), EMAN::EMData::do_ift(), EMAN::EMData::get_data(), EMAN::EMData::get_ysize(), in, EMAN::EMData::rotate(), EMAN::EMData::set_attr(), EMAN::EMData::set_complex(), EMAN::EMData::set_fftodd(), EMAN::EMData::set_ri(), EMAN::EMData::set_size(), EMAN::Transform::set_trans(), and EMAN::EMData::translate().

Referenced by EMAN::FRM2DAligner::align().

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

void refalidf const gsl_vector *  v,
void *  params,
gsl_vector *  df
[static]
 

Definition at line 1601 of file aligner.cpp.

References refalifn(), v, and vp.

01601                                                                          {
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 }

void refalifdf const gsl_vector *  v,
void *  params,
double *  f,
gsl_vector *  df
[static]
 

Definition at line 1624 of file aligner.cpp.

References refalifn(), v, and vp.

01624                                                                                        {
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 }

double refalifn const gsl_vector *  v,
void *  params
[static]
 

Definition at line 1569 of file aligner.cpp.

References EMAN::Cmp::cmp(), EMAN::Dict::has_key(), EMAN::EMData::mult(), EMAN::EMData::process(), EMAN::Transform::set_mirror(), EMAN::Transform::set_scale(), EMAN::Transform::set_trans(), t, v, x, and y.

Referenced by refalidf(), and refalifdf().

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 }

double refalifn3dquat const gsl_vector *  v,
void *  params
[static]
 

Definition at line 1980 of file aligner.cpp.

References EMAN::Cmp::cmp(), EMAN::EMData::process(), refalin3d_perturbquat(), t, v, x, and y.

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 }

double refalifnfast const gsl_vector *  v,
void *  params
[static]
 

Definition at line 1648 of file aligner.cpp.

References EMAN::EMData::dot_rotate_translate(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), v, x, and y.

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 }

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
[static]
 

Definition at line 1930 of file aligner.cpp.

References EMAN::Vec3< Type >::normalize(), q, EMAN::Transform::set_trans(), sqrt(), EMAN::Vec3f, x, and y.

Referenced by EMAN::Refine3DAlignerQuaternion::align(), EMAN::SymAlignProcessorQuat::align(), refalifn3dquat(), and symquat().

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 }

double symquat const gsl_vector *  v,
void *  params
[static]
 

Definition at line 1952 of file aligner.cpp.

References EMAN::Cmp::cmp(), EMAN::EMData::process(), refalin3d_perturbquat(), t, v, x, and y.

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 }


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