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:

[legend]
Collaboration diagram for EMAN::PhaseCmp:
[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 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 759 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(), nx, ny, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::update(), x, and y.

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

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.

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

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


The documentation for this class was generated from the following files:
Generated on Fri Apr 30 15:39:14 2010 for EMAN2 by  doxygen 1.3.9.1