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

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 418 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 762 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.

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

string EMAN::PhaseCmp::get_desc  )  const [inline, virtual]
 

Implements EMAN::Cmp.

Definition at line 428 of file cmp.h.

00429                 {
00430                         return "Mean phase difference";
00431                 }

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 423 of file cmp.h.

References NAME.

00424                 {
00425                         return NAME;
00426                 }

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 438 of file cmp.h.

References EMAN::EMObject::FLOAT, EMAN::EMObject::INT, and 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                 }

static Cmp* EMAN::PhaseCmp::NEW  )  [inline, static]
 

Definition at line 433 of file cmp.h.

00434                 {
00435                         return new PhaseCmp();
00436                 }


Member Data Documentation

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

Definition at line 450 of file cmp.h.

Referenced by get_name().


The documentation for this class was generated from the following files:
Generated on Tue May 25 17:36:23 2010 for EMAN2 by  doxygen 1.4.4