EMAN2
Public Member Functions | Static Public Member Functions | Static Public Attributes
EMAN::TomoFscCmp Class Reference

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

static CmpNEW ()

Static Public Attributes

static 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:
sigmasThe number of times the standard deviation of Fourier amplitudes to accept
minresThe minimum resolution to accept
maxesThe maximum resolution to accept
apixThe 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:
imageThe first image to be compared.
withThe second image to be comppared.
Returns:
The comparison result. Smaller better by default

Implements EMAN::Cmp.

Definition at line 683 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().

{
        ENTERFUNC;
        bool usecpu = 1;
        bool del_imagefft = 0;
        bool del_withfft = 0;
        float score = 1.0f;
        
        //get parameters
        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");
        // BAD DESIGN!!!! The fact that I have to load attrs into a variable before I can do operations on them is silly
        
        //get threshold information
        float image_meanwedgeamp = image->get_attr("spt_wedge_mean");
        float image_sigmawedgeamp = image->get_attr("spt_wedge_sigma");
        
        // Use with amp stats if avaliable otherwise effectivly ignore
        float with_meanwedgeamp = image_meanwedgeamp;
        float with_sigmawedgeamp = image_sigmawedgeamp;
        if (with->has_attr("spt_wedge_mean") && with->has_attr("spt_wedge_sigma"))
        {
                with_meanwedgeamp = with->get_attr("spt_wedge_mean");
                with_sigmawedgeamp = with->get_attr("spt_wedge_sigma");
        }
        
        // Find threshold
        float sigmas = params.set_default("sigmas",5.0f);
        float img_amp_thres = pow(image_meanwedgeamp + sigmas*image_sigmawedgeamp, 2.0f);
        float with_amp_thres = pow(with_meanwedgeamp + sigmas*with_sigmawedgeamp, 2.0f);
        
        // take negative of score
        float negative = (float)params.set_default("negative", 1.0f);
        if (negative) negative=-1.0; else negative=1.0;
        //get apix, use param apix, if not specified use apix_x, if this is not specified then apix=1.0
        float apix = params.set_default("apix",image->get_attr_default("apix_x", 1.0f));
        //get min and max res
        float minres = params.set_default("minres",std::numeric_limits<float>::max());
        float maxres = params.set_default("maxres", 0.0f);
        
        //Check to ensure that images are complex
        EMData* image_fft = image;
        EMData* with_fft = with;
        
#ifdef EMAN2_USING_CUDA 
        //do CUDA FFT, does not allow minres, maxres yet
        if(EMData::usecuda == 1 && image->getcudarwdata() && with->getcudarwdata()) {
                if(!image->is_complex()){
                        del_imagefft = 1;
                        image_fft = image->do_fft_cuda();
                }
                if(!with->is_complex()){
                        del_withfft = 1;
                        with_fft = with->do_fft_cuda();
                }
                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());
                usecpu = 0;
        }
#endif  
        if(usecpu){
                if(!image->is_complex()){
                        del_imagefft = 1;
                        image_fft = image->do_fft();
                }
                if(!with->is_complex()){
                        del_withfft = 1;
                        with_fft = with->do_fft();
                }
                
                //loop over all voxels
                int count = 0;
                double sum_imgamp_sq = 0.0;
                double sum_withamp_sq = 0.0;
                double cong = 0.0;
                float* img_data = image_fft->get_data();
                float* with_data = with_fft->get_data();
        
                int nx  = image_fft->get_xsize();
                int ny  = image_fft->get_ysize();
                int nz  = image_fft->get_zsize();
                int ny2 = ny/2;
                int nz2 = nz/2;
        
                //compute FSC
                int ii, kz, ky;
                for (int iz = 0; iz < nz; iz++) {
                        if(iz > nz2) kz = nz-iz; else kz=iz;
                        for (int iy = 0; iy < ny; iy++) {
                                if(iy > ny2) ky = ny-iy; else ky=iy;
                                for (int ix = 0; ix < nx; ix+=2) {
                                        //compute spatial frequency and convert to resolution stat
                                        float freq = std::sqrt(kz*kz + ky*ky + ix*ix/4.0f)/float(nz);
                                        float reso = apix/freq;
                                        
                                        //only look within a resolution domain
                                        if(reso < minres && reso > maxres){
                                                ii = ix + (iy  + iz * ny)* nx;
                                                float img_r = img_data[ii];
                                                float img_i = img_data[ii+1];
                                                float with_r = with_data[ii];
                                                float with_i = with_data[ii+1];
                                                double img_amp_sq = img_r*img_r + img_i*img_i;
                                                double with_amp_sq = with_r*with_r + with_i*with_i;

                                                if((img_amp_sq >  img_amp_thres) && (with_amp_sq >  with_amp_thres)){
                                                        count ++;
                                                        sum_imgamp_sq += img_amp_sq;
                                                        sum_withamp_sq += with_amp_sq;
                                                        cong += img_r*with_r + img_i*with_i;
                                                }
                                        }
                                }
                        }
                }
                
                if(count > 0){ 
                        score = (float)(cong/sqrt(sum_imgamp_sq*sum_withamp_sq));
                }else{
                        score = 1.0f;
                }
        }
        
        //avoid mem leaks
        if(del_imagefft) delete image_fft;
        if(del_withfft) delete with_fft;
        
        return negative*score;
        EXITFUNC;
}
virtual string EMAN::TomoFscCmp::get_desc ( ) const [inline, virtual]

Implements EMAN::Cmp.

Definition at line 372 of file cmp.h.

                {
                        return "EXPERIMENTAL - Fsc with consideration given for the missing wedge. Not ready for routine use.";
                }
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.

References NAME.

                {
                        return NAME;
                }
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::EMObject::EMDATA, EMAN::EMObject::FLOAT, and EMAN::TypeDict::put().

                {
                        TypeDict d;
                        d.put("normalize", EMObject::EMDATA,"Return the negative value (which is EMAN2 convention), Defalut is true(1)");
                        d.put("sigmas", EMObject::FLOAT, "The number of times the standard deviation of Fourier amplitudes to accept");
                        d.put("minres", EMObject::FLOAT, "The minimum resolution to accept (1/A) Default is inf");
                        d.put("maxres", EMObject::FLOAT, "The maximum resolution to accept (1/A) Default=0.0");
                        d.put("apix", EMObject::FLOAT, "The angstroms per pixel to use. Default = apix_x(1.0 if not present)");
                        return d;
                }
static Cmp* EMAN::TomoFscCmp::NEW ( ) [inline, static]

Definition at line 377 of file cmp.h.

                {
                        return new TomoFscCmp();
                }

Member Data Documentation

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

Definition at line 393 of file cmp.h.

Referenced by get_name().


The documentation for this class was generated from the following files: