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