#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:
Go to the source code of this file.
Defines | |
#define | EMAN2_ALIGNER_DEBUG 0 |
| |
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) |
|
Definition at line 55 of file aligner.cpp. |
|
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
|
|
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 }
|
|
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 }
|
|
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 }
|
|
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 }
|
|
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 }
|
|
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 }
|
|
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 }
|