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