#include <cmp.h>
Inheritance diagram for EMAN::DotCmp:
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 Cmp * | NEW () |
Static Public Attributes | |
static const string | NAME = "dot" |
// Added mask option PAP 04/23/06 For complex images, it does not check r/i vs a/p.
negative | Returns -1 * dot product, default true | |
normalize | Returns normalized dot product -1.0 - 1.0 |
Definition at line 272 of file cmp.h.
To compare 'image' with another image passed in through its parameters.
An optional transformation may be used to transform the 2 images.
image | The first image to be compared. | |
with | The second image to be comppared. |
Implements EMAN::Cmp.
Definition at line 406 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().
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] |
string EMAN::DotCmp::get_name | ( | ) | const [inline, virtual] |
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.
Implements EMAN::Cmp.
Definition at line 292 of file cmp.h.
References EMAN::EMObject::EMDATA, EMAN::EMObject::INT, and 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 }
static Cmp* EMAN::DotCmp::NEW | ( | ) | [inline, static] |
const string DotCmp::NAME = "dot" [static] |