Amplitude weighted mean phase difference (radians) with optional SNR weight. More...

#include <cmp.h>

Inheritance diagram for EMAN::PhaseCmp:

Collaboration diagram for EMAN::PhaseCmp:

[legend]

Public Member Functions
float	cmp (EMData image, EMData with) const
	To compare 'image' with another image passed in through its parameters. More...

string	get_name () const
	Get the Cmp's name. More...

string	get_desc () const

TypeDict	get_param_types () const
	Get Cmp parameter information in a dictionary. More...

Public Member Functions inherited from EMAN::Cmp
virtual	~Cmp ()

virtual Dict	get_params () const
	Get the Cmp parameters in a key/value dictionary. More...

virtual void	set_params (const Dict &new_params)
	Set the Cmp parameters using a key/value dictionary. More...

Static Public Member Functions
static Cmp *	NEW ()

Static Public Attributes
static const string	NAME = "phase"

Additional Inherited Members
Protected Member Functions inherited from EMAN::Cmp
void	validate_input_args (const EMData image, const EMData with) const

Protected Attributes inherited from EMAN::Cmp
Dict	params

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

Member Function Documentation

◆ cmp()

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

{
        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;
        
        float pmin = params.set_default("pmin",0.0f);
        float pmax = params.set_default("pmax",0.0f);
        if (minres>0 && pmin==0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;             //cutoff in pixels, assume square
//      else pmin=0;
        if (maxres>0 && pmax==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);
}

References EMAN::Util::angle_err_ri(), EMAN::Ctf::apix, EMAN::Ctf::compute_1d(), EMAN::Ctf::CTF_SNR, ENTERFUNC, EXITFUNC, EMAN::Util::hypot3(), EMAN::Util::hypot_fast_int(), InvalidCallException, EMAN::Cmp::params, EMAN::Dict::set_default(), x, and y.

◆ get_desc()

string EMAN::PhaseCmp::get_desc ( ) const

inlinevirtual

Implements EMAN::Cmp.

Definition at line 646 of file cmp.h.

                {
                        return "Mean phase difference";
                }

◆ get_name()

string EMAN::PhaseCmp::get_name ( ) const

inlinevirtual

Get the Cmp's name.

Each Cmp is identified by a unique name.

Returns: The Cmp's name.

Implements EMAN::Cmp.

Definition at line 641 of file cmp.h.

                {
                        return NAME;
                }

References NAME.

◆ get_param_types()

TypeDict EMAN::PhaseCmp::get_param_types ( ) const

inlinevirtual

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

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

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

◆ NEW()

static Cmp * EMAN::PhaseCmp::NEW ( )

inlinestatic

Definition at line 651 of file cmp.h.

                {
                        return new PhaseCmp();
                }

Member Data Documentation

◆ NAME

const string PhaseCmp::NAME = "phase"

static

Definition at line 668 of file cmp.h.

Referenced by get_name().

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

Public Member Functions

Static Public Member Functions

Static Public Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ cmp()

◆ get_desc()

◆ get_name()

◆ get_param_types()

◆ NEW()

Member Data Documentation

◆ NAME