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

Segment a volume into ~n subvolumes using K-means classification. More...

#include <processor.h>

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

List of all members.

Public Member Functions

string get_name () const
 Get the processor's name.
virtual EMDataprocess (const EMData *const image)
 To proccess an image out-of-place.
void process_inplace (EMData *image)
 To process an image in-place.
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 = "segment.kmeans"

Detailed Description

Segment a volume into ~n subvolumes using K-means classification.

Author:
Steve Ludtke
Date:
2008/11/03
Parameters:
ctf[in]A Ctf object to use

Definition at line 788 of file processor.h.


Member Function Documentation

string EMAN::KmeansSegmentProcessor::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 825 of file processor.h.

                {
                        return "Performs K-means segmentation on a volume. Note that this method uses random seeds, and thus will return different results each time it is run. Returned map contains number of segment for each voxel (or 0 for unsegmented voxels). Segmentation centers are stored in 'segmentcenters' attribute, consisting of a list of 3n floats in x,y,z triples.";
                }
string EMAN::KmeansSegmentProcessor::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 791 of file processor.h.

References NAME.

                {
                        return NAME;
                }
TypeDict EMAN::KmeansSegmentProcessor::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.

An option for psudoatom generation in pathwalker. Instead of random seeding, seed on the gird initially.

Author:
Muyuan Chen
Date:
2014/06/05

Reimplemented from EMAN::Processor.

Definition at line 799 of file processor.h.

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

                {
                        TypeDict d ;
                        d.put("nseg", EMObject::INT, "Number of segments to divide the image into. default=12" );
                        d.put("thr",EMObject::FLOAT,"Isosurface threshold value. Pixels below this will not be segmented");
                        d.put("ampweight",EMObject::INT,"If set, will weight centers by voxel amplitude. default = 1");
                        d.put("maxsegsize",EMObject::FLOAT,"Maximum radial distance from segment center to member voxel. Default=10000");
                        d.put("minsegsep",EMObject::FLOAT,"Minimum segment separation. Segments too close will trigger a reseed");
                        d.put("maxiter",EMObject::FLOAT,"Maximum number of iterations to run before stopping. Default=100");
                        d.put("maxvoxmove",EMObject::FLOAT,"Maximum number of voxels that can move before quitting. Default=25");
                        d.put("verbose",EMObject::INT,"Be verbose while running");
                        d.put("psudoatom",EMObject::BOOL,"Doing psudoatom generation");
                        d.put("sep",EMObject::FLOAT,"Separation distance, used only in psudoatom generation. Default=3.78");
                        return d;
                }
static Processor* EMAN::KmeansSegmentProcessor::NEW ( ) [inline, static]

Definition at line 820 of file processor.h.

                {
                        return new KmeansSegmentProcessor();
                }
EMData * KmeansSegmentProcessor::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 1010 of file processor.cpp.

References abs, EMAN::EMData::copy(), EMAN::EMData::get_attr(), EMAN::Util::get_frand(), EMAN::EMData::get_value_at(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Util::hypot3(), nx, ny, EMAN::Processor::params, q, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), EMAN::EMData::set_value_at(), v, x, and y.

{
        EMData * result = image->copy();

        int nseg = params.set_default("nseg",12);       // Muyuan, you cannot make major changes to default values without thorough testing! I returned this to 12
        float thr = params.set_default("thr",-1.0e30f);
        int ampweight = params.set_default("ampweight",1);
        float maxsegsize = params.set_default("maxsegsize",10000.0f);
        float minsegsep = params.set_default("minsegsep",0.0f);
        int maxiter = params.set_default("maxiter",100);
        int maxvoxmove = params.set_default("maxvoxmove",25);
        int verbose = params.set_default("verbose",0);
        bool psudoatom = params.set_default("psudoatom",0);
        float sep = params.set_default("sep",3.78f);

        int nx=image->get_xsize();
        int ny=image->get_ysize();
        int nz=image->get_zsize();
//      int nxy=nx*ny;

        // seed
        vector<float> centers(nseg*3);
        vector<float> count(nseg);
        // Alternative seeding method for paudoatom generation. Seed on the gird.
        if (psudoatom){
                float ax=image->get_attr("apix_x");
                sep/=ax;
                if (verbose) printf("Seeding .....\n");
                int sx=int(nx/sep)+1,sy=int(ny/sep)+1,sz=int(nz/sep)+1;
                EMData m(sx,sy,sz);
                EMData mcount(sx,sy,sz);

                for (int i=0; i<nx; i++){
                        for (int j=0; j<ny; j++){
                                for (int k=0; k<nz; k++){
                                        int ni=(i/sep),nj=(j/sep),nk=(k/sep);
                                        float v=image->get_value_at(i,j,k);
                                        if (v>thr){
                                                m.set_value_at(ni,nj,nk,(m.get_value_at(ni,nj,nk)+v));
                                                mcount.set_value_at(ni,nj,nk,(mcount.get_value_at(ni,nj,nk)+1));
                                        }
                                }
                        }
                }
                int nsum=0;
                float l=0,r=2000,th=5;
                while (abs(nsum-nseg)>0){
                        th=(l+r)/2;
                        nsum=0;
                        for (int i=0; i<sx; i++){
                                for (int j=0; j<sy; j++){
                                        for (int k=0; k<sz; k++){
                                                if (m.get_value_at(i,j,k)>th)  nsum+=1;
                                        }
                                }
                        }
                        if (verbose) printf("%3f\t %3f\t %3f,\t %4d\t %4d\n", l,th,r,nsum,nseg);
                        if (nsum>nseg) l=th;
                        if (nsum<nseg) r=th;
                        if ((r-l)<.01) break;
                }
//              nseg=nsum;
                int q=0;
                for (int i=0; i<sx; i++){
                        for (int j=0; j<sy; j++){
                                for (int k=0; k<sz; k++){
                                        if (m.get_value_at(i,j,k)>th){
                                                if(q<nseg*3){
                                                        centers[q]=  float(i+.5)*sep;
                                                        centers[q+1]=float(j+.5)*sep;
                                                        centers[q+2]=float(k+.5)*sep;
                                                        q+=3;
                                                }
                                        }
                                }
                        }
                }
        }
        // Default: random seeding.
        else{
                for (int i=0; i<nseg*3; i+=3) {
                        centers[i]=  Util::get_frand(0.0f,(float)nx);
                        centers[i+1]=Util::get_frand(0.0f,(float)ny);
                        centers[i+2]=Util::get_frand(0.0f,(float)nz);
                }
        }

        for (int iter=0; iter<maxiter; iter++) {
                // **** classify
                size_t pixmov=0;                // count of moved pixels
                for (int z=0; z<nz; z++) {
                        for (int y=0; y<ny; y++) {
                                for (int x=0; x<nx; 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 (int c=0; c<nseg; 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; }
                                        }
                                        if ((int)result->get_value_at(x,y,z)!=bcls) pixmov++;
                                        if (bdist>maxsegsize) result->set_value_at(x,y,z,-1);           // pixel is too far from any center
                                        else result->set_value_at(x,y,z,(float)bcls);           // set the pixel to the class number
                                }
                        }
                }

                // **** adjust centers
                for (int i=0; i<nseg*3; i++) centers[i]=0;
                for (int i=0; i<nseg; i++) count[i]=0;
                // weighted sums
                for (int z=0; z<nz; z++) {
                        for (int y=0; y<ny; y++) {
                                for (int x=0; x<nx; x++) {
                                        int cls = (int)result->get_value_at(x,y,z);
                                        if (cls==-1) continue;
                                        float w=1.0;
                                        if (ampweight) w=image->get_value_at(x,y,z);

                                        centers[cls*3]+=x*w;
                                        centers[cls*3+1]+=y*w;
                                        centers[cls*3+2]+=z*w;
                                        count[cls]+=w;
                                }
                        }
                }
                // now each becomes center of mass, or gets randomly reseeded
                int nreseed=0;
                for (int c=0; c<nseg; c++) {
                        // reseed
                        if (count[c]==0) {
                                nreseed++;
                                do {
                                        centers[c*3]=  Util::get_frand(0.0f,(float)nx);
                                        centers[c*3+1]=Util::get_frand(0.0f,(float)ny);
                                        centers[c*3+2]=Util::get_frand(0.0f,(float)nz);
                                } while (image->get_value_at((int)centers[c*3],(int)centers[c*3+1],(int)centers[c*3+2])<thr);           // This makes sure the new point is inside density
                        }
                        // center of mass
                        else {
                                centers[c*3]/=count[c];
                                centers[c*3+1]/=count[c];
                                centers[c*3+2]/=count[c];
                        }
                }

                // with minsegsep, check separation
                if (minsegsep>0) {
                        for (int c1=0; c1<nseg-1; c1++) {
                                for (int c2=c1+1; c2<nseg; c2++) {
                                        if (Util::hypot3(centers[c1*3]-centers[c2*3],centers[c1*3+1]-centers[c2*3+1],centers[c1*3+2]-centers[c2*3+2])<=minsegsep) {
                                                nreseed++;
                                                do {
                                                        centers[c1*3]=  Util::get_frand(0.0f,(float)nx);
                                                        centers[c1*3+1]=Util::get_frand(0.0f,(float)ny);
                                                        centers[c1*3+2]=Util::get_frand(0.0f,(float)nz);
                                                } while (image->get_value_at((int)centers[c1*3],(int)centers[c1*3+1],(int)centers[c1*3+2])<thr);
                                        }
                                }
                        }
                }


                if (verbose) printf("Iteration %3d: %6ld voxels moved, %3d classes reseeded\n",iter,pixmov,nreseed);
                if (nreseed==0 && pixmov<(size_t)maxvoxmove) break;             // termination conditions met
        }

        result->set_attr("segment_centers",centers);

        return result;
}
void KmeansSegmentProcessor::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 1185 of file processor.cpp.

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

Member Data Documentation

const string KmeansSegmentProcessor::NAME = "segment.kmeans" [static]

Definition at line 830 of file processor.h.

Referenced by get_name().


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