#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 418 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 759 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(), nx, ny, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::update(), x, and y. 00760 { 00761 ENTERFUNC; 00762 00763 int snrweight = params.set_default("snrweight", 0); 00764 int snrfn = params.set_default("snrfn",0); 00765 int ampweight = params.set_default("ampweight",0); 00766 int zeromask = params.set_default("zeromask",0); 00767 float minres = params.set_default("minres",500.0f); 00768 float maxres = params.set_default("maxres",10.0f); 00769 00770 if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator"); 00771 00772 #ifdef EMAN2_USING_CUDA 00773 if (image->gpu_operation_preferred()) { 00774 // cout << "Cuda cmp" << endl; 00775 EXITFUNC; 00776 return cuda_cmp(image,with); 00777 } 00778 #endif 00779 00780 EMData *image_fft = NULL; 00781 EMData *with_fft = NULL; 00782 00783 int ny = image->get_ysize(); 00784 // int np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2; 00785 int np = 0; 00786 vector<float> snr; 00787 00788 // weighting based on SNR estimate from CTF 00789 if (snrweight) { 00790 Ctf *ctf = NULL; 00791 if (!image->has_attr("ctf")) { 00792 if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters"); 00793 ctf=with->get_attr("ctf"); 00794 } 00795 else ctf=image->get_attr("ctf"); 00796 00797 float ds=1.0f/(ctf->apix*ny); 00798 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR); // note that this returns ny/2 values 00799 if(ctf) {delete ctf; ctf=0;} 00800 np=snr.size(); 00801 } 00802 // weighting based on empirical SNR function (not really good) 00803 else if (snrfn==1) { 00804 np=ny/2; 00805 snr.resize(np); 00806 00807 for (int i = 0; i < np; i++) { 00808 float x2 = 10.0f*i/np; 00809 snr[i] = x2 * exp(-x2); 00810 } 00811 } 00812 // no weighting 00813 else { 00814 np=ny/2; 00815 snr.resize(np); 00816 00817 for (int i = 0; i < np; i++) snr[i]=1.0; 00818 } 00819 00820 // Min/max modifications to weighting 00821 float pmin,pmax; 00822 if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres; //cutoff in pixels, assume square 00823 else pmin=0; 00824 if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres; 00825 else pmax=0; 00826 00827 // printf("%f\t%f\n",pmin,pmax); 00828 00829 // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation 00830 for (int i = 0; i < np; i++) { 00831 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f; 00832 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f; 00833 // printf("%d\t%f\n",i,snr[i]); 00834 } 00835 00836 if (zeromask) { 00837 image_fft=image->copy(); 00838 with_fft=with->copy(); 00839 00840 if (image_fft->is_complex()) image_fft->do_ift_inplace(); 00841 if (with_fft->is_complex()) with_fft->do_ift_inplace(); 00842 00843 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize(); 00844 float *d1=image_fft->get_data(); 00845 float *d2=with_fft->get_data(); 00846 00847 for (int i=0; i<sz; i++) { 00848 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; } 00849 } 00850 00851 image_fft->update(); 00852 with_fft->update(); 00853 image_fft->do_fft_inplace(); 00854 with_fft->do_fft_inplace(); 00855 image_fft->set_attr("free_me",1); 00856 with_fft->set_attr("free_me",1); 00857 } 00858 else { 00859 if (image->is_complex()) image_fft=image; 00860 else { 00861 image_fft=image->do_fft(); 00862 image_fft->set_attr("free_me",1); 00863 } 00864 00865 if (with->is_complex()) with_fft=with; 00866 else { 00867 with_fft=with->do_fft(); 00868 with_fft->set_attr("free_me",1); 00869 } 00870 } 00871 // image_fft->ri2ap(); 00872 // with_fft->ri2ap(); 00873 00874 const float *const image_fft_data = image_fft->get_const_data(); 00875 const float *const with_fft_data = with_fft->get_const_data(); 00876 double sum = 0; 00877 double norm = FLT_MIN; 00878 size_t i = 0; 00879 int nx=image_fft->get_xsize(); 00880 ny=image_fft->get_ysize(); 00881 int nz=image_fft->get_zsize(); 00882 int nx2=image_fft->get_xsize()/2; 00883 int ny2=image_fft->get_ysize()/2; 00884 int nz2=image_fft->get_zsize()/2; 00885 00886 // This can never happen any more... 00887 if (np==0) { 00888 for (int z = 0; z < nz; z++){ 00889 for (int y = 0; y < ny; y++) { 00890 for (int x = 0; x < nx2; x ++) { 00891 float a; 00892 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 00893 else a=1.0f; 00894 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 00895 norm += a; 00896 i += 2; 00897 } 00898 } 00899 } 00900 00901 } 00902 else if (nz==1) { 00903 for (int y = 0; y < ny; y++) { 00904 for (int x = 0; x < nx2; x ++) { 00905 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y); 00906 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 00907 00908 float a; 00909 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 00910 else a=1.0f; 00911 a*=snr[r]; 00912 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 00913 norm += a; 00914 i += 2; 00915 } 00916 } 00917 } 00918 else { 00919 for (int z = 0; z < nz; z++){ 00920 for (int y = 0; y < ny; y++) { 00921 for (int x = 0; x < nx2; x ++) { 00922 int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z); 00923 if (r>=ny2) { i+=2; continue; } // we only have snr values to the box radius 00924 00925 float a; 00926 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]); 00927 else a=1.0f; 00928 a*=snr[r]; 00929 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a; 00930 norm += a; 00931 i += 2; 00932 } 00933 } 00934 } 00935 00936 } 00937 00938 EXITFUNC; 00939 00940 if( image_fft->has_attr("free_me") ) 00941 { 00942 delete image_fft; 00943 image_fft = 0; 00944 } 00945 if( with_fft->has_attr("free_me") ) 00946 { 00947 delete with_fft; 00948 with_fft = 0; 00949 } 00950 #if 0 00951 return (1.0f - sum / norm); 00952 #endif 00953 return (float)(sum / norm); 00954 }
|
|
Implements EMAN::Cmp. Definition at line 428 of file cmp.h. 00429 { 00430 return "Mean phase difference"; 00431 }
|
|
Get the Cmp's name. Each Cmp is identified by a unique name.
Implements EMAN::Cmp. Definition at line 423 of file cmp.h. 00424 {
00425 return NAME;
00426 }
|
|
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 438 of file cmp.h. References EMAN::TypeDict::put(). 00439 { 00440 TypeDict d; 00441 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)"); 00442 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)"); 00443 d.put("ampweight", EMObject::INT, "If set, the amplitude of 'with' will be used as a weight in the averaging'. (default=0)"); 00444 d.put("zeromask", EMObject::INT, "Treat regions in either image that are zero as a mask"); 00445 d.put("minres", EMObject::FLOAT, "Lowest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=500"); 00446 d.put("maxres", EMObject::FLOAT, "Highest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=10"); 00447 return d; 00448 }
|
|
Definition at line 433 of file cmp.h. 00434 { 00435 return new PhaseCmp(); 00436 }
|
|
|