00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "cmp.h"
00037 #include "emdata.h"
00038 #include "ctf.h"
00039 #include "plugins/cmp_template.h"
00040
00041 #ifdef EMAN2_USING_CUDA
00042
00043 #include "cuda/cuda_processor.h"
00044 #include "cuda/cuda_cmp.h"
00045 #endif // EMAN2_USING_CUDA
00046
00047 using namespace EMAN;
00048
00049 const string CccCmp::NAME = "ccc";
00050 const string LodCmp::NAME = "lod";
00051 const string SqEuclideanCmp::NAME = "sqeuclidean";
00052 const string DotCmp::NAME = "dot";
00053 const string TomoCccCmp::NAME = "ccc.tomo";
00054 const string QuadMinDotCmp::NAME = "quadmindot";
00055 const string OptVarianceCmp::NAME = "optvariance";
00056 const string PhaseCmp::NAME = "phase";
00057 const string FRCCmp::NAME = "frc";
00058
00059 template <> Factory < Cmp >::Factory()
00060 {
00061 force_add<CccCmp>();
00062 force_add<LodCmp>();
00063 force_add<SqEuclideanCmp>();
00064 force_add<DotCmp>();
00065 force_add<TomoCccCmp>();
00066 force_add<QuadMinDotCmp>();
00067 force_add<OptVarianceCmp>();
00068 force_add<PhaseCmp>();
00069 force_add<FRCCmp>();
00070
00071 }
00072
00073 void Cmp::validate_input_args(const EMData * image, const EMData *with) const
00074 {
00075
00076 #ifdef EMAN2_USING_CUDA
00077 if (image->getcudarwdata() && with->getcudarwdata()) {
00078
00079 return;
00080 }
00081 #endif
00082 if (!image) {
00083 throw NullPointerException("compared image");
00084 }
00085 if (!with) {
00086 throw NullPointerException("compare-with image");
00087 }
00088
00089 if (!EMUtil::is_same_size(image, with)) {
00090 throw ImageFormatException( "images not same size");
00091 }
00092
00093 float *d1 = image->get_data();
00094 if (!d1) {
00095 throw NullPointerException("image contains no data");
00096 }
00097
00098 float *d2 = with->get_data();
00099 if (!d2) {
00100 throw NullPointerException("compare-with image data");
00101 }
00102 }
00103
00104
00105 float CccCmp::cmp(EMData * image, EMData *with) const
00106 {
00107 ENTERFUNC;
00108 if (image->is_complex() || with->is_complex())
00109 throw ImageFormatException( "Complex images not supported by CMP::CccCmp");
00110 validate_input_args(image, with);
00111
00112 const float *const d1 = image->get_const_data();
00113 const float *const d2 = with->get_const_data();
00114
00115 float negative = (float)params.set_default("negative", 1);
00116 if (negative) negative=-1.0; else negative=1.0;
00117
00118 double avg1 = 0.0, var1 = 0.0, avg2 = 0.0, var2 = 0.0, ccc = 0.0;
00119 long n = 0;
00120 size_t totsize = image->get_xsize()*image->get_ysize()*image->get_zsize();
00121
00122 bool has_mask = false;
00123 EMData* mask = 0;
00124 if (params.has_key("mask")) {
00125 mask = params["mask"];
00126 if(mask!=0) {has_mask=true;}
00127 }
00128 #ifdef EMAN2_USING_CUDA
00129 if (image->getcudarwdata() && with->getcudarwdata()) {
00130
00131 float* maskdata = 0;
00132 if(has_mask && !mask->getcudarwdata()){
00133 mask->copy_to_cuda();
00134 maskdata = mask->getcudarwdata();
00135 }
00136 float ccc = ccc_cmp_cuda(image->getcudarwdata(), with->getcudarwdata(), maskdata, image->get_xsize(), image->get_ysize(), image->get_zsize());
00137 ccc *= negative;
00138
00139 return ccc;
00140 }
00141 #endif
00142 if (has_mask) {
00143 const float *const dm = mask->get_const_data();
00144 for (size_t i = 0; i < totsize; ++i) {
00145 if (dm[i] > 0.5) {
00146 avg1 += double(d1[i]);
00147 var1 += d1[i]*double(d1[i]);
00148 avg2 += double(d2[i]);
00149 var2 += d2[i]*double(d2[i]);
00150 ccc += d1[i]*double(d2[i]);
00151 n++;
00152 }
00153 }
00154 } else {
00155 for (size_t i = 0; i < totsize; ++i) {
00156 avg1 += double(d1[i]);
00157 var1 += d1[i]*double(d1[i]);
00158 avg2 += double(d2[i]);
00159 var2 += d2[i]*double(d2[i]);
00160 ccc += d1[i]*double(d2[i]);
00161 }
00162 n = totsize;
00163 }
00164
00165 avg1 /= double(n);
00166 var1 = var1/double(n) - avg1*avg1;
00167 avg2 /= double(n);
00168 var2 = var2/double(n) - avg2*avg2;
00169 ccc = ccc/double(n) - avg1*avg2;
00170 ccc /= sqrt(var1*var2);
00171 ccc *= negative;
00172 return static_cast<float>(ccc);
00173 EXITFUNC;
00174 }
00175
00176
00177
00178
00179 float LodCmp::cmp(EMData * image, EMData *with) const
00180 {
00181 ENTERFUNC;
00182 if (image->is_complex() || with->is_complex())
00183 throw ImageFormatException( "Complex images not yet supported by CMP::LodCmp");
00184 validate_input_args(image, with);
00185
00186 const float *const d1 = image->get_const_data();
00187 const float *const d2 = with->get_const_data();
00188
00189 float negative = (float)params.set_default("negative", 1);
00190 if (negative) negative=-1.0; else negative=1.0;
00191
00192 double avg1 = 0.0, sig1 = 0.0, avg2 = 0.0, sig2 = 0.0, lod = 0.0;
00193 long n = 0;
00194 size_t totsize = image->get_xsize()*image->get_ysize()*image->get_zsize();
00195 size_t i;
00196
00197 bool has_mask = false;
00198 EMData* mask = 0;
00199 if (params.has_key("mask")) {
00200 mask = params["mask"];
00201 if(mask!=0) {has_mask=true;}
00202 }
00203
00204 int normalize = 0;
00205 if (params.has_key("normalize")) {
00206 normalize = params["normalize"];
00207 }
00208
00209 if (normalize) {
00210 if (has_mask) {
00211 const float *const dm = mask->get_const_data();
00212 for (i = 0; i < totsize; ++i) {
00213 if (dm[i] > 0.5) {
00214 avg1 += double(d1[i]);
00215 avg2 += double(d2[i]);
00216 n++;
00217 }
00218 }
00219 } else {
00220 for (i = 0; i < totsize; ++i) {
00221 avg1 += double(d1[i]);
00222 avg2 += double(d2[i]);
00223 }
00224 n = totsize;
00225 }
00226
00227 avg1 /= double(n);
00228 avg2 /= double(n);
00229
00230 if (has_mask) {
00231 const float *const dm = mask->get_const_data();
00232 for (i = 0; i < totsize; ++i) {
00233 if (dm[i] > 0.5) {
00234 sig1 += fabs(double(d1[i])-avg1);
00235 sig2 += fabs(double(d2[i])-avg2);
00236 }
00237 }
00238 } else {
00239 for (i = 0; i < totsize; ++i) {
00240 sig1 += fabs(double(d1[i])-avg1);
00241 sig2 += fabs(double(d2[i])-avg2);
00242 }
00243 }
00244 } else {
00245 avg1 = 0.0; avg2 = 0.0;
00246 sig1 = 1.0; sig2 = 1.0;
00247 }
00248
00249 if (has_mask) {
00250 const float *const dm = mask->get_const_data();
00251 for (i = 0; i < totsize; ++i) {
00252 if (dm[i] > 0.5) {
00253 lod += fabs((double(d1[i])-avg1)/sig1 - (double(d2[i])-avg2)/sig2);
00254 }
00255 }
00256 } else {
00257 for (i = 0; i < totsize; ++i) {
00258 lod += fabs((double(d1[i])-avg1)/sig1 - (double(d2[i])-avg2)/sig2);
00259 }
00260 }
00261
00262 lod *= (-0.5);
00263 lod *= negative;
00264 return static_cast<float>(lod);
00265 EXITFUNC;
00266 }
00267
00268
00269
00270 float SqEuclideanCmp::cmp(EMData *image,EMData * withorig ) const
00271 {
00272 ENTERFUNC;
00273 EMData *with = withorig;
00274 validate_input_args(image, with);
00275
00276 int zeromask = params.set_default("zeromask",0);
00277 int normto = params.set_default("normto",0);
00278
00279 if (normto) {
00280 if (zeromask) with = withorig->process("normalize.toimage",Dict("to",image));
00281 else with = withorig->process("normalize.toimage",Dict("to",image,"ignore_zero",0));
00282 with->set_attr("deleteme",1);
00283 if ((float)(with->get_attr("norm_mult"))<=0) {
00284 delete with;
00285 with=withorig;
00286 }
00287 }
00288
00289 const float *const y_data = with->get_const_data();
00290 const float *const x_data = image->get_const_data();
00291 double result = 0.;
00292 float n = 0.0f;
00293 if(image->is_complex() && with->is_complex()) {
00294
00295 int nx = with->get_xsize();
00296 int ny = with->get_ysize();
00297 int nz = with->get_zsize();
00298 nx = (nx - 2 + with->is_fftodd());
00299 int lsd2 = (nx + 2 - nx%2) ;
00300
00301 int ixb = 2*((nx+1)%2);
00302 int iyb = ny%2;
00303
00304 if(nz == 1) {
00305
00306 for ( int iz = 0; iz <= nz-1; iz++) {
00307 double part = 0.;
00308 for ( int iy = 0; iy <= ny-1; iy++) {
00309 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ix++) {
00310 size_t ii = ix + (iy + iz * ny)* lsd2;
00311 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00312 }
00313 }
00314 for ( int iy = 1; iy <= ny/2-1 + iyb; iy++) {
00315 size_t ii = (iy + iz * ny)* lsd2;
00316 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00317 part += (x_data[ii+1] - y_data[ii+1])*double(x_data[ii+1] - y_data[ii+1]);
00318 }
00319 if(nx%2 == 0) {
00320 for ( int iy = 1; iy <= ny/2-1 + iyb; iy++) {
00321 size_t ii = lsd2 - 2 + (iy + iz * ny)* lsd2;
00322 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00323 part += (x_data[ii+1] - y_data[ii+1])*double(x_data[ii+1] - y_data[ii+1]);
00324 }
00325
00326 }
00327 part *= 2;
00328 part += (x_data[0] - y_data[0])*double(x_data[0] - y_data[0]);
00329 if(ny%2 == 0) {
00330 int ii = (ny/2 + iz * ny)* lsd2;
00331 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00332 }
00333 if(nx%2 == 0) {
00334 int ii = lsd2 - 2 + (0 + iz * ny)* lsd2;
00335 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00336 if(ny%2 == 0) {
00337 int ii = lsd2 - 2 +(ny/2 + iz * ny)* lsd2;
00338 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00339 }
00340 }
00341 result += part;
00342 }
00343 n = (float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz;
00344
00345 }else{
00346 int ky, kz;
00347 int ny2 = ny/2; int nz2 = nz/2;
00348 for ( int iz = 0; iz <= nz-1; iz++) {
00349 if(iz>nz2) kz=iz-nz; else kz=iz;
00350 for ( int iy = 0; iy <= ny-1; iy++) {
00351 if(iy>ny2) ky=iy-ny; else ky=iy;
00352 for ( int ix = 0; ix <= lsd2-1; ix++) {
00353
00354 if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00355 size_t ii = ix + (iy + iz * ny)* (size_t)lsd2;
00356 result += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00357 }
00358 }
00359 }
00360 }
00361 n = ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz)/2.0f;
00362 }
00363 } else {
00364 size_t totsize = (size_t)image->get_xsize()*image->get_ysize()*image->get_zsize();
00365 if (params.has_key("mask")) {
00366 EMData* mask;
00367 mask = params["mask"];
00368 const float *const dm = mask->get_const_data();
00369 for (size_t i = 0; i < totsize; i++) {
00370 if (dm[i] > 0.5) {
00371 double temp = x_data[i]- y_data[i];
00372 result += temp*temp;
00373 n++;
00374 }
00375 }
00376 }
00377 else if (zeromask) {
00378 n=0;
00379 for (size_t i = 0; i < totsize; i++) {
00380 if (x_data[i]==0 || y_data[i]==0) continue;
00381 double temp = x_data[i]- y_data[i];
00382 result += temp*temp;
00383 n++;
00384 }
00385
00386 }
00387 else {
00388 for (size_t i = 0; i < totsize; i++) {
00389 double temp = x_data[i]- y_data[i];
00390 result += temp*temp;
00391 }
00392 n = (float)totsize;
00393 }
00394 }
00395 result/=n;
00396
00397 EXITFUNC;
00398 if (with->has_attr("deleteme")) delete with;
00399 float ret = (float)result;
00400 if (!Util::goodf(&ret)) return FLT_MAX;
00401 return ret;
00402 }
00403
00404
00405
00406
00407 float DotCmp::cmp(EMData* image, EMData* with) const
00408 {
00409 ENTERFUNC;
00410
00411 validate_input_args(image, with);
00412
00413 int normalize = params.set_default("normalize", 0);
00414 float negative = (float)params.set_default("negative", 1);
00415 if (negative) negative=-1.0; else negative=1.0;
00416 #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!!!!)
00417 if(image->is_complex() && with->is_complex()) {
00418 } else {
00419 if (image->getcudarwdata() && with->getcudarwdata()) {
00420
00421 float* maskdata = 0;
00422 bool has_mask = false;
00423 EMData* mask = 0;
00424 if (params.has_key("mask")) {
00425 mask = params["mask"];
00426 if(mask!=0) {has_mask=true;}
00427 }
00428 if(has_mask && !mask->getcudarwdata()){
00429 mask->copy_to_cuda();
00430 maskdata = mask->getcudarwdata();
00431 }
00432
00433 float result = dot_cmp_cuda(image->getcudarwdata(), with->getcudarwdata(), maskdata, image->get_xsize(), image->get_ysize(), image->get_zsize());
00434 result *= negative;
00435
00436 return result;
00437
00438 }
00439 }
00440 #endif
00441 const float *const x_data = image->get_const_data();
00442 const float *const y_data = with->get_const_data();
00443
00444 double result = 0.;
00445 long n = 0;
00446 if(image->is_complex() && with->is_complex()) {
00447
00448 int nx = with->get_xsize();
00449 int ny = with->get_ysize();
00450 int nz = with->get_zsize();
00451 nx = (nx - 2 + with->is_fftodd());
00452 int lsd2 = (nx + 2 - nx%2) ;
00453
00454 int ixb = 2*((nx+1)%2);
00455 int iyb = ny%2;
00456
00457 if(nz == 1) {
00458
00459 for ( int iz = 0; iz <= nz-1; ++iz) {
00460 double part = 0.;
00461 for ( int iy = 0; iy <= ny-1; ++iy) {
00462 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00463 size_t ii = ix + (iy + iz * ny)* lsd2;
00464 part += x_data[ii] * double(y_data[ii]);
00465 }
00466 }
00467 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00468 size_t ii = (iy + iz * ny)* lsd2;
00469 part += x_data[ii] * double(y_data[ii]);
00470 part += x_data[ii+1] * double(y_data[ii+1]);
00471 }
00472 if(nx%2 == 0) {
00473 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00474 size_t ii = lsd2 - 2 + (iy + iz * ny)* lsd2;
00475 part += x_data[ii] * double(y_data[ii]);
00476 part += x_data[ii+1] * double(y_data[ii+1]);
00477 }
00478
00479 }
00480 part *= 2;
00481 part += x_data[0] * double(y_data[0]);
00482 if(ny%2 == 0) {
00483 size_t ii = (ny/2 + iz * ny)* lsd2;
00484 part += x_data[ii] * double(y_data[ii]);
00485 }
00486 if(nx%2 == 0) {
00487 size_t ii = lsd2 - 2 + (0 + iz * ny)* lsd2;
00488 part += x_data[ii] * double(y_data[ii]);
00489 if(ny%2 == 0) {
00490 int ii = lsd2 - 2 +(ny/2 + iz * ny)* lsd2;
00491 part += x_data[ii] * double(y_data[ii]);
00492 }
00493 }
00494 result += part;
00495 }
00496 if( normalize ) {
00497
00498 double square_sum1 = 0., square_sum2 = 0.;
00499 for ( int iz = 0; iz <= nz-1; ++iz) {
00500 for ( int iy = 0; iy <= ny-1; ++iy) {
00501 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00502 size_t ii = ix + (iy + iz * ny)* lsd2;
00503 square_sum1 += x_data[ii] * double(x_data[ii]);
00504 square_sum2 += y_data[ii] * double(y_data[ii]);
00505 }
00506 }
00507 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00508 size_t ii = (iy + iz * ny)* lsd2;
00509 square_sum1 += x_data[ii] * double(x_data[ii]);
00510 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00511 square_sum2 += y_data[ii] * double(y_data[ii]);
00512 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00513 }
00514 if(nx%2 == 0) {
00515 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00516 size_t ii = lsd2 - 2 + (iy + iz * ny)* lsd2;
00517 square_sum1 += x_data[ii] * double(x_data[ii]);
00518 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00519 square_sum2 += y_data[ii] * double(y_data[ii]);
00520 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00521 }
00522
00523 }
00524 square_sum1 *= 2;
00525 square_sum1 += x_data[0] * double(x_data[0]);
00526 square_sum2 *= 2;
00527 square_sum2 += y_data[0] * double(y_data[0]);
00528 if(ny%2 == 0) {
00529 int ii = (ny/2 + iz * ny)* lsd2;
00530 square_sum1 += x_data[ii] * double(x_data[ii]);
00531 square_sum2 += y_data[ii] * double(y_data[ii]);
00532 }
00533 if(nx%2 == 0) {
00534 int ii = lsd2 - 2 + (0 + iz * ny)* lsd2;
00535 square_sum1 += x_data[ii] * double(x_data[ii]);
00536 square_sum2 += y_data[ii] * double(y_data[ii]);
00537 if(ny%2 == 0) {
00538 int ii = lsd2 - 2 +(ny/2 + iz * ny)* lsd2;
00539 square_sum1 += x_data[ii] * double(x_data[ii]);
00540 square_sum2 += y_data[ii] * double(y_data[ii]);
00541 }
00542 }
00543 }
00544 result /= sqrt(square_sum1*square_sum2);
00545 } else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
00546
00547 } else {
00548 int ky, kz;
00549 int ny2 = ny/2; int nz2 = nz/2;
00550 for ( int iz = 0; iz <= nz-1; ++iz) {
00551 if(iz>nz2) kz=iz-nz; else kz=iz;
00552 for ( int iy = 0; iy <= ny-1; ++iy) {
00553 if(iy>ny2) ky=iy-ny; else ky=iy;
00554 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00555
00556 if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00557 size_t ii = ix + (iy + iz * ny)* (size_t)lsd2;
00558 result += x_data[ii] * double(y_data[ii]);
00559 }
00560 }
00561 }
00562 }
00563 if( normalize ) {
00564
00565 double square_sum1 = 0., square_sum2 = 0.;
00566 int ky, kz;
00567 int ny2 = ny/2; int nz2 = nz/2;
00568 for ( int iz = 0; iz <= nz-1; ++iz) {
00569 if(iz>nz2) kz=iz-nz; else kz=iz;
00570 for ( int iy = 0; iy <= ny-1; ++iy) {
00571 if(iy>ny2) ky=iy-ny; else ky=iy;
00572 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00573
00574 if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00575 size_t ii = ix + (iy + iz * ny)* (size_t)lsd2;
00576 square_sum1 += x_data[ii] * double(x_data[ii]);
00577 square_sum2 += y_data[ii] * double(y_data[ii]);
00578 }
00579 }
00580 }
00581 }
00582 result /= sqrt(square_sum1*square_sum2);
00583 } else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz/2);
00584 }
00585 } else {
00586
00587 size_t totsize = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00588
00589 double square_sum1 = 0., square_sum2 = 0.;
00590
00591 if (params.has_key("mask")) {
00592 EMData* mask;
00593 mask = params["mask"];
00594 const float *const dm = mask->get_const_data();
00595 if (normalize) {
00596 for (size_t i = 0; i < totsize; i++) {
00597 if (dm[i] > 0.5) {
00598 square_sum1 += x_data[i]*double(x_data[i]);
00599 square_sum2 += y_data[i]*double(y_data[i]);
00600 result += x_data[i]*double(y_data[i]);
00601 }
00602 }
00603 } else {
00604 for (size_t i = 0; i < totsize; i++) {
00605 if (dm[i] > 0.5) {
00606 result += x_data[i]*double(y_data[i]);
00607 n++;
00608 }
00609 }
00610 }
00611 } else {
00612
00613 for (size_t i=0; i<totsize; i++) result+=x_data[i]*y_data[i];
00614
00615 if (normalize) {
00616 square_sum1 = image->get_attr("square_sum");
00617 square_sum2 = with->get_attr("square_sum");
00618 } else n = totsize;
00619 }
00620 if (normalize) result /= (sqrt(square_sum1*square_sum2)); else result /= n;
00621 }
00622
00623
00624 EXITFUNC;
00625 return (float) (negative*result);
00626 }
00627
00628
00629
00630
00631
00632
00633
00634 float TomoCccCmp::cmp(EMData * image, EMData *with) const
00635 {
00636 ENTERFUNC;
00637 EMData* ccf = params.set_default("ccf",(EMData*) NULL);
00638 bool ccf_ownership = false;
00639 bool norm = params.set_default("norm",true);
00640 float negative = (float)params.set_default("negative", 1);
00641 if (negative) negative=-1.0; else negative=1.0;
00642 int searchx, searchy, searchz;
00643
00644
00645 searchx = params.set_default("searchx",-1);
00646 searchy = params.set_default("searchy",-1);
00647 searchz = params.set_default("searchz",-1);
00648
00649 #ifdef EMAN2_USING_CUDA
00650 if(image->getcudarwdata() && with->getcudarwdata()){
00651 if (!ccf) {
00652 ccf = image->calc_ccf(with);
00653 ccf_ownership = true;
00654 }
00655
00656 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
00657 float best_score = get_value_at_wrap_cuda(ccf->getcudarwdata(), 0, 0, 0, ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
00658 if(norm) {
00659 best_score = negative*(best_score - stats.x)/sqrt(stats.y);
00660 } else {
00661 best_score = negative*best_score;
00662 }
00663
00664 if (ccf_ownership) delete ccf; ccf = 0;
00665
00666 return best_score;
00667
00668 }
00669 #endif
00670
00671 if (!ccf) {
00672 ccf = image->calc_ccf(with);
00673 ccf_ownership = true;
00674 }
00675 if (norm) ccf->process_inplace("normalize");
00676
00677 float best_score = ccf->get_value_at_wrap(0,0,0);
00678 if (ccf_ownership) delete ccf; ccf = 0;
00679
00680 return negative*best_score;
00681
00682 }
00683
00684
00685
00686 float QuadMinDotCmp::cmp(EMData * image, EMData *with) const
00687 {
00688 ENTERFUNC;
00689 validate_input_args(image, with);
00690
00691 if (image->get_zsize()!=1) throw InvalidValueException(0, "QuadMinDotCmp supports 2D only");
00692
00693 int nx=image->get_xsize();
00694 int ny=image->get_ysize();
00695
00696 int normalize = params.set_default("normalize", 0);
00697 float negative = (float)params.set_default("negative", 1);
00698
00699 if (negative) negative=-1.0; else negative=1.0;
00700
00701 double result[4] = { 0,0,0,0 }, sq1[4] = { 0,0,0,0 }, sq2[4] = { 0,0,0,0 } ;
00702
00703 vector<int> image_saved_offsets = image->get_array_offsets();
00704 vector<int> with_saved_offsets = with->get_array_offsets();
00705 image->set_array_offsets(-nx/2,-ny/2);
00706 with->set_array_offsets(-nx/2,-ny/2);
00707 int i,x,y;
00708 for (y=-ny/2; y<ny/2; y++) {
00709 for (x=-nx/2; x<nx/2; x++) {
00710 int quad=(x<0?0:1) + (y<0?0:2);
00711 result[quad]+=(*image)(x,y)*(*with)(x,y);
00712 if (normalize) {
00713 sq1[quad]+=(*image)(x,y)*(*image)(x,y);
00714 sq2[quad]+=(*with)(x,y)*(*with)(x,y);
00715 }
00716 }
00717 }
00718 image->set_array_offsets(image_saved_offsets);
00719 with->set_array_offsets(with_saved_offsets);
00720
00721 if (normalize) {
00722 for (i=0; i<4; i++) result[i]/=sqrt(sq1[i]*sq2[i]);
00723 } else {
00724 for (i=0; i<4; i++) result[i]/=nx*ny/4;
00725 }
00726
00727 float worst=static_cast<float>(result[0]);
00728 for (i=1; i<4; i++) if (static_cast<float>(result[i])<worst) worst=static_cast<float>(result[i]);
00729
00730 EXITFUNC;
00731 return (float) (negative*worst);
00732 }
00733
00734 float OptVarianceCmp::cmp(EMData * image, EMData *with) const
00735 {
00736 ENTERFUNC;
00737 validate_input_args(image, with);
00738
00739 int keepzero = params.set_default("keepzero", 1);
00740 int invert = params.set_default("invert",0);
00741 int matchfilt = params.set_default("matchfilt",1);
00742 int matchamp = params.set_default("matchamp",0);
00743 int radweight = params.set_default("radweight",0);
00744 int dbug = params.set_default("debug",0);
00745
00746 size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00747
00748
00749 EMData *with2=NULL;
00750 if (matchfilt) {
00751 EMData *a = image->do_fft();
00752 EMData *b = with->do_fft();
00753
00754 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1);
00755 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1);
00756
00757 float avg=0;
00758 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00759 rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i]));
00760 avg+=rfa[i];
00761 }
00762
00763 avg/=a->get_ysize()/2.0f;
00764 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00765 if (rfa[i]>avg*10.0) rfa[i]=10.0;
00766 }
00767 rfa[0]=0.0;
00768
00769 if (dbug) b->write_image("a.hdf",-1);
00770
00771 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa);
00772 with2=b->do_ift();
00773
00774 if (dbug) b->write_image("a.hdf",-1);
00775 if (dbug) a->write_image("a.hdf",-1);
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788 delete a;
00789 delete b;
00790
00791 if (dbug) {
00792 with2->write_image("a.hdf",-1);
00793 image->write_image("a.hdf",-1);
00794 }
00795
00796
00797
00798 }
00799
00800
00801
00802 if (matchamp) {
00803 EMData *a = image->do_fft();
00804 EMData *b = with->do_fft();
00805 size_t size2 = (size_t)a->get_xsize() * a->get_ysize() * a->get_zsize();
00806
00807 a->ri2ap();
00808 b->ri2ap();
00809
00810 const float *const ad=a->get_const_data();
00811 float * bd=b->get_data();
00812
00813 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i];
00814 b->update();
00815
00816 b->ap2ri();
00817 with2=b->do_ift();
00818
00819 delete a;
00820 delete b;
00821 }
00822
00823 const float * x_data;
00824 if (with2) x_data=with2->get_const_data();
00825 else x_data = with->get_const_data();
00826 const float *const y_data = image->get_const_data();
00827
00828 size_t nx = image->get_xsize();
00829 float m = 0;
00830 float b = 0;
00831
00832
00833
00834 if (dbug) {
00835 FILE *out=fopen("dbug.optvar.txt","w");
00836 if (out) {
00837 for (size_t i=0; i<size; i++) {
00838 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]);
00839 }
00840 fclose(out);
00841 }
00842 }
00843
00844
00845 Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero);
00846 if (m == 0) {
00847 m = FLT_MIN;
00848 }
00849 b = -b / m;
00850 m = 1.0f / m;
00851
00852
00853
00854
00855
00856
00857
00858
00859 double result = 0;
00860 int count = 0;
00861
00862 if (radweight) {
00863 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only");
00864 if (keepzero) {
00865 for (size_t i = 0,y=0; i < size; y++) {
00866 for (size_t x=0; x<nx; i++,x++) {
00867 if (y_data[i] && x_data[i]) {
00868 #ifdef _WIN32
00869 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00870 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00871 #else
00872 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00873 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00874 #endif
00875 count++;
00876 }
00877 }
00878 }
00879 result/=count;
00880 }
00881 else {
00882 for (size_t i = 0,y=0; i < size; y++) {
00883 for (size_t x=0; x<nx; i++,x++) {
00884 #ifdef _WIN32
00885 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00886 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00887 #else
00888 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00889 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00890 #endif
00891 }
00892 }
00893 result = result / size;
00894 }
00895 }
00896 else {
00897 if (keepzero) {
00898 for (size_t i = 0; i < size; i++) {
00899 if (y_data[i] && x_data[i]) {
00900 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00901 else result += Util::square((x_data[i] * m) + b - y_data[i]);
00902 count++;
00903 }
00904 }
00905 result/=count;
00906 }
00907 else {
00908 for (size_t i = 0; i < size; i++) {
00909 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00910 else result += Util::square((x_data[i] * m) + b - y_data[i]);
00911 }
00912 result = result / size;
00913 }
00914 }
00915 scale = m;
00916 shift = b;
00917
00918 image->set_attr("ovcmp_m",m);
00919 image->set_attr("ovcmp_b",b);
00920 if (with2) delete with2;
00921 EXITFUNC;
00922
00923 #if 0
00924 return (1 - result);
00925 #endif
00926
00927 return static_cast<float>(result);
00928 }
00929
00930 float PhaseCmp::cmp(EMData * image, EMData *with) const
00931 {
00932 ENTERFUNC;
00933
00934 int snrweight = params.set_default("snrweight", 0);
00935 int snrfn = params.set_default("snrfn",0);
00936 int ampweight = params.set_default("ampweight",0);
00937 int zeromask = params.set_default("zeromask",0);
00938 float minres = params.set_default("minres",500.0f);
00939 float maxres = params.set_default("maxres",2.0f);
00940
00941 if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator");
00942
00943 EMData *image_fft = NULL;
00944 EMData *with_fft = NULL;
00945
00946 int ny = image->get_ysize();
00947
00948 int np = 0;
00949 vector<float> snr;
00950
00951
00952 if (snrweight) {
00953 Ctf *ctf = NULL;
00954 if (!image->has_attr("ctf")) {
00955 if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
00956 ctf=with->get_attr("ctf");
00957 }
00958 else ctf=image->get_attr("ctf");
00959
00960 float ds=1.0f/(ctf->apix*ny);
00961 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR);
00962 if(ctf) {delete ctf; ctf=0;}
00963 np=snr.size();
00964 }
00965
00966 else if (snrfn==1) {
00967 np=ny/2;
00968 snr.resize(np);
00969
00970 for (int i = 0; i < np; i++) {
00971 float x2 = 10.0f*i/np;
00972 snr[i] = x2 * exp(-x2);
00973 }
00974 }
00975
00976 else {
00977 np=ny/2;
00978 snr.resize(np);
00979
00980 for (int i = 0; i < np; i++) snr[i]=1.0;
00981 }
00982
00983
00984 float pmin,pmax;
00985 if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;
00986 else pmin=0;
00987 if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
00988 else pmax=0;
00989
00990
00991
00992
00993 for (int i = 0; i < np; i++) {
00994 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f;
00995 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f;
00996
00997 }
00998
00999 if (zeromask) {
01000 image_fft=image->copy();
01001 with_fft=with->copy();
01002
01003 if (image_fft->is_complex()) image_fft->do_ift_inplace();
01004 if (with_fft->is_complex()) with_fft->do_ift_inplace();
01005
01006 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize();
01007 float *d1=image_fft->get_data();
01008 float *d2=with_fft->get_data();
01009
01010 for (int i=0; i<sz; i++) {
01011 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
01012 }
01013
01014 image_fft->update();
01015 with_fft->update();
01016 image_fft->do_fft_inplace();
01017 with_fft->do_fft_inplace();
01018 image_fft->set_attr("free_me",1);
01019 with_fft->set_attr("free_me",1);
01020 }
01021 else {
01022 if (image->is_complex()) image_fft=image;
01023 else {
01024 image_fft=image->do_fft();
01025 image_fft->set_attr("free_me",1);
01026 }
01027
01028 if (with->is_complex()) with_fft=with;
01029 else {
01030 with_fft=with->do_fft();
01031 with_fft->set_attr("free_me",1);
01032 }
01033 }
01034
01035
01036
01037 const float *const image_fft_data = image_fft->get_const_data();
01038 const float *const with_fft_data = with_fft->get_const_data();
01039 double sum = 0;
01040 double norm = FLT_MIN;
01041 size_t i = 0;
01042
01043 ny=image_fft->get_ysize();
01044 int nz=image_fft->get_zsize();
01045 int nx2=image_fft->get_xsize()/2;
01046 int ny2=image_fft->get_ysize()/2;
01047
01048
01049
01050 if (np==0) {
01051 for (int z = 0; z < nz; z++){
01052 for (int y = 0; y < ny; y++) {
01053 for (int x = 0; x < nx2; x ++) {
01054 float a;
01055 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01056 else a=1.0f;
01057 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01058 norm += a;
01059 i += 2;
01060 }
01061 }
01062 }
01063
01064 }
01065 else if (nz==1) {
01066 for (int y = 0; y < ny; y++) {
01067 for (int x = 0; x < nx2; x ++) {
01068 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y);
01069 if (r>=ny2) { i+=2; continue; }
01070
01071 float a;
01072 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01073 else a=1.0f;
01074 a*=snr[r];
01075 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01076 norm += a;
01077 i += 2;
01078 }
01079 }
01080 }
01081 else {
01082 for (int z = 0; z < nz; z++){
01083 for (int y = 0; y < ny; y++) {
01084 for (int x = 0; x < nx2; x ++) {
01085 int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z);
01086 if (r>=ny2) { i+=2; continue; }
01087
01088 float a;
01089 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01090 else a=1.0f;
01091 a*=snr[r];
01092 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01093 norm += a;
01094 i += 2;
01095 }
01096 }
01097 }
01098
01099 }
01100
01101 EXITFUNC;
01102
01103 if( image_fft->has_attr("free_me") )
01104 {
01105 delete image_fft;
01106 image_fft = 0;
01107 }
01108 if( with_fft->has_attr("free_me") )
01109 {
01110 delete with_fft;
01111 with_fft = 0;
01112 }
01113 #if 0
01114 return (1.0f - sum / norm);
01115 #endif
01116 return (float)(sum / norm);
01117 }
01118
01119 float FRCCmp::cmp(EMData * image, EMData * with) const
01120 {
01121 ENTERFUNC;
01122 validate_input_args(image, with);
01123
01124 int snrweight = params.set_default("snrweight", 0);
01125 int ampweight = params.set_default("ampweight", 0);
01126 int sweight = params.set_default("sweight", 1);
01127 int nweight = params.set_default("nweight", 0);
01128 int zeromask = params.set_default("zeromask",0);
01129 float minres = params.set_default("minres",500.0f);
01130 float maxres = params.set_default("maxres",2.0f);
01131
01132 vector < float >fsc;
01133 bool use_cpu = true;
01134
01135 if (use_cpu) {
01136 if (zeromask) {
01137 image=image->copy();
01138 with=with->copy();
01139
01140 int sz=image->get_xsize()*image->get_ysize()*image->get_zsize();
01141 float *d1=image->get_data();
01142 float *d2=with->get_data();
01143
01144 for (int i=0; i<sz; i++) {
01145 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
01146 }
01147
01148 image->update();
01149 with->update();
01150 image->do_fft_inplace();
01151 with->do_fft_inplace();
01152 image->set_attr("free_me",1);
01153 with->set_attr("free_me",1);
01154 }
01155
01156
01157 if (!image->is_complex()) {
01158 image=image->do_fft();
01159 image->set_attr("free_me",1);
01160 }
01161 if (!with->is_complex()) {
01162 with=with->do_fft();
01163 with->set_attr("free_me",1);
01164 }
01165
01166 fsc = image->calc_fourier_shell_correlation(with,1);
01167 }
01168
01169 int ny = image->get_ysize();
01170 int ny2=ny/2+1;
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205 vector<float> snr;
01206 if (snrweight) {
01207 Ctf *ctf = NULL;
01208 if (!image->has_attr("ctf")) {
01209 if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
01210 ctf=with->get_attr("ctf");
01211 }
01212 else ctf=image->get_attr("ctf");
01213
01214 float ds=1.0f/(ctf->apix*ny);
01215 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR);
01216 if(ctf) {delete ctf; ctf=0;}
01217 }
01218
01219 vector<float> amp;
01220 if (ampweight) amp=image->calc_radial_dist(ny/2,0,1,0);
01221
01222
01223 float pmin,pmax;
01224 if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;
01225 else pmin=0;
01226 if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
01227 else pmax=0;
01228
01229 double sum=0.0, norm=0.0;
01230
01231 for (int i=0; i<ny/2; i++) {
01232 double weight=1.0;
01233 if (sweight) weight*=fsc[(ny2)*2+i];
01234 if (ampweight) weight*=amp[i];
01235 if (snrweight) weight*=snr[i];
01236
01237
01238
01239
01240
01241 if (pmin>0) weight*=(tanh(5.0*(i-pmin)/pmin)+1.0)/2.0;
01242 if (pmax>0) weight*=(1.0-tanh(i-pmax))/2.0;
01243
01244 sum+=weight*fsc[ny2+i];
01245 norm+=weight;
01246
01247 }
01248
01249
01250 sum/=norm;
01251 if (nweight && with->get_attr_default("ptcl_repr",0) && sum>=0 && sum<1.0) {
01252 sum=sum/(1.0-sum);
01253 sum/=(float)with->get_attr_default("ptcl_repr",0);
01254 sum=sum/(1.0+sum);
01255 }
01256
01257 if (image->has_attr("free_me")) delete image;
01258 if (with->has_attr("free_me")) delete with;
01259
01260 EXITFUNC;
01261
01262
01263 if (!Util::goodf(&sum)) sum=-3.0e38;
01264
01265
01266
01267 return (float)-sum;
01268 }
01269
01270 void EMAN::dump_cmps()
01271 {
01272 dump_factory < Cmp > ();
01273 }
01274
01275 map<string, vector<string> > EMAN::dump_cmps_list()
01276 {
01277 return dump_factory_list < Cmp > ();
01278 }
01279
01280