EMAN::FRCCmp Class Reference

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

#include <cmp.h>

Inheritance diagram for EMAN::FRCCmp:

Inheritance graph
[legend]
Collaboration diagram for EMAN::FRCCmp:

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

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


Member Function Documentation

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

References EMAN::Ctf::apix, EMAN::EMData::calc_fourier_shell_correlation(), EMAN::EMData::calc_radial_dist(), EMAN::Ctf::compute_1d(), EMAN::EMData::copy(), EMAN::Ctf::CTF_SNR, EMAN::EMData::do_fft(), EMAN::EMData::do_fft_inplace(), ENTERFUNC, EXITFUNC, EMAN::EMObject::f, EMAN::EMData::get_attr(), EMAN::EMData::get_attr_default(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Util::goodf(), EMAN::EMData::has_attr(), InvalidCallException, EMAN::EMData::is_complex(), norm(), ny, EMAN::Cmp::params, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::update(), EMAN::Cmp::validate_input_args(), and weight.

01249 {
01250         ENTERFUNC;
01251         validate_input_args(image, with);
01252 
01253         int snrweight = params.set_default("snrweight", 0);
01254         int ampweight = params.set_default("ampweight", 0);
01255         int sweight = params.set_default("sweight", 1);
01256         int nweight = params.set_default("nweight", 0);
01257         int zeromask = params.set_default("zeromask",0);
01258         float minres = params.set_default("minres",500.0f);
01259         float maxres = params.set_default("maxres",10.0f);
01260 
01261         vector < float >fsc;
01262         bool use_cpu = true;
01263 
01264         if (use_cpu) {
01265                 if (zeromask) {
01266                         image=image->copy();
01267                         with=with->copy();
01268                 
01269                         int sz=image->get_xsize()*image->get_ysize()*image->get_zsize();
01270                         float *d1=image->get_data();
01271                         float *d2=with->get_data();
01272                 
01273                         for (int i=0; i<sz; i++) {
01274                                 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
01275                         }
01276                 
01277                         image->update();
01278                         with->update();
01279                         image->do_fft_inplace();
01280                         with->do_fft_inplace();
01281                         image->set_attr("free_me",1); 
01282                         with->set_attr("free_me",1); 
01283                 }
01284 
01285 
01286                 if (!image->is_complex()) {
01287                         image=image->do_fft(); 
01288                         image->set_attr("free_me",1); 
01289                 }
01290                 if (!with->is_complex()) { 
01291                         with=with->do_fft(); 
01292                         with->set_attr("free_me",1); 
01293                 }
01294 
01295                 fsc = image->calc_fourier_shell_correlation(with,1);
01296         }
01297         
01298         int ny = image->get_ysize();
01299         int ny2=ny/2+1;
01300 
01301         // The fast hypot here was supposed to speed things up. Little effect
01302 //      if (image->get_zsize()>1) fsc = image->calc_fourier_shell_correlation(with,1);
01303 //      else {
01304 //              double *sxy = (double *)malloc(ny2*sizeof(double)*4);
01305 //              double *sxx = sxy+ny2;
01306 //              double *syy = sxy+2*ny2;
01307 //              double *norm= sxy+3*ny2;
01308 //
01309 //              float *df1=image->get_data();
01310 //              float *df2=with->get_data();
01311 //              int nx2=image->get_xsize();
01312 //
01313 //              for (int y=-ny/2; y<ny/2; y++) {
01314 //                      for (int x=0; x<nx2/2; x++) {
01315 //                              if (x==0 && y<0) continue;      // skip Friedel pair
01316 //                              short r=Util::hypot_fast_int(x,y);
01317 //                              if (r>ny2-1) continue;
01318 //                              int l=x*2+(y<0?ny+y:y)*nx2;
01319 //                              sxy[r]+=df1[l]*df2[l]+df1[l+1]*df2[l+1];
01320 //                              sxx[r]+=df1[l]*df1[l];
01321 //                              syy[r]+=df2[l]*df2[l];
01322 //                              norm[r]+=1.0;
01323 //                      }
01324 //              }
01325 //              fsc.resize(ny2*3);
01326 //              for (int r=0; r<ny2; r++) {
01327 //                      fsc[r]=r*0.5/ny2;
01328 //                      fsc[ny2+r]=sxy[r]/(sqrt(sxx[r])*sqrt(syy[r]));
01329 //                      fsc[ny2*2+r]=norm[r];
01330 //              }
01331 //              free(sxy);
01332 //      }
01333 
01334         vector<float> snr;
01335         if (snrweight) {
01336                 Ctf *ctf = NULL;
01337                 if (!image->has_attr("ctf")) {
01338                         if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
01339                         ctf=with->get_attr("ctf");
01340                 }
01341                 else ctf=image->get_attr("ctf");
01342 
01343                 float ds=1.0f/(ctf->apix*ny);
01344                 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR);
01345                 for (int i=0; i<snr.size(); i++) {
01346                         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
01347                 }
01348                 if(ctf) {delete ctf; ctf=0;}
01349         }
01350 
01351         vector<float> amp;
01352         if (ampweight) amp=image->calc_radial_dist(ny/2,0,1,0);
01353 
01354         // Min/max modifications to weighting
01355         float pmin,pmax;
01356         if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;                //cutoff in pixels, assume square
01357         else pmin=0;
01358         if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
01359         else pmax=0;
01360 
01361         double sum=0.0, norm=0.0;
01362 
01363         for (int i=0; i<ny/2; i++) {
01364                 double weight=1.0;
01365                 if (sweight) weight*=fsc[(ny2)*2+i];
01366                 if (ampweight) weight*=amp[i];
01367                 if (snrweight) weight*=snr[i];
01368 //              if (snrweight)  {
01369 //                      if (snr[i]>0) weight*=sqrt(snr[i]);
01370 //                      else weight=0;
01371 //              }
01372 //if(snr[i]<0) printf("snr[%d] = %1.5g\n",i,snr[i]);
01373                 if (pmin>0) weight*=(tanh(5.0*(i-pmin)/pmin)+1.0)/2.0;
01374                 if (pmax>0) weight*=(1.0-tanh(i-pmax))/2.0;
01375                 
01376                 sum+=weight*fsc[ny2+i];
01377                 norm+=weight;
01378 //              printf("%d\t%f\t%f\n",i,weight,fsc[ny/2+1+i]);
01379         }
01380 
01381         // This performs a weighting that tries to normalize FRC by correcting from the number of particles represented by the average
01382         sum/=norm;
01383         if (nweight && with->get_attr_default("ptcl_repr",0) && sum>=0 && sum<1.0) {
01384                 sum=sum/(1.0-sum);                                                      // convert to SNR
01385                 sum/=(float)with->get_attr_default("ptcl_repr",0);      // divide by ptcl represented
01386                 sum=sum/(1.0+sum);                                                      // convert back to correlation
01387         }
01388 
01389         if (image->has_attr("free_me")) delete image;
01390         if (with->has_attr("free_me")) delete with;
01391 
01392         EXITFUNC;
01393 
01394         if (!Util::goodf(&sum)) sum=-3.0e38;
01395 
01396         //.Note the negative! This is because EMAN2 follows the convention that
01397         // smaller return values from comparitors indicate higher similarity -
01398         // this enables comparitors to be used in a generic fashion.
01399         return (float)-sum;
01400 }

string EMAN::FRCCmp::get_desc (  )  const [inline, virtual]

Implements EMAN::Cmp.

Definition at line 562 of file cmp.h.

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

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

References NAME.

00558                 {
00559                         return NAME;
00560                 }

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

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

00573                 {
00574                         TypeDict d;
00575                         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)");
00576                         d.put("ampweight", EMObject::INT, "If set, the amplitude of 'this' will be used to weight the result (default=0)");
00577                         d.put("sweight", EMObject::INT, "If set, weight the (1-D) average by the number of pixels in each ring (default=1)");
00578                         d.put("nweight", EMObject::INT, "Downweight similarity based on number of particles in reference (default=0)");
00579                         d.put("zeromask", EMObject::INT, "Treat regions in either image that are zero as a mask");
00580                         d.put("minres", EMObject::FLOAT, "Lowest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=500");
00581                         d.put("maxres", EMObject::FLOAT, "Highest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables.  Default=10");
00582                         return d;
00583                 }

static Cmp* EMAN::FRCCmp::NEW (  )  [inline, static]

Definition at line 567 of file cmp.h.

00568                 {
00569                         return new FRCCmp();
00570                 }


Member Data Documentation

const string FRCCmp::NAME = "frc" [static]

Definition at line 585 of file cmp.h.

Referenced by get_name().


The documentation for this class was generated from the following files:
Generated on Tue Jun 11 12:43:00 2013 for EMAN2 by  doxygen 1.4.7