#include <cmp.h>


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