EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Attributes
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 490 of file cmp.h.


Constructor & Destructor Documentation

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

Definition at line 493 of file cmp.h.

Referenced by NEW().

: 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:
imageThe first image to be compared.
withThe second image to be comppared.
Returns:
The comparison result. Smaller better by default

Implements EMAN::Cmp.

Definition at line 862 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(), ENTERFUNC, EXITFUNC, EMAN::EMData::get_const_data(), EMAN::EMData::get_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::EMData::update(), EMAN::Cmp::validate_input_args(), EMAN::EMData::write_image(), x, and y.

{
        ENTERFUNC;
        validate_input_args(image, with);

        int keepzero = params.set_default("keepzero", 1);
        int invert = params.set_default("invert",0);
        int matchfilt = params.set_default("matchfilt",1);
        int matchamp = params.set_default("matchamp",0);
        int radweight = params.set_default("radweight",0);
        int dbug = params.set_default("debug",0);

        size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();


        EMData *with2=NULL;
        if (matchfilt) {
                EMData *a = image->do_fft();
                EMData *b = with->do_fft();

                vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1);
                vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1);

                float avg=0;
                for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
                        rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i]));
                        avg+=rfa[i];
                }

                avg/=a->get_ysize()/2.0f;
                for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
                        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
                }
                rfa[0]=0.0;

                if (dbug) b->write_image("a.hdf",-1);

                b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa);
                with2=b->do_ift();

                if (dbug) b->write_image("a.hdf",-1);
                if (dbug) a->write_image("a.hdf",-1);

/*              if (dbug) {
                        FILE *out=fopen("a.txt","w");
                        for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]);
                        fclose(out);

                        out=fopen("b.txt","w");
                        for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]);
                        fclose(out);
                }*/


                delete a;
                delete b;

                if (dbug) {
                        with2->write_image("a.hdf",-1);
                        image->write_image("a.hdf",-1);
                }

//              with2->process_inplace("matchfilt",Dict("to",this));
//              x_data = with2->get_data();
        }

        // This applies the individual Fourier amplitudes from 'image' and
        // applies them to 'with'
        if (matchamp) {
                EMData *a = image->do_fft();
                EMData *b = with->do_fft();
                size_t size2 = (size_t)a->get_xsize() * a->get_ysize() * a->get_zsize();

                a->ri2ap();
                b->ri2ap();

                const float *const ad=a->get_const_data();
                float * bd=b->get_data();

                for (size_t i=0; i<size2; i+=2) bd[i]=ad[i];
                b->update();

                b->ap2ri();
                with2=b->do_ift();
//with2->write_image("a.hdf",-1);
                delete a;
                delete b;
        }

        const float * x_data;
        if (with2) x_data=with2->get_const_data();
        else x_data = with->get_const_data();
        const float *const y_data = image->get_const_data();

        size_t nx = image->get_xsize();
        float m = 0;
        float b = 0;

        // This will write the x vs y file used to calculate the density
        // optimization. This behavior may change in the future
        if (dbug) {
                FILE *out=fopen("dbug.optvar.txt","w");
                if (out) {
                        for (size_t i=0; i<size; i++) {
                                if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]);
                        }
                        fclose(out);
                }
        }


        Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero);
        if (m == 0) {
                m = FLT_MIN;
        }
        b = -b / m;
        m = 1.0f / m;

        // While negative slopes are really not a valid comparison in most cases, we
        // still want to detect these instances, so this if is removed
/*      if (m < 0) {
                b = 0;
                m = 1000.0;
        }*/

        double  result = 0;
        int count = 0;

        if (radweight) {
                if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only");
                if (keepzero) {
                        for (size_t i = 0,y=0; i < size; y++) {
                                for (size_t x=0; x<nx; i++,x++) {
                                        if (y_data[i] && x_data[i]) {
#ifdef  _WIN32
                                                if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
                                                else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
#else
                                                if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
                                                else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
#endif
                                                count++;
                                        }
                                }
                        }
                        result/=count;
                }
                else {
                        for (size_t i = 0,y=0; i < size; y++) {
                                for (size_t x=0; x<nx; i++,x++) {
#ifdef  _WIN32
                                        if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
                                        else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
#else
                                        if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
                                        else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
#endif
                                }
                        }
                        result = result / size;
                }
        }
        else {
                if (keepzero) {
                        for (size_t i = 0; i < size; i++) {
                                if (y_data[i] && x_data[i]) {
                                        if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
                                        else result += Util::square((x_data[i] * m) + b - y_data[i]);
                                        count++;
                                }
                        }
                        result/=count;
                }
                else {
                        for (size_t i = 0; i < size; i++) {
                                if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
                                else result += Util::square((x_data[i] * m) + b - y_data[i]);
                        }
                        result = result / size;
                }
        }
        scale = m;
        shift = b;

        image->set_attr("ovcmp_m",m);
        image->set_attr("ovcmp_b",b);
        if (with2) delete with2;
        EXITFUNC;

#if 0
        return (1 - result);
#endif

        return static_cast<float>(result);
}
string EMAN::OptVarianceCmp::get_desc ( ) const [inline, virtual]

Implements EMAN::Cmp.

Definition at line 502 of file cmp.h.

                {
                        return "Real-space variance after density optimization, self should be noisy and target less noisy. Linear transform applied to density to minimize variance.";
                }
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 497 of file cmp.h.

References NAME.

                {
                        return NAME;
                }
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 512 of file cmp.h.

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

                {
                        TypeDict d;
                        d.put("invert", EMObject::INT, "If set, 'with' is rescaled rather than 'this'. 'this' should still be the noisier image. (default=0)");
                        d.put("keepzero", EMObject::INT, "If set, zero pixels will not be adjusted in the linear density optimization. (default=1)");
                        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)");
                        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)");
                        d.put("radweight", EMObject::INT, "Upweight variances closer to the edge of the image. (default=0)");
                        d.put("debug", EMObject::INT, "Performs various debugging actions if set.");
                        return d;
                }
float EMAN::OptVarianceCmp::get_scale ( ) const [inline]

Definition at line 524 of file cmp.h.

References scale.

                {
                        return scale;
                }
float EMAN::OptVarianceCmp::get_shift ( ) const [inline]

Definition at line 529 of file cmp.h.

References shift.

                {
                        return shift;
                }
static Cmp* EMAN::OptVarianceCmp::NEW ( ) [inline, static]

Definition at line 507 of file cmp.h.

References OptVarianceCmp().

                {
                        return new OptVarianceCmp();
                }

Member Data Documentation

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

Definition at line 534 of file cmp.h.

Referenced by get_name().

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

Definition at line 537 of file cmp.h.

Referenced by cmp(), and get_scale().

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

Definition at line 538 of file cmp.h.

Referenced by cmp(), and get_shift().


The documentation for this class was generated from the following files: