#include <cmp.h>
Inheritance diagram for EMAN::PhaseCmp:
Public Member Functions | |
float | cmp (EMData *image, EMData *with) const |
To compare 'image' with another image passed in through its parameters. | |
string | get_name () const |
Get the Cmp's name. | |
string | get_desc () const |
TypeDict | get_param_types () const |
Get Cmp parameter information in a dictionary. | |
Static Public Member Functions | |
static Cmp * | NEW () |
Static Public Attributes | |
static const string | NAME = "phase" |
SNR should be an array as returned by ctfcurve() 'with' should be the less noisy image, since it's amplitudes will be used to weight the phase residual. 2D only.
Use Phase Residual as a measure of similarity Differential phase residual (DPR) is a measure of statistical dependency between two averages, computed over rings in Fourier space as a function of ring radius (= spatial frequency, or resolution)
Definition at line 463 of file cmp.h.
To compare 'image' with another image passed in through its parameters.
An optional transformation may be used to transform the 2 images.
image | The first image to be compared. | |
with | The second image to be comppared. |
Implements EMAN::Cmp.
Definition at line 929 of file cmp.cpp.
References EMAN::Util::angle_err_ri(), EMAN::Ctf::apix, EMAN::Ctf::compute_1d(), EMAN::EMData::copy(), EMAN::Ctf::CTF_SNR, EMAN::EMData::do_fft(), EMAN::EMData::do_fft_inplace(), EMAN::EMData::do_ift_inplace(), ENTERFUNC, EXITFUNC, EMAN::EMObject::f, EMAN::EMData::get_attr(), EMAN::EMData::get_const_data(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::EMData::has_attr(), EMAN::Util::hypot3(), EMAN::Util::hypot_fast_int(), InvalidCallException, EMAN::EMData::is_complex(), norm(), ny, EMAN::Cmp::params, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::update(), x, and y.
00930 { 00931 ENTERFUNC; 00932 00933 int snrweight = params.set_default("snrweight", 0); 00934 int snrfn = params.set_default("snrfn",0); 00935 int ampweight = params.set_default("ampweight",0); 00936 int zeromask = params.set_default("zeromask",0); 00937 float minres = params.set_default("minres",500.0f); 00938 float maxres = params.set_default("maxres",2.0f); 00939 00940 if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator"); 00941 00942 //#ifdef EMAN2_USING_CUDA 00943 // if (image->gpu_operation_preferred()) { 00944 // cout << "Cuda cmp" << endl; 00945 // EXITFUNC; 00946 // return cuda_cmp(image,with); 00947 // } 00948 //#endif 00949 00950 EMData *image_fft = NULL; 00951 EMData *with_fft = NULL; 00952 00953 int ny = image->get_ysize(); 00954 // int np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2; 00955 int np = 0; 00956 vector<float> snr; 00957 00958 // weighting based on SNR estimate from CTF 00959 if (snrweight) { 00960 Ctf *ctf = NULL; 00961 if (!image->has_attr("ctf")) { 00962 if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters"); 00963 ctf=with->get_attr("ctf"); 00964 } 00965 else ctf=image->get_attr("ctf"); 00966 00967 float ds=1.0f/(ctf->apix*ny); 00968 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR); // note that this returns ny/2 values 00969 if(ctf) {delete ctf; ctf=0;} 00970 np=snr.size(); 00971 } 00972 // weighting based on empirical SNR function (not really good) 00973 else if (snrfn==1) { 00974 np=ny/2; 00975 snr.resize(np); 00976 00977 for (int i = 0; i < np; i++) { 00978 float x2 = 10.0f*i/np; 00979 snr[i] = x2 * exp(-x2); 00980 } 00981 } 00982 // no weighting 00983 else { 00984 np=ny/2; 00985 snr.resize(np); 00986 00987 for (int i = 0; i < np; i++) snr[i]=1.0; 00988 } 00989 00990 // Min/max modifications to weighting 00991 float pmin,pmax; 00992 if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres; //cutoff in pixels, assume square 00993 else pmin=0; 00994 if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres; 00995 else pmax=0; 00996 00997 // printf("%f\t%f\n",pmin,pmax); 00998 00999 // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation 01000 for (int i = 0; i < np; i++) { 01001 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f; 01002 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f; 01003 // printf("%d\t%f\n",i,snr[i]); 01004 } 01005 01006 if (zeromask) { 01007 image_fft=image->copy(); 01008 with_fft=with->copy(); 01009 01010 if (image_fft->is_complex()) image_fft->do_ift_inplace(); 01011 if (with_fft->is_complex()) with_fft->do_ift_inplace(); 01012 01013 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize(); 01014 float *d1=image_fft->get_data(); 01015 float *d2=with_fft->get_data(); 01016 01017 for (int i=0; i<sz; i++) { 01018 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; } 01019 } 01020 01021 image_fft->update(); 01022 with_fft->update(); 01023 image_fft->do_fft_inplace(); 01024 with_fft->do_fft_inplace(); 01025 image_fft->set_attr("free_me",1); 01026 with_fft->set_attr("free_me",1); 01027 } 01028 else { 01029 if (image->is_complex()) image_fft=image; 01030 else { 01031 image_fft=image->do_fft(); 01032 image_fft->set_attr("free_me",1); 01033 } 01034 01035 if (with->is_complex()) with_fft=with; 01036 else { 01037 with_fft=with->do_fft(); 01038 with_fft->set_attr("free_me",1); 01039 } 01040 } 01041 // image_fft->ri2ap(); 01042 // with_fft->ri2ap(); 01043 01044 const float *const image_fft_data = image_fft->get_const_data(); 01045 const float *const with_fft_data = with_fft->get_const_data(); 01046 double sum = 0; 01047 double norm = FLT_MIN; 01048 size_t i = 0; 01049 // int nx=image_fft->get_xsize(); 01050 ny=image_fft->get_ysize(); 01051 int nz=image_fft->get_zsize(); 01052 int nx2=image_fft->get_xsize()/2; 01053 int ny2=image_fft->get_ysize()/2; 01054 // int nz2=image_fft->get_zsize()/2; 01055 01056 // This can never happen any more... 01057 if (np==0) { 01058 for (int z = 0; z < nz; z++){ 01059 for (int y = 0; y < ny; y++) { 01060 for (int x = 0; x < nx2; x ++) { 01061 float a; 01062 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01063 else a=1.0f; 01064 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01065 norm += a; 01066 i += 2; 01067 } 01068 } 01069 } 01070 01071 } 01072 else if (nz==1) { 01073 for (int y = 0; y < ny; y++) { 01074 for (int x = 0; x < nx2; x ++) { 01075 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y); 01076 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 01077 01078 float a; 01079 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01080 else a=1.0f; 01081 a*=snr[r]; 01082 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01083 norm += a; 01084 i += 2; 01085 } 01086 } 01087 } 01088 else { 01089 for (int z = 0; z < nz; z++){ 01090 for (int y = 0; y < ny; y++) { 01091 for (int x = 0; x < nx2; x ++) { 01092 int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z); 01093 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 01094 01095 float a; 01096 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01097 else a=1.0f; 01098 a*=snr[r]; 01099 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01100 norm += a; 01101 i += 2; 01102 } 01103 } 01104 } 01105 01106 } 01107 01108 EXITFUNC; 01109 01110 if( image_fft->has_attr("free_me") ) 01111 { 01112 delete image_fft; 01113 image_fft = 0; 01114 } 01115 if( with_fft->has_attr("free_me") ) 01116 { 01117 delete with_fft; 01118 with_fft = 0; 01119 } 01120 #if 0 01121 return (1.0f - sum / norm); 01122 #endif 01123 return (float)(sum / norm); 01124 }
string EMAN::PhaseCmp::get_desc | ( | ) | const [inline, virtual] |
string EMAN::PhaseCmp::get_name | ( | ) | const [inline, virtual] |
TypeDict EMAN::PhaseCmp::get_param_types | ( | ) | const [inline, virtual] |
Get Cmp parameter information in a dictionary.
Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.
Implements EMAN::Cmp.
Definition at line 483 of file cmp.h.
References EMAN::EMObject::FLOAT, EMAN::EMObject::INT, and EMAN::TypeDict::put().
00484 { 00485 TypeDict d; 00486 d.put("snrweight", EMObject::INT, "If set, the SNR of 'this' will be used to weight the result. If 'this' lacks CTF info, it will check 'with'. (default=0)"); 00487 d.put("snrfn", EMObject::INT, "If nonzero, an empirical function will be used as a radial weight rather than the true SNR. (1 - exp decay)'. (default=0)"); 00488 d.put("ampweight", EMObject::INT, "If set, the amplitude of 'with' will be used as a weight in the averaging'. (default=0)"); 00489 d.put("zeromask", EMObject::INT, "Treat regions in either image that are zero as a mask"); 00490 d.put("minres", EMObject::FLOAT, "Lowest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=500"); 00491 d.put("maxres", EMObject::FLOAT, "Highest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=10"); 00492 return d; 00493 }
static Cmp* EMAN::PhaseCmp::NEW | ( | ) | [inline, static] |
const string PhaseCmp::NAME = "phase" [static] |