Segment a volume by sequentially finding the highest peak and subtracting a Gaussian at that point from the density after strongly filtering the map to a specified resolvability. More...

#include <processor.h>

Inheritance diagram for EMAN::GaussSegmentProcessor:

[legend]

Collaboration diagram for EMAN::GaussSegmentProcessor:

[legend]

Public Member Functions
string	get_name () const
	Get the processor's name. More...

virtual EMData *	process (const EMData *const image)
	To proccess an image out-of-place. More...

void	process_inplace (EMData *image)
	To process an image in-place. More...

TypeDict	get_param_types () const
	Get processor parameter information in a dictionary. More...

string	get_desc () const
	Get the descrition of this specific processor. More...

Public Member Functions inherited from EMAN::Processor
virtual	~Processor ()

virtual void	process_list_inplace (vector< EMData * > &images)
	To process multiple images using the same algorithm. More...

virtual Dict	get_params () const
	Get the processor parameters in a key/value dictionary. More...

virtual void	set_params (const Dict &new_params)
	Set the processor parameters using a key/value dictionary. More...

Static Public Member Functions
static Processor *	NEW ()

Static Public Member Functions inherited from EMAN::Processor
static string	get_group_desc ()
	Get the description of this group of processors. More...

static void	EMFourierFilterInPlace (EMData *fimage, Dict params)
	Compute a Fourier-filter processed image in place. More...

static EMData *	EMFourierFilter (EMData *fimage, Dict params)
	Compute a Fourier-processor processed image without altering the original image. More...

Static Public Attributes
static const string	NAME = "segment.gauss"

Additional Inherited Members
Public Types inherited from EMAN::Processor
enum	fourier_filter_types { TOP_HAT_LOW_PASS , TOP_HAT_HIGH_PASS , TOP_HAT_BAND_PASS , TOP_HOMOMORPHIC , GAUSS_LOW_PASS , GAUSS_HIGH_PASS , GAUSS_BAND_PASS , GAUSS_INVERSE , GAUSS_HOMOMORPHIC , BUTTERWORTH_LOW_PASS , BUTTERWORTH_HIGH_PASS , BUTTERWORTH_HOMOMORPHIC , KAISER_I0 , KAISER_SINH , KAISER_I0_INVERSE , KAISER_SINH_INVERSE , SHIFT , TANH_LOW_PASS , TANH_HIGH_PASS , TANH_HOMOMORPHIC , TANH_BAND_PASS , RADIAL_TABLE , CTF_ }
	Fourier filter Processor type enum. More...

Protected Attributes inherited from EMAN::Processor
Dict	params

Detailed Description

Segment a volume by sequentially finding the highest peak and subtracting a Gaussian at that point from the density after strongly filtering the map to a specified resolvability.

Author: Steve Ludtke

Date: 2021/04/18

Definition at line 1660 of file processor.h.

Member Function Documentation

◆ get_desc()

string EMAN::GaussSegmentProcessor::get_desc ( ) const

inlinevirtual

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 1688 of file processor.h.

                {
                        return "Segments a volume by sequentially finding and subtracting Gaussians at a specified resolvability.";
                }

◆ get_name()

string EMAN::GaussSegmentProcessor::get_name ( ) const

inlinevirtual

Get the processor's name.

Each processor is identified by a unique name.

Returns: The processor's name.

Implements EMAN::Processor.

Definition at line 1663 of file processor.h.

                {
                        return NAME;
                }

References NAME.

◆ get_param_types()

TypeDict EMAN::GaussSegmentProcessor::get_param_types ( ) const

inlinevirtual

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 1671 of file processor.h.

                {
                        TypeDict d ;
                        d.put("minratio",EMObject::FLOAT,"The ratio of the smallest amplitude segment to locate relative to the strongest peak (default=0.5)");
                        d.put("maxnseg",EMObject::INT,"Maximum number of segments to return (default = unlimited)");
                        d.put("width",EMObject::FLOAT,"Required: full width of Gaussians in A at 1/e (FWHM). Also used to determine map prefiltration.");
                        d.put("mask",EMObject::EMDATA,"Optional: mask to apply to map after filtration to limit where centers are placed");
                        d.put("skipseg",EMObject::INT,"Normally the returned map is a segmentation map, but this is unnecessary if only the center coordinates are needed. If 1, the returned map will be a residual volume, not a segmentation map. If 2, it will be the filtered map before segmentation");
                        d.put("verbose",EMObject::INT,"Be verbose while running");
                        return d;
                }

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

◆ NEW()

static Processor * EMAN::GaussSegmentProcessor::NEW ( )

inlinestatic

Definition at line 1683 of file processor.h.

                {
                        return new GaussSegmentProcessor();
                }

◆ process()

EMData * GaussSegmentProcessor::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

image The 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 1448 of file processor.cpp.

{
 
        float minratio = params.set_default("minratio",0.5f);
        int maxnseg = params.set_default("maxnseg",0);
        int skipseg = params.set_default("skipseg",0);
        float width = params.set_default("width",10.0f);
        int verbose = params.set_default("verbose",0);
        EMData *mask= (EMData *)params.set_default("mask",(EMData *)NULL);
        float apix=image->get_attr("apix_x");
        int nx=image->get_xsize();
        int ny=image->get_ysize();
        int nz=image->get_zsize();
 
        EMData * result = image->process("threshold.belowtozero",Dict("minval",0.0f));
        // The intent of these filters is to insure that the map effectively 
//      result->process_inplace("filter.lowpass.gauss",Dict("cutoff_freq",0.45/width)); // 0.45 = sqrt(2)/pi, resolvability threshold
//      result->process_inplace("filter.lowpass.gauss",Dict("cutoff_freq",0.637/width));        // 0.637 = 2/pi, (may be Rayleigh criterion?)
        if (mask!=NULL) {
                result->mult(*mask);
                result->process_inplace("normalize.local",Dict("radius",nx*apix/(3.0f*width),"threshold",(float)result->get_attr("sigma_nonzero")*1.2));
        }
        
        XYData gsf;
        gsf.make_gauss(nx*4,1.0f/apix,0.45/width);
//      gsf.make_gauss(nx*4,1.0f/apix,0.637/width);
        result->process_inplace("filter.setstrucfac",Dict("apix",apix,"strucfac",&gsf));
                
        // Generate a Gaussian we can subtract out of the volume quickly
        int gs=2*width/(float)image->get_attr("apix_x");
        EMData gsub(gs,gs,gs);
        gsub.to_one();
        gsub.process_inplace("mask.gaussian",Dict("outer_radius",int(width/(2.0*apix))));
 
        if (verbose>2) {
                result->set_attr("render_bits",12);
                result->set_attr("render_min",(float)result->get_attr("minimum"));
                result->set_attr("render_max",(float)result->get_attr("maximum"));
                result->write_image("segfilt.hdf",0,EMUtil::IMAGE_UNKNOWN,false,0,EMUtil::EM_COMPRESSED);
        }
        
        EMData *cache = NULL;
        if (skipseg==2) cache=result->copy();
        
        // original version of the code. A bit slow but works very well
//      vector<float> centers;
//      vector<float> amps;
//      if (maxnseg==0) maxnseg=2000000000;
//      float maxval=0.0f;
//      
//      while (amps.size()<maxnseg) {
//              IntPoint p = result->calc_max_location();
//              FloatPoint fp(p[0],p[1],p[2]);
//              
//              float amp=result->get_value_at(p[0],p[1],p[2]);
//              if (amp<maxval*minratio) break;
//              amps.push_back(amp);
//              centers.push_back(p[0]);
//              centers.push_back(p[1]);
//              centers.push_back(p[2]);
//              if (amp>maxval) maxval=amp;             // really this should only happen once...
//              
//              result->insert_scaled_sum(&gsub,fp,1.0,-amp);
//      }
        
        // This was an alternate version, to improve speed. While it is faster, the approximations its making
        // produce slightly less optimal results.
        vector<float> centers;
        vector<float> amps;
        if (maxnseg==0) maxnseg=2000000000;
        
        float thr=(float)result->get_attr("maximum")*minratio;
        while (amps.size()<maxnseg) {
                if ((float)result->get_attr("maximum")<=thr) break;
                // We find all peak values greater than our threshold
                vector<Pixel> pixels=result->calc_highest_locations(thr);
                if (pixels.size()==0) break;
                if (verbose>1) printf("%d %f %f %f %d\n",pixels.size(),pixels[0].get_value(),pixels[1].get_value(),pixels[2].get_value(),amps.size());
                
                for (int i=0; i<pixels.size(); i++) {
                        IntPoint p = pixels[i].get_point();
                        FloatPoint fp(p[0],p[1],p[2]);
                        
                        float amp=result->get_value_at(p[0],p[1],p[2]);
                        if (amp!=pixels[i].get_value()) continue;       // skip any values near a value we've just changed, should come back in the next round
                        //                      printf("%d %d %d     %f  %f   %f\n",p[0],p[1],p[2],amp,pixels[i].get_value(),pixels[i+1].get_value());
//                      if (i<pixels.size()-1 && amp<pixels[i+1].get_value()) continue; // if one of the identified peaks has been reduced by a nearby peak too much
                        if (amp<thr) continue;
//                      if (amp<thr) { maxnseg=amps.size(); break; }
                        amps.push_back(amp);
                        centers.push_back(p[0]);
                        centers.push_back(p[1]);
                        centers.push_back(p[2]);
                        
                        result->insert_scaled_sum(&gsub,fp,1.0,-amp);
                }
        }
 
        if (verbose) printf("Found %d centers\n",amps.size());
        
        if (verbose>2) {
                result->set_attr("render_bits",12);
                result->set_attr("render_min",(float)result->get_attr("minimum"));
                result->set_attr("render_max",(float)result->get_attr("maximum"));
                result->write_image("segfilt.hdf",1,EMUtil::IMAGE_UNKNOWN,false,0,EMUtil::EM_COMPRESSED);
        }
        if (!skipseg) {
        // after we have our list of centers classify pixels
        for (int z=0; z<nz; z++) {
                for (int y=0; y<ny; y++) {
                        for (int x=0; x<nz; x++) {
//                              if (image->get_value_at(x,y,z)<thr) {
//                                      result->set_value_at(x,y,z,-1.0);               //below threshold -> -1 (unclassified)
//                                      continue;
//                              }
                                int bcls=-1;                    // best matching class
                                float bdist=(float)(nx+ny+nz);  // distance for best class
                                for (unsigned int c=0; c<centers.size()/3; c++) {
                                        float d=Util::hypot3(x-centers[c*3],y-centers[c*3+1],z-centers[c*3+2]);
                                        if (d<bdist) { bdist=d; bcls=c; }
                                }
                                result->set_value_at(x,y,z,(float)bcls);                // set the pixel to the class number
                        }
                }
        }
        if (verbose) printf("segmented\n");
        }
 
        //if skipseg is set to 2, we return the filtered, unsegmented map (with valid centers)
        if (skipseg==2) {
                delete result;
                result=cache;
        }
        
        result->set_attr("segment_centers",centers);
        result->set_attr("segment_amps",amps);
        
        return result;
}

References EMAN::EMUtil::EM_COMPRESSED, EMAN::Util::hypot3(), EMAN::EMUtil::IMAGE_UNKNOWN, EMAN::XYData::make_gauss(), EMAN::Processor::params, EMAN::Dict::set_default(), x, and y.

◆ process_inplace()

void GaussSegmentProcessor::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

image The image to be processed.

Implements EMAN::Processor.

Definition at line 1441 of file processor.cpp.

{
        printf("Process inplace not implemented. Please use process.\n");
        return;
}

Member Data Documentation

◆ NAME

const string GaussSegmentProcessor::NAME = "segment.gauss"

static

Definition at line 1693 of file processor.h.

Referenced by get_name().

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

Public Member Functions

Static Public Member Functions

Static Public Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ get_desc()

◆ get_name()

◆ get_param_types()

◆ NEW()

◆ process()

◆ process_inplace()

Member Data Documentation

◆ NAME