EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all 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]

Public Member Functions

float cmp (EMData *image, EMData *with) const
 To compare 'image' with another image passed in through its parameters. More...
 
string get_name () const
 Get the Cmp's name. More...
 
string get_desc () const
 
TypeDict get_param_types () const
 Get Cmp parameter information in a dictionary. More...
 
- Public Member Functions inherited from EMAN::Cmp
virtual ~Cmp ()
 
virtual Dict get_params () const
 Get the Cmp parameters in a key/value dictionary. More...
 
virtual void set_params (const Dict &new_params)
 Set the Cmp parameters using a key/value dictionary. More...
 

Static Public Member Functions

static Cmp * NEW ()
 

Static Public Attributes

static const string NAME = "dot"
 

Additional Inherited Members

- Protected Member Functions inherited from EMAN::Cmp
void validate_input_args (const EMData *image, const EMData *with) const
 
- Protected Attributes inherited from EMAN::Cmp
Dict params
 

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 (sludt.nosp@m.ke@b.nosp@m.cm.tm.nosp@m.c.ed.nosp@m.u)
Date
2005-07-13
Parameters
negativeReturns -1 * dot product, default true
normalizeReturns normalized dot product -1.0 - 1.0

Definition at line 268 of file cmp.h.

Member Function Documentation

◆ cmp()

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

418{
419 ENTERFUNC;
420
421 validate_input_args(image, with);
422
423 int normalize = params.set_default("normalize", 0);
424 float negative = (float)params.set_default("negative", 1);
425 if (negative) negative=-1.0; else negative=1.0;
426#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!!!!)
427 if(image->is_complex() && with->is_complex()) {
428 } else {
429 if (image->getcudarwdata() && with->getcudarwdata()) {
430 //cout << "CUDA dot cmp" << endl;
431 float* maskdata = 0;
432 bool has_mask = false;
433 EMData* mask = 0;
434 if (params.has_key("mask")) {
435 mask = params["mask"];
436 if(mask!=0) {has_mask=true;}
437 }
438 if(has_mask && !mask->getcudarwdata()){
439 mask->copy_to_cuda();
440 maskdata = mask->getcudarwdata();
441 }
442
443 float result = dot_cmp_cuda(image->getcudarwdata(), with->getcudarwdata(), maskdata, image->get_xsize(), image->get_ysize(), image->get_zsize());
444 result *= negative;
445
446 return result;
447
448 }
449 }
450#endif
451 const float *const x_data = image->get_const_data();
452 const float *const y_data = with->get_const_data();
453
454 double result = 0.;
455 long n = 0;
456 if(image->is_complex() && with->is_complex()) {
457 // Implemented by PAP 01/09/06 - please do not change. If in doubts, write/call me.
458 int nx = with->get_xsize();
459 int ny = with->get_ysize();
460 int nz = with->get_zsize();
461 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
462 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
463
464 int ixb = 2*((nx+1)%2);
465 int iyb = ny%2;
466 //
467 if(nz == 1) {
468 // it looks like it could work in 3D, but does not
469 for ( int iz = 0; iz <= nz-1; ++iz) {
470 double part = 0.;
471 for ( int iy = 0; iy <= ny-1; ++iy) {
472 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
473 size_t ii = ix + (iy + iz * ny)* lsd2;
474 part += x_data[ii] * double(y_data[ii]);
475 }
476 }
477 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
478 size_t ii = (iy + iz * ny)* lsd2;
479 part += x_data[ii] * double(y_data[ii]);
480 part += x_data[ii+1] * double(y_data[ii+1]);
481 }
482 if(nx%2 == 0) {
483 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
484 size_t ii = lsd2 - 2 + (iy + iz * ny)* lsd2;
485 part += x_data[ii] * double(y_data[ii]);
486 part += x_data[ii+1] * double(y_data[ii+1]);
487 }
488
489 }
490 part *= 2;
491 part += x_data[0] * double(y_data[0]);
492 if(ny%2 == 0) {
493 size_t ii = (ny/2 + iz * ny)* lsd2;
494 part += x_data[ii] * double(y_data[ii]);
495 }
496 if(nx%2 == 0) {
497 size_t ii = lsd2 - 2 + (0 + iz * ny)* lsd2;
498 part += x_data[ii] * double(y_data[ii]);
499 if(ny%2 == 0) {
500 int ii = lsd2 - 2 +(ny/2 + iz * ny)* lsd2;
501 part += x_data[ii] * double(y_data[ii]);
502 }
503 }
504 result += part;
505 }
506 } else {
507 // The 3D code
508 int ky, kz;
509 int ny2 = ny/2; int nz2 = nz/2;
510 for ( int iz = 0; iz <= nz-1; ++iz) {
511 if(iz>nz2) kz=iz-nz; else kz=iz;
512 for ( int iy = 0; iy <= ny-1; ++iy) {
513 if(iy>ny2) ky=iy-ny; else ky=iy;
514 for (int ix = 0; ix <= lsd2-1; ++ix) {
515 int p = 2;
516 if (ix <= 1 || (ix >= lsd2-2 && nx%2 == 0) ) p = 1;
517 size_t ii = ix + (iy + iz * ny)* (size_t)lsd2;
518 result += p * x_data[ii] * double(y_data[ii]);
519 }
520 }
521 }
522 if( normalize ) {
523 double square_sum1 = 0., square_sum2 = 0.;
524 int ky, kz;
525 int ny2 = ny/2; int nz2 = nz/2;
526 for ( int iz = 0; iz <= nz-1; ++iz) {
527 if(iz>nz2) kz=iz-nz; else kz=iz;
528 for ( int iy = 0; iy <= ny-1; ++iy) {
529 if(iy>ny2) ky=iy-ny; else ky=iy;
530 for (int ix = 0; ix <= lsd2-1; ++ix) {
531 int p = 2;
532 if (ix <= 1 || (ix >= lsd2-2 && nx%2 == 0) ) p = 1;
533 size_t ii = ix + (iy + iz * ny)* (size_t)lsd2;
534 square_sum1 += p * x_data[ii] * double(x_data[ii]);
535 square_sum2 += p * y_data[ii] * double(y_data[ii]);
536 }
537 }
538 }
539 result /= sqrt(square_sum1*square_sum2);
540 } else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
541 if( normalize ) {
542 // it looks like it could work in 3D, but does not
543 double square_sum1 = 0., square_sum2 = 0.;
544 for ( int iz = 0; iz <= nz-1; ++iz) {
545 for ( int iy = 0; iy <= ny-1; ++iy) {
546 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
547 size_t ii = ix + (iy + iz * ny)* lsd2;
548 square_sum1 += x_data[ii] * double(x_data[ii]);
549 square_sum2 += y_data[ii] * double(y_data[ii]);
550 }
551 }
552 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
553 size_t ii = (iy + iz * ny)* lsd2;
554 square_sum1 += x_data[ii] * double(x_data[ii]);
555 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
556 square_sum2 += y_data[ii] * double(y_data[ii]);
557 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
558 }
559 if(nx%2 == 0) {
560 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
561 size_t ii = lsd2 - 2 + (iy + iz * ny)* lsd2;
562 square_sum1 += x_data[ii] * double(x_data[ii]);
563 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
564 square_sum2 += y_data[ii] * double(y_data[ii]);
565 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
566 }
567
568 }
569 square_sum1 *= 2;
570 square_sum1 += x_data[0] * double(x_data[0]);
571 square_sum2 *= 2;
572 square_sum2 += y_data[0] * double(y_data[0]);
573 if(ny%2 == 0) {
574 int ii = (ny/2 + iz * ny)* lsd2;
575 square_sum1 += x_data[ii] * double(x_data[ii]);
576 square_sum2 += y_data[ii] * double(y_data[ii]);
577 }
578 if(nx%2 == 0) {
579 int ii = lsd2 - 2 + (0 + iz * ny)* lsd2;
580 square_sum1 += x_data[ii] * double(x_data[ii]);
581 square_sum2 += y_data[ii] * double(y_data[ii]);
582 if(ny%2 == 0) {
583 int ii = lsd2 - 2 +(ny/2 + iz * ny)* lsd2;
584 square_sum1 += x_data[ii] * double(x_data[ii]);
585 square_sum2 += y_data[ii] * double(y_data[ii]);
586 }
587 }
588 }
589 result /= sqrt(square_sum1*square_sum2);
590 } else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
591 // This concludes the 2D part.
592 }
593 } else {
594 // This part is for when two images are real, which is much easier because 2-D or 3-D
595 // is the same code.
596 size_t totsize = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
597
598 double square_sum1 = 0., square_sum2 = 0.;
599
600 if (params.has_key("mask")) {
601 EMData* mask;
602 mask = params["mask"];
603 const float *const dm = mask->get_const_data();
604 if (normalize) {
605 for (size_t i = 0; i < totsize; i++) {
606 if (dm[i] > 0.5) {
607 square_sum1 += x_data[i]*double(x_data[i]);
608 square_sum2 += y_data[i]*double(y_data[i]);
609 result += x_data[i]*double(y_data[i]);
610 }
611 }
612 } else {
613 for (size_t i = 0; i < totsize; i++) {
614 if (dm[i] > 0.5) {
615 result += x_data[i]*double(y_data[i]);
616 n++;
617 }
618 }
619 }
620 } else {
621 // this little bit of manual loop unrolling makes the dot product as fast as sqeuclidean with -O2
622 for (size_t i=0; i<totsize; i++) result += x_data[i]*y_data[i];
623
624 if (normalize) {
625 square_sum1 = image->get_attr("square_sum");
626 square_sum2 = with->get_attr("square_sum");
627 } else n = totsize;
628 }
629 if (normalize) result /= (sqrt(square_sum1*square_sum2)); else result /= n;
630 }
631
632
633 EXITFUNC;
634 return (float) (negative*result);
635}
EMAN::Cmp::validate_input_args
void validate_input_args(const EMData *image, const EMData *with) const
Definition: cmp.cpp:81
EMAN::Cmp::params
Dict params
Definition: cmp.h:132
EMAN::Dict::set_default
type set_default(const string &key, type val)
Default setting behavior This can be achieved using a template - d.woolford Jan 2008 (before there wa...
Definition: emobject.h:569
EMAN::Dict::has_key
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
dot_cmp_cuda
float dot_cmp_cuda(float *data1, float *data2, const float *dm, const int &nx, const int &ny, const int &nz)
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
ENTERFUNC
#define ENTERFUNC
Definition: log.h:48
EXITFUNC
#define EXITFUNC
Definition: log.h:49
dm
#define dm(i)
Definition: projector.cpp:1606

References dm, dot_cmp_cuda(), ENTERFUNC, EXITFUNC, EMAN::Dict::has_key(), EMAN::Cmp::params, EMAN::Dict::set_default(), sqrt(), and EMAN::Cmp::validate_input_args().

◆ get_desc()

string EMAN::DotCmp::get_desc ( ) const
inlinevirtual

Implements EMAN::Cmp.

Definition at line 278 of file cmp.h.

279 {
280 return "Dot product (default -1 * dot product)";
281 }

◆ get_name()

string EMAN::DotCmp::get_name ( ) const
inlinevirtual

Get the Cmp's name.

Each Cmp is identified by a unique name.

Returns
The Cmp's name.

Implements EMAN::Cmp.

Definition at line 273 of file cmp.h.

274 {
275 return NAME;
276 }
EMAN::DotCmp::NAME
static const string NAME
Definition: cmp.h:297

References NAME.

◆ get_param_types()

TypeDict EMAN::DotCmp::get_param_types ( ) const
inlinevirtual

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

289 {
290 TypeDict d;
291 d.put("negative", EMObject::INT, "If set, returns -1 * dot product. Set by default so smaller is better");
292 d.put("normalize", EMObject::INT, "If set, returns normalized dot product (cosine of the angle) -1.0 - 1.0.");
293 d.put("mask", EMObject::EMDATA, "image mask");
294 return d;
295 }

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

◆ NEW()

static Cmp * EMAN::DotCmp::NEW ( )
inlinestatic

Definition at line 283 of file cmp.h.

284 {
285 return new DotCmp();
286 }

Member Data Documentation

◆ NAME

const string DotCmp::NAME = "dot"
static

Definition at line 297 of file cmp.h.

Referenced by get_name().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.2