Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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

CmpNEW ()

Static Public Attributes

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::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 }

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.

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::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                 }

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 55 of file cmp.cpp.


The documentation for this class was generated from the following files:
Generated on Tue Jun 11 13:47:58 2013 for EMAN2 by  doxygen 1.3.9.1