Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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)
 * Thin the current volume while preserving voxels with values > highthr or <= lowthr in grayvol
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)
 * Adding ends *// Move its neighboring unvisited node to queue2
void erodeHelix ()
void erodeHelix (int disthr)
int erodeSheet ()
int erodeSheet (int disthr)
void surfaceSkeletonPres (float thr, Volume *preserve)
 * Commented for debugging/ Find complex nodes in its 3x3 neighborhood/ move them to queue2
void threshold (double thr)
 Normalize to a given range.
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

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().

00010                 {
00011                         this->volData = new VolumeData(em);
00012                 }

Volume::Volume int  x,
int  y,
int  z
 

Definition at line 13 of file volume.cpp.

References volData, x, and y.

00013                                                   {
00014                         volData = new VolumeData(x, y, z);
00015                 }

Volume::Volume int  x,
int  y,
int  z,
float  val
 

Definition at line 17 of file volume.cpp.

References volData, x, and y.

00017                                                              {
00018                         volData = new VolumeData(x, y, z, val);
00019                 }

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(), volData, x, and y.

00021                                                                                               {
00022                         volData = new VolumeData(x, y, z, offx, offy, offz, vol->getVolumeData());
00023                 }

Volume::~Volume  ) 
 

Definition at line 26 of file volume.cpp.

00027                 {
00028                         delete volData;
00029                 }


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().

02768                 {
02769                         // Stupid flood fill
02770                         int tot = 0 ;
02771                         int queue[27][3] ;
02772                         int vis[3][3][3] ;
02773                         int head = 0, tail = 1 ;
02774                         int i, j, k ;
02775                         for ( i = 0 ; i < 3 ; i ++ )
02776                                 for ( j = 0 ; j < 3 ; j ++ )
02777                                         for ( k = 0 ; k < 3 ; k ++ )
02778                                         {
02779                                                 vis[i][j][k] = 0 ;
02780                                                 if ( vox[i][j][k] )
02781                                                 {
02782                                                         if ( tot == 0 )
02783                                                         {
02784                                                                 queue[0][0] = i ;
02785                                                                 queue[0][1] = j ;
02786                                                                 queue[0][2] = k ;
02787                                                                 vis[i][j][k] = 1 ;
02788                                                         }
02789                                                         tot ++ ;
02790                                                 }
02791                                         }
02792                         if ( tot == 0 )
02793                         {
02794                                 return 0 ;
02795                         }
02796 
02797                         int ct = 1 ;
02798                         while ( head != tail )
02799                         {
02800                                 int x = queue[head][0] ;
02801                                 int y = queue[head][1] ;
02802                                 int z = queue[head][2] ;
02803                                 head ++ ;
02804 
02805                                 for ( i = -1 ; i < 2 ; i ++ )
02806                                         for ( j = -1 ; j < 2 ; j ++ )
02807                                                 for ( k = -1 ; k < 2 ; k ++ )
02808                                                 {
02809                                                         int nx = x + i ;
02810                                                         int ny = y + j ;
02811                                                         int nz = z + k ;
02812                                                         if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
02813                                                         {
02814                                                                 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
02815                                                                 {
02816                                                                         queue[tail][0] = nx ;
02817                                                                         queue[tail][1] = ny ;
02818                                                                         queue[tail][2] = nz ;
02819                                                                         tail ++ ;
02820                                                                         vis[nx][ny][nz] = 1 ;
02821                                                                         ct ++ ;
02822                                                                 }
02823                                                         }
02824                                                 }
02825                         }
02826 
02827                         if ( ct == tot )
02828                         {
02829                                 return 1 ;
02830                         }
02831                         else
02832                         {
02833                                 return 2 ;
02834                         }
02835 
02836                 }

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().

02699                 {
02700                         // Stupid flood fill
02701                         int tot = 0 ;
02702                         int queue[27][3] ;
02703                         int vis[3][3][3] ;
02704                         int head = 0, tail = 1 ;
02705                         int i, j, k ;
02706                         for ( i = 0 ; i < 3 ; i ++ )
02707                                 for ( j = 0 ; j < 3 ; j ++ )
02708                                         for ( k = 0 ; k < 3 ; k ++ )
02709                                         {
02710                                                 vis[i][j][k] = 0 ;
02711                                                 if ( vox[i][j][k] )
02712                                                 {
02713                                                         if ( tot == 0 )
02714                                                         {
02715                                                                 queue[0][0] = i ;
02716                                                                 queue[0][1] = j ;
02717                                                                 queue[0][2] = k ;
02718                                                                 vis[i][j][k] = 1 ;
02719                                                         }
02720                                                         tot ++ ;
02721                                                 }
02722                                         }
02723                         if ( tot == 0 )
02724                         {
02725                                 return 0 ;
02726                         }
02727                         // printf("total: %d\n", tot) ;
02728 
02729                         int ct = 1 ;
02730                         while ( head != tail )
02731                         {
02732                                 int x = queue[head][0] ;
02733                                 int y = queue[head][1] ;
02734                                 int z = queue[head][2] ;
02735                                 head ++ ;
02736 
02737                                 for ( i = 0 ; i < 6 ; i ++ )
02738                                 {
02739                                         int nx = x + neighbor6[i][0] ;
02740                                         int ny = y + neighbor6[i][1] ;
02741                                         int nz = z + neighbor6[i][2] ;
02742                                         if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
02743                                         {
02744                                                 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
02745                                                 {
02746                                                         queue[tail][0] = nx ;
02747                                                         queue[tail][1] = ny ;
02748                                                         queue[tail][2] = nz ;
02749                                                         tail ++ ;
02750                                                         vis[nx][ny][nz] = 1 ;
02751                                                         ct ++ ;
02752                                                 }
02753                                         }
02754                                 }
02755                         }
02756 
02757                         if ( ct == tot )
02758                         {
02759                                 return 1 ;
02760                         }
02761                         else
02762                         {
02763                                 return 2 ;
02764                         }
02765 
02766                 }

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

Definition at line 2838 of file volume.cpp.

References components26().

Referenced by isPiercable(), and isSimple().

02839                 {
02840                         int tvox[3][3][3] ;
02841 
02842                         for ( int i = 0 ; i < 3 ; i ++ )
02843                                 for ( int j = 0 ; j < 3 ; j ++ )
02844                                         for ( int k = 0 ; k < 3 ; k ++ )
02845                                         {
02846                                                 if ( vox[i][j][k] < 0 )
02847                                                 {
02848                                                         tvox[i][j][k] = 1 ;
02849                                                 }
02850                                                 else
02851                                                 {
02852                                                         tvox[i][j][k] = 0 ;
02853                                                 }
02854                                         }
02855 
02856                         return components26( tvox ) ;
02857                 }

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().

02860                 {
02861                         int i, j, k ;
02862                         int tvox[3][3][3] ;
02863 
02864                         for ( i = 0 ; i < 3 ; i ++ )
02865                                 for ( j = 0 ; j < 3 ; j ++ )
02866                                         for ( k = 0 ; k < 3 ; k ++ )
02867                                         {
02868                                                 tvox[i][j][k] = 0 ;
02869                                         }
02870 
02871                         for ( i = 0 ; i < 6 ; i ++ )
02872                         {
02873                                 int nx = 1 + neighbor6[i][0] ;
02874                                 int ny = 1 + neighbor6[i][1] ;
02875                                 int nz = 1 + neighbor6[i][2] ;
02876                                 if ( vox[ nx ][ ny ][ nz ] >= 0 )
02877                                 {
02878                                         tvox[ nx ][ ny ][ nz ] = 1 ;
02879                                         for ( j = 0 ; j < 4 ; j ++ )
02880                                         {
02881                                                 int nnx = neighbor64[i][j][0] + nx ;
02882                                                 int nny = neighbor64[i][j][1] + ny ;
02883                                                 int nnz = neighbor64[i][j][2] + nz ;
02884                                                 if ( vox[ nnx ][ nny ][ nnz ] >= 0 )
02885                                                 {
02886                                                         tvox[ nnx ][ nny ][ nnz ] = 1 ;
02887                                                 }
02888                                         }
02889                                 }
02890                         }
02891 
02892                         return components6( tvox ) ;
02893                 }

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().

02896                 {
02897                         int nv = 0 , ne = 0, nc = 0 ;
02898 
02899                         int i, j, k ;
02900                         int tvox[3][3][3] ;
02901                         double vox[3][3][3] ;
02902 
02903                         for ( i = 0 ; i < 3 ; i ++ )
02904                                 for ( j = 0 ; j < 3 ; j ++ )
02905                                         for ( k = 0 ; k < 3 ; k ++ )
02906                                         {
02907                                                 vox[i][j][k] = getDataAt( ox - 1 + i, oy - 1 + j, oz - 1 + k ) ;
02908                                                 tvox[i][j][k] = 0 ;
02909                                         }
02910 
02911                         for ( i = 0 ; i < 6 ; i ++ )
02912                         {
02913                                 int nx = 1 + neighbor6[i][0] ;
02914                                 int ny = 1 + neighbor6[i][1] ;
02915                                 int nz = 1 + neighbor6[i][2] ;
02916                                 if ( vox[nx][ny][nz] >= 0 )
02917                                 {
02918                                         tvox[ nx ][ ny ][ nz ] = 1 ;
02919 
02920                                         nv ++ ;
02921 
02922                                         for ( j = 0 ; j < 4 ; j ++ )
02923                                         {
02924                                                 int nnx = neighbor64[i][j][0] + nx ;
02925                                                 int nny = neighbor64[i][j][1] + ny ;
02926                                                 int nnz = neighbor64[i][j][2] + nz ;
02927                                                 if ( vox[nnx][nny][nnz] >= 0 )
02928                                                 {
02929                                                         if ( tvox[nnx][nny][nnz] == 0 )
02930                                                         {
02931                                                                 tvox[nnx][nny][nnz] = 1 ;
02932                                                                 nv ++ ;
02933                                                         }
02934 
02935                                                         ne ++ ;
02936                                                 }
02937                                         }
02938                                 }
02939                         }
02940 
02941                         nc = components6( tvox ) ;
02942 
02943                         return ( nc - ( nv - ne ) ) ;
02944                 }

void Volume::curveSkeleton float  thr,
Volume svol
 

Definition at line 4268 of file volume.cpp.

References EMAN::PriorityQueue< ValueT, KeyT >::add(), flag, getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isSimple(), MAX_ERODE, 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.

04269                 {
04270                         int i, j, k ;
04271                         // First, threshold the volume
04272                         #ifdef VERBOSE
04273                         printf("Thresholding the volume to -1/0...\n") ;
04274                         #endif
04275                         threshold( thr, -1, 0 ) ;
04276 
04277                         // Next, apply convergent erosion
04278                         // by preserving: complex nodes, curve end-points, and sheet points
04279 
04280                         // Next, initialize the linked queue
04281                         #ifdef VERBOSE
04282                         printf("Initializing queue...\n") ;
04283                         #endif
04284                         GridQueue2* queue2 = new GridQueue2( ) ;
04285                         GridQueue2* queue3 = new GridQueue2( ) ;
04286                         GridQueue2* queue4 = new GridQueue2( ) ;
04287                         PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
04288 
04289                         for ( i = 0 ; i < getSizeX() ; i ++ )
04290                                 for ( j = 0 ; j < getSizeY() ; j ++ )
04291                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
04292                                         {
04293                                                 if ( getDataAt( i, j, k ) >= 0 )
04294                                                 {
04295                                                         if ( svol->getDataAt(i,j,k) > 0 )
04296                                                         {
04297                                                                 setDataAt( i, j, k, MAX_ERODE ) ;
04298                                                         }
04299                                                         else
04300                                                         {
04301                                                                 for ( int m = 0 ; m < 6 ; m ++ )
04302                                                                 {
04303                                                                         if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
04304                                                                         {
04305                                                                                 // setDataAt( i, j, k, 1 ) ;
04306                                                                                 queue2->prepend( i, j, k ) ;
04307                                                                                 break ;
04308                                                                         }
04309                                                                 }
04310                                                         }
04311                                                 }
04312                                         }
04313 
04314                         int wid = MAX_ERODE ;
04315                         #ifdef VERBOSE
04316                         printf("Total %d nodes\n", queue2->getNumElements() ) ;
04317                         printf("Start erosion to %d...\n", wid) ;
04318                         #endif
04319 
04320 
04321                         // Perform erosion
04322                         gridQueueEle* ele ;
04323                         gridPoint* gp ;
04324                         int ox, oy, oz ;
04325                         int score ;
04326                         Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
04327                         for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
04328                         {
04329                                 scrvol->setDataAt( i, -1 ) ;
04330                         }
04331 
04332         #ifdef  NOISE_DIS_HELIX
04333                         Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
04334         #endif
04335 
04336                         for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
04337                         {
04338                                 // At the start of each iteration,
04339                                 // queue2 holds all the nodes for this layer
04340                                 // queue3 and queue are empty
04341 
04342                                 int numComplex = 0, numSimple = 0 ;
04343                                 #ifdef VERBOSE
04344                                 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
04345                                 #endif
04346 
04347                                 /*
04348                                 We first need to assign curwid + 1 to every node in this layer
04349                                 */
04350                                 queue2->reset() ;
04351                                 ele = queue2->getNext() ;
04352                                 while ( ele != NULL )
04353                                 {
04354                                         ox = ele->x ;
04355                                         oy = ele->y ;
04356                                         oz = ele->z ;
04357 
04358                                         if ( getDataAt(ox,oy,oz) == curwid )
04359                                         {
04360                                                 ele = queue2->remove() ;
04361                                         }
04362                                         else
04363                                         {
04364                                                 setDataAt(ox,oy,oz, curwid) ;
04365                                                 ele = queue2->getNext() ;
04366                                         }
04367                                 }
04368                                 queue4->reset() ;
04369                                 ele = queue4->getNext() ;
04370                                 while ( ele != NULL )
04371                                 {
04372                                         ox = ele->x ;
04373                                         oy = ele->y ;
04374                                         oz = ele->z ;
04375 
04376                                         queue2->prepend(ox,oy,oz) ;
04377                                         ele = queue4->remove() ;
04378                                 }
04379 
04380                                 // Now queue2 holds all the nodes for this layer
04381 
04382         #ifdef NOISE_DIS_HELIX
04383                                 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
04384                                 queue2->reset() ;
04385 
04386                                 // First run
04387                                 int flag = 0 ;
04388                                 while ( ( ele = queue2->getNext() ) != NULL )
04389                                 {
04390                                         ox = ele->x ;
04391                                         oy = ele->y ;
04392                                         oz = ele->z ;
04393                                         if ( NOISE_DIS_HELIX <= 1 )
04394                                         {
04395                                                 noisevol->setDataAt( ox, oy, oz, 0 ) ;
04396                                         }
04397                                         else
04398                                         {
04399                                                 flag = 0 ;
04400                                                 for ( int m = 0 ; m < 6 ; m ++ )
04401                                                 {
04402                                                         int nx = ox + neighbor6[m][0] ;
04403                                                         int ny = oy + neighbor6[m][1] ;
04404                                                         int nz = oz + neighbor6[m][2] ;
04405                                                         if ( getDataAt( nx, ny, nz ) == 0 )
04406                                                         {
04407                                                                 noisevol->setDataAt( ox, oy, oz, 1 ) ;
04408                                                                 flag = 1 ;
04409                                                                 break ;
04410                                                         }
04411                                                 }
04412                                                 if ( ! flag )
04413                                                 {
04414                                                         noisevol->setDataAt( ox, oy, oz, 0 ) ;
04415                                                 }
04416                                         }
04417                                 }
04418 
04419                                 int cur, visited ;
04420                                 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
04421                                 {
04422                                         queue2->reset() ;
04423                                         int count = 0 ;
04424                                         visited = 0 ;
04425 
04426                                         while ( ( ele = queue2->getNext() ) != NULL )
04427                                         {
04428                                                 ox = ele->x ;
04429                                                 oy = ele->y ;
04430                                                 oz = ele->z ;
04431 
04432                                                 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
04433                                                 {
04434                                                         visited ++ ;
04435                                                         continue ;
04436                                                 }
04437 
04438                                                 flag = 0 ;
04439                                                 for ( int m = 0 ; m < 6 ; m ++ )
04440                                                 {
04441                                                         int nx = ox + neighbor6[m][0] ;
04442                                                         int ny = oy + neighbor6[m][1] ;
04443                                                         int nz = oz + neighbor6[m][2] ;
04444                                                         if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
04445                                                         {
04446                                                                 noisevol->setDataAt( ox, oy, oz, 1 ) ;
04447                                                                 visited ++ ;
04448                                                                 count ++ ;
04449                                                                 break ;
04450                                                         }
04451                                                 }
04452                                         }
04453 
04454                                         if ( count == 0 )
04455                                         {
04456                                                 break ;
04457                                         }
04458                                 }
04459                                 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
04460 
04461 
04462         #endif
04463                                 /* Commented out for debugging
04464 
04465                                 // First,
04466                                 // check for complex nodes in queue2
04467                                 // move them from queue2 to queue3
04468                                 queue2->reset() ;
04469                                 ele = queue2->getNext() ;
04470                                 while ( ele != NULL )
04471                                 {
04472                                         ox = ele->x ;
04473                                         oy = ele->y ;
04474                                         oz = ele->z ;
04475 
04476                                         // Check simple
04477         #ifndef NOISE_DIS_HELIX
04478                                         if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04479         #else
04480                                         if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
04481         #endif
04482                                         {
04483                                                 // Complex, set to next layer
04484                                                 setDataAt( ox, oy, oz, curwid + 1 ) ;
04485                                                 queue3->prepend( ox, oy, oz ) ;
04486                                                 ele = queue2->remove() ;
04487 
04488                                                 numComplex ++ ;
04489                                         }
04490                                         else
04491                                         {
04492                                                 ele = queue2->getNext() ;
04493                                         }
04494                                 }
04495                                 */
04496 
04497                                 // Next,
04498                                 // Compute score for each node left in queue2
04499                                 // move them into priority queue
04500                                 queue2->reset() ;
04501                                 ele = queue2->getNext() ;
04502                                 while ( ele != NULL )
04503                                 {
04504                                         ox = ele->x ;
04505                                         oy = ele->y ;
04506                                         oz = ele->z ;
04507 
04508                                         // Compute score
04509                                         score = getNumPotComplex2( ox, oy, oz ) ;
04510                                         scrvol->setDataAt( ox, oy, oz, score ) ;
04511 
04512                                         // Push to queue
04513                                         gp = new gridPoint ;
04514                                         gp->x = ox ;
04515                                         gp->y = oy ;
04516                                         gp->z = oz ;
04517                                         // queue->add( gp, -score ) ;
04518                                         queue->add( gp, score ) ;
04519 
04520                                         ele = queue2->remove() ;
04521                                 }
04522 
04523                                 // Rename queue3 to be queue2,
04524                                 // Clear queue3
04525                                 // From now on, queue2 holds nodes of next level
04526                                 delete queue2 ;
04527                                 queue2 = queue3 ;
04528                                 queue3 = new GridQueue2( ) ;
04529 
04530                                 // Next, start priority queue iteration
04531                                 while ( ! queue->isEmpty() )
04532                                 {
04533                                         // Retrieve the node with the highest score
04534                                         queue->remove( gp, score ) ;
04535                                         ox = gp->x ;
04536                                         oy = gp->y ;
04537                                         oz = gp->z ;
04538                                         delete gp ;
04539         //                              score = -score ;
04540 
04541                                         // Ignore the node
04542                                         // if it has been processed before
04543                                         // or it has an updated score
04544                                         if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
04545                                         {
04546                                                 continue ;
04547                                         }
04548 
04549                                         /* Commented out for debugging
04550 
04551                                         // Remove this simple node
04552                                         setDataAt( ox, oy, oz, -1 ) ;
04553                                         numSimple ++ ;
04554                                         // printf("Highest score: %d\n", score) ;
04555                                         */
04556 
04557                                         /* Added for debugging */
04558                                         // Check simple
04559         #ifndef NOISE_DIS_HELIX
04560                                         // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
04561                                         if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04562         #else
04563                                         if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04564         #endif
04565                                         {
04566                                                 // Complex, set to next layer
04567                                                 setDataAt( ox, oy, oz, curwid + 1 ) ;
04568                                                 queue4->prepend( ox, oy, oz ) ;
04569                                                 numComplex ++ ;
04570                                         }
04571                                         else
04572                                         {
04573                                                 setDataAt( ox, oy, oz, -1 ) ;
04574                                                 numSimple ++ ;
04575                                         }
04576                                         /* Adding ends */
04577 
04578                                         // Move its neighboring unvisited node to queue2
04579                                         for ( int m = 0 ; m < 6 ; m ++ )
04580                                         {
04581                                                 int nx = ox + neighbor6[m][0] ;
04582                                                 int ny = oy + neighbor6[m][1] ;
04583                                                 int nz = oz + neighbor6[m][2] ;
04584                                                 if ( getDataAt( nx, ny, nz ) == 0 )
04585                                                 {
04586                                                         // setDataAt( nx, ny, nz, curwid + 1 ) ;
04587                                                         queue2->prepend( nx, ny, nz ) ;
04588                                                 }
04589                                         }
04590 
04591                                         /* Commented out for debugging
04592 
04593                                         // Find complex nodes in its 3x3 neighborhood
04594                                         // move them to queue2
04595                                         for ( i = -1 ; i < 2 ; i ++ )
04596                                                 for ( j = -1 ; j < 2 ; j ++ )
04597                                                         for ( k = -1 ; k < 2 ; k ++ )
04598                                                         {
04599                                                                 int nx = ox + i ;
04600                                                                 int ny = oy + j ;
04601                                                                 int nz = oz + k ;
04602 
04603                                                                 // Check simple
04604                                                                 if ( getDataAt( nx, ny, nz ) == curwid &&
04605                                                                         // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
04606         #ifndef NOISE_DIS_HELIX
04607                                                                         ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
04608         #else
04609                                                                         ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
04610         #endif
04611 
04612                                                                 {
04613                                                                         // Complex, set to next layer
04614                                                                         setDataAt( nx, ny, nz, curwid + 1 ) ;
04615                                                                         queue2->prepend( nx, ny, nz ) ;
04616                                                                         numComplex ++ ;
04617                                                                 }
04618                                                         }
04619                                         */
04620 
04621                                         // Update scores for nodes in its 5x5 neighborhood
04622                                         // insert them back into priority queue
04623 
04624                                         for ( i = -2 ; i < 3 ;i ++ )
04625                                                 for ( j = -2 ; j < 3 ; j ++ )
04626                                                         for ( k = -2 ; k < 3 ; k ++ )
04627                                                         {
04628                                                                 int nx = ox + i ;
04629                                                                 int ny = oy + j ;
04630                                                                 int nz = oz + k ;
04631 
04632                                                                 if ( getDataAt( nx, ny, nz ) == curwid )
04633                                                                 {
04634                                                                         // Compute score
04635                                                                         score = getNumPotComplex2( nx, ny, nz ) ;
04636 
04637                                                                         if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
04638                                                                         {
04639                                                                                 // printf("Update\n") ;
04640                                                                                 scrvol->setDataAt( nx, ny, nz, score ) ;
04641                                                                                 // Push to queue
04642                                                                                 gp = new gridPoint ;
04643                                                                                 gp->x = nx ;
04644                                                                                 gp->y = ny ;
04645                                                                                 gp->z = nz ;
04646                                                                                 // queue->add( gp, -score ) ;
04647                                                                                 queue->add( gp, score ) ;
04648                                                                         }
04649                                                                 }
04650                                                         }
04651 
04652 
04653                                 }
04654 
04655                                 #ifdef VERBOSE
04656                                 printf("%d complex, %d simple\n", numComplex, numSimple) ;
04657                                 #endif
04658 
04659                                 if ( numSimple == 0 )
04660                                 {
04661                                                 break ;
04662                                 }
04663                         }
04664 
04665                         // Finally, clean up
04666                         delete scrvol;
04667                         delete queue;
04668                         delete queue2;
04669                         delete queue3;
04670                         delete queue4;
04671                         #ifdef VERBOSE
04672                         printf("Thresholding the volume to 0/1...\n") ;
04673                         #endif
04674                         threshold( 0, 0, 1 ) ;
04675                 }

void Volume::curveSkeleton Volume grayvol,
float  lowthr,
float  highthr,
Volume svol
 

* 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(), flag, 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, 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().

03813                 {
03814                         int i, j, k ;
03815                         // First, threshold the volume
03816                         #ifdef VERBOSE
03817                         printf("Thresholding the volume to -1/0...\n") ;
03818                         #endif
03819                         threshold( 0.5f, -1, 0 ) ;
03820 
03821                         // Next, apply convergent erosion
03822                         // by preserving: complex nodes, curve end-points, and sheet points
03823 
03824                         // Next, initialize the linked queue
03825                         #ifdef VERBOSE
03826                         printf("Initializing queue...\n") ;
03827                         #endif
03828                         GridQueue2* queue2 = new GridQueue2( ) ;
03829                         GridQueue2* queue3 = new GridQueue2( ) ;
03830                         GridQueue2* queue4 = new GridQueue2( ) ;
03831                         PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
03832 
03833                         for ( i = 0 ; i < getSizeX() ; i ++ )
03834                                 for ( j = 0 ; j < getSizeY() ; j ++ )
03835                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
03836                                         {
03837                                                 if ( getDataAt( i, j, k ) >= 0 )
03838                                                 {
03839                                                         float v = (float)grayvol->getDataAt(i,j,k) ;
03840                                                         if ( v <= lowthr || v > highthr || svol->getDataAt(i,j,k) > 0 )
03841                                                         {
03842                                                                 setDataAt( i, j, k, MAX_ERODE ) ;
03843                                                         }
03844                                                         else
03845                                                         {
03846                                                                 for ( int m = 0 ; m < 6 ; m ++ )
03847                                                                 {
03848                                                                         if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
03849                                                                         {
03850                                                                                 // setDataAt( i, j, k, 1 ) ;
03851                                                                                 queue2->prepend( i, j, k ) ;
03852                                                                                 break ;
03853                                                                         }
03854                                                                 }
03855                                                         }
03856                                                 }
03857                                         }
03858                         int wid = MAX_ERODE ;
03859                         #ifdef VERBOSE
03860                         printf("Total %d nodes\n", queue2->getNumElements() ) ;
03861                         printf("Start erosion to %d...\n", wid) ;
03862                         #endif
03863 
03864 
03865                         // Perform erosion
03866                         gridQueueEle* ele ;
03867                         gridPoint* gp ;
03868                         int ox, oy, oz ;
03869                         int score ;
03870                         Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
03871                         for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
03872                         {
03873                                 scrvol->setDataAt( i, -1 ) ;
03874                         }
03875 
03876         #ifdef  NOISE_DIS_HELIX
03877                         Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
03878         #endif
03879 
03880                         for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
03881                         {
03882                                 // At the start of each iteration,
03883                                 // queue2 holds all the nodes for this layer
03884                                 // queue3 and queue are empty
03885 
03886                                 int numComplex = 0, numSimple = 0 ;
03887                                 #ifdef VERBOSE
03888                                 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
03889                                 #endif
03890 
03891                                 /*
03892                                 We first need to assign curwid + 1 to every node in this layer
03893                                 */
03894                                 queue2->reset() ;
03895                                 ele = queue2->getNext() ;
03896                                 while ( ele != NULL )
03897                                 {
03898                                         ox = ele->x ;
03899                                         oy = ele->y ;
03900                                         oz = ele->z ;
03901 
03902                                         if ( getDataAt(ox,oy,oz) == curwid )
03903                                         {
03904                                                 ele = queue2->remove() ;
03905                                         }
03906                                         else
03907                                         {
03908                                                 setDataAt(ox,oy,oz, curwid) ;
03909                                                 ele = queue2->getNext() ;
03910                                         }
03911                                 }
03912                                 queue4->reset() ;
03913                                 ele = queue4->getNext() ;
03914                                 while ( ele != NULL )
03915                                 {
03916                                         ox = ele->x ;
03917                                         oy = ele->y ;
03918                                         oz = ele->z ;
03919 
03920                                         queue2->prepend(ox,oy,oz) ;
03921                                         ele = queue4->remove() ;
03922                                 }
03923 
03924                                 // Now queue2 holds all the nodes for this layer
03925 
03926         #ifdef NOISE_DIS_HELIX
03927                                 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
03928                                 queue2->reset() ;
03929 
03930                                 // First run
03931                                 int flag = 0 ;
03932                                 while ( ( ele = queue2->getNext() ) != NULL )
03933                                 {
03934                                         ox = ele->x ;
03935                                         oy = ele->y ;
03936                                         oz = ele->z ;
03937                                         if ( NOISE_DIS_HELIX <= 1 )
03938                                         {
03939                                                 noisevol->setDataAt( ox, oy, oz, 0 ) ;
03940                                         }
03941                                         else
03942                                         {
03943                                                 flag = 0 ;
03944                                                 for ( int m = 0 ; m < 6 ; m ++ )
03945                                                 {
03946                                                         int nx = ox + neighbor6[m][0] ;
03947                                                         int ny = oy + neighbor6[m][1] ;
03948                                                         int nz = oz + neighbor6[m][2] ;
03949                                                         if ( getDataAt( nx, ny, nz ) == 0 )
03950                                                         {
03951                                                                 noisevol->setDataAt( ox, oy, oz, 1 ) ;
03952                                                                 flag = 1 ;
03953                                                                 break ;
03954                                                         }
03955                                                 }
03956                                                 if ( ! flag )
03957                                                 {
03958                                                         noisevol->setDataAt( ox, oy, oz, 0 ) ;
03959                                                 }
03960                                         }
03961                                 }
03962 
03963                                 int cur, visited ;
03964                                 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
03965                                 {
03966                                         queue2->reset() ;
03967                                         int count = 0 ;
03968                                         visited = 0 ;
03969 
03970                                         while ( ( ele = queue2->getNext() ) != NULL )
03971                                         {
03972                                                 ox = ele->x ;
03973                                                 oy = ele->y ;
03974                                                 oz = ele->z ;
03975 
03976                                                 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
03977                                                 {
03978                                                         visited ++ ;
03979                                                         continue ;
03980                                                 }
03981 
03982                                                 flag = 0 ;
03983                                                 for ( int m = 0 ; m < 6 ; m ++ )
03984                                                 {
03985                                                         int nx = ox + neighbor6[m][0] ;
03986                                                         int ny = oy + neighbor6[m][1] ;
03987                                                         int nz = oz + neighbor6[m][2] ;
03988                                                         if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
03989                                                         {
03990                                                                 noisevol->setDataAt( ox, oy, oz, 1 ) ;
03991                                                                 visited ++ ;
03992                                                                 count ++ ;
03993                                                                 break ;
03994                                                         }
03995                                                 }
03996                                         }
03997 
03998                                         if ( count == 0 )
03999                                         {
04000                                                 break ;
04001                                         }
04002                                 }
04003                                 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
04004 
04005 
04006         #endif
04007                                 /* Commented out for debugging
04008 
04009                                 // First,
04010                                 // check for complex nodes in queue2
04011                                 // move them from queue2 to queue3
04012                                 queue2->reset() ;
04013                                 ele = queue2->getNext() ;
04014                                 while ( ele != NULL )
04015                                 {
04016                                         ox = ele->x ;
04017                                         oy = ele->y ;
04018                                         oz = ele->z ;
04019 
04020                                         // Check simple
04021         #ifndef NOISE_DIS_HELIX
04022                                         if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04023         #else
04024                                         if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
04025         #endif
04026                                         {
04027                                                 // Complex, set to next layer
04028                                                 setDataAt( ox, oy, oz, curwid + 1 ) ;
04029                                                 queue3->prepend( ox, oy, oz ) ;
04030                                                 ele = queue2->remove() ;
04031 
04032                                                 numComplex ++ ;
04033                                         }
04034                                         else
04035                                         {
04036                                                 ele = queue2->getNext() ;
04037                                         }
04038                                 }
04039                                 */
04040 
04041                                 // Next,
04042                                 // Compute score for each node left in queue2
04043                                 // move them into priority queue
04044                                 queue2->reset() ;
04045                                 ele = queue2->getNext() ;
04046                                 while ( ele != NULL )
04047                                 {
04048                                         ox = ele->x ;
04049                                         oy = ele->y ;
04050                                         oz = ele->z ;
04051 
04052                                         // Compute score
04053                                         score = getNumPotComplex2( ox, oy, oz ) ;
04054                                         scrvol->setDataAt( ox, oy, oz, score ) ;
04055 
04056                                         // Push to queue
04057                                         gp = new gridPoint ;
04058                                         gp->x = ox ;
04059                                         gp->y = oy ;
04060                                         gp->z = oz ;
04061                                         // queue->add( gp, -score ) ;
04062                                         queue->add( gp, score ) ;
04063 
04064                                         ele = queue2->remove() ;
04065                                 }
04066 
04067                                 // Rename queue3 to be queue2,
04068                                 // Clear queue3
04069                                 // From now on, queue2 holds nodes of next level
04070                                 delete queue2 ;
04071                                 queue2 = queue3 ;
04072                                 queue3 = new GridQueue2( ) ;
04073 
04074                                 // Next, start priority queue iteration
04075                                 while ( ! queue->isEmpty() )
04076                                 {
04077                                         // Retrieve the node with the highest score
04078                                         queue->remove( gp, score ) ;
04079                                         ox = gp->x ;
04080                                         oy = gp->y ;
04081                                         oz = gp->z ;
04082                                         delete gp ;
04083         //                              score = -score ;
04084 
04085                                         // Ignore the node
04086                                         // if it has been processed before
04087                                         // or it has an updated score
04088                                         if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
04089                                         {
04090                                                 continue ;
04091                                         }
04092 
04093                                         /* Commented out for debugging
04094 
04095                                         // Remove this simple node
04096                                         setDataAt( ox, oy, oz, -1 ) ;
04097                                         numSimple ++ ;
04098                                         // printf("Highest score: %d\n", score) ;
04099                                         */
04100 
04101                                         /* Added for debugging */
04102                                         // Check simple
04103         #ifndef NOISE_DIS_HELIX
04104                                         // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
04105                                         if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04106         #else
04107                                         if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04108         #endif
04109                                         {
04110                                                 // Complex, set to next layer
04111                                                 setDataAt( ox, oy, oz, curwid + 1 ) ;
04112                                                 queue4->prepend( ox, oy, oz ) ;
04113                                                 numComplex ++ ;
04114                                         }
04115                                         else
04116                                         {
04117                                                 setDataAt( ox, oy, oz, -1 ) ;
04118                                                 numSimple ++ ;
04119 
04120 
04121                                         /* Adding ends */
04122                                                 // Move its neighboring unvisited node to queue2
04123                                                 for ( int m = 0 ; m < 6 ; m ++ )
04124                                                 {
04125                                                         int nx = ox + neighbor6[m][0] ;
04126                                                         int ny = oy + neighbor6[m][1] ;
04127                                                         int nz = oz + neighbor6[m][2] ;
04128                                                         if ( getDataAt( nx, ny, nz ) == 0 )
04129                                                         {
04130                                                                 // setDataAt( nx, ny, nz, curwid + 1 ) ;
04131                                                                 queue2->prepend( nx, ny, nz ) ;
04132                                                         }
04133                                                 }
04134 
04135                                         }
04136 
04137 
04138                                         /* Commented out for debugging
04139 
04140                                         // Find complex nodes in its 3x3 neighborhood
04141                                         // move them to queue2
04142                                         for ( i = -1 ; i < 2 ; i ++ )
04143                                                 for ( j = -1 ; j < 2 ; j ++ )
04144                                                         for ( k = -1 ; k < 2 ; k ++ )
04145                                                         {
04146                                                                 int nx = ox + i ;
04147                                                                 int ny = oy + j ;
04148                                                                 int nz = oz + k ;
04149 
04150                                                                 // Check simple
04151                                                                 if ( getDataAt( nx, ny, nz ) == curwid &&
04152                                                                         // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
04153         #ifndef NOISE_DIS_HELIX
04154                                                                         ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
04155         #else
04156                                                                         ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
04157         #endif
04158 
04159                                                                 {
04160                                                                         // Complex, set to next layer
04161                                                                         setDataAt( nx, ny, nz, curwid + 1 ) ;
04162                                                                         queue2->prepend( nx, ny, nz ) ;
04163                                                                         numComplex ++ ;
04164                                                                 }
04165                                                         }
04166                                         */
04167 
04168                                         // Update scores for nodes in its 5x5 neighborhood
04169                                         // insert them back into priority queue
04170 
04171                                         for ( i = -2 ; i < 3 ;i ++ )
04172                                                 for ( j = -2 ; j < 3 ; j ++ )
04173                                                         for ( k = -2 ; k < 3 ; k ++ )
04174                                                         {
04175                                                                 int nx = ox + i ;
04176                                                                 int ny = oy + j ;
04177                                                                 int nz = oz + k ;
04178 
04179                                                                 if ( getDataAt( nx, ny, nz ) == curwid )
04180                                                                 {
04181                                                                         // Compute score
04182                                                                         score = getNumPotComplex2( nx, ny, nz ) ;
04183 
04184                                                                         if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
04185                                                                         {
04186                                                                                 // printf("Update\n") ;
04187                                                                                 scrvol->setDataAt( nx, ny, nz, score ) ;
04188                                                                                 // Push to queue
04189                                                                                 gp = new gridPoint ;
04190                                                                                 gp->x = nx ;
04191                                                                                 gp->y = ny ;
04192                                                                                 gp->z = nz ;
04193                                                                                 // queue->add( gp, -score ) ;
04194                                                                                 queue->add( gp, score ) ;
04195                                                                         }
04196                                                                 }
04197                                                         }
04198 
04199 
04200                                 }
04201 
04202                                 #ifdef VERBOSE
04203                                 printf("%d complex, %d simple\n", numComplex, numSimple) ;
04204                                 #endif
04205 
04206                                 if ( numSimple == 0 )
04207                                 {
04208                                         if ( queue2->getNumElements() > 0 )
04209                                         {
04210                                                 printf("*************************wierd here*************************\n");
04211                                         }
04212                                                 break ;
04213                                 }
04214                         }
04215 
04216                         // Remove all internal voxels (contained in manifold surfaces)
04217                         queue2->reset() ;
04218                         queue4->reset() ;
04219                         ele = queue4->getNext() ;
04220                         while ( ele != NULL )
04221                         {
04222                                 ox = ele->x ;
04223                                 oy = ele->y ;
04224                                 oz = ele->z ;
04225 
04226                                 if ( isPiercable(ox,oy,oz) == 1 )  // hasCompleteSheet( ox, oy, oz ) == 1 ) //
04227                                 {
04228                                         queue2->prepend(ox,oy,oz) ;
04229                                 //      setDataAt( ox, oy, oz, -1 ) ;
04230                                 }
04231                                 ele = queue4->remove() ;
04232                         }
04233 
04234                         for ( i = 0 ; i < getSizeX() ; i ++ )
04235                                 for ( j = 0 ; j < getSizeY() ; j ++ )
04236                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
04237                                         {
04238                                                 if ( getDataAt( i, j, k ) == 0 && isPiercable(i,j,k) ) //hasCompleteSheet(i,j,k) == 1) //
04239                                                 {
04240                                                         queue2->prepend( i, j, k ) ;
04241                                                 }
04242                                         }
04243                         queue2->reset() ;
04244                         ele = queue2->getNext() ;
04245                         while ( ele != NULL )
04246                         {
04247                                 ox = ele->x ;
04248                                 oy = ele->y ;
04249                                 oz = ele->z ;
04250                                 setDataAt( ox, oy, oz, -1 ) ;
04251                                 ele = queue2->remove() ;
04252                         }
04253 
04254 
04255                         // Finally, clean up
04256                         delete scrvol;
04257                         delete queue;
04258                         delete queue2;
04259                         delete queue3;
04260                         delete queue4;
04261                         #ifdef VERBOSE
04262                         printf("Thresholding the volume to 0/1...\n") ;
04263                         #endif
04264                         threshold( 0, 0, 1 ) ;
04265                 }

void Volume::curveSkeleton2D float  thr,
Volume svol
 

Definition at line 4678 of file volume.cpp.

References EMAN::PriorityQueue< ValueT, KeyT >::add(), flag, getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isSimple(), MAX_ERODE, 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().

04679                 {
04680                         int i, j, k ;
04681                         // First, threshold the volume
04682                         #ifdef VERBOSE
04683                         printf("Thresholding the volume to -1/0...\n") ;
04684                         #endif
04685                         threshold( thr, -1, 0 ) ;
04686 
04687                         // Next, apply convergent erosion
04688                         // by preserving: complex nodes, curve end-points, and sheet points
04689 
04690                         // Next, initialize the linked queue
04691                         #ifdef VERBOSE
04692                         printf("Initializing queue...\n") ;
04693                         #endif
04694                         GridQueue2* queue2 = new GridQueue2( ) ;
04695                         GridQueue2* queue3 = new GridQueue2( ) ;
04696                         GridQueue2* queue4 = new GridQueue2( ) ;
04697                         PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
04698 
04699                         for ( i = 0 ; i < getSizeX() ; i ++ )
04700                                 for ( j = 0 ; j < getSizeY() ; j ++ )
04701                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
04702                                         {
04703                                                 if ( getDataAt( i, j, k ) >= 0 )
04704                                                 {
04705                                                         if ( svol->getDataAt(i,j,k) > 0 )
04706                                                         {
04707                                                                 setDataAt( i, j, k, MAX_ERODE ) ;
04708                                                         }
04709                                                         else
04710                                                         {
04711                                                                 for ( int m = 0 ; m < 4 ; m ++ )
04712                                                                 {
04713                                                                         if ( getDataAt( i + neighbor4[m][0], j + neighbor4[m][1], k ) < 0 )
04714                                                                         {
04715                                                                                 // setDataAt( i, j, k, 1 ) ;
04716                                                                                 queue2->prepend( i, j, k ) ;
04717                                                                                 break ;
04718                                                                         }
04719                                                                 }
04720                                                         }
04721                                                 }
04722                                         }
04723                         int wid = MAX_ERODE ;
04724                         #ifdef VERBOSE
04725                         printf("Total %d nodes\n", queue2->getNumElements() ) ;
04726                         printf("Start erosion to %d...\n", wid) ;
04727                         #endif
04728 
04729 
04730                         // Perform erosion
04731                         gridQueueEle* ele ;
04732                         gridPoint* gp ;
04733                         int ox, oy, oz ;
04734                         int score ;
04735                         Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
04736                         for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
04737                         {
04738                                 scrvol->setDataAt( i, -1 ) ;
04739                         }
04740 
04741         #ifdef  NOISE_DIS_HELIX
04742                         Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
04743         #endif
04744 
04745                         for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
04746                         {
04747                                 // At the start of each iteration,
04748                                 // queue2 holds all the nodes for this layer
04749                                 // queue3 and queue are empty
04750 
04751                                 int numComplex = 0, numSimple = 0 ;
04752                                 #ifdef VERBOSE
04753                                 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
04754                                 #endif
04755 
04756                                 /*
04757                                 We first need to assign curwid + 1 to every node in this layer
04758                                 */
04759                                 queue2->reset() ;
04760                                 ele = queue2->getNext() ;
04761                                 while ( ele != NULL )
04762                                 {
04763                                         ox = ele->x ;
04764                                         oy = ele->y ;
04765                                         oz = ele->z ;
04766 
04767                                         if ( getDataAt(ox,oy,oz) == curwid )
04768                                         {
04769                                                 ele = queue2->remove() ;
04770                                         }
04771                                         else
04772                                         {
04773                                                 setDataAt(ox,oy,oz, curwid) ;
04774                                                 ele = queue2->getNext() ;
04775                                         }
04776                                 }
04777                                 queue4->reset() ;
04778                                 ele = queue4->getNext() ;
04779                                 while ( ele != NULL )
04780                                 {
04781                                         ox = ele->x ;
04782                                         oy = ele->y ;
04783                                         oz = ele->z ;
04784 
04785                                         queue2->prepend(ox,oy,oz) ;
04786                                         ele = queue4->remove() ;
04787                                 }
04788 
04789                                 // Now queue2 holds all the nodes for this layer
04790 
04791                                 #ifdef NOISE_DIS_HELIX
04792                                 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
04793                                 queue2->reset() ;
04794 
04795                                 // First run
04796                                 int flag = 0 ;
04797                                 while ( ( ele = queue2->getNext() ) != NULL )
04798                                 {
04799                                         ox = ele->x ;
04800                                         oy = ele->y ;
04801                                         oz = ele->z ;
04802                                         if ( NOISE_DIS_HELIX <= 1 )
04803                                         {
04804                                                 noisevol->setDataAt( ox, oy, oz, 0 ) ;
04805                                         }
04806                                         else
04807                                         {
04808                                                 flag = 0 ;
04809                                                 for ( int m = 0 ; m < 6 ; m ++ )
04810                                                 {
04811                                                         int nx = ox + neighbor6[m][0] ;
04812                                                         int ny = oy + neighbor6[m][1] ;
04813                                                         int nz = oz + neighbor6[m][2] ;
04814                                                         if ( getDataAt( nx, ny, nz ) == 0 )
04815                                                         {
04816                                                                 noisevol->setDataAt( ox, oy, oz, 1 ) ;
04817                                                                 flag = 1 ;
04818                                                                 break ;
04819                                                         }
04820                                                 }
04821                                                 if ( ! flag )
04822                                                 {
04823                                                         noisevol->setDataAt( ox, oy, oz, 0 ) ;
04824                                                 }
04825                                         }
04826                                 }
04827 
04828                                 int cur, visited ;
04829                                 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
04830                                 {
04831                                         queue2->reset() ;
04832                                         int count = 0 ;
04833                                         visited = 0 ;
04834 
04835                                         while ( ( ele = queue2->getNext() ) != NULL )
04836                                         {
04837                                                 ox = ele->x ;
04838                                                 oy = ele->y ;
04839                                                 oz = ele->z ;
04840 
04841                                                 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
04842                                                 {
04843                                                         visited ++ ;
04844                                                         continue ;
04845                                                 }
04846 
04847                                                 flag = 0 ;
04848                                                 for ( int m = 0 ; m < 6 ; m ++ )
04849                                                 {
04850                                                         int nx = ox + neighbor6[m][0] ;
04851                                                         int ny = oy + neighbor6[m][1] ;
04852                                                         int nz = oz + neighbor6[m][2] ;
04853                                                         if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
04854                                                         {
04855                                                                 noisevol->setDataAt( ox, oy, oz, 1 ) ;
04856                                                                 visited ++ ;
04857                                                                 count ++ ;
04858                                                                 break ;
04859                                                         }
04860                                                 }
04861                                         }
04862 
04863                                         if ( count == 0 )
04864                                         {
04865                                                 break ;
04866                                         }
04867                                 }
04868                                 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
04869 
04870 
04871                                 #endif
04872                                 /* Commented out for debugging
04873 
04874                                 // First,
04875                                 // check for complex nodes in queue2
04876                                 // move them from queue2 to queue3
04877                                 queue2->reset() ;
04878                                 ele = queue2->getNext() ;
04879                                 while ( ele != NULL )
04880                                 {
04881                                         ox = ele->x ;
04882                                         oy = ele->y ;
04883                                         oz = ele->z ;
04884 
04885                                         // Check simple
04886         #ifndef NOISE_DIS_HELIX
04887                                         if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04888         #else
04889                                         if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
04890         #endif
04891                                         {
04892                                                 // Complex, set to next layer
04893                                                 setDataAt( ox, oy, oz, curwid + 1 ) ;
04894                                                 queue3->prepend( ox, oy, oz ) ;
04895                                                 ele = queue2->remove() ;
04896 
04897                                                 numComplex ++ ;
04898                                         }
04899                                         else
04900                                         {
04901                                                 ele = queue2->getNext() ;
04902                                         }
04903                                 }
04904                                 */
04905 
04906                                 // Next,
04907                                 // Compute score for each node left in queue2
04908                                 // move them into priority queue
04909                                 queue2->reset() ;
04910                                 ele = queue2->getNext() ;
04911                                 while ( ele != NULL )
04912                                 {
04913                                         ox = ele->x ;
04914                                         oy = ele->y ;
04915                                         oz = ele->z ;
04916 
04917                                         // Compute score
04918                                         score = getNumPotComplex2( ox, oy, oz ) ;
04919                                         //score = getNumNeighbor6( ox, oy, oz ) ;
04920                                         scrvol->setDataAt( ox, oy, oz, score ) ;
04921 
04922                                         // Push to queue
04923                                         gp = new gridPoint ;
04924                                         gp->x = ox ;
04925                                         gp->y = oy ;
04926                                         gp->z = oz ;
04927                                         // queue->add( gp, -score ) ;
04928                                         queue->add( gp, score ) ;
04929 
04930                                         ele = queue2->remove() ;
04931                                 }
04932 
04933                                 // Rename queue3 to be queue2,
04934                                 // Clear queue3
04935                                 // From now on, queue2 holds nodes of next level
04936                                 delete queue2 ;
04937                                 queue2 = queue3 ;
04938                                 queue3 = new GridQueue2( ) ;
04939 
04940                                 // Next, start priority queue iteration
04941                                 while ( ! queue->isEmpty() )
04942                                 {
04943                                         // Retrieve the node with the highest score
04944                                         queue->remove( gp, score ) ;
04945                                         ox = gp->x ;
04946                                         oy = gp->y ;
04947                                         oz = gp->z ;
04948                                         delete gp ;
04949         //                              score = -score ;
04950 
04951                                         // Ignore the node
04952                                         // if it has been processed before
04953                                         // or it has an updated score
04954                                         if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
04955                                         {
04956                                                 continue ;
04957                                         }
04958 
04959                                         /* Commented out for debugging
04960 
04961                                         // Remove this simple node
04962                                         setDataAt( ox, oy, oz, -1 ) ;
04963                                         numSimple ++ ;
04964                                         // printf("Highest score: %d\n", score) ;
04965                                         */
04966 
04967                                         /* Added for debugging */
04968                                         // Check simple
04969                                         #ifndef NOISE_DIS_HELIX
04970                                         // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
04971                                         if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04972                                         #else
04973                                         if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04974                                         #endif
04975                                         {
04976                                                 // Complex, set to next layer
04977                                                 setDataAt( ox, oy, oz, curwid + 1 ) ;
04978                                                 queue4->prepend( ox, oy, oz ) ;
04979                                                 numComplex ++ ;
04980                                         }
04981                                         else
04982                                         {
04983                                                 setDataAt( ox, oy, oz, -1 ) ;
04984                                                 numSimple ++ ;
04985                                         }
04986                                         /* Adding ends */
04987 
04988                                         // Move its neighboring unvisited node to queue2
04989                                         for ( int m = 0 ; m < 4 ; m ++ )
04990                                         {
04991                                                 int nx = ox + neighbor4[m][0] ;
04992                                                 int ny = oy + neighbor4[m][1] ;
04993                                                 int nz = oz ;
04994                                                 if ( getDataAt( nx, ny, nz ) == 0 )
04995                                                 {
04996                                                         // setDataAt( nx, ny, nz, curwid + 1 ) ;
04997                                                         queue2->prepend( nx, ny, nz ) ;
04998                                                 }
04999                                         }
05000 
05001                                         /* Commented out for debugging
05002 
05003                                         // Find complex nodes in its 3x3 neighborhood
05004                                         // move them to queue2
05005                                         for ( i = -1 ; i < 2 ; i ++ )
05006                                                 for ( j = -1 ; j < 2 ; j ++ )
05007                                                         for ( k = -1 ; k < 2 ; k ++ )
05008                                                         {
05009                                                                 int nx = ox + i ;
05010                                                                 int ny = oy + j ;
05011                                                                 int nz = oz + k ;
05012 
05013                                                                 // Check simple
05014                                                                 if ( getDataAt( nx, ny, nz ) == curwid &&
05015                                                                         // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
05016         #ifndef NOISE_DIS_HELIX
05017                                                                         ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
05018         #else
05019                                                                         ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
05020         #endif
05021 
05022                                                                 {
05023                                                                         // Complex, set to next layer
05024                                                                         setDataAt( nx, ny, nz, curwid + 1 ) ;
05025                                                                         queue2->prepend( nx, ny, nz ) ;
05026                                                                         numComplex ++ ;
05027                                                                 }
05028                                                         }
05029                                         */
05030 
05031                                         // Update scores for nodes in its 5x5 neighborhood
05032                                         // insert them back into priority queue
05033 
05034                                         for ( i = -2 ; i < 3 ;i ++ )
05035                                                 for ( j = -2 ; j < 3 ; j ++ )
05036                                                         for ( k = -2 ; k < 3 ; k ++ )
05037                                                         {
05038                                                                 int nx = ox + i ;
05039                                                                 int ny = oy + j ;
05040                                                                 int nz = oz + k ;
05041 
05042                                                                 if ( getDataAt( nx, ny, nz ) == curwid )
05043                                                                 {
05044                                                                         // Compute score
05045                                                                         score = getNumPotComplex2( nx, ny, nz ) ;
05046                                                                         //score = getNumNeighbor6( nx, ny, nz ) ;
05047 
05048                                                                         if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
05049                                                                         {
05050                                                                                 // printf("Update\n") ;
05051                                                                                 scrvol->setDataAt( nx, ny, nz, score ) ;
05052                                                                                 // Push to queue
05053                                                                                 gp = new gridPoint ;
05054                                                                                 gp->x = nx ;
05055                                                                                 gp->y = ny ;
05056                                                                                 gp->z = nz ;
05057                                                                                 // queue->add( gp, -score ) ;
05058                                                                                 queue->add( gp, score ) ;
05059                                                                         }
05060                                                                 }
05061                                                         }
05062 
05063 
05064                                 }
05065 
05066                                 #ifdef VERBOSE
05067                                 printf("%d complex, %d simple\n", numComplex, numSimple) ;
05068                                 #endif
05069 
05070                                 if ( numSimple == 0 )
05071                                 {
05072                                                 break ;
05073                                 }
05074                         }
05075 
05076                         // Finally, clean up
05077                         #ifdef VERBOSE
05078                         printf("Thresholding the volume to 0/1...\n") ;
05079                         #endif
05080                         threshold( 0, 0, 1 ) ;
05081                         delete scrvol;
05082                         delete queue;
05083                         delete queue2;
05084                         delete queue3;
05085                         delete queue4;
05086                 }

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.

05912                 {
05913                         int i, j, k ;
05914                         // First, threshold the volume
05915                         //printf("Thresholding the volume to -1/0...\n") ;
05916                         threshold( 0.1f, -1, 0 ) ;
05917 
05918                         /* Debug: remove faces */
05919                         //Volume* facevol = markFaceEdge() ;
05920                         /* End debugging */
05921 
05922                         // Next, initialize the linked queue
05923                         //printf("Initializing queue...\n") ;
05924                         Volume * fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
05925                         GridQueue2* queue2 = new GridQueue2( ) ;
05926                         GridQueue2* queue3 = new GridQueue2( ) ;
05927                         GridQueue2** queues = new GridQueue2* [ 10000 ] ;
05928 
05929                         for ( i = 0 ; i < getSizeX() ; i ++ )
05930                                 for ( j = 0 ; j < getSizeY() ; j ++ )
05931                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
05932                                         {
05933                                                 if ( getDataAt( i, j, k ) >= 0 )
05934                                                 {
05935                                                         if ( ! hasCompleteHelix( i, j, k ) )
05936                                                         // if ( ! hasCompleteHelix( i, j, k, facevol ) )
05937                                                         {
05938                                                                 queue2->prepend( i, j, k ) ;
05939                                                                 fvol->setDataAt( i, j, k, -1 ) ;
05940                                                         }
05941                                                 }
05942                                         }
05943                         //printf("Total %d nodes\n", queue2->getNumElements() ) ;
05944 
05945                         // Start erosion
05946                         gridQueueEle* ele ;
05947                         int dis = -1 ;
05948                         while ( queue2->getNumElements() > 0 )
05949                         {
05950                                 // First, set distance
05951                                 dis -- ;
05952                                 queues[ - dis ] = new GridQueue2( ) ;
05953                                 //printf("Distance transform to %d...", dis) ;
05954                                 queue2->reset() ;
05955                                 while( ( ele = queue2->getNext() ) != NULL )
05956                                 {
05957                                         setDataAt( ele->x, ele->y, ele->z, dis ) ;
05958                                         queues[ -dis ]->prepend( ele->x, ele->y, ele->z ) ;
05959                                 }
05960                                 //printf("%d nodes\n", queues[-dis]->getNumElements()) ;
05961 
05962                                 // Next, find next layer
05963                                 queue2->reset() ;
05964                                 ele = queue2->getNext() ;
05965                                 while ( ele != NULL )
05966                                 {
05967                                         for ( int m = 0 ; m < 6 ; m ++ )
05968                                         {
05969                                                 int nx = ele->x + neighbor6[m][0] ;
05970                                                 int ny = ele->y + neighbor6[m][1] ;
05971                                                 int nz = ele->z + neighbor6[m][2] ;
05972                                                 if ( getDataAt( nx, ny, nz ) == 0 )
05973                                                 {
05974                                                         fvol->setDataAt( nx, ny, nz, dis ) ;
05975 
05976                                                         if ( ! hasCompleteHelix( nx, ny, nz ) )
05977                                                         // if ( ! hasCompleteHelix( nx, ny, nz, facevol ) )
05978                                                         {
05979                                                                 setDataAt( nx, ny, nz, 1 ) ;
05980                                                                 queue3->prepend( nx, ny, nz ) ;
05981                                                         }
05982                                                 }
05983                                         }
05984 
05985                                         ele = queue2->remove() ;
05986                                 }
05987 
05988                                 // Next, swap queues
05989                                 GridQueue2* temp = queue2 ;
05990                                 queue2 = queue3 ;
05991                                 queue3 = temp ;
05992                         }
05993 
05994                         // Deal with closed rings
05995                         dis -- ;
05996                         queues[ - dis ] = new GridQueue2( ) ;
05997                         #ifdef VERBOSE
05998                         printf("Detecting closed rings %d...", dis) ;
05999                         #endif
06000                         int ftot = 0 ;
06001                         for ( i = 0 ; i < getSizeX() ; i ++ )
06002                                 for ( j = 0 ; j < getSizeY() ; j ++ )
06003                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
06004                                         {
06005                                                 if ( getDataAt( i, j, k ) == 0 )
06006                                                 {
06007                                                         setDataAt( i, j, k, dis ) ;
06008                                                         queues[ -dis ]->prepend( i, j, k ) ;
06009         /*
06010                                                         int fval = (int) fvol->getDataAt( i, j, k ) ;
06011                                                         if ( fval == 0)
06012                                                         {
06013                                                                 // queues[ -dis ]->prepend( i, j, k ) ;
06014                                                         }
06015                                                         else
06016                                                         {
06017                                                                 setDataAt( i, j, k, fval - 1 ) ;
06018                                                                 // queues[ -fval + 1 ]->prepend( i, j, k ) ;
06019                                                         }
06020         */
06021                                                         ftot ++ ;
06022                                                 }
06023                                         }
06024                         #ifdef VERBOSE
06025                         printf("%d nodes\n", ftot) ;
06026                         #endif
06027 
06028 
06029                         // return ;
06030 
06031                         /* Find local minimum: to help determine erosion level
06032                         int cts[ 64 ] ;
06033                         for ( i = 0 ; i <= - dis ; i ++ )
06034                         {
06035                                 cts[ i ] = 0 ;
06036                         }
06037                         for ( i = 0 ; i < getSizeX() ; i ++ )
06038                                 for ( j = 0 ; j < getSizeY() ; j ++ )
06039                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
06040                                         {
06041                                                 double val = getDataAt( i, j, k ) ;
06042                                                 if ( val < -1 )
06043                                                 {
06044                                                         int min = 1 ;
06045                                                         for ( int m = 0 ; m < 6 ; m ++ )
06046                                                         {
06047                                                                 int nx = i + neighbor6[m][0] ;
06048                                                                 int ny = j + neighbor6[m][1] ;
06049                                                                 int nz = k + neighbor6[m][2] ;
06050                                                                 if ( getDataAt( nx, ny, nz ) < val )
06051                                                                 {
06052                                                                         min = 0 ;
06053                                                                         break ;
06054                                                                 }
06055                                                         }
06056 
06057                                                         if ( min )
06058                                                         {
06059                                                                 cts[ (int) - val ] ++ ;
06060                                                         }
06061                                                 }
06062                                         }
06063 
06064                         for ( i = 2 ; i <= - dis ; i ++ )
06065                         {
06066                                 printf("Local minima: %d with %d nodes.\n", -i, cts[ i ] ) ;
06067                         }
06068                         */
06069 
06070                         // Dilate back
06071                         // Starting from nodes with distance - 2 - disthr
06072 
06073                         if ( disthr + 2 > - dis )
06074                         {
06075                                 disthr = - dis - 2 ;
06076                         }
06077                         int d;
06078 
06079                         for ( d = - dis ; d > disthr + 1 ; d -- )
06080                         {
06081                                 queues[ d ]->reset() ;
06082                                 while ( (ele = queues[ d ]->getNext() ) != NULL )
06083                                 {
06084                                         setDataAt( ele->x, ele->y, ele->z, d ) ;
06085                                 }
06086                         }
06087 
06088 
06089                         for ( d = disthr + 1 ; d >= 2 ; d -- )
06090                         {
06091                                 //delete queue3 ;
06092                                 //queue3 = new GridQueue2( ) ;
06093                                 queues[ d ]->reset() ;
06094                                 ele = queues[ d ]->getNext() ;
06095                                 while ( ele != NULL )
06096                                 {
06097                                         int dilatable = 0 ;
06098                                         for ( int m = 0 ; m < 6 ; m ++ )
06099                                                         {
06100                                                                 int nx = ele->x + neighbor6[m][0] ;
06101                                                                 int ny = ele->y + neighbor6[m][1] ;
06102                                                                 int nz = ele->z + neighbor6[m][2] ;
06103                                                                 if ( getDataAt( nx, ny, nz ) == d + 1 )
06104                                                                 {
06105                                                                         dilatable = 1 ;
06106                                                                         break ;
06107                                                                 }
06108                                                         }
06109 
06110 
06111                                         if ( ! dilatable )
06112                                         {
06113                                                 /*
06114                                                 setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
06115                                                 */
06116 
06117                                                 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
06118                                                 if ( d > 2 )
06119                                                 {
06120                                                         queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
06121                                                 }
06122                                                 ele = queues[ d ]->remove() ;
06123                                         }
06124                                         else
06125                                         {
06126                                                 setDataAt( ele->x, ele->y, ele->z, d ) ;
06127                                                 ele = queues[ d ]->getNext() ;
06128                                         }
06129 
06130                                 }
06131 
06132                                 /* Debug: extract minimal set */
06133                                 while ( 1 )
06134                                 {
06135                                         int num = 0 ;
06136                                         queues[ d ]->reset() ;
06137                                         ele = queues[ d ]->getNext() ;
06138                                         while ( ele != NULL )
06139                                         {
06140                                                 int critical = 0 ;
06141                                                 setDataAt( ele->x, ele->y, ele->z, -1 ) ;
06142 
06143                                                 for ( i = -1 ; i < 2 ; i ++ )
06144                                                 {
06145                                                         for ( j = -1 ; j < 2 ; j ++ )
06146                                                         {
06147                                                                 for ( k = -1 ; k < 2 ; k ++ )
06148                                                                 {
06149                                                                         if ( i != 0 && j != 0 && k != 0 )
06150                                                                         {
06151                                                                                 continue ;
06152                                                                         }
06153                                                                         int nx = ele->x + i ;
06154                                                                         int ny = ele->y + j ;
06155                                                                         int nz = ele->z + k ;
06156                                                                         if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteHelix( nx,ny,nz ) ) //, facevol ) )
06157                                                                         {
06158                                                                                 critical = 1 ;
06159                                                                                 break ;
06160                                                                         }
06161                                                                 }
06162                                                                 if ( critical )
06163                                                                 {
06164                                                                         break ;
06165                                                                 }
06166                                                         }
06167                                                         if ( critical )
06168                                                         {
06169                                                                 break ;
06170                                                         }
06171                                                 }
06172 
06173                                                 if ( critical )
06174                                                 {
06175                                                         setDataAt( ele->x, ele->y, ele->z, d ) ;
06176                                                         ele = queues[ d ]->getNext() ;
06177                                                 }
06178                                                 else
06179                                                 {
06180                                                         setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
06181                                                         if ( d > 2 )
06182                                                         {
06183                                                                 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
06184                                                         }
06185                                                         ele = queues[ d ]->remove() ;
06186                                                         num ++ ;
06187                                                 }
06188 
06189                                         }
06190 
06191                                         #ifdef VERBOSE
06192                                         printf("Non-minimal: %d\n", num) ;
06193                                         #endif
06194 
06195                                         if ( num == 0 )
06196                                         {
06197                                                 break ;
06198                                         }
06199                                 }
06200 
06201 
06202                                 /* End of debugging */
06203 
06204                                 /*
06205                                 queue3->reset() ;
06206                                 ele = queue3->getNext() ;
06207                                 while ( ele != NULL )
06208                                 {
06209                                         setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
06210                                         ele = queue3->remove() ;
06211                                 }
06212                                 */
06213                         }
06214 
06215                         // Finally, threshold the volume
06216                         #ifdef VERBOSE
06217                         //printf("Thresholding the volume to 0/1...\n") ;
06218                         #endif
06219                         //threshold( -1, 1, 0, 0 ) ;
06220                         threshold( 0, 0, 1 ) ;
06221                         delete fvol ;
06222                         delete queue2;
06223                         delete queue3;
06224                         for ( d = -dis ; d >= 2 ; d -- ) {
06225                                 delete queues[d];
06226                         }
06227                         delete [] queues;
06228 
06229                 }

void Volume::erodeHelix  ) 
 

Definition at line 5905 of file volume.cpp.

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

05906                 {
05907                         erodeHelix( 3 ) ;
05908                 }

int Volume::erodeSheet int  disthr  ) 
 

Definition at line 6240 of file volume.cpp.

References flag, 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(), Volume(), wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridQueueEle::y, and wustl_mm::SkeletonMaker::gridQueueEle::z.

06241                 {
06242                         int i, j, k ;
06243                         // First, threshold the volume
06244                         //printf("Thresholding the volume to -1/0...\n") ;
06245                         threshold( 0.1f, -1, 0 ) ;
06246 
06247                         /* Debug: remove cells */
06248                         Volume* facevol = markCellFace() ;
06249                         /* End debugging */
06250 
06251                         // Next, initialize the linked queue
06252                         //printf("Initializing queue...\n") ;
06253                         Volume * fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
06254                         GridQueue2* queue2 = new GridQueue2( ) ;
06255                         GridQueue2* queue3 = new GridQueue2( ) ;
06256                         GridQueue2** queues = new GridQueue2* [ 10000 ] ;
06257 
06258                         for ( i = 0 ; i < getSizeX() ; i ++ )
06259                                 for ( j = 0 ; j < getSizeY() ; j ++ )
06260                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
06261                                         {
06262                                                 if ( getDataAt( i, j, k ) >= 0 )
06263                                                 {
06264                                                         if ( ! hasCompleteSheet( i, j, k ) )
06265                                                         //if ( ! hasCompleteSheet( i, j, k, facevol ) )
06266                                                         {
06267                                                                 queue2->prepend( i, j, k ) ;
06268                                                                 fvol->setDataAt( i, j, k, -1 ) ;
06269                                                         }
06270                                                 }
06271                                         }
06272                         #ifdef VERBOSE
06273                         printf("Total %d nodes\n", queue2->getNumElements() ) ;
06274                         #endif
06275 
06276                         // Start erosion
06277                         gridQueueEle* ele ;
06278                         int dis = -1 ;
06279                         while ( queue2->getNumElements() > 0 )
06280                         {
06281                                 // First, set distance
06282                                 dis -- ;
06283                                 queues[ - dis ] = new GridQueue2( ) ;
06284                                 //printf("Distance transform to %d...", dis) ;
06285                                 queue2->reset() ;
06286                                 while( ( ele = queue2->getNext() ) != NULL )
06287                                 {
06288                                         setDataAt( ele->x, ele->y, ele->z, dis ) ;
06289                                         queues[ -dis ]->prepend( ele->x, ele->y, ele->z ) ;
06290                                 }
06291                                 //printf("%d nodes\n", queues[-dis]->getNumElements()) ;
06292 
06293                                 // Next, find next layer
06294                                 queue2->reset() ;
06295                                 ele = queue2->getNext() ;
06296                                 while ( ele != NULL )
06297                                 {
06298                                         // for ( int m = 0 ; m < 6 ; m ++ )
06299                                         for ( int mx = -1 ; mx < 2 ; mx ++ )
06300                                                 for ( int my = -1 ; my < 2 ; my ++ )
06301                                                         for ( int mz = -1 ; mz < 2 ; mz ++ )
06302                                                         {
06303                                                                 if ( mx != 0 && my != 0 && mz != 0 )
06304                                                                 {
06305                                                                         continue ;
06306                                                                 }
06307                                                                 //int nx = ele->x + neighbor6[m][0] ;
06308                                                                 //int ny = ele->y + neighbor6[m][1] ;
06309                                                                 //int nz = ele->z + neighbor6[m][2] ;
06310                                                                 int nx = ele->x + mx ;
06311                                                                 int ny = ele->y + my ;
06312                                                                 int nz = ele->z + mz ;
06313 
06314                                                                 if ( getDataAt( nx, ny, nz ) == 0 )
06315                                                                 {
06316                                                                         fvol->setDataAt( nx, ny, nz, dis ) ;
06317 
06318                                                                         if  ( ! hasCompleteSheet( nx, ny, nz ) )
06319                                                                         // if  ( ! hasCompleteSheet( nx, ny, nz, facevol ) )
06320                                                                         {
06321                                                                                 setDataAt( nx, ny, nz, 1 ) ;
06322                                                                                 queue3->prepend( nx, ny, nz ) ;
06323                                                                         }
06324                                                                 }
06325                                                         }
06326 
06327                                         ele = queue2->remove() ;
06328                                 }
06329 
06330                                 // Next, swap queues
06331                                 GridQueue2* temp = queue2 ;
06332                                 queue2 = queue3 ;
06333                                 queue3 = temp ;
06334                         }
06335 
06336                         /* Deal with closed rings */
06337 
06338                         dis -- ;
06339                         queues[ - dis ] = new GridQueue2( ) ;
06340                         #ifdef VERBOSE
06341                         printf("Detecting closed rings %d...", dis) ;
06342                         #endif
06343                         int ftot = 0 ;
06344                         for ( i = 0 ; i < getSizeX() ; i ++ )
06345                                 for ( j = 0 ; j < getSizeY() ; j ++ )
06346                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
06347                                         {
06348                                                 if ( getDataAt( i, j, k ) == 0 )
06349                                                 {
06350                                                         /*
06351                                                         int fval = (int) fvol->getDataAt( i, j, k ) ;
06352                                                         if ( fval == 0)
06353                                                         {
06354                                                                 setDataAt( i, j, k, dis - 2 ) ;
06355                                                                 // queues[ -dis ]->prepend( i, j, k ) ;
06356                                                         }
06357                                                         else
06358                                                         {
06359                                                                 setDataAt( i, j, k, fval - 1 ) ;
06360                                                                 queues[ -fval + 1 ]->prepend( i, j, k ) ;
06361                                                         }
06362                                                         */
06363                                                         setDataAt( i, j, k, dis ) ;
06364                                                         queues[ -dis ]->prepend( i, j, k ) ;
06365 
06366                                                         ftot ++ ;
06367                                                 }
06368                                         }
06369                         #ifdef VERBOSE
06370                         printf("%d nodes\n", ftot) ;
06371                         #endif
06372 
06373 
06374                         /* Find local minimum: to help determine erosion level
06375                         int cts[ 64 ] ;
06376                         for ( i = 0 ; i <= - dis ; i ++ )
06377                         {
06378                                 cts[ i ] = 0 ;
06379                         }
06380                         for ( i = 0 ; i < getSizeX() ; i ++ )
06381                                 for ( j = 0 ; j < getSizeY() ; j ++ )
06382                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
06383                                         {
06384                                                 double val = getDataAt( i, j, k ) ;
06385                                                 if ( val < -1 )
06386                                                 {
06387                                                         int min = 1 ;
06388                                                         for ( int m = 0 ; m < 6 ; m ++ )
06389                                                         {
06390                                                                 int nx = i + neighbor6[m][0] ;
06391                                                                 int ny = j + neighbor6[m][1] ;
06392                                                                 int nz = k + neighbor6[m][2] ;
06393                                                                 if ( getDataAt( nx, ny, nz ) < val )
06394                                                                 {
06395                                                                         min = 0 ;
06396                                                                         break ;
06397                                                                 }
06398                                                         }
06399 
06400                                                         if ( min )
06401                                                         {
06402                                                                 cts[ (int) - val ] ++ ;
06403                                                         }
06404                                                 }
06405                                         }
06406 
06407                         for ( i = 2 ; i <= - dis ; i ++ )
06408                         {
06409                                 printf("Local minima: %d with %d nodes.\n", -i, cts[ i ] ) ;
06410                         }
06411                         */
06412 
06413                         // return ;
06414 
06415                         // Dilate back
06416                         // Starting from nodes with distance - 2 - disthr
06417                         int d ;
06418                         if ( disthr + 2 > - dis )
06419                         {
06420                                 disthr = - dis - 2 ;
06421 
06422                         }
06423                         for ( d = - dis ; d > disthr + 1 ; d -- )
06424                         {
06425                                 queues[ d ]->reset() ;
06426                                 while ( (ele = queues[ d ]->getNext() ) != NULL )
06427                                 {
06428                                         setDataAt( ele->x, ele->y, ele->z, d ) ;
06429                                 }
06430                         }
06431 
06432                         for (d = disthr + 1 ; d >= 2 ; d -- )
06433                         {
06434 
06435                                 //delete queue3 ;
06436                                 //queue3 = new GridQueue2( ) ;
06437                                 queues[ d ]->reset() ;
06438                                 ele = queues[ d ]->getNext() ;
06439                                 while ( ele != NULL )
06440                                 {
06441                                         int dilatable = 0 ;
06442                                         // for ( int m = 0 ; m < 6 ; m ++ )
06443                                         /*
06444                                         for ( int mx = -1 ; mx < 2 ; mx ++ )
06445                                         for ( int my = -1 ; my < 2 ; my ++ )
06446                                         for ( int mz = -1 ; mz < 2 ; mz ++ )
06447                                         {
06448                                         if ( mx == 0 || my == 0 || mz == 0 )
06449                                         {
06450                                         int nx = ele->x + mx ; // neighbor6[m][0] ;
06451                                         int ny = ele->y + my ; // neighbor6[m][1] ;
06452                                         int nz = ele->z + mz ; // neighbor6[m][2] ;
06453                                         if ( getDataAt( nx, ny, nz ) == - d - 1 )
06454                                         {
06455                                         dilatable = 1 ;
06456                                         break ;
06457                                         }
06458                                         }
06459                                         }
06460                                         */
06461                                         for ( i = 0 ; i < 12 ; i ++ )
06462                                         {
06463                                                 int flag = 1, flag2 = 0 ;
06464                                                 for ( j = 0 ; j < 4 ; j ++ )
06465                                                 {
06466                                                         int nx = ele->x + sheetNeighbor[i][j][0] ;
06467                                                         int ny = ele->y + sheetNeighbor[i][j][1] ;
06468                                                         int nz = ele->z + sheetNeighbor[i][j][2] ;
06469 
06470                                                         double val = getDataAt( nx, ny, nz ) ;
06471 
06472                                                         if ( val > - d && val < 0)
06473                                                         {
06474                                                                 flag = 0 ;
06475                                                                 break ;
06476                                                         }
06477                                                         if ( val == d + 1 )
06478                                                         {
06479                                                                 flag2 ++ ;
06480                                                         }
06481                                                 }
06482 
06483                                                 if ( flag && flag2 )
06484                                                 {
06485                                                         dilatable = 1 ;
06486                                                         break ;
06487                                                 }
06488                                         }
06489 
06490                                         if ( ! dilatable )
06491                                         {
06492                                                 // setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
06493                                                 // queue3->prepend( ele->x, ele->y, ele->z ) ;
06494 
06495                                                 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
06496                                                 if ( d > 2 )
06497                                                 {
06498                                                         queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
06499                                                 }
06500                                                 ele = queues[ d ]->remove() ;
06501                                         }
06502                                         else
06503                                         {
06504                                                 setDataAt( ele->x, ele->y, ele->z, d ) ;
06505                                                 ele = queues[ d ]->getNext() ;
06506                                         }
06507                                 }
06508 
06509                                 /* Debug: extract minimal set */
06510                                 while ( 1 )
06511                                 {
06512                                         int num = 0 ;
06513                                         queues[ d ]->reset() ;
06514                                         ele = queues[ d ]->getNext() ;
06515                                         while ( ele != NULL )
06516                                         {
06517                                                 int critical = 0 ;
06518                                                 setDataAt( ele->x, ele->y, ele->z, -1 ) ;
06519 
06520                                                 for ( i = -1 ; i < 2 ; i ++ )
06521                                                 {
06522                                                         for ( j = -1 ; j < 2 ; j ++ )
06523                                                         {
06524                                                                 for ( k = -1 ; k < 2 ; k ++ )
06525                                                                 {
06526                                                                         if ( i != 0 && j != 0 && k != 0 )
06527                                                                         {
06528                                                                                 continue ;
06529                                                                         }
06530                                                                         int nx = ele->x + i ;
06531                                                                         int ny = ele->y + j ;
06532                                                                         int nz = ele->z + k ;
06533                                                                         // if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteSheet( nx,ny,nz, facevol ) )
06534                                                                         if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteSheet( nx,ny,nz ) )
06535                                                                         {
06536                                                                                 critical = 1 ;
06537                                                                                 break ;
06538                                                                         }
06539                                                                 }
06540                                                                 if ( critical )
06541                                                                 {
06542                                                                         break ;
06543                                                                 }
06544                                                         }
06545                                                         if ( critical )
06546                                                         {
06547                                                                 break ;
06548                                                         }
06549                                                 }
06550 
06551                                                 if ( critical )
06552                                                 {
06553                                                         setDataAt( ele->x, ele->y, ele->z, d ) ;
06554                                                         ele = queues[ d ]->getNext() ;
06555                                                 }
06556                                                 else
06557                                                 {
06558                                                         setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
06559                                                         if ( d > 2 )
06560                                                         {
06561                                                                 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
06562                                                         }
06563                                                         ele = queues[ d ]->remove() ;
06564                                                         num ++ ;
06565                                                 }
06566 
06567                                         }
06568                                         #ifdef VERBOSE
06569                                         printf("Non-minimal: %d\n", num) ;
06570                                         #endif
06571 
06572                                         if ( num == 0 )
06573                                         {
06574                                                 break ;
06575                                         }
06576                                 }
06577 
06578 
06579                                 /* End of debugging */
06580 
06581                                 /*
06582                                 queue3->reset() ;
06583                                 ele = queue3->getNext() ;
06584                                 while ( ele != NULL )
06585                                 {
06586                                         setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
06587                                         ele = queue3->remove() ;
06588                                 }
06589                                 */
06590                         }
06591 
06592 
06593                         // Finally, threshold the volume
06594                         #ifdef VERBOSE
06595                         //printf("Thresholding the volume to 0/1...\n") ;
06596                         #endif
06597                         //threshold( -1, 1, 0, 0 ) ;
06598                         threshold( 0, 0, 1 ) ;
06599 
06600                         delete facevol ;
06601                         delete fvol ;
06602                         delete queue2;
06603                         delete queue3;
06604                         for (d = -dis ; d >= 2 ; d -- ) {
06605                                 delete queues[d];
06606                         }
06607                         delete [] queues;
06608 
06609                         return - dis - 1 ;
06610                 }

int Volume::erodeSheet  ) 
 

Definition at line 6234 of file volume.cpp.

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

06235                 {
06236                         return erodeSheet( 3 ) ;
06237                 }

EMData * Volume::get_emdata  ) 
 

Definition at line 31 of file volume.cpp.

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

Referenced by EMAN::BinarySkeletonizerProcessor::process().

00032                 {
00033                         return this->getVolumeData()->get_emdata();
00034                 }

double Volume::getDataAt int  index  ) 
 

Definition at line 65 of file volume.cpp.

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

00065                                                     {
00066                         return volData->GetDataAt(index);
00067                 }

double Volume::getDataAt int  x,
int  y,
int  z
 

Definition at line 60 of file volume.cpp.

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

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::CleanUpSkeleton(), countIntEuler(), curveSkeleton(), curveSkeleton2D(), erodeHelix(), erodeSheet(), getNumNeighbor6(), getNumPotComplex(), hasCell(), hasCompleteHelix(), hasCompleteSheet(), isFeatureFace(), isHelixEnd(), isPiercable(), isSheetEnd(), isSimple(), markCellFace(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::MarkSurfaces(), skeleton(), surfaceSkeletonPres(), threshold(), and wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::VoxelOr().

00061                 {
00062                         return volData->GetDataAt(x, y, z);
00063                 }

int Volume::getIndex int  x,
int  y,
int  z
 

Definition at line 48 of file volume.cpp.

References wustl_mm::SkeletonMaker::VolumeData::GetIndex(), volData, x, and y.

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

00048                                                         {
00049                         return volData->GetIndex(x, y, z);
00050                 }

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().

00661                                                                     {
00662                         int rvalue = 0 ;
00663                         for ( int i = 0 ; i < 6 ; i ++ ) {
00664                                 int nx = ox + neighbor6[i][0] ;
00665                                 int ny = oy + neighbor6[i][1] ;
00666                                 int nz = oz + neighbor6[i][2] ;
00667                                 if ( getDataAt( nx, ny, nz ) >= 0 ) {
00668                                         rvalue ++ ;
00669                                 }
00670                         }
00671 
00672                         return rvalue ;
00673                 }

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().

02549                 {
02550                         //return 0 ;
02551 
02552                         int i;
02553                         double val = getDataAt( ox, oy, oz ) ;
02554                         if ( val <= 0 )
02555                         {
02556                 //              return 70 ;
02557                         }
02558 
02559                         // return getNumNeighbor6( ox, oy, oz ) ;
02560 
02561                         // int v = ((getNumNeighbor6( ox, oy, oz ) & 255) << 24) ;
02562                         //int v = 0  ;
02563 
02564                         int rvalue = 0, nx, ny, nz ;
02565                         setDataAt( ox, oy, oz, -val ) ;
02566 
02567                         /*
02568                         for ( i = -1 ; i < 2 ; i ++ )
02569                                 for ( j = -1 ; j < 2 ; j ++ )
02570                                         for ( k = -1 ; k < 2 ; k ++ )
02571                                         {
02572                                                 nx = ox + i ;
02573                                                 ny = oy + j ;
02574                                                 nz = oz + k ;
02575                                                 if ( getDataAt( nx, ny, nz ) == val )
02576                                                 {
02577                                                         if ( isSheetEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
02578                                                         {
02579                                                                 rvalue ++ ;
02580                                                         }
02581                                                 }
02582                                         }
02583                         */
02584 
02585                         for ( i = 0 ; i < 6 ; i ++ )
02586                         {
02587                                 nx = ox + neighbor6[i][0] ;
02588                                 ny = oy + neighbor6[i][1] ;
02589                                 nz = oz + neighbor6[i][2] ;
02590 
02591                                 if ( getDataAt(nx,ny,nz) >= 0 )
02592                                 {
02593                                         int num = getNumNeighbor6( nx, ny, nz ) ;
02594                                         if ( num > rvalue )
02595                                         {
02596                                                 rvalue = num ;
02597                                         }
02598                                 }
02599                         }
02600 
02601 
02602                         setDataAt( ox, oy, oz, val ) ;
02603 
02604                         return rvalue + getNumNeighbor6( ox, oy, oz ) * 10 ;
02605                         /*
02606                         int v = (((rvalue + getNumNeighbor6( ox, oy, oz ) * 10) & 255) << 24) ;
02607                         v |= ( ( ox & 255 ) << 16 )  ;
02608                         v |= ( ( oy & 255 ) << 8 ) ;
02609                         v |= ( ( oz & 255 ) ) ;
02610                         return v ;
02611                         */
02612 
02613                 }

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().

02616                 {
02617                         return getNumPotComplex( ox, oy, oz ) ;
02618 
02619                         //int i, j, k ;
02620                         //double val = getDataAt( ox, oy, oz ) ;
02621                         //if ( val <= 0 )
02622                         //{
02623                         //      return 0 ;
02624                         //}
02625 
02626                         //int rvalue = 0, nx, ny, nz ;
02627                         //setDataAt( ox, oy, oz, -val ) ;
02628 
02629                         //for ( i = -1 ; i < 2 ; i ++ )
02630                         //      for ( j = -1 ; j < 2 ; j ++ )
02631                         //              for ( k = -1 ; k < 2 ; k ++ )
02632                         //              {
02633                         //                      nx = ox + i ;
02634                         //                      ny = oy + j ;
02635                         //                      nz = oz + k ;
02636                         //                      if ( getDataAt( nx, ny, nz ) == val )
02637                         //                      {
02638                         //                              if ( isHelixEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
02639                         //                              {
02640                         //                                      rvalue ++ ;
02641                         //                              }
02642                         //                      }
02643                         //              }
02644 
02645                         //setDataAt( ox, oy, oz, val ) ;
02646 
02647                         //return rvalue + getNumNeighbor6( ox, oy, oz ) * 30 ;
02648                 }

float Volume::getOriginX  ) 
 

Definition at line 93 of file volume.cpp.

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

00093                                          {
00094                         return volData->GetOriginX();
00095                 }

float Volume::getOriginY  ) 
 

Definition at line 97 of file volume.cpp.

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

00097                                          {
00098                         return volData->GetOriginY();
00099                 }

float Volume::getOriginZ  ) 
 

Definition at line 101 of file volume.cpp.

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

00101                                          {
00102                         return volData->GetOriginZ();
00103                 }

int Volume::getSizeX  ) 
 

Definition at line 36 of file volume.cpp.

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

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::CleanUpSkeleton(), curveSkeleton(), curveSkeleton2D(), erodeHelix(), erodeSheet(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning(), markCellFace(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::MarkSurfaces(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PerformPureJuSkeletonization(), skeleton(), surfaceSkeletonPres(), threshold(), and wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::VoxelOr().

00036                                      {
00037                         return volData->GetSizeX();
00038                 }

int Volume::getSizeY  ) 
 

Definition at line 40 of file volume.cpp.

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

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::CleanUpSkeleton(), curveSkeleton(), curveSkeleton2D(), erodeHelix(), erodeSheet(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning(), markCellFace(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::MarkSurfaces(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PerformPureJuSkeletonization(), skeleton(), surfaceSkeletonPres(), threshold(), and wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::VoxelOr().

00040                                      {
00041                         return volData->GetSizeY();
00042                 }

int Volume::getSizeZ  ) 
 

Definition at line 44 of file volume.cpp.

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

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::CleanUpSkeleton(), curveSkeleton(), curveSkeleton2D(), erodeHelix(), erodeSheet(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning(), markCellFace(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::MarkSurfaces(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PerformPureJuSkeletonization(), skeleton(), surfaceSkeletonPres(), threshold(), and wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::VoxelOr().

00044                                      {
00045                         return volData->GetSizeZ();
00046                 }

float Volume::getSpacingX  ) 
 

Definition at line 81 of file volume.cpp.

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

00081                                           {
00082                         return volData->GetSpacingX();
00083                 }

float Volume::getSpacingY  ) 
 

Definition at line 85 of file volume.cpp.

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

00085                                           {
00086                         return volData->GetSpacingY();
00087                 }

float Volume::getSpacingZ  ) 
 

Definition at line 89 of file volume.cpp.

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

00089                                           {
00090                         return volData->GetSpacingZ();
00091                 }

VolumeData * Volume::getVolumeData  ) 
 

Definition at line 69 of file volume.cpp.

Referenced by get_emdata(), EMAN::BinarySkeletonizerProcessor::process(), and Volume().

00069                                                    {
00070                         return volData;
00071                 }

int Volume::hasCell int  ox,
int  oy,
int  oz
 

Definition at line 1016 of file volume.cpp.

References getDataAt().

Referenced by hasCompleteSheet(), and markCellFace().

01016                                                             {
01017                         for ( int i = 0 ; i < 2 ; i ++ )
01018                                 for ( int j = 0 ; j < 2 ; j ++ )
01019                                         for ( int k = 0 ; k < 2 ; k ++ ) {
01020                                                 if ( getDataAt( ox + i, oy + j, oz + k ) < 0 ) {
01021                                                         return 0 ;
01022                                                 }
01023                                         }
01024                         return 1 ;
01025                 }

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.

01513                 {
01514 
01515                         int i ;
01516                         int c1 = 0;
01517                         int nx, ny, nz ;
01518                         int j ;
01519 
01520                         for ( i = 0 ; i < 6 ; i ++ )
01521                         {
01522                                 nx = ox + neighbor6[i][0] ;
01523                                 ny = oy + neighbor6[i][1] ;
01524                                 nz = oz + neighbor6[i][2] ;
01525                                 if ( getDataAt( nx, ny, nz ) >= 0 )
01526                                 {
01527                                         if ( i % 2 == 0 )
01528                                         {
01529                                                 nx = ox ;
01530                                                 ny = oy ;
01531                                                 nz = oz ;
01532                                         }
01533 
01534                                         int val = (int)fvol->getDataAt( nx, ny, nz) ;
01535                                         if ( (( val >> ( 2 * ( i / 2 ) ) ) & 1) == 0 )
01536                                         {
01537                                                 c1 ++ ;
01538                                                 j = i ;
01539                                         }
01540                                 }
01541 
01542                         }
01543 
01544                         if ( c1 > 1 )
01545                         {
01546                                 return 1 ;
01547                         }
01548 
01549                         return 0 ;
01550 
01551                         /*
01552                         ox = ox + neighbor6[j][0] ;
01553                         oy = oy + neighbor6[j][1] ;
01554                         oz = oz + neighbor6[j][2] ;
01555                         c1 = 0 ;
01556                         for ( i = 0 ; i < 6 ; i ++ )
01557                         {
01558                                 nx = ox + neighbor6[i][0] ;
01559                                 ny = oy + neighbor6[i][1] ;
01560                                 nz = oz + neighbor6[i][2] ;
01561                                 if ( getDataAt( nx, ny, nz ) >= 0 )
01562                                 {
01563                                         c1 ++ ;
01564                                 }
01565 
01566                         }
01567 
01568                         if ( c1 > 1 )
01569                         {
01570                                 return 0 ;
01571                         }
01572                         else
01573                         {
01574                                 return 1 ;
01575                         }
01576                         */
01577                 }

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().

01457                 {
01458                         // Returns 1 if it has a complete helix
01459                         int i ;
01460                         int c1 = 0;
01461                         int nx, ny, nz ;
01462                         int j ;
01463 
01464                         for ( i = 0 ; i < 6 ; i ++ )
01465                         {
01466                                 nx = ox + neighbor6[i][0] ;
01467                                 ny = oy + neighbor6[i][1] ;
01468                                 nz = oz + neighbor6[i][2] ;
01469                                 if ( getDataAt( nx, ny, nz ) >= 0 )
01470                                 {
01471                                         c1 ++ ;
01472                                         j = i ;
01473                                 }
01474 
01475                         }
01476 
01477                         if ( c1 > 1 ) // || c1 == 0 )
01478                         {
01479                                 return 1 ;
01480                         }
01481 
01482                         return 0 ;
01483 
01484                         /*
01485                         ox = ox + neighbor6[j][0] ;
01486                         oy = oy + neighbor6[j][1] ;
01487                         oz = oz + neighbor6[j][2] ;
01488                         c1 = 0 ;
01489                         for ( i = 0 ; i < 6 ; i ++ )
01490                         {
01491                                 nx = ox + neighbor6[i][0] ;
01492                                 ny = oy + neighbor6[i][1] ;
01493                                 nz = oz + neighbor6[i][2] ;
01494                                 if ( getDataAt( nx, ny, nz ) >= 0 )
01495                                 {
01496                                         c1 ++ ;
01497                                 }
01498 
01499                         }
01500 
01501                         if ( c1 > 1 )
01502                         {
01503                                 return 0 ;
01504                         }
01505                         else
01506                         {
01507                                 return 1 ;
01508                         }
01509                         */
01510                 }

int Volume::hasCompleteSheet int  ox,
int  oy,
int  oz
 

Definition at line 1322 of file volume.cpp.

References countIntEuler().

01322                                                                      {
01323                         // Returns 1 if it lies in the middle of a sheet
01324                         int temp = countIntEuler( ox, oy, oz ) ;
01325                         if ( temp > 0 )
01326                         {
01327                                 return 1 ;
01328                         }
01329                         else
01330                         {
01331                                 return 0 ;
01332                         }
01333                 }

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().

01174                                                                                    {
01175                         int i, j, k ;
01176                         int nx, ny, nz ;
01177 
01178                         int edge[6] = { 0,0,0,0,0,0 } ;
01179                         int faceflag[ 12 ] ;
01180                         int tot = 0 ;
01181                         int cellflag[ 8 ] ;
01182 
01183                         int ct = 0 ;
01184                         for (  i = -1 ; i < 1 ; i ++ )
01185                                 for (  j = -1 ; j < 1 ; j ++ )
01186                                         for (  k = -1 ; k < 1 ; k ++ )
01187                                         {
01188                                                 if ( hasCell( ox + i, oy + j, oz + k ) )
01189                                                 {
01190                                                         cellflag[ ct ] = 1 ;
01191                                                 }
01192                                                 else
01193                                                 {
01194                                                         cellflag[ ct ] = 0 ;
01195                                                 }
01196                                                 ct ++ ;
01197                                         }
01198 
01199                         for ( i = 0 ; i < 12 ; i ++ )
01200                         {
01201                                 faceflag[ i ] = 1 ;
01202                                 for ( j = 0 ; j < 4 ; j ++ )
01203                                 {
01204                                         nx = ox + sheetNeighbor[i][j][0] ;
01205                                         ny = oy + sheetNeighbor[i][j][1] ;
01206                                         nz = oz + sheetNeighbor[i][j][2] ;
01207 
01208                                         if ( getDataAt( nx, ny, nz ) < 0 )
01209                                         {
01210                                                 faceflag[ i ] = 0 ;
01211                                                 break ;
01212                                         }
01213                                 }
01214 
01215                                 if ( faceflag[ i ] )
01216                                 {
01217                                         if ( cellflag[ faceCells[i][0] ] ^ cellflag[ faceCells[i][1] ] )
01218                                         {
01219                                                 int v1 = (int)( fvol->getDataAt(
01220                                                         ox - 1 + (( faceCells[i][0] >> 2 ) & 1 ),
01221                                                         oy - 1 + (( faceCells[i][0] >> 1 ) & 1 ),
01222                                                         oz - 1 + (( faceCells[i][0] ) & 1)) ) ;
01223                                                 int v2 = (int)( fvol->getDataAt(
01224                                                         ox - 1 + (( faceCells[i][1] >> 2 ) & 1 ),
01225                                                         oy - 1 + (( faceCells[i][1] >> 1 ) & 1 ),
01226                                                         oz - 1 + (( faceCells[i][1] ) & 1)) ) ;
01227                                                 if ( ((v1 >> (2 * (2 - i/4))) & 1) ||
01228                                                          ((v2 >> (2 * (2 - i/4) + 1 )) & 1) )
01229                                                 {
01230                                                         faceflag[ i ] = 0 ;
01231                                                 }
01232                                         }
01233                                 }
01234 
01235                                 if ( faceflag[ i ] )
01236                                 {
01237                                         int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
01238                                         edge[ e0 ] ++ ;
01239                                         edge[ e1 ] ++ ;
01240                                         tot ++ ;
01241                                 }
01242                         }
01243 
01244                         // Removing 1s
01245                         int numones = 0 ;
01246                         for ( i = 0 ; i < 6 ; i ++ )
01247                         {
01248                                 if ( edge[ i ] == 1 )
01249                                 {
01250                                         numones ++ ;
01251                                 }
01252                         }
01253                         while( numones > 0 )
01254                         {
01255                                 int e = 0;
01256                                 for ( i = 0 ; i < 6 ; i ++ )
01257                                 {
01258                                         if ( edge[ i ] == 1 )
01259                                         {
01260                                                 e = i ;
01261                                                 break ;
01262                                         }
01263                                 }
01264                                 /*
01265                                 if ( edge[ e ] != 1 )
01266                                 {
01267                                         printf("Wrong Again!********\n") ;
01268                                 }
01269                                 */
01270 
01271                                 int f, e2 ;
01272                                 for ( j = 0 ; j < 4 ; j ++ )
01273                                 {
01274                                         f = edgeFaces[ e ][ j ] ;
01275                                         if ( faceflag[ f ] )
01276                                         {
01277                                                 break ;
01278                                         }
01279                                 }
01280 
01281                                 /*
01282                                 if ( faceflag[ f ] == 0 )
01283                                 {
01284                                         printf("Wrong!********\n") ;
01285                                 }
01286                                 */
01287 
01288                                 if ( faceEdges[ f ][ 0 ] == e )
01289                                 {
01290                                         e2 = faceEdges[ f ][ 1 ] ;
01291                                 }
01292                                 else
01293                                 {
01294                                         e2 = faceEdges[ f ][ 0 ] ;
01295                                 }
01296 
01297 
01298                                 edge[ e ] -- ;
01299                                 numones -- ;
01300                                 edge[ e2 ] -- ;
01301                                 faceflag[ f ] = 0 ;
01302                                 tot -- ;
01303 
01304                                 if ( edge[ e2 ] == 1 )
01305                                 {
01306                                         numones ++ ;
01307                                 }
01308                                 else if ( edge[ e2 ] == 0 )
01309                                 {
01310                                         numones -- ;
01311                                 }
01312                         }
01313 
01314                         if ( tot > 0 )
01315                         {
01316                                 return 1 ;
01317                         }
01318 
01319                         return 0 ;
01320                 }

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().

02136                 {
02137                         // return 1 ;
02138 
02139                         int i, j ;
02140                         int nx, ny, nz ;
02141 
02142                         int faces = 12 ;
02143                         for ( i = 0 ; i < 12 ; i ++ )
02144                         {
02145                                 int ct = 0 ;
02146                                 for ( j = 0 ; j < 4 ; j ++ )
02147                                 {
02148                                         nx = ox + sheetNeighbor[i][j][0] ;
02149                                         ny = oy + sheetNeighbor[i][j][1] ;
02150                                         nz = oz + sheetNeighbor[i][j][2] ;
02151 
02152                                         if ( getDataAt( nx, ny, nz ) < 0 )
02153                                         {
02154                                                 ct = -1 ;
02155                                                 break ;
02156                                         }
02157                                         else if ( getNumNeighbor6( nx, ny, nz ) == 6 )
02158                                         {
02159                                                 ct = -1 ;
02160                                                 break ;
02161                                         }
02162         //                              else if ( getDataAt( nx, ny, nz ) == 0 )
02163         //                              {
02164         //                                      ct ++ ;
02165         //                              }
02166 
02167 
02168                                 }
02169                                 if ( ct == -1 || ct >= 1 )
02170                                 {
02171                                         faces -- ;
02172                                 }
02173                         }
02174 
02175                         if ( faces > 0 )
02176                         {
02177                                 return 1 ;
02178                         }
02179                         return 0 ;
02180 
02181                 }

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.

01789                                                                {
01790 
01791                         int i ;
01792                         int c1 = 0 , c2 = 0 ;
01793                         int nx, ny, nz ;
01794 
01795                         for ( i = 0 ; i < 6 ; i ++ )
01796                         {
01797                                 nx = ox + neighbor6[i][0] ;
01798                                 ny = oy + neighbor6[i][1] ;
01799                                 nz = oz + neighbor6[i][2] ;
01800 
01801                                 double val = getDataAt( nx, ny, nz ) ;
01802 
01803                                 if ( val >= 0 )
01804                                 {
01805                                         c1 ++ ;
01806                                         if ( getNumNeighbor6(nx,ny,nz) < 6 ) // if ( val > 0 && val < MAX_ERODE )
01807                                         {
01808                                                 c2 ++ ;
01809                                         }
01810                                 }
01811 
01812                         }
01813 
01814                         if ( c1 == 1 && c2 == 1 )
01815                         {
01816                                 return 1 ;
01817                         }
01818 
01819                         return 0 ;
01820                 }

int Volume::isHelixEnd int  ox,
int  oy,
int  oz,
Volume nvol
 

Definition at line 1579 of file volume.cpp.

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

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

01579                                                                              {
01580                         // Returns 1 if it is a curve endpoint
01581                         int i ;
01582                         int c1 = 0 , c2 = 0 ;
01583                         int nx, ny, nz ;
01584 
01585                         for ( i = 0 ; i < 6 ; i ++ )
01586                         {
01587                                 nx = ox + neighbor6[i][0] ;
01588                                 ny = oy + neighbor6[i][1] ;
01589                                 nz = oz + neighbor6[i][2] ;
01590 
01591                                 double val = getDataAt( nx, ny, nz ) ;
01592 
01593                                 if ( val >= 0 )
01594                                 {
01595                                         c1 ++ ;
01596                                         if ( val > 0 && val < MAX_ERODE && nvol->getDataAt( nx, ny, nz ) == 0 )
01597                                         {
01598                                                 c2 ++ ;
01599                                         }
01600                                 }
01601 
01602                         }
01603 
01604                         if ( c1 == 1 && c2 == 1 )
01605                         {
01606                                 return 1 ;
01607                         }
01608 
01609                         return 0 ;
01610                 }

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().

02382                 {
02383                         /* Test if this is a simple voxel */
02384                         // int flag = 0 ;
02385                         double vox[3][3][3] ;
02386 
02387                         int i, j, k ;
02388                         for ( i = -1 ; i < 2 ; i ++ )
02389                                 for ( j = -1 ; j < 2 ; j ++ )
02390                                         for ( k = -1 ; k < 2 ; k ++ )
02391                                         {
02392                                                 double tval = getDataAt( ox + i, oy + j, oz + k ) ;
02393 
02394                                                 /*
02395                                                 if ( tval == 0 || tval > (va )
02396                                                 {
02397                                                         flag = 1 ;
02398                                                 }
02399                                                 */
02400                                                 /*
02401                                                 if ( tval < 0 && tval == - curwid )
02402                                                 {
02403                                                 printf("Here %d", (int)-tval) ;
02404                                                 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
02405                                                 }
02406                                                 else
02407                                                 */
02408                                                 {
02409                                                         vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
02410                                                 }
02411                                         }
02412 
02413                                 /* Debugging
02414                                 printf("{") ;
02415                                 for ( i = 0 ; i < 3 ; i ++ )
02416                                 {
02417                                         if ( i ) printf(",") ;
02418                                         printf("{") ;
02419                                         for ( j = 0 ; j < 3 ; j ++ )
02420                                         {
02421                                                 if ( j ) printf(",") ;
02422                                                 printf("{") ;
02423                                                 for ( k = 0 ; k < 3 ; k ++ )
02424                                                 {
02425                                                         if ( k ) printf(",") ;
02426                                                         printf("%d", (vox[i][j][k] >=0 ? 1: 0));
02427                                                 }
02428                                                 printf("}") ;
02429                                         }
02430                                         printf("}") ;
02431                                 }
02432                                 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
02433                                 */
02434 
02435                         if ( countInt( vox ) == 1 && countExt( vox ) != 1 )
02436                         {
02437                                 return 1 ;
02438                         }
02439                         else
02440                         {
02441                                 return 0 ;
02442                         }
02443                 }

int Volume::isSheetEnd int  ox,
int  oy,
int  oz
 

*

Definition at line 2238 of file volume.cpp.

References hasCompleteSheet(), and isFeatureFace().

02239                 {
02240                         return ( ( hasCompleteSheet(ox,oy,oz) == 0 ) && isFeatureFace(ox,oy,oz) ) ;
02241 
02243                         //
02244                         //int i, j, k ;
02245                         //int nx, ny, nz ;
02246 
02247                         //double vox[3][3][3] ;
02248                         //for ( i = -1 ; i < 2 ; i ++ )
02249                         //      for ( j = -1 ; j < 2 ; j ++ )
02250                         //              for ( k = -1 ; k < 2 ; k ++ )
02251                         //              {
02252                         //                      vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ;
02253                         //              }
02254 
02255                         //int edge[6] = { 4,4,4,4,4,4 } ;
02256                         //int edge2[6] = { 4,4,4,4,4,4 } ;
02257 
02258                         //for ( i = 0 ; i < 12 ; i ++ )
02259                         //{
02260                         //      for ( j = 0 ; j < 4 ; j ++ )
02261                         //      {
02262                         //              nx = 1 + sheetNeighbor[i][j][0] ;
02263                         //              ny = 1 + sheetNeighbor[i][j][1] ;
02264                         //              nz = 1 + sheetNeighbor[i][j][2] ;
02265 
02266                         //              if ( vox[nx][ny][nz] < 0 )
02267                         //              {
02268                         //                      edge[ faceEdges[ i ][ 0 ] ] -- ;
02269                         //                      edge[ faceEdges[ i ][ 1 ] ] -- ;
02270                         //                      break ;
02271                         //              }
02272                         //      }
02273 
02274                         //      for ( j = 0 ; j < 4 ; j ++ )
02275                         //      {
02276                         //              nx = 1 + sheetNeighbor[i][j][0] ;
02277                         //              ny = 1 + sheetNeighbor[i][j][1] ;
02278                         //              nz = 1 + sheetNeighbor[i][j][2] ;
02279 
02280                         //              if ( vox[nx][ny][nz] <= 0 )
02281                         //              {
02282                         //                      edge2[ faceEdges[ i ][ 0 ] ] -- ;
02283                         //                      edge2[ faceEdges[ i ][ 1 ] ] -- ;
02284                         //                      break ;
02285                         //              }
02286                         //      }
02287                         //}
02288 
02289                         //
02291                         //for ( i = 0 ; i < 6 ; i ++ )
02292                         //{
02293                         //      nx = 1 + neighbor6[i][0] ;
02294                         //      ny = 1 + neighbor6[i][1] ;
02295                         //      nz = 1 + neighbor6[i][2] ;
02296                         //      if ( edge[i] == 0 && vox[nx][ny][nz] >= 0 )
02297                         //      {
02298                         //              return 0 ;
02299                         //      }
02300                         //}
02301                         //*/
02302                         //
02303 
02304 
02305                         //for ( i = 0 ; i < 6 ; i ++ )
02306                         //{
02307                         //      if ( edge[ i ] == 1 && edge2[ i ] == 1 )
02308                         //      {
02309                         //              return 1 ;
02310                         //      }
02311                         //}
02312 
02313                         //return 0 ;
02314                 }

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().

01822                 {
01823                         // Returns 1 if it contains a sheet boundary. Noise-resistant
01824                         int i, j ;
01825                         int nx, ny, nz ;
01826 
01827                         int edge[6] = { 0,0,0,0,0,0 } ;
01828                         int faceflag[ 12 ] ;
01829                         int hasFeatureFace = 0 ;
01830                         int tot = 0 ;
01831 
01832                         for ( i = 0 ; i < 12 ; i ++ )
01833                         {
01834                                 faceflag[ i ] = 1 ;
01835                                 int hasFeature = 1 ;
01836 
01837                                 for ( j = 0 ; j < 4 ; j ++ )
01838                                 {
01839                                         nx = ox + sheetNeighbor[i][j][0] ;
01840                                         ny = oy + sheetNeighbor[i][j][1] ;
01841                                         nz = oz + sheetNeighbor[i][j][2] ;
01842 
01843                                         if ( getDataAt( nx, ny, nz ) == 0 || nvol->getDataAt( nx, ny, nz ) == 1 )
01844                                         {
01845                                                 hasFeature = 0 ;
01846                                         }
01847                                         if ( getDataAt( nx, ny, nz ) < 0 )
01848                                         {
01849                                                 faceflag[ i ] = 0 ;
01850                                                 break ;
01851                                         }
01852                                 }
01853                                 if ( faceflag[ i ] == 1 && hasFeature )
01854                                 {
01855                                         hasFeatureFace ++ ;
01856                                         // return 0 ;
01857                                 }
01858 
01859                                 if ( faceflag[ i ] )
01860                                 {
01861                                         int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
01862                                         edge[ e0 ] ++ ;
01863                                         edge[ e1 ] ++ ;
01864                                         tot ++ ;
01865                                 }
01866                         }
01867 
01868                         if ( tot == 0 || hasFeatureFace == 0 )
01869                         {
01870                                 return 0 ;
01871                         }
01872 
01873                         // Removing 1s
01874                         int numones = 0 ;
01875                         for ( i = 0 ; i < 6 ; i ++ )
01876                         {
01877                                 if ( edge[ i ] == 1 )
01878                                 {
01879                                         numones ++ ;
01880                                 }
01881                         }
01882                         while( numones > 0 )
01883                         {
01884                                 int e = 0;
01885                                 for ( i = 0 ; i < 6 ; i ++ )
01886                                 {
01887                                         if ( edge[ i ] == 1 )
01888                                         {
01889                                                 e = i ;
01890                                                 break ;
01891                                         }
01892                                 }
01893                                 /*
01894                                 if ( edge[ e ] != 1 )
01895                                 {
01896                                         printf("Wrong Again!********\n") ;
01897                                 }
01898                                 */
01899 
01900                                 int f, e2 ;
01901                                 for ( j = 0 ; j < 4 ; j ++ )
01902                                 {
01903                                         f = edgeFaces[ e ][ j ] ;
01904                                         if ( faceflag[ f ] )
01905                                         {
01906                                                 break ;
01907                                         }
01908                                 }
01909 
01910                                 /*
01911                                 if ( faceflag[ f ] == 0 )
01912                                 {
01913                                         printf("Wrong!********\n") ;
01914                                 }
01915                                 */
01916 
01917                                 if ( faceEdges[ f ][ 0 ] == e )
01918                                 {
01919                                         e2 = faceEdges[ f ][ 1 ] ;
01920                                 }
01921                                 else
01922                                 {
01923                                         e2 = faceEdges[ f ][ 0 ] ;
01924                                 }
01925 
01926 
01927                                 edge[ e ] -- ;
01928                                 numones -- ;
01929                                 edge[ e2 ] -- ;
01930                                 faceflag[ f ] = 0 ;
01931                                 tot -- ;
01932 
01933                                 if ( edge[ e2 ] == 1 )
01934                                 {
01935                                         numones ++ ;
01936                                 }
01937                                 else if ( edge[ e2 ] == 0 )
01938                                 {
01939                                         numones -- ;
01940                                 }
01941                         }
01942 
01943                         if ( tot > 0 )
01944                         {
01945                                 return 0 ;
01946                         }
01947 
01948                         return 1 ;
01949                 }

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().

02317                 {
02318                         /* Test if this is a simple voxel */
02319                         // int flag = 0 ;
02320                         double vox[3][3][3] ;
02321 
02322                         int i, j, k ;
02323                         for ( i = -1 ; i < 2 ; i ++ )
02324                                 for ( j = -1 ; j < 2 ; j ++ )
02325                                         for ( k = -1 ; k < 2 ; k ++ )
02326                                         {
02327                                                 double tval = getDataAt( ox + i, oy + j, oz + k ) ;
02328 
02329                                                 /*
02330                                                 if ( tval == 0 || tval > (va )
02331                                                 {
02332                                                         flag = 1 ;
02333                                                 }
02334                                                 */
02335                                                 /*
02336                                                 if ( tval < 0 && tval == - curwid )
02337                                                 {
02338                                                 printf("Here %d", (int)-tval) ;
02339                                                 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
02340                                                 }
02341                                                 else
02342                                                 */
02343                                                 {
02344                                                         vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
02345                                                 }
02346                                         }
02347 
02348                                 /* Debugging
02349                                 printf("{") ;
02350                                 for ( i = 0 ; i < 3 ; i ++ )
02351                                 {
02352                                         if ( i ) printf(",") ;
02353                                         printf("{") ;
02354                                         for ( j = 0 ; j < 3 ; j ++ )
02355                                         {
02356                                                 if ( j ) printf(",") ;
02357                                                 printf("{") ;
02358                                                 for ( k = 0 ; k < 3 ; k ++ )
02359                                                 {
02360                                                         if ( k ) printf(",") ;
02361                                                         printf("%d", (vox[i][j][k] >=0 ? 1: 0));
02362                                                 }
02363                                                 printf("}") ;
02364                                         }
02365                                         printf("}") ;
02366                                 }
02367                                 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
02368                                 */
02369 
02370                         if ( countInt( vox ) == 1 && countExt( vox ) == 1 )
02371                         {
02372                                 return 1 ;
02373                         }
02374                         else
02375                         {
02376                                 return 0 ;
02377                         }
02378                 }

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().

01027                                                {
01028                         int i, j, k ;
01029                         Volume* fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
01030 
01031                         //return fvol ;
01032 
01033                         for ( i = 0 ; i < getSizeX() ; i ++ )
01034                                 for ( j = 0 ; j < getSizeY() ; j ++ )
01035                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
01036                                         {
01037                                                 if( getDataAt(i,j,k) >= 0 )
01038                                                 {
01039                                                         if ( hasCell(i,j,k) )
01040                                                         {
01041                                                                 for ( int m = 0 ; m < 6 ; m ++ )
01042                                                                 {
01043                                                                         int nx = i + neighbor6[m][0] ;
01044                                                                         int ny = j + neighbor6[m][1] ;
01045                                                                         int nz = k + neighbor6[m][2] ;
01046                                                                         if ( ! hasCell(nx,ny,nz) )
01047                                                                         {
01048                                                                                 fvol->setDataAt(i,j,k,(double)(1<<m)) ;
01049                                                                                 break ;
01050                                                                         }
01051                                                                 }
01052                                                         }
01053                                                 }
01054                                         }
01055 
01056 
01057                         return fvol ;
01058                 }

void Volume::pad int  padBy,
double  padValue
 

* 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().

00273                                                             {
00274                         volData->Pad(padBy, padValue);
00275                 }

void Volume::setDataAt int  index,
double  d
 

Definition at line 56 of file volume.cpp.

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

00056                                                             {
00057                         volData->SetDataAt(index, (float)d);
00058                 }

void Volume::setDataAt int  x,
int  y,
int  z,
double  d
 

Definition at line 52 of file volume.cpp.

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

Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::CleanUpSkeleton(), curveSkeleton(), curveSkeleton2D(), erodeHelix(), erodeSheet(), getNumPotComplex(), markCellFace(), wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::MarkSurfaces(), skeleton(), surfaceSkeletonPres(), threshold(), and wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::VoxelOr().

00052                                                                       {
00053                         volData->SetDataAt(x, y, z, (float)d);
00054                 }

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.

00077                                                                          {
00078                         volData->SetOrigin(orgX, orgY, orgZ);
00079                 }

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.

00073                                                                         {
00074                         volData->SetSpacing(spx, spy, spz);
00075                 }

void Volume::skeleton float  thr,
Volume svol,
Volume hvol
 

* Adding ends *// Move its neighboring unvisited node to queue2

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, 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.

05693                 {
05694                         int i, j, k ;
05695                         // First, threshold the volume
05696                         #ifdef VERBOSE
05697                         printf("Thresholding the volume to -1/0...\n") ;
05698                         #endif
05699                         threshold( thr, -1, 0 ) ;
05700 
05701                         // Next, apply convergent erosion
05702                         // by preserving: complex nodes, curve end-points, and sheet points
05703 
05704                         // Next, initialize the linked queue
05705                         #ifdef VERBOSE
05706                         printf("Initializing queue...\n") ;
05707                         #endif
05708                         GridQueue2* queue2 = new GridQueue2( ) ;
05709                         GridQueue2* queue3 = new GridQueue2( ) ;
05710                         PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
05711 
05712                         for ( i = 0 ; i < getSizeX() ; i ++ )
05713                                 for ( j = 0 ; j < getSizeY() ; j ++ )
05714                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
05715                                         {
05716                                                 if ( getDataAt( i, j, k ) >= 0 )
05717                                                 {
05718                                                         if ( svol->getDataAt(i,j,k) > 0 || hvol->getDataAt(i,j,k) > 0 )
05719                                                         {
05720                                                                 setDataAt( i, j, k, MAX_ERODE ) ;
05721                                                         }
05722                                                         else
05723                                                         {
05724                                                                 for ( int m = 0 ; m < 6 ; m ++ )
05725                                                                 {
05726                                                                         if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
05727                                                                         {
05728                                                                                 setDataAt( i, j, k, 1 ) ;
05729                                                                                 queue2->prepend( i, j, k ) ;
05730                                                                                 break ;
05731                                                                         }
05732                                                                 }
05733                                                         }
05734                                                 }
05735                                         }
05736                         int wid = MAX_ERODE ;
05737                         #ifdef VERBOSE
05738                         printf("Total %d nodes\n", queue2->getNumElements() ) ;
05739 
05740 
05741                         // Perform erosion
05742                         printf("Start erosion to %d...\n", wid) ;
05743                         #endif
05744                         gridQueueEle* ele ;
05745                         gridPoint* gp ;
05746                         int ox, oy, oz ;
05747                         int score ;
05748                         Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
05749                         for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
05750                         {
05751                                 scrvol->setDataAt( i, -1 ) ;
05752                         }
05753 
05754 
05755                         for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
05756                         {
05757                                 // At the start of each iteration,
05758                                 // queue2 holds all the nodes for this layer
05759                                 // queue3 and queue are empty
05760 
05761                                 int numComplex = 0, numSimple = 0 ;
05762                                 #ifdef VERBOSE
05763                                 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
05764                                 #endif
05765 
05766 
05767                                 // Next,
05768                                 // Compute score for each node left in queue2
05769                                 // move them into priority queue
05770                                 queue2->reset() ;
05771                                 ele = queue2->getNext() ;
05772                                 while ( ele != NULL )
05773                                 {
05774                                         ox = ele->x ;
05775                                         oy = ele->y ;
05776                                         oz = ele->z ;
05777 
05778                                         // Compute score
05779                                         score = getNumPotComplex2( ox, oy, oz ) ;
05780                                         scrvol->setDataAt( ox, oy, oz, score ) ;
05781 
05782                                         // Push to queue
05783                                         gp = new gridPoint ;
05784                                         gp->x = ox ;
05785                                         gp->y = oy ;
05786                                         gp->z = oz ;
05787                                         // queue->add( gp, -score ) ;
05788                                         queue->add( gp, score ) ;
05789 
05790                                         ele = queue2->remove() ;
05791                                 }
05792 
05793                                 // Rename queue3 to be queue2,
05794                                 // Clear queue3
05795                                 // From now on, queue2 holds nodes of next level
05796                                 delete queue2 ;
05797                                 queue2 = queue3 ;
05798                                 queue3 = new GridQueue2( ) ;
05799 
05800                                 // Next, start priority queue iteration
05801                                 while ( ! queue->isEmpty() )
05802                                 {
05803                                         // Retrieve the node with the highest score
05804                                         queue->remove( gp, score ) ;
05805                                         ox = gp->x ;
05806                                         oy = gp->y ;
05807                                         oz = gp->z ;
05808                                         delete gp ;
05809         //                              score = -score ;
05810 
05811                                         // Ignore the node
05812                                         // if it has been processed before
05813                                         // or it has an updated score
05814                                         if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
05815                                         {
05816                                                 continue ;
05817                                         }
05818 
05819                                         /* Added for debugging */
05820                                         // Check simple
05821                                         if ( ! isSimple( ox, oy, oz ) )
05822                                         {
05823                                                 // Complex, set to next layer
05824                                                 setDataAt( ox, oy, oz, curwid + 1 ) ;
05825                                                 queue2->prepend( ox, oy, oz ) ;
05826                                                 numComplex ++ ;
05827                                         }
05828                                         else
05829                                         {
05830                                                 setDataAt( ox, oy, oz, -1 ) ;
05831                                                 numSimple ++ ;
05832                                         }
05833                                         /* Adding ends */
05834 
05835                                         // Move its neighboring unvisited node to queue2
05836                                         for ( int m = 0 ; m < 6 ; m ++ )
05837                                         {
05838                                                 int nx = ox + neighbor6[m][0] ;
05839                                                 int ny = oy + neighbor6[m][1] ;
05840                                                 int nz = oz + neighbor6[m][2] ;
05841                                                 if ( getDataAt( nx, ny, nz ) == 0 )
05842                                                 {
05843                                                         setDataAt( nx, ny, nz, curwid + 1 ) ;
05844                                                         queue2->prepend( nx, ny, nz ) ;
05845                                                 }
05846                                         }
05847 
05848                                         // Update scores for nodes in its 5x5 neighborhood
05849                                         // insert them back into priority queue
05850 
05851                                         for ( i = -2 ; i < 3 ;i ++ )
05852                                                 for ( j = -2 ; j < 3 ; j ++ )
05853                                                         for ( k = -2 ; k < 3 ; k ++ )
05854                                                         {
05855                                                                 int nx = ox + i ;
05856                                                                 int ny = oy + j ;
05857                                                                 int nz = oz + k ;
05858 
05859                                                                 if ( getDataAt( nx, ny, nz ) == curwid )
05860                                                                 {
05861                                                                         // Compute score
05862                                                                         score = getNumPotComplex2( nx, ny, nz ) ;
05863 
05864                                                                         if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
05865                                                                         {
05866                                                                                 // printf("Update\n") ;
05867                                                                                 scrvol->setDataAt( nx, ny, nz, score ) ;
05868                                                                                 // Push to queue
05869                                                                                 gp = new gridPoint ;
05870                                                                                 gp->x = nx ;
05871                                                                                 gp->y = ny ;
05872                                                                                 gp->z = nz ;
05873                                                                                 // queue->add( gp, -score ) ;
05874                                                                                 queue->add( gp, score ) ;
05875                                                                         }
05876                                                                 }
05877                                                         }
05878 
05879 
05880                                 }
05881 
05882                                 #ifdef VERBOSE
05883                                 printf("%d complex, %d simple\n", numComplex, numSimple) ;
05884                                 #endif
05885 
05886                                 if ( numSimple == 0 )
05887                                 {
05888                                         delete queue2 ;
05889                                         break ;
05890                                 }
05891                         }
05892 
05893                         // Finally, clean up
05894                         #ifdef VERBOSE
05895                         printf("Thresholding the volume to 0/1...\n") ;
05896                         #endif
05897                         threshold( 0, 0, 1 ) ;
05898                         delete scrvol;
05899                         delete queue;
05900                         delete queue3;
05901                 }

void Volume::skeleton float  thr,
int  off
 

Definition at line 5089 of file volume.cpp.

References getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getSizeX(), getSizeY(), getSizeZ(), isSimple(), wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), 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().

05090                 {
05091                         int i, j, k ;
05092                         // First, threshold the volume
05093                         #ifdef VERBOSE
05094                         printf("Thresholding the volume to -1/0...\n") ;
05095                         #endif
05096                         threshold( thr, -1, 0 ) ;
05097 
05098                         // Next, apply convergent erosion
05099                         // by preserving: complex nodes, curve end-points, and sheet points
05100 
05101                         // Next, initialize the linked queue
05102                         #ifdef VERBOSE
05103                         printf("Initializing queue...\n") ;
05104                         #endif
05105                         GridQueue2* queue2 = new GridQueue2( ) ;
05106                         GridQueue2* queue = new GridQueue2( ) ;
05107 
05108                         for ( i = 0 ; i < getSizeX() ; i ++ )
05109                                 for ( j = 0 ; j < getSizeY() ; j ++ )
05110                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
05111                                         {
05112                                                 if ( getDataAt( i, j, k ) >= 0 )
05113                                                 {
05114                                                         {
05115                                                                 for ( int m = 0 ; m < 6 ; m ++ )
05116                                                                 {
05117                                                                         if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
05118                                                                         {
05119                                                                                 setDataAt( i, j, k, 1 ) ;
05120                                                                                 queue2->prepend( i, j, k ) ;
05121                                                                                 break ;
05122                                                                         }
05123                                                                 }
05124                                                         }
05125                                                 }
05126                                         }
05127                         int wid = 0 ;
05128                         #ifdef VERBOSE
05129                         printf("Total %d nodes\n", queue2->getNumElements() ) ;
05130                         printf("Start erosion to %d...\n", wid) ;
05131                         #endif
05132 
05133                         // Perform erosion
05134                         gridQueueEle* ele ;
05135                         int ox, oy, oz ;
05136 
05137                         while( 1 )
05138                         {
05139                                 wid ++ ;
05140 
05141                                 // At the start of each iteration,
05142                                 // queue2 holds all the nodes for this layer
05143                                 // queue is empty
05144 
05145                                 int numComplex = 0, numSimple = 0 ;
05146                                 #ifdef VERBOSE
05147                                 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), wid) ;
05148                                 #endif
05149 
05150                                 // Rename queue2 to be queue,
05151                                 // Clear queue2
05152                                 // From now on, queue2 holds nodes of next level
05153                                 delete queue ;
05154                                 queue = queue2 ;
05155                                 queue2 = new GridQueue2( ) ;
05156 
05157                                 // Next, start queue iteration
05158                                 queue->reset() ;
05159                                 ele = queue->getNext();
05160                                 while ( ele != NULL )
05161                                 {
05162                                         ox = ele->x ;
05163                                         oy = ele->y ;
05164                                         oz = ele->z ;
05165         //                              delete ele ;
05166 
05167                                         // Check simple
05168                                         if ( ! isSimple( ox, oy, oz ) )
05169                                         {
05170                                                 // Complex, set to next layer
05171                                                 queue2->prepend( ox, oy, oz ) ;
05172                                                 numComplex ++ ;
05173                                         }
05174                                         /*
05175                                         else if ( ox == off || oy == off || oz == off ||
05176                                                 ox == getSizeX() - off - 1 || oy == getSizeY() - off - 1 || oz == getSizeZ() - off - 1 )
05177                                         {
05178                                                 // Wall, don't erode, set to next layer
05179                                                 queue2->prepend( ox, oy, oz ) ;
05180                                                 numComplex ++ ;
05181                                         }
05182                                         */
05183                                         else
05184                                         {
05185                                                 setDataAt( ox, oy, oz, -1 ) ;
05186                                                 numSimple ++ ;
05187 
05188                                                 for ( int m = 0 ; m < 6 ; m ++ )
05189                                                 {
05190                                                         int nx = ox + neighbor6[m][0] ;
05191                                                         int ny = oy + neighbor6[m][1] ;
05192                                                         int nz = oz + neighbor6[m][2] ;
05193                                                         if ( getDataAt( nx, ny, nz ) == 0 )
05194                                                         {
05195                                                                 setDataAt( nx, ny, nz, 1 ) ;
05196                                                                 queue2->prepend( nx, ny, nz ) ;
05197                                                         }
05198                                                 }
05199 
05200                                         }
05201 
05202                                         ele = queue->remove() ;
05203                                 }
05204                                 #ifdef VERBOSE
05205                                 printf("Level %d: %d complex, %d simple\n", wid, numComplex, numSimple) ;
05206                                 #endif
05207 
05208                                 if ( numSimple == 0 )
05209                                 {
05210                                         break ;
05211                                 }
05212                         }
05213 
05214                         // Finally, clean up
05215                         delete queue;
05216                         delete queue2;
05217                         #ifdef VERBOSE
05218                         printf("Thresholding the volume to 0/1...\n") ;
05219                         #endif
05220                         threshold( 0, 0, 1 ) ;
05221                 }

*void Volume::surfaceSkeletonPres float  thr,
Volume preserve
 

* Commented for debugging/ Find complex nodes in its 3x3 neighborhood/ move them to queue2

Parameters:
noise: 0 for no-noise-reduction, n for level-n noise reduction

Definition at line 8672 of file volume.cpp.

References EMAN::PriorityQueue< ValueT, KeyT >::add(), flag, 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, 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().

08673                 {
08674                         int i, j, k ;
08675                         // First, threshold the volume
08676                         #ifdef VERBOSE
08677                         printf("Thresholding the volume to -MAX_ERODE/0...\n") ;
08678                         #endif
08679                         threshold( thr, -MAX_ERODE, 0 ) ;
08680 
08681                         // Next, initialize the linked queue
08682                         #ifdef VERBOSE
08683                         printf("Initializing queue...\n") ;
08684                         #endif
08685                         GridQueue2* queue2 = new GridQueue2( ) ;
08686                         GridQueue2* queue3 = new GridQueue2( ) ;
08687                         GridQueue2* queue4 = new GridQueue2( ) ;
08688 
08689                         PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
08690 
08691                         for ( i = 0 ; i < getSizeX() ; i ++ )
08692                                 for ( j = 0 ; j < getSizeY() ; j ++ )
08693                                         for ( k = 0 ; k < getSizeZ() ; k ++ )
08694                                         {
08695                                                 if ( getDataAt( i, j, k ) >= 0 ) {
08696                                                         if(preserve->getDataAt(i, j, k) > 0) {
08697                                                                 setDataAt(i, j, k, MAX_ERODE);
08698                                                         } else {
08699                                                                 for ( int m = 0 ; m < 6 ; m ++ )
08700                                                                 {
08701                                                                         if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
08702                                                                         {
08703                                                                                 // setDataAt( i, j, k, 1 ) ;
08704                                                                                 queue2->prepend( i, j, k ) ;
08705                                                                                 break ;
08706                                                                         }
08707                                                                 }
08708                                                         }
08709                                                 }
08710                                         }
08711                         int wid = MAX_ERODE ;
08712                         #ifdef VERBOSE
08713                         printf("Total %d nodes\n", queue2->getNumElements() ) ;
08714                         printf("Start erosion to %d...\n", wid) ;
08715                         #endif
08716 
08717 
08718                         // Perform erosion
08719                         gridQueueEle* ele ;
08720                         gridPoint* gp ;
08721                         int ox, oy, oz ;
08722                         int score ;
08723                         Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
08724                         for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
08725                         {
08726                                 scrvol->setDataAt( i, -1 ) ;
08727                         }
08728 
08729         #ifdef  NOISE_DIS_SHEET
08730                         Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
08731         #endif
08732 
08733                         for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
08734                         {
08735                                 // At the start of each iteration,
08736                                 // queue2 and queue4 holds all the nodes for this layer
08737                                 // queue3 and queue are empty
08738 
08739                                 int numComplex = 0, numSimple = 0 ;
08740                                 #ifdef VERBOSE
08741                                 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
08742                                 #endif
08743 
08744                                 /*
08745                                 We first need to assign curwid + 1 to every node in this layer
08746                                 */
08747                                 queue2->reset() ;
08748                                 ele = queue2->getNext() ;
08749                                 while ( ele != NULL )
08750                                 {
08751                                         ox = ele->x ;
08752                                         oy = ele->y ;
08753                                         oz = ele->z ;
08754 
08755                                         if ( getDataAt(ox,oy,oz) == curwid )
08756                                         {
08757                                                 ele = queue2->remove() ;
08758                                         }
08759                                         else
08760                                         {
08761                                                 setDataAt(ox,oy,oz, curwid) ;
08762                                                 ele = queue2->getNext() ;
08763                                         }
08764                                 }
08765                                 queue4->reset() ;
08766                                 ele = queue4->getNext() ;
08767                                 while ( ele != NULL )
08768                                 {
08769                                         ox = ele->x ;
08770                                         oy = ele->y ;
08771                                         oz = ele->z ;
08772 
08773                                         queue2->prepend(ox,oy,oz) ;
08774                                         ele = queue4->remove() ;
08775                                 }
08776 
08777                                 // Now queue2 holds all the nodes for this layer
08778 
08779         #ifdef NOISE_DIS_SHEET
08780                                 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
08781                                 queue2->reset() ;
08782 
08783                                 // First run
08784                                 int flag = 0 ;
08785                                 while ( ( ele = queue2->getNext() ) != NULL )
08786                                 {
08787                                         ox = ele->x ;
08788                                         oy = ele->y ;
08789                                         oz = ele->z ;
08790                                         if ( NOISE_DIS_SHEET <= 1 )
08791                                         {
08792                                                 noisevol->setDataAt( ox, oy, oz, 0 ) ;
08793                                         }
08794                                         else
08795                                         {
08796                                                 flag = 0 ;
08797                                                 for ( int m = 0 ; m < 6 ; m ++ )
08798                                                 {
08799                                                         int nx = ox + neighbor6[m][0] ;
08800                                                         int ny = oy + neighbor6[m][1] ;
08801                                                         int nz = oz + neighbor6[m][2] ;
08802                                                         if ( getDataAt( nx, ny, nz ) == 0 )
08803                                                         {
08804                                                                 noisevol->setDataAt( ox, oy, oz, 1 ) ;
08805                                                                 flag = 1 ;
08806                                                                 break ;
08807                                                         }
08808                                                 }
08809                                                 if ( ! flag )
08810                                                 {
08811                                                         noisevol->setDataAt( ox, oy, oz, 0 ) ;
08812                                                 }
08813                                         }
08814                                 }
08815 
08816                                 for ( int cur = 1 ; cur < NOISE_DIS_SHEET ; cur ++ )
08817                                 {
08818                                         queue2->reset() ;
08819                                         int count = 0 ;
08820 
08821                                         while ( ( ele = queue2->getNext() ) != NULL )
08822                                         {
08823                                                 ox = ele->x ;
08824                                                 oy = ele->y ;
08825                                                 oz = ele->z ;
08826 
08827                                                 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
08828                                                 {
08829                                                         continue ;
08830                                                 }
08831 
08832                                                 flag = 0 ;
08833                                                 for ( int m = 0 ; m < 6 ; m ++ )
08834                                                 {
08835                                                         int nx = ox + neighbor6[m][0] ;
08836                                                         int ny = oy + neighbor6[m][1] ;
08837                                                         int nz = oz + neighbor6[m][2] ;
08838                                                         if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
08839                                                         {
08840                                                                 noisevol->setDataAt( ox, oy, oz, 1 ) ;
08841                                                                 count ++ ;
08842                                                                 break ;
08843                                                         }
08844                                                 }
08845                                         }
08846 
08847                                         if ( count == 0 )
08848                                         {
08849                                                 break ;
08850                                         }
08851                                 }
08852 
08853 
08854         #endif
08855 
08856                                 /* Commented for debugging
08857 
08858                                 // First,
08859                                 // check for complex nodes in queue2
08860                                 // move them from queue2 to queue3
08861                                 queue2->reset() ;
08862                                 ele = queue2->getNext() ;
08863                                 while ( ele != NULL )
08864                                 {
08865                                         ox = ele->x ;
08866                                         oy = ele->y ;
08867                                         oz = ele->z ;
08868 
08869                                         // Check simple
08870         #ifndef NOISE_DIS_SHEET
08871                                         if ( isSheetEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
08872         #else
08873                                         if ( isSheetEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
08874         #endif
08875                                         {
08876                                                 // Complex, set to next layer
08877                                                 setDataAt( ox, oy, oz, curwid + 1 ) ;
08878                                                 queue3->prepend( ox, oy, oz ) ;
08879                                                 ele = queue2->remove() ;
08880 
08881                                                 numComplex ++ ;
08882                                         }
08883                                         else
08884                                         {
08885                                                 ele = queue2->getNext() ;
08886                                         }
08887                                 }
08888                                 */
08889 
08890 
08891                                 // Next,
08892                                 // Compute score for each node left in queue2
08893                                 // move them into priority queue
08894                                 queue2->reset() ;
08895                                 ele = queue2->getNext() ;
08896                                 while ( ele != NULL )
08897                                 {
08898                                         ox = ele->x ;
08899                                         oy = ele->y ;
08900                                         oz = ele->z ;
08901 
08902                                         // Compute score
08903                                         score = getNumPotComplex( ox, oy, oz ) ;
08904                                         scrvol->setDataAt( ox, oy, oz, score ) ;
08905 
08906                                         // Push to queue
08907                                         gp = new gridPoint ;
08908                                         gp->x = ox ;
08909                                         gp->y = oy ;
08910                                         gp->z = oz ;
08911                                         // queue->add( gp, -score ) ;
08912                                         queue->add( gp, score ) ;
08913 
08914                                         ele = queue2->remove() ;
08915                                 }
08916 
08917                                 // Rename queue3 to be queue2,
08918                                 // Clear queue3
08919                                 // From now on, queue2 holds nodes of next level
08920                                 delete queue2 ;
08921                                 queue2 = queue3 ;
08922                                 queue3 = new GridQueue2( ) ;
08923 
08924 
08925                                 // Next, start priority queue iteration
08926                                 while ( ! queue->isEmpty() )
08927                                 {
08928                                         // Retrieve the node with the highest score
08929                                         queue->remove( gp, score ) ;
08930                                         ox = gp->x ;
08931                                         oy = gp->y ;
08932                                         oz = gp->z ;
08933                                         delete gp ;
08934                                         // printf("%d\n", score);
08935         //                              score = -score ;
08936 
08937                                         // Ignore the node
08938                                         // if it has been processed before
08939                                         // or it has an updated score
08940                                         if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
08941                                         {
08942                                                 continue ;
08943                                         }
08944 
08945                                         /* Commented for debugging
08946 
08947                                         // Remove this simple node
08948                                         setDataAt( ox, oy, oz, -1 ) ;
08949                                         numSimple ++ ;
08950                                         // printf("Highest score: %d\n", score) ;
08951                                         */
08952 
08953                                         /* Added for debugging */
08954                                         // Check simple
08955         #ifndef NOISE_DIS_SHEET
08956                                         // if ( hasFeatureFace( ox, oy, oz ) )
08957                                         if ( (! isSimple( ox, oy, oz )) || isSheetEnd( ox, oy, oz ) )
08958                                         // if ( hasIsolatedFace(ox,oy,oz)  && (! isNoiseSheetEnd(ox,oy,oz)))
08959         #else
08960                                         // if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz, noisevol ) )
08961                                         if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz, noisevol ) || isHelixEnd( ox, oy, oz, noisevol ))
08962                                         // if ( isBertrandEndPoint( ox, oy, oz ) )
08963         #endif
08964                                         {
08965                                                 // Complex, set to next layer
08966                                                 setDataAt( ox, oy, oz, curwid + 1 ) ;
08967                                                 queue4->prepend( ox, oy, oz ) ;
08968                                                 numComplex ++ ;
08969 
08970                                         }
08971                                         else
08972                                         {
08973                                                 setDataAt( ox, oy, oz, -1 ) ;
08974                                                 numSimple ++ ;
08975 
08976                                         }
08977                                         /* Adding ends */
08978 
08979                                         // Move its neighboring unvisited node to queue2
08980                                         for ( int m = 0 ; m < 6 ; m ++ )
08981                                         {
08982                                                 int nx = ox + neighbor6[m][0] ;
08983                                                 int ny = oy + neighbor6[m][1] ;
08984                                                 int nz = oz + neighbor6[m][2] ;
08985                                                 if ( getDataAt( nx, ny, nz ) == 0 )
08986                                                 {
08987                                                         // setDataAt( nx, ny, nz, curwid + 1 ) ;
08988                                                         queue2->prepend( nx, ny, nz ) ;
08989                                                 }
08990                                         }
08991 
08992                                         /* Commented for debugging
08993 
08994                                         // Find complex nodes in its 3x3 neighborhood
08995                                         // move them to queue2
08996                                         for ( i = -1 ; i < 2 ; i ++ )
08997                                                 for ( j = -1 ; j < 2 ; j ++ )
08998                                                         for ( k = -1 ; k < 2 ; k ++ )
08999                                                         {
09000                                                                 int nx = ox + i ;
09001                                                                 int ny = oy + j ;
09002                                                                 int nz = oz + k ;
09003 
09004                                                                 // Check simple
09005                                                                 if ( getDataAt( nx, ny, nz ) == curwid &&
09006                                                                         // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
09007         #ifndef NOISE_DIS_SHEET
09008                                                                         ( isSheetEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
09009         #else
09010                                                                         ( isSheetEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
09011         #endif
09012 
09013                                                                 {
09014                                                                         // Complex, set to next layer
09015                                                                         setDataAt( nx, ny, nz, curwid + 1 ) ;
09016                                                                         queue2->prepend( nx, ny, nz ) ;
09017                                                                         numComplex ++ ;
09018                                                                 }
09019                                                         }
09020                                         */
09021 
09022                                         // Update scores for nodes in its 5x5 neighborhood
09023                                         // insert them back into priority queue
09024 
09025                                         for ( i = -2 ; i < 3 ;i ++ )
09026                                                 for ( j = -2 ; j < 3 ; j ++ )
09027                                                         for ( k = -2 ; k < 3 ; k ++ )
09028                                                         {
09029                                                                 int nx = ox + i ;
09030                                                                 int ny = oy + j ;
09031                                                                 int nz = oz + k ;
09032 
09033                                                                 if ( getDataAt( nx, ny, nz ) == curwid )
09034                                                                 {
09035                                                                         // Compute score
09036                                                                         score = getNumPotComplex( nx, ny, nz ) ;
09037 
09038                                                                         if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
09039                                                                         {
09040                                                                                 // printf("Update\n") ;
09041                                                                                 scrvol->setDataAt( nx, ny, nz, score ) ;
09042                                                                                 // Push to queue
09043                                                                                 gp = new gridPoint ;
09044                                                                                 gp->x = nx ;
09045                                                                                 gp->y = ny ;
09046                                                                                 gp->z = nz ;
09047                                                                                 // queue->add( gp, -score ) ;
09048                                                                                 queue->add( gp, score ) ;
09049                                                                         }
09050                                                                 }
09051                                                         }
09052 
09053 
09054                                 }
09055 
09056                                 #ifdef VERBOSE
09057                                 printf("%d complex, %d simple\n", numComplex, numSimple) ;
09058                                 #endif
09059 
09060                                 if ( numSimple == 0 )
09061                                 {
09062                                                 break ;
09063                                 }
09064                         }
09065 
09066                         // Finally, clean up
09067                         #ifdef VERBOSE
09068                         printf("Thresholding the volume to 0/1...\n") ;
09069                         #endif
09070                         threshold( 0, 0, 1 ) ;
09071 
09072                         delete scrvol;
09073                         delete queue;
09074                         delete queue2;
09075                         delete queue3;
09076                         delete queue4;
09077 
09078                 }

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(), in, and setDataAt().

09531                 {
09532                         float val;
09533                         for ( int i = 0 ; i < getSizeX() ; i ++ )
09534                                 for ( int j = 0 ; j < getSizeY() ; j ++ )
09535                                         for ( int k = 0 ; k < getSizeZ() ; k ++ )
09536                                         {
09537                                                 val = (float)getDataAt(i, j, k);
09538                                                 if(markBoundary) {
09539                                                         if ( i > 1 && i < getSizeX() - 2 && j > 1 && j < getSizeY() - 2 && k > 1 && k < getSizeZ() - 2 ) {
09540                                                                 if(val < thr) {
09541                                                                         setDataAt(i, j, k, out);
09542                                                                 } else {
09543                                                                         setDataAt(i, j, k, in);
09544                                                                 }
09545                                                         }
09546                                                         else
09547                                                         {
09548                                                                 setDataAt(i, j, k, boundary);
09549                                                         }
09550                                                 } else {
09551                                                         if(val < thr) {
09552                                                                 setDataAt(i, j, k, out);
09553                                                         } else {
09554                                                                 setDataAt(i, j, k, in);
09555                                                         }
09556                                                 }
09557                                         }
09558                 }

void Volume::threshold double  thr,
int  out,
int  in,
int  boundary
 

Definition at line 9525 of file volume.cpp.

References in, and threshold().

09526                 {
09527                         threshold(thr, out, in, boundary, true);
09528                 }

void Volume::threshold double  thr,
int  out,
int  in
 

Definition at line 9520 of file volume.cpp.

References in, and threshold().

09521                 {
09522                         threshold( thr, out, in, out, true) ;
09523                 }

void Volume::threshold double  thr  ) 
 

Normalize to a given range.

Definition at line 9515 of file volume.cpp.

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

09516                 {
09517                         threshold( thr, 0, 1, 0, true) ;
09518                 }


Member Data Documentation

VolumeData* wustl_mm::SkeletonMaker::Volume::volData [private]
 

Definition at line 222 of file volume.h.

Referenced by getDataAt(), getIndex(), getOriginX(), getOriginY(), getOriginZ(), getSizeX(), getSizeY(), getSizeZ(), getSpacingX(), getSpacingY(), getSpacingZ(), pad(), setDataAt(), setOrigin(), setSpacing(), and Volume().


The documentation for this class was generated from the following files:
Generated on Tue Jun 11 13:49:55 2013 for EMAN2 by  doxygen 1.3.9.1