EMAN2
Public Member Functions | Private Attributes
wustl_mm::SkeletonMaker::Volume Class Reference

#include <volume.h>

Collaboration diagram for wustl_mm::SkeletonMaker::Volume:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 Volume (EMData *em)
 Volume (int x, int y, int z)
 Volume (int x, int y, int z, float val)
 Volume (int x, int y, int z, int offx, int offy, int offz, Volume *vol)
 ~Volume ()
EMDataget_emdata ()
float getSpacingX ()
float getSpacingY ()
float getSpacingZ ()
float getOriginX ()
float getOriginY ()
float getOriginZ ()
int getSizeX ()
int getSizeY ()
int getSizeZ ()
int getIndex (int x, int y, int z)
double getDataAt (int x, int y, int z)
double getDataAt (int index)
void setSpacing (float spx, float spy, float spz)
void setOrigin (float orgX, float orgY, float orgZ)
void setDataAt (int x, int y, int z, double d)
void setDataAt (int index, double d)
void pad (int padBy, double padValue)
 * Next, smooth
int getNumNeighbor6 (int ox, int oy, int oz)
int hasCell (int ox, int oy, int oz)
VolumemarkCellFace ()
int hasCompleteSheet (int ox, int oy, int oz, Volume *fvol)
int hasCompleteSheet (int ox, int oy, int oz)
int hasCompleteHelix (int ox, int oy, int oz)
int hasCompleteHelix (int ox, int oy, int oz, Volume *fvol)
int isHelixEnd (int ox, int oy, int oz, Volume *nvol)
int isHelixEnd (int ox, int oy, int oz)
int isSheetEnd (int ox, int oy, int oz, Volume *nvol)
int isFeatureFace (int ox, int oy, int oz)
int isSheetEnd (int ox, int oy, int oz)
int isSimple (int ox, int oy, int oz)
int isPiercable (int ox, int oy, int oz)
int getNumPotComplex (int ox, int oy, int oz)
int getNumPotComplex2 (int ox, int oy, int oz)
int components6 (int vox[3][3][3])
 *
int components26 (int vox[3][3][3])
int countExt (double vox[3][3][3])
int countInt (double vox[3][3][3])
int countIntEuler (int ox, int oy, int oz)
void curveSkeleton (Volume *grayvol, float lowthr, float highthr, Volume *svol)
 * No topology check */
void curveSkeleton (float thr, Volume *svol)
void curveSkeleton2D (float thr, Volume *svol)
void skeleton (float thr, int off)
void skeleton (float thr, Volume *svol, Volume *hvol)
 * Added for debugging */ / Check simple
void erodeHelix ()
void erodeHelix (int disthr)
int erodeSheet ()
int erodeSheet (int disthr)
void surfaceSkeletonPres (float thr, Volume *preserve)
 * Deal with closed rings
void threshold (double thr)
 * Bertrand's parallel 6-connected surface thinning */
void threshold (double thr, int out, int in)
void threshold (double thr, int out, int in, int boundary)
void threshold (double thr, int out, int in, int boundary, bool markBoundary)
VolumeDatagetVolumeData ()

Private Attributes

VolumeDatavolData

Detailed Description

Definition at line 70 of file volume.h.


Constructor & Destructor Documentation

Volume::Volume ( EMData em)

Definition at line 9 of file volume.cpp.

References volData.

Referenced by curveSkeleton(), curveSkeleton2D(), erodeHelix(), erodeSheet(), markCellFace(), skeleton(), and surfaceSkeletonPres().

                {
                        this->volData = new VolumeData(em);
                }
Volume::Volume ( int  x,
int  y,
int  z 
)

Definition at line 13 of file volume.cpp.

References volData.

                                                  {
                        volData = new VolumeData(x, y, z);
                }
Volume::Volume ( int  x,
int  y,
int  z,
float  val 
)

Definition at line 17 of file volume.cpp.

References volData.

                                                             {
                        volData = new VolumeData(x, y, z, val);
                }
Volume::Volume ( int  x,
int  y,
int  z,
int  offx,
int  offy,
int  offz,
Volume vol 
)

Definition at line 21 of file volume.cpp.

References getVolumeData(), and volData.

                                                                                              {
                        volData = new VolumeData(x, y, z, offx, offy, offz, vol->getVolumeData());
                }
Volume::~Volume ( )

Definition at line 26 of file volume.cpp.

References volData.

                {
                        delete volData;
                }

Member Function Documentation

int Volume::components26 ( int  vox[3][3][3])

Definition at line 2767 of file volume.cpp.

References nx, ny, tot, x, and y.

Referenced by countExt().

                {
                        // Stupid flood fill
                        int tot = 0 ;
                        int queue[27][3] ;
                        int vis[3][3][3] ;
                        int head = 0, tail = 1 ;
                        int i, j, k ;
                        for ( i = 0 ; i < 3 ; i ++ )
                                for ( j = 0 ; j < 3 ; j ++ )
                                        for ( k = 0 ; k < 3 ; k ++ )
                                        {
                                                vis[i][j][k] = 0 ;
                                                if ( vox[i][j][k] )
                                                {
                                                        if ( tot == 0 )
                                                        {
                                                                queue[0][0] = i ;
                                                                queue[0][1] = j ;
                                                                queue[0][2] = k ;
                                                                vis[i][j][k] = 1 ;
                                                        }
                                                        tot ++ ;
                                                }
                                        }
                        if ( tot == 0 )
                        {
                                return 0 ;
                        }

                        int ct = 1 ;
                        while ( head != tail )
                        {
                                int x = queue[head][0] ;
                                int y = queue[head][1] ;
                                int z = queue[head][2] ;
                                head ++ ;

                                for ( i = -1 ; i < 2 ; i ++ )
                                        for ( j = -1 ; j < 2 ; j ++ )
                                                for ( k = -1 ; k < 2 ; k ++ )
                                                {
                                                        int nx = x + i ;
                                                        int ny = y + j ;
                                                        int nz = z + k ;
                                                        if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
                                                        {
                                                                if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
                                                                {
                                                                        queue[tail][0] = nx ;
                                                                        queue[tail][1] = ny ;
                                                                        queue[tail][2] = nz ;
                                                                        tail ++ ;
                                                                        vis[nx][ny][nz] = 1 ;
                                                                        ct ++ ;
                                                                }
                                                        }
                                                }
                        }

                        if ( ct == tot )
                        {
                                return 1 ;
                        }
                        else
                        {
                                return 2 ;
                        }

                }
int Volume::components6 ( int  vox[3][3][3])

*

Definition at line 2698 of file volume.cpp.

References wustl_mm::SkeletonMaker::neighbor6, nx, ny, tot, x, and y.

Referenced by countInt(), and countIntEuler().

                {
                        // Stupid flood fill
                        int tot = 0 ;
                        int queue[27][3] ;
                        int vis[3][3][3] ;
                        int head = 0, tail = 1 ;
                        int i, j, k ;
                        for ( i = 0 ; i < 3 ; i ++ )
                                for ( j = 0 ; j < 3 ; j ++ )
                                        for ( k = 0 ; k < 3 ; k ++ )
                                        {
                                                vis[i][j][k] = 0 ;
                                                if ( vox[i][j][k] )
                                                {
                                                        if ( tot == 0 )
                                                        {
                                                                queue[0][0] = i ;
                                                                queue[0][1] = j ;
                                                                queue[0][2] = k ;
                                                                vis[i][j][k] = 1 ;
                                                        }
                                                        tot ++ ;
                                                }
                                        }
                        if ( tot == 0 )
                        {
                                return 0 ;
                        }
                        // printf("total: %d\n", tot) ;

                        int ct = 1 ;
                        while ( head != tail )
                        {
                                int x = queue[head][0] ;
                                int y = queue[head][1] ;
                                int z = queue[head][2] ;
                                head ++ ;

                                for ( i = 0 ; i < 6 ; i ++ )
                                {
                                        int nx = x + neighbor6[i][0] ;
                                        int ny = y + neighbor6[i][1] ;
                                        int nz = z + neighbor6[i][2] ;
                                        if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
                                        {
                                                if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
                                                {
                                                        queue[tail][0] = nx ;
                                                        queue[tail][1] = ny ;
                                                        queue[tail][2] = nz ;
                                                        tail ++ ;
                                                        vis[nx][ny][nz] = 1 ;
                                                        ct ++ ;
                                                }
                                        }
                                }
                        }

                        if ( ct == tot )
                        {
                                return 1 ;
                        }
                        else
                        {
                                return 2 ;
                        }

                }
int Volume::countExt ( double  vox[3][3][3])

Definition at line 2838 of file volume.cpp.

References components26().

Referenced by isPiercable(), and isSimple().

                {
                        int tvox[3][3][3] ;

                        for ( int i = 0 ; i < 3 ; i ++ )
                                for ( int j = 0 ; j < 3 ; j ++ )
                                        for ( int k = 0 ; k < 3 ; k ++ )
                                        {
                                                if ( vox[i][j][k] < 0 )
                                                {
                                                        tvox[i][j][k] = 1 ;
                                                }
                                                else
                                                {
                                                        tvox[i][j][k] = 0 ;
                                                }
                                        }

                        return components26( tvox ) ;
                }
int Volume::countInt ( double  vox[3][3][3])

Definition at line 2859 of file volume.cpp.

References components6(), wustl_mm::SkeletonMaker::neighbor6, wustl_mm::SkeletonMaker::neighbor64, nnz, nx, and ny.

Referenced by isPiercable(), and isSimple().

                {
                        int i, j, k ;
                        int tvox[3][3][3] ;

                        for ( i = 0 ; i < 3 ; i ++ )
                                for ( j = 0 ; j < 3 ; j ++ )
                                        for ( k = 0 ; k < 3 ; k ++ )
                                        {
                                                tvox[i][j][k] = 0 ;
                                        }

                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                int nx = 1 + neighbor6[i][0] ;
                                int ny = 1 + neighbor6[i][1] ;
                                int nz = 1 + neighbor6[i][2] ;
                                if ( vox[ nx ][ ny ][ nz ] >= 0 )
                                {
                                        tvox[ nx ][ ny ][ nz ] = 1 ;
                                        for ( j = 0 ; j < 4 ; j ++ )
                                        {
                                                int nnx = neighbor64[i][j][0] + nx ;
                                                int nny = neighbor64[i][j][1] + ny ;
                                                int nnz = neighbor64[i][j][2] + nz ;
                                                if ( vox[ nnx ][ nny ][ nnz ] >= 0 )
                                                {
                                                        tvox[ nnx ][ nny ][ nnz ] = 1 ;
                                                }
                                        }
                                }
                        }

                        return components6( tvox ) ;
                }
int Volume::countIntEuler ( int  ox,
int  oy,
int  oz 
)

Definition at line 2895 of file volume.cpp.

References components6(), getDataAt(), wustl_mm::SkeletonMaker::neighbor6, wustl_mm::SkeletonMaker::neighbor64, nnz, nx, and ny.

Referenced by hasCompleteSheet().

                {
                        int nv = 0 , ne = 0, nc = 0 ;

                        int i, j, k ;
                        int tvox[3][3][3] ;
                        double vox[3][3][3] ;

                        for ( i = 0 ; i < 3 ; i ++ )
                                for ( j = 0 ; j < 3 ; j ++ )
                                        for ( k = 0 ; k < 3 ; k ++ )
                                        {
                                                vox[i][j][k] = getDataAt( ox - 1 + i, oy - 1 + j, oz - 1 + k ) ;
                                                tvox[i][j][k] = 0 ;
                                        }

                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                int nx = 1 + neighbor6[i][0] ;
                                int ny = 1 + neighbor6[i][1] ;
                                int nz = 1 + neighbor6[i][2] ;
                                if ( vox[nx][ny][nz] >= 0 )
                                {
                                        tvox[ nx ][ ny ][ nz ] = 1 ;

                                        nv ++ ;

                                        for ( j = 0 ; j < 4 ; j ++ )
                                        {
                                                int nnx = neighbor64[i][j][0] + nx ;
                                                int nny = neighbor64[i][j][1] + ny ;
                                                int nnz = neighbor64[i][j][2] + nz ;
                                                if ( vox[nnx][nny][nnz] >= 0 )
                                                {
                                                        if ( tvox[nnx][nny][nnz] == 0 )
                                                        {
                                                                tvox[nnx][nny][nnz] = 1 ;
                                                                nv ++ ;
                                                        }

                                                        ne ++ ;
                                                }
                                        }
                                }
                        }

                        nc = components6( tvox ) ;

                        return ( nc - ( nv - ne ) ) ;
                }
void Volume::curveSkeleton ( Volume grayvol,
float  lowthr,
float  highthr,
Volume svol 
)

* No topology check */

* Minimal topology check */ * set nodes for next layer * Faster version of erode2 using priority queue */ insert them back into priority queue * * Erode to atoms */ insert them back into priority queue * * Thin the current volume while preserving voxels with values > highthr or <= lowthr in grayvol

Definition at line 3812 of file volume.cpp.

References EMAN::PriorityQueue< ValueT, KeyT >::add(), getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isPiercable(), isSimple(), MAX_ERODE, MAX_QUEUELEN, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), v, Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z.

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning().

                {
                        int i, j, k ;
                        // First, threshold the volume
                        #ifdef VERBOSE
                        printf("Thresholding the volume to -1/0...\n") ;
                        #endif
                        threshold( 0.5f, -1, 0 ) ;

                        // Next, apply convergent erosion
                        // by preserving: complex nodes, curve end-points, and sheet points

                        // Next, initialize the linked queue
                        #ifdef VERBOSE
                        printf("Initializing queue...\n") ;
                        #endif
                        GridQueue2* queue2 = new GridQueue2( ) ;
                        GridQueue2* queue3 = new GridQueue2( ) ;
                        GridQueue2* queue4 = new GridQueue2( ) ;
                        PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) >= 0 )
                                                {
                                                        float v = (float)grayvol->getDataAt(i,j,k) ;
                                                        if ( v <= lowthr || v > highthr || svol->getDataAt(i,j,k) > 0 )
                                                        {
                                                                setDataAt( i, j, k, MAX_ERODE ) ;
                                                        }
                                                        else
                                                        {
                                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                                {
                                                                        if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
                                                                        {
                                                                                // setDataAt( i, j, k, 1 ) ;
                                                                                queue2->prepend( i, j, k ) ;
                                                                                break ;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                        int wid = MAX_ERODE ;
                        #ifdef VERBOSE
                        printf("Total %d nodes\n", queue2->getNumElements() ) ;
                        printf("Start erosion to %d...\n", wid) ;
                        #endif


                        // Perform erosion
                        gridQueueEle* ele ;
                        gridPoint* gp ;
                        int ox, oy, oz ;
                        int score ;
                        Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
                        for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
                        {
                                scrvol->setDataAt( i, -1 ) ;
                        }

        #ifdef  NOISE_DIS_HELIX
                        Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
        #endif

                        for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
                        {
                                // At the start of each iteration,
                                // queue2 holds all the nodes for this layer
                                // queue3 and queue are empty

                                int numComplex = 0, numSimple = 0 ;
                                #ifdef VERBOSE
                                printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
                                #endif

                                /*
                                We first need to assign curwid + 1 to every node in this layer
                                */
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        if ( getDataAt(ox,oy,oz) == curwid )
                                        {
                                                ele = queue2->remove() ;
                                        }
                                        else
                                        {
                                                setDataAt(ox,oy,oz, curwid) ;
                                                ele = queue2->getNext() ;
                                        }
                                }
                                queue4->reset() ;
                                ele = queue4->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        queue2->prepend(ox,oy,oz) ;
                                        ele = queue4->remove() ;
                                }

                                // Now queue2 holds all the nodes for this layer

        #ifdef NOISE_DIS_HELIX
                                /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
                                queue2->reset() ;

                                // First run
                                int flag = 0 ;
                                while ( ( ele = queue2->getNext() ) != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;
                                        if ( NOISE_DIS_HELIX <= 1 )
                                        {
                                                noisevol->setDataAt( ox, oy, oz, 0 ) ;
                                        }
                                        else
                                        {
                                                flag = 0 ;
                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) == 0 )
                                                        {
                                                                noisevol->setDataAt( ox, oy, oz, 1 ) ;
                                                                flag = 1 ;
                                                                break ;
                                                        }
                                                }
                                                if ( ! flag )
                                                {
                                                        noisevol->setDataAt( ox, oy, oz, 0 ) ;
                                                }
                                        }
                                }

                                int cur, visited ;
                                for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
                                {
                                        queue2->reset() ;
                                        int count = 0 ;
                                        visited = 0 ;

                                        while ( ( ele = queue2->getNext() ) != NULL )
                                        {
                                                ox = ele->x ;
                                                oy = ele->y ;
                                                oz = ele->z ;

                                                if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
                                                {
                                                        visited ++ ;
                                                        continue ;
                                                }

                                                flag = 0 ;
                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
                                                        {
                                                                noisevol->setDataAt( ox, oy, oz, 1 ) ;
                                                                visited ++ ;
                                                                count ++ ;
                                                                break ;
                                                        }
                                                }
                                        }

                                        if ( count == 0 )
                                        {
                                                break ;
                                        }
                                }
                                printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;


        #endif
                                /* Commented out for debugging

                                // First,
                                // check for complex nodes in queue2
                                // move them from queue2 to queue3
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        // Check simple
        #ifndef NOISE_DIS_HELIX
                                        if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
        #else
                                        if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
        #endif
                                        {
                                                // Complex, set to next layer
                                                setDataAt( ox, oy, oz, curwid + 1 ) ;
                                                queue3->prepend( ox, oy, oz ) ;
                                                ele = queue2->remove() ;

                                                numComplex ++ ;
                                        }
                                        else
                                        {
                                                ele = queue2->getNext() ;
                                        }
                                }
                                */

                                // Next,
                                // Compute score for each node left in queue2
                                // move them into priority queue
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        // Compute score
                                        score = getNumPotComplex2( ox, oy, oz ) ;
                                        scrvol->setDataAt( ox, oy, oz, score ) ;

                                        // Push to queue
                                        gp = new gridPoint ;
                                        gp->x = ox ;
                                        gp->y = oy ;
                                        gp->z = oz ;
                                        // queue->add( gp, -score ) ;
                                        queue->add( gp, score ) ;

                                        ele = queue2->remove() ;
                                }

                                // Rename queue3 to be queue2,
                                // Clear queue3
                                // From now on, queue2 holds nodes of next level
                                delete queue2 ;
                                queue2 = queue3 ;
                                queue3 = new GridQueue2( ) ;

                                // Next, start priority queue iteration
                                while ( ! queue->isEmpty() )
                                {
                                        // Retrieve the node with the highest score
                                        queue->remove( gp, score ) ;
                                        ox = gp->x ;
                                        oy = gp->y ;
                                        oz = gp->z ;
                                        delete gp ;
        //                              score = -score ;

                                        // Ignore the node
                                        // if it has been processed before
                                        // or it has an updated score
                                        if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
                                        {
                                                continue ;
                                        }

                                        /* Commented out for debugging

                                        // Remove this simple node
                                        setDataAt( ox, oy, oz, -1 ) ;
                                        numSimple ++ ;
                                        // printf("Highest score: %d\n", score) ;
                                        */

                                        /* Added for debugging */
                                        // Check simple
        #ifndef NOISE_DIS_HELIX
                                        // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
                                        if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
        #else
                                        if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
        #endif
                                        {
                                                // Complex, set to next layer
                                                setDataAt( ox, oy, oz, curwid + 1 ) ;
                                                queue4->prepend( ox, oy, oz ) ;
                                                numComplex ++ ;
                                        }
                                        else
                                        {
                                                setDataAt( ox, oy, oz, -1 ) ;
                                                numSimple ++ ;


                                        /* Adding ends */
                                                // Move its neighboring unvisited node to queue2
                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) == 0 )
                                                        {
                                                                // setDataAt( nx, ny, nz, curwid + 1 ) ;
                                                                queue2->prepend( nx, ny, nz ) ;
                                                        }
                                                }

                                        }


                                        /* Commented out for debugging

                                        // Find complex nodes in its 3x3 neighborhood
                                        // move them to queue2
                                        for ( i = -1 ; i < 2 ; i ++ )
                                                for ( j = -1 ; j < 2 ; j ++ )
                                                        for ( k = -1 ; k < 2 ; k ++ )
                                                        {
                                                                int nx = ox + i ;
                                                                int ny = oy + j ;
                                                                int nz = oz + k ;

                                                                // Check simple
                                                                if ( getDataAt( nx, ny, nz ) == curwid &&
                                                                        // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
        #ifndef NOISE_DIS_HELIX
                                                                        ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
        #else
                                                                        ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
        #endif

                                                                {
                                                                        // Complex, set to next layer
                                                                        setDataAt( nx, ny, nz, curwid + 1 ) ;
                                                                        queue2->prepend( nx, ny, nz ) ;
                                                                        numComplex ++ ;
                                                                }
                                                        }
                                        */

                                        // Update scores for nodes in its 5x5 neighborhood
                                        // insert them back into priority queue

                                        for ( i = -2 ; i < 3 ;i ++ )
                                                for ( j = -2 ; j < 3 ; j ++ )
                                                        for ( k = -2 ; k < 3 ; k ++ )
                                                        {
                                                                int nx = ox + i ;
                                                                int ny = oy + j ;
                                                                int nz = oz + k ;

                                                                if ( getDataAt( nx, ny, nz ) == curwid )
                                                                {
                                                                        // Compute score
                                                                        score = getNumPotComplex2( nx, ny, nz ) ;

                                                                        if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
                                                                        {
                                                                                // printf("Update\n") ;
                                                                                scrvol->setDataAt( nx, ny, nz, score ) ;
                                                                                // Push to queue
                                                                                gp = new gridPoint ;
                                                                                gp->x = nx ;
                                                                                gp->y = ny ;
                                                                                gp->z = nz ;
                                                                                // queue->add( gp, -score ) ;
                                                                                queue->add( gp, score ) ;
                                                                        }
                                                                }
                                                        }


                                }

                                #ifdef VERBOSE
                                printf("%d complex, %d simple\n", numComplex, numSimple) ;
                                #endif

                                if ( numSimple == 0 )
                                {
                                        if ( queue2->getNumElements() > 0 )
                                        {
                                                printf("*************************wierd here*************************\n");
                                        }
                                                break ;
                                }
                        }

                        // Remove all internal voxels (contained in manifold surfaces)
                        queue2->reset() ;
                        queue4->reset() ;
                        ele = queue4->getNext() ;
                        while ( ele != NULL )
                        {
                                ox = ele->x ;
                                oy = ele->y ;
                                oz = ele->z ;

                                if ( isPiercable(ox,oy,oz) == 1 )  // hasCompleteSheet( ox, oy, oz ) == 1 ) //
                                {
                                        queue2->prepend(ox,oy,oz) ;
                                //      setDataAt( ox, oy, oz, -1 ) ;
                                }
                                ele = queue4->remove() ;
                        }

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) == 0 && isPiercable(i,j,k) ) //hasCompleteSheet(i,j,k) == 1) //
                                                {
                                                        queue2->prepend( i, j, k ) ;
                                                }
                                        }
                        queue2->reset() ;
                        ele = queue2->getNext() ;
                        while ( ele != NULL )
                        {
                                ox = ele->x ;
                                oy = ele->y ;
                                oz = ele->z ;
                                setDataAt( ox, oy, oz, -1 ) ;
                                ele = queue2->remove() ;
                        }


                        // Finally, clean up
                        delete scrvol;
                        delete queue;
                        delete queue2;
                        delete queue3;
                        delete queue4;
                        #ifdef VERBOSE
                        printf("Thresholding the volume to 0/1...\n") ;
                        #endif
                        threshold( 0, 0, 1 ) ;
                }
void Volume::curveSkeleton ( float  thr,
Volume svol 
)

Definition at line 4268 of file volume.cpp.

References EMAN::PriorityQueue< ValueT, KeyT >::add(), getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isSimple(), MAX_ERODE, MAX_QUEUELEN, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z.

                {
                        int i, j, k ;
                        // First, threshold the volume
                        #ifdef VERBOSE
                        printf("Thresholding the volume to -1/0...\n") ;
                        #endif
                        threshold( thr, -1, 0 ) ;

                        // Next, apply convergent erosion
                        // by preserving: complex nodes, curve end-points, and sheet points

                        // Next, initialize the linked queue
                        #ifdef VERBOSE
                        printf("Initializing queue...\n") ;
                        #endif
                        GridQueue2* queue2 = new GridQueue2( ) ;
                        GridQueue2* queue3 = new GridQueue2( ) ;
                        GridQueue2* queue4 = new GridQueue2( ) ;
                        PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) >= 0 )
                                                {
                                                        if ( svol->getDataAt(i,j,k) > 0 )
                                                        {
                                                                setDataAt( i, j, k, MAX_ERODE ) ;
                                                        }
                                                        else
                                                        {
                                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                                {
                                                                        if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
                                                                        {
                                                                                // setDataAt( i, j, k, 1 ) ;
                                                                                queue2->prepend( i, j, k ) ;
                                                                                break ;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }

                        int wid = MAX_ERODE ;
                        #ifdef VERBOSE
                        printf("Total %d nodes\n", queue2->getNumElements() ) ;
                        printf("Start erosion to %d...\n", wid) ;
                        #endif


                        // Perform erosion
                        gridQueueEle* ele ;
                        gridPoint* gp ;
                        int ox, oy, oz ;
                        int score ;
                        Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
                        for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
                        {
                                scrvol->setDataAt( i, -1 ) ;
                        }

        #ifdef  NOISE_DIS_HELIX
                        Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
        #endif

                        for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
                        {
                                // At the start of each iteration,
                                // queue2 holds all the nodes for this layer
                                // queue3 and queue are empty

                                int numComplex = 0, numSimple = 0 ;
                                #ifdef VERBOSE
                                printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
                                #endif

                                /*
                                We first need to assign curwid + 1 to every node in this layer
                                */
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        if ( getDataAt(ox,oy,oz) == curwid )
                                        {
                                                ele = queue2->remove() ;
                                        }
                                        else
                                        {
                                                setDataAt(ox,oy,oz, curwid) ;
                                                ele = queue2->getNext() ;
                                        }
                                }
                                queue4->reset() ;
                                ele = queue4->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        queue2->prepend(ox,oy,oz) ;
                                        ele = queue4->remove() ;
                                }

                                // Now queue2 holds all the nodes for this layer

        #ifdef NOISE_DIS_HELIX
                                /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
                                queue2->reset() ;

                                // First run
                                int flag = 0 ;
                                while ( ( ele = queue2->getNext() ) != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;
                                        if ( NOISE_DIS_HELIX <= 1 )
                                        {
                                                noisevol->setDataAt( ox, oy, oz, 0 ) ;
                                        }
                                        else
                                        {
                                                flag = 0 ;
                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) == 0 )
                                                        {
                                                                noisevol->setDataAt( ox, oy, oz, 1 ) ;
                                                                flag = 1 ;
                                                                break ;
                                                        }
                                                }
                                                if ( ! flag )
                                                {
                                                        noisevol->setDataAt( ox, oy, oz, 0 ) ;
                                                }
                                        }
                                }

                                int cur, visited ;
                                for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
                                {
                                        queue2->reset() ;
                                        int count = 0 ;
                                        visited = 0 ;

                                        while ( ( ele = queue2->getNext() ) != NULL )
                                        {
                                                ox = ele->x ;
                                                oy = ele->y ;
                                                oz = ele->z ;

                                                if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
                                                {
                                                        visited ++ ;
                                                        continue ;
                                                }

                                                flag = 0 ;
                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
                                                        {
                                                                noisevol->setDataAt( ox, oy, oz, 1 ) ;
                                                                visited ++ ;
                                                                count ++ ;
                                                                break ;
                                                        }
                                                }
                                        }

                                        if ( count == 0 )
                                        {
                                                break ;
                                        }
                                }
                                printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;


        #endif
                                /* Commented out for debugging

                                // First,
                                // check for complex nodes in queue2
                                // move them from queue2 to queue3
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        // Check simple
        #ifndef NOISE_DIS_HELIX
                                        if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
        #else
                                        if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
        #endif
                                        {
                                                // Complex, set to next layer
                                                setDataAt( ox, oy, oz, curwid + 1 ) ;
                                                queue3->prepend( ox, oy, oz ) ;
                                                ele = queue2->remove() ;

                                                numComplex ++ ;
                                        }
                                        else
                                        {
                                                ele = queue2->getNext() ;
                                        }
                                }
                                */

                                // Next,
                                // Compute score for each node left in queue2
                                // move them into priority queue
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        // Compute score
                                        score = getNumPotComplex2( ox, oy, oz ) ;
                                        scrvol->setDataAt( ox, oy, oz, score ) ;

                                        // Push to queue
                                        gp = new gridPoint ;
                                        gp->x = ox ;
                                        gp->y = oy ;
                                        gp->z = oz ;
                                        // queue->add( gp, -score ) ;
                                        queue->add( gp, score ) ;

                                        ele = queue2->remove() ;
                                }

                                // Rename queue3 to be queue2,
                                // Clear queue3
                                // From now on, queue2 holds nodes of next level
                                delete queue2 ;
                                queue2 = queue3 ;
                                queue3 = new GridQueue2( ) ;

                                // Next, start priority queue iteration
                                while ( ! queue->isEmpty() )
                                {
                                        // Retrieve the node with the highest score
                                        queue->remove( gp, score ) ;
                                        ox = gp->x ;
                                        oy = gp->y ;
                                        oz = gp->z ;
                                        delete gp ;
        //                              score = -score ;

                                        // Ignore the node
                                        // if it has been processed before
                                        // or it has an updated score
                                        if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
                                        {
                                                continue ;
                                        }

                                        /* Commented out for debugging

                                        // Remove this simple node
                                        setDataAt( ox, oy, oz, -1 ) ;
                                        numSimple ++ ;
                                        // printf("Highest score: %d\n", score) ;
                                        */

                                        /* Added for debugging */
                                        // Check simple
        #ifndef NOISE_DIS_HELIX
                                        // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
                                        if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
        #else
                                        if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
        #endif
                                        {
                                                // Complex, set to next layer
                                                setDataAt( ox, oy, oz, curwid + 1 ) ;
                                                queue4->prepend( ox, oy, oz ) ;
                                                numComplex ++ ;
                                        }
                                        else
                                        {
                                                setDataAt( ox, oy, oz, -1 ) ;
                                                numSimple ++ ;
                                        }
                                        /* Adding ends */

                                        // Move its neighboring unvisited node to queue2
                                        for ( int m = 0 ; m < 6 ; m ++ )
                                        {
                                                int nx = ox + neighbor6[m][0] ;
                                                int ny = oy + neighbor6[m][1] ;
                                                int nz = oz + neighbor6[m][2] ;
                                                if ( getDataAt( nx, ny, nz ) == 0 )
                                                {
                                                        // setDataAt( nx, ny, nz, curwid + 1 ) ;
                                                        queue2->prepend( nx, ny, nz ) ;
                                                }
                                        }

                                        /* Commented out for debugging

                                        // Find complex nodes in its 3x3 neighborhood
                                        // move them to queue2
                                        for ( i = -1 ; i < 2 ; i ++ )
                                                for ( j = -1 ; j < 2 ; j ++ )
                                                        for ( k = -1 ; k < 2 ; k ++ )
                                                        {
                                                                int nx = ox + i ;
                                                                int ny = oy + j ;
                                                                int nz = oz + k ;

                                                                // Check simple
                                                                if ( getDataAt( nx, ny, nz ) == curwid &&
                                                                        // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
        #ifndef NOISE_DIS_HELIX
                                                                        ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
        #else
                                                                        ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
        #endif

                                                                {
                                                                        // Complex, set to next layer
                                                                        setDataAt( nx, ny, nz, curwid + 1 ) ;
                                                                        queue2->prepend( nx, ny, nz ) ;
                                                                        numComplex ++ ;
                                                                }
                                                        }
                                        */

                                        // Update scores for nodes in its 5x5 neighborhood
                                        // insert them back into priority queue

                                        for ( i = -2 ; i < 3 ;i ++ )
                                                for ( j = -2 ; j < 3 ; j ++ )
                                                        for ( k = -2 ; k < 3 ; k ++ )
                                                        {
                                                                int nx = ox + i ;
                                                                int ny = oy + j ;
                                                                int nz = oz + k ;

                                                                if ( getDataAt( nx, ny, nz ) == curwid )
                                                                {
                                                                        // Compute score
                                                                        score = getNumPotComplex2( nx, ny, nz ) ;

                                                                        if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
                                                                        {
                                                                                // printf("Update\n") ;
                                                                                scrvol->setDataAt( nx, ny, nz, score ) ;
                                                                                // Push to queue
                                                                                gp = new gridPoint ;
                                                                                gp->x = nx ;
                                                                                gp->y = ny ;
                                                                                gp->z = nz ;
                                                                                // queue->add( gp, -score ) ;
                                                                                queue->add( gp, score ) ;
                                                                        }
                                                                }
                                                        }


                                }

                                #ifdef VERBOSE
                                printf("%d complex, %d simple\n", numComplex, numSimple) ;
                                #endif

                                if ( numSimple == 0 )
                                {
                                                break ;
                                }
                        }

                        // Finally, clean up
                        delete scrvol;
                        delete queue;
                        delete queue2;
                        delete queue3;
                        delete queue4;
                        #ifdef VERBOSE
                        printf("Thresholding the volume to 0/1...\n") ;
                        #endif
                        threshold( 0, 0, 1 ) ;
                }
void Volume::curveSkeleton2D ( float  thr,
Volume svol 
)

Definition at line 4678 of file volume.cpp.

References EMAN::PriorityQueue< ValueT, KeyT >::add(), getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isSimple(), MAX_ERODE, MAX_QUEUELEN, wustl_mm::SkeletonMaker::neighbor4, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z.

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning().

                {
                        int i, j, k ;
                        // First, threshold the volume
                        #ifdef VERBOSE
                        printf("Thresholding the volume to -1/0...\n") ;
                        #endif
                        threshold( thr, -1, 0 ) ;

                        // Next, apply convergent erosion
                        // by preserving: complex nodes, curve end-points, and sheet points

                        // Next, initialize the linked queue
                        #ifdef VERBOSE
                        printf("Initializing queue...\n") ;
                        #endif
                        GridQueue2* queue2 = new GridQueue2( ) ;
                        GridQueue2* queue3 = new GridQueue2( ) ;
                        GridQueue2* queue4 = new GridQueue2( ) ;
                        PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) >= 0 )
                                                {
                                                        if ( svol->getDataAt(i,j,k) > 0 )
                                                        {
                                                                setDataAt( i, j, k, MAX_ERODE ) ;
                                                        }
                                                        else
                                                        {
                                                                for ( int m = 0 ; m < 4 ; m ++ )
                                                                {
                                                                        if ( getDataAt( i + neighbor4[m][0], j + neighbor4[m][1], k ) < 0 )
                                                                        {
                                                                                // setDataAt( i, j, k, 1 ) ;
                                                                                queue2->prepend( i, j, k ) ;
                                                                                break ;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                        int wid = MAX_ERODE ;
                        #ifdef VERBOSE
                        printf("Total %d nodes\n", queue2->getNumElements() ) ;
                        printf("Start erosion to %d...\n", wid) ;
                        #endif


                        // Perform erosion
                        gridQueueEle* ele ;
                        gridPoint* gp ;
                        int ox, oy, oz ;
                        int score ;
                        Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
                        for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
                        {
                                scrvol->setDataAt( i, -1 ) ;
                        }

        #ifdef  NOISE_DIS_HELIX
                        Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
        #endif

                        for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
                        {
                                // At the start of each iteration,
                                // queue2 holds all the nodes for this layer
                                // queue3 and queue are empty

                                int numComplex = 0, numSimple = 0 ;
                                #ifdef VERBOSE
                                printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
                                #endif

                                /*
                                We first need to assign curwid + 1 to every node in this layer
                                */
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        if ( getDataAt(ox,oy,oz) == curwid )
                                        {
                                                ele = queue2->remove() ;
                                        }
                                        else
                                        {
                                                setDataAt(ox,oy,oz, curwid) ;
                                                ele = queue2->getNext() ;
                                        }
                                }
                                queue4->reset() ;
                                ele = queue4->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        queue2->prepend(ox,oy,oz) ;
                                        ele = queue4->remove() ;
                                }

                                // Now queue2 holds all the nodes for this layer

                                #ifdef NOISE_DIS_HELIX
                                /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
                                queue2->reset() ;

                                // First run
                                int flag = 0 ;
                                while ( ( ele = queue2->getNext() ) != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;
                                        if ( NOISE_DIS_HELIX <= 1 )
                                        {
                                                noisevol->setDataAt( ox, oy, oz, 0 ) ;
                                        }
                                        else
                                        {
                                                flag = 0 ;
                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) == 0 )
                                                        {
                                                                noisevol->setDataAt( ox, oy, oz, 1 ) ;
                                                                flag = 1 ;
                                                                break ;
                                                        }
                                                }
                                                if ( ! flag )
                                                {
                                                        noisevol->setDataAt( ox, oy, oz, 0 ) ;
                                                }
                                        }
                                }

                                int cur, visited ;
                                for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
                                {
                                        queue2->reset() ;
                                        int count = 0 ;
                                        visited = 0 ;

                                        while ( ( ele = queue2->getNext() ) != NULL )
                                        {
                                                ox = ele->x ;
                                                oy = ele->y ;
                                                oz = ele->z ;

                                                if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
                                                {
                                                        visited ++ ;
                                                        continue ;
                                                }

                                                flag = 0 ;
                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
                                                        {
                                                                noisevol->setDataAt( ox, oy, oz, 1 ) ;
                                                                visited ++ ;
                                                                count ++ ;
                                                                break ;
                                                        }
                                                }
                                        }

                                        if ( count == 0 )
                                        {
                                                break ;
                                        }
                                }
                                printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;


                                #endif
                                /* Commented out for debugging

                                // First,
                                // check for complex nodes in queue2
                                // move them from queue2 to queue3
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        // Check simple
        #ifndef NOISE_DIS_HELIX
                                        if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
        #else
                                        if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
        #endif
                                        {
                                                // Complex, set to next layer
                                                setDataAt( ox, oy, oz, curwid + 1 ) ;
                                                queue3->prepend( ox, oy, oz ) ;
                                                ele = queue2->remove() ;

                                                numComplex ++ ;
                                        }
                                        else
                                        {
                                                ele = queue2->getNext() ;
                                        }
                                }
                                */

                                // Next,
                                // Compute score for each node left in queue2
                                // move them into priority queue
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        // Compute score
                                        score = getNumPotComplex2( ox, oy, oz ) ;
                                        //score = getNumNeighbor6( ox, oy, oz ) ;
                                        scrvol->setDataAt( ox, oy, oz, score ) ;

                                        // Push to queue
                                        gp = new gridPoint ;
                                        gp->x = ox ;
                                        gp->y = oy ;
                                        gp->z = oz ;
                                        // queue->add( gp, -score ) ;
                                        queue->add( gp, score ) ;

                                        ele = queue2->remove() ;
                                }

                                // Rename queue3 to be queue2,
                                // Clear queue3
                                // From now on, queue2 holds nodes of next level
                                delete queue2 ;
                                queue2 = queue3 ;
                                queue3 = new GridQueue2( ) ;

                                // Next, start priority queue iteration
                                while ( ! queue->isEmpty() )
                                {
                                        // Retrieve the node with the highest score
                                        queue->remove( gp, score ) ;
                                        ox = gp->x ;
                                        oy = gp->y ;
                                        oz = gp->z ;
                                        delete gp ;
        //                              score = -score ;

                                        // Ignore the node
                                        // if it has been processed before
                                        // or it has an updated score
                                        if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
                                        {
                                                continue ;
                                        }

                                        /* Commented out for debugging

                                        // Remove this simple node
                                        setDataAt( ox, oy, oz, -1 ) ;
                                        numSimple ++ ;
                                        // printf("Highest score: %d\n", score) ;
                                        */

                                        /* Added for debugging */
                                        // Check simple
                                        #ifndef NOISE_DIS_HELIX
                                        // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
                                        if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
                                        #else
                                        if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
                                        #endif
                                        {
                                                // Complex, set to next layer
                                                setDataAt( ox, oy, oz, curwid + 1 ) ;
                                                queue4->prepend( ox, oy, oz ) ;
                                                numComplex ++ ;
                                        }
                                        else
                                        {
                                                setDataAt( ox, oy, oz, -1 ) ;
                                                numSimple ++ ;
                                        }
                                        /* Adding ends */

                                        // Move its neighboring unvisited node to queue2
                                        for ( int m = 0 ; m < 4 ; m ++ )
                                        {
                                                int nx = ox + neighbor4[m][0] ;
                                                int ny = oy + neighbor4[m][1] ;
                                                int nz = oz ;
                                                if ( getDataAt( nx, ny, nz ) == 0 )
                                                {
                                                        // setDataAt( nx, ny, nz, curwid + 1 ) ;
                                                        queue2->prepend( nx, ny, nz ) ;
                                                }
                                        }

                                        /* Commented out for debugging

                                        // Find complex nodes in its 3x3 neighborhood
                                        // move them to queue2
                                        for ( i = -1 ; i < 2 ; i ++ )
                                                for ( j = -1 ; j < 2 ; j ++ )
                                                        for ( k = -1 ; k < 2 ; k ++ )
                                                        {
                                                                int nx = ox + i ;
                                                                int ny = oy + j ;
                                                                int nz = oz + k ;

                                                                // Check simple
                                                                if ( getDataAt( nx, ny, nz ) == curwid &&
                                                                        // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
        #ifndef NOISE_DIS_HELIX
                                                                        ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
        #else
                                                                        ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
        #endif

                                                                {
                                                                        // Complex, set to next layer
                                                                        setDataAt( nx, ny, nz, curwid + 1 ) ;
                                                                        queue2->prepend( nx, ny, nz ) ;
                                                                        numComplex ++ ;
                                                                }
                                                        }
                                        */

                                        // Update scores for nodes in its 5x5 neighborhood
                                        // insert them back into priority queue

                                        for ( i = -2 ; i < 3 ;i ++ )
                                                for ( j = -2 ; j < 3 ; j ++ )
                                                        for ( k = -2 ; k < 3 ; k ++ )
                                                        {
                                                                int nx = ox + i ;
                                                                int ny = oy + j ;
                                                                int nz = oz + k ;

                                                                if ( getDataAt( nx, ny, nz ) == curwid )
                                                                {
                                                                        // Compute score
                                                                        score = getNumPotComplex2( nx, ny, nz ) ;
                                                                        //score = getNumNeighbor6( nx, ny, nz ) ;

                                                                        if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
                                                                        {
                                                                                // printf("Update\n") ;
                                                                                scrvol->setDataAt( nx, ny, nz, score ) ;
                                                                                // Push to queue
                                                                                gp = new gridPoint ;
                                                                                gp->x = nx ;
                                                                                gp->y = ny ;
                                                                                gp->z = nz ;
                                                                                // queue->add( gp, -score ) ;
                                                                                queue->add( gp, score ) ;
                                                                        }
                                                                }
                                                        }


                                }

                                #ifdef VERBOSE
                                printf("%d complex, %d simple\n", numComplex, numSimple) ;
                                #endif

                                if ( numSimple == 0 )
                                {
                                                break ;
                                }
                        }

                        // Finally, clean up
                        #ifdef VERBOSE
                        printf("Thresholding the volume to 0/1...\n") ;
                        #endif
                        threshold( 0, 0, 1 ) ;
                        delete scrvol;
                        delete queue;
                        delete queue2;
                        delete queue3;
                        delete queue4;
                }
void Volume::erodeHelix ( )

Definition at line 5905 of file volume.cpp.

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PruneCurves().

                {
                        erodeHelix( 3 ) ;
                }
void Volume::erodeHelix ( int  disthr)

Definition at line 5911 of file volume.cpp.

References getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getSizeX(), getSizeY(), getSizeZ(), hasCompleteHelix(), wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridQueueEle::y, and wustl_mm::SkeletonMaker::gridQueueEle::z.

                {
                        int i, j, k ;
                        // First, threshold the volume
                        //printf("Thresholding the volume to -1/0...\n") ;
                        threshold( 0.1f, -1, 0 ) ;

                        /* Debug: remove faces */
                        //Volume* facevol = markFaceEdge() ;
                        /* End debugging */

                        // Next, initialize the linked queue
                        //printf("Initializing queue...\n") ;
                        Volume * fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
                        GridQueue2* queue2 = new GridQueue2( ) ;
                        GridQueue2* queue3 = new GridQueue2( ) ;
                        GridQueue2** queues = new GridQueue2* [ 10000 ] ;

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) >= 0 )
                                                {
                                                        if ( ! hasCompleteHelix( i, j, k ) )
                                                        // if ( ! hasCompleteHelix( i, j, k, facevol ) )
                                                        {
                                                                queue2->prepend( i, j, k ) ;
                                                                fvol->setDataAt( i, j, k, -1 ) ;
                                                        }
                                                }
                                        }
                        //printf("Total %d nodes\n", queue2->getNumElements() ) ;

                        // Start erosion
                        gridQueueEle* ele ;
                        int dis = -1 ;
                        while ( queue2->getNumElements() > 0 )
                        {
                                // First, set distance
                                dis -- ;
                                queues[ - dis ] = new GridQueue2( ) ;
                                //printf("Distance transform to %d...", dis) ;
                                queue2->reset() ;
                                while( ( ele = queue2->getNext() ) != NULL )
                                {
                                        setDataAt( ele->x, ele->y, ele->z, dis ) ;
                                        queues[ -dis ]->prepend( ele->x, ele->y, ele->z ) ;
                                }
                                //printf("%d nodes\n", queues[-dis]->getNumElements()) ;

                                // Next, find next layer
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        for ( int m = 0 ; m < 6 ; m ++ )
                                        {
                                                int nx = ele->x + neighbor6[m][0] ;
                                                int ny = ele->y + neighbor6[m][1] ;
                                                int nz = ele->z + neighbor6[m][2] ;
                                                if ( getDataAt( nx, ny, nz ) == 0 )
                                                {
                                                        fvol->setDataAt( nx, ny, nz, dis ) ;

                                                        if ( ! hasCompleteHelix( nx, ny, nz ) )
                                                        // if ( ! hasCompleteHelix( nx, ny, nz, facevol ) )
                                                        {
                                                                setDataAt( nx, ny, nz, 1 ) ;
                                                                queue3->prepend( nx, ny, nz ) ;
                                                        }
                                                }
                                        }

                                        ele = queue2->remove() ;
                                }

                                // Next, swap queues
                                GridQueue2* temp = queue2 ;
                                queue2 = queue3 ;
                                queue3 = temp ;
                        }

                        // Deal with closed rings
                        dis -- ;
                        queues[ - dis ] = new GridQueue2( ) ;
                        #ifdef VERBOSE
                        printf("Detecting closed rings %d...", dis) ;
                        #endif
                        int ftot = 0 ;
                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) == 0 )
                                                {
                                                        setDataAt( i, j, k, dis ) ;
                                                        queues[ -dis ]->prepend( i, j, k ) ;
        /*
                                                        int fval = (int) fvol->getDataAt( i, j, k ) ;
                                                        if ( fval == 0)
                                                        {
                                                                // queues[ -dis ]->prepend( i, j, k ) ;
                                                        }
                                                        else
                                                        {
                                                                setDataAt( i, j, k, fval - 1 ) ;
                                                                // queues[ -fval + 1 ]->prepend( i, j, k ) ;
                                                        }
        */
                                                        ftot ++ ;
                                                }
                                        }
                        #ifdef VERBOSE
                        printf("%d nodes\n", ftot) ;
                        #endif


                        // return ;

                        /* Find local minimum: to help determine erosion level
                        int cts[ 64 ] ;
                        for ( i = 0 ; i <= - dis ; i ++ )
                        {
                                cts[ i ] = 0 ;
                        }
                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                double val = getDataAt( i, j, k ) ;
                                                if ( val < -1 )
                                                {
                                                        int min = 1 ;
                                                        for ( int m = 0 ; m < 6 ; m ++ )
                                                        {
                                                                int nx = i + neighbor6[m][0] ;
                                                                int ny = j + neighbor6[m][1] ;
                                                                int nz = k + neighbor6[m][2] ;
                                                                if ( getDataAt( nx, ny, nz ) < val )
                                                                {
                                                                        min = 0 ;
                                                                        break ;
                                                                }
                                                        }

                                                        if ( min )
                                                        {
                                                                cts[ (int) - val ] ++ ;
                                                        }
                                                }
                                        }

                        for ( i = 2 ; i <= - dis ; i ++ )
                        {
                                printf("Local minima: %d with %d nodes.\n", -i, cts[ i ] ) ;
                        }
                        */

                        // Dilate back
                        // Starting from nodes with distance - 2 - disthr

                        if ( disthr + 2 > - dis )
                        {
                                disthr = - dis - 2 ;
                        }
                        int d;

                        for ( d = - dis ; d > disthr + 1 ; d -- )
                        {
                                queues[ d ]->reset() ;
                                while ( (ele = queues[ d ]->getNext() ) != NULL )
                                {
                                        setDataAt( ele->x, ele->y, ele->z, d ) ;
                                }
                        }


                        for ( d = disthr + 1 ; d >= 2 ; d -- )
                        {
                                //delete queue3 ;
                                //queue3 = new GridQueue2( ) ;
                                queues[ d ]->reset() ;
                                ele = queues[ d ]->getNext() ;
                                while ( ele != NULL )
                                {
                                        int dilatable = 0 ;
                                        for ( int m = 0 ; m < 6 ; m ++ )
                                                        {
                                                                int nx = ele->x + neighbor6[m][0] ;
                                                                int ny = ele->y + neighbor6[m][1] ;
                                                                int nz = ele->z + neighbor6[m][2] ;
                                                                if ( getDataAt( nx, ny, nz ) == d + 1 )
                                                                {
                                                                        dilatable = 1 ;
                                                                        break ;
                                                                }
                                                        }


                                        if ( ! dilatable )
                                        {
                                                /*
                                                setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
                                                */

                                                setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
                                                if ( d > 2 )
                                                {
                                                        queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
                                                }
                                                ele = queues[ d ]->remove() ;
                                        }
                                        else
                                        {
                                                setDataAt( ele->x, ele->y, ele->z, d ) ;
                                                ele = queues[ d ]->getNext() ;
                                        }

                                }

                                /* Debug: extract minimal set */
                                while ( 1 )
                                {
                                        int num = 0 ;
                                        queues[ d ]->reset() ;
                                        ele = queues[ d ]->getNext() ;
                                        while ( ele != NULL )
                                        {
                                                int critical = 0 ;
                                                setDataAt( ele->x, ele->y, ele->z, -1 ) ;

                                                for ( i = -1 ; i < 2 ; i ++ )
                                                {
                                                        for ( j = -1 ; j < 2 ; j ++ )
                                                        {
                                                                for ( k = -1 ; k < 2 ; k ++ )
                                                                {
                                                                        if ( i != 0 && j != 0 && k != 0 )
                                                                        {
                                                                                continue ;
                                                                        }
                                                                        int nx = ele->x + i ;
                                                                        int ny = ele->y + j ;
                                                                        int nz = ele->z + k ;
                                                                        if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteHelix( nx,ny,nz ) ) //, facevol ) )
                                                                        {
                                                                                critical = 1 ;
                                                                                break ;
                                                                        }
                                                                }
                                                                if ( critical )
                                                                {
                                                                        break ;
                                                                }
                                                        }
                                                        if ( critical )
                                                        {
                                                                break ;
                                                        }
                                                }

                                                if ( critical )
                                                {
                                                        setDataAt( ele->x, ele->y, ele->z, d ) ;
                                                        ele = queues[ d ]->getNext() ;
                                                }
                                                else
                                                {
                                                        setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
                                                        if ( d > 2 )
                                                        {
                                                                queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
                                                        }
                                                        ele = queues[ d ]->remove() ;
                                                        num ++ ;
                                                }

                                        }

                                        #ifdef VERBOSE
                                        printf("Non-minimal: %d\n", num) ;
                                        #endif

                                        if ( num == 0 )
                                        {
                                                break ;
                                        }
                                }


                                /* End of debugging */

                                /*
                                queue3->reset() ;
                                ele = queue3->getNext() ;
                                while ( ele != NULL )
                                {
                                        setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
                                        ele = queue3->remove() ;
                                }
                                */
                        }

                        // Finally, threshold the volume
                        #ifdef VERBOSE
                        //printf("Thresholding the volume to 0/1...\n") ;
                        #endif
                        //threshold( -1, 1, 0, 0 ) ;
                        threshold( 0, 0, 1 ) ;
                        delete fvol ;
                        delete queue2;
                        delete queue3;
                        for ( d = -dis ; d >= 2 ; d -- ) {
                                delete queues[d];
                        }
                        delete [] queues;

                }
int Volume::erodeSheet ( )

Definition at line 6234 of file volume.cpp.

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PruneSurfaces().

                {
                        return erodeSheet( 3 ) ;
                }
int Volume::erodeSheet ( int  disthr)

Definition at line 6240 of file volume.cpp.

References getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getSizeX(), getSizeY(), getSizeZ(), hasCompleteSheet(), markCellFace(), nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), wustl_mm::SkeletonMaker::sheetNeighbor, threshold(), and Volume().

                {
                        int i, j, k ;
                        // First, threshold the volume
                        //printf("Thresholding the volume to -1/0...\n") ;
                        threshold( 0.1f, -1, 0 ) ;

                        /* Debug: remove cells */
                        Volume* facevol = markCellFace() ;
                        /* End debugging */

                        // Next, initialize the linked queue
                        //printf("Initializing queue...\n") ;
                        Volume * fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
                        GridQueue2* queue2 = new GridQueue2( ) ;
                        GridQueue2* queue3 = new GridQueue2( ) ;
                        GridQueue2** queues = new GridQueue2* [ 10000 ] ;

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) >= 0 )
                                                {
                                                        if ( ! hasCompleteSheet( i, j, k ) )
                                                        //if ( ! hasCompleteSheet( i, j, k, facevol ) )
                                                        {
                                                                queue2->prepend( i, j, k ) ;
                                                                fvol->setDataAt( i, j, k, -1 ) ;
                                                        }
                                                }
                                        }
                        #ifdef VERBOSE
                        printf("Total %d nodes\n", queue2->getNumElements() ) ;
                        #endif

                        // Start erosion
                        gridQueueEle* ele ;
                        int dis = -1 ;
                        while ( queue2->getNumElements() > 0 )
                        {
                                // First, set distance
                                dis -- ;
                                queues[ - dis ] = new GridQueue2( ) ;
                                //printf("Distance transform to %d...", dis) ;
                                queue2->reset() ;
                                while( ( ele = queue2->getNext() ) != NULL )
                                {
                                        setDataAt( ele->x, ele->y, ele->z, dis ) ;
                                        queues[ -dis ]->prepend( ele->x, ele->y, ele->z ) ;
                                }
                                //printf("%d nodes\n", queues[-dis]->getNumElements()) ;

                                // Next, find next layer
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        // for ( int m = 0 ; m < 6 ; m ++ )
                                        for ( int mx = -1 ; mx < 2 ; mx ++ )
                                                for ( int my = -1 ; my < 2 ; my ++ )
                                                        for ( int mz = -1 ; mz < 2 ; mz ++ )
                                                        {
                                                                if ( mx != 0 && my != 0 && mz != 0 )
                                                                {
                                                                        continue ;
                                                                }
                                                                //int nx = ele->x + neighbor6[m][0] ;
                                                                //int ny = ele->y + neighbor6[m][1] ;
                                                                //int nz = ele->z + neighbor6[m][2] ;
                                                                int nx = ele->x + mx ;
                                                                int ny = ele->y + my ;
                                                                int nz = ele->z + mz ;

                                                                if ( getDataAt( nx, ny, nz ) == 0 )
                                                                {
                                                                        fvol->setDataAt( nx, ny, nz, dis ) ;

                                                                        if  ( ! hasCompleteSheet( nx, ny, nz ) )
                                                                        // if  ( ! hasCompleteSheet( nx, ny, nz, facevol ) )
                                                                        {
                                                                                setDataAt( nx, ny, nz, 1 ) ;
                                                                                queue3->prepend( nx, ny, nz ) ;
                                                                        }
                                                                }
                                                        }

                                        ele = queue2->remove() ;
                                }

                                // Next, swap queues
                                GridQueue2* temp = queue2 ;
                                queue2 = queue3 ;
                                queue3 = temp ;
                        }

                        /* Deal with closed rings */

                        dis -- ;
                        queues[ - dis ] = new GridQueue2( ) ;
                        #ifdef VERBOSE
                        printf("Detecting closed rings %d...", dis) ;
                        #endif
                        int ftot = 0 ;
                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) == 0 )
                                                {
                                                        /*
                                                        int fval = (int) fvol->getDataAt( i, j, k ) ;
                                                        if ( fval == 0)
                                                        {
                                                                setDataAt( i, j, k, dis - 2 ) ;
                                                                // queues[ -dis ]->prepend( i, j, k ) ;
                                                        }
                                                        else
                                                        {
                                                                setDataAt( i, j, k, fval - 1 ) ;
                                                                queues[ -fval + 1 ]->prepend( i, j, k ) ;
                                                        }
                                                        */
                                                        setDataAt( i, j, k, dis ) ;
                                                        queues[ -dis ]->prepend( i, j, k ) ;

                                                        ftot ++ ;
                                                }
                                        }
                        #ifdef VERBOSE
                        printf("%d nodes\n", ftot) ;
                        #endif


                        /* Find local minimum: to help determine erosion level
                        int cts[ 64 ] ;
                        for ( i = 0 ; i <= - dis ; i ++ )
                        {
                                cts[ i ] = 0 ;
                        }
                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                double val = getDataAt( i, j, k ) ;
                                                if ( val < -1 )
                                                {
                                                        int min = 1 ;
                                                        for ( int m = 0 ; m < 6 ; m ++ )
                                                        {
                                                                int nx = i + neighbor6[m][0] ;
                                                                int ny = j + neighbor6[m][1] ;
                                                                int nz = k + neighbor6[m][2] ;
                                                                if ( getDataAt( nx, ny, nz ) < val )
                                                                {
                                                                        min = 0 ;
                                                                        break ;
                                                                }
                                                        }

                                                        if ( min )
                                                        {
                                                                cts[ (int) - val ] ++ ;
                                                        }
                                                }
                                        }

                        for ( i = 2 ; i <= - dis ; i ++ )
                        {
                                printf("Local minima: %d with %d nodes.\n", -i, cts[ i ] ) ;
                        }
                        */

                        // return ;

                        // Dilate back
                        // Starting from nodes with distance - 2 - disthr
                        int d ;
                        if ( disthr + 2 > - dis )
                        {
                                disthr = - dis - 2 ;

                        }
                        for ( d = - dis ; d > disthr + 1 ; d -- )
                        {
                                queues[ d ]->reset() ;
                                while ( (ele = queues[ d ]->getNext() ) != NULL )
                                {
                                        setDataAt( ele->x, ele->y, ele->z, d ) ;
                                }
                        }

                        for (d = disthr + 1 ; d >= 2 ; d -- )
                        {

                                //delete queue3 ;
                                //queue3 = new GridQueue2( ) ;
                                queues[ d ]->reset() ;
                                ele = queues[ d ]->getNext() ;
                                while ( ele != NULL )
                                {
                                        int dilatable = 0 ;
                                        // for ( int m = 0 ; m < 6 ; m ++ )
                                        /*
                                        for ( int mx = -1 ; mx < 2 ; mx ++ )
                                        for ( int my = -1 ; my < 2 ; my ++ )
                                        for ( int mz = -1 ; mz < 2 ; mz ++ )
                                        {
                                        if ( mx == 0 || my == 0 || mz == 0 )
                                        {
                                        int nx = ele->x + mx ; // neighbor6[m][0] ;
                                        int ny = ele->y + my ; // neighbor6[m][1] ;
                                        int nz = ele->z + mz ; // neighbor6[m][2] ;
                                        if ( getDataAt( nx, ny, nz ) == - d - 1 )
                                        {
                                        dilatable = 1 ;
                                        break ;
                                        }
                                        }
                                        }
                                        */
                                        for ( i = 0 ; i < 12 ; i ++ )
                                        {
                                                int flag = 1, flag2 = 0 ;
                                                for ( j = 0 ; j < 4 ; j ++ )
                                                {
                                                        int nx = ele->x + sheetNeighbor[i][j][0] ;
                                                        int ny = ele->y + sheetNeighbor[i][j][1] ;
                                                        int nz = ele->z + sheetNeighbor[i][j][2] ;

                                                        double val = getDataAt( nx, ny, nz ) ;

                                                        if ( val > - d && val < 0)
                                                        {
                                                                flag = 0 ;
                                                                break ;
                                                        }
                                                        if ( val == d + 1 )
                                                        {
                                                                flag2 ++ ;
                                                        }
                                                }

                                                if ( flag && flag2 )
                                                {
                                                        dilatable = 1 ;
                                                        break ;
                                                }
                                        }

                                        if ( ! dilatable )
                                        {
                                                // setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
                                                // queue3->prepend( ele->x, ele->y, ele->z ) ;

                                                setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
                                                if ( d > 2 )
                                                {
                                                        queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
                                                }
                                                ele = queues[ d ]->remove() ;
                                        }
                                        else
                                        {
                                                setDataAt( ele->x, ele->y, ele->z, d ) ;
                                                ele = queues[ d ]->getNext() ;
                                        }
                                }

                                /* Debug: extract minimal set */
                                while ( 1 )
                                {
                                        int num = 0 ;
                                        queues[ d ]->reset() ;
                                        ele = queues[ d ]->getNext() ;
                                        while ( ele != NULL )
                                        {
                                                int critical = 0 ;
                                                setDataAt( ele->x, ele->y, ele->z, -1 ) ;

                                                for ( i = -1 ; i < 2 ; i ++ )
                                                {
                                                        for ( j = -1 ; j < 2 ; j ++ )
                                                        {
                                                                for ( k = -1 ; k < 2 ; k ++ )
                                                                {
                                                                        if ( i != 0 && j != 0 && k != 0 )
                                                                        {
                                                                                continue ;
                                                                        }
                                                                        int nx = ele->x + i ;
                                                                        int ny = ele->y + j ;
                                                                        int nz = ele->z + k ;
                                                                        // if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteSheet( nx,ny,nz, facevol ) )
                                                                        if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteSheet( nx,ny,nz ) )
                                                                        {
                                                                                critical = 1 ;
                                                                                break ;
                                                                        }
                                                                }
                                                                if ( critical )
                                                                {
                                                                        break ;
                                                                }
                                                        }
                                                        if ( critical )
                                                        {
                                                                break ;
                                                        }
                                                }

                                                if ( critical )
                                                {
                                                        setDataAt( ele->x, ele->y, ele->z, d ) ;
                                                        ele = queues[ d ]->getNext() ;
                                                }
                                                else
                                                {
                                                        setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
                                                        if ( d > 2 )
                                                        {
                                                                queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
                                                        }
                                                        ele = queues[ d ]->remove() ;
                                                        num ++ ;
                                                }

                                        }
                                        #ifdef VERBOSE
                                        printf("Non-minimal: %d\n", num) ;
                                        #endif

                                        if ( num == 0 )
                                        {
                                                break ;
                                        }
                                }


                                /* End of debugging */

                                /*
                                queue3->reset() ;
                                ele = queue3->getNext() ;
                                while ( ele != NULL )
                                {
                                        setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
                                        ele = queue3->remove() ;
                                }
                                */
                        }


                        // Finally, threshold the volume
                        #ifdef VERBOSE
                        //printf("Thresholding the volume to 0/1...\n") ;
                        #endif
                        //threshold( -1, 1, 0, 0 ) ;
                        threshold( 0, 0, 1 ) ;

                        delete facevol ;
                        delete fvol ;
                        delete queue2;
                        delete queue3;
                        for (d = -dis ; d >= 2 ; d -- ) {
                                delete queues[d];
                        }
                        delete [] queues;

                        return - dis - 1 ;
                }
EMData * Volume::get_emdata ( )

Definition at line 31 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::get_emdata(), and getVolumeData().

                {
                        return this->getVolumeData()->get_emdata();
                }
double Volume::getDataAt ( int  x,
int  y,
int  z 
)
double Volume::getDataAt ( int  index)

Definition at line 65 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::GetDataAt(), and volData.

                                                    {
                        return volData->GetDataAt(index);
                }
int Volume::getIndex ( int  x,
int  y,
int  z 
)
int Volume::getNumNeighbor6 ( int  ox,
int  oy,
int  oz 
)

Definition at line 661 of file volume.cpp.

References getDataAt(), wustl_mm::SkeletonMaker::neighbor6, nx, and ny.

Referenced by getNumPotComplex(), isFeatureFace(), and isHelixEnd().

                                                                    {
                        int rvalue = 0 ;
                        for ( int i = 0 ; i < 6 ; i ++ ) {
                                int nx = ox + neighbor6[i][0] ;
                                int ny = oy + neighbor6[i][1] ;
                                int nz = oz + neighbor6[i][2] ;
                                if ( getDataAt( nx, ny, nz ) >= 0 ) {
                                        rvalue ++ ;
                                }
                        }

                        return rvalue ;
                }
int Volume::getNumPotComplex ( int  ox,
int  oy,
int  oz 
)

Definition at line 2548 of file volume.cpp.

References getDataAt(), getNumNeighbor6(), wustl_mm::SkeletonMaker::neighbor6, nx, ny, and setDataAt().

Referenced by getNumPotComplex2(), and surfaceSkeletonPres().

                {
                        //return 0 ;

                        int i;
                        double val = getDataAt( ox, oy, oz ) ;
                        if ( val <= 0 )
                        {
                //              return 70 ;
                        }

                        // return getNumNeighbor6( ox, oy, oz ) ;

                        // int v = ((getNumNeighbor6( ox, oy, oz ) & 255) << 24) ;
                        //int v = 0  ;

                        int rvalue = 0, nx, ny, nz ;
                        setDataAt( ox, oy, oz, -val ) ;

                        /*
                        for ( i = -1 ; i < 2 ; i ++ )
                                for ( j = -1 ; j < 2 ; j ++ )
                                        for ( k = -1 ; k < 2 ; k ++ )
                                        {
                                                nx = ox + i ;
                                                ny = oy + j ;
                                                nz = oz + k ;
                                                if ( getDataAt( nx, ny, nz ) == val )
                                                {
                                                        if ( isSheetEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
                                                        {
                                                                rvalue ++ ;
                                                        }
                                                }
                                        }
                        */

                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                nx = ox + neighbor6[i][0] ;
                                ny = oy + neighbor6[i][1] ;
                                nz = oz + neighbor6[i][2] ;

                                if ( getDataAt(nx,ny,nz) >= 0 )
                                {
                                        int num = getNumNeighbor6( nx, ny, nz ) ;
                                        if ( num > rvalue )
                                        {
                                                rvalue = num ;
                                        }
                                }
                        }


                        setDataAt( ox, oy, oz, val ) ;

                        return rvalue + getNumNeighbor6( ox, oy, oz ) * 10 ;
                        /*
                        int v = (((rvalue + getNumNeighbor6( ox, oy, oz ) * 10) & 255) << 24) ;
                        v |= ( ( ox & 255 ) << 16 )  ;
                        v |= ( ( oy & 255 ) << 8 ) ;
                        v |= ( ( oz & 255 ) ) ;
                        return v ;
                        */

                }
int Volume::getNumPotComplex2 ( int  ox,
int  oy,
int  oz 
)

Definition at line 2615 of file volume.cpp.

References getNumPotComplex().

Referenced by curveSkeleton(), curveSkeleton2D(), and skeleton().

                {
                        return getNumPotComplex( ox, oy, oz ) ;

                        //int i, j, k ;
                        //double val = getDataAt( ox, oy, oz ) ;
                        //if ( val <= 0 )
                        //{
                        //      return 0 ;
                        //}

                        //int rvalue = 0, nx, ny, nz ;
                        //setDataAt( ox, oy, oz, -val ) ;

                        //for ( i = -1 ; i < 2 ; i ++ )
                        //      for ( j = -1 ; j < 2 ; j ++ )
                        //              for ( k = -1 ; k < 2 ; k ++ )
                        //              {
                        //                      nx = ox + i ;
                        //                      ny = oy + j ;
                        //                      nz = oz + k ;
                        //                      if ( getDataAt( nx, ny, nz ) == val )
                        //                      {
                        //                              if ( isHelixEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
                        //                              {
                        //                                      rvalue ++ ;
                        //                              }
                        //                      }
                        //              }

                        //setDataAt( ox, oy, oz, val ) ;

                        //return rvalue + getNumNeighbor6( ox, oy, oz ) * 30 ;
                }
float Volume::getOriginX ( )

Definition at line 93 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::GetOriginX(), and volData.

                                         {
                        return volData->GetOriginX();
                }
float Volume::getOriginY ( )

Definition at line 97 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::GetOriginY(), and volData.

                                         {
                        return volData->GetOriginY();
                }
float Volume::getOriginZ ( )

Definition at line 101 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::GetOriginZ(), and volData.

                                         {
                        return volData->GetOriginZ();
                }
int Volume::getSizeX ( )
int Volume::getSizeY ( )
int Volume::getSizeZ ( )
float Volume::getSpacingX ( )

Definition at line 81 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::GetSpacingX(), and volData.

                                          {
                        return volData->GetSpacingX();
                }
float Volume::getSpacingY ( )

Definition at line 85 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::GetSpacingY(), and volData.

                                          {
                        return volData->GetSpacingY();
                }
float Volume::getSpacingZ ( )

Definition at line 89 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::GetSpacingZ(), and volData.

                                          {
                        return volData->GetSpacingZ();
                }
VolumeData * Volume::getVolumeData ( )

Definition at line 69 of file volume.cpp.

References volData.

Referenced by get_emdata(), and Volume().

                                                   {
                        return volData;
                }
int Volume::hasCell ( int  ox,
int  oy,
int  oz 
)

Definition at line 1016 of file volume.cpp.

References getDataAt().

Referenced by hasCompleteSheet(), and markCellFace().

                                                            {
                        for ( int i = 0 ; i < 2 ; i ++ )
                                for ( int j = 0 ; j < 2 ; j ++ )
                                        for ( int k = 0 ; k < 2 ; k ++ ) {
                                                if ( getDataAt( ox + i, oy + j, oz + k ) < 0 ) {
                                                        return 0 ;
                                                }
                                        }
                        return 1 ;
                }
int Volume::hasCompleteHelix ( int  ox,
int  oy,
int  oz 
)

*

Definition at line 1456 of file volume.cpp.

References getDataAt(), wustl_mm::SkeletonMaker::neighbor6, nx, and ny.

Referenced by erodeHelix().

                {
                        // Returns 1 if it has a complete helix
                        int i ;
                        int c1 = 0;
                        int nx, ny, nz ;
                        int j ;

                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                nx = ox + neighbor6[i][0] ;
                                ny = oy + neighbor6[i][1] ;
                                nz = oz + neighbor6[i][2] ;
                                if ( getDataAt( nx, ny, nz ) >= 0 )
                                {
                                        c1 ++ ;
                                        j = i ;
                                }

                        }

                        if ( c1 > 1 ) // || c1 == 0 )
                        {
                                return 1 ;
                        }

                        return 0 ;

                        /*
                        ox = ox + neighbor6[j][0] ;
                        oy = oy + neighbor6[j][1] ;
                        oz = oz + neighbor6[j][2] ;
                        c1 = 0 ;
                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                nx = ox + neighbor6[i][0] ;
                                ny = oy + neighbor6[i][1] ;
                                nz = oz + neighbor6[i][2] ;
                                if ( getDataAt( nx, ny, nz ) >= 0 )
                                {
                                        c1 ++ ;
                                }

                        }

                        if ( c1 > 1 )
                        {
                                return 0 ;
                        }
                        else
                        {
                                return 1 ;
                        }
                        */
                }
int Volume::hasCompleteHelix ( int  ox,
int  oy,
int  oz,
Volume fvol 
)

Definition at line 1512 of file volume.cpp.

References getDataAt(), wustl_mm::SkeletonMaker::neighbor6, nx, and ny.

                {

                        int i ;
                        int c1 = 0;
                        int nx, ny, nz ;
                        int j ;

                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                nx = ox + neighbor6[i][0] ;
                                ny = oy + neighbor6[i][1] ;
                                nz = oz + neighbor6[i][2] ;
                                if ( getDataAt( nx, ny, nz ) >= 0 )
                                {
                                        if ( i % 2 == 0 )
                                        {
                                                nx = ox ;
                                                ny = oy ;
                                                nz = oz ;
                                        }

                                        int val = (int)fvol->getDataAt( nx, ny, nz) ;
                                        if ( (( val >> ( 2 * ( i / 2 ) ) ) & 1) == 0 )
                                        {
                                                c1 ++ ;
                                                j = i ;
                                        }
                                }

                        }

                        if ( c1 > 1 )
                        {
                                return 1 ;
                        }

                        return 0 ;

                        /*
                        ox = ox + neighbor6[j][0] ;
                        oy = oy + neighbor6[j][1] ;
                        oz = oz + neighbor6[j][2] ;
                        c1 = 0 ;
                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                nx = ox + neighbor6[i][0] ;
                                ny = oy + neighbor6[i][1] ;
                                nz = oz + neighbor6[i][2] ;
                                if ( getDataAt( nx, ny, nz ) >= 0 )
                                {
                                        c1 ++ ;
                                }

                        }

                        if ( c1 > 1 )
                        {
                                return 0 ;
                        }
                        else
                        {
                                return 1 ;
                        }
                        */
                }
int Volume::hasCompleteSheet ( int  ox,
int  oy,
int  oz 
)

Definition at line 1322 of file volume.cpp.

References countIntEuler().

                                                                     {
                        // Returns 1 if it lies in the middle of a sheet
                        int temp = countIntEuler( ox, oy, oz ) ;
                        if ( temp > 0 )
                        {
                                return 1 ;
                        }
                        else
                        {
                                return 0 ;
                        }
                }
int Volume::hasCompleteSheet ( int  ox,
int  oy,
int  oz,
Volume fvol 
)

Definition at line 1174 of file volume.cpp.

References wustl_mm::SkeletonMaker::edgeFaces, wustl_mm::SkeletonMaker::faceCells, wustl_mm::SkeletonMaker::faceEdges, getDataAt(), hasCell(), nx, ny, wustl_mm::SkeletonMaker::sheetNeighbor, and tot.

Referenced by erodeSheet(), and isSheetEnd().

                                                                                   {
                        int i, j, k ;
                        int nx, ny, nz ;

                        int edge[6] = { 0,0,0,0,0,0 } ;
                        int faceflag[ 12 ] ;
                        int tot = 0 ;
                        int cellflag[ 8 ] ;

                        int ct = 0 ;
                        for (  i = -1 ; i < 1 ; i ++ )
                                for (  j = -1 ; j < 1 ; j ++ )
                                        for (  k = -1 ; k < 1 ; k ++ )
                                        {
                                                if ( hasCell( ox + i, oy + j, oz + k ) )
                                                {
                                                        cellflag[ ct ] = 1 ;
                                                }
                                                else
                                                {
                                                        cellflag[ ct ] = 0 ;
                                                }
                                                ct ++ ;
                                        }

                        for ( i = 0 ; i < 12 ; i ++ )
                        {
                                faceflag[ i ] = 1 ;
                                for ( j = 0 ; j < 4 ; j ++ )
                                {
                                        nx = ox + sheetNeighbor[i][j][0] ;
                                        ny = oy + sheetNeighbor[i][j][1] ;
                                        nz = oz + sheetNeighbor[i][j][2] ;

                                        if ( getDataAt( nx, ny, nz ) < 0 )
                                        {
                                                faceflag[ i ] = 0 ;
                                                break ;
                                        }
                                }

                                if ( faceflag[ i ] )
                                {
                                        if ( cellflag[ faceCells[i][0] ] ^ cellflag[ faceCells[i][1] ] )
                                        {
                                                int v1 = (int)( fvol->getDataAt(
                                                        ox - 1 + (( faceCells[i][0] >> 2 ) & 1 ),
                                                        oy - 1 + (( faceCells[i][0] >> 1 ) & 1 ),
                                                        oz - 1 + (( faceCells[i][0] ) & 1)) ) ;
                                                int v2 = (int)( fvol->getDataAt(
                                                        ox - 1 + (( faceCells[i][1] >> 2 ) & 1 ),
                                                        oy - 1 + (( faceCells[i][1] >> 1 ) & 1 ),
                                                        oz - 1 + (( faceCells[i][1] ) & 1)) ) ;
                                                if ( ((v1 >> (2 * (2 - i/4))) & 1) ||
                                                         ((v2 >> (2 * (2 - i/4) + 1 )) & 1) )
                                                {
                                                        faceflag[ i ] = 0 ;
                                                }
                                        }
                                }

                                if ( faceflag[ i ] )
                                {
                                        int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
                                        edge[ e0 ] ++ ;
                                        edge[ e1 ] ++ ;
                                        tot ++ ;
                                }
                        }

                        // Removing 1s
                        int numones = 0 ;
                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                if ( edge[ i ] == 1 )
                                {
                                        numones ++ ;
                                }
                        }
                        while( numones > 0 )
                        {
                                int e = 0;
                                for ( i = 0 ; i < 6 ; i ++ )
                                {
                                        if ( edge[ i ] == 1 )
                                        {
                                                e = i ;
                                                break ;
                                        }
                                }
                                /*
                                if ( edge[ e ] != 1 )
                                {
                                        printf("Wrong Again!********\n") ;
                                }
                                */

                                int f, e2 ;
                                for ( j = 0 ; j < 4 ; j ++ )
                                {
                                        f = edgeFaces[ e ][ j ] ;
                                        if ( faceflag[ f ] )
                                        {
                                                break ;
                                        }
                                }

                                /*
                                if ( faceflag[ f ] == 0 )
                                {
                                        printf("Wrong!********\n") ;
                                }
                                */

                                if ( faceEdges[ f ][ 0 ] == e )
                                {
                                        e2 = faceEdges[ f ][ 1 ] ;
                                }
                                else
                                {
                                        e2 = faceEdges[ f ][ 0 ] ;
                                }


                                edge[ e ] -- ;
                                numones -- ;
                                edge[ e2 ] -- ;
                                faceflag[ f ] = 0 ;
                                tot -- ;

                                if ( edge[ e2 ] == 1 )
                                {
                                        numones ++ ;
                                }
                                else if ( edge[ e2 ] == 0 )
                                {
                                        numones -- ;
                                }
                        }

                        if ( tot > 0 )
                        {
                                return 1 ;
                        }

                        return 0 ;
                }
int Volume::isFeatureFace ( int  ox,
int  oy,
int  oz 
)

Definition at line 2135 of file volume.cpp.

References getDataAt(), getNumNeighbor6(), nx, ny, and wustl_mm::SkeletonMaker::sheetNeighbor.

Referenced by isSheetEnd().

                {
                        // return 1 ;

                        int i, j ;
                        int nx, ny, nz ;

                        int faces = 12 ;
                        for ( i = 0 ; i < 12 ; i ++ )
                        {
                                int ct = 0 ;
                                for ( j = 0 ; j < 4 ; j ++ )
                                {
                                        nx = ox + sheetNeighbor[i][j][0] ;
                                        ny = oy + sheetNeighbor[i][j][1] ;
                                        nz = oz + sheetNeighbor[i][j][2] ;

                                        if ( getDataAt( nx, ny, nz ) < 0 )
                                        {
                                                ct = -1 ;
                                                break ;
                                        }
                                        else if ( getNumNeighbor6( nx, ny, nz ) == 6 )
                                        {
                                                ct = -1 ;
                                                break ;
                                        }
        //                              else if ( getDataAt( nx, ny, nz ) == 0 )
        //                              {
        //                                      ct ++ ;
        //                              }


                                }
                                if ( ct == -1 || ct >= 1 )
                                {
                                        faces -- ;
                                }
                        }

                        if ( faces > 0 )
                        {
                                return 1 ;
                        }
                        return 0 ;

                }
int Volume::isHelixEnd ( int  ox,
int  oy,
int  oz 
)

Definition at line 1789 of file volume.cpp.

References getDataAt(), getNumNeighbor6(), wustl_mm::SkeletonMaker::neighbor6, nx, and ny.

                                                               {

                        int i ;
                        int c1 = 0 , c2 = 0 ;
                        int nx, ny, nz ;

                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                nx = ox + neighbor6[i][0] ;
                                ny = oy + neighbor6[i][1] ;
                                nz = oz + neighbor6[i][2] ;

                                double val = getDataAt( nx, ny, nz ) ;

                                if ( val >= 0 )
                                {
                                        c1 ++ ;
                                        if ( getNumNeighbor6(nx,ny,nz) < 6 ) // if ( val > 0 && val < MAX_ERODE )
                                        {
                                                c2 ++ ;
                                        }
                                }

                        }

                        if ( c1 == 1 && c2 == 1 )
                        {
                                return 1 ;
                        }

                        return 0 ;
                }
int Volume::isHelixEnd ( int  ox,
int  oy,
int  oz,
Volume nvol 
)

Definition at line 1579 of file volume.cpp.

References getDataAt(), wustl_mm::SkeletonMaker::neighbor6, nx, and ny.

Referenced by curveSkeleton(), curveSkeleton2D(), and surfaceSkeletonPres().

                                                                             {
                        // Returns 1 if it is a curve endpoint
                        int i ;
                        int c1 = 0 , c2 = 0 ;
                        int nx, ny, nz ;

                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                nx = ox + neighbor6[i][0] ;
                                ny = oy + neighbor6[i][1] ;
                                nz = oz + neighbor6[i][2] ;

                                double val = getDataAt( nx, ny, nz ) ;

                                if ( val >= 0 )
                                {
                                        c1 ++ ;
                                        if ( val > 0 && val < MAX_ERODE && nvol->getDataAt( nx, ny, nz ) == 0 )
                                        {
                                                c2 ++ ;
                                        }
                                }

                        }

                        if ( c1 == 1 && c2 == 1 )
                        {
                                return 1 ;
                        }

                        return 0 ;
                }
int Volume::isPiercable ( int  ox,
int  oy,
int  oz 
)

Definition at line 2381 of file volume.cpp.

References countExt(), countInt(), and getDataAt().

Referenced by curveSkeleton().

                {
                        /* Test if this is a simple voxel */
                        // int flag = 0 ;
                        double vox[3][3][3] ;

                        int i, j, k ;
                        for ( i = -1 ; i < 2 ; i ++ )
                                for ( j = -1 ; j < 2 ; j ++ )
                                        for ( k = -1 ; k < 2 ; k ++ )
                                        {
                                                double tval = getDataAt( ox + i, oy + j, oz + k ) ;

                                                /*
                                                if ( tval == 0 || tval > (va )
                                                {
                                                        flag = 1 ;
                                                }
                                                */
                                                /*
                                                if ( tval < 0 && tval == - curwid )
                                                {
                                                printf("Here %d", (int)-tval) ;
                                                vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
                                                }
                                                else
                                                */
                                                {
                                                        vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
                                                }
                                        }

                                /* Debugging
                                printf("{") ;
                                for ( i = 0 ; i < 3 ; i ++ )
                                {
                                        if ( i ) printf(",") ;
                                        printf("{") ;
                                        for ( j = 0 ; j < 3 ; j ++ )
                                        {
                                                if ( j ) printf(",") ;
                                                printf("{") ;
                                                for ( k = 0 ; k < 3 ; k ++ )
                                                {
                                                        if ( k ) printf(",") ;
                                                        printf("%d", (vox[i][j][k] >=0 ? 1: 0));
                                                }
                                                printf("}") ;
                                        }
                                        printf("}") ;
                                }
                                printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
                                */

                        if ( countInt( vox ) == 1 && countExt( vox ) != 1 )
                        {
                                return 1 ;
                        }
                        else
                        {
                                return 0 ;
                        }
                }
int Volume::isSheetEnd ( int  ox,
int  oy,
int  oz,
Volume nvol 
)

Definition at line 1821 of file volume.cpp.

References wustl_mm::SkeletonMaker::edgeFaces, wustl_mm::SkeletonMaker::faceEdges, getDataAt(), nx, ny, wustl_mm::SkeletonMaker::sheetNeighbor, and tot.

Referenced by surfaceSkeletonPres().

                {
                        // Returns 1 if it contains a sheet boundary. Noise-resistant
                        int i, j ;
                        int nx, ny, nz ;

                        int edge[6] = { 0,0,0,0,0,0 } ;
                        int faceflag[ 12 ] ;
                        int hasFeatureFace = 0 ;
                        int tot = 0 ;

                        for ( i = 0 ; i < 12 ; i ++ )
                        {
                                faceflag[ i ] = 1 ;
                                int hasFeature = 1 ;

                                for ( j = 0 ; j < 4 ; j ++ )
                                {
                                        nx = ox + sheetNeighbor[i][j][0] ;
                                        ny = oy + sheetNeighbor[i][j][1] ;
                                        nz = oz + sheetNeighbor[i][j][2] ;

                                        if ( getDataAt( nx, ny, nz ) == 0 || nvol->getDataAt( nx, ny, nz ) == 1 )
                                        {
                                                hasFeature = 0 ;
                                        }
                                        if ( getDataAt( nx, ny, nz ) < 0 )
                                        {
                                                faceflag[ i ] = 0 ;
                                                break ;
                                        }
                                }
                                if ( faceflag[ i ] == 1 && hasFeature )
                                {
                                        hasFeatureFace ++ ;
                                        // return 0 ;
                                }

                                if ( faceflag[ i ] )
                                {
                                        int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
                                        edge[ e0 ] ++ ;
                                        edge[ e1 ] ++ ;
                                        tot ++ ;
                                }
                        }

                        if ( tot == 0 || hasFeatureFace == 0 )
                        {
                                return 0 ;
                        }

                        // Removing 1s
                        int numones = 0 ;
                        for ( i = 0 ; i < 6 ; i ++ )
                        {
                                if ( edge[ i ] == 1 )
                                {
                                        numones ++ ;
                                }
                        }
                        while( numones > 0 )
                        {
                                int e = 0;
                                for ( i = 0 ; i < 6 ; i ++ )
                                {
                                        if ( edge[ i ] == 1 )
                                        {
                                                e = i ;
                                                break ;
                                        }
                                }
                                /*
                                if ( edge[ e ] != 1 )
                                {
                                        printf("Wrong Again!********\n") ;
                                }
                                */

                                int f, e2 ;
                                for ( j = 0 ; j < 4 ; j ++ )
                                {
                                        f = edgeFaces[ e ][ j ] ;
                                        if ( faceflag[ f ] )
                                        {
                                                break ;
                                        }
                                }

                                /*
                                if ( faceflag[ f ] == 0 )
                                {
                                        printf("Wrong!********\n") ;
                                }
                                */

                                if ( faceEdges[ f ][ 0 ] == e )
                                {
                                        e2 = faceEdges[ f ][ 1 ] ;
                                }
                                else
                                {
                                        e2 = faceEdges[ f ][ 0 ] ;
                                }


                                edge[ e ] -- ;
                                numones -- ;
                                edge[ e2 ] -- ;
                                faceflag[ f ] = 0 ;
                                tot -- ;

                                if ( edge[ e2 ] == 1 )
                                {
                                        numones ++ ;
                                }
                                else if ( edge[ e2 ] == 0 )
                                {
                                        numones -- ;
                                }
                        }

                        if ( tot > 0 )
                        {
                                return 0 ;
                        }

                        return 1 ;
                }
int Volume::isSheetEnd ( int  ox,
int  oy,
int  oz 
)

*

Definition at line 2238 of file volume.cpp.

References hasCompleteSheet(), and isFeatureFace().

                {
                        return ( ( hasCompleteSheet(ox,oy,oz) == 0 ) && isFeatureFace(ox,oy,oz) ) ;

                        //
                        //int i, j, k ;
                        //int nx, ny, nz ;

                        //double vox[3][3][3] ;
                        //for ( i = -1 ; i < 2 ; i ++ )
                        //      for ( j = -1 ; j < 2 ; j ++ )
                        //              for ( k = -1 ; k < 2 ; k ++ )
                        //              {
                        //                      vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ;
                        //              }

                        //int edge[6] = { 4,4,4,4,4,4 } ;
                        //int edge2[6] = { 4,4,4,4,4,4 } ;

                        //for ( i = 0 ; i < 12 ; i ++ )
                        //{
                        //      for ( j = 0 ; j < 4 ; j ++ )
                        //      {
                        //              nx = 1 + sheetNeighbor[i][j][0] ;
                        //              ny = 1 + sheetNeighbor[i][j][1] ;
                        //              nz = 1 + sheetNeighbor[i][j][2] ;

                        //              if ( vox[nx][ny][nz] < 0 )
                        //              {
                        //                      edge[ faceEdges[ i ][ 0 ] ] -- ;
                        //                      edge[ faceEdges[ i ][ 1 ] ] -- ;
                        //                      break ;
                        //              }
                        //      }

                        //      for ( j = 0 ; j < 4 ; j ++ )
                        //      {
                        //              nx = 1 + sheetNeighbor[i][j][0] ;
                        //              ny = 1 + sheetNeighbor[i][j][1] ;
                        //              nz = 1 + sheetNeighbor[i][j][2] ;

                        //              if ( vox[nx][ny][nz] <= 0 )
                        //              {
                        //                      edge2[ faceEdges[ i ][ 0 ] ] -- ;
                        //                      edge2[ faceEdges[ i ][ 1 ] ] -- ;
                        //                      break ;
                        //              }
                        //      }
                        //}

                        //
                        //for ( i = 0 ; i < 6 ; i ++ )
                        //{
                        //      nx = 1 + neighbor6[i][0] ;
                        //      ny = 1 + neighbor6[i][1] ;
                        //      nz = 1 + neighbor6[i][2] ;
                        //      if ( edge[i] == 0 && vox[nx][ny][nz] >= 0 )
                        //      {
                        //              return 0 ;
                        //      }
                        //}
                        //*/
                        //


                        //for ( i = 0 ; i < 6 ; i ++ )
                        //{
                        //      if ( edge[ i ] == 1 && edge2[ i ] == 1 )
                        //      {
                        //              return 1 ;
                        //      }
                        //}

                        //return 0 ;
                }
int Volume::isSimple ( int  ox,
int  oy,
int  oz 
)

Definition at line 2316 of file volume.cpp.

References countExt(), countInt(), and getDataAt().

Referenced by curveSkeleton(), curveSkeleton2D(), skeleton(), and surfaceSkeletonPres().

                {
                        /* Test if this is a simple voxel */
                        // int flag = 0 ;
                        double vox[3][3][3] ;

                        int i, j, k ;
                        for ( i = -1 ; i < 2 ; i ++ )
                                for ( j = -1 ; j < 2 ; j ++ )
                                        for ( k = -1 ; k < 2 ; k ++ )
                                        {
                                                double tval = getDataAt( ox + i, oy + j, oz + k ) ;

                                                /*
                                                if ( tval == 0 || tval > (va )
                                                {
                                                        flag = 1 ;
                                                }
                                                */
                                                /*
                                                if ( tval < 0 && tval == - curwid )
                                                {
                                                printf("Here %d", (int)-tval) ;
                                                vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
                                                }
                                                else
                                                */
                                                {
                                                        vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
                                                }
                                        }

                                /* Debugging
                                printf("{") ;
                                for ( i = 0 ; i < 3 ; i ++ )
                                {
                                        if ( i ) printf(",") ;
                                        printf("{") ;
                                        for ( j = 0 ; j < 3 ; j ++ )
                                        {
                                                if ( j ) printf(",") ;
                                                printf("{") ;
                                                for ( k = 0 ; k < 3 ; k ++ )
                                                {
                                                        if ( k ) printf(",") ;
                                                        printf("%d", (vox[i][j][k] >=0 ? 1: 0));
                                                }
                                                printf("}") ;
                                        }
                                        printf("}") ;
                                }
                                printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
                                */

                        if ( countInt( vox ) == 1 && countExt( vox ) == 1 )
                        {
                                return 1 ;
                        }
                        else
                        {
                                return 0 ;
                        }
                }
Volume * Volume::markCellFace ( )

Definition at line 1027 of file volume.cpp.

References getDataAt(), getSizeX(), getSizeY(), getSizeZ(), hasCell(), wustl_mm::SkeletonMaker::neighbor6, nx, ny, setDataAt(), and Volume().

Referenced by erodeSheet().

                                               {
                        int i, j, k ;
                        Volume* fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;

                        //return fvol ;

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if( getDataAt(i,j,k) >= 0 )
                                                {
                                                        if ( hasCell(i,j,k) )
                                                        {
                                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                                {
                                                                        int nx = i + neighbor6[m][0] ;
                                                                        int ny = j + neighbor6[m][1] ;
                                                                        int nz = k + neighbor6[m][2] ;
                                                                        if ( ! hasCell(nx,ny,nz) )
                                                                        {
                                                                                fvol->setDataAt(i,j,k,(double)(1<<m)) ;
                                                                                break ;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }


                        return fvol ;
                }
void Volume::pad ( int  padBy,
double  padValue 
)

* Next, smooth

* Next, smooth */

Definition at line 273 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::Pad(), and volData.

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PerformPureJuSkeletonization().

                                                            {
                        volData->Pad(padBy, padValue);
                }
void Volume::setDataAt ( int  x,
int  y,
int  z,
double  d 
)
void Volume::setDataAt ( int  index,
double  d 
)

Definition at line 56 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::SetDataAt(), and volData.

                                                            {
                        volData->SetDataAt(index, (float)d);
                }
void Volume::setOrigin ( float  orgX,
float  orgY,
float  orgZ 
)

Definition at line 77 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::SetOrigin(), and volData.

                                                                         {
                        volData->SetOrigin(orgX, orgY, orgZ);
                }
void Volume::setSpacing ( float  spx,
float  spy,
float  spz 
)

Definition at line 73 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::SetSpacing(), and volData.

                                                                        {
                        volData->SetSpacing(spx, spy, spz);
                }
void Volume::skeleton ( float  thr,
int  off 
)

Definition at line 5089 of file volume.cpp.

References getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getSizeX(), getSizeY(), getSizeZ(), isSimple(), wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), setDataAt(), threshold(), wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridQueueEle::y, and wustl_mm::SkeletonMaker::gridQueueEle::z.

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning().

                {
                        int i, j, k ;
                        // First, threshold the volume
                        #ifdef VERBOSE
                        printf("Thresholding the volume to -1/0...\n") ;
                        #endif
                        threshold( thr, -1, 0 ) ;

                        // Next, apply convergent erosion
                        // by preserving: complex nodes, curve end-points, and sheet points

                        // Next, initialize the linked queue
                        #ifdef VERBOSE
                        printf("Initializing queue...\n") ;
                        #endif
                        GridQueue2* queue2 = new GridQueue2( ) ;
                        GridQueue2* queue = new GridQueue2( ) ;

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) >= 0 )
                                                {
                                                        {
                                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                                {
                                                                        if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
                                                                        {
                                                                                setDataAt( i, j, k, 1 ) ;
                                                                                queue2->prepend( i, j, k ) ;
                                                                                break ;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                        int wid = 0 ;
                        #ifdef VERBOSE
                        printf("Total %d nodes\n", queue2->getNumElements() ) ;
                        printf("Start erosion to %d...\n", wid) ;
                        #endif

                        // Perform erosion
                        gridQueueEle* ele ;
                        int ox, oy, oz ;

                        while( 1 )
                        {
                                wid ++ ;

                                // At the start of each iteration,
                                // queue2 holds all the nodes for this layer
                                // queue is empty

                                int numComplex = 0, numSimple = 0 ;
                                #ifdef VERBOSE
                                printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), wid) ;
                                #endif

                                // Rename queue2 to be queue,
                                // Clear queue2
                                // From now on, queue2 holds nodes of next level
                                delete queue ;
                                queue = queue2 ;
                                queue2 = new GridQueue2( ) ;

                                // Next, start queue iteration
                                queue->reset() ;
                                ele = queue->getNext();
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;
        //                              delete ele ;

                                        // Check simple
                                        if ( ! isSimple( ox, oy, oz ) )
                                        {
                                                // Complex, set to next layer
                                                queue2->prepend( ox, oy, oz ) ;
                                                numComplex ++ ;
                                        }
                                        /*
                                        else if ( ox == off || oy == off || oz == off ||
                                                ox == getSizeX() - off - 1 || oy == getSizeY() - off - 1 || oz == getSizeZ() - off - 1 )
                                        {
                                                // Wall, don't erode, set to next layer
                                                queue2->prepend( ox, oy, oz ) ;
                                                numComplex ++ ;
                                        }
                                        */
                                        else
                                        {
                                                setDataAt( ox, oy, oz, -1 ) ;
                                                numSimple ++ ;

                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) == 0 )
                                                        {
                                                                setDataAt( nx, ny, nz, 1 ) ;
                                                                queue2->prepend( nx, ny, nz ) ;
                                                        }
                                                }

                                        }

                                        ele = queue->remove() ;
                                }
                                #ifdef VERBOSE
                                printf("Level %d: %d complex, %d simple\n", wid, numComplex, numSimple) ;
                                #endif

                                if ( numSimple == 0 )
                                {
                                        break ;
                                }
                        }

                        // Finally, clean up
                        delete queue;
                        delete queue2;
                        #ifdef VERBOSE
                        printf("Thresholding the volume to 0/1...\n") ;
                        #endif
                        threshold( 0, 0, 1 ) ;
                }
void Volume::skeleton ( float  thr,
Volume svol,
Volume hvol 
)

* Added for debugging */ / Check simple

* Adding ends */ insert them back into priority queue * * Thin the current volume while preserving voxels with values > highthr or <= lowthr in grayvol * Added for debugging */ / Check simple * Adding ends */

Definition at line 5692 of file volume.cpp.

References EMAN::PriorityQueue< ValueT, KeyT >::add(), getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isSimple(), MAX_ERODE, MAX_QUEUELEN, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z.

                {
                        int i, j, k ;
                        // First, threshold the volume
                        #ifdef VERBOSE
                        printf("Thresholding the volume to -1/0...\n") ;
                        #endif
                        threshold( thr, -1, 0 ) ;

                        // Next, apply convergent erosion
                        // by preserving: complex nodes, curve end-points, and sheet points

                        // Next, initialize the linked queue
                        #ifdef VERBOSE
                        printf("Initializing queue...\n") ;
                        #endif
                        GridQueue2* queue2 = new GridQueue2( ) ;
                        GridQueue2* queue3 = new GridQueue2( ) ;
                        PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) >= 0 )
                                                {
                                                        if ( svol->getDataAt(i,j,k) > 0 || hvol->getDataAt(i,j,k) > 0 )
                                                        {
                                                                setDataAt( i, j, k, MAX_ERODE ) ;
                                                        }
                                                        else
                                                        {
                                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                                {
                                                                        if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
                                                                        {
                                                                                setDataAt( i, j, k, 1 ) ;
                                                                                queue2->prepend( i, j, k ) ;
                                                                                break ;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                        int wid = MAX_ERODE ;
                        #ifdef VERBOSE
                        printf("Total %d nodes\n", queue2->getNumElements() ) ;


                        // Perform erosion
                        printf("Start erosion to %d...\n", wid) ;
                        #endif
                        gridQueueEle* ele ;
                        gridPoint* gp ;
                        int ox, oy, oz ;
                        int score ;
                        Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
                        for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
                        {
                                scrvol->setDataAt( i, -1 ) ;
                        }


                        for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
                        {
                                // At the start of each iteration,
                                // queue2 holds all the nodes for this layer
                                // queue3 and queue are empty

                                int numComplex = 0, numSimple = 0 ;
                                #ifdef VERBOSE
                                printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
                                #endif


                                // Next,
                                // Compute score for each node left in queue2
                                // move them into priority queue
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        // Compute score
                                        score = getNumPotComplex2( ox, oy, oz ) ;
                                        scrvol->setDataAt( ox, oy, oz, score ) ;

                                        // Push to queue
                                        gp = new gridPoint ;
                                        gp->x = ox ;
                                        gp->y = oy ;
                                        gp->z = oz ;
                                        // queue->add( gp, -score ) ;
                                        queue->add( gp, score ) ;

                                        ele = queue2->remove() ;
                                }

                                // Rename queue3 to be queue2,
                                // Clear queue3
                                // From now on, queue2 holds nodes of next level
                                delete queue2 ;
                                queue2 = queue3 ;
                                queue3 = new GridQueue2( ) ;

                                // Next, start priority queue iteration
                                while ( ! queue->isEmpty() )
                                {
                                        // Retrieve the node with the highest score
                                        queue->remove( gp, score ) ;
                                        ox = gp->x ;
                                        oy = gp->y ;
                                        oz = gp->z ;
                                        delete gp ;
        //                              score = -score ;

                                        // Ignore the node
                                        // if it has been processed before
                                        // or it has an updated score
                                        if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
                                        {
                                                continue ;
                                        }

                                        /* Added for debugging */
                                        // Check simple
                                        if ( ! isSimple( ox, oy, oz ) )
                                        {
                                                // Complex, set to next layer
                                                setDataAt( ox, oy, oz, curwid + 1 ) ;
                                                queue2->prepend( ox, oy, oz ) ;
                                                numComplex ++ ;
                                        }
                                        else
                                        {
                                                setDataAt( ox, oy, oz, -1 ) ;
                                                numSimple ++ ;
                                        }
                                        /* Adding ends */

                                        // Move its neighboring unvisited node to queue2
                                        for ( int m = 0 ; m < 6 ; m ++ )
                                        {
                                                int nx = ox + neighbor6[m][0] ;
                                                int ny = oy + neighbor6[m][1] ;
                                                int nz = oz + neighbor6[m][2] ;
                                                if ( getDataAt( nx, ny, nz ) == 0 )
                                                {
                                                        setDataAt( nx, ny, nz, curwid + 1 ) ;
                                                        queue2->prepend( nx, ny, nz ) ;
                                                }
                                        }

                                        // Update scores for nodes in its 5x5 neighborhood
                                        // insert them back into priority queue

                                        for ( i = -2 ; i < 3 ;i ++ )
                                                for ( j = -2 ; j < 3 ; j ++ )
                                                        for ( k = -2 ; k < 3 ; k ++ )
                                                        {
                                                                int nx = ox + i ;
                                                                int ny = oy + j ;
                                                                int nz = oz + k ;

                                                                if ( getDataAt( nx, ny, nz ) == curwid )
                                                                {
                                                                        // Compute score
                                                                        score = getNumPotComplex2( nx, ny, nz ) ;

                                                                        if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
                                                                        {
                                                                                // printf("Update\n") ;
                                                                                scrvol->setDataAt( nx, ny, nz, score ) ;
                                                                                // Push to queue
                                                                                gp = new gridPoint ;
                                                                                gp->x = nx ;
                                                                                gp->y = ny ;
                                                                                gp->z = nz ;
                                                                                // queue->add( gp, -score ) ;
                                                                                queue->add( gp, score ) ;
                                                                        }
                                                                }
                                                        }


                                }

                                #ifdef VERBOSE
                                printf("%d complex, %d simple\n", numComplex, numSimple) ;
                                #endif

                                if ( numSimple == 0 )
                                {
                                        delete queue2 ;
                                        break ;
                                }
                        }

                        // Finally, clean up
                        #ifdef VERBOSE
                        printf("Thresholding the volume to 0/1...\n") ;
                        #endif
                        threshold( 0, 0, 1 ) ;
                        delete scrvol;
                        delete queue;
                        delete queue3;
                }
*void Volume::surfaceSkeletonPres ( float  thr,
Volume preserve 
)

* Deal with closed rings

* Find local minimum: to help determine erosion level */ for ( int m = 0 ; m < 6 ; m ++ )* ************************************************************************/ / A sequential thinning method for extracting 6-connected skeletons / Properties: medial, topology preserving, shape preserving, noise resistance! /

Parameters:
thr,:threshold /
type,:0 for curve preserving, 1 for surface preserving
noise,:0 for no-noise-reduction, n for level-n noise reduction* Thin the current volume while preserving voxels with values > highthr or <= lowthr in grayvol * Extra step: classify nodes in queue2 into noise and non-noise nodes */ * * Extra step: classify nodes in queue2 into noise and non-noise nodes */ * Commented for debugging * Commented for debugging * Added for debugging */ / Check simple * Adding ends */ * Commented for debugging * * Extra step: classify nodes in queue2 into noise and non-noise nodes */ * Commented for debugging * Commented for debugging * Added for debugging */ / Check simple * Adding ends */ * Commented for debugging * * Extra step: classify nodes in queue2 into noise and non-noise nodes */ * Commented for debugging * Commented for debugging * Added for debugging */ / Check simple * Adding ends */ * Commented for debugging

Definition at line 8672 of file volume.cpp.

References EMAN::PriorityQueue< ValueT, KeyT >::add(), getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isSheetEnd(), isSimple(), MAX_ERODE, MAX_QUEUELEN, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z.

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning().

                {
                        int i, j, k ;
                        // First, threshold the volume
                        #ifdef VERBOSE
                        printf("Thresholding the volume to -MAX_ERODE/0...\n") ;
                        #endif
                        threshold( thr, -MAX_ERODE, 0 ) ;

                        // Next, initialize the linked queue
                        #ifdef VERBOSE
                        printf("Initializing queue...\n") ;
                        #endif
                        GridQueue2* queue2 = new GridQueue2( ) ;
                        GridQueue2* queue3 = new GridQueue2( ) ;
                        GridQueue2* queue4 = new GridQueue2( ) ;

                        PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );

                        for ( i = 0 ; i < getSizeX() ; i ++ )
                                for ( j = 0 ; j < getSizeY() ; j ++ )
                                        for ( k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                if ( getDataAt( i, j, k ) >= 0 ) {
                                                        if(preserve->getDataAt(i, j, k) > 0) {
                                                                setDataAt(i, j, k, MAX_ERODE);
                                                        } else {
                                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                                {
                                                                        if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
                                                                        {
                                                                                // setDataAt( i, j, k, 1 ) ;
                                                                                queue2->prepend( i, j, k ) ;
                                                                                break ;
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                        int wid = MAX_ERODE ;
                        #ifdef VERBOSE
                        printf("Total %d nodes\n", queue2->getNumElements() ) ;
                        printf("Start erosion to %d...\n", wid) ;
                        #endif


                        // Perform erosion
                        gridQueueEle* ele ;
                        gridPoint* gp ;
                        int ox, oy, oz ;
                        int score ;
                        Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
                        for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
                        {
                                scrvol->setDataAt( i, -1 ) ;
                        }

        #ifdef  NOISE_DIS_SHEET
                        Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
        #endif

                        for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
                        {
                                // At the start of each iteration,
                                // queue2 and queue4 holds all the nodes for this layer
                                // queue3 and queue are empty

                                int numComplex = 0, numSimple = 0 ;
                                #ifdef VERBOSE
                                printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
                                #endif

                                /*
                                We first need to assign curwid + 1 to every node in this layer
                                */
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        if ( getDataAt(ox,oy,oz) == curwid )
                                        {
                                                ele = queue2->remove() ;
                                        }
                                        else
                                        {
                                                setDataAt(ox,oy,oz, curwid) ;
                                                ele = queue2->getNext() ;
                                        }
                                }
                                queue4->reset() ;
                                ele = queue4->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        queue2->prepend(ox,oy,oz) ;
                                        ele = queue4->remove() ;
                                }

                                // Now queue2 holds all the nodes for this layer

        #ifdef NOISE_DIS_SHEET
                                /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
                                queue2->reset() ;

                                // First run
                                int flag = 0 ;
                                while ( ( ele = queue2->getNext() ) != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;
                                        if ( NOISE_DIS_SHEET <= 1 )
                                        {
                                                noisevol->setDataAt( ox, oy, oz, 0 ) ;
                                        }
                                        else
                                        {
                                                flag = 0 ;
                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) == 0 )
                                                        {
                                                                noisevol->setDataAt( ox, oy, oz, 1 ) ;
                                                                flag = 1 ;
                                                                break ;
                                                        }
                                                }
                                                if ( ! flag )
                                                {
                                                        noisevol->setDataAt( ox, oy, oz, 0 ) ;
                                                }
                                        }
                                }

                                for ( int cur = 1 ; cur < NOISE_DIS_SHEET ; cur ++ )
                                {
                                        queue2->reset() ;
                                        int count = 0 ;

                                        while ( ( ele = queue2->getNext() ) != NULL )
                                        {
                                                ox = ele->x ;
                                                oy = ele->y ;
                                                oz = ele->z ;

                                                if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
                                                {
                                                        continue ;
                                                }

                                                flag = 0 ;
                                                for ( int m = 0 ; m < 6 ; m ++ )
                                                {
                                                        int nx = ox + neighbor6[m][0] ;
                                                        int ny = oy + neighbor6[m][1] ;
                                                        int nz = oz + neighbor6[m][2] ;
                                                        if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
                                                        {
                                                                noisevol->setDataAt( ox, oy, oz, 1 ) ;
                                                                count ++ ;
                                                                break ;
                                                        }
                                                }
                                        }

                                        if ( count == 0 )
                                        {
                                                break ;
                                        }
                                }


        #endif

                                /* Commented for debugging

                                // First,
                                // check for complex nodes in queue2
                                // move them from queue2 to queue3
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        // Check simple
        #ifndef NOISE_DIS_SHEET
                                        if ( isSheetEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
        #else
                                        if ( isSheetEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
        #endif
                                        {
                                                // Complex, set to next layer
                                                setDataAt( ox, oy, oz, curwid + 1 ) ;
                                                queue3->prepend( ox, oy, oz ) ;
                                                ele = queue2->remove() ;

                                                numComplex ++ ;
                                        }
                                        else
                                        {
                                                ele = queue2->getNext() ;
                                        }
                                }
                                */


                                // Next,
                                // Compute score for each node left in queue2
                                // move them into priority queue
                                queue2->reset() ;
                                ele = queue2->getNext() ;
                                while ( ele != NULL )
                                {
                                        ox = ele->x ;
                                        oy = ele->y ;
                                        oz = ele->z ;

                                        // Compute score
                                        score = getNumPotComplex( ox, oy, oz ) ;
                                        scrvol->setDataAt( ox, oy, oz, score ) ;

                                        // Push to queue
                                        gp = new gridPoint ;
                                        gp->x = ox ;
                                        gp->y = oy ;
                                        gp->z = oz ;
                                        // queue->add( gp, -score ) ;
                                        queue->add( gp, score ) ;

                                        ele = queue2->remove() ;
                                }

                                // Rename queue3 to be queue2,
                                // Clear queue3
                                // From now on, queue2 holds nodes of next level
                                delete queue2 ;
                                queue2 = queue3 ;
                                queue3 = new GridQueue2( ) ;


                                // Next, start priority queue iteration
                                while ( ! queue->isEmpty() )
                                {
                                        // Retrieve the node with the highest score
                                        queue->remove( gp, score ) ;
                                        ox = gp->x ;
                                        oy = gp->y ;
                                        oz = gp->z ;
                                        delete gp ;
                                        // printf("%d\n", score);
        //                              score = -score ;

                                        // Ignore the node
                                        // if it has been processed before
                                        // or it has an updated score
                                        if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
                                        {
                                                continue ;
                                        }

                                        /* Commented for debugging

                                        // Remove this simple node
                                        setDataAt( ox, oy, oz, -1 ) ;
                                        numSimple ++ ;
                                        // printf("Highest score: %d\n", score) ;
                                        */

                                        /* Added for debugging */
                                        // Check simple
        #ifndef NOISE_DIS_SHEET
                                        // if ( hasFeatureFace( ox, oy, oz ) )
                                        if ( (! isSimple( ox, oy, oz )) || isSheetEnd( ox, oy, oz ) )
                                        // if ( hasIsolatedFace(ox,oy,oz)  && (! isNoiseSheetEnd(ox,oy,oz)))
        #else
                                        // if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz, noisevol ) )
                                        if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz, noisevol ) || isHelixEnd( ox, oy, oz, noisevol ))
                                        // if ( isBertrandEndPoint( ox, oy, oz ) )
        #endif
                                        {
                                                // Complex, set to next layer
                                                setDataAt( ox, oy, oz, curwid + 1 ) ;
                                                queue4->prepend( ox, oy, oz ) ;
                                                numComplex ++ ;

                                        }
                                        else
                                        {
                                                setDataAt( ox, oy, oz, -1 ) ;
                                                numSimple ++ ;

                                        }
                                        /* Adding ends */

                                        // Move its neighboring unvisited node to queue2
                                        for ( int m = 0 ; m < 6 ; m ++ )
                                        {
                                                int nx = ox + neighbor6[m][0] ;
                                                int ny = oy + neighbor6[m][1] ;
                                                int nz = oz + neighbor6[m][2] ;
                                                if ( getDataAt( nx, ny, nz ) == 0 )
                                                {
                                                        // setDataAt( nx, ny, nz, curwid + 1 ) ;
                                                        queue2->prepend( nx, ny, nz ) ;
                                                }
                                        }

                                        /* Commented for debugging

                                        // Find complex nodes in its 3x3 neighborhood
                                        // move them to queue2
                                        for ( i = -1 ; i < 2 ; i ++ )
                                                for ( j = -1 ; j < 2 ; j ++ )
                                                        for ( k = -1 ; k < 2 ; k ++ )
                                                        {
                                                                int nx = ox + i ;
                                                                int ny = oy + j ;
                                                                int nz = oz + k ;

                                                                // Check simple
                                                                if ( getDataAt( nx, ny, nz ) == curwid &&
                                                                        // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
        #ifndef NOISE_DIS_SHEET
                                                                        ( isSheetEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
        #else
                                                                        ( isSheetEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
        #endif

                                                                {
                                                                        // Complex, set to next layer
                                                                        setDataAt( nx, ny, nz, curwid + 1 ) ;
                                                                        queue2->prepend( nx, ny, nz ) ;
                                                                        numComplex ++ ;
                                                                }
                                                        }
                                        */

                                        // Update scores for nodes in its 5x5 neighborhood
                                        // insert them back into priority queue

                                        for ( i = -2 ; i < 3 ;i ++ )
                                                for ( j = -2 ; j < 3 ; j ++ )
                                                        for ( k = -2 ; k < 3 ; k ++ )
                                                        {
                                                                int nx = ox + i ;
                                                                int ny = oy + j ;
                                                                int nz = oz + k ;

                                                                if ( getDataAt( nx, ny, nz ) == curwid )
                                                                {
                                                                        // Compute score
                                                                        score = getNumPotComplex( nx, ny, nz ) ;

                                                                        if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
                                                                        {
                                                                                // printf("Update\n") ;
                                                                                scrvol->setDataAt( nx, ny, nz, score ) ;
                                                                                // Push to queue
                                                                                gp = new gridPoint ;
                                                                                gp->x = nx ;
                                                                                gp->y = ny ;
                                                                                gp->z = nz ;
                                                                                // queue->add( gp, -score ) ;
                                                                                queue->add( gp, score ) ;
                                                                        }
                                                                }
                                                        }


                                }

                                #ifdef VERBOSE
                                printf("%d complex, %d simple\n", numComplex, numSimple) ;
                                #endif

                                if ( numSimple == 0 )
                                {
                                                break ;
                                }
                        }

                        // Finally, clean up
                        #ifdef VERBOSE
                        printf("Thresholding the volume to 0/1...\n") ;
                        #endif
                        threshold( 0, 0, 1 ) ;

                        delete scrvol;
                        delete queue;
                        delete queue2;
                        delete queue3;
                        delete queue4;

                }
void Volume::threshold ( double  thr)

* Bertrand's parallel 6-connected surface thinning */

* Palagyi's parallel surface thinning */ * * Debug */ Normalize to a given range

Definition at line 9515 of file volume.cpp.

Referenced by curveSkeleton(), curveSkeleton2D(), erodeHelix(), erodeSheet(), skeleton(), surfaceSkeletonPres(), and threshold().

                {
                        threshold( thr, 0, 1, 0, true) ;
                }
void Volume::threshold ( double  thr,
int  out,
int  in 
)

Definition at line 9520 of file volume.cpp.

References threshold().

                {
                        threshold( thr, out, in, out, true) ;
                }
void Volume::threshold ( double  thr,
int  out,
int  in,
int  boundary 
)

Definition at line 9525 of file volume.cpp.

References threshold().

                {
                        threshold(thr, out, in, boundary, true);
                }
void Volume::threshold ( double  thr,
int  out,
int  in,
int  boundary,
bool  markBoundary 
)

Definition at line 9530 of file volume.cpp.

References getDataAt(), getSizeX(), getSizeY(), getSizeZ(), and setDataAt().

                {
                        float val;
                        for ( int i = 0 ; i < getSizeX() ; i ++ )
                                for ( int j = 0 ; j < getSizeY() ; j ++ )
                                        for ( int k = 0 ; k < getSizeZ() ; k ++ )
                                        {
                                                val = (float)getDataAt(i, j, k);
                                                if(markBoundary) {
                                                        if ( i > 1 && i < getSizeX() - 2 && j > 1 && j < getSizeY() - 2 && k > 1 && k < getSizeZ() - 2 ) {
                                                                if(val < thr) {
                                                                        setDataAt(i, j, k, out);
                                                                } else {
                                                                        setDataAt(i, j, k, in);
                                                                }
                                                        }
                                                        else
                                                        {
                                                                setDataAt(i, j, k, boundary);
                                                        }
                                                } else {
                                                        if(val < thr) {
                                                                setDataAt(i, j, k, out);
                                                        } else {
                                                                setDataAt(i, j, k, in);
                                                        }
                                                }
                                        }
                }

Member Data Documentation


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