#include <cmp.h>
Inheritance diagram for EMAN::TomoFscCmp:
Public Member Functions | |
virtual float | cmp (EMData *image, EMData *with) const |
To compare 'image' with another image passed in through its parameters. | |
virtual string | get_name () const |
Get the Cmp's name. | |
virtual 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 = "fsc.tomo" |
I only counts voxels above a threshold, which is obtained from subtomogram metadata
sigmas | The number of times the standard deviation of Fourier amplitudes to accept | |
minres | The minimum resolution to accept | |
maxes | The maximum resolution to accept | |
apix | The angstroms per pixel to use |
Definition at line 362 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 681 of file cmp.cpp.
References EMAN::EMData::do_fft(), ENTERFUNC, EXITFUNC, fsc_tomo_cmp_cuda(), 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::EMData::has_attr(), InvalidCallException, EMAN::EMData::is_complex(), max, nx, ny, EMAN::Cmp::params, EMAN::Dict::set_default(), and sqrt().
00682 { 00683 ENTERFUNC; 00684 bool usecpu = 1; 00685 bool del_imagefft = 0; 00686 bool del_withfft = 0; 00687 float score = 1.0f; 00688 00689 //get parameters 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 // BAD DESIGN!!!! The fact that I have to load attrs into a variable before I can do operations on them is silly 00692 00693 //get threshold information 00694 float image_meanwedgeamp = image->get_attr("spt_wedge_mean"); 00695 float image_sigmawedgeamp = image->get_attr("spt_wedge_sigma"); 00696 00697 // Use with amp stats if avaliable otherwise effectivly ignore 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 // Find threshold 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 // take negative of score 00712 float negative = (float)params.set_default("negative", 1.0f); 00713 if (negative) negative=-1.0; else negative=1.0; 00714 //get apix, use param apix, if not specified use apix_x, if this is not specified then apix=1.0 00715 float apix = params.set_default("apix",image->get_attr_default("apix_x", 1.0f)); 00716 //get min and max res 00717 float minres = params.set_default("minres",std::numeric_limits<float>::max()); 00718 float maxres = params.set_default("maxres", 0.0f); 00719 00720 //Check to ensure that images are complex 00721 EMData* image_fft = image; 00722 EMData* with_fft = with; 00723 00724 #ifdef EMAN2_USING_CUDA 00725 //do CUDA FFT, does not allow minres, maxres yet 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 //loop over all voxels 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 //compute FSC 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 //compute spatial frequency and convert to resolution stat 00771 float freq = std::sqrt(kz*kz + ky*ky + ix*ix/4.0f)/float(nz); 00772 float reso = apix/freq; 00773 00774 //only look within a resolution domain 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 //avoid mem leaks 00803 if(del_imagefft) delete image_fft; 00804 if(del_withfft) delete with_fft; 00805 00806 return negative*score; 00807 EXITFUNC; 00808 }
virtual string EMAN::TomoFscCmp::get_desc | ( | ) | const [inline, virtual] |
virtual string EMAN::TomoFscCmp::get_name | ( | ) | const [inline, virtual] |
TypeDict EMAN::TomoFscCmp::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 382 of file cmp.h.
References EMAN::EMObject::EMDATA, EMAN::EMObject::FLOAT, and EMAN::TypeDict::put().
00383 { 00384 TypeDict d; 00385 d.put("normalize", EMObject::EMDATA,"Return the negative value (which is EMAN2 convention), Defalut is true(1)"); 00386 d.put("sigmas", EMObject::FLOAT, "The number of times the standard deviation of Fourier amplitudes to accept"); 00387 d.put("minres", EMObject::FLOAT, "The minimum resolution to accept (1/A) Default is inf"); 00388 d.put("maxres", EMObject::FLOAT, "The maximum resolution to accept (1/A) Default=0.0"); 00389 d.put("apix", EMObject::FLOAT, "The angstroms per pixel to use. Default = apix_x(1.0 if not present)"); 00390 return d; 00391 }
static Cmp* EMAN::TomoFscCmp::NEW | ( | ) | [inline, static] |
const string TomoFscCmp::NAME = "fsc.tomo" [static] |