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
static Cmp *	NEW ()
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 417 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 733 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.

{
        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");

#ifdef EMAN2_USING_CUDA
        if (image->gpu_operation_preferred()) {
//              cout << "Cuda cmp" << endl;
                EXITFUNC;
                return cuda_cmp(image,with);
        }
#endif

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

References NAME.

                {
                        return NAME;
                }

string EMAN::PhaseCmp::get_desc

(

)

const [inline, virtual]

Implements EMAN::Cmp.

Definition at line 427 of file cmp.h.

                {
                        return "Mean phase difference";
                }

static Cmp* EMAN::PhaseCmp::NEW

(

)

[inline, static]

Definition at line 432 of file cmp.h.

                {
                        return new PhaseCmp();
                }

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 437 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;
                }

---

Member Data Documentation

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

Definition at line 449 of file cmp.h.

Referenced by get_name().

---

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

• cmp.h
• cmp.cpp