#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 | |
Cmp * | NEW () |
Static Public Attributes | |
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.
Implements EMAN::Cmp. Definition at line 930 of file cmp.cpp. References EMAN::Util::angle_err_ri(), EMAN::Ctf::apix, EMAN::Ctf::compute_1d(), EMAN::EMData::copy(), EMAN::EMData::do_fft(), EMAN::EMData::do_fft_inplace(), EMAN::EMData::do_ift_inplace(), 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::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::update(), x, and y. 00931 { 00932 ENTERFUNC; 00933 00934 int snrweight = params.set_default("snrweight", 0); 00935 int snrfn = params.set_default("snrfn",0); 00936 int ampweight = params.set_default("ampweight",0); 00937 int zeromask = params.set_default("zeromask",0); 00938 float minres = params.set_default("minres",500.0f); 00939 float maxres = params.set_default("maxres",2.0f); 00940 00941 if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator"); 00942 00943 EMData *image_fft = NULL; 00944 EMData *with_fft = NULL; 00945 00946 int ny = image->get_ysize(); 00947 // int np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2; 00948 int np = 0; 00949 vector<float> snr; 00950 00951 // weighting based on SNR estimate from CTF 00952 if (snrweight) { 00953 Ctf *ctf = NULL; 00954 if (!image->has_attr("ctf")) { 00955 if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters"); 00956 ctf=with->get_attr("ctf"); 00957 } 00958 else ctf=image->get_attr("ctf"); 00959 00960 float ds=1.0f/(ctf->apix*ny); 00961 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR); // note that this returns ny/2 values 00962 if(ctf) {delete ctf; ctf=0;} 00963 np=snr.size(); 00964 } 00965 // weighting based on empirical SNR function (not really good) 00966 else if (snrfn==1) { 00967 np=ny/2; 00968 snr.resize(np); 00969 00970 for (int i = 0; i < np; i++) { 00971 float x2 = 10.0f*i/np; 00972 snr[i] = x2 * exp(-x2); 00973 } 00974 } 00975 // no weighting 00976 else { 00977 np=ny/2; 00978 snr.resize(np); 00979 00980 for (int i = 0; i < np; i++) snr[i]=1.0; 00981 } 00982 00983 // Min/max modifications to weighting 00984 float pmin,pmax; 00985 if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres; //cutoff in pixels, assume square 00986 else pmin=0; 00987 if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres; 00988 else pmax=0; 00989 00990 // printf("%f\t%f\n",pmin,pmax); 00991 00992 // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation 00993 for (int i = 0; i < np; i++) { 00994 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f; 00995 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f; 00996 // printf("%d\t%f\n",i,snr[i]); 00997 } 00998 00999 if (zeromask) { 01000 image_fft=image->copy(); 01001 with_fft=with->copy(); 01002 01003 if (image_fft->is_complex()) image_fft->do_ift_inplace(); 01004 if (with_fft->is_complex()) with_fft->do_ift_inplace(); 01005 01006 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize(); 01007 float *d1=image_fft->get_data(); 01008 float *d2=with_fft->get_data(); 01009 01010 for (int i=0; i<sz; i++) { 01011 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; } 01012 } 01013 01014 image_fft->update(); 01015 with_fft->update(); 01016 image_fft->do_fft_inplace(); 01017 with_fft->do_fft_inplace(); 01018 image_fft->set_attr("free_me",1); 01019 with_fft->set_attr("free_me",1); 01020 } 01021 else { 01022 if (image->is_complex()) image_fft=image; 01023 else { 01024 image_fft=image->do_fft(); 01025 image_fft->set_attr("free_me",1); 01026 } 01027 01028 if (with->is_complex()) with_fft=with; 01029 else { 01030 with_fft=with->do_fft(); 01031 with_fft->set_attr("free_me",1); 01032 } 01033 } 01034 // image_fft->ri2ap(); 01035 // with_fft->ri2ap(); 01036 01037 const float *const image_fft_data = image_fft->get_const_data(); 01038 const float *const with_fft_data = with_fft->get_const_data(); 01039 double sum = 0; 01040 double norm = FLT_MIN; 01041 size_t i = 0; 01042 // int nx=image_fft->get_xsize(); 01043 ny=image_fft->get_ysize(); 01044 int nz=image_fft->get_zsize(); 01045 int nx2=image_fft->get_xsize()/2; 01046 int ny2=image_fft->get_ysize()/2; 01047 // int nz2=image_fft->get_zsize()/2; 01048 01049 // This can never happen any more... 01050 if (np==0) { 01051 for (int z = 0; z < nz; z++){ 01052 for (int y = 0; y < ny; y++) { 01053 for (int x = 0; x < nx2; x ++) { 01054 float a; 01055 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01056 else a=1.0f; 01057 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01058 norm += a; 01059 i += 2; 01060 } 01061 } 01062 } 01063 01064 } 01065 else if (nz==1) { 01066 for (int y = 0; y < ny; y++) { 01067 for (int x = 0; x < nx2; x ++) { 01068 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y); 01069 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 01070 01071 float a; 01072 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01073 else a=1.0f; 01074 a*=snr[r]; 01075 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01076 norm += a; 01077 i += 2; 01078 } 01079 } 01080 } 01081 else { 01082 for (int z = 0; z < nz; z++){ 01083 for (int y = 0; y < ny; y++) { 01084 for (int x = 0; x < nx2; x ++) { 01085 int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z); 01086 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 01087 01088 float a; 01089 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01090 else a=1.0f; 01091 a*=snr[r]; 01092 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01093 norm += a; 01094 i += 2; 01095 } 01096 } 01097 } 01098 01099 } 01100 01101 EXITFUNC; 01102 01103 if( image_fft->has_attr("free_me") ) 01104 { 01105 delete image_fft; 01106 image_fft = 0; 01107 } 01108 if( with_fft->has_attr("free_me") ) 01109 { 01110 delete with_fft; 01111 with_fft = 0; 01112 } 01113 #if 0 01114 return (1.0f - sum / norm); 01115 #endif 01116 return (float)(sum / norm); 01117 }
|
|
Implements EMAN::Cmp. Definition at line 473 of file cmp.h. 00474 { 00475 return "Mean phase difference"; 00476 }
|
|
Get the Cmp's name. Each Cmp is identified by a unique name.
Implements EMAN::Cmp. Definition at line 468 of file cmp.h. 00469 {
00470 return NAME;
00471 }
|
|
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::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 }
|
|
Definition at line 478 of file cmp.h. 00479 { 00480 return new PhaseCmp(); 00481 }
|
|
|