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 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 930 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.

00931 {
00932         ENTERFUNC;
00933 
00934         int snrweight = params.set_default("snrweight", 0);
00935         int snrfn = params.set_default("snrfn",0);
00936         int ampweight = params.set_default("ampweight",0);
00937         int zeromask = params.set_default("zeromask",0);
00938         float minres = params.set_default("minres",500.0f);
00939         float maxres = params.set_default("maxres",2.0f);
00940 
00941         if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator");
00942 
00943         EMData *image_fft = NULL;
00944         EMData *with_fft = NULL;
00945 
00946         int ny = image->get_ysize();
00947 //      int np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2;
00948         int np = 0;
00949         vector<float> snr;
00950 
00951         // weighting based on SNR estimate from CTF
00952         if (snrweight) {
00953                 Ctf *ctf = NULL;
00954                 if (!image->has_attr("ctf")) {
00955                         if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
00956                         ctf=with->get_attr("ctf");
00957                 }
00958                 else ctf=image->get_attr("ctf");
00959 
00960                 float ds=1.0f/(ctf->apix*ny);
00961                 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR); // note that this returns ny/2 values
00962                 if(ctf) {delete ctf; ctf=0;}
00963                 np=snr.size();
00964         }
00965         // weighting based on empirical SNR function (not really good)
00966         else if (snrfn==1) {
00967                 np=ny/2;
00968                 snr.resize(np);
00969 
00970                 for (int i = 0; i < np; i++) {
00971                         float x2 = 10.0f*i/np;
00972                         snr[i] = x2 * exp(-x2);
00973                 }
00974         }
00975         // no weighting
00976         else {
00977                 np=ny/2;
00978                 snr.resize(np);
00979 
00980                 for (int i = 0; i < np; i++) snr[i]=1.0;
00981         }
00982 
00983         // Min/max modifications to weighting
00984         float pmin,pmax;
00985         if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;                //cutoff in pixels, assume square
00986         else pmin=0;
00987         if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
00988         else pmax=0;
00989 
00990 //      printf("%f\t%f\n",pmin,pmax);
00991 
00992         // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation
00993         for (int i = 0; i < np; i++) {
00994                 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f;
00995                 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f;
00996 //              printf("%d\t%f\n",i,snr[i]);
00997         }
00998 
00999         if (zeromask) {
01000                 image_fft=image->copy();
01001                 with_fft=with->copy();
01002                 
01003                 if (image_fft->is_complex()) image_fft->do_ift_inplace();
01004                 if (with_fft->is_complex()) with_fft->do_ift_inplace();
01005                 
01006                 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize();
01007                 float *d1=image_fft->get_data();
01008                 float *d2=with_fft->get_data();
01009                 
01010                 for (int i=0; i<sz; i++) {
01011                         if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
01012                 }
01013                 
01014                 image_fft->update();
01015                 with_fft->update();
01016                 image_fft->do_fft_inplace();
01017                 with_fft->do_fft_inplace();
01018                 image_fft->set_attr("free_me",1); 
01019                 with_fft->set_attr("free_me",1); 
01020         }
01021         else {
01022                 if (image->is_complex()) image_fft=image;
01023                 else {
01024                         image_fft=image->do_fft();
01025                         image_fft->set_attr("free_me",1);
01026                 }
01027                 
01028                 if (with->is_complex()) with_fft=with;
01029                 else {
01030                         with_fft=with->do_fft();
01031                         with_fft->set_attr("free_me",1);
01032                 }
01033         }
01034 //      image_fft->ri2ap();
01035 //      with_fft->ri2ap();
01036 
01037         const float *const image_fft_data = image_fft->get_const_data();
01038         const float *const with_fft_data = with_fft->get_const_data();
01039         double sum = 0;
01040         double norm = FLT_MIN;
01041         size_t i = 0;
01042 //      int nx=image_fft->get_xsize();
01043             ny=image_fft->get_ysize();
01044         int nz=image_fft->get_zsize();
01045         int nx2=image_fft->get_xsize()/2;
01046         int ny2=image_fft->get_ysize()/2;
01047 //      int nz2=image_fft->get_zsize()/2;
01048 
01049         // This can never happen any more...
01050         if (np==0) {
01051                 for (int z = 0; z < nz; z++){
01052                         for (int y = 0; y < ny; y++) {
01053                                 for (int x = 0; x < nx2; x ++) {
01054                                         float a;
01055                                         if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01056                                         else a=1.0f;
01057                                         sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01058                                         norm += a;
01059                                         i += 2;
01060                                 }
01061                         }
01062                 }
01063                 
01064         }
01065         else if (nz==1) {
01066                 for (int y = 0; y < ny; y++) {
01067                         for (int x = 0; x < nx2; x ++) {
01068                                 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y);
01069                                 if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
01070                                 
01071                                 float a;
01072                                 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01073                                 else a=1.0f;
01074                                 a*=snr[r];
01075                                 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01076                                 norm += a;
01077                                 i += 2;
01078                         }
01079                 }
01080         }
01081         else {
01082                 for (int z = 0; z < nz; z++){
01083                         for (int y = 0; y < ny; y++) {
01084                                 for (int x = 0; x < nx2; x ++) {
01085                                         int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z);
01086                                         if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
01087                                         
01088                                         float a;
01089                                         if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01090                                         else a=1.0f;
01091                                         a*=snr[r];
01092                                         sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01093                                         norm += a;
01094                                         i += 2;
01095                                 } 
01096                         }
01097                 }
01098                 
01099         }
01100 
01101         EXITFUNC;
01102 
01103         if( image_fft->has_attr("free_me") )
01104         {
01105                 delete image_fft;
01106                 image_fft = 0;
01107         }
01108         if( with_fft->has_attr("free_me") )
01109         {
01110                 delete with_fft;
01111                 with_fft = 0;
01112         }
01113 #if 0
01114         return (1.0f - sum / norm);
01115 #endif
01116         return (float)(sum / norm);
01117 }

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.

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

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


The documentation for this class was generated from the following files:
Generated on Thu Nov 17 12:45:27 2011 for EMAN2 by  doxygen 1.3.9.1