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

Sorry for the pun. More...

#include <processor.h>

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

List of all members.

Public Member Functions

virtual void process_inplace (EMData *image)
 To process an image in-place.
virtual EMDataprocess (const EMData *const image)
 To proccess an image out-of-place.
string get_name () const
 Get the processor's name.
TypeDict get_param_types () const
 Get processor parameter information in a dictionary.
string get_desc () const
 Get the descrition of this specific processor.

Static Public Member Functions

static ProcessorNEW ()

Static Public Attributes

static const string NAME = "math.sub.optimal"

Detailed Description

Sorry for the pun.

This processor will take a second image and try to filter/scale it to optimally subtract it from the original image. The idea here is that if you have an image with noise plus a linear-filter modified projection, that a good measure of the similarity of the image to the projection would be to try and remove the projection from the image as optimally as possible, then compute the standard deviation of what's left.

Now you might say that if the total energy in the noisy image is normalized then this should be equivalent to just integrating the FSC, which is what we use to do the optimal subtraction in the first place. This would be true, but this "optimal subtraction" has other purposes as well, such as the e2extractsubparticles program.

Parameters:
refReference image to subtract
return_radialWill return the radial filter function applied to ref as filter_curve

Definition at line 4657 of file processor.h.


Member Function Documentation

string EMAN::SubtractOptProcessor::get_desc ( ) const [inline, virtual]

Get the descrition of this specific processor.

This function must be overwritten by a subclass.

Returns:
The description of this processor.

Implements EMAN::Processor.

Definition at line 4687 of file processor.h.

                {
                        return "This will filter/scale 'ref' optimally and subtract it from image using ring dot products in Fourier space for normalization. Cutoff frequencies apply a bandpass tophat filter to the output.";
                }
string EMAN::SubtractOptProcessor::get_name ( ) const [inline, virtual]

Get the processor's name.

Each processor is identified by a unique name.

Returns:
The processor's name.

Implements EMAN::Processor.

Definition at line 4663 of file processor.h.

References NAME.

                {
                        return NAME;
                }
TypeDict EMAN::SubtractOptProcessor::get_param_types ( ) const [inline, virtual]

Get processor 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.

Reimplemented from EMAN::Processor.

Definition at line 4673 of file processor.h.

References EMAN::EMObject::BOOL, EMAN::EMObject::EMDATA, EMAN::EMObject::FLOAT, and EMAN::TypeDict::put().

                {
                        TypeDict d;
                        d.put("ref", EMObject::EMDATA, "Reference image to subtract");
                        d.put("actual", EMObject::EMDATA, "If specified, ref is used for normalization, but actual is subtracted.");
                        d.put("low_cutoff_frequency", EMObject::FLOAT, "Absolute [0,0.5] low cut-off frequency.");
                        d.put("high_cutoff_frequency", EMObject::FLOAT, "Absolute [0,0.5] high cut-off frequency.");
                        d.put("ctfweight",EMObject::BOOL, "Filter the image by CTF before subtraction");
                        d.put("return_fft",EMObject::BOOL, "Skips the final IFT, and returns the FFT of the subtracted image");
                        d.put("return_radial", EMObject::BOOL, "Return the radial filter function as an attribute (filter_curve)");
                        d.put("return_presigma", EMObject::BOOL, "Return the sigma of the pre-subtracted image in real-space with the specified filter applied as sigma_presub. This is an expensive option.");
                        return d;
                }
static Processor* EMAN::SubtractOptProcessor::NEW ( ) [inline, static]

Definition at line 4668 of file processor.h.

                {
                        return new SubtractOptProcessor();
                }
EMData * SubtractOptProcessor::process ( const EMData *const  image) [virtual]

To proccess an image out-of-place.

For those processors which can only be processed out-of-place, override this function to give the right behavior.

Parameters:
imageThe image will be copied, actual process happen on copy of image.
Returns:
the image processing result, may or may not be the same size of the input image

Reimplemented from EMAN::Processor.

Definition at line 3971 of file processor.cpp.

References EMAN::Ctf::compute_2d_complex(), EMAN::EMData::copy(), EMAN::EMData::copy_head(), EMAN::Ctf::CTF_INTEN, EMAN::EMData::do_fft(), EMAN::EMData::do_ift(), EMAN::EMData::get_attr(), EMAN::EMData::get_complex_at(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::Util::hypot_fast(), InvalidCallException, EMAN::EMData::is_complex(), LOGWARN, EMAN::EMData::mult(), norm(), EMAN::Processor::params, EMAN::EMData::set_attr(), EMAN::EMData::set_complex_at(), EMAN::Dict::set_default(), sqrt(), and x.

{
        if (!image) {
                LOGWARN("NULL Image");
                return NULL;
        }

        EMData *refr = params["ref"];
        EMData *actual = params.set_default("actual",(EMData*)NULL);
        EMData *ref;
        bool return_radial = params.set_default("return_radial",false);
        bool return_fft = params.set_default("return_fft",false);
        bool ctfweight = params.set_default("ctfweight",false);
        bool return_presigma = params.set_default("return_presigma",false);
        int si0=(int)floor(params.set_default("low_cutoff_frequency",0.0f)*image->get_ysize());
        int si1=(int)ceil(params.set_default("high_cutoff_frequency",0.7071f)*image->get_ysize());              // include the corners unless explicitly excluded
        
        // We will be modifying imf, so it needs to be a copy
        EMData *imf;
        if (image->is_complex()) imf=image->copy();
        else imf=image->do_fft();

        if (ctfweight) {
                EMData *ctfi=imf->copy_head();
                Ctf *ctf;
//              if (image->has_attr("ctf")) 
                        ctf=(Ctf *)(image->get_attr("ctf"));
//              else ctf=(Ctf *)(ref->get_attr("ctf"));
                ctf->compute_2d_complex(ctfi,Ctf::CTF_INTEN);
                imf->mult(*ctfi);
                delete ctfi;
        }
        
        // Make sure ref is complex
        if (refr->is_complex()) ref=refr;
        else ref=refr->do_fft();

        EMData *actf;
        if (actual==NULL) actf=ref;
        else {
                if (ctfweight) throw InvalidCallException("math.sub.optimal: Sorry, cannot use ctfweight in combination with actual");
                if (actual->is_complex()) actf=actual;
                else actf=actual->do_fft();
        }
                
        int ny2=(int)(image->get_ysize()*sqrt(2.0)/2);
        vector <double>rad(ny2+1);
        vector <double>norm(ny2+1);

        // We are essentially computing an FSC here, but while the reference (the image
        // we plan to subtract) is normalized, the other image is not. This gives us a filter
        // to apply to the reference to optimally eliminate its contents from 'image'.
        for (int y=-ny2; y<ny2; y++) {
                for (int x=0; x<ny2; x++) {
                        int r=int(Util::hypot_fast(x,y));
                        if (r>ny2) continue;
                        std::complex<float> v1=imf->get_complex_at(x,y);
                        std::complex<float> v2=ref->get_complex_at(x,y);
                        rad[r]+=(double)(v1.real()*v2.real()+v1.imag()*v2.imag());
//                      norm[r]+=v2.real()*v2.real()+v2.imag()*v2.imag()+v1.real()*v1.real()+v1.imag()*v1.imag();
                        norm[r]+=(double)(v2.real()*v2.real()+v2.imag()*v2.imag());
                }
        }
        for (int i=0; i<ny2; i++) rad[i]/=norm[i];

//      FILE *out=fopen("dbug.txt","w");
//      for (int i=0; i<ny2; i++) fprintf(out,"%lf\t%lf\t%lf\n",(float)i,rad[i],norm[i]);
//      fclose(out);
        
        float oldsig=-1.0;
        // This option computes the real-space sigma on the input-image after the specified filter
        // This is an expensive option, but more efficient than computing the same using other means
        if (return_presigma) {
                for (int y=-ny2; y<ny2; y++) {
                        for (int x=0; x<imf->get_xsize()/2; x++) {
                                int r=int(Util::hypot_fast(x,y));
                                if (r>=ny2 || r>=si1 || r<si0) {
                                        imf->set_complex_at(x,y,0);
                                        continue;
                                }
                                std::complex<float> v1=imf->get_complex_at(x,y);
                                imf->set_complex_at(x,y,v1);
                        }
                }
                EMData *tmp=imf->do_ift();
                oldsig=(float)tmp->get_attr("sigma");
                delete tmp;
        }
        
        for (int y=-ny2; y<ny2; y++) {
                for (int x=0; x<imf->get_xsize()/2; x++) {
                        int r=int(Util::hypot_fast(x,y));
                        if (r>=ny2 || r>=si1 || r<si0) {
                                imf->set_complex_at(x,y,0);
                                continue;
                        }
                        std::complex<float> v1=imf->get_complex_at(x,y);
                        std::complex<float> v2=actf->get_complex_at(x,y);
                        v2*=(float)rad[r];
                        imf->set_complex_at(x,y,v1-v2);
                }
        }
        
        if (!refr->is_complex()) delete ref;
        if (actual!=NULL && !actual->is_complex()) delete actf;

        vector <float>radf;
        if (return_radial) {
                radf.resize(ny2);
                for (int i=0; i<ny2; i++) radf[i]=(float)rad[i];
        }
                
        if (!return_fft) {
                EMData *ret=imf->do_ift();
                delete imf;
                if (return_radial) ret->set_attr("filter_curve",radf);
                if (return_presigma) {
                        ret->set_attr("sigma_presub",oldsig);
//                      printf("Set %f\n",(float)image->get_attr("sigma_presub"));
                }
                return ret;
        }
        if (return_radial) imf->set_attr("filter_curve",radf);
        if (return_presigma) imf->set_attr("sigma_presub",oldsig);
        return imf;
}
void SubtractOptProcessor::process_inplace ( EMData image) [virtual]

To process an image in-place.

For those processors which can only be processed out-of-place, override this function to just print out some error message to remind user call the out-of-place version.

Parameters:
imageThe image to be processed.

Implements EMAN::Processor.

Definition at line 4098 of file processor.cpp.

References EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), LOGWARN, EMAN::Processor::process(), and EMAN::EMData::update().

{
        if (!image) {
                LOGWARN("NULL image");
                return;
        }

        EMData *tmp=process(image);
        memcpy(image->get_data(),tmp->get_data(),(size_t)image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));
        delete tmp;
        image->update();
        return;
}

Member Data Documentation

const string SubtractOptProcessor::NAME = "math.sub.optimal" [static]

Definition at line 4692 of file processor.h.

Referenced by get_name().


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