FRCCmp returns a quality factor based on FRC between images. More...

#include <cmp.h>

Inheritance diagram for EMAN::FRCCmp:

Collaboration diagram for EMAN::FRCCmp:

[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 = "frc"

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

FRCCmp returns a quality factor based on FRC between images.

Fourier ring correlation (FRC) is a measure of statistical dependency between two averages, computed by comparison of rings in Fourier space. 1 means prefect agreement. 0 means no correlation.

Definition at line 681 of file cmp.h.

Member Function Documentation

◆ cmp()

float FRCCmp::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 1417 of file cmp.cpp.

{
        ENTERFUNC;
        validate_input_args(image, with);
 
        int snrweight = params.set_default("snrweight", 0);
        int ampweight = params.set_default("ampweight", 0);
        int sweight = params.set_default("sweight", 1);
        int nweight = params.set_default("nweight", 0);
        int zeromask = params.set_default("zeromask",0);
        float minres = params.set_default("minres",200.0f);
        float maxres = params.set_default("maxres",8.0f);
 
        vector < float >fsc;
        bool use_cpu = true;
 
        if (use_cpu) {
                if (zeromask) {
                        image=image->copy();
                        with=with->copy();
                
                        int sz=image->get_xsize()*image->get_ysize()*image->get_zsize();
                        float *d1=image->get_data();
                        float *d2=with->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->update();
                        with->update();
                        image->do_fft_inplace();
                        with->do_fft_inplace();
                        image->set_attr("free_me",1); 
                        with->set_attr("free_me",1); 
                }
 
 
                if (!image->is_complex()) {
                        image=image->do_fft(); 
                        image->set_attr("free_me",1); 
                }
                if (!with->is_complex()) { 
                        with=with->do_fft(); 
                        with->set_attr("free_me",1); 
                }
 
                fsc = image->calc_fourier_shell_correlation(with,1);
        }
        
        int ny = image->get_ysize();
        int ny2=ny/2+1;
 
        // The fast hypot here was supposed to speed things up. Little effect
//      if (image->get_zsize()>1) fsc = image->calc_fourier_shell_correlation(with,1);
//      else {
//              double *sxy = (double *)malloc(ny2*sizeof(double)*4);
//              double *sxx = sxy+ny2;
//              double *syy = sxy+2*ny2;
//              double *norm= sxy+3*ny2;
//
//              float *df1=image->get_data();
//              float *df2=with->get_data();
//              int nx2=image->get_xsize();
//
//              for (int y=-ny/2; y<ny/2; y++) {
//                      for (int x=0; x<nx2/2; x++) {
//                              if (x==0 && y<0) continue;      // skip Friedel pair
//                              short r=Util::hypot_fast_int(x,y);
//                              if (r>ny2-1) continue;
//                              int l=x*2+(y<0?ny+y:y)*nx2;
//                              sxy[r]+=df1[l]*df2[l]+df1[l+1]*df2[l+1];
//                              sxx[r]+=df1[l]*df1[l];
//                              syy[r]+=df2[l]*df2[l];
//                              norm[r]+=1.0;
//                      }
//              }
//              fsc.resize(ny2*3);
//              for (int r=0; r<ny2; r++) {
//                      fsc[r]=r*0.5/ny2;
//                      fsc[ny2+r]=sxy[r]/(sqrt(sxx[r])*sqrt(syy[r]));
//                      fsc[ny2*2+r]=norm[r];
//              }
//              free(sxy);
//      }
 
        vector<float> snr;
        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);
                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;}
        }
 
        vector<float> amp;
        if (ampweight) amp=image->calc_radial_dist(ny/2,0,1,0);
 
        // Min/max modifications to weighting
        float pmin = params.set_default("pmin",0);
        float pmax = params.set_default("pmax", 0);
        
        if (pmin==0 && minres>0)
                pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;              //cutoff in pixels, assume square
        
        if (pmax==0 && maxres>0) 
                pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
        
 
        double sum=0.0, norm=0.0;
 
        for (int i=0; i<ny/2; i++) {
                double weight=1.0;
                if (sweight) weight*=fsc[(ny2)*2+i];
                if (ampweight) weight*=amp[i];
                if (snrweight) weight*=snr[i];
//              if (snrweight)  {
//                      if (snr[i]>0) weight*=sqrt(snr[i]);
//                      else weight=0;
//              }
//if(snr[i]<0) printf("snr[%d] = %1.5g\n",i,snr[i]);
                if (pmin>0) weight*=(tanh(5.0*(i-pmin)/pmin)+1.0)/2.0;
                if (pmax>0) weight*=(1.0-tanh(i-pmax))/2.0;
                
                sum+=weight*fsc[ny2+i];
                norm+=weight;
//              printf("%d\t%f\t%f\n",i,weight,fsc[ny/2+1+i]);
        }
 
        // This performs a weighting that tries to normalize FRC by correcting from the number of particles represented by the average
        sum/=norm;
        if (nweight && with->get_attr_default("ptcl_repr",0) && sum>=0 && sum<1.0) {
                sum=sum/(1.0-sum);                                                      // convert to SNR
                sum/=(float)with->get_attr_default("ptcl_repr",0);      // divide by ptcl represented
                sum=sum/(1.0+sum);                                                      // convert back to correlation
        }
 
        if (image->has_attr("free_me")) delete image;
        if (with->has_attr("free_me")) delete with;
 
        EXITFUNC;
 
        if (!Util::goodf(&sum)) sum=-2.0;       // normally should be >-1.0
 
        //.Note the negative! This is because EMAN2 follows the convention that
        // smaller return values from comparitors indicate higher similarity -
        // this enables comparitors to be used in a generic fashion.
        return (float)-sum;
}

References EMAN::Ctf::apix, EMAN::EMData::calc_radial_dist(), EMAN::Ctf::compute_1d(), EMAN::Ctf::CTF_SNR, ENTERFUNC, EXITFUNC, EMAN::Util::goodf(), InvalidCallException, EMAN::Cmp::params, EMAN::Dict::set_default(), and EMAN::Cmp::validate_input_args().

◆ get_desc()

string EMAN::FRCCmp::get_desc ( ) const

inlinevirtual

Implements EMAN::Cmp.

Definition at line 691 of file cmp.h.

                {
                        return "Computes the mean Fourier Ring Correlation between the image and reference (with optional weighting factors).";
                }

◆ get_name()

string EMAN::FRCCmp::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 686 of file cmp.h.

                {
                        return NAME;
                }

References NAME.

◆ get_param_types()

TypeDict EMAN::FRCCmp::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 701 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("ampweight", EMObject::INT, "If set, the amplitude of 'this' will be used to weight the result (default=0)");
                        d.put("sweight", EMObject::INT, "If set, weight the (1-D) average by the number of pixels in each ring (default=1)");
                        d.put("nweight", EMObject::INT, "Downweight similarity based on number of particles in reference (default=0)");
                        d.put("zeromask", EMObject::INT, "Treat regions in either image that are zero as a mask (default=0)");
                        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");
                        d.put("pmin", EMObject::FLOAT, "The minimum resolution in pixels.");
                        d.put("pmax", EMObject::FLOAT, "The maximum resolution in pixels.");
                        return d;
                }

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

◆ NEW()

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

inlinestatic

Definition at line 696 of file cmp.h.

                {
                        return new FRCCmp();
                }

Member Data Documentation

◆ NAME

const string FRCCmp::NAME = "frc"

static

Definition at line 716 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