EMAN::PhaseCmp Class Reference

Amplitude weighted mean phase difference (radians) with optional SNR weight. More...

#include <cmp.h>

Inheritance diagram for EMAN::PhaseCmp:

Inheritance graph
[legend]
Collaboration diagram for EMAN::PhaseCmp:

Collaboration graph
[legend]
List of all members.

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 CmpNEW ()

Static Public Attributes

static const string NAME = "phase"

Detailed Description

Amplitude weighted mean phase difference (radians) with optional SNR weight.

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.


Member Function Documentation

float PhaseCmp::cmp ( EMData image,
EMData with 
) const [virtual]

To compare 'image' with another image passed in through its parameters.

An optional transformation may be used to transform the 2 images.

Parameters:
image The first image to be compared.
with The second image to be comppared.
Returns:
The comparison result. Smaller better by default

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]

Implements EMAN::Cmp.

Definition at line 517 of file cmp.h.

00518                 {
00519                         return "Mean phase difference";
00520                 }

string EMAN::PhaseCmp::get_name (  )  const [inline, virtual]

Get the Cmp's name.

Each Cmp is identified by a unique name.

Returns:
The Cmp's name.

Implements EMAN::Cmp.

Definition at line 512 of file cmp.h.

References NAME.

00513                 {
00514                         return NAME;
00515                 }

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.

Returns:
A dictionary containing the parameter info.

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]

Definition at line 522 of file cmp.h.

00523                 {
00524                         return new PhaseCmp();
00525                 }


Member Data Documentation

const string PhaseCmp::NAME = "phase" [static]

Definition at line 539 of file cmp.h.

Referenced by get_name().


The documentation for this class was generated from the following files:
Generated on Thu May 3 10:09:04 2012 for EMAN2 by  doxygen 1.4.7