#include <cmp.h>


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 () |
// 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 231 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 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] |
| string EMAN::DotCmp::get_desc | ( | ) | const [inline, virtual] |
| static Cmp* EMAN::DotCmp::NEW | ( | ) | [inline, static] |
| 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 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 }
1.5.6