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

EMAN::DotCmp Class Reference
[a function or class that is CUDA enabled]

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

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 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:
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 406 of file cmp.cpp.

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

00407 {
00408         ENTERFUNC;
00409         
00410         validate_input_args(image, with);
00411 
00412         int normalize = params.set_default("normalize", 0);
00413         float negative = (float)params.set_default("negative", 1);
00414         if (negative) negative=-1.0; else negative=1.0;
00415 #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!!!!)
00416         if(image->is_complex() && with->is_complex()) {
00417         } else {
00418                 if (image->getcudarwdata() && with->getcudarwdata()) {
00419                         //cout << "CUDA dot cmp" << endl;
00420                         float* maskdata = 0;
00421                         bool has_mask = false;
00422                         EMData* mask = 0;
00423                         if (params.has_key("mask")) {
00424                                 mask = params["mask"];
00425                                 if(mask!=0) {has_mask=true;}
00426                         }
00427                         if(has_mask && !mask->getcudarwdata()){
00428                                 mask->copy_to_cuda();
00429                                 maskdata = mask->getcudarwdata();
00430                         }
00431 
00432                         float result = dot_cmp_cuda(image->getcudarwdata(), with->getcudarwdata(), maskdata, image->get_xsize(), image->get_ysize(), image->get_zsize());
00433                         result *= negative;
00434 
00435                         return result;
00436                         
00437                 }
00438         }
00439 #endif
00440         const float *const x_data = image->get_const_data();
00441         const float *const y_data = with->get_const_data();
00442 
00443         double result = 0.;
00444         long n = 0;
00445         if(image->is_complex() && with->is_complex()) {
00446         // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
00447                 int nx  = with->get_xsize();
00448                 int ny  = with->get_ysize();
00449                 int nz  = with->get_zsize();
00450                 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
00451                 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
00452 
00453                 int ixb = 2*((nx+1)%2);
00454                 int iyb = ny%2;
00455                 //
00456                 if(nz == 1) {
00457                 //  it looks like it could work in 3D, but does not
00458                 for ( int iz = 0; iz <= nz-1; ++iz) {
00459                         double part = 0.;
00460                         for ( int iy = 0; iy <= ny-1; ++iy) {
00461                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00462                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00463                                         part += x_data[ii] * double(y_data[ii]);
00464                                 }
00465                         }
00466                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00467                                 size_t ii = (iy  + iz * ny)* lsd2;
00468                                 part += x_data[ii] * double(y_data[ii]);
00469                                 part += x_data[ii+1] * double(y_data[ii+1]);
00470                         }
00471                         if(nx%2 == 0) {
00472                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00473                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00474                                         part += x_data[ii] * double(y_data[ii]);
00475                                         part += x_data[ii+1] * double(y_data[ii+1]);
00476                                 }
00477 
00478                         }
00479                         part *= 2;
00480                         part += x_data[0] * double(y_data[0]);
00481                         if(ny%2 == 0) {
00482                                 size_t ii = (ny/2  + iz * ny)* lsd2;
00483                                 part += x_data[ii] * double(y_data[ii]);
00484                         }
00485                         if(nx%2 == 0) {
00486                                 size_t ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00487                                 part += x_data[ii] * double(y_data[ii]);
00488                                 if(ny%2 == 0) {
00489                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00490                                         part += x_data[ii] * double(y_data[ii]);
00491                                 }
00492                         }
00493                         result += part;
00494                 }
00495                 if( normalize ) {
00496                 //  it looks like it could work in 3D, but does not
00497                 double square_sum1 = 0., square_sum2 = 0.;
00498                 for ( int iz = 0; iz <= nz-1; ++iz) {
00499                         for ( int iy = 0; iy <= ny-1; ++iy) {
00500                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00501                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00502                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00503                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00504                                 }
00505                         }
00506                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00507                                 size_t ii = (iy  + iz * ny)* lsd2;
00508                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00509                                 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00510                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00511                                 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00512                         }
00513                         if(nx%2 == 0) {
00514                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00515                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00516                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00517                                         square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00518                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00519                                         square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00520                                 }
00521 
00522                         }
00523                         square_sum1 *= 2;
00524                         square_sum1 += x_data[0] * double(x_data[0]);
00525                         square_sum2 *= 2;
00526                         square_sum2 += y_data[0] * double(y_data[0]);
00527                         if(ny%2 == 0) {
00528                                 int ii = (ny/2  + iz * ny)* lsd2;
00529                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00530                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00531                         }
00532                         if(nx%2 == 0) {
00533                                 int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00534                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00535                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00536                                 if(ny%2 == 0) {
00537                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00538                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00539                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00540                                 }
00541                         }
00542                 }
00543                 result /= sqrt(square_sum1*square_sum2);
00544                 } else  result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
00545 
00546                 } else { //This 3D code is incorrect, but it is the best I can do now 01/09/06 PAP
00547                 int ky, kz;
00548                 int ny2 = ny/2; int nz2 = nz/2;
00549                 for ( int iz = 0; iz <= nz-1; ++iz) {
00550                         if(iz>nz2) kz=iz-nz; else kz=iz;
00551                         for ( int iy = 0; iy <= ny-1; ++iy) {
00552                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00553                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00554                                         // Skip Friedel related values
00555                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00556                                                 size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
00557                                                 result += x_data[ii] * double(y_data[ii]);
00558                                         }
00559                                 }
00560                         }
00561                 }
00562                 if( normalize ) {
00563                 //  still incorrect
00564                 double square_sum1 = 0., square_sum2 = 0.;
00565                 int ky, kz;
00566                 int ny2 = ny/2; int nz2 = nz/2;
00567                 for ( int iz = 0; iz <= nz-1; ++iz) {
00568                         if(iz>nz2) kz=iz-nz; else kz=iz;
00569                         for ( int iy = 0; iy <= ny-1; ++iy) {
00570                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00571                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00572                                         // Skip Friedel related values
00573                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00574                                                 size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
00575                                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00576                                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00577                                         }
00578                                 }
00579                         }
00580                 }
00581                 result /= sqrt(square_sum1*square_sum2);
00582                 } else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz/2);
00583                 }
00584         } else {
00585                 
00586                 size_t totsize = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00587 
00588                 double square_sum1 = 0., square_sum2 = 0.;
00589 
00590                 if (params.has_key("mask")) {
00591                         EMData* mask;
00592                         mask = params["mask"];
00593                         const float *const dm = mask->get_const_data();
00594                         if (normalize) {
00595                                 for (size_t i = 0; i < totsize; i++) {
00596                                         if (dm[i] > 0.5) {
00597                                                 square_sum1 += x_data[i]*double(x_data[i]);
00598                                                 square_sum2 += y_data[i]*double(y_data[i]);
00599                                                 result += x_data[i]*double(y_data[i]);
00600                                         }
00601                                 }
00602                         } else {
00603                                 for (size_t i = 0; i < totsize; i++) {
00604                                         if (dm[i] > 0.5) {
00605                                                 result += x_data[i]*double(y_data[i]);
00606                                                 n++;
00607                                         }
00608                                 }
00609                         }
00610                 } else {
00611                         // this little bit of manual loop unrolling makes the dot product as fast as sqeuclidean with -O2
00612                         for (size_t i=0; i<totsize; i++) result+=x_data[i]*y_data[i];
00613 
00614                         if (normalize) {
00615                                 square_sum1 = image->get_attr("square_sum");
00616                                 square_sum2 = with->get_attr("square_sum");
00617                         } else n = totsize;
00618                 }
00619                 if (normalize) result /= (sqrt(square_sum1*square_sum2)); else result /= n;
00620         }
00621 
00622 
00623         EXITFUNC;
00624         return (float) (negative*result);
00625 }

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

Implements EMAN::Cmp.

Definition at line 282 of file cmp.h.

00283                 {
00284                         return "Dot product (default -1 * dot product)";
00285                 }

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.

00278                 {
00279                         return NAME;
00280                 }

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::TypeDict::put().

00293                 {
00294                         TypeDict d;
00295                         d.put("negative", EMObject::INT, "If set, returns -1 * dot product. Set by default so smaller is better");
00296                         d.put("normalize", EMObject::INT, "If set, returns normalized dot product (cosine of the angle) -1.0 - 1.0.");
00297                         d.put("mask", EMObject::EMDATA, "image mask");
00298                         return d;
00299                 }

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

Definition at line 287 of file cmp.h.

00288                 {
00289                         return new DotCmp();
00290                 }


Member Data Documentation

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

Definition at line 51 of file cmp.cpp.


The documentation for this class was generated from the following files:
Generated on Tue Jul 12 13:51:18 2011 for EMAN2 by  doxygen 1.3.9.1