EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes
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 551 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:
imageThe first image to be compared.
withThe second image to be comppared.
Returns:
The comparison result. Smaller better by default

Implements EMAN::Cmp.

Definition at line 1058 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.

{
        ENTERFUNC;

        int snrweight = params.set_default("snrweight", 0);
        int snrfn = params.set_default("snrfn",0);
        int ampweight = params.set_default("ampweight",0);
        int zeromask = params.set_default("zeromask",0);
        float minres = params.set_default("minres",500.0f);
        float maxres = params.set_default("maxres",10.0f);

        if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator");

        EMData *image_fft = NULL;
        EMData *with_fft = NULL;

        int ny = image->get_ysize();
//      int np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2;
        int np = 0;
        vector<float> snr;

        // weighting based on SNR estimate from CTF
        if (snrweight) {
                Ctf *ctf = NULL;
                if (!image->has_attr("ctf")) {
                        if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
                        ctf=with->get_attr("ctf");
                }
                else ctf=image->get_attr("ctf");

                float ds=1.0f/(ctf->apix*ny);
                snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR); // note that this returns ny/2 values
                for (int i=0; i<snr.size(); i++) {
                        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
                }
                if(ctf) {delete ctf; ctf=0;}
                np=snr.size();
        }
        // weighting based on empirical SNR function (not really good)
        else if (snrfn==1) {
                np=ny/2;
                snr.resize(np);

                for (int i = 0; i < np; i++) {
                        float x2 = 10.0f*i/np;
                        snr[i] = x2 * exp(-x2);
                }
        }
        // no weighting
        else {
                np=ny/2;
                snr.resize(np);

                for (int i = 0; i < np; i++) snr[i]=1.0;
        }

        // Min/max modifications to weighting
        float pmin,pmax;
        if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;                //cutoff in pixels, assume square
        else pmin=0;
        if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
        else pmax=0;

//      printf("%f\t%f\n",pmin,pmax);

        // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation
        for (int i = 0; i < np; i++) {
                if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f;
                if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f;
//              printf("%d\t%f\n",i,snr[i]);
        }

        if (zeromask) {
                image_fft=image->copy();
                with_fft=with->copy();
                
                if (image_fft->is_complex()) image_fft->do_ift_inplace();
                if (with_fft->is_complex()) with_fft->do_ift_inplace();
                
                int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize();
                float *d1=image_fft->get_data();
                float *d2=with_fft->get_data();
                
                for (int i=0; i<sz; i++) {
                        if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
                }
                
                image_fft->update();
                with_fft->update();
                image_fft->do_fft_inplace();
                with_fft->do_fft_inplace();
                image_fft->set_attr("free_me",1); 
                with_fft->set_attr("free_me",1); 
        }
        else {
                if (image->is_complex()) image_fft=image;
                else {
                        image_fft=image->do_fft();
                        image_fft->set_attr("free_me",1);
                }
                
                if (with->is_complex()) with_fft=with;
                else {
                        with_fft=with->do_fft();
                        with_fft->set_attr("free_me",1);
                }
        }
//      image_fft->ri2ap();
//      with_fft->ri2ap();

        const float *const image_fft_data = image_fft->get_const_data();
        const float *const with_fft_data = with_fft->get_const_data();
        double sum = 0;
        double norm = FLT_MIN;
        size_t i = 0;
//      int nx=image_fft->get_xsize();
            ny=image_fft->get_ysize();
        int nz=image_fft->get_zsize();
        int nx2=image_fft->get_xsize()/2;
        int ny2=image_fft->get_ysize()/2;
//      int nz2=image_fft->get_zsize()/2;

        // This can never happen any more...
        if (np==0) {
                for (int z = 0; z < nz; z++){
                        for (int y = 0; y < ny; y++) {
                                for (int x = 0; x < nx2; x ++) {
                                        float a;
                                        if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
                                        else a=1.0f;
                                        sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
                                        norm += a;
                                        i += 2;
                                }
                        }
                }
                
        }
        else if (nz==1) {
                for (int y = 0; y < ny; y++) {
                        for (int x = 0; x < nx2; x ++) {
                                int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y);
                                if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
                                
                                float a;
                                if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
                                else a=1.0f;
                                a*=snr[r];
                                sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
                                norm += a;
                                i += 2;
                        }
                }
        }
        else {
                for (int z = 0; z < nz; z++){
                        for (int y = 0; y < ny; y++) {
                                for (int x = 0; x < nx2; x ++) {
                                        int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z);
                                        if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
                                        
                                        float a;
                                        if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
                                        else a=1.0f;
                                        a*=snr[r];
                                        sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
                                        norm += a;
                                        i += 2;
                                } 
                        }
                }
                
        }

        EXITFUNC;

        if( image_fft->has_attr("free_me") )
        {
                delete image_fft;
                image_fft = 0;
        }
        if( with_fft->has_attr("free_me") )
        {
                delete with_fft;
                with_fft = 0;
        }
#if 0
        return (1.0f - sum / norm);
#endif
        return (float)(sum / norm);
}
string EMAN::PhaseCmp::get_desc ( ) const [inline, virtual]

Implements EMAN::Cmp.

Definition at line 561 of file cmp.h.

                {
                        return "Mean phase difference";
                }
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 556 of file cmp.h.

References NAME.

                {
                        return NAME;
                }
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 571 of file cmp.h.

References EMAN::EMObject::FLOAT, EMAN::EMObject::INT, and EMAN::TypeDict::put().

                {
                        TypeDict d;
                        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)");
                        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)");
                        d.put("ampweight", EMObject::INT, "If set, the amplitude of 'with' will be used as a weight in the averaging'. (default=0)");
                        d.put("zeromask", EMObject::INT, "Treat regions in either image that are zero as a mask");
                        d.put("minres", EMObject::FLOAT, "Lowest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=500");
                        d.put("maxres", EMObject::FLOAT, "Highest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables.  Default=10");
                        return d;
                }
static Cmp* EMAN::PhaseCmp::NEW ( ) [inline, static]

Definition at line 566 of file cmp.h.

                {
                        return new PhaseCmp();
                }

Member Data Documentation

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

Definition at line 583 of file cmp.h.

Referenced by get_name().


The documentation for this class was generated from the following files: