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
static Cmp *	NEW ()
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 351 of file cmp.h.

---

Constructor & Destructor Documentation

EMAN::OptVarianceCmp::OptVarianceCmp

(

)

[inline]

Definition at line 354 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:

	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 526 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 = 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 = 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_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 358 of file cmp.h.

                {
                        return "optvariance";
                }

string EMAN::OptVarianceCmp::get_desc

(

)

const [inline, virtual]

Implements EMAN::Cmp.

Definition at line 363 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.";
                }

static Cmp* EMAN::OptVarianceCmp::NEW

(

)

[inline, static]

Definition at line 368 of file cmp.h.

References OptVarianceCmp().

                {
                        return new OptVarianceCmp();
                }

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

References scale.

                {
                        return scale;
                }

float EMAN::OptVarianceCmp::get_shift

(

)

const [inline]

Definition at line 390 of file cmp.h.

References shift.

                {
                        return shift;
                }

---

Member Data Documentation

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

Definition at line 396 of file cmp.h.

Referenced by cmp(), and get_scale().

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

Definition at line 397 of file cmp.h.

Referenced by cmp(), and get_shift().

---

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

• cmp.h
• cmp.cpp