#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 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.
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::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. 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",10.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 for (int i=0; i<snr.size(); i++) { 01089 if (snr[i]<=0) snr[i]=0.001; // make sure that points don't get completely excluded due to SNR estimation issues, or worse, contribute with a negative weight 01090 } 01091 if(ctf) {delete ctf; ctf=0;} 01092 np=snr.size(); 01093 } 01094 // weighting based on empirical SNR function (not really good) 01095 else if (snrfn==1) { 01096 np=ny/2; 01097 snr.resize(np); 01098 01099 for (int i = 0; i < np; i++) { 01100 float x2 = 10.0f*i/np; 01101 snr[i] = x2 * exp(-x2); 01102 } 01103 } 01104 // no weighting 01105 else { 01106 np=ny/2; 01107 snr.resize(np); 01108 01109 for (int i = 0; i < np; i++) snr[i]=1.0; 01110 } 01111 01112 // Min/max modifications to weighting 01113 float pmin,pmax; 01114 if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres; //cutoff in pixels, assume square 01115 else pmin=0; 01116 if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres; 01117 else pmax=0; 01118 01119 // printf("%f\t%f\n",pmin,pmax); 01120 01121 // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation 01122 for (int i = 0; i < np; i++) { 01123 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f; 01124 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f; 01125 // printf("%d\t%f\n",i,snr[i]); 01126 } 01127 01128 if (zeromask) { 01129 image_fft=image->copy(); 01130 with_fft=with->copy(); 01131 01132 if (image_fft->is_complex()) image_fft->do_ift_inplace(); 01133 if (with_fft->is_complex()) with_fft->do_ift_inplace(); 01134 01135 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize(); 01136 float *d1=image_fft->get_data(); 01137 float *d2=with_fft->get_data(); 01138 01139 for (int i=0; i<sz; i++) { 01140 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; } 01141 } 01142 01143 image_fft->update(); 01144 with_fft->update(); 01145 image_fft->do_fft_inplace(); 01146 with_fft->do_fft_inplace(); 01147 image_fft->set_attr("free_me",1); 01148 with_fft->set_attr("free_me",1); 01149 } 01150 else { 01151 if (image->is_complex()) image_fft=image; 01152 else { 01153 image_fft=image->do_fft(); 01154 image_fft->set_attr("free_me",1); 01155 } 01156 01157 if (with->is_complex()) with_fft=with; 01158 else { 01159 with_fft=with->do_fft(); 01160 with_fft->set_attr("free_me",1); 01161 } 01162 } 01163 // image_fft->ri2ap(); 01164 // with_fft->ri2ap(); 01165 01166 const float *const image_fft_data = image_fft->get_const_data(); 01167 const float *const with_fft_data = with_fft->get_const_data(); 01168 double sum = 0; 01169 double norm = FLT_MIN; 01170 size_t i = 0; 01171 // int nx=image_fft->get_xsize(); 01172 ny=image_fft->get_ysize(); 01173 int nz=image_fft->get_zsize(); 01174 int nx2=image_fft->get_xsize()/2; 01175 int ny2=image_fft->get_ysize()/2; 01176 // int nz2=image_fft->get_zsize()/2; 01177 01178 // This can never happen any more... 01179 if (np==0) { 01180 for (int z = 0; z < nz; z++){ 01181 for (int y = 0; y < ny; y++) { 01182 for (int x = 0; x < nx2; x ++) { 01183 float a; 01184 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01185 else a=1.0f; 01186 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01187 norm += a; 01188 i += 2; 01189 } 01190 } 01191 } 01192 01193 } 01194 else if (nz==1) { 01195 for (int y = 0; y < ny; y++) { 01196 for (int x = 0; x < nx2; x ++) { 01197 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y); 01198 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 01199 01200 float a; 01201 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01202 else a=1.0f; 01203 a*=snr[r]; 01204 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01205 norm += a; 01206 i += 2; 01207 } 01208 } 01209 } 01210 else { 01211 for (int z = 0; z < nz; z++){ 01212 for (int y = 0; y < ny; y++) { 01213 for (int x = 0; x < nx2; x ++) { 01214 int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z); 01215 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 01216 01217 float a; 01218 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 01219 else a=1.0f; 01220 a*=snr[r]; 01221 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 01222 norm += a; 01223 i += 2; 01224 } 01225 } 01226 } 01227 01228 } 01229 01230 EXITFUNC; 01231 01232 if( image_fft->has_attr("free_me") ) 01233 { 01234 delete image_fft; 01235 image_fft = 0; 01236 } 01237 if( with_fft->has_attr("free_me") ) 01238 { 01239 delete with_fft; 01240 with_fft = 0; 01241 } 01242 #if 0 01243 return (1.0f - sum / norm); 01244 #endif 01245 return (float)(sum / norm); 01246 }
|
|
Implements EMAN::Cmp. Definition at line 517 of file cmp.h. 00518 { 00519 return "Mean phase difference"; 00520 }
|
|
Get the Cmp's name. Each Cmp is identified by a unique name.
Implements EMAN::Cmp. Definition at line 512 of file cmp.h. 00513 {
00514 return NAME;
00515 }
|
|
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::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 }
|
|
Definition at line 522 of file cmp.h. 00523 { 00524 return new PhaseCmp(); 00525 }
|
|
|