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:

[legend]
Collaboration diagram for EMAN::OptVarianceCmp:
[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 358 of file cmp.h.


Constructor & Destructor Documentation

EMAN::OptVarianceCmp::OptVarianceCmp  )  [inline]
 

Definition at line 361 of file cmp.h.

00361 : 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 563 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.

00564 {
00565         ENTERFUNC;
00566         validate_input_args(image, with);
00567 
00568         int keepzero = params.set_default("keepzero", 1);
00569         int invert = params.set_default("invert",0);
00570         int matchfilt = params.set_default("matchfilt",1);
00571         int matchamp = params.set_default("matchamp",0);
00572         int radweight = params.set_default("radweight",0);
00573         int dbug = params.set_default("debug",0);
00574 
00575         size_t size = image->get_xsize() * image->get_ysize() * image->get_zsize();
00576 
00577 
00578         EMData *with2=NULL;
00579         if (matchfilt) {
00580                 EMData *a = image->do_fft();
00581                 EMData *b = with->do_fft();
00582 
00583                 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1);
00584                 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1);
00585 
00586                 float avg=0;
00587                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00588                         rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i]));
00589                         avg+=rfa[i];
00590                 }
00591 
00592                 avg/=a->get_ysize()/2.0f;
00593                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00594                         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
00595                 }
00596                 rfa[0]=0.0;
00597 
00598                 if (dbug) b->write_image("a.hdf",-1);
00599 
00600                 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa);
00601                 with2=b->do_ift();
00602 
00603                 if (dbug) b->write_image("a.hdf",-1);
00604                 if (dbug) a->write_image("a.hdf",-1);
00605 
00606 /*              if (dbug) {
00607                         FILE *out=fopen("a.txt","w");
00608                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]);
00609                         fclose(out);
00610 
00611                         out=fopen("b.txt","w");
00612                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]);
00613                         fclose(out);
00614                 }*/
00615 
00616 
00617                 delete a;
00618                 delete b;
00619 
00620                 if (dbug) {
00621                         with2->write_image("a.hdf",-1);
00622                         image->write_image("a.hdf",-1);
00623                 }
00624 
00625 //              with2->process_inplace("matchfilt",Dict("to",this));
00626 //              x_data = with2->get_data();
00627         }
00628 
00629         // This applies the individual Fourier amplitudes from 'image' and
00630         // applies them to 'with'
00631         if (matchamp) {
00632                 EMData *a = image->do_fft();
00633                 EMData *b = with->do_fft();
00634                 size_t size2 = a->get_xsize() * a->get_ysize() * a->get_zsize();
00635 
00636                 a->ri2ap();
00637                 b->ri2ap();
00638 
00639                 const float *const ad=a->get_const_data();
00640                 float * bd=b->get_data();
00641 
00642                 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i];
00643                 b->update();
00644 
00645                 b->ap2ri();
00646                 with2=b->do_ift();
00647 //with2->write_image("a.hdf",-1);
00648                 delete a;
00649                 delete b;
00650         }
00651 
00652         const float * x_data;
00653         if (with2) x_data=with2->get_const_data();
00654         else x_data = with->get_const_data();
00655         const float *const y_data = image->get_const_data();
00656 
00657         size_t nx = image->get_xsize();
00658         float m = 0;
00659         float b = 0;
00660 
00661         // This will write the x vs y file used to calculate the density
00662         // optimization. This behavior may change in the future
00663         if (dbug) {
00664                 FILE *out=fopen("dbug.optvar.txt","w");
00665                 if (out) {
00666                         for (size_t i=0; i<size; i++) {
00667                                 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]);
00668                         }
00669                         fclose(out);
00670                 }
00671         }
00672 
00673 
00674         Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero);
00675         if (m == 0) {
00676                 m = FLT_MIN;
00677         }
00678         b = -b / m;
00679         m = 1.0f / m;
00680 
00681         // While negative slopes are really not a valid comparison in most cases, we
00682         // still want to detect these instances, so this if is removed
00683 /*      if (m < 0) {
00684                 b = 0;
00685                 m = 1000.0;
00686         }*/
00687 
00688         double  result = 0;
00689         int count = 0;
00690 
00691         if (radweight) {
00692                 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only");
00693                 if (keepzero) {
00694                         for (size_t i = 0,y=0; i < size; y++) {
00695                                 for (size_t x=0; x<nx; i++,x++) {
00696                                         if (y_data[i] && x_data[i]) {
00697 #ifdef  _WIN32
00698                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00699                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00700 #else
00701                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00702                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00703 #endif
00704                                                 count++;
00705                                         }
00706                                 }
00707                         }
00708                         result/=count;
00709                 }
00710                 else {
00711                         for (size_t i = 0,y=0; i < size; y++) {
00712                                 for (size_t x=0; x<nx; i++,x++) {
00713 #ifdef  _WIN32
00714                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00715                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00716 #else
00717                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00718                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00719 #endif
00720                                 }
00721                         }
00722                         result = result / size;
00723                 }
00724         }
00725         else {
00726                 if (keepzero) {
00727                         for (size_t i = 0; i < size; i++) {
00728                                 if (y_data[i] && x_data[i]) {
00729                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00730                                         else result += Util::square((x_data[i] * m) + b - y_data[i]);
00731                                         count++;
00732                                 }
00733                         }
00734                         result/=count;
00735                 }
00736                 else {
00737                         for (size_t i = 0; i < size; i++) {
00738                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00739                                 else result += Util::square((x_data[i] * m) + b - y_data[i]);
00740                         }
00741                         result = result / size;
00742                 }
00743         }
00744         scale = m;
00745         shift = b;
00746 
00747         image->set_attr("ovcmp_m",m);
00748         image->set_attr("ovcmp_b",b);
00749         if (with2) delete with2;
00750         EXITFUNC;
00751 
00752 #if 0
00753         return (1 - result);
00754 #endif
00755 
00756         return static_cast<float>(result);
00757 }

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

Implements EMAN::Cmp.

Definition at line 370 of file cmp.h.

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

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

00366                 {
00367                         return NAME;
00368                 }

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

References EMAN::TypeDict::put().

00381                 {
00382                         TypeDict d;
00383                         d.put("invert", EMObject::INT, "If set, 'with' is rescaled rather than 'this'. 'this' should still be the noisier image. (default=0)");
00384                         d.put("keepzero", EMObject::INT, "If set, zero pixels will not be adjusted in the linear density optimization. (default=1)");
00385                         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)");
00386                         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)");
00387                         d.put("radweight", EMObject::INT, "Upweight variances closer to the edge of the image. (default=0)");
00388                         d.put("debug", EMObject::INT, "Performs various debugging actions if set.");
00389                         return d;
00390                 }

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

Definition at line 392 of file cmp.h.

00393                 {
00394                         return scale;
00395                 }

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

Definition at line 397 of file cmp.h.

00398                 {
00399                         return shift;
00400                 }

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

Definition at line 375 of file cmp.h.

00376                 {
00377                         return new OptVarianceCmp();
00378                 }


Member Data Documentation

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

Definition at line 48 of file cmp.cpp.

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

Definition at line 405 of file cmp.h.

Referenced by cmp().

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

Definition at line 406 of file cmp.h.

Referenced by cmp().


The documentation for this class was generated from the following files:
Generated on Fri Apr 30 15:39:13 2010 for EMAN2 by  doxygen 1.3.9.1