#include <cmp.h>
Inheritance diagram for EMAN::FRCCmp:
Public Member Functions | |
float | cmp (EMData *image, EMData *with) const |
To compare 'image' with another image passed in through its parameters. | |
string | get_name () const |
Get the Cmp's name. | |
string | get_desc () const |
TypeDict | get_param_types () const |
Get Cmp parameter information in a dictionary. | |
Static Public Member Functions | |
static Cmp * | NEW () |
Static Public Attributes | |
static const string | NAME = "frc" |
Fourier ring correlation (FRC) is a measure of statistical dependency between two averages, computed by comparison of rings in Fourier space. 1 means prefect agreement. 0 means no correlation.
Definition at line 552 of file cmp.h.
To compare 'image' with another image passed in through its parameters.
An optional transformation may be used to transform the 2 images.
image | The first image to be compared. | |
with | The second image to be comppared. |
Implements EMAN::Cmp.
Definition at line 1248 of file cmp.cpp.
References EMAN::Ctf::apix, EMAN::EMData::calc_fourier_shell_correlation(), EMAN::EMData::calc_radial_dist(), EMAN::Ctf::compute_1d(), EMAN::EMData::copy(), EMAN::Ctf::CTF_SNR, EMAN::EMData::do_fft(), EMAN::EMData::do_fft_inplace(), ENTERFUNC, EXITFUNC, EMAN::EMObject::f, EMAN::EMData::get_attr(), EMAN::EMData::get_attr_default(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Util::goodf(), EMAN::EMData::has_attr(), InvalidCallException, EMAN::EMData::is_complex(), norm(), ny, EMAN::Cmp::params, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::update(), EMAN::Cmp::validate_input_args(), and weight.
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 // The fast hypot here was supposed to speed things up. Little effect 01302 // if (image->get_zsize()>1) fsc = image->calc_fourier_shell_correlation(with,1); 01303 // else { 01304 // double *sxy = (double *)malloc(ny2*sizeof(double)*4); 01305 // double *sxx = sxy+ny2; 01306 // double *syy = sxy+2*ny2; 01307 // double *norm= sxy+3*ny2; 01308 // 01309 // float *df1=image->get_data(); 01310 // float *df2=with->get_data(); 01311 // int nx2=image->get_xsize(); 01312 // 01313 // for (int y=-ny/2; y<ny/2; y++) { 01314 // for (int x=0; x<nx2/2; x++) { 01315 // if (x==0 && y<0) continue; // skip Friedel pair 01316 // short r=Util::hypot_fast_int(x,y); 01317 // if (r>ny2-1) continue; 01318 // int l=x*2+(y<0?ny+y:y)*nx2; 01319 // sxy[r]+=df1[l]*df2[l]+df1[l+1]*df2[l+1]; 01320 // sxx[r]+=df1[l]*df1[l]; 01321 // syy[r]+=df2[l]*df2[l]; 01322 // norm[r]+=1.0; 01323 // } 01324 // } 01325 // fsc.resize(ny2*3); 01326 // for (int r=0; r<ny2; r++) { 01327 // fsc[r]=r*0.5/ny2; 01328 // fsc[ny2+r]=sxy[r]/(sqrt(sxx[r])*sqrt(syy[r])); 01329 // fsc[ny2*2+r]=norm[r]; 01330 // } 01331 // free(sxy); 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; // make sure that points don't get completely excluded due to SNR estimation issues, or worse, contribute with a negative weight 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 // Min/max modifications to weighting 01355 float pmin,pmax; 01356 if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres; //cutoff in pixels, assume square 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 // if (snrweight) { 01369 // if (snr[i]>0) weight*=sqrt(snr[i]); 01370 // else weight=0; 01371 // } 01372 //if(snr[i]<0) printf("snr[%d] = %1.5g\n",i,snr[i]); 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 // printf("%d\t%f\t%f\n",i,weight,fsc[ny/2+1+i]); 01379 } 01380 01381 // This performs a weighting that tries to normalize FRC by correcting from the number of particles represented by the average 01382 sum/=norm; 01383 if (nweight && with->get_attr_default("ptcl_repr",0) && sum>=0 && sum<1.0) { 01384 sum=sum/(1.0-sum); // convert to SNR 01385 sum/=(float)with->get_attr_default("ptcl_repr",0); // divide by ptcl represented 01386 sum=sum/(1.0+sum); // convert back to correlation 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 //.Note the negative! This is because EMAN2 follows the convention that 01397 // smaller return values from comparitors indicate higher similarity - 01398 // this enables comparitors to be used in a generic fashion. 01399 return (float)-sum; 01400 }
string EMAN::FRCCmp::get_desc | ( | ) | const [inline, virtual] |
string EMAN::FRCCmp::get_name | ( | ) | const [inline, virtual] |
TypeDict EMAN::FRCCmp::get_param_types | ( | ) | const [inline, virtual] |
Get Cmp parameter information in a dictionary.
Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.
Implements EMAN::Cmp.
Definition at line 572 of file cmp.h.
References EMAN::EMObject::FLOAT, EMAN::EMObject::INT, and EMAN::TypeDict::put().
00573 { 00574 TypeDict d; 00575 d.put("snrweight", EMObject::INT, "If set, the SNR of 'this' will be used to weight the result. If 'this' lacks CTF info, it will check 'with'. (default=0)"); 00576 d.put("ampweight", EMObject::INT, "If set, the amplitude of 'this' will be used to weight the result (default=0)"); 00577 d.put("sweight", EMObject::INT, "If set, weight the (1-D) average by the number of pixels in each ring (default=1)"); 00578 d.put("nweight", EMObject::INT, "Downweight similarity based on number of particles in reference (default=0)"); 00579 d.put("zeromask", EMObject::INT, "Treat regions in either image that are zero as a mask"); 00580 d.put("minres", EMObject::FLOAT, "Lowest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=500"); 00581 d.put("maxres", EMObject::FLOAT, "Highest resolution to use in comparison (soft cutoff). Requires accurate A/pix in image. <0 disables. Default=10"); 00582 return d; 00583 }
static Cmp* EMAN::FRCCmp::NEW | ( | ) | [inline, static] |
const string FRCCmp::NAME = "frc" [static] |