EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes
EMAN::DotCmp Class Reference

Use dot product of 2 same-size images to do the comparison. More...

#include <cmp.h>

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

List of all members.

Public Member Functions

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.

Static Public Member Functions

static CmpNEW ()

Static Public Attributes

static const string NAME = "dot"

Detailed Description

Use dot product of 2 same-size images to do the comparison.

// Added mask option PAP 04/23/06 For complex images, it does not check r/i vs a/p.

Author:
Steve Ludtke (sludtke@bcm.tmc.edu)
Date:
2005-07-13
Parameters:
negativeReturns -1 * dot product, default true
normalizeReturns normalized dot product -1.0 - 1.0

Definition at line 272 of file cmp.h.


Member Function Documentation

float DotCmp::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 410 of file cmp.cpp.

References dm, dot_cmp_cuda(), ENTERFUNC, EXITFUNC, EMAN::EMData::get_attr(), EMAN::EMData::get_const_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Dict::has_key(), EMAN::EMData::is_complex(), EMAN::EMData::is_fftodd(), nx, ny, EMAN::Cmp::params, EMAN::Dict::set_default(), sqrt(), and EMAN::Cmp::validate_input_args().

Referenced by EMAN::EMData::dot().

{
        ENTERFUNC;
        
        validate_input_args(image, with);

        int normalize = params.set_default("normalize", 0);
        float negative = (float)params.set_default("negative", 1);
        if (negative) negative=-1.0; else negative=1.0;
#ifdef EMAN2_USING_CUDA // SO far only works for real images I put CUDA first to avoid running non CUDA overhead (calls to getdata are expensive!!!!)
        if(image->is_complex() && with->is_complex()) {
        } else {
                if (image->getcudarwdata() && with->getcudarwdata()) {
                        //cout << "CUDA dot cmp" << endl;
                        float* maskdata = 0;
                        bool has_mask = false;
                        EMData* mask = 0;
                        if (params.has_key("mask")) {
                                mask = params["mask"];
                                if(mask!=0) {has_mask=true;}
                        }
                        if(has_mask && !mask->getcudarwdata()){
                                mask->copy_to_cuda();
                                maskdata = mask->getcudarwdata();
                        }

                        float result = dot_cmp_cuda(image->getcudarwdata(), with->getcudarwdata(), maskdata, image->get_xsize(), image->get_ysize(), image->get_zsize());
                        result *= negative;

                        return result;
                        
                }
        }
#endif
        const float *const x_data = image->get_const_data();
        const float *const y_data = with->get_const_data();

        double result = 0.;
        long n = 0;
        if(image->is_complex() && with->is_complex()) {
        // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
                int nx  = with->get_xsize();
                int ny  = with->get_ysize();
                int nz  = with->get_zsize();
                nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
                int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image

                int ixb = 2*((nx+1)%2);
                int iyb = ny%2;
                //
                if(nz == 1) {
                //  it looks like it could work in 3D, but does not
                for ( int iz = 0; iz <= nz-1; ++iz) {
                        double part = 0.;
                        for ( int iy = 0; iy <= ny-1; ++iy) {
                                for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
                                        size_t ii = ix + (iy  + iz * ny)* lsd2;
                                        part += x_data[ii] * double(y_data[ii]);
                                }
                        }
                        for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
                                size_t ii = (iy  + iz * ny)* lsd2;
                                part += x_data[ii] * double(y_data[ii]);
                                part += x_data[ii+1] * double(y_data[ii+1]);
                        }
                        if(nx%2 == 0) {
                                for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
                                        size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
                                        part += x_data[ii] * double(y_data[ii]);
                                        part += x_data[ii+1] * double(y_data[ii+1]);
                                }

                        }
                        part *= 2;
                        part += x_data[0] * double(y_data[0]);
                        if(ny%2 == 0) {
                                size_t ii = (ny/2  + iz * ny)* lsd2;
                                part += x_data[ii] * double(y_data[ii]);
                        }
                        if(nx%2 == 0) {
                                size_t ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
                                part += x_data[ii] * double(y_data[ii]);
                                if(ny%2 == 0) {
                                        int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
                                        part += x_data[ii] * double(y_data[ii]);
                                }
                        }
                        result += part;
                }
                if( normalize ) {
                //  it looks like it could work in 3D, but does not
                double square_sum1 = 0., square_sum2 = 0.;
                for ( int iz = 0; iz <= nz-1; ++iz) {
                        for ( int iy = 0; iy <= ny-1; ++iy) {
                                for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
                                        size_t ii = ix + (iy  + iz * ny)* lsd2;
                                        square_sum1 += x_data[ii] * double(x_data[ii]);
                                        square_sum2 += y_data[ii] * double(y_data[ii]);
                                }
                        }
                        for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
                                size_t ii = (iy  + iz * ny)* lsd2;
                                square_sum1 += x_data[ii] * double(x_data[ii]);
                                square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
                                square_sum2 += y_data[ii] * double(y_data[ii]);
                                square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
                        }
                        if(nx%2 == 0) {
                                for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
                                        size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
                                        square_sum1 += x_data[ii] * double(x_data[ii]);
                                        square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
                                        square_sum2 += y_data[ii] * double(y_data[ii]);
                                        square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
                                }

                        }
                        square_sum1 *= 2;
                        square_sum1 += x_data[0] * double(x_data[0]);
                        square_sum2 *= 2;
                        square_sum2 += y_data[0] * double(y_data[0]);
                        if(ny%2 == 0) {
                                int ii = (ny/2  + iz * ny)* lsd2;
                                square_sum1 += x_data[ii] * double(x_data[ii]);
                                square_sum2 += y_data[ii] * double(y_data[ii]);
                        }
                        if(nx%2 == 0) {
                                int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
                                square_sum1 += x_data[ii] * double(x_data[ii]);
                                square_sum2 += y_data[ii] * double(y_data[ii]);
                                if(ny%2 == 0) {
                                        int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
                                        square_sum1 += x_data[ii] * double(x_data[ii]);
                                        square_sum2 += y_data[ii] * double(y_data[ii]);
                                }
                        }
                }
                result /= sqrt(square_sum1*square_sum2);
                } 
                        else  result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
                        // This concludes the 2D part.
                } else {
                // The 3D code
                int ky, kz;
                int ny2 = ny/2; int nz2 = nz/2;
                for ( int iz = 0; iz <= nz-1; ++iz) {
                        if(iz>nz2) kz=iz-nz; else kz=iz;
                        for ( int iy = 0; iy <= ny-1; ++iy) {
                                if(iy>ny2) ky=iy-ny; else ky=iy;
                                for (int ix = 0; ix <= lsd2-1; ++ix) {
                                                int p = 2;
                                                if  (ix <= 1 || ix >= lsd2-2 && nx%2 == 0)   p = 1;
                                                size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
                                                result += p * x_data[ii] * double(y_data[ii]);
                                }
                        }
                }
                if( normalize ) {
                double square_sum1 = 0., square_sum2 = 0.;
                int ky, kz;
                int ny2 = ny/2; int nz2 = nz/2;
                for ( int iz = 0; iz <= nz-1; ++iz) {
                        if(iz>nz2) kz=iz-nz; else kz=iz;
                        for ( int iy = 0; iy <= ny-1; ++iy) {
                                if(iy>ny2) ky=iy-ny; else ky=iy;
                                for (int ix = 0; ix <= lsd2-1; ++ix) {
                                                int p = 2;
                                                if  (ix <= 1 || ix >= lsd2-2 && nx%2 == 0)   p = 1;
                                                size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
                                                square_sum1 += p * x_data[ii] * double(x_data[ii]);
                                                square_sum2 += p * y_data[ii] * double(y_data[ii]);
                                }
                        }
                }
                result /= sqrt(square_sum1*square_sum2);
                } 
                        else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
                }
        } else {
                // This part is for when two images are real, which is much easier because 2-D or 3-D
                // is the same code.
                size_t totsize = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();

                double square_sum1 = 0., square_sum2 = 0.;

                if (params.has_key("mask")) {
                        EMData* mask;
                        mask = params["mask"];
                        const float *const dm = mask->get_const_data();
                        if (normalize) {
                                for (size_t i = 0; i < totsize; i++) {
                                        if (dm[i] > 0.5) {
                                                square_sum1 += x_data[i]*double(x_data[i]);
                                                square_sum2 += y_data[i]*double(y_data[i]);
                                                result += x_data[i]*double(y_data[i]);
                                        }
                                }
                        } else {
                                for (size_t i = 0; i < totsize; i++) {
                                        if (dm[i] > 0.5) {
                                                result += x_data[i]*double(y_data[i]);
                                                n++;
                                        }
                                }
                        }
                } else {
                        // this little bit of manual loop unrolling makes the dot product as fast as sqeuclidean with -O2
                        for (size_t i=0; i<totsize; i++) result+=x_data[i]*y_data[i];

                        if (normalize) {
                                square_sum1 = image->get_attr("square_sum");
                                square_sum2 = with->get_attr("square_sum");
                        } else n = totsize;
                }
                if (normalize) result /= (sqrt(square_sum1*square_sum2)); else result /= n;
        }


        EXITFUNC;
        return (float) (negative*result);
}
string EMAN::DotCmp::get_desc ( ) const [inline, virtual]

Implements EMAN::Cmp.

Definition at line 282 of file cmp.h.

                {
                        return "Dot product (default -1 * dot product)";
                }
string EMAN::DotCmp::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 277 of file cmp.h.

References NAME.

                {
                        return NAME;
                }
TypeDict EMAN::DotCmp::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 292 of file cmp.h.

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

                {
                        TypeDict d;
                        d.put("negative", EMObject::INT, "If set, returns -1 * dot product. Set by default so smaller is better");
                        d.put("normalize", EMObject::INT, "If set, returns normalized dot product (cosine of the angle) -1.0 - 1.0.");
                        d.put("mask", EMObject::EMDATA, "image mask");
                        return d;
                }
static Cmp* EMAN::DotCmp::NEW ( ) [inline, static]

Definition at line 287 of file cmp.h.

                {
                        return new DotCmp();
                }

Member Data Documentation

const string DotCmp::NAME = "dot" [static]

Definition at line 301 of file cmp.h.

Referenced by get_name().


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