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