#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 507 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 1056 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.
01057 { 01058 ENTERFUNC; 01059 01060 int snrweight = params.set_default("snrweight", 0); 01061 int snrfn = params.set_default("snrfn",0); 01062 int ampweight = params.set_default("ampweight",0); 01063 int zeromask = params.set_default("zeromask",0); 01064 float minres = params.set_default("minres",500.0f); 01065 float maxres = params.set_default("maxres",2.0f); 01066 01067 if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator"); 01068 01069 EMData *image_fft = NULL; 01070 EMData *with_fft = NULL; 01071 01072 int ny = image->get_ysize(); 01073 // int np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2; 01074 int np = 0; 01075 vector<float> snr; 01076 01077 // weighting based on SNR estimate from CTF 01078 if (snrweight) { 01079 Ctf *ctf = NULL; 01080 if (!image->has_attr("ctf")) { 01081 if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters"); 01082 ctf=with->get_attr("ctf"); 01083 } 01084 else ctf=image->get_attr("ctf"); 01085 01086 float ds=1.0f/(ctf->apix*ny); 01087 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR); // note that this returns ny/2 values 01088 if(ctf) {delete ctf; ctf=0;} 01089 np=snr.size(); 01090 } 01091 // weighting based on empirical SNR function (not really good) 01092 else if (snrfn==1) { 01093 np=ny/2; 01094 snr.resize(np); 01095 01096 for (int i = 0; i < np; i++) { 01097 float x2 = 10.0f*i/np; 01098 snr[i] = x2 * exp(-x2); 01099 } 01100 } 01101 // no weighting 01102 else { 01103 np=ny/2; 01104 snr.resize(np); 01105 01106 for (int i = 0; i < np; i++) snr[i]=1.0; 01107 } 01108 01109 // Min/max modifications to weighting 01110 float pmin,pmax; 01111 if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres; //cutoff in pixels, assume square 01112 else pmin=0; 01113 if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres; 01114 else pmax=0; 01115 01116 // printf("%f\t%f\n",pmin,pmax); 01117 01118 // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation 01119 for (int i = 0; i < np; i++) { 01120 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f; 01121 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f; 01122 // printf("%d\t%f\n",i,snr[i]); 01123 } 01124 01125 if (zeromask) { 01126 image_fft=image->copy(); 01127 with_fft=with->copy(); 01128 01129 if (image_fft->is_complex()) image_fft->do_ift_inplace(); 01130 if (with_fft->is_complex()) with_fft->do_ift_inplace(); 01131 01132 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize(); 01133 float *d1=image_fft->get_data(); 01134 float *d2=with_fft->get_data(); 01135 01136 for (int i=0; i<sz; i++) { 01137 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; } 01138 } 01139 01140 image_fft->update(); 01141 with_fft->update(); 01142 image_fft->do_fft_inplace(); 01143 with_fft->do_fft_inplace(); 01144 image_fft->set_attr("free_me",1); 01145 with_fft->set_attr("free_me",1); 01146 } 01147 else { 01148 if (image->is_complex()) image_fft=image; 01149 else { 01150 image_fft=image->do_fft(); 01151 image_fft->set_attr("free_me",1); 01152 } 01153 01154 if (with->is_complex()) with_fft=with; 01155 else { 01156 with_fft=with->do_fft(); 01157 with_fft->set_attr("free_me",1); 01158 } 01159 } 01160 // image_fft->ri2ap(); 01161 // with_fft->ri2ap(); 01162 01163 const float *const image_fft_data = image_fft->get_const_data(); 01164 const float *const with_fft_data = with_fft->get_const_data(); 01165 double sum = 0; 01166 double norm = FLT_MIN; 01167 size_t i = 0; 01168 // int nx=image_fft->get_xsize(); 01169 ny=image_fft->get_ysize(); 01170 int nz=image_fft->get_zsize(); 01171 int nx2=image_fft->get_xsize()/2; 01172 int ny2=image_fft->get_ysize()/2; 01173 // int nz2=image_fft->get_zsize()/2; 01174 01175 // This can never happen any more... 01176 if (np==0) { 01177 for (int z = 0; z < nz; z++){ 01178 for (int y = 0; y < ny; y++) { 01179 for (int x = 0; x < nx2; x ++) { 01180 float a; 01181 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01182 else a=1.0f; 01183 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01184 norm += a; 01185 i += 2; 01186 } 01187 } 01188 } 01189 01190 } 01191 else if (nz==1) { 01192 for (int y = 0; y < ny; y++) { 01193 for (int x = 0; x < nx2; x ++) { 01194 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y); 01195 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 01196 01197 float a; 01198 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01199 else a=1.0f; 01200 a*=snr[r]; 01201 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01202 norm += a; 01203 i += 2; 01204 } 01205 } 01206 } 01207 else { 01208 for (int z = 0; z < nz; z++){ 01209 for (int y = 0; y < ny; y++) { 01210 for (int x = 0; x < nx2; x ++) { 01211 int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z); 01212 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 01213 01214 float a; 01215 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01216 else a=1.0f; 01217 a*=snr[r]; 01218 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01219 norm += a; 01220 i += 2; 01221 } 01222 } 01223 } 01224 01225 } 01226 01227 EXITFUNC; 01228 01229 if( image_fft->has_attr("free_me") ) 01230 { 01231 delete image_fft; 01232 image_fft = 0; 01233 } 01234 if( with_fft->has_attr("free_me") ) 01235 { 01236 delete with_fft; 01237 with_fft = 0; 01238 } 01239 #if 0 01240 return (1.0f - sum / norm); 01241 #endif 01242 return (float)(sum / norm); 01243 }
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 527 of file cmp.h.
References EMAN::EMObject::FLOAT, EMAN::EMObject::INT, and EMAN::TypeDict::put().
00528 { 00529 TypeDict d; 00530 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)"); 00531 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)"); 00532 d.put("ampweight", EMObject::INT, "If set, the amplitude of 'with' will be used as a weight in the averaging'. (default=0)"); 00533 d.put("zeromask", EMObject::INT, "Treat regions in either image that are zero as a mask"); 00534 d.put("minres", EMObject::FLOAT, "Lowest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=500"); 00535 d.put("maxres", EMObject::FLOAT, "Highest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=10"); 00536 return d; 00537 }
static Cmp* EMAN::PhaseCmp::NEW | ( | ) | [inline, static] |
const string PhaseCmp::NAME = "phase" [static] |