'paints' a circle into the image at x,y,z with values inside r1 set to v1, values between r1 and r2 will be set to a value between v1 and v2, and values outside r2 will be unchanged More...

#include <processor.h>

Inheritance diagram for EMAN::WatershedProcessor:

[legend]

Collaboration diagram for EMAN::WatershedProcessor:

[legend]

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

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

virtual string	get_name () const
	Get the processor's name. More...

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

virtual TypeDict	get_param_types () const
	Get processor parameter information in a dictionary. 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.watershed"

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

'paints' a circle into the image at x,y,z with values inside r1 set to v1, values between r1 and r2 will be set to a value between v1 and v2, and values outside r2 will be unchanged

Parameters

xpoints	x coordinates
ypoints	y coordinates
zpoints	z coordinates
minval	min value

Definition at line 7695 of file processor.h.

Member Function Documentation

◆ get_desc()

virtual string EMAN::WatershedProcessor::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 7711 of file processor.h.

                {
                        return "Watershed segmentation. Warning: uses up to 2.5x the map size in RAM. This will segment all voxels above threshold except for a 1-voxel wide border on all edges.";
                }

◆ get_name()

virtual string EMAN::WatershedProcessor::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 7701 of file processor.h.

                {
                        return NAME;
                }

References NAME.

◆ get_param_types()

virtual TypeDict EMAN::WatershedProcessor::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 7716 of file processor.h.

                {
                        TypeDict d;
                        d.put("nseg", EMObject::INT, "Number of segments to (attempt) to produce. The actual number may be fewer. (default=12)" );
                        d.put("thr",EMObject::FLOAT,"Isosurface threshold value. Pixels below this value will not be segmented. All voxels above this value will be segmented. (default=0.5)");
                        d.put("segbymerge", EMObject::INT, "If set, will achieve the specified number of segments by progressively merging the most connected segments. Can produce very different results." );
                        d.put("verbose", EMObject::INT, "If set, will print console output while running" );
                        return d;
                }

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

◆ NEW()

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

inlinestatic

Definition at line 7706 of file processor.h.

                {
                        return new WatershedProcessor();
                }

◆ process()

EMData * WatershedProcessor::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 11885 of file processor.cpp.

                                                             {
        //unsigned int nseg = params.set_default("nseg",12);
        int nseg = params.set_default("nseg",12);
        float thr = params.set_default("thr",0.5f);
        int segbymerge = params.set_default("segbymerge",0);
        vector<float> centers;
        int verbose = params.set_default("verbose",0);
        if (nseg<=1) throw InvalidValueException(nseg,"nseg must be greater than 1");
 
        if (segbymerge) { segbymerge=nseg; nseg=4096; }         // set max number of segments to a large (but not infinite) value. Too many will make connectivity matrix too big
 
        EMData *ret=new EMData(image->get_xsize(),image->get_ysize(),image->get_zsize());
        ret->to_zero();
 
 
        int nx=image->get_xsize();
        int ny=image->get_ysize();
        int nz=image->get_zsize();
        if (nz==1) throw ImageDimensionException("Only 3-D data supported");
 
        // Count the number of above threshold pixels
        size_t n2seg = 0;
        for (int z=1; z<nz-1; z++) {
                for (int y=1; y<ny-1; y++) {
                        for (int x=1; x<nx-1; x++) {
                                if (image->get_value_at(x,y,z)>=thr) n2seg++;
                        }
                }
        }
        if (verbose) printf("%ld voxels above threshold\n",n2seg);
 
        // Extract the pixels for sorting
        vector<WSsortlist> srt(n2seg);
        size_t i=0;
        for (int z=1; z<nz-1; z++) {
                for (int y=1; y<ny-1; y++) {
                        for (int x=1; x<nx-1; x++) {
                                if (image->get_value_at(x,y,z)>=thr) {
                                        srt[i].pix=image->get_value_at(x,y,z);
                                        srt[i].x=x;
                                        srt[i].y=y;
                                        srt[i].z=z;
                                        i++;
                                }
                        }
                }
        }
        if (verbose) printf("Voxels extracted, sorting\n");
 
        // actual sort
        sort(srt.begin(), srt.end());
        if (verbose) printf("Voxels sorted (%1.4g max), starting watershed\n",srt[0].pix);
 
        // now we start with the highest value and fill in the segments
        float cseg=1.0;
        int start=n2seg;
        for (i=0; i<n2seg; i++) {
                int x=srt[i].x;
                int y=srt[i].y;
                int z=srt[i].z;
                float lvl=0;
                for (int zz=z-1; zz<=z+1; zz++) {
                        for (int yy=y-1; yy<=y+1; yy++) {
                                for (int xx=x-1; xx<=x+1; xx++) {
                                        float v=ret->get_value_at(xx,yy,zz);
                                        if (v>lvl) lvl=v;                               // use the highest numbered border segment (arbitrary)
                                }
                        }
                }
                if (lvl==0) {
                        if (verbose) printf("%d %d %d\t%1.0f\t%1.3g\n",x,y,z,cseg,srt[i].pix);
                        lvl=cseg;
                        centers.push_back(x);
                        centers.push_back(y);
                        centers.push_back(z);
                        cseg+=1.0;
                }
                if (lvl>nseg) {
                        start=i;
                        if (verbose) printf("Requested number of segments achieved at density %1.4g\n",srt[i].pix);
                        break;
                }               // This means we've made as many segments as we need, so we switch to flood-filling
                ret->set_value_at_fast(x,y,z,lvl);
        }
 
        // We have as many segments as we'll get, but not all voxels have been segmented, so we do a progressive flood fill in density order
        size_t chg=1;
        while (chg) {
                chg=0;
                for (i=start; i<n2seg; i++) {
                        int x=srt[i].x;
                        int y=srt[i].y;
                        int z=srt[i].z;
                        if (ret->get_value_at(x,y,z)!=0) continue;      // This voxel is already done
 
                        float lvl=0;
                        for (int zz=z-1; zz<=z+1; zz++) {
                                for (int yy=y-1; yy<=y+1; yy++) {
                                        for (int xx=x-1; xx<=x+1; xx++) {
                                        float v=ret->get_value_at(xx,yy,zz);
                                        if (v>lvl) lvl=v;                               // use the highest numbered border segment (arbitrary)
                                        }
                                }
                        }
                        if (lvl==0) continue;                                   // we just skip voxels without any segmented neighbors
                        ret->set_value_at_fast(x,y,z,lvl);
                        chg+=1;
                }
                if (verbose) printf("%ld voxels changed\n",chg);
        }
        ret->set_attr("segment_centers",centers);
 
        if (segbymerge) {
                if (cseg<segbymerge) return ret;
        }
        else if (cseg<=nseg) return ret;                // We don't have too many segments, so we just return now
 
        if (verbose) printf("Merging segments\n");
        // If requested, we now merge segments with the most surface contact until we have the correct final number
        if (segbymerge) {
                int nsegstart=(int)cseg;        // Number of segments we actually generated
                nseg=(int)cseg;
                EMData *mx=new EMData(nsegstart,nsegstart,1);           // This will be a "contact matrix" among segments
                float *mxd=mx->get_data();
 
                // each cycle of the while loop, we eliminate one segment by merging
                int sub1=-1,sub2=-1;            // sub2 will be merged into sub1
                nseg++;                                         // since we don't actually remove one on the first pass, but decrement the counter
                while (segbymerge<nseg) {
                        mx->to_zero();
 
                        for (i=0; i<n2seg; i++) {
                                int x=srt[i].x;
                                int y=srt[i].y;
                                int z=srt[i].z;
 
                                int v1=(int)ret->get_value_at(x,y,z);
                                if (v1==sub2) { ret->set_value_at_fast(x,y,z,sub1); v1=sub1; }
                                mxd[v1+v1*nsegstart]++;                                 // the diagonal is a count of the number of voxels in the segment
                                for (int zz=z-1; zz<=z+1; zz++) {
                                        for (int yy=y-1; yy<=y+1; yy++) {
                                                for (int xx=x-1; xx<=x+1; xx++) {
                                                        int v2=(int)ret->get_value_at(xx,yy,zz);
                                                        if (v2==sub2) v2=sub1;          // pretend that any sub2 values are actually sub1
                                                        if (v1==v2) continue;
                                                        mxd[v1+v2*nsegstart]+=image->get_value_at(xx,yy,zz);            // We weight the connectivity by the image value
                                                }
                                        }
                                }
                        }
                        mx->update();
                        nseg--;                                 // number of segments left
                        if (verbose && sub1==-1) { mx->write_image("contactmx.hdf",0); }                // for debugging
 
                        sub1=-1;
                        sub2=-1;
                        // contact matrix complete, now figure out which 2 segments to merge
                        // diagonal of matrix is a count of the 'volume' of the segment. off-diagonal elements are surface area of contact region (roughly)
                        // we want to normalize the surface area elements so they are roughly proportional to the size of the segment, so we don't merge
                        // based on total contact area, but contact area as a fraction of the total area.
                        float bestv=-1.0;
                        for (int s1=1; s1<nsegstart; s1++) {
                                for (int s2=1; s2<nsegstart; s2++) {
                                        if (s1==s2) continue;                           // ignore the diagonal
                                        float v=mxd[s1+s2*nsegstart];
                                        if (v==0) continue;                                     // empty segment
//                                      v/=(pow(mxd[s1+s1*nsegstart],0.6667f)+pow(mxd[s2+s2*nsegstart],0.6667f));       // normalize by the sum of the estimated surface areas (no shape effects)
                                        v/=max(mxd[s1+s1*nsegstart],mxd[s2+s2*nsegstart]);      // normalize by the sum of the estimated surface areas (no shape effects)
                                        if (v>bestv) { bestv=v; sub1=s1; sub2=s2; }
                                }
                        }
                        float mv=0;
                        int mvl=0;
                        for (i=nsegstart+1; i<nsegstart*nsegstart; i+=nsegstart+1)
                                if (mxd[i]>mv) { mv=mxd[i]; mvl=i/nsegstart; }
                        if (verbose) printf("Merging %d to %d (%1.0f, %d)\n",sub2,sub1,mv,mvl);
                        if (sub1==-1) {
                                if (verbose) printf("Unable to find segments to merge, aborting\n");
                                break;
                        }
                }
 
        }
 
        return ret;
}

References ImageDimensionException, InvalidValueException, EMAN::Processor::params, EMAN::Dict::set_default(), x, and y.

Referenced by process_inplace().

◆ process_inplace()

void WatershedProcessor::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 11879 of file processor.cpp.

                                                      {
        EMData *tmp=process(image);
        memcpy(image->get_data(),tmp->get_data(),image->get_size()*sizeof(float));
        delete tmp;
}

References process().

Member Data Documentation

◆ NAME

const string WatershedProcessor::NAME = "segment.watershed"

static

Definition at line 7726 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