Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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

CmpNEW ()

Static Public Attributes

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


Constructor & Destructor Documentation

EMAN::OptVarianceCmp::OptVarianceCmp  )  [inline]
 

Definition at line 406 of file cmp.h.

00406 : 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 733 of file cmp.cpp.

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

00734 {
00735         ENTERFUNC;
00736         validate_input_args(image, with);
00737 
00738         int keepzero = params.set_default("keepzero", 1);
00739         int invert = params.set_default("invert",0);
00740         int matchfilt = params.set_default("matchfilt",1);
00741         int matchamp = params.set_default("matchamp",0);
00742         int radweight = params.set_default("radweight",0);
00743         int dbug = params.set_default("debug",0);
00744 
00745         size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00746 
00747 
00748         EMData *with2=NULL;
00749         if (matchfilt) {
00750                 EMData *a = image->do_fft();
00751                 EMData *b = with->do_fft();
00752 
00753                 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1);
00754                 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1);
00755 
00756                 float avg=0;
00757                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00758                         rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i]));
00759                         avg+=rfa[i];
00760                 }
00761 
00762                 avg/=a->get_ysize()/2.0f;
00763                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00764                         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
00765                 }
00766                 rfa[0]=0.0;
00767 
00768                 if (dbug) b->write_image("a.hdf",-1);
00769 
00770                 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa);
00771                 with2=b->do_ift();
00772 
00773                 if (dbug) b->write_image("a.hdf",-1);
00774                 if (dbug) a->write_image("a.hdf",-1);
00775 
00776 /*              if (dbug) {
00777                         FILE *out=fopen("a.txt","w");
00778                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]);
00779                         fclose(out);
00780 
00781                         out=fopen("b.txt","w");
00782                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]);
00783                         fclose(out);
00784                 }*/
00785 
00786 
00787                 delete a;
00788                 delete b;
00789 
00790                 if (dbug) {
00791                         with2->write_image("a.hdf",-1);
00792                         image->write_image("a.hdf",-1);
00793                 }
00794 
00795 //              with2->process_inplace("matchfilt",Dict("to",this));
00796 //              x_data = with2->get_data();
00797         }
00798 
00799         // This applies the individual Fourier amplitudes from 'image' and
00800         // applies them to 'with'
00801         if (matchamp) {
00802                 EMData *a = image->do_fft();
00803                 EMData *b = with->do_fft();
00804                 size_t size2 = (size_t)a->get_xsize() * a->get_ysize() * a->get_zsize();
00805 
00806                 a->ri2ap();
00807                 b->ri2ap();
00808 
00809                 const float *const ad=a->get_const_data();
00810                 float * bd=b->get_data();
00811 
00812                 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i];
00813                 b->update();
00814 
00815                 b->ap2ri();
00816                 with2=b->do_ift();
00817 //with2->write_image("a.hdf",-1);
00818                 delete a;
00819                 delete b;
00820         }
00821 
00822         const float * x_data;
00823         if (with2) x_data=with2->get_const_data();
00824         else x_data = with->get_const_data();
00825         const float *const y_data = image->get_const_data();
00826 
00827         size_t nx = image->get_xsize();
00828         float m = 0;
00829         float b = 0;
00830 
00831         // This will write the x vs y file used to calculate the density
00832         // optimization. This behavior may change in the future
00833         if (dbug) {
00834                 FILE *out=fopen("dbug.optvar.txt","w");
00835                 if (out) {
00836                         for (size_t i=0; i<size; i++) {
00837                                 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]);
00838                         }
00839                         fclose(out);
00840                 }
00841         }
00842 
00843 
00844         Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero);
00845         if (m == 0) {
00846                 m = FLT_MIN;
00847         }
00848         b = -b / m;
00849         m = 1.0f / m;
00850 
00851         // While negative slopes are really not a valid comparison in most cases, we
00852         // still want to detect these instances, so this if is removed
00853 /*      if (m < 0) {
00854                 b = 0;
00855                 m = 1000.0;
00856         }*/
00857 
00858         double  result = 0;
00859         int count = 0;
00860 
00861         if (radweight) {
00862                 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only");
00863                 if (keepzero) {
00864                         for (size_t i = 0,y=0; i < size; y++) {
00865                                 for (size_t x=0; x<nx; i++,x++) {
00866                                         if (y_data[i] && x_data[i]) {
00867 #ifdef  _WIN32
00868                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00869                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00870 #else
00871                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00872                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00873 #endif
00874                                                 count++;
00875                                         }
00876                                 }
00877                         }
00878                         result/=count;
00879                 }
00880                 else {
00881                         for (size_t i = 0,y=0; i < size; y++) {
00882                                 for (size_t x=0; x<nx; i++,x++) {
00883 #ifdef  _WIN32
00884                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00885                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00886 #else
00887                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00888                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00889 #endif
00890                                 }
00891                         }
00892                         result = result / size;
00893                 }
00894         }
00895         else {
00896                 if (keepzero) {
00897                         for (size_t i = 0; i < size; i++) {
00898                                 if (y_data[i] && x_data[i]) {
00899                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00900                                         else result += Util::square((x_data[i] * m) + b - y_data[i]);
00901                                         count++;
00902                                 }
00903                         }
00904                         result/=count;
00905                 }
00906                 else {
00907                         for (size_t i = 0; i < size; i++) {
00908                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00909                                 else result += Util::square((x_data[i] * m) + b - y_data[i]);
00910                         }
00911                         result = result / size;
00912                 }
00913         }
00914         scale = m;
00915         shift = b;
00916 
00917         image->set_attr("ovcmp_m",m);
00918         image->set_attr("ovcmp_b",b);
00919         if (with2) delete with2;
00920         EXITFUNC;
00921 
00922 #if 0
00923         return (1 - result);
00924 #endif
00925 
00926         return static_cast<float>(result);
00927 }

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

Implements EMAN::Cmp.

Definition at line 415 of file cmp.h.

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

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

00411                 {
00412                         return NAME;
00413                 }

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

References EMAN::TypeDict::put().

00426                 {
00427                         TypeDict d;
00428                         d.put("invert", EMObject::INT, "If set, 'with' is rescaled rather than 'this'. 'this' should still be the noisier image. (default=0)");
00429                         d.put("keepzero", EMObject::INT, "If set, zero pixels will not be adjusted in the linear density optimization. (default=1)");
00430                         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)");
00431                         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)");
00432                         d.put("radweight", EMObject::INT, "Upweight variances closer to the edge of the image. (default=0)");
00433                         d.put("debug", EMObject::INT, "Performs various debugging actions if set.");
00434                         return d;
00435                 }

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

Definition at line 437 of file cmp.h.

00438                 {
00439                         return scale;
00440                 }

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

Definition at line 442 of file cmp.h.

00443                 {
00444                         return shift;
00445                 }

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

Definition at line 420 of file cmp.h.

00421                 {
00422                         return new OptVarianceCmp();
00423                 }


Member Data Documentation

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

Definition at line 54 of file cmp.cpp.

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

Definition at line 450 of file cmp.h.

Referenced by cmp().

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

Definition at line 451 of file cmp.h.

Referenced by cmp().


The documentation for this class was generated from the following files:
Generated on Tue Jul 12 13:50:30 2011 for EMAN2 by  doxygen 1.3.9.1