Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

EMAN::TomoFscCmp Class Reference
[a function or class that is CUDA enabled]

This is a FSC comparitor for tomography. More...

#include <cmp.h>

Inheritance diagram for EMAN::TomoFscCmp:

Inheritance graph
[legend]
Collaboration diagram for EMAN::TomoFscCmp:

Collaboration graph
[legend]
List of all members.

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

CmpNEW ()

Static Public Attributes

const string NAME = "fsc.tomo"

Detailed Description

This is a FSC comparitor for tomography.

I only counts voxels above a threshold, which is obtained from subtomogram metadata

Author:
John Flanagan
Date:
2012-2-02
Parameters:
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.


Member Function Documentation

float TomoFscCmp::cmp EMData image,
EMData with
const [virtual]
 

To compare 'image' with another image passed in through its parameters.

An optional transformation may be used to transform the 2 images.

Parameters:
image The first image to be compared.
with The second image to be comppared.
Returns:
The comparison result. Smaller better by default

Implements EMAN::Cmp.

Definition at line 681 of file cmp.cpp.

References EMAN::EMData::do_fft(), 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::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]
 

Implements EMAN::Cmp.

Definition at line 372 of file cmp.h.

00373                 {
00374                         return "Fsc with consideration given for the missing wedge";
00375                 }

virtual string EMAN::TomoFscCmp::get_name  )  const [inline, virtual]
 

Get the Cmp's name.

Each Cmp is identified by a unique name.

Returns:
The Cmp's name.

Implements EMAN::Cmp.

Definition at line 367 of file cmp.h.

00368                 {
00369                         return NAME;
00370                 }

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.

Returns:
A dictionary containing the parameter info.

Implements EMAN::Cmp.

Definition at line 382 of file cmp.h.

References 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                 }

Cmp* EMAN::TomoFscCmp::NEW  )  [inline, static]
 

Definition at line 377 of file cmp.h.

00378                 {
00379                         return new TomoFscCmp();
00380                 }


Member Data Documentation

const string TomoFscCmp::NAME = "fsc.tomo" [static]
 

Definition at line 52 of file cmp.cpp.


The documentation for this class was generated from the following files:
Generated on Tue Jun 11 13:41:56 2013 for EMAN2 by  doxygen 1.3.9.1