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 463 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 929 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.

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 }

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

Implements EMAN::Cmp.

Definition at line 473 of file cmp.h.

00474                 {
00475                         return "Mean phase difference";
00476                 }

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

References NAME.

00469                 {
00470                         return NAME;
00471                 }

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

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

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

Definition at line 478 of file cmp.h.

00479                 {
00480                         return new PhaseCmp();
00481                 }


Member Data Documentation

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

Definition at line 495 of file cmp.h.

Referenced by get_name().


The documentation for this class was generated from the following files:
Generated on Tue Jul 12 13:48:09 2011 for EMAN2 by  doxygen 1.4.7