EMAN::OptVarianceCmp Class Reference

Variance between two data sets after various modifications. More...

#include <cmp.h>

Inheritance diagram for EMAN::OptVarianceCmp:

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

Collaboration graph
[legend]
List of all members.

Public Member Functions

 OptVarianceCmp ()
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.
float get_scale () const
float get_shift () const

Static Public Member Functions

static CmpNEW ()

Static Public Attributes

static const string NAME = "optvariance"

Private Attributes

float scale
float shift

Detailed Description

Variance between two data sets after various modifications.

Generally, 'this' should be noisy and 'with' should be less noisy. linear scaling (mx + b) of the densities in 'this' is performed to produce the smallest possible variance between images.

If keepzero is set, then zero pixels are left at zero (b is not added). if matchfilt is set, then 'with' is filtered so its radial power spectrum matches 'this' If invert is set, then (y-b)/m is applied to the second image rather than mx+b to the first.

To modify 'this' to match the operation performed here, scale should be applied first, then b should be added

Definition at line 447 of file cmp.h.


Constructor & Destructor Documentation

EMAN::OptVarianceCmp::OptVarianceCmp (  )  [inline]

Definition at line 450 of file cmp.h.

Referenced by NEW().

00450 : scale(0), shift(0) {}


Member Function Documentation

float OptVarianceCmp::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 860 of file cmp.cpp.

References b, EMAN::Util::calc_least_square_fit(), EMAN::EMData::calc_radial_dist(), EMAN::EMData::do_fft(), ENTERFUNC, EXITFUNC, EMAN::EMData::get_const_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ImageDimensionException, nx, EMAN::Cmp::params, EMAN::EMData::ri2ap(), scale, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), shift, square, EMAN::Util::square(), EMAN::Cmp::validate_input_args(), EMAN::EMData::write_image(), x, and y.

00861 {
00862         ENTERFUNC;
00863         validate_input_args(image, with);
00864 
00865         int keepzero = params.set_default("keepzero", 1);
00866         int invert = params.set_default("invert",0);
00867         int matchfilt = params.set_default("matchfilt",1);
00868         int matchamp = params.set_default("matchamp",0);
00869         int radweight = params.set_default("radweight",0);
00870         int dbug = params.set_default("debug",0);
00871 
00872         size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00873 
00874 
00875         EMData *with2=NULL;
00876         if (matchfilt) {
00877                 EMData *a = image->do_fft();
00878                 EMData *b = with->do_fft();
00879 
00880                 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1);
00881                 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1);
00882 
00883                 float avg=0;
00884                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00885                         rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i]));
00886                         avg+=rfa[i];
00887                 }
00888 
00889                 avg/=a->get_ysize()/2.0f;
00890                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00891                         if (rfa[i]>avg*10.0) rfa[i]=10.0;                       // If some particular location has a small but non-zero value, we don't want to overcorrect it
00892                 }
00893                 rfa[0]=0.0;
00894 
00895                 if (dbug) b->write_image("a.hdf",-1);
00896 
00897                 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa);
00898                 with2=b->do_ift();
00899 
00900                 if (dbug) b->write_image("a.hdf",-1);
00901                 if (dbug) a->write_image("a.hdf",-1);
00902 
00903 /*              if (dbug) {
00904                         FILE *out=fopen("a.txt","w");
00905                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]);
00906                         fclose(out);
00907 
00908                         out=fopen("b.txt","w");
00909                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]);
00910                         fclose(out);
00911                 }*/
00912 
00913 
00914                 delete a;
00915                 delete b;
00916 
00917                 if (dbug) {
00918                         with2->write_image("a.hdf",-1);
00919                         image->write_image("a.hdf",-1);
00920                 }
00921 
00922 //              with2->process_inplace("matchfilt",Dict("to",this));
00923 //              x_data = with2->get_data();
00924         }
00925 
00926         // This applies the individual Fourier amplitudes from 'image' and
00927         // applies them to 'with'
00928         if (matchamp) {
00929                 EMData *a = image->do_fft();
00930                 EMData *b = with->do_fft();
00931                 size_t size2 = (size_t)a->get_xsize() * a->get_ysize() * a->get_zsize();
00932 
00933                 a->ri2ap();
00934                 b->ri2ap();
00935 
00936                 const float *const ad=a->get_const_data();
00937                 float * bd=b->get_data();
00938 
00939                 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i];
00940                 b->update();
00941 
00942                 b->ap2ri();
00943                 with2=b->do_ift();
00944 //with2->write_image("a.hdf",-1);
00945                 delete a;
00946                 delete b;
00947         }
00948 
00949         const float * x_data;
00950         if (with2) x_data=with2->get_const_data();
00951         else x_data = with->get_const_data();
00952         const float *const y_data = image->get_const_data();
00953 
00954         size_t nx = image->get_xsize();
00955         float m = 0;
00956         float b = 0;
00957 
00958         // This will write the x vs y file used to calculate the density
00959         // optimization. This behavior may change in the future
00960         if (dbug) {
00961                 FILE *out=fopen("dbug.optvar.txt","w");
00962                 if (out) {
00963                         for (size_t i=0; i<size; i++) {
00964                                 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]);
00965                         }
00966                         fclose(out);
00967                 }
00968         }
00969 
00970 
00971         Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero);
00972         if (m == 0) {
00973                 m = FLT_MIN;
00974         }
00975         b = -b / m;
00976         m = 1.0f / m;
00977 
00978         // While negative slopes are really not a valid comparison in most cases, we
00979         // still want to detect these instances, so this if is removed
00980 /*      if (m < 0) {
00981                 b = 0;
00982                 m = 1000.0;
00983         }*/
00984 
00985         double  result = 0;
00986         int count = 0;
00987 
00988         if (radweight) {
00989                 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only");
00990                 if (keepzero) {
00991                         for (size_t i = 0,y=0; i < size; y++) {
00992                                 for (size_t x=0; x<nx; i++,x++) {
00993                                         if (y_data[i] && x_data[i]) {
00994 #ifdef  _WIN32
00995                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00996                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00997 #else
00998                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00999                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
01000 #endif
01001                                                 count++;
01002                                         }
01003                                 }
01004                         }
01005                         result/=count;
01006                 }
01007                 else {
01008                         for (size_t i = 0,y=0; i < size; y++) {
01009                                 for (size_t x=0; x<nx; i++,x++) {
01010 #ifdef  _WIN32
01011                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
01012                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
01013 #else
01014                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
01015                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
01016 #endif
01017                                 }
01018                         }
01019                         result = result / size;
01020                 }
01021         }
01022         else {
01023                 if (keepzero) {
01024                         for (size_t i = 0; i < size; i++) {
01025                                 if (y_data[i] && x_data[i]) {
01026                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
01027                                         else result += Util::square((x_data[i] * m) + b - y_data[i]);
01028                                         count++;
01029                                 }
01030                         }
01031                         result/=count;
01032                 }
01033                 else {
01034                         for (size_t i = 0; i < size; i++) {
01035                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
01036                                 else result += Util::square((x_data[i] * m) + b - y_data[i]);
01037                         }
01038                         result = result / size;
01039                 }
01040         }
01041         scale = m;
01042         shift = b;
01043 
01044         image->set_attr("ovcmp_m",m);
01045         image->set_attr("ovcmp_b",b);
01046         if (with2) delete with2;
01047         EXITFUNC;
01048 
01049 #if 0
01050         return (1 - result);
01051 #endif
01052 
01053         return static_cast<float>(result);
01054 }

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

Implements EMAN::Cmp.

Definition at line 459 of file cmp.h.

00460                 {
00461                         return "Real-space variance after density optimization, self should be noisy and target less noisy. Linear transform applied to density to minimize variance.";
00462                 }

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

References NAME.

00455                 {
00456                         return NAME;
00457                 }

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

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

00470                 {
00471                         TypeDict d;
00472                         d.put("invert", EMObject::INT, "If set, 'with' is rescaled rather than 'this'. 'this' should still be the noisier image. (default=0)");
00473                         d.put("keepzero", EMObject::INT, "If set, zero pixels will not be adjusted in the linear density optimization. (default=1)");
00474                         d.put("matchfilt", EMObject::INT, "If set, with will be filtered so its radial power spectrum matches 'this' before density optimization of this. (default=1)");
00475                         d.put("matchamp", EMObject::INT, "Takes per-pixel Fourier amplitudes from self and imposes them on the target, but leaves the phases alone. (default=0)");
00476                         d.put("radweight", EMObject::INT, "Upweight variances closer to the edge of the image. (default=0)");
00477                         d.put("debug", EMObject::INT, "Performs various debugging actions if set.");
00478                         return d;
00479                 }

float EMAN::OptVarianceCmp::get_scale (  )  const [inline]

Definition at line 481 of file cmp.h.

References scale.

00482                 {
00483                         return scale;
00484                 }

float EMAN::OptVarianceCmp::get_shift (  )  const [inline]

Definition at line 486 of file cmp.h.

References shift.

00487                 {
00488                         return shift;
00489                 }

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

Definition at line 464 of file cmp.h.

References OptVarianceCmp().

00465                 {
00466                         return new OptVarianceCmp();
00467                 }


Member Data Documentation

const string OptVarianceCmp::NAME = "optvariance" [static]

Definition at line 491 of file cmp.h.

Referenced by get_name().

float EMAN::OptVarianceCmp::scale [mutable, private]

Definition at line 494 of file cmp.h.

Referenced by cmp(), and get_scale().

float EMAN::OptVarianceCmp::shift [mutable, private]

Definition at line 495 of file cmp.h.

Referenced by cmp(), and get_shift().


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