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 ()


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 231 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 248 of file cmp.cpp.

References dm, 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().

00249 {
00250         ENTERFUNC;
00251         validate_input_args(image, with);
00252 
00253         const float *const x_data = image->get_const_data();
00254         const float *const y_data = with->get_const_data();
00255 
00256         int normalize = params.set_default("normalize", 0);
00257         float negative = (float)params.set_default("negative", 1);
00258 
00259         if (negative) negative=-1.0; else negative=1.0;
00260         double result = 0.;
00261         long n = 0;
00262         if(image->is_complex() && with->is_complex()) {
00263         // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
00264                 int nx  = with->get_xsize();
00265                 int ny  = with->get_ysize();
00266                 int nz  = with->get_zsize();
00267                 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
00268                 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
00269 
00270                 int ixb = 2*((nx+1)%2);
00271                 int iyb = ny%2;
00272                 //
00273                 if(nz == 1) {
00274                 //  it looks like it could work in 3D, but does not
00275                 for ( int iz = 0; iz <= nz-1; ++iz) {
00276                         double part = 0.;
00277                         for ( int iy = 0; iy <= ny-1; ++iy) {
00278                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00279                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00280                                         part += x_data[ii] * double(y_data[ii]);
00281                                 }
00282                         }
00283                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00284                                 size_t ii = (iy  + iz * ny)* lsd2;
00285                                 part += x_data[ii] * double(y_data[ii]);
00286                                 part += x_data[ii+1] * double(y_data[ii+1]);
00287                         }
00288                         if(nx%2 == 0) {
00289                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00290                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00291                                         part += x_data[ii] * double(y_data[ii]);
00292                                         part += x_data[ii+1] * double(y_data[ii+1]);
00293                                 }
00294 
00295                         }
00296                         part *= 2;
00297                         part += x_data[0] * double(y_data[0]);
00298                         if(ny%2 == 0) {
00299                                 size_t ii = (ny/2  + iz * ny)* lsd2;
00300                                 part += x_data[ii] * double(y_data[ii]);
00301                         }
00302                         if(nx%2 == 0) {
00303                                 size_t ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00304                                 part += x_data[ii] * double(y_data[ii]);
00305                                 if(ny%2 == 0) {
00306                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00307                                         part += x_data[ii] * double(y_data[ii]);
00308                                 }
00309                         }
00310                         result += part;
00311                 }
00312                 if( normalize ) {
00313                 //  it looks like it could work in 3D, but does not
00314                 double square_sum1 = 0., square_sum2 = 0.;
00315                 for ( int iz = 0; iz <= nz-1; ++iz) {
00316                         for ( int iy = 0; iy <= ny-1; ++iy) {
00317                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00318                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00319                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00320                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00321                                 }
00322                         }
00323                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00324                                 size_t ii = (iy  + iz * ny)* lsd2;
00325                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00326                                 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00327                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00328                                 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00329                         }
00330                         if(nx%2 == 0) {
00331                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00332                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00333                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00334                                         square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00335                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00336                                         square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00337                                 }
00338 
00339                         }
00340                         square_sum1 *= 2;
00341                         square_sum1 += x_data[0] * double(x_data[0]);
00342                         square_sum2 *= 2;
00343                         square_sum2 += y_data[0] * double(y_data[0]);
00344                         if(ny%2 == 0) {
00345                                 int ii = (ny/2  + iz * ny)* lsd2;
00346                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00347                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00348                         }
00349                         if(nx%2 == 0) {
00350                                 int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00351                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00352                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00353                                 if(ny%2 == 0) {
00354                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00355                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00356                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00357                                 }
00358                         }
00359                 }
00360                 result /= sqrt(square_sum1*square_sum2);
00361                 } else  result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
00362 
00363                 } else { //This 3D code is incorrect, but it is the best I can do now 01/09/06 PAP
00364                 int ky, kz;
00365                 int ny2 = ny/2; int nz2 = nz/2;
00366                 for ( int iz = 0; iz <= nz-1; ++iz) {
00367                         if(iz>nz2) kz=iz-nz; else kz=iz;
00368                         for ( int iy = 0; iy <= ny-1; ++iy) {
00369                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00370                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00371                                         // Skip Friedel related values
00372                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00373                                                 size_t ii = ix + (iy  + iz * ny)* lsd2;
00374                                                 result += x_data[ii] * double(y_data[ii]);
00375                                         }
00376                                 }
00377                         }
00378                 }
00379                 if( normalize ) {
00380                 //  still incorrect
00381                 double square_sum1 = 0., square_sum2 = 0.;
00382                 int ky, kz;
00383                 int ny2 = ny/2; int nz2 = nz/2;
00384                 for ( int iz = 0; iz <= nz-1; ++iz) {
00385                         if(iz>nz2) kz=iz-nz; else kz=iz;
00386                         for ( int iy = 0; iy <= ny-1; ++iy) {
00387                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00388                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00389                                         // Skip Friedel related values
00390                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00391                                                 size_t ii = ix + (iy  + 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                 }
00398                 result /= sqrt(square_sum1*square_sum2);
00399                 } else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz/2);
00400                 }
00401         } else {
00402                 size_t totsize = image->get_xsize() * image->get_ysize() * image->get_zsize();
00403 
00404                 double square_sum1 = 0., square_sum2 = 0.;
00405 
00406                 if (params.has_key("mask")) {
00407                         EMData* mask;
00408                         mask = params["mask"];
00409                         const float *const dm = mask->get_const_data();
00410                         if (normalize) {
00411                                 for (size_t i = 0; i < totsize; i++) {
00412                                         if (dm[i] > 0.5) {
00413                                                 square_sum1 += x_data[i]*double(x_data[i]);
00414                                                 square_sum2 += y_data[i]*double(y_data[i]);
00415                                                 result += x_data[i]*double(y_data[i]);
00416                                         }
00417                                 }
00418                         } else {
00419                                 for (size_t i = 0; i < totsize; i++) {
00420                                         if (dm[i] > 0.5) {
00421                                                 result += x_data[i]*double(y_data[i]);
00422                                                 n++;
00423                                         }
00424                                 }
00425                         }
00426                 } else {
00427                         // this little bit of manual loop unrolling makes the dot product as fast as sqeuclidean with -O2
00428                         for (size_t i=0; i<totsize; i++) result+=x_data[i]*y_data[i];
00429 
00430                         if (normalize) {
00431                                 square_sum1 = image->get_attr("square_sum");
00432                                 square_sum2 = with->get_attr("square_sum");
00433                         } else n = totsize;
00434                 }
00435                 if (normalize) result /= (sqrt(square_sum1*square_sum2)); else result /= n;
00436         }
00437 
00438 
00439         EXITFUNC;
00440         return (float) (negative*result);
00441 }

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

00237                 {
00238                         return "dot";
00239                 }

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

Implements EMAN::Cmp.

Definition at line 241 of file cmp.h.

00242                 {
00243                         return "Dot product (default -1 * dot product)";
00244                 }

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

Definition at line 246 of file cmp.h.

00247                 {
00248                         return new DotCmp();
00249                 }

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

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

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


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

Generated on Sat Nov 21 02:20:17 2009 for EMAN2 by  doxygen 1.5.6