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

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:

[legend]
Collaboration diagram for EMAN::DotCmp:
[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

CmpNEW ()

Static Public Attributes

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:
negative Returns -1 * dot product, default true
normalize Returns normalized dot product -1.0 - 1.0

Definition at line 234 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:
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 285 of file cmp.cpp.

References dm, 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::Dict::set_default(), sqrt(), and EMAN::Cmp::validate_input_args().

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

00286 {
00287         ENTERFUNC;
00288         validate_input_args(image, with);
00289 
00290         const float *const x_data = image->get_const_data();
00291         const float *const y_data = with->get_const_data();
00292 
00293         int normalize = params.set_default("normalize", 0);
00294         float negative = (float)params.set_default("negative", 1);
00295 
00296         if (negative) negative=-1.0; else negative=1.0;
00297         double result = 0.;
00298         long n = 0;
00299         if(image->is_complex() && with->is_complex()) {
00300         // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
00301                 int nx  = with->get_xsize();
00302                 int ny  = with->get_ysize();
00303                 int nz  = with->get_zsize();
00304                 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
00305                 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
00306 
00307                 int ixb = 2*((nx+1)%2);
00308                 int iyb = ny%2;
00309                 //
00310                 if(nz == 1) {
00311                 //  it looks like it could work in 3D, but does not
00312                 for ( int iz = 0; iz <= nz-1; ++iz) {
00313                         double part = 0.;
00314                         for ( int iy = 0; iy <= ny-1; ++iy) {
00315                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00316                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00317                                         part += x_data[ii] * double(y_data[ii]);
00318                                 }
00319                         }
00320                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00321                                 size_t ii = (iy  + iz * ny)* lsd2;
00322                                 part += x_data[ii] * double(y_data[ii]);
00323                                 part += x_data[ii+1] * double(y_data[ii+1]);
00324                         }
00325                         if(nx%2 == 0) {
00326                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00327                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00328                                         part += x_data[ii] * double(y_data[ii]);
00329                                         part += x_data[ii+1] * double(y_data[ii+1]);
00330                                 }
00331 
00332                         }
00333                         part *= 2;
00334                         part += x_data[0] * double(y_data[0]);
00335                         if(ny%2 == 0) {
00336                                 size_t ii = (ny/2  + iz * ny)* lsd2;
00337                                 part += x_data[ii] * double(y_data[ii]);
00338                         }
00339                         if(nx%2 == 0) {
00340                                 size_t ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00341                                 part += x_data[ii] * double(y_data[ii]);
00342                                 if(ny%2 == 0) {
00343                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00344                                         part += x_data[ii] * double(y_data[ii]);
00345                                 }
00346                         }
00347                         result += part;
00348                 }
00349                 if( normalize ) {
00350                 //  it looks like it could work in 3D, but does not
00351                 double square_sum1 = 0., square_sum2 = 0.;
00352                 for ( int iz = 0; iz <= nz-1; ++iz) {
00353                         for ( int iy = 0; iy <= ny-1; ++iy) {
00354                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00355                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00356                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00357                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00358                                 }
00359                         }
00360                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00361                                 size_t ii = (iy  + iz * ny)* lsd2;
00362                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00363                                 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00364                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00365                                 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00366                         }
00367                         if(nx%2 == 0) {
00368                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00369                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00370                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00371                                         square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00372                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00373                                         square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00374                                 }
00375 
00376                         }
00377                         square_sum1 *= 2;
00378                         square_sum1 += x_data[0] * double(x_data[0]);
00379                         square_sum2 *= 2;
00380                         square_sum2 += y_data[0] * double(y_data[0]);
00381                         if(ny%2 == 0) {
00382                                 int ii = (ny/2  + iz * ny)* lsd2;
00383                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00384                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00385                         }
00386                         if(nx%2 == 0) {
00387                                 int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00388                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00389                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00390                                 if(ny%2 == 0) {
00391                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00392                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00393                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00394                                 }
00395                         }
00396                 }
00397                 result /= sqrt(square_sum1*square_sum2);
00398                 } else  result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
00399 
00400                 } else { //This 3D code is incorrect, but it is the best I can do now 01/09/06 PAP
00401                 int ky, kz;
00402                 int ny2 = ny/2; int nz2 = nz/2;
00403                 for ( int iz = 0; iz <= nz-1; ++iz) {
00404                         if(iz>nz2) kz=iz-nz; else kz=iz;
00405                         for ( int iy = 0; iy <= ny-1; ++iy) {
00406                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00407                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00408                                         // Skip Friedel related values
00409                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00410                                                 size_t ii = ix + (iy  + iz * ny)* lsd2;
00411                                                 result += x_data[ii] * double(y_data[ii]);
00412                                         }
00413                                 }
00414                         }
00415                 }
00416                 if( normalize ) {
00417                 //  still incorrect
00418                 double square_sum1 = 0., square_sum2 = 0.;
00419                 int ky, kz;
00420                 int ny2 = ny/2; int nz2 = nz/2;
00421                 for ( int iz = 0; iz <= nz-1; ++iz) {
00422                         if(iz>nz2) kz=iz-nz; else kz=iz;
00423                         for ( int iy = 0; iy <= ny-1; ++iy) {
00424                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00425                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00426                                         // Skip Friedel related values
00427                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00428                                                 size_t ii = ix + (iy  + iz * ny)* lsd2;
00429                                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00430                                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00431                                         }
00432                                 }
00433                         }
00434                 }
00435                 result /= sqrt(square_sum1*square_sum2);
00436                 } else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz/2);
00437                 }
00438         } else {
00439                 size_t totsize = image->get_xsize() * image->get_ysize() * image->get_zsize();
00440 
00441                 double square_sum1 = 0., square_sum2 = 0.;
00442 
00443                 if (params.has_key("mask")) {
00444                         EMData* mask;
00445                         mask = params["mask"];
00446                         const float *const dm = mask->get_const_data();
00447                         if (normalize) {
00448                                 for (size_t i = 0; i < totsize; i++) {
00449                                         if (dm[i] > 0.5) {
00450                                                 square_sum1 += x_data[i]*double(x_data[i]);
00451                                                 square_sum2 += y_data[i]*double(y_data[i]);
00452                                                 result += x_data[i]*double(y_data[i]);
00453                                         }
00454                                 }
00455                         } else {
00456                                 for (size_t i = 0; i < totsize; i++) {
00457                                         if (dm[i] > 0.5) {
00458                                                 result += x_data[i]*double(y_data[i]);
00459                                                 n++;
00460                                         }
00461                                 }
00462                         }
00463                 } else {
00464                         // this little bit of manual loop unrolling makes the dot product as fast as sqeuclidean with -O2
00465                         for (size_t i=0; i<totsize; i++) result+=x_data[i]*y_data[i];
00466 
00467                         if (normalize) {
00468                                 square_sum1 = image->get_attr("square_sum");
00469                                 square_sum2 = with->get_attr("square_sum");
00470                         } else n = totsize;
00471                 }
00472                 if (normalize) result /= (sqrt(square_sum1*square_sum2)); else result /= n;
00473         }
00474 
00475 
00476         EXITFUNC;
00477         return (float) (negative*result);
00478 }

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

Implements EMAN::Cmp.

Definition at line 244 of file cmp.h.

00245                 {
00246                         return "Dot product (default -1 * dot product)";
00247                 }

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

00240                 {
00241                         return NAME;
00242                 }

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

References EMAN::TypeDict::put().

00255                 {
00256                         TypeDict d;
00257                         d.put("negative", EMObject::INT, "If set, returns -1 * dot product. Set by default so smaller is better");
00258                         d.put("normalize", EMObject::INT, "If set, returns normalized dot product (cosine of the angle) -1.0 - 1.0.");
00259                         d.put("mask", EMObject::EMDATA, "image mask");
00260                         return d;
00261                 }

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

Definition at line 249 of file cmp.h.

00250                 {
00251                         return new DotCmp();
00252                 }


Member Data Documentation

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

Definition at line 45 of file cmp.cpp.


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