#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 | |
static double | refalifn (const gsl_vector *v, void *params) |
static double | refalifnfast (const gsl_vector *v, void *params) |
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) |
static double | symquat (const gsl_vector *v, void *params) |
static 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) |
void | EMAN::dump_aligners () |
map< string, vector< string > > | EMAN::dump_aligners_list () |
#define EMAN2_ALIGNER_DEBUG 0 |
float @0::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 2712 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::Transform::set_trans(), and EMAN::EMData::translate().
Referenced by EMAN::FRM2DAligner::align().
02713 { 02714 int size=rsize; 02715 float dx,dy; 02716 int bw=size/2; 02717 int MAXR=this_img->get_ysize()/2; 02718 //int MAXR=size; 02719 unsigned long tsize=2*size; 02720 unsigned long ind1=0, ind2=0, ind3=0, ind4=0, ind41=0; 02721 unsigned long index0=0; 02722 int i=0, j=0, m=0, n=0, r=0; 02723 int loop_rho=0, rho_best=0; 02724 02725 float* gnr2 = new float[size*2]; 02726 float* maxcor = new float[size+1]; // MAXR need change 02727 02728 int p_max=p_max_input; 02729 float* result = new float[5*(p_max+1)]; 02730 float* cr=new float[size*(bw+1)]; 02731 float* ci=new float[size*(bw+1)]; 02732 EMData *data_in=new EMData; 02733 data_in->set_complex(true); 02734 data_in->set_fftodd(false); 02735 data_in->set_ri(true); 02736 data_in->set_size(size+2,size,1); 02737 float *in=data_in->get_data(); 02738 02739 float *self_sampl_fft = selfpcsfft->get_data(); // ming f(r) 02740 02741 float maxcor_sofar=0.0f; 02742 int p=0; 02743 02744 for(p=0; p<=p_max; ++p){ 02745 ind1=p*size*bw; 02746 for (i=0;i<size;++i) 02747 for (j=0;j<bw+1;++j){ 02748 cr[i*(bw+1)+j]=0.0; 02749 ci[i*(bw+1)+j]=0.0; 02750 } 02751 for(n=0;n<bw;++n){ // loop for n 02752 ind2=(ind1+n); 02753 index0=n*(bw+1); 02754 for(r=0;r<=MAXR;++r) { 02755 ind3=(ind2+r*bw)*size; 02756 for(m=0;m<size;m++){ // take back hat{h(n,r,p)}(m) 02757 ind4=(ind3+m)*2; 02758 ind41=ind4+1; 02759 gnr2[2*m]=frm2dhhat[ind4]; 02760 gnr2[2*m+1]=frm2dhhat[ind41]; 02761 } 02762 for(m=0;m<bw;++m){ 02763 float tempr=self_sampl_fft[2*m+r*(size+2)]*r; 02764 float tempi=self_sampl_fft[2*m+1+r*(size+2)]*r; 02765 float gnr2_r=gnr2[2*m]; 02766 float gnr2_i=gnr2[2*m+1]; 02767 cr[n*(bw+1)+m]+=gnr2_r*tempr+gnr2_i*tempi; 02768 ci[n*(bw+1)+m]+=gnr2_i*tempr-gnr2_r*tempi; 02769 if(n!=0){ // m,-n 02770 if(m!= 0){ 02771 int ssize=tsize-2*m; // ssize = 2*size-2m 02772 int ssize1=ssize+1; 02773 float gnr2_r=gnr2[ssize]; 02774 float gnr2_i=gnr2[ssize1]; 02775 cr[(size-n)*(bw+1)+m]+=gnr2_r*tempr-gnr2_i*tempi; 02776 ci[(size-n)*(bw+1)+m]-=gnr2_i*tempr+gnr2_r*tempi; 02777 } 02778 else{ 02779 cr[(size-n)*(bw+1)+m]+=*(gnr2)*tempr-*(gnr2+1)*tempi; 02780 ci[(size-n)*(bw+1)+m]-=*(gnr2+1)*tempr+*(gnr2)*tempi; 02781 } 02782 } 02783 } 02784 } 02785 } 02786 for (int cii=0; cii<size*(bw+1);++cii){ 02787 in[2*cii]=cr[cii]; 02788 in[2*cii+1]=ci[cii]; 02789 //printf("cii=%d,in[2i+1]=%f\n",cii, cr[cii]); 02790 } 02791 02792 EMData *data_out; 02793 data_out=data_in->do_ift(); 02794 float *c=data_out->get_data(); 02795 float tempr=0.0f, corre_fcs=999.0f; 02796 02797 int n_best=0, m_best=0; 02798 float temp=-100.0f; 02799 for(n=0;n<size;++n){// move Tri_2D to Tri = c(phi,phi';rho) 02800 for(m=0;m<size;++m){ 02801 temp=c[n*size+m]; 02802 if(temp>tempr) { 02803 tempr=temp; 02804 n_best=n; 02805 m_best=m; 02806 } 02807 } 02808 } 02809 delete data_out; 02810 02811 float corre,Phi2,Phi,Tx,Ty,Vx, Vy; 02812 02813 //for (n_best=0;n_best<bw;n_best++) 02814 // for (m_best=0;m_best<2*bw;m_best++){ 02815 //n_best=0; 02816 //m_best=70; 02817 Phi2=n_best*M_PI/bw; // ming this is reference image rotation angle 02818 Phi=m_best*M_PI/bw; // ming this is particle image rotation angle 02819 Vx=p*cos(Phi);//should use the angle of the centered one 02820 Vy=-p*sin(Phi); 02821 Tx=Vx+(floor(com_this_x+0.5f)-floor(com_with_x+0.5f)); 02822 Ty=Vy+(floor(com_this_y+0.5f)-floor(com_with_y+0.5f)); 02823 02824 dx=-Tx; // the Rota & Trans value (Tx,Ty, ang_keep) are for the projection image, 02825 dy=-Ty; // need to convert to raw image 02826 02827 EMData *this_tmp=this_img->copy();//ming change to to 02828 this_tmp->rotate(-(Phi2-Phi)*180.0f,0.0f,0.0f); 02829 this_tmp->translate(dx,dy,0.0); 02830 02831 corre=this_tmp->cmp(cmp_name,to,cmp_params); 02832 //printf("corre=%f\n",corre); 02833 delete this_tmp; 02834 if(corre<=corre_fcs) { //ming, cmp use smaller value stands for more similarity 02835 corre_fcs=corre; 02836 result[0+5*p] = float(p); // rho 02837 result[1+5*p] = corre; // correlation_fcs 02838 result[2+5*p] = (Phi2-Phi)*180.0f; // rotation angle 02839 result[3+5*p] = Tx; // Translation_x 02840 result[4+5*p] = Ty; // Translation_y 02841 } 02842 maxcor[p]=corre_fcs; // maximum correlation for current rho 02843 if(corre_fcs<maxcor_sofar) { 02844 maxcor_sofar=corre_fcs; // max correlation up to current rho 02845 rho_best=p; // the rho value with maxinum correlation value 02846 } 02847 if(p>=4){ 02848 if(maxcor[p] < maxcor[p-1] && maxcor[p-1] < maxcor[p-2]&& maxcor[p-2] < maxcor[p-3] && maxcor[p-3] < maxcor[p-4]){ 02849 loop_rho=1; 02850 break; //exit p loop 02851 } 02852 } 02853 } // end for p 02854 //}//test my method 02855 if(loop_rho == 1) 02856 p=p+1; 02857 int rb5=5*rho_best; 02858 float fsc = result[1+rb5]; 02859 float ang_keep = result[2+rb5]; 02860 float Tx = result[3+rb5]; 02861 float Ty = result[4+rb5]; 02862 delete[] gnr2; 02863 delete[] maxcor; 02864 delete[] result; 02865 delete cr; 02866 cr=0; 02867 delete ci; 02868 ci=0; 02869 delete data_in; // ming add 02870 dx = -Tx; // the Rota & Trans value (Tx,Ty, ang_keep) are for the projection image, 02871 dy = -Ty; // need to convert to raw image 02872 this_img->rotate(-ang_keep,0,0); // ming change this to this_img?? 02873 this_img->translate(dx,dy,0.0); // ming change this to this_img 02874 02875 02876 Transform tsoln(Dict("type","2d","alpha",ang_keep)); 02877 tsoln.set_trans(dx,dy); 02878 this_img->set_attr("xform.align2d",&tsoln); 02879 #ifdef DEBUG 02880 float fsc_best=this_img->cmp(cmp_name,to,cmp_params); 02881 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); 02882 #endif 02883 return fsc; // return the fsc coefficients 02884 } // FRM2D aligner sub_class
static double refalifn | ( | const gsl_vector * | v, | |
void * | params | |||
) | [static] |
Definition at line 1564 of file aligner.cpp.
References EMAN::Cmp::cmp(), EMAN::EMData::process(), t, x, and y.
Referenced by EMAN::RefineAligner::align().
01565 { 01566 Dict *dict = (Dict *) params; 01567 01568 double x = gsl_vector_get(v, 0); 01569 double y = gsl_vector_get(v, 1); 01570 double a = gsl_vector_get(v, 2); 01571 01572 EMData *this_img = (*dict)["this"]; 01573 EMData *with = (*dict)["with"]; 01574 bool mirror = (*dict)["mirror"]; 01575 01576 // float mean = (float)this_img->get_attr("mean"); 01577 // if ( Util::goodf(&mean) ) { 01578 // //cout << "tmps mean is nan even before rotation" << endl; 01579 // } 01580 01581 Transform t(Dict("type","2d","alpha",static_cast<float>(a))); 01582 // Transform3D t3d(Transform3D::EMAN, (float)a, 0.0f, 0.0f); 01583 // t3d.set_posttrans( (float) x, (float) y); 01584 // tmp->rotate_translate(t3d); 01585 t.set_trans((float)x,(float)y); 01586 t.set_mirror(mirror); 01587 if (v->size>3) { 01588 float sca=(float)gsl_vector_get(v, 3); 01589 if (sca<.7 || sca>1.3) return 1.0e20; 01590 t.set_scale((float)gsl_vector_get(v, 3)); 01591 } 01592 EMData *tmp = this_img->process("xform",Dict("transform",&t)); 01593 01594 // printf("GSL %f %f %f %d %f\n",x,y,a,mirror,(float)gsl_vector_get(v, 3)); 01595 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]); 01596 double result = c->cmp(tmp,with); 01597 01598 // DELETE AT SOME STAGE, USEFUL FOR PRERELEASE STUFF 01599 // float test_result = (float)result; 01600 // if ( Util::goodf(&test_result) ) { 01601 // cout << "result " << result << " " << x << " " << y << " " << a << endl; 01602 // cout << (float)this_img->get_attr("mean") << " " << (float)tmp->get_attr("mean") << " " << (float)with->get_attr("mean") << endl; 01603 // tmp->write_image("tmp.hdf"); 01604 // with->write_image("with.hdf"); 01605 // this_img->write_image("this_img.hdf"); 01606 // EMData* t = this_img->copy(); 01607 // cout << (float)t->get_attr("mean") << endl; 01608 // t->rotate_translate( t3d ); 01609 // cout << (float)t->get_attr("mean") << endl; 01610 // cout << "exit" << endl; 01612 // cout << (float)t->get_attr("mean") << endl; 01613 // cout << "now exit" << endl; 01614 // delete t; 01615 // } 01616 01617 01618 if (tmp != 0) delete tmp; 01619 01620 return result; 01621 }
static double refalifn3dquat | ( | const gsl_vector * | v, | |
void * | params | |||
) | [static] |
Definition at line 1826 of file aligner.cpp.
References EMAN::Cmp::cmp(), EMAN::EMData::process(), refalin3d_perturbquat(), t, x, and y.
Referenced by EMAN::Refine3DAlignerQuaternion::align().
01827 { 01828 Dict *dict = (Dict *) params; 01829 01830 double n0 = gsl_vector_get(v, 0); 01831 double n1 = gsl_vector_get(v, 1); 01832 double n2 = gsl_vector_get(v, 2); 01833 double x = gsl_vector_get(v, 3); 01834 double y = gsl_vector_get(v, 4); 01835 double z = gsl_vector_get(v, 5); 01836 01837 EMData *this_img = (*dict)["this"]; 01838 EMData *with = (*dict)["with"]; 01839 01840 Transform* t = (*dict)["transform"]; 01841 float spincoeff = (*dict)["spincoeff"]; 01842 01843 Transform soln = refalin3d_perturbquat(t,spincoeff,(float)n0,(float)n1,(float)n2,(float)x,(float)y,(float)z); 01844 01845 EMData *tmp = this_img->process("xform",Dict("transform",&soln)); 01846 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]); 01847 double result = c->cmp(tmp,with); 01848 if ( tmp != 0 ) delete tmp; 01849 01850 //cout << result << endl; 01851 return result; 01852 }
static double refalifnfast | ( | const gsl_vector * | v, | |
void * | params | |||
) | [static] |
Definition at line 1623 of file aligner.cpp.
References EMAN::EMData::dot_rotate_translate(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), x, and y.
Referenced by EMAN::RefineAligner::align().
01624 { 01625 Dict *dict = (Dict *) params; 01626 EMData *this_img = (*dict)["this"]; 01627 EMData *img_to = (*dict)["with"]; 01628 bool mirror = (*dict)["mirror"]; 01629 01630 double x = gsl_vector_get(v, 0); 01631 double y = gsl_vector_get(v, 1); 01632 double a = gsl_vector_get(v, 2); 01633 01634 double r = this_img->dot_rotate_translate(img_to, (float)x, (float)y, (float)a, mirror); 01635 int nsec = this_img->get_xsize() * this_img->get_ysize(); 01636 double result = 1.0 - r / nsec; 01637 01638 // cout << result << " x " << x << " y " << y << " az " << a << endl; 01639 return result; 01640 }
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 | |||
) | [static] |
Definition at line 1776 of file aligner.cpp.
References EMAN::Vec3< Type >::normalize(), q, and sqrt().
Referenced by EMAN::Refine3DAlignerQuaternion::align(), EMAN::SymAlignProcessorQuat::align(), refalifn3dquat(), and symquat().
01777 { 01778 Vec3f normal(n0,n1,n2); 01779 normal.normalize(); 01780 01781 float omega = spincoeff*sqrt(n0*n0 + n1*n1 + n2*n2); // Here we compute the spin by the rotation axis vector length 01782 Dict d; 01783 d["type"] = "spin"; 01784 d["Omega"] = omega; 01785 d["n1"] = normal[0]; 01786 d["n2"] = normal[1]; 01787 d["n3"] = normal[2]; 01788 //cout << omega << " " << normal[0] << " " << normal[1] << " " << normal[2] << " " << n0 << " " << n1 << " " << n2 << endl; 01789 01790 Transform q(d); 01791 q.set_trans((float)x,(float)y,(float)z); 01792 01793 q = q*(*t); //compose transforms 01794 01795 return q; 01796 }
static double symquat | ( | const gsl_vector * | v, | |
void * | params | |||
) | [static] |
Definition at line 1798 of file aligner.cpp.
References EMAN::Cmp::cmp(), EMAN::EMData::process(), refalin3d_perturbquat(), t, x, and y.
Referenced by EMAN::SymAlignProcessorQuat::align().
01799 { 01800 Dict *dict = (Dict *) params; 01801 01802 double n0 = gsl_vector_get(v, 0); 01803 double n1 = gsl_vector_get(v, 1); 01804 double n2 = gsl_vector_get(v, 2); 01805 double x = gsl_vector_get(v, 3); 01806 double y = gsl_vector_get(v, 4); 01807 double z = gsl_vector_get(v, 5); 01808 01809 EMData* volume = (*dict)["volume"]; 01810 float spincoeff = (*dict)["spincoeff"]; 01811 Transform* t = (*dict)["transform"]; 01812 01813 Transform soln = refalin3d_perturbquat(t,spincoeff,(float)n0,(float)n1,(float)n2,(float)x,(float)y,(float)z); 01814 01815 EMData *tmp = volume->process("xform",Dict("transform",&soln)); 01816 EMData *symtmp = tmp->process("xform.applysym",Dict("sym",(*dict)["sym"])); 01817 Cmp* c = (Cmp*) ((void*)(*dict)["cmp"]); 01818 double result = c->cmp(symtmp,tmp); 01819 delete tmp; 01820 delete symtmp; 01821 01822 //cout << result << endl; 01823 return result; 01824 }