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.274 2013/02/06 21:02:20 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.274 2013/02/06 21:02:20 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 2887 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().

02888 {
02889         int size=rsize;
02890         float dx,dy;
02891         int bw=size/2;
02892         int MAXR=this_img->get_ysize()/2;
02893         //int MAXR=size;
02894         unsigned long tsize=2*size;
02895         unsigned long ind1=0, ind2=0, ind3=0, ind4=0, ind41=0;
02896         unsigned long index0=0;
02897         int i=0, j=0, m=0, n=0, r=0;
02898         int loop_rho=0, rho_best=0;
02899 
02900         float* gnr2   = new float[size*2];
02901         float* maxcor = new float[size+1];                  // MAXR need change
02902 
02903         int p_max=p_max_input;
02904         float* result = new float[5*(p_max+1)];
02905         float* cr=new float[size*(bw+1)];
02906         float* ci=new float[size*(bw+1)];
02907         EMData *data_in=new EMData;
02908         data_in->set_complex(true);
02909         data_in->set_fftodd(false);
02910         data_in->set_ri(true);
02911         data_in->set_size(size+2,size,1);
02912         float *in=data_in->get_data();
02913 
02914         float *self_sampl_fft = selfpcsfft->get_data(); // ming f(r)
02915 
02916         float maxcor_sofar=0.0f;
02917         int p=0;
02918 
02919         for(p=0; p<=p_max; ++p){
02920                 ind1=p*size*bw;
02921                 for (i=0;i<size;++i)
02922                         for (j=0;j<bw+1;++j){
02923                                 cr[i*(bw+1)+j]=0.0;
02924                                 ci[i*(bw+1)+j]=0.0;
02925                         }
02926         for(n=0;n<bw;++n){                                // loop for n
02927                 ind2=(ind1+n);
02928                 index0=n*(bw+1);
02929                         for(r=0;r<=MAXR;++r) {
02930                         ind3=(ind2+r*bw)*size;
02931                         for(m=0;m<size;m++){              // take back hat{h(n,r,p)}(m)
02932                                 ind4=(ind3+m)*2;
02933                                     ind41=ind4+1;
02934                                     gnr2[2*m]=frm2dhhat[ind4];
02935                                     gnr2[2*m+1]=frm2dhhat[ind41];
02936                                 }
02937                         for(m=0;m<bw;++m){
02938                                         float tempr=self_sampl_fft[2*m+r*(size+2)]*r;
02939                                 float tempi=self_sampl_fft[2*m+1+r*(size+2)]*r;
02940                                 float gnr2_r=gnr2[2*m];
02941                                 float gnr2_i=gnr2[2*m+1];
02942                                 cr[n*(bw+1)+m]+=gnr2_r*tempr+gnr2_i*tempi;
02943                                         ci[n*(bw+1)+m]+=gnr2_i*tempr-gnr2_r*tempi;
02944                                         if(n!=0){                                       // m,-n
02945                                         if(m!= 0){
02946                                                 int ssize=tsize-2*m;    // ssize = 2*size-2m
02947                                                 int ssize1=ssize+1;
02948                                                 float gnr2_r=gnr2[ssize];
02949                                                 float gnr2_i=gnr2[ssize1];
02950                                                         cr[(size-n)*(bw+1)+m]+=gnr2_r*tempr-gnr2_i*tempi;
02951                                                 ci[(size-n)*(bw+1)+m]-=gnr2_i*tempr+gnr2_r*tempi;
02952                                         }
02953                                                 else{
02954                                                         cr[(size-n)*(bw+1)+m]+=*(gnr2)*tempr-*(gnr2+1)*tempi;
02955                                                         ci[(size-n)*(bw+1)+m]-=*(gnr2+1)*tempr+*(gnr2)*tempi;
02956                                                 }
02957                                 }
02958                                 }
02959                         }
02960         }
02961         for (int cii=0; cii<size*(bw+1);++cii){
02962                         in[2*cii]=cr[cii];
02963                         in[2*cii+1]=ci[cii];
02964                         //printf("cii=%d,in[2i+1]=%f\n",cii, cr[cii]);
02965         }
02966 
02967         EMData *data_out;
02968                 data_out=data_in->do_ift();
02969                 float *c=data_out->get_data();
02970                 float tempr=0.0f, corre_fcs=999.0f;
02971 
02972             int n_best=0, m_best=0;
02973         float temp=-100.0f;
02974                 for(n=0;n<size;++n){// move Tri_2D to Tri = c(phi,phi';rho)
02975                         for(m=0;m<size;++m){
02976                                 temp=c[n*size+m];
02977                                 if(temp>tempr) {
02978                                         tempr=temp;
02979                                         n_best=n;
02980                                         m_best=m;
02981                                 }
02982                         }
02983                 }
02984                 delete data_out;
02985 
02986                 float corre,Phi2,Phi,Tx,Ty,Vx, Vy;
02987 
02988                 //for (n_best=0;n_best<bw;n_best++)
02989                   //  for (m_best=0;m_best<2*bw;m_best++){
02990                 //n_best=0;
02991                 //m_best=70;
02992                 Phi2=n_best*M_PI/bw;  // ming this is reference image rotation angle
02993                 Phi=m_best*M_PI/bw;   // ming this is particle image rotation angle
02994                 Vx=p*cos(Phi);//should use the angle of the centered one
02995                 Vy=-p*sin(Phi);
02996                 Tx=Vx+(floor(com_this_x+0.5f)-floor(com_with_x+0.5f));
02997                 Ty=Vy+(floor(com_this_y+0.5f)-floor(com_with_y+0.5f));
02998 
02999                 dx=-Tx; // the Rota & Trans value (Tx,Ty, ang_keep) are for the projection image,
03000                 dy=-Ty; // need to convert to raw image
03001 
03002                 EMData *this_tmp=this_img->copy();//ming change to to
03003                 this_tmp->rotate(-(Phi2-Phi)*180.0f,0.0f,0.0f);
03004                 this_tmp->translate(dx,dy,0.0);
03005 
03006                 corre=this_tmp->cmp(cmp_name,to,cmp_params);
03007                 //printf("corre=%f\n",corre);
03008                 delete this_tmp;
03009                 if(corre<=corre_fcs) { //ming, cmp use smaller value stands for more similarity
03010                         corre_fcs=corre;
03011                         result[0+5*p] = float(p);       // rho
03012                         result[1+5*p] = corre;          // correlation_fcs
03013                         result[2+5*p] = (Phi2-Phi)*180.0f;      // rotation angle
03014                         result[3+5*p] = Tx;             // Translation_x
03015                         result[4+5*p] = Ty;             // Translation_y
03016                 }
03017                 maxcor[p]=corre_fcs;                            //  maximum correlation for current rho
03018                 if(corre_fcs<maxcor_sofar) {
03019                         maxcor_sofar=corre_fcs;                 // max correlation up to current rho
03020                     rho_best=p;                         // the rho value with maxinum correlation value
03021                 }
03022                 if(p>=4){
03023                         if(maxcor[p] < maxcor[p-1] && maxcor[p-1] < maxcor[p-2]&& maxcor[p-2] < maxcor[p-3] && maxcor[p-3] < maxcor[p-4]){
03024                                 loop_rho=1;
03025                                 break; //exit p loop
03026                         }
03027                 }
03028         } // end for p
03029         //}//test my method
03030         if(loop_rho == 1)
03031                 p=p+1;
03032         int rb5=5*rho_best;
03033         float fsc      = result[1+rb5];
03034         float ang_keep = result[2+rb5];
03035         float Tx       = result[3+rb5];
03036         float Ty       = result[4+rb5];
03037         delete[] gnr2;
03038         delete[] maxcor;
03039         delete[] result;
03040         delete cr;
03041         cr=0;
03042         delete ci;
03043         ci=0;
03044         delete data_in; // ming add
03045         dx = -Tx;               // the Rota & Trans value (Tx,Ty, ang_keep) are for the projection image,
03046         dy = -Ty;               // need to convert to raw image
03047         this_img->rotate(-ang_keep,0,0); // ming change this to this_img??
03048         this_img->translate(dx,dy,0.0); // ming change this to this_img
03049 
03050 
03051         Transform  tsoln(Dict("type","2d","alpha",ang_keep));
03052         tsoln.set_trans(dx,dy);
03053         this_img->set_attr("xform.align2d",&tsoln);
03054 #ifdef DEBUG
03055         float fsc_best=this_img->cmp(cmp_name,to,cmp_params);
03056         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);
03057 #endif
03058         return fsc;     // return the fsc coefficients
03059 } // FRM2D aligner sub_class

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

Definition at line 1617 of file aligner.cpp.

References refalifn(), v, and vp.

01617                                                                          {
01618         // we do this using a simple local difference estimate due to the expense of the calculation. 
01619         // The step has to be large enough for the similarity metric
01620         // To provide an accurate change in value. 
01621         static double lstep[4] = { 0.05, 0.05, 0.1, 0.01 }; 
01622         
01623         gsl_vector *vc = gsl_vector_alloc(v->size);
01624         gsl_vector_memcpy(vc,v);
01625         
01626         double f = refalifn(v,params);
01627         for (unsigned int i=0; i<v->size; i++) {
01628                 double *vp = gsl_vector_ptr(vc,i);
01629                 *vp+=lstep[i];
01630                 double f2 = refalifn(vc,params);
01631                 *vp-=lstep[i];
01632                 
01633                 gsl_vector_set(df,i,(f2-f)/lstep[i]);
01634         }
01635         
01636         gsl_vector_free(vc);
01637         return;
01638 }

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

Definition at line 1640 of file aligner.cpp.

References refalifn(), v, and vp.

01640                                                                                        {
01641         // we do this using a simple local difference estimate due to the expense of the calculation. 
01642         // The step has to be large enough for the similarity metric
01643         // To provide an accurate change in value. 
01644         static double lstep[4] = { 0.05, 0.05, 0.1, 0.01 }; 
01645         
01646         gsl_vector *vc = gsl_vector_alloc(v->size);
01647         gsl_vector_memcpy(vc,v);
01648         
01649         *f = refalifn(v,params);
01650         for (unsigned int i=0; i<v->size; i++) {
01651                 double *vp = gsl_vector_ptr(vc,i);
01652                 *vp+=lstep[i];
01653                 double f2 = refalifn(vc,params);
01654                 *vp-=lstep[i];
01655                 
01656                 gsl_vector_set(df,i,(f2-*f)/lstep[i]);
01657         }
01658         
01659         gsl_vector_free(vc);
01660         return;
01661 
01662 }

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

Definition at line 1585 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().

01586 {
01587         Dict *dict = (Dict *) params;
01588 
01589         double x = gsl_vector_get(v, 0);
01590         double y = gsl_vector_get(v, 1);
01591         double a = gsl_vector_get(v, 2);
01592 
01593         EMData *this_img = (*dict)["this"];
01594         EMData *with = (*dict)["with"];
01595         bool mirror = (*dict)["mirror"];
01596 
01597         Transform t(Dict("type","2d","alpha",static_cast<float>(a)));
01598         t.set_trans((float)x,(float)y);
01599         t.set_mirror(mirror);
01600         if (v->size>3) {
01601                 float sca=(float)gsl_vector_get(v, 3);
01602                 if (sca<.7 || sca>1.3) return 1.0e20;
01603                 t.set_scale((float)gsl_vector_get(v, 3));
01604         }
01605         EMData *tmp = this_img->process("xform",Dict("transform",&t));
01606         if (dict->has_key("mask")) tmp->mult(*(EMData *)((*dict)["mask"]));
01607 
01608 //      printf("GSL %f %f %f %d %f\n",x,y,a,mirror,(float)gsl_vector_get(v, 3));
01609         Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
01610         double result = c->cmp(tmp,with);
01611 
01612         if (tmp != 0) delete tmp;
01613         
01614         return result;
01615 }

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

Definition at line 1996 of file aligner.cpp.

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

01997 {
01998         Dict *dict = (Dict *) params;
01999 
02000         double n0 = gsl_vector_get(v, 0);
02001         double n1 = gsl_vector_get(v, 1);
02002         double n2 = gsl_vector_get(v, 2);
02003         double x = gsl_vector_get(v, 3);
02004         double y = gsl_vector_get(v, 4);
02005         double z = gsl_vector_get(v, 5);
02006 
02007         EMData *this_img = (*dict)["this"];
02008         EMData *with = (*dict)["with"];
02009 
02010         Transform* t = (*dict)["transform"];
02011         float spincoeff = (*dict)["spincoeff"];
02012 
02013         Transform soln = refalin3d_perturbquat(t,spincoeff,(float)n0,(float)n1,(float)n2,(float)x,(float)y,(float)z);
02014 
02015         EMData *tmp = this_img->process("xform",Dict("transform",&soln));
02016         Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
02017         double result = c->cmp(tmp,with);
02018         if ( tmp != 0 ) delete tmp;
02019 
02020         //cout << result << endl;
02021         return result;
02022 }

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

Definition at line 1664 of file aligner.cpp.

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

01665 {
01666         Dict *dict = (Dict *) params;
01667         EMData *this_img = (*dict)["this"];
01668         EMData *img_to = (*dict)["with"];
01669         bool mirror = (*dict)["mirror"];
01670 
01671         double x = gsl_vector_get(v, 0);
01672         double y = gsl_vector_get(v, 1);
01673         double a = gsl_vector_get(v, 2);
01674 
01675         double r = this_img->dot_rotate_translate(img_to, (float)x, (float)y, (float)a, mirror);
01676         int nsec = this_img->get_xsize() * this_img->get_ysize();
01677         double result = 1.0 - r / nsec;
01678 
01679 //      cout << result << " x " << x << " y " << y << " az " << a <<  endl;
01680         return result;
01681 }

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 1946 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().

01947 {
01948         Vec3f normal(n0,n1,n2);
01949         normal.normalize();
01950         
01951         float omega = spincoeff*sqrt(n0*n0 + n1*n1 + n2*n2); // Here we compute the spin by the rotation axis vector length
01952         Dict d;
01953         d["type"] = "spin";
01954         d["omega"] = omega;
01955         d["n1"] = normal[0];
01956         d["n2"] = normal[1];
01957         d["n3"] = normal[2];
01958         //cout << omega << " " << normal[0] << " " << normal[1] << " " << normal[2] << " " << n0 << " " << n1 << " " << n2 << endl;
01959         
01960         Transform q(d);
01961         q.set_trans((float)x,(float)y,(float)z);
01962         
01963         q = q*(*t); //compose transforms        
01964         
01965         return q;
01966 }

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

Definition at line 1968 of file aligner.cpp.

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

01969 {
01970         Dict *dict = (Dict *) params;
01971 
01972         double n0 = gsl_vector_get(v, 0);
01973         double n1 = gsl_vector_get(v, 1);
01974         double n2 = gsl_vector_get(v, 2);
01975         double x = gsl_vector_get(v, 3);
01976         double y = gsl_vector_get(v, 4);
01977         double z = gsl_vector_get(v, 5);
01978 
01979         EMData* volume = (*dict)["volume"];
01980         float spincoeff = (*dict)["spincoeff"];
01981         Transform* t = (*dict)["transform"];
01982 
01983         Transform soln = refalin3d_perturbquat(t,spincoeff,(float)n0,(float)n1,(float)n2,(float)x,(float)y,(float)z);
01984 
01985         EMData *tmp = volume->process("xform",Dict("transform",&soln));
01986         EMData *symtmp = tmp->process("xform.applysym",Dict("sym",(*dict)["sym"]));
01987         Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]);
01988         double result = c->cmp(symtmp,tmp);
01989         delete tmp;
01990         delete symtmp;
01991 
01992         //cout << result << endl;
01993         return result;
01994 }


Generated on Tue Jun 11 13:40:46 2013 for EMAN2 by  doxygen 1.3.9.1