EMAN2
Public Member Functions | Private Attributes | List of all members
wustl_mm::SkeletonMaker::Volume Class Reference

#include <volume.h>

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

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 ()
 
EMData * get_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)
 
int getNumNeighbor6 (int ox, int oy, int oz)
 
int hasCell (int ox, int oy, int oz)
 
Volume * markCellFace ()
 
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)
 
More...
 
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])
 
More...
 
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)
 insert them back into priority queue More...
 
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)
 insert them back into priority queue More...
 
void erodeHelix ()
 
void erodeHelix (int disthr)
 
int erodeSheet ()
 
int erodeSheet (int disthr)
 
void surfaceSkeletonPres (float thr, Volume *preserve)
 for ( int m = 0 ; m < 6 ; m ++ ) More...
 
void threshold (double thr)
 
More...
 
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)
 
VolumeData * getVolumeData ()
 

Private Attributes

VolumeData * volData
 

Detailed Description

Definition at line 70 of file volume.h.

Constructor & Destructor Documentation

◆ Volume() [1/4]

Volume::Volume ( EMData *  em)

Definition at line 9 of file volume.cpp.

10 {
11 this->volData = new VolumeData(em);
12 }
wustl_mm::SkeletonMaker::VolumeData
Definition: volume_data.h:12
wustl_mm::SkeletonMaker::Volume::volData
VolumeData * volData
Definition: volume.h:222

References volData.

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

◆ Volume() [2/4]

Volume::Volume ( int  x,
int  y,
int  z 
)

Definition at line 13 of file volume.cpp.

13 {
14 volData = new VolumeData(x, y, z);
15 }
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517

References volData, x, and y.

◆ Volume() [3/4]

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

Definition at line 17 of file volume.cpp.

17 {
18 volData = new VolumeData(x, y, z, val);
19 }

References volData, x, and y.

◆ Volume() [4/4]

Volume::Volume ( int  x,
int  y,
int  z,
int  offx,
int  offy,
int  offz,
Volume *  vol 
)

Definition at line 21 of file volume.cpp.

21 {
22 volData = new VolumeData(x, y, z, offx, offy, offz, vol->getVolumeData());
23 }
wustl_mm::SkeletonMaker::Volume::getVolumeData
VolumeData * getVolumeData()
Definition: volume.cpp:69

References getVolumeData(), volData, x, and y.

◆ ~Volume()

Volume::~Volume ( )

Definition at line 26 of file volume.cpp.

27 {
28 delete volData;
29 }

References volData.

Member Function Documentation

◆ components26()

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

Definition at line 2767 of file volume.cpp.

2768 {
2769 // Stupid flood fill
2770 int tot = 0 ;
2771 int queue[27][3] ;
2772 int vis[3][3][3] ;
2773 int head = 0, tail = 1 ;
2774 int i, j, k ;
2775 for ( i = 0 ; i < 3 ; i ++ )
2776 for ( j = 0 ; j < 3 ; j ++ )
2777 for ( k = 0 ; k < 3 ; k ++ )
2778 {
2779 vis[i][j][k] = 0 ;
2780 if ( vox[i][j][k] )
2781 {
2782 if ( tot == 0 )
2783 {
2784 queue[0][0] = i ;
2785 queue[0][1] = j ;
2786 queue[0][2] = k ;
2787 vis[i][j][k] = 1 ;
2788 }
2789 tot ++ ;
2790 }
2791 }
2792 if ( tot == 0 )
2793 {
2794 return 0 ;
2795 }
2796
2797 int ct = 1 ;
2798 while ( head != tail )
2799 {
2800 int x = queue[head][0] ;
2801 int y = queue[head][1] ;
2802 int z = queue[head][2] ;
2803 head ++ ;
2804
2805 for ( i = -1 ; i < 2 ; i ++ )
2806 for ( j = -1 ; j < 2 ; j ++ )
2807 for ( k = -1 ; k < 2 ; k ++ )
2808 {
2809 int nx = x + i ;
2810 int ny = y + j ;
2811 int nz = z + k ;
2812 if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
2813 {
2814 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
2815 {
2816 queue[tail][0] = nx ;
2817 queue[tail][1] = ny ;
2818 queue[tail][2] = nz ;
2819 tail ++ ;
2820 vis[nx][ny][nz] = 1 ;
2821 ct ++ ;
2822 }
2823 }
2824 }
2825 }
2826
2827 if ( ct == tot )
2828 {
2829 return 1 ;
2830 }
2831 else
2832 {
2833 return 2 ;
2834 }
2835
2836 }

References x, and y.

Referenced by countExt().

◆ components6()

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

Definition at line 2698 of file volume.cpp.

2699 {
2700 // Stupid flood fill
2701 int tot = 0 ;
2702 int queue[27][3] ;
2703 int vis[3][3][3] ;
2704 int head = 0, tail = 1 ;
2705 int i, j, k ;
2706 for ( i = 0 ; i < 3 ; i ++ )
2707 for ( j = 0 ; j < 3 ; j ++ )
2708 for ( k = 0 ; k < 3 ; k ++ )
2709 {
2710 vis[i][j][k] = 0 ;
2711 if ( vox[i][j][k] )
2712 {
2713 if ( tot == 0 )
2714 {
2715 queue[0][0] = i ;
2716 queue[0][1] = j ;
2717 queue[0][2] = k ;
2718 vis[i][j][k] = 1 ;
2719 }
2720 tot ++ ;
2721 }
2722 }
2723 if ( tot == 0 )
2724 {
2725 return 0 ;
2726 }
2727 // printf("total: %d\n", tot) ;
2728
2729 int ct = 1 ;
2730 while ( head != tail )
2731 {
2732 int x = queue[head][0] ;
2733 int y = queue[head][1] ;
2734 int z = queue[head][2] ;
2735 head ++ ;
2736
2737 for ( i = 0 ; i < 6 ; i ++ )
2738 {
2739 int nx = x + neighbor6[i][0] ;
2740 int ny = y + neighbor6[i][1] ;
2741 int nz = z + neighbor6[i][2] ;
2742 if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
2743 {
2744 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
2745 {
2746 queue[tail][0] = nx ;
2747 queue[tail][1] = ny ;
2748 queue[tail][2] = nz ;
2749 tail ++ ;
2750 vis[nx][ny][nz] = 1 ;
2751 ct ++ ;
2752 }
2753 }
2754 }
2755 }
2756
2757 if ( ct == tot )
2758 {
2759 return 1 ;
2760 }
2761 else
2762 {
2763 return 2 ;
2764 }
2765
2766 }
wustl_mm::SkeletonMaker::neighbor6
const int neighbor6[6][3]
Definition: volume.h:27

References wustl_mm::SkeletonMaker::neighbor6, x, and y.

Referenced by countInt(), and countIntEuler().

◆ countExt()

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

Definition at line 2838 of file volume.cpp.

2839 {
2840 int tvox[3][3][3] ;
2841
2842 for ( int i = 0 ; i < 3 ; i ++ )
2843 for ( int j = 0 ; j < 3 ; j ++ )
2844 for ( int k = 0 ; k < 3 ; k ++ )
2845 {
2846 if ( vox[i][j][k] < 0 )
2847 {
2848 tvox[i][j][k] = 1 ;
2849 }
2850 else
2851 {
2852 tvox[i][j][k] = 0 ;
2853 }
2854 }
2855
2856 return components26( tvox ) ;
2857 }
wustl_mm::SkeletonMaker::Volume::components26
int components26(int vox[3][3][3])
Definition: volume.cpp:2767

References components26().

Referenced by isPiercable(), and isSimple().

◆ countInt()

int Volume::countInt ( double  vox[3][3][3])

Definition at line 2859 of file volume.cpp.

2860 {
2861 int i, j, k ;
2862 int tvox[3][3][3] ;
2863
2864 for ( i = 0 ; i < 3 ; i ++ )
2865 for ( j = 0 ; j < 3 ; j ++ )
2866 for ( k = 0 ; k < 3 ; k ++ )
2867 {
2868 tvox[i][j][k] = 0 ;
2869 }
2870
2871 for ( i = 0 ; i < 6 ; i ++ )
2872 {
2873 int nx = 1 + neighbor6[i][0] ;
2874 int ny = 1 + neighbor6[i][1] ;
2875 int nz = 1 + neighbor6[i][2] ;
2876 if ( vox[ nx ][ ny ][ nz ] >= 0 )
2877 {
2878 tvox[ nx ][ ny ][ nz ] = 1 ;
2879 for ( j = 0 ; j < 4 ; j ++ )
2880 {
2881 int nnx = neighbor64[i][j][0] + nx ;
2882 int nny = neighbor64[i][j][1] + ny ;
2883 int nnz = neighbor64[i][j][2] + nz ;
2884 if ( vox[ nnx ][ nny ][ nnz ] >= 0 )
2885 {
2886 tvox[ nnx ][ nny ][ nnz ] = 1 ;
2887 }
2888 }
2889 }
2890 }
2891
2892 return components6( tvox ) ;
2893 }
wustl_mm::SkeletonMaker::Volume::components6
int components6(int vox[3][3][3])
Definition: volume.cpp:2698
wustl_mm::SkeletonMaker::neighbor64
const int neighbor64[6][4][3]
Definition: volume.h:29

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

Referenced by isPiercable(), and isSimple().

◆ countIntEuler()

int Volume::countIntEuler ( int  ox,
int  oy,
int  oz 
)

Definition at line 2895 of file volume.cpp.

2896 {
2897 int nv = 0 , ne = 0, nc = 0 ;
2898
2899 int i, j, k ;
2900 int tvox[3][3][3] ;
2901 double vox[3][3][3] ;
2902
2903 for ( i = 0 ; i < 3 ; i ++ )
2904 for ( j = 0 ; j < 3 ; j ++ )
2905 for ( k = 0 ; k < 3 ; k ++ )
2906 {
2907 vox[i][j][k] = getDataAt( ox - 1 + i, oy - 1 + j, oz - 1 + k ) ;
2908 tvox[i][j][k] = 0 ;
2909 }
2910
2911 for ( i = 0 ; i < 6 ; i ++ )
2912 {
2913 int nx = 1 + neighbor6[i][0] ;
2914 int ny = 1 + neighbor6[i][1] ;
2915 int nz = 1 + neighbor6[i][2] ;
2916 if ( vox[nx][ny][nz] >= 0 )
2917 {
2918 tvox[ nx ][ ny ][ nz ] = 1 ;
2919
2920 nv ++ ;
2921
2922 for ( j = 0 ; j < 4 ; j ++ )
2923 {
2924 int nnx = neighbor64[i][j][0] + nx ;
2925 int nny = neighbor64[i][j][1] + ny ;
2926 int nnz = neighbor64[i][j][2] + nz ;
2927 if ( vox[nnx][nny][nnz] >= 0 )
2928 {
2929 if ( tvox[nnx][nny][nnz] == 0 )
2930 {
2931 tvox[nnx][nny][nnz] = 1 ;
2932 nv ++ ;
2933 }
2934
2935 ne ++ ;
2936 }
2937 }
2938 }
2939 }
2940
2941 nc = components6( tvox ) ;
2942
2943 return ( nc - ( nv - ne ) ) ;
2944 }
wustl_mm::SkeletonMaker::Volume::getDataAt
double getDataAt(int x, int y, int z)
Definition: volume.cpp:60

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

Referenced by hasCompleteSheet().

◆ curveSkeleton() [1/2]

void Volume::curveSkeleton ( float  thr,
Volume *  svol 
)

Definition at line 4268 of file volume.cpp.

4269 {
4270 int i, j, k ;
4271 // First, threshold the volume
4272 #ifdef VERBOSE
4273 printf("Thresholding the volume to -1/0...\n") ;
4274 #endif
4275 threshold( thr, -1, 0 ) ;
4276
4277 // Next, apply convergent erosion
4278 // by preserving: complex nodes, curve end-points, and sheet points
4279
4280 // Next, initialize the linked queue
4281 #ifdef VERBOSE
4282 printf("Initializing queue...\n") ;
4283 #endif
4284 GridQueue2* queue2 = new GridQueue2( ) ;
4285 GridQueue2* queue3 = new GridQueue2( ) ;
4286 GridQueue2* queue4 = new GridQueue2( ) ;
4287 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
4288
4289 for ( i = 0 ; i < getSizeX() ; i ++ )
4290 for ( j = 0 ; j < getSizeY() ; j ++ )
4291 for ( k = 0 ; k < getSizeZ() ; k ++ )
4292 {
4293 if ( getDataAt( i, j, k ) >= 0 )
4294 {
4295 if ( svol->getDataAt(i,j,k) > 0 )
4296 {
4297 setDataAt( i, j, k, MAX_ERODE ) ;
4298 }
4299 else
4300 {
4301 for ( int m = 0 ; m < 6 ; m ++ )
4302 {
4303 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
4304 {
4305 // setDataAt( i, j, k, 1 ) ;
4306 queue2->prepend( i, j, k ) ;
4307 break ;
4308 }
4309 }
4310 }
4311 }
4312 }
4313
4314 int wid = MAX_ERODE ;
4315 #ifdef VERBOSE
4316 printf("Total %d nodes\n", queue2->getNumElements() ) ;
4317 printf("Start erosion to %d...\n", wid) ;
4318 #endif
4319
4320
4321 // Perform erosion
4322 gridQueueEle* ele ;
4323 gridPoint* gp ;
4324 int ox, oy, oz ;
4325 int score ;
4326 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
4327 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
4328 {
4329 scrvol->setDataAt( i, -1 ) ;
4330 }
4331
4332 #ifdef NOISE_DIS_HELIX
4333 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
4334 #endif
4335
4336 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
4337 {
4338 // At the start of each iteration,
4339 // queue2 holds all the nodes for this layer
4340 // queue3 and queue are empty
4341
4342 int numComplex = 0, numSimple = 0 ;
4343 #ifdef VERBOSE
4344 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
4345 #endif
4346
4347 /*
4348 We first need to assign curwid + 1 to every node in this layer
4349 */
4350 queue2->reset() ;
4351 ele = queue2->getNext() ;
4352 while ( ele != NULL )
4353 {
4354 ox = ele->x ;
4355 oy = ele->y ;
4356 oz = ele->z ;
4357
4358 if ( getDataAt(ox,oy,oz) == curwid )
4359 {
4360 ele = queue2->remove() ;
4361 }
4362 else
4363 {
4364 setDataAt(ox,oy,oz, curwid) ;
4365 ele = queue2->getNext() ;
4366 }
4367 }
4368 queue4->reset() ;
4369 ele = queue4->getNext() ;
4370 while ( ele != NULL )
4371 {
4372 ox = ele->x ;
4373 oy = ele->y ;
4374 oz = ele->z ;
4375
4376 queue2->prepend(ox,oy,oz) ;
4377 ele = queue4->remove() ;
4378 }
4379
4380 // Now queue2 holds all the nodes for this layer
4381
4382 #ifdef NOISE_DIS_HELIX
4383 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
4384 queue2->reset() ;
4385
4386 // First run
4387 int flag = 0 ;
4388 while ( ( ele = queue2->getNext() ) != NULL )
4389 {
4390 ox = ele->x ;
4391 oy = ele->y ;
4392 oz = ele->z ;
4393 if ( NOISE_DIS_HELIX <= 1 )
4394 {
4395 noisevol->setDataAt( ox, oy, oz, 0 ) ;
4396 }
4397 else
4398 {
4399 flag = 0 ;
4400 for ( int m = 0 ; m < 6 ; m ++ )
4401 {
4402 int nx = ox + neighbor6[m][0] ;
4403 int ny = oy + neighbor6[m][1] ;
4404 int nz = oz + neighbor6[m][2] ;
4405 if ( getDataAt( nx, ny, nz ) == 0 )
4406 {
4407 noisevol->setDataAt( ox, oy, oz, 1 ) ;
4408 flag = 1 ;
4409 break ;
4410 }
4411 }
4412 if ( ! flag )
4413 {
4414 noisevol->setDataAt( ox, oy, oz, 0 ) ;
4415 }
4416 }
4417 }
4418
4419 int cur, visited ;
4420 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
4421 {
4422 queue2->reset() ;
4423 int count = 0 ;
4424 visited = 0 ;
4425
4426 while ( ( ele = queue2->getNext() ) != NULL )
4427 {
4428 ox = ele->x ;
4429 oy = ele->y ;
4430 oz = ele->z ;
4431
4432 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
4433 {
4434 visited ++ ;
4435 continue ;
4436 }
4437
4438 flag = 0 ;
4439 for ( int m = 0 ; m < 6 ; m ++ )
4440 {
4441 int nx = ox + neighbor6[m][0] ;
4442 int ny = oy + neighbor6[m][1] ;
4443 int nz = oz + neighbor6[m][2] ;
4444 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
4445 {
4446 noisevol->setDataAt( ox, oy, oz, 1 ) ;
4447 visited ++ ;
4448 count ++ ;
4449 break ;
4450 }
4451 }
4452 }
4453
4454 if ( count == 0 )
4455 {
4456 break ;
4457 }
4458 }
4459 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
4460
4461
4462 #endif
4463 /* Commented out for debugging
4464
4465 // First,
4466 // check for complex nodes in queue2
4467 // move them from queue2 to queue3
4468 queue2->reset() ;
4469 ele = queue2->getNext() ;
4470 while ( ele != NULL )
4471 {
4472 ox = ele->x ;
4473 oy = ele->y ;
4474 oz = ele->z ;
4475
4476 // Check simple
4477 #ifndef NOISE_DIS_HELIX
4478 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4479 #else
4480 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
4481 #endif
4482 {
4483 // Complex, set to next layer
4484 setDataAt( ox, oy, oz, curwid + 1 ) ;
4485 queue3->prepend( ox, oy, oz ) ;
4486 ele = queue2->remove() ;
4487
4488 numComplex ++ ;
4489 }
4490 else
4491 {
4492 ele = queue2->getNext() ;
4493 }
4494 }
4495 */
4496
4497 // Next,
4498 // Compute score for each node left in queue2
4499 // move them into priority queue
4500 queue2->reset() ;
4501 ele = queue2->getNext() ;
4502 while ( ele != NULL )
4503 {
4504 ox = ele->x ;
4505 oy = ele->y ;
4506 oz = ele->z ;
4507
4508 // Compute score
4509 score = getNumPotComplex2( ox, oy, oz ) ;
4510 scrvol->setDataAt( ox, oy, oz, score ) ;
4511
4512 // Push to queue
4513 gp = new gridPoint ;
4514 gp->x = ox ;
4515 gp->y = oy ;
4516 gp->z = oz ;
4517 // queue->add( gp, -score ) ;
4518 queue->add( gp, score ) ;
4519
4520 ele = queue2->remove() ;
4521 }
4522
4523 // Rename queue3 to be queue2,
4524 // Clear queue3
4525 // From now on, queue2 holds nodes of next level
4526 delete queue2 ;
4527 queue2 = queue3 ;
4528 queue3 = new GridQueue2( ) ;
4529
4530 // Next, start priority queue iteration
4531 while ( ! queue->isEmpty() )
4532 {
4533 // Retrieve the node with the highest score
4534 queue->remove( gp, score ) ;
4535 ox = gp->x ;
4536 oy = gp->y ;
4537 oz = gp->z ;
4538 delete gp ;
4539 // score = -score ;
4540
4541 // Ignore the node
4542 // if it has been processed before
4543 // or it has an updated score
4544 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
4545 {
4546 continue ;
4547 }
4548
4549 /* Commented out for debugging
4550
4551 // Remove this simple node
4552 setDataAt( ox, oy, oz, -1 ) ;
4553 numSimple ++ ;
4554 // printf("Highest score: %d\n", score) ;
4555 */
4556
4557 /* Added for debugging */
4558 // Check simple
4559 #ifndef NOISE_DIS_HELIX
4560 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
4561 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4562 #else
4563 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4564 #endif
4565 {
4566 // Complex, set to next layer
4567 setDataAt( ox, oy, oz, curwid + 1 ) ;
4568 queue4->prepend( ox, oy, oz ) ;
4569 numComplex ++ ;
4570 }
4571 else
4572 {
4573 setDataAt( ox, oy, oz, -1 ) ;
4574 numSimple ++ ;
4575 }
4576 /* Adding ends */
4577
4578 // Move its neighboring unvisited node to queue2
4579 for ( int m = 0 ; m < 6 ; m ++ )
4580 {
4581 int nx = ox + neighbor6[m][0] ;
4582 int ny = oy + neighbor6[m][1] ;
4583 int nz = oz + neighbor6[m][2] ;
4584 if ( getDataAt( nx, ny, nz ) == 0 )
4585 {
4586 // setDataAt( nx, ny, nz, curwid + 1 ) ;
4587 queue2->prepend( nx, ny, nz ) ;
4588 }
4589 }
4590
4591 /* Commented out for debugging
4592
4593 // Find complex nodes in its 3x3 neighborhood
4594 // move them to queue2
4595 for ( i = -1 ; i < 2 ; i ++ )
4596 for ( j = -1 ; j < 2 ; j ++ )
4597 for ( k = -1 ; k < 2 ; k ++ )
4598 {
4599 int nx = ox + i ;
4600 int ny = oy + j ;
4601 int nz = oz + k ;
4602
4603 // Check simple
4604 if ( getDataAt( nx, ny, nz ) == curwid &&
4605 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
4606 #ifndef NOISE_DIS_HELIX
4607 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
4608 #else
4609 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
4610 #endif
4611
4612 {
4613 // Complex, set to next layer
4614 setDataAt( nx, ny, nz, curwid + 1 ) ;
4615 queue2->prepend( nx, ny, nz ) ;
4616 numComplex ++ ;
4617 }
4618 }
4619 */
4620
4621 // Update scores for nodes in its 5x5 neighborhood
4622 // insert them back into priority queue
4623
4624 for ( i = -2 ; i < 3 ;i ++ )
4625 for ( j = -2 ; j < 3 ; j ++ )
4626 for ( k = -2 ; k < 3 ; k ++ )
4627 {
4628 int nx = ox + i ;
4629 int ny = oy + j ;
4630 int nz = oz + k ;
4631
4632 if ( getDataAt( nx, ny, nz ) == curwid )
4633 {
4634 // Compute score
4635 score = getNumPotComplex2( nx, ny, nz ) ;
4636
4637 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
4638 {
4639 // printf("Update\n") ;
4640 scrvol->setDataAt( nx, ny, nz, score ) ;
4641 // Push to queue
4642 gp = new gridPoint ;
4643 gp->x = nx ;
4644 gp->y = ny ;
4645 gp->z = nz ;
4646 // queue->add( gp, -score ) ;
4647 queue->add( gp, score ) ;
4648 }
4649 }
4650 }
4651
4652
4653 }
4654
4655 #ifdef VERBOSE
4656 printf("%d complex, %d simple\n", numComplex, numSimple) ;
4657 #endif
4658
4659 if ( numSimple == 0 )
4660 {
4661 break ;
4662 }
4663 }
4664
4665 // Finally, clean up
4666 delete scrvol;
4667 delete queue;
4668 delete queue2;
4669 delete queue3;
4670 delete queue4;
4671 #ifdef VERBOSE
4672 printf("Thresholding the volume to 0/1...\n") ;
4673 #endif
4674 threshold( 0, 0, 1 ) ;
4675 }
EMAN::PriorityQueue
Template class for a priority queue.
Definition: priority_queue.h:18
EMAN::PriorityQueue::remove
void remove(ValueT *&v, KeyT &k)
Remove an element.
Definition: priority_queue.h:120
EMAN::PriorityQueue::add
void add(ValueT *v, KeyT k)
Add an element.
Definition: priority_queue.h:86
EMAN::PriorityQueue::isEmpty
bool isEmpty()
Test whether empty.
Definition: priority_queue.h:70
wustl_mm::SkeletonMaker::GridQueue2
Definition: grid_queue2.h:14
wustl_mm::SkeletonMaker::GridQueue2::getNumElements
int getNumElements()
Definition: grid_queue2.cpp:50
wustl_mm::SkeletonMaker::GridQueue2::reset
void reset()
Definition: grid_queue2.cpp:43
wustl_mm::SkeletonMaker::GridQueue2::remove
gridQueueEle * remove()
Definition: grid_queue2.cpp:70
wustl_mm::SkeletonMaker::GridQueue2::getNext
gridQueueEle * getNext()
Definition: grid_queue2.cpp:25
wustl_mm::SkeletonMaker::GridQueue2::prepend
void prepend(int xx, int yy, int zz)
Definition: grid_queue2.cpp:55
wustl_mm::SkeletonMaker::Volume
Definition: volume.h:70
wustl_mm::SkeletonMaker::Volume::isHelixEnd
int isHelixEnd(int ox, int oy, int oz, Volume *nvol)
Definition: volume.cpp:1579
wustl_mm::SkeletonMaker::Volume::isSimple
int isSimple(int ox, int oy, int oz)
Definition: volume.cpp:2316
wustl_mm::SkeletonMaker::Volume::getSizeZ
int getSizeZ()
Definition: volume.cpp:44
wustl_mm::SkeletonMaker::Volume::getSizeX
int getSizeX()
Definition: volume.cpp:36
wustl_mm::SkeletonMaker::Volume::Volume
Volume(EMData *em)
Definition: volume.cpp:9
wustl_mm::SkeletonMaker::Volume::getNumPotComplex2
int getNumPotComplex2(int ox, int oy, int oz)
Definition: volume.cpp:2615
wustl_mm::SkeletonMaker::Volume::getSizeY
int getSizeY()
Definition: volume.cpp:40
wustl_mm::SkeletonMaker::Volume::threshold
void threshold(double thr)
Definition: volume.cpp:9515
wustl_mm::SkeletonMaker::Volume::setDataAt
void setDataAt(int x, int y, int z, double d)
Definition: volume.cpp:52
wustl_mm::SkeletonMaker::gridPoint
Definition: volume.h:66
wustl_mm::SkeletonMaker::gridPoint::y
int y
Definition: volume.h:67
wustl_mm::SkeletonMaker::gridPoint::z
int z
Definition: volume.h:67
wustl_mm::SkeletonMaker::gridPoint::x
int x
Definition: volume.h:67
wustl_mm::SkeletonMaker::gridQueueEle
Definition: grid_queue.h:15
wustl_mm::SkeletonMaker::gridQueueEle::y
int y
Definition: grid_queue.h:16
wustl_mm::SkeletonMaker::gridQueueEle::x
int x
Definition: grid_queue.h:16
wustl_mm::SkeletonMaker::gridQueueEle::z
int z
Definition: grid_queue.h:16
MAX_ERODE
#define MAX_ERODE
Definition: volume.h:20
MAX_QUEUELEN
#define MAX_QUEUELEN
Definition: volume.h:19

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

◆ curveSkeleton() [2/2]

void Volume::curveSkeleton ( Volume *  grayvol,
float  lowthr,
float  highthr,
Volume *  svol 
)

insert them back into priority queue

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

Definition at line 3812 of file volume.cpp.

3813 {
3814 int i, j, k ;
3815 // First, threshold the volume
3816 #ifdef VERBOSE
3817 printf("Thresholding the volume to -1/0...\n") ;
3818 #endif
3819 threshold( 0.5f, -1, 0 ) ;
3820
3821 // Next, apply convergent erosion
3822 // by preserving: complex nodes, curve end-points, and sheet points
3823
3824 // Next, initialize the linked queue
3825 #ifdef VERBOSE
3826 printf("Initializing queue...\n") ;
3827 #endif
3828 GridQueue2* queue2 = new GridQueue2( ) ;
3829 GridQueue2* queue3 = new GridQueue2( ) ;
3830 GridQueue2* queue4 = new GridQueue2( ) ;
3831 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
3832
3833 for ( i = 0 ; i < getSizeX() ; i ++ )
3834 for ( j = 0 ; j < getSizeY() ; j ++ )
3835 for ( k = 0 ; k < getSizeZ() ; k ++ )
3836 {
3837 if ( getDataAt( i, j, k ) >= 0 )
3838 {
3839 float v = (float)grayvol->getDataAt(i,j,k) ;
3840 if ( v <= lowthr || v > highthr || svol->getDataAt(i,j,k) > 0 )
3841 {
3842 setDataAt( i, j, k, MAX_ERODE ) ;
3843 }
3844 else
3845 {
3846 for ( int m = 0 ; m < 6 ; m ++ )
3847 {
3848 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
3849 {
3850 // setDataAt( i, j, k, 1 ) ;
3851 queue2->prepend( i, j, k ) ;
3852 break ;
3853 }
3854 }
3855 }
3856 }
3857 }
3858 int wid = MAX_ERODE ;
3859 #ifdef VERBOSE
3860 printf("Total %d nodes\n", queue2->getNumElements() ) ;
3861 printf("Start erosion to %d...\n", wid) ;
3862 #endif
3863
3864
3865 // Perform erosion
3866 gridQueueEle* ele ;
3867 gridPoint* gp ;
3868 int ox, oy, oz ;
3869 int score ;
3870 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
3871 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
3872 {
3873 scrvol->setDataAt( i, -1 ) ;
3874 }
3875
3876 #ifdef NOISE_DIS_HELIX
3877 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
3878 #endif
3879
3880 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
3881 {
3882 // At the start of each iteration,
3883 // queue2 holds all the nodes for this layer
3884 // queue3 and queue are empty
3885
3886 int numComplex = 0, numSimple = 0 ;
3887 #ifdef VERBOSE
3888 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
3889 #endif
3890
3891 /*
3892 We first need to assign curwid + 1 to every node in this layer
3893 */
3894 queue2->reset() ;
3895 ele = queue2->getNext() ;
3896 while ( ele != NULL )
3897 {
3898 ox = ele->x ;
3899 oy = ele->y ;
3900 oz = ele->z ;
3901
3902 if ( getDataAt(ox,oy,oz) == curwid )
3903 {
3904 ele = queue2->remove() ;
3905 }
3906 else
3907 {
3908 setDataAt(ox,oy,oz, curwid) ;
3909 ele = queue2->getNext() ;
3910 }
3911 }
3912 queue4->reset() ;
3913 ele = queue4->getNext() ;
3914 while ( ele != NULL )
3915 {
3916 ox = ele->x ;
3917 oy = ele->y ;
3918 oz = ele->z ;
3919
3920 queue2->prepend(ox,oy,oz) ;
3921 ele = queue4->remove() ;
3922 }
3923
3924 // Now queue2 holds all the nodes for this layer
3925
3926 #ifdef NOISE_DIS_HELIX
3927 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
3928 queue2->reset() ;
3929
3930 // First run
3931 int flag = 0 ;
3932 while ( ( ele = queue2->getNext() ) != NULL )
3933 {
3934 ox = ele->x ;
3935 oy = ele->y ;
3936 oz = ele->z ;
3937 if ( NOISE_DIS_HELIX <= 1 )
3938 {
3939 noisevol->setDataAt( ox, oy, oz, 0 ) ;
3940 }
3941 else
3942 {
3943 flag = 0 ;
3944 for ( int m = 0 ; m < 6 ; m ++ )
3945 {
3946 int nx = ox + neighbor6[m][0] ;
3947 int ny = oy + neighbor6[m][1] ;
3948 int nz = oz + neighbor6[m][2] ;
3949 if ( getDataAt( nx, ny, nz ) == 0 )
3950 {
3951 noisevol->setDataAt( ox, oy, oz, 1 ) ;
3952 flag = 1 ;
3953 break ;
3954 }
3955 }
3956 if ( ! flag )
3957 {
3958 noisevol->setDataAt( ox, oy, oz, 0 ) ;
3959 }
3960 }
3961 }
3962
3963 int cur, visited ;
3964 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
3965 {
3966 queue2->reset() ;
3967 int count = 0 ;
3968 visited = 0 ;
3969
3970 while ( ( ele = queue2->getNext() ) != NULL )
3971 {
3972 ox = ele->x ;
3973 oy = ele->y ;
3974 oz = ele->z ;
3975
3976 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
3977 {
3978 visited ++ ;
3979 continue ;
3980 }
3981
3982 flag = 0 ;
3983 for ( int m = 0 ; m < 6 ; m ++ )
3984 {
3985 int nx = ox + neighbor6[m][0] ;
3986 int ny = oy + neighbor6[m][1] ;
3987 int nz = oz + neighbor6[m][2] ;
3988 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
3989 {
3990 noisevol->setDataAt( ox, oy, oz, 1 ) ;
3991 visited ++ ;
3992 count ++ ;
3993 break ;
3994 }
3995 }
3996 }
3997
3998 if ( count == 0 )
3999 {
4000 break ;
4001 }
4002 }
4003 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
4004
4005
4006 #endif
4007 /* Commented out for debugging
4008
4009 // First,
4010 // check for complex nodes in queue2
4011 // move them from queue2 to queue3
4012 queue2->reset() ;
4013 ele = queue2->getNext() ;
4014 while ( ele != NULL )
4015 {
4016 ox = ele->x ;
4017 oy = ele->y ;
4018 oz = ele->z ;
4019
4020 // Check simple
4021 #ifndef NOISE_DIS_HELIX
4022 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4023 #else
4024 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
4025 #endif
4026 {
4027 // Complex, set to next layer
4028 setDataAt( ox, oy, oz, curwid + 1 ) ;
4029 queue3->prepend( ox, oy, oz ) ;
4030 ele = queue2->remove() ;
4031
4032 numComplex ++ ;
4033 }
4034 else
4035 {
4036 ele = queue2->getNext() ;
4037 }
4038 }
4039 */
4040
4041 // Next,
4042 // Compute score for each node left in queue2
4043 // move them into priority queue
4044 queue2->reset() ;
4045 ele = queue2->getNext() ;
4046 while ( ele != NULL )
4047 {
4048 ox = ele->x ;
4049 oy = ele->y ;
4050 oz = ele->z ;
4051
4052 // Compute score
4053 score = getNumPotComplex2( ox, oy, oz ) ;
4054 scrvol->setDataAt( ox, oy, oz, score ) ;
4055
4056 // Push to queue
4057 gp = new gridPoint ;
4058 gp->x = ox ;
4059 gp->y = oy ;
4060 gp->z = oz ;
4061 // queue->add( gp, -score ) ;
4062 queue->add( gp, score ) ;
4063
4064 ele = queue2->remove() ;
4065 }
4066
4067 // Rename queue3 to be queue2,
4068 // Clear queue3
4069 // From now on, queue2 holds nodes of next level
4070 delete queue2 ;
4071 queue2 = queue3 ;
4072 queue3 = new GridQueue2( ) ;
4073
4074 // Next, start priority queue iteration
4075 while ( ! queue->isEmpty() )
4076 {
4077 // Retrieve the node with the highest score
4078 queue->remove( gp, score ) ;
4079 ox = gp->x ;
4080 oy = gp->y ;
4081 oz = gp->z ;
4082 delete gp ;
4083 // score = -score ;
4084
4085 // Ignore the node
4086 // if it has been processed before
4087 // or it has an updated score
4088 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
4089 {
4090 continue ;
4091 }
4092
4093 /* Commented out for debugging
4094
4095 // Remove this simple node
4096 setDataAt( ox, oy, oz, -1 ) ;
4097 numSimple ++ ;
4098 // printf("Highest score: %d\n", score) ;
4099 */
4100
4101 /* Added for debugging */
4102 // Check simple
4103 #ifndef NOISE_DIS_HELIX
4104 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
4105 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4106 #else
4107 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4108 #endif
4109 {
4110 // Complex, set to next layer
4111 setDataAt( ox, oy, oz, curwid + 1 ) ;
4112 queue4->prepend( ox, oy, oz ) ;
4113 numComplex ++ ;
4114 }
4115 else
4116 {
4117 setDataAt( ox, oy, oz, -1 ) ;
4118 numSimple ++ ;
4119
4120
4121 /* Adding ends */
4122 // Move its neighboring unvisited node to queue2
4123 for ( int m = 0 ; m < 6 ; m ++ )
4124 {
4125 int nx = ox + neighbor6[m][0] ;
4126 int ny = oy + neighbor6[m][1] ;
4127 int nz = oz + neighbor6[m][2] ;
4128 if ( getDataAt( nx, ny, nz ) == 0 )
4129 {
4130 // setDataAt( nx, ny, nz, curwid + 1 ) ;
4131 queue2->prepend( nx, ny, nz ) ;
4132 }
4133 }
4134
4135 }
4136
4137
4138 /* Commented out for debugging
4139
4140 // Find complex nodes in its 3x3 neighborhood
4141 // move them to queue2
4142 for ( i = -1 ; i < 2 ; i ++ )
4143 for ( j = -1 ; j < 2 ; j ++ )
4144 for ( k = -1 ; k < 2 ; k ++ )
4145 {
4146 int nx = ox + i ;
4147 int ny = oy + j ;
4148 int nz = oz + k ;
4149
4150 // Check simple
4151 if ( getDataAt( nx, ny, nz ) == curwid &&
4152 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
4153 #ifndef NOISE_DIS_HELIX
4154 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
4155 #else
4156 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
4157 #endif
4158
4159 {
4160 // Complex, set to next layer
4161 setDataAt( nx, ny, nz, curwid + 1 ) ;
4162 queue2->prepend( nx, ny, nz ) ;
4163 numComplex ++ ;
4164 }
4165 }
4166 */
4167
4168 // Update scores for nodes in its 5x5 neighborhood
4169 // insert them back into priority queue
4170
4171 for ( i = -2 ; i < 3 ;i ++ )
4172 for ( j = -2 ; j < 3 ; j ++ )
4173 for ( k = -2 ; k < 3 ; k ++ )
4174 {
4175 int nx = ox + i ;
4176 int ny = oy + j ;
4177 int nz = oz + k ;
4178
4179 if ( getDataAt( nx, ny, nz ) == curwid )
4180 {
4181 // Compute score
4182 score = getNumPotComplex2( nx, ny, nz ) ;
4183
4184 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
4185 {
4186 // printf("Update\n") ;
4187 scrvol->setDataAt( nx, ny, nz, score ) ;
4188 // Push to queue
4189 gp = new gridPoint ;
4190 gp->x = nx ;
4191 gp->y = ny ;
4192 gp->z = nz ;
4193 // queue->add( gp, -score ) ;
4194 queue->add( gp, score ) ;
4195 }
4196 }
4197 }
4198
4199
4200 }
4201
4202 #ifdef VERBOSE
4203 printf("%d complex, %d simple\n", numComplex, numSimple) ;
4204 #endif
4205
4206 if ( numSimple == 0 )
4207 {
4208 if ( queue2->getNumElements() > 0 )
4209 {
4210 printf("*************************wierd here*************************\n");
4211 }
4212 break ;
4213 }
4214 }
4215
4216 // Remove all internal voxels (contained in manifold surfaces)
4217 queue2->reset() ;
4218 queue4->reset() ;
4219 ele = queue4->getNext() ;
4220 while ( ele != NULL )
4221 {
4222 ox = ele->x ;
4223 oy = ele->y ;
4224 oz = ele->z ;
4225
4226 if ( isPiercable(ox,oy,oz) == 1 ) // hasCompleteSheet( ox, oy, oz ) == 1 ) //
4227 {
4228 queue2->prepend(ox,oy,oz) ;
4229 // setDataAt( ox, oy, oz, -1 ) ;
4230 }
4231 ele = queue4->remove() ;
4232 }
4233
4234 for ( i = 0 ; i < getSizeX() ; i ++ )
4235 for ( j = 0 ; j < getSizeY() ; j ++ )
4236 for ( k = 0 ; k < getSizeZ() ; k ++ )
4237 {
4238 if ( getDataAt( i, j, k ) == 0 && isPiercable(i,j,k) ) //hasCompleteSheet(i,j,k) == 1) //
4239 {
4240 queue2->prepend( i, j, k ) ;
4241 }
4242 }
4243 queue2->reset() ;
4244 ele = queue2->getNext() ;
4245 while ( ele != NULL )
4246 {
4247 ox = ele->x ;
4248 oy = ele->y ;
4249 oz = ele->z ;
4250 setDataAt( ox, oy, oz, -1 ) ;
4251 ele = queue2->remove() ;
4252 }
4253
4254
4255 // Finally, clean up
4256 delete scrvol;
4257 delete queue;
4258 delete queue2;
4259 delete queue3;
4260 delete queue4;
4261 #ifdef VERBOSE
4262 printf("Thresholding the volume to 0/1...\n") ;
4263 #endif
4264 threshold( 0, 0, 1 ) ;
4265 }
wustl_mm::SkeletonMaker::Volume::isPiercable
int isPiercable(int ox, int oy, int oz)
Definition: volume.cpp:2381

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

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

◆ curveSkeleton2D()

void Volume::curveSkeleton2D ( float  thr,
Volume *  svol 
)

Definition at line 4678 of file volume.cpp.

4679 {
4680 int i, j, k ;
4681 // First, threshold the volume
4682 #ifdef VERBOSE
4683 printf("Thresholding the volume to -1/0...\n") ;
4684 #endif
4685 threshold( thr, -1, 0 ) ;
4686
4687 // Next, apply convergent erosion
4688 // by preserving: complex nodes, curve end-points, and sheet points
4689
4690 // Next, initialize the linked queue
4691 #ifdef VERBOSE
4692 printf("Initializing queue...\n") ;
4693 #endif
4694 GridQueue2* queue2 = new GridQueue2( ) ;
4695 GridQueue2* queue3 = new GridQueue2( ) ;
4696 GridQueue2* queue4 = new GridQueue2( ) ;
4697 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
4698
4699 for ( i = 0 ; i < getSizeX() ; i ++ )
4700 for ( j = 0 ; j < getSizeY() ; j ++ )
4701 for ( k = 0 ; k < getSizeZ() ; k ++ )
4702 {
4703 if ( getDataAt( i, j, k ) >= 0 )
4704 {
4705 if ( svol->getDataAt(i,j,k) > 0 )
4706 {
4707 setDataAt( i, j, k, MAX_ERODE ) ;
4708 }
4709 else
4710 {
4711 for ( int m = 0 ; m < 4 ; m ++ )
4712 {
4713 if ( getDataAt( i + neighbor4[m][0], j + neighbor4[m][1], k ) < 0 )
4714 {
4715 // setDataAt( i, j, k, 1 ) ;
4716 queue2->prepend( i, j, k ) ;
4717 break ;
4718 }
4719 }
4720 }
4721 }
4722 }
4723 int wid = MAX_ERODE ;
4724 #ifdef VERBOSE
4725 printf("Total %d nodes\n", queue2->getNumElements() ) ;
4726 printf("Start erosion to %d...\n", wid) ;
4727 #endif
4728
4729
4730 // Perform erosion
4731 gridQueueEle* ele ;
4732 gridPoint* gp ;
4733 int ox, oy, oz ;
4734 int score ;
4735 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
4736 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
4737 {
4738 scrvol->setDataAt( i, -1 ) ;
4739 }
4740
4741 #ifdef NOISE_DIS_HELIX
4742 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
4743 #endif
4744
4745 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
4746 {
4747 // At the start of each iteration,
4748 // queue2 holds all the nodes for this layer
4749 // queue3 and queue are empty
4750
4751 int numComplex = 0, numSimple = 0 ;
4752 #ifdef VERBOSE
4753 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
4754 #endif
4755
4756 /*
4757 We first need to assign curwid + 1 to every node in this layer
4758 */
4759 queue2->reset() ;
4760 ele = queue2->getNext() ;
4761 while ( ele != NULL )
4762 {
4763 ox = ele->x ;
4764 oy = ele->y ;
4765 oz = ele->z ;
4766
4767 if ( getDataAt(ox,oy,oz) == curwid )
4768 {
4769 ele = queue2->remove() ;
4770 }
4771 else
4772 {
4773 setDataAt(ox,oy,oz, curwid) ;
4774 ele = queue2->getNext() ;
4775 }
4776 }
4777 queue4->reset() ;
4778 ele = queue4->getNext() ;
4779 while ( ele != NULL )
4780 {
4781 ox = ele->x ;
4782 oy = ele->y ;
4783 oz = ele->z ;
4784
4785 queue2->prepend(ox,oy,oz) ;
4786 ele = queue4->remove() ;
4787 }
4788
4789 // Now queue2 holds all the nodes for this layer
4790
4791 #ifdef NOISE_DIS_HELIX
4792 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
4793 queue2->reset() ;
4794
4795 // First run
4796 int flag = 0 ;
4797 while ( ( ele = queue2->getNext() ) != NULL )
4798 {
4799 ox = ele->x ;
4800 oy = ele->y ;
4801 oz = ele->z ;
4802 if ( NOISE_DIS_HELIX <= 1 )
4803 {
4804 noisevol->setDataAt( ox, oy, oz, 0 ) ;
4805 }
4806 else
4807 {
4808 flag = 0 ;
4809 for ( int m = 0 ; m < 6 ; m ++ )
4810 {
4811 int nx = ox + neighbor6[m][0] ;
4812 int ny = oy + neighbor6[m][1] ;
4813 int nz = oz + neighbor6[m][2] ;
4814 if ( getDataAt( nx, ny, nz ) == 0 )
4815 {
4816 noisevol->setDataAt( ox, oy, oz, 1 ) ;
4817 flag = 1 ;
4818 break ;
4819 }
4820 }
4821 if ( ! flag )
4822 {
4823 noisevol->setDataAt( ox, oy, oz, 0 ) ;
4824 }
4825 }
4826 }
4827
4828 int cur, visited ;
4829 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
4830 {
4831 queue2->reset() ;
4832 int count = 0 ;
4833 visited = 0 ;
4834
4835 while ( ( ele = queue2->getNext() ) != NULL )
4836 {
4837 ox = ele->x ;
4838 oy = ele->y ;
4839 oz = ele->z ;
4840
4841 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
4842 {
4843 visited ++ ;
4844 continue ;
4845 }
4846
4847 flag = 0 ;
4848 for ( int m = 0 ; m < 6 ; m ++ )
4849 {
4850 int nx = ox + neighbor6[m][0] ;
4851 int ny = oy + neighbor6[m][1] ;
4852 int nz = oz + neighbor6[m][2] ;
4853 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
4854 {
4855 noisevol->setDataAt( ox, oy, oz, 1 ) ;
4856 visited ++ ;
4857 count ++ ;
4858 break ;
4859 }
4860 }
4861 }
4862
4863 if ( count == 0 )
4864 {
4865 break ;
4866 }
4867 }
4868 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
4869
4870
4871 #endif
4872 /* Commented out for debugging
4873
4874 // First,
4875 // check for complex nodes in queue2
4876 // move them from queue2 to queue3
4877 queue2->reset() ;
4878 ele = queue2->getNext() ;
4879 while ( ele != NULL )
4880 {
4881 ox = ele->x ;
4882 oy = ele->y ;
4883 oz = ele->z ;
4884
4885 // Check simple
4886 #ifndef NOISE_DIS_HELIX
4887 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4888 #else
4889 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
4890 #endif
4891 {
4892 // Complex, set to next layer
4893 setDataAt( ox, oy, oz, curwid + 1 ) ;
4894 queue3->prepend( ox, oy, oz ) ;
4895 ele = queue2->remove() ;
4896
4897 numComplex ++ ;
4898 }
4899 else
4900 {
4901 ele = queue2->getNext() ;
4902 }
4903 }
4904 */
4905
4906 // Next,
4907 // Compute score for each node left in queue2
4908 // move them into priority queue
4909 queue2->reset() ;
4910 ele = queue2->getNext() ;
4911 while ( ele != NULL )
4912 {
4913 ox = ele->x ;
4914 oy = ele->y ;
4915 oz = ele->z ;
4916
4917 // Compute score
4918 score = getNumPotComplex2( ox, oy, oz ) ;
4919 //score = getNumNeighbor6( ox, oy, oz ) ;
4920 scrvol->setDataAt( ox, oy, oz, score ) ;
4921
4922 // Push to queue
4923 gp = new gridPoint ;
4924 gp->x = ox ;
4925 gp->y = oy ;
4926 gp->z = oz ;
4927 // queue->add( gp, -score ) ;
4928 queue->add( gp, score ) ;
4929
4930 ele = queue2->remove() ;
4931 }
4932
4933 // Rename queue3 to be queue2,
4934 // Clear queue3
4935 // From now on, queue2 holds nodes of next level
4936 delete queue2 ;
4937 queue2 = queue3 ;
4938 queue3 = new GridQueue2( ) ;
4939
4940 // Next, start priority queue iteration
4941 while ( ! queue->isEmpty() )
4942 {
4943 // Retrieve the node with the highest score
4944 queue->remove( gp, score ) ;
4945 ox = gp->x ;
4946 oy = gp->y ;
4947 oz = gp->z ;
4948 delete gp ;
4949 // score = -score ;
4950
4951 // Ignore the node
4952 // if it has been processed before
4953 // or it has an updated score
4954 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
4955 {
4956 continue ;
4957 }
4958
4959 /* Commented out for debugging
4960
4961 // Remove this simple node
4962 setDataAt( ox, oy, oz, -1 ) ;
4963 numSimple ++ ;
4964 // printf("Highest score: %d\n", score) ;
4965 */
4966
4967 /* Added for debugging */
4968 // Check simple
4969 #ifndef NOISE_DIS_HELIX
4970 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
4971 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4972 #else
4973 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4974 #endif
4975 {
4976 // Complex, set to next layer
4977 setDataAt( ox, oy, oz, curwid + 1 ) ;
4978 queue4->prepend( ox, oy, oz ) ;
4979 numComplex ++ ;
4980 }
4981 else
4982 {
4983 setDataAt( ox, oy, oz, -1 ) ;
4984 numSimple ++ ;
4985 }
4986 /* Adding ends */
4987
4988 // Move its neighboring unvisited node to queue2
4989 for ( int m = 0 ; m < 4 ; m ++ )
4990 {
4991 int nx = ox + neighbor4[m][0] ;
4992 int ny = oy + neighbor4[m][1] ;
4993 int nz = oz ;
4994 if ( getDataAt( nx, ny, nz ) == 0 )
4995 {
4996 // setDataAt( nx, ny, nz, curwid + 1 ) ;
4997 queue2->prepend( nx, ny, nz ) ;
4998 }
4999 }
5000
5001 /* Commented out for debugging
5002
5003 // Find complex nodes in its 3x3 neighborhood
5004 // move them to queue2
5005 for ( i = -1 ; i < 2 ; i ++ )
5006 for ( j = -1 ; j < 2 ; j ++ )
5007 for ( k = -1 ; k < 2 ; k ++ )
5008 {
5009 int nx = ox + i ;
5010 int ny = oy + j ;
5011 int nz = oz + k ;
5012
5013 // Check simple
5014 if ( getDataAt( nx, ny, nz ) == curwid &&
5015 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
5016 #ifndef NOISE_DIS_HELIX
5017 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
5018 #else
5019 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
5020 #endif
5021
5022 {
5023 // Complex, set to next layer
5024 setDataAt( nx, ny, nz, curwid + 1 ) ;
5025 queue2->prepend( nx, ny, nz ) ;
5026 numComplex ++ ;
5027 }
5028 }
5029 */
5030
5031 // Update scores for nodes in its 5x5 neighborhood
5032 // insert them back into priority queue
5033
5034 for ( i = -2 ; i < 3 ;i ++ )
5035 for ( j = -2 ; j < 3 ; j ++ )
5036 for ( k = -2 ; k < 3 ; k ++ )
5037 {
5038 int nx = ox + i ;
5039 int ny = oy + j ;
5040 int nz = oz + k ;
5041
5042 if ( getDataAt( nx, ny, nz ) == curwid )
5043 {
5044 // Compute score
5045 score = getNumPotComplex2( nx, ny, nz ) ;
5046 //score = getNumNeighbor6( nx, ny, nz ) ;
5047
5048 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
5049 {
5050 // printf("Update\n") ;
5051 scrvol->setDataAt( nx, ny, nz, score ) ;
5052 // Push to queue
5053 gp = new gridPoint ;
5054 gp->x = nx ;
5055 gp->y = ny ;
5056 gp->z = nz ;
5057 // queue->add( gp, -score ) ;
5058 queue->add( gp, score ) ;
5059 }
5060 }
5061 }
5062
5063
5064 }
5065
5066 #ifdef VERBOSE
5067 printf("%d complex, %d simple\n", numComplex, numSimple) ;
5068 #endif
5069
5070 if ( numSimple == 0 )
5071 {
5072 break ;
5073 }
5074 }
5075
5076 // Finally, clean up
5077 #ifdef VERBOSE
5078 printf("Thresholding the volume to 0/1...\n") ;
5079 #endif
5080 threshold( 0, 0, 1 ) ;
5081 delete scrvol;
5082 delete queue;
5083 delete queue2;
5084 delete queue3;
5085 delete queue4;
5086 }
wustl_mm::SkeletonMaker::neighbor4
const int neighbor4[4][2]
Definition: volume.h:28

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

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

◆ erodeHelix() [1/2]

void Volume::erodeHelix ( )

Definition at line 5905 of file volume.cpp.

5906 {
5907 erodeHelix( 3 ) ;
5908 }
wustl_mm::SkeletonMaker::Volume::erodeHelix
void erodeHelix()
Definition: volume.cpp:5905

References erodeHelix().

Referenced by erodeHelix(), and wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PruneCurves().

◆ erodeHelix() [2/2]

void Volume::erodeHelix ( int  disthr)

Definition at line 5911 of file volume.cpp.

5912 {
5913 int i, j, k ;
5914 // First, threshold the volume
5915 //printf("Thresholding the volume to -1/0...\n") ;
5916 threshold( 0.1f, -1, 0 ) ;
5917
5918 /* Debug: remove faces */
5919 //Volume* facevol = markFaceEdge() ;
5920 /* End debugging */
5921
5922 // Next, initialize the linked queue
5923 //printf("Initializing queue...\n") ;
5924 Volume * fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
5925 GridQueue2* queue2 = new GridQueue2( ) ;
5926 GridQueue2* queue3 = new GridQueue2( ) ;
5927 GridQueue2** queues = new GridQueue2* [ 10000 ] ;
5928
5929 for ( i = 0 ; i < getSizeX() ; i ++ )
5930 for ( j = 0 ; j < getSizeY() ; j ++ )
5931 for ( k = 0 ; k < getSizeZ() ; k ++ )
5932 {
5933 if ( getDataAt( i, j, k ) >= 0 )
5934 {
5935 if ( ! hasCompleteHelix( i, j, k ) )
5936 // if ( ! hasCompleteHelix( i, j, k, facevol ) )
5937 {
5938 queue2->prepend( i, j, k ) ;
5939 fvol->setDataAt( i, j, k, -1 ) ;
5940 }
5941 }
5942 }
5943 //printf("Total %d nodes\n", queue2->getNumElements() ) ;
5944
5945 // Start erosion
5946 gridQueueEle* ele ;
5947 int dis = -1 ;
5948 while ( queue2->getNumElements() > 0 )
5949 {
5950 // First, set distance
5951 dis -- ;
5952 queues[ - dis ] = new GridQueue2( ) ;
5953 //printf("Distance transform to %d...", dis) ;
5954 queue2->reset() ;
5955 while( ( ele = queue2->getNext() ) != NULL )
5956 {
5957 setDataAt( ele->x, ele->y, ele->z, dis ) ;
5958 queues[ -dis ]->prepend( ele->x, ele->y, ele->z ) ;
5959 }
5960 //printf("%d nodes\n", queues[-dis]->getNumElements()) ;
5961
5962 // Next, find next layer
5963 queue2->reset() ;
5964 ele = queue2->getNext() ;
5965 while ( ele != NULL )
5966 {
5967 for ( int m = 0 ; m < 6 ; m ++ )
5968 {
5969 int nx = ele->x + neighbor6[m][0] ;
5970 int ny = ele->y + neighbor6[m][1] ;
5971 int nz = ele->z + neighbor6[m][2] ;
5972 if ( getDataAt( nx, ny, nz ) == 0 )
5973 {
5974 fvol->setDataAt( nx, ny, nz, dis ) ;
5975
5976 if ( ! hasCompleteHelix( nx, ny, nz ) )
5977 // if ( ! hasCompleteHelix( nx, ny, nz, facevol ) )
5978 {
5979 setDataAt( nx, ny, nz, 1 ) ;
5980 queue3->prepend( nx, ny, nz ) ;
5981 }
5982 }
5983 }
5984
5985 ele = queue2->remove() ;
5986 }
5987
5988 // Next, swap queues
5989 GridQueue2* temp = queue2 ;
5990 queue2 = queue3 ;
5991 queue3 = temp ;
5992 }
5993
5994 // Deal with closed rings
5995 dis -- ;
5996 queues[ - dis ] = new GridQueue2( ) ;
5997 #ifdef VERBOSE
5998 printf("Detecting closed rings %d...", dis) ;
5999 #endif
6000 int ftot = 0 ;
6001 for ( i = 0 ; i < getSizeX() ; i ++ )
6002 for ( j = 0 ; j < getSizeY() ; j ++ )
6003 for ( k = 0 ; k < getSizeZ() ; k ++ )
6004 {
6005 if ( getDataAt( i, j, k ) == 0 )
6006 {
6007 setDataAt( i, j, k, dis ) ;
6008 queues[ -dis ]->prepend( i, j, k ) ;
6009 /*
6010 int fval = (int) fvol->getDataAt( i, j, k ) ;
6011 if ( fval == 0)
6012 {
6013 // queues[ -dis ]->prepend( i, j, k ) ;
6014 }
6015 else
6016 {
6017 setDataAt( i, j, k, fval - 1 ) ;
6018 // queues[ -fval + 1 ]->prepend( i, j, k ) ;
6019 }
6020 */
6021 ftot ++ ;
6022 }
6023 }
6024 #ifdef VERBOSE
6025 printf("%d nodes\n", ftot) ;
6026 #endif
6027
6028
6029 // return ;
6030
6031 /* Find local minimum: to help determine erosion level
6032 int cts[ 64 ] ;
6033 for ( i = 0 ; i <= - dis ; i ++ )
6034 {
6035 cts[ i ] = 0 ;
6036 }
6037 for ( i = 0 ; i < getSizeX() ; i ++ )
6038 for ( j = 0 ; j < getSizeY() ; j ++ )
6039 for ( k = 0 ; k < getSizeZ() ; k ++ )
6040 {
6041 double val = getDataAt( i, j, k ) ;
6042 if ( val < -1 )
6043 {
6044 int min = 1 ;
6045 for ( int m = 0 ; m < 6 ; m ++ )
6046 {
6047 int nx = i + neighbor6[m][0] ;
6048 int ny = j + neighbor6[m][1] ;
6049 int nz = k + neighbor6[m][2] ;
6050 if ( getDataAt( nx, ny, nz ) < val )
6051 {
6052 min = 0 ;
6053 break ;
6054 }
6055 }
6056
6057 if ( min )
6058 {
6059 cts[ (int) - val ] ++ ;
6060 }
6061 }
6062 }
6063
6064 for ( i = 2 ; i <= - dis ; i ++ )
6065 {
6066 printf("Local minima: %d with %d nodes.\n", -i, cts[ i ] ) ;
6067 }
6068 */
6069
6070 // Dilate back
6071 // Starting from nodes with distance - 2 - disthr
6072
6073 if ( disthr + 2 > - dis )
6074 {
6075 disthr = - dis - 2 ;
6076 }
6077 int d;
6078
6079 for ( d = - dis ; d > disthr + 1 ; d -- )
6080 {
6081 queues[ d ]->reset() ;
6082 while ( (ele = queues[ d ]->getNext() ) != NULL )
6083 {
6084 setDataAt( ele->x, ele->y, ele->z, d ) ;
6085 }
6086 }
6087
6088
6089 for ( d = disthr + 1 ; d >= 2 ; d -- )
6090 {
6091 //delete queue3 ;
6092 //queue3 = new GridQueue2( ) ;
6093 queues[ d ]->reset() ;
6094 ele = queues[ d ]->getNext() ;
6095 while ( ele != NULL )
6096 {
6097 int dilatable = 0 ;
6098 for ( int m = 0 ; m < 6 ; m ++ )
6099 {
6100 int nx = ele->x + neighbor6[m][0] ;
6101 int ny = ele->y + neighbor6[m][1] ;
6102 int nz = ele->z + neighbor6[m][2] ;
6103 if ( getDataAt( nx, ny, nz ) == d + 1 )
6104 {
6105 dilatable = 1 ;
6106 break ;
6107 }
6108 }
6109
6110
6111 if ( ! dilatable )
6112 {
6113 /*
6114 setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
6115 */
6116
6117 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
6118 if ( d > 2 )
6119 {
6120 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
6121 }
6122 ele = queues[ d ]->remove() ;
6123 }
6124 else
6125 {
6126 setDataAt( ele->x, ele->y, ele->z, d ) ;
6127 ele = queues[ d ]->getNext() ;
6128 }
6129
6130 }
6131
6132 /* Debug: extract minimal set */
6133 while ( 1 )
6134 {
6135 int num = 0 ;
6136 queues[ d ]->reset() ;
6137 ele = queues[ d ]->getNext() ;
6138 while ( ele != NULL )
6139 {
6140 int critical = 0 ;
6141 setDataAt( ele->x, ele->y, ele->z, -1 ) ;
6142
6143 for ( i = -1 ; i < 2 ; i ++ )
6144 {
6145 for ( j = -1 ; j < 2 ; j ++ )
6146 {
6147 for ( k = -1 ; k < 2 ; k ++ )
6148 {
6149 if ( i != 0 && j != 0 && k != 0 )
6150 {
6151 continue ;
6152 }
6153 int nx = ele->x + i ;
6154 int ny = ele->y + j ;
6155 int nz = ele->z + k ;
6156 if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteHelix( nx,ny,nz ) ) //, facevol ) )
6157 {
6158 critical = 1 ;
6159 break ;
6160 }
6161 }
6162 if ( critical )
6163 {
6164 break ;
6165 }
6166 }
6167 if ( critical )
6168 {
6169 break ;
6170 }
6171 }
6172
6173 if ( critical )
6174 {
6175 setDataAt( ele->x, ele->y, ele->z, d ) ;
6176 ele = queues[ d ]->getNext() ;
6177 }
6178 else
6179 {
6180 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
6181 if ( d > 2 )
6182 {
6183 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
6184 }
6185 ele = queues[ d ]->remove() ;
6186 num ++ ;
6187 }
6188
6189 }
6190
6191 #ifdef VERBOSE
6192 printf("Non-minimal: %d\n", num) ;
6193 #endif
6194
6195 if ( num == 0 )
6196 {
6197 break ;
6198 }
6199 }
6200
6201
6202 /* End of debugging */
6203
6204 /*
6205 queue3->reset() ;
6206 ele = queue3->getNext() ;
6207 while ( ele != NULL )
6208 {
6209 setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
6210 ele = queue3->remove() ;
6211 }
6212 */
6213 }
6214
6215 // Finally, threshold the volume
6216 #ifdef VERBOSE
6217 //printf("Thresholding the volume to 0/1...\n") ;
6218 #endif
6219 //threshold( -1, 1, 0, 0 ) ;
6220 threshold( 0, 0, 1 ) ;
6221 delete fvol ;
6222 delete queue2;
6223 delete queue3;
6224 for ( d = -dis ; d >= 2 ; d -- ) {
6225 delete queues[d];
6226 }
6227 delete [] queues;
6228
6229 }
wustl_mm::SkeletonMaker::Volume::hasCompleteHelix
int hasCompleteHelix(int ox, int oy, int oz)
Definition: volume.cpp:1456

References getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getSizeX(), getSizeY(), getSizeZ(), hasCompleteHelix(), wustl_mm::SkeletonMaker::neighbor6, 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.

◆ erodeSheet() [1/2]

int Volume::erodeSheet ( )

Definition at line 6234 of file volume.cpp.

6235 {
6236 return erodeSheet( 3 ) ;
6237 }
wustl_mm::SkeletonMaker::Volume::erodeSheet
int erodeSheet()
Definition: volume.cpp:6234

References erodeSheet().

Referenced by erodeSheet(), and wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PruneSurfaces().

◆ erodeSheet() [2/2]

int Volume::erodeSheet ( int  disthr)

Definition at line 6240 of file volume.cpp.

6241 {
6242 int i, j, k ;
6243 // First, threshold the volume
6244 //printf("Thresholding the volume to -1/0...\n") ;
6245 threshold( 0.1f, -1, 0 ) ;
6246
6247 /* Debug: remove cells */
6248 Volume* facevol = markCellFace() ;
6249 /* End debugging */
6250
6251 // Next, initialize the linked queue
6252 //printf("Initializing queue...\n") ;
6253 Volume * fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
6254 GridQueue2* queue2 = new GridQueue2( ) ;
6255 GridQueue2* queue3 = new GridQueue2( ) ;
6256 GridQueue2** queues = new GridQueue2* [ 10000 ] ;
6257
6258 for ( i = 0 ; i < getSizeX() ; i ++ )
6259 for ( j = 0 ; j < getSizeY() ; j ++ )
6260 for ( k = 0 ; k < getSizeZ() ; k ++ )
6261 {
6262 if ( getDataAt( i, j, k ) >= 0 )
6263 {
6264 if ( ! hasCompleteSheet( i, j, k ) )
6265 //if ( ! hasCompleteSheet( i, j, k, facevol ) )
6266 {
6267 queue2->prepend( i, j, k ) ;
6268 fvol->setDataAt( i, j, k, -1 ) ;
6269 }
6270 }
6271 }
6272 #ifdef VERBOSE
6273 printf("Total %d nodes\n", queue2->getNumElements() ) ;
6274 #endif
6275
6276 // Start erosion
6277 gridQueueEle* ele ;
6278 int dis = -1 ;
6279 while ( queue2->getNumElements() > 0 )
6280 {
6281 // First, set distance
6282 dis -- ;
6283 queues[ - dis ] = new GridQueue2( ) ;
6284 //printf("Distance transform to %d...", dis) ;
6285 queue2->reset() ;
6286 while( ( ele = queue2->getNext() ) != NULL )
6287 {
6288 setDataAt( ele->x, ele->y, ele->z, dis ) ;
6289 queues[ -dis ]->prepend( ele->x, ele->y, ele->z ) ;
6290 }
6291 //printf("%d nodes\n", queues[-dis]->getNumElements()) ;
6292
6293 // Next, find next layer
6294 queue2->reset() ;
6295 ele = queue2->getNext() ;
6296 while ( ele != NULL )
6297 {
6298 // for ( int m = 0 ; m < 6 ; m ++ )
6299 for ( int mx = -1 ; mx < 2 ; mx ++ )
6300 for ( int my = -1 ; my < 2 ; my ++ )
6301 for ( int mz = -1 ; mz < 2 ; mz ++ )
6302 {
6303 if ( mx != 0 && my != 0 && mz != 0 )
6304 {
6305 continue ;
6306 }
6307 //int nx = ele->x + neighbor6[m][0] ;
6308 //int ny = ele->y + neighbor6[m][1] ;
6309 //int nz = ele->z + neighbor6[m][2] ;
6310 int nx = ele->x + mx ;
6311 int ny = ele->y + my ;
6312 int nz = ele->z + mz ;
6313
6314 if ( getDataAt( nx, ny, nz ) == 0 )
6315 {
6316 fvol->setDataAt( nx, ny, nz, dis ) ;
6317
6318 if ( ! hasCompleteSheet( nx, ny, nz ) )
6319 // if ( ! hasCompleteSheet( nx, ny, nz, facevol ) )
6320 {
6321 setDataAt( nx, ny, nz, 1 ) ;
6322 queue3->prepend( nx, ny, nz ) ;
6323 }
6324 }
6325 }
6326
6327 ele = queue2->remove() ;
6328 }
6329
6330 // Next, swap queues
6331 GridQueue2* temp = queue2 ;
6332 queue2 = queue3 ;
6333 queue3 = temp ;
6334 }
6335
6336 /* Deal with closed rings */
6337
6338 dis -- ;
6339 queues[ - dis ] = new GridQueue2( ) ;
6340 #ifdef VERBOSE
6341 printf("Detecting closed rings %d...", dis) ;
6342 #endif
6343 int ftot = 0 ;
6344 for ( i = 0 ; i < getSizeX() ; i ++ )
6345 for ( j = 0 ; j < getSizeY() ; j ++ )
6346 for ( k = 0 ; k < getSizeZ() ; k ++ )
6347 {
6348 if ( getDataAt( i, j, k ) == 0 )
6349 {
6350 /*
6351 int fval = (int) fvol->getDataAt( i, j, k ) ;
6352 if ( fval == 0)
6353 {
6354 setDataAt( i, j, k, dis - 2 ) ;
6355 // queues[ -dis ]->prepend( i, j, k ) ;
6356 }
6357 else
6358 {
6359 setDataAt( i, j, k, fval - 1 ) ;
6360 queues[ -fval + 1 ]->prepend( i, j, k ) ;
6361 }
6362 */
6363 setDataAt( i, j, k, dis ) ;
6364 queues[ -dis ]->prepend( i, j, k ) ;
6365
6366 ftot ++ ;
6367 }
6368 }
6369 #ifdef VERBOSE
6370 printf("%d nodes\n", ftot) ;
6371 #endif
6372
6373
6374 /* Find local minimum: to help determine erosion level
6375 int cts[ 64 ] ;
6376 for ( i = 0 ; i <= - dis ; i ++ )
6377 {
6378 cts[ i ] = 0 ;
6379 }
6380 for ( i = 0 ; i < getSizeX() ; i ++ )
6381 for ( j = 0 ; j < getSizeY() ; j ++ )
6382 for ( k = 0 ; k < getSizeZ() ; k ++ )
6383 {
6384 double val = getDataAt( i, j, k ) ;
6385 if ( val < -1 )
6386 {
6387 int min = 1 ;
6388 for ( int m = 0 ; m < 6 ; m ++ )
6389 {
6390 int nx = i + neighbor6[m][0] ;
6391 int ny = j + neighbor6[m][1] ;
6392 int nz = k + neighbor6[m][2] ;
6393 if ( getDataAt( nx, ny, nz ) < val )
6394 {
6395 min = 0 ;
6396 break ;
6397 }
6398 }
6399
6400 if ( min )
6401 {
6402 cts[ (int) - val ] ++ ;
6403 }
6404 }
6405 }
6406
6407 for ( i = 2 ; i <= - dis ; i ++ )
6408 {
6409 printf("Local minima: %d with %d nodes.\n", -i, cts[ i ] ) ;
6410 }
6411 */
6412
6413 // return ;
6414
6415 // Dilate back
6416 // Starting from nodes with distance - 2 - disthr
6417 int d ;
6418 if ( disthr + 2 > - dis )
6419 {
6420 disthr = - dis - 2 ;
6421
6422 }
6423 for ( d = - dis ; d > disthr + 1 ; d -- )
6424 {
6425 queues[ d ]->reset() ;
6426 while ( (ele = queues[ d ]->getNext() ) != NULL )
6427 {
6428 setDataAt( ele->x, ele->y, ele->z, d ) ;
6429 }
6430 }
6431
6432 for (d = disthr + 1 ; d >= 2 ; d -- )
6433 {
6434
6435 //delete queue3 ;
6436 //queue3 = new GridQueue2( ) ;
6437 queues[ d ]->reset() ;
6438 ele = queues[ d ]->getNext() ;
6439 while ( ele != NULL )
6440 {
6441 int dilatable = 0 ;
6442 // for ( int m = 0 ; m < 6 ; m ++ )
6443 /*
6444 for ( int mx = -1 ; mx < 2 ; mx ++ )
6445 for ( int my = -1 ; my < 2 ; my ++ )
6446 for ( int mz = -1 ; mz < 2 ; mz ++ )
6447 {
6448 if ( mx == 0 || my == 0 || mz == 0 )
6449 {
6450 int nx = ele->x + mx ; // neighbor6[m][0] ;
6451 int ny = ele->y + my ; // neighbor6[m][1] ;
6452 int nz = ele->z + mz ; // neighbor6[m][2] ;
6453 if ( getDataAt( nx, ny, nz ) == - d - 1 )
6454 {
6455 dilatable = 1 ;
6456 break ;
6457 }
6458 }
6459 }
6460 */
6461 for ( i = 0 ; i < 12 ; i ++ )
6462 {
6463 int flag = 1, flag2 = 0 ;
6464 for ( j = 0 ; j < 4 ; j ++ )
6465 {
6466 int nx = ele->x + sheetNeighbor[i][j][0] ;
6467 int ny = ele->y + sheetNeighbor[i][j][1] ;
6468 int nz = ele->z + sheetNeighbor[i][j][2] ;
6469
6470 double val = getDataAt( nx, ny, nz ) ;
6471
6472 if ( val > - d && val < 0)
6473 {
6474 flag = 0 ;
6475 break ;
6476 }
6477 if ( val == d + 1 )
6478 {
6479 flag2 ++ ;
6480 }
6481 }
6482
6483 if ( flag && flag2 )
6484 {
6485 dilatable = 1 ;
6486 break ;
6487 }
6488 }
6489
6490 if ( ! dilatable )
6491 {
6492 // setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
6493 // queue3->prepend( ele->x, ele->y, ele->z ) ;
6494
6495 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
6496 if ( d > 2 )
6497 {
6498 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
6499 }
6500 ele = queues[ d ]->remove() ;
6501 }
6502 else
6503 {
6504 setDataAt( ele->x, ele->y, ele->z, d ) ;
6505 ele = queues[ d ]->getNext() ;
6506 }
6507 }
6508
6509 /* Debug: extract minimal set */
6510 while ( 1 )
6511 {
6512 int num = 0 ;
6513 queues[ d ]->reset() ;
6514 ele = queues[ d ]->getNext() ;
6515 while ( ele != NULL )
6516 {
6517 int critical = 0 ;
6518 setDataAt( ele->x, ele->y, ele->z, -1 ) ;
6519
6520 for ( i = -1 ; i < 2 ; i ++ )
6521 {
6522 for ( j = -1 ; j < 2 ; j ++ )
6523 {
6524 for ( k = -1 ; k < 2 ; k ++ )
6525 {
6526 if ( i != 0 && j != 0 && k != 0 )
6527 {
6528 continue ;
6529 }
6530 int nx = ele->x + i ;
6531 int ny = ele->y + j ;
6532 int nz = ele->z + k ;
6533 // if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteSheet( nx,ny,nz, facevol ) )
6534 if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteSheet( nx,ny,nz ) )
6535 {
6536 critical = 1 ;
6537 break ;
6538 }
6539 }
6540 if ( critical )
6541 {
6542 break ;
6543 }
6544 }
6545 if ( critical )
6546 {
6547 break ;
6548 }
6549 }
6550
6551 if ( critical )
6552 {
6553 setDataAt( ele->x, ele->y, ele->z, d ) ;
6554 ele = queues[ d ]->getNext() ;
6555 }
6556 else
6557 {
6558 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
6559 if ( d > 2 )
6560 {
6561 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
6562 }
6563 ele = queues[ d ]->remove() ;
6564 num ++ ;
6565 }
6566
6567 }
6568 #ifdef VERBOSE
6569 printf("Non-minimal: %d\n", num) ;
6570 #endif
6571
6572 if ( num == 0 )
6573 {
6574 break ;
6575 }
6576 }
6577
6578
6579 /* End of debugging */
6580
6581 /*
6582 queue3->reset() ;
6583 ele = queue3->getNext() ;
6584 while ( ele != NULL )
6585 {
6586 setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
6587 ele = queue3->remove() ;
6588 }
6589 */
6590 }
6591
6592
6593 // Finally, threshold the volume
6594 #ifdef VERBOSE
6595 //printf("Thresholding the volume to 0/1...\n") ;
6596 #endif
6597 //threshold( -1, 1, 0, 0 ) ;
6598 threshold( 0, 0, 1 ) ;
6599
6600 delete facevol ;
6601 delete fvol ;
6602 delete queue2;
6603 delete queue3;
6604 for (d = -dis ; d >= 2 ; d -- ) {
6605 delete queues[d];
6606 }
6607 delete [] queues;
6608
6609 return - dis - 1 ;
6610 }
wustl_mm::SkeletonMaker::Volume::hasCompleteSheet
int hasCompleteSheet(int ox, int oy, int oz, Volume *fvol)
Definition: volume.cpp:1174
wustl_mm::SkeletonMaker::Volume::markCellFace
Volume * markCellFace()
Definition: volume.cpp:1027
wustl_mm::SkeletonMaker::sheetNeighbor
const int sheetNeighbor[12][4][3]
Definition: volume.h:37

References getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getSizeX(), getSizeY(), getSizeZ(), hasCompleteSheet(), markCellFace(), 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.

◆ get_emdata()

EMData * Volume::get_emdata ( )

Definition at line 31 of file volume.cpp.

32 {
33 return this->getVolumeData()->get_emdata();
34 }
wustl_mm::SkeletonMaker::VolumeData::get_emdata
EMData * get_emdata()
Definition: volume_data.cpp:166

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

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

◆ getDataAt() [1/2]

double Volume::getDataAt ( int  index)

Definition at line 65 of file volume.cpp.

65 {
66 return volData->GetDataAt(index);
67 }
wustl_mm::SkeletonMaker::VolumeData::GetDataAt
float GetDataAt(int x, int y, int z)
Definition: volume_data.cpp:149

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

◆ getDataAt() [2/2]

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

Definition at line 60 of file volume.cpp.

61 {
62 return volData->GetDataAt(x, y, z);
63 }

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

◆ getIndex()

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

Definition at line 48 of file volume.cpp.

48 {
49 return volData->GetIndex(x, y, z);
50 }
wustl_mm::SkeletonMaker::VolumeData::GetIndex
int GetIndex(int x, int y, int z)
Definition: volume_data.cpp:158

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

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

◆ getNumNeighbor6()

int Volume::getNumNeighbor6 ( int  ox,
int  oy,
int  oz 
)

Definition at line 661 of file volume.cpp.

661 {
662 int rvalue = 0 ;
663 for ( int i = 0 ; i < 6 ; i ++ ) {
664 int nx = ox + neighbor6[i][0] ;
665 int ny = oy + neighbor6[i][1] ;
666 int nz = oz + neighbor6[i][2] ;
667 if ( getDataAt( nx, ny, nz ) >= 0 ) {
668 rvalue ++ ;
669 }
670 }
671
672 return rvalue ;
673 }

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

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

◆ getNumPotComplex()

int Volume::getNumPotComplex ( int  ox,
int  oy,
int  oz 
)

Definition at line 2548 of file volume.cpp.

2549 {
2550 //return 0 ;
2551
2552 int i;
2553 double val = getDataAt( ox, oy, oz ) ;
2554 if ( val <= 0 )
2555 {
2556 // return 70 ;
2557 }
2558
2559 // return getNumNeighbor6( ox, oy, oz ) ;
2560
2561 // int v = ((getNumNeighbor6( ox, oy, oz ) & 255) << 24) ;
2562 //int v = 0 ;
2563
2564 int rvalue = 0, nx, ny, nz ;
2565 setDataAt( ox, oy, oz, -val ) ;
2566
2567 /*
2568 for ( i = -1 ; i < 2 ; i ++ )
2569 for ( j = -1 ; j < 2 ; j ++ )
2570 for ( k = -1 ; k < 2 ; k ++ )
2571 {
2572 nx = ox + i ;
2573 ny = oy + j ;
2574 nz = oz + k ;
2575 if ( getDataAt( nx, ny, nz ) == val )
2576 {
2577 if ( isSheetEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
2578 {
2579 rvalue ++ ;
2580 }
2581 }
2582 }
2583 */
2584
2585 for ( i = 0 ; i < 6 ; i ++ )
2586 {
2587 nx = ox + neighbor6[i][0] ;
2588 ny = oy + neighbor6[i][1] ;
2589 nz = oz + neighbor6[i][2] ;
2590
2591 if ( getDataAt(nx,ny,nz) >= 0 )
2592 {
2593 int num = getNumNeighbor6( nx, ny, nz ) ;
2594 if ( num > rvalue )
2595 {
2596 rvalue = num ;
2597 }
2598 }
2599 }
2600
2601
2602 setDataAt( ox, oy, oz, val ) ;
2603
2604 return rvalue + getNumNeighbor6( ox, oy, oz ) * 10 ;
2605 /*
2606 int v = (((rvalue + getNumNeighbor6( ox, oy, oz ) * 10) & 255) << 24) ;
2607 v |= ( ( ox & 255 ) << 16 ) ;
2608 v |= ( ( oy & 255 ) << 8 ) ;
2609 v |= ( ( oz & 255 ) ) ;
2610 return v ;
2611 */
2612
2613 }
wustl_mm::SkeletonMaker::Volume::getNumNeighbor6
int getNumNeighbor6(int ox, int oy, int oz)
Definition: volume.cpp:661

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

Referenced by getNumPotComplex2(), and surfaceSkeletonPres().

◆ getNumPotComplex2()

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

Definition at line 2615 of file volume.cpp.

2616 {
2617 return getNumPotComplex( ox, oy, oz ) ;
2618
2619 //int i, j, k ;
2620 //double val = getDataAt( ox, oy, oz ) ;
2621 //if ( val <= 0 )
2622 //{
2623 // return 0 ;
2624 //}
2625
2626 //int rvalue = 0, nx, ny, nz ;
2627 //setDataAt( ox, oy, oz, -val ) ;
2628
2629 //for ( i = -1 ; i < 2 ; i ++ )
2630 // for ( j = -1 ; j < 2 ; j ++ )
2631 // for ( k = -1 ; k < 2 ; k ++ )
2632 // {
2633 // nx = ox + i ;
2634 // ny = oy + j ;
2635 // nz = oz + k ;
2636 // if ( getDataAt( nx, ny, nz ) == val )
2637 // {
2638 // if ( isHelixEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
2639 // {
2640 // rvalue ++ ;
2641 // }
2642 // }
2643 // }
2644
2645 //setDataAt( ox, oy, oz, val ) ;
2646
2647 //return rvalue + getNumNeighbor6( ox, oy, oz ) * 30 ;
2648 }
wustl_mm::SkeletonMaker::Volume::getNumPotComplex
int getNumPotComplex(int ox, int oy, int oz)
Definition: volume.cpp:2548

References getNumPotComplex().

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

◆ getOriginX()

float Volume::getOriginX ( )

Definition at line 93 of file volume.cpp.

93 {
94 return volData->GetOriginX();
95 }
wustl_mm::SkeletonMaker::VolumeData::GetOriginX
float GetOriginX()
Definition: volume_data.cpp:136

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

◆ getOriginY()

float Volume::getOriginY ( )

Definition at line 97 of file volume.cpp.

97 {
98 return volData->GetOriginY();
99 }
wustl_mm::SkeletonMaker::VolumeData::GetOriginY
float GetOriginY()
Definition: volume_data.cpp:140

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

◆ getOriginZ()

float Volume::getOriginZ ( )

Definition at line 101 of file volume.cpp.

101 {
102 return volData->GetOriginZ();
103 }
wustl_mm::SkeletonMaker::VolumeData::GetOriginZ
float GetOriginZ()
Definition: volume_data.cpp:144

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

◆ getSizeX()

int Volume::getSizeX ( )

Definition at line 36 of file volume.cpp.

36 {
37 return volData->GetSizeX();
38 }
wustl_mm::SkeletonMaker::VolumeData::GetSizeX
int GetSizeX()
Definition: volume_data.cpp:78

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

◆ getSizeY()

int Volume::getSizeY ( )

Definition at line 40 of file volume.cpp.

40 {
41 return volData->GetSizeY();
42 }
wustl_mm::SkeletonMaker::VolumeData::GetSizeY
int GetSizeY()
Definition: volume_data.cpp:82

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

◆ getSizeZ()

int Volume::getSizeZ ( )

Definition at line 44 of file volume.cpp.

44 {
45 return volData->GetSizeZ();
46 }
wustl_mm::SkeletonMaker::VolumeData::GetSizeZ
int GetSizeZ()
Definition: volume_data.cpp:86

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

◆ getSpacingX()

float Volume::getSpacingX ( )

Definition at line 81 of file volume.cpp.

81 {
82 return volData->GetSpacingX();
83 }
wustl_mm::SkeletonMaker::VolumeData::GetSpacingX
float GetSpacingX()
Definition: volume_data.cpp:107

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

◆ getSpacingY()

float Volume::getSpacingY ( )

Definition at line 85 of file volume.cpp.

85 {
86 return volData->GetSpacingY();
87 }
wustl_mm::SkeletonMaker::VolumeData::GetSpacingY
float GetSpacingY()
Definition: volume_data.cpp:111

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

◆ getSpacingZ()

float Volume::getSpacingZ ( )

Definition at line 89 of file volume.cpp.

89 {
90 return volData->GetSpacingZ();
91 }
wustl_mm::SkeletonMaker::VolumeData::GetSpacingZ
float GetSpacingZ()
Definition: volume_data.cpp:115

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

◆ getVolumeData()

VolumeData * Volume::getVolumeData ( )

Definition at line 69 of file volume.cpp.

69 {
70 return volData;
71 }

References volData.

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

◆ hasCell()

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

Definition at line 1016 of file volume.cpp.

1016 {
1017 for ( int i = 0 ; i < 2 ; i ++ )
1018 for ( int j = 0 ; j < 2 ; j ++ )
1019 for ( int k = 0 ; k < 2 ; k ++ ) {
1020 if ( getDataAt( ox + i, oy + j, oz + k ) < 0 ) {
1021 return 0 ;
1022 }
1023 }
1024 return 1 ;
1025 }

References getDataAt().

Referenced by hasCompleteSheet(), and markCellFace().

◆ hasCompleteHelix() [1/2]

int Volume::hasCompleteHelix ( int  ox,
int  oy,
int  oz 
)

Definition at line 1456 of file volume.cpp.

1457 {
1458 // Returns 1 if it has a complete helix
1459 int i ;
1460 int c1 = 0;
1461 int nx, ny, nz ;
1462 int j ;
1463
1464 for ( i = 0 ; i < 6 ; i ++ )
1465 {
1466 nx = ox + neighbor6[i][0] ;
1467 ny = oy + neighbor6[i][1] ;
1468 nz = oz + neighbor6[i][2] ;
1469 if ( getDataAt( nx, ny, nz ) >= 0 )
1470 {
1471 c1 ++ ;
1472 j = i ;
1473 }
1474
1475 }
1476
1477 if ( c1 > 1 ) // || c1 == 0 )
1478 {
1479 return 1 ;
1480 }
1481
1482 return 0 ;
1483
1484 /*
1485 ox = ox + neighbor6[j][0] ;
1486 oy = oy + neighbor6[j][1] ;
1487 oz = oz + neighbor6[j][2] ;
1488 c1 = 0 ;
1489 for ( i = 0 ; i < 6 ; i ++ )
1490 {
1491 nx = ox + neighbor6[i][0] ;
1492 ny = oy + neighbor6[i][1] ;
1493 nz = oz + neighbor6[i][2] ;
1494 if ( getDataAt( nx, ny, nz ) >= 0 )
1495 {
1496 c1 ++ ;
1497 }
1498
1499 }
1500
1501 if ( c1 > 1 )
1502 {
1503 return 0 ;
1504 }
1505 else
1506 {
1507 return 1 ;
1508 }
1509 */
1510 }

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

Referenced by erodeHelix().

◆ hasCompleteHelix() [2/2]

int Volume::hasCompleteHelix ( int  ox,
int  oy,
int  oz,
Volume *  fvol 
)

Definition at line 1512 of file volume.cpp.

1513 {
1514
1515 int i ;
1516 int c1 = 0;
1517 int nx, ny, nz ;
1518 int j ;
1519
1520 for ( i = 0 ; i < 6 ; i ++ )
1521 {
1522 nx = ox + neighbor6[i][0] ;
1523 ny = oy + neighbor6[i][1] ;
1524 nz = oz + neighbor6[i][2] ;
1525 if ( getDataAt( nx, ny, nz ) >= 0 )
1526 {
1527 if ( i % 2 == 0 )
1528 {
1529 nx = ox ;
1530 ny = oy ;
1531 nz = oz ;
1532 }
1533
1534 int val = (int)fvol->getDataAt( nx, ny, nz) ;
1535 if ( (( val >> ( 2 * ( i / 2 ) ) ) & 1) == 0 )
1536 {
1537 c1 ++ ;
1538 j = i ;
1539 }
1540 }
1541
1542 }
1543
1544 if ( c1 > 1 )
1545 {
1546 return 1 ;
1547 }
1548
1549 return 0 ;
1550
1551 /*
1552 ox = ox + neighbor6[j][0] ;
1553 oy = oy + neighbor6[j][1] ;
1554 oz = oz + neighbor6[j][2] ;
1555 c1 = 0 ;
1556 for ( i = 0 ; i < 6 ; i ++ )
1557 {
1558 nx = ox + neighbor6[i][0] ;
1559 ny = oy + neighbor6[i][1] ;
1560 nz = oz + neighbor6[i][2] ;
1561 if ( getDataAt( nx, ny, nz ) >= 0 )
1562 {
1563 c1 ++ ;
1564 }
1565
1566 }
1567
1568 if ( c1 > 1 )
1569 {
1570 return 0 ;
1571 }
1572 else
1573 {
1574 return 1 ;
1575 }
1576 */
1577 }

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

◆ hasCompleteSheet() [1/2]

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

Definition at line 1322 of file volume.cpp.

1322 {
1323 // Returns 1 if it lies in the middle of a sheet
1324 int temp = countIntEuler( ox, oy, oz ) ;
1325 if ( temp > 0 )
1326 {
1327 return 1 ;
1328 }
1329 else
1330 {
1331 return 0 ;
1332 }
1333 }
wustl_mm::SkeletonMaker::Volume::countIntEuler
int countIntEuler(int ox, int oy, int oz)
Definition: volume.cpp:2895

References countIntEuler().

◆ hasCompleteSheet() [2/2]

int Volume::hasCompleteSheet ( int  ox,
int  oy,
int  oz,
Volume *  fvol 
)

Definition at line 1174 of file volume.cpp.

1174 {
1175 int i, j, k ;
1176 int nx, ny, nz ;
1177
1178 int edge[6] = { 0,0,0,0,0,0 } ;
1179 int faceflag[ 12 ] ;
1180 int tot = 0 ;
1181 int cellflag[ 8 ] ;
1182
1183 int ct = 0 ;
1184 for ( i = -1 ; i < 1 ; i ++ )
1185 for ( j = -1 ; j < 1 ; j ++ )
1186 for ( k = -1 ; k < 1 ; k ++ )
1187 {
1188 if ( hasCell( ox + i, oy + j, oz + k ) )
1189 {
1190 cellflag[ ct ] = 1 ;
1191 }
1192 else
1193 {
1194 cellflag[ ct ] = 0 ;
1195 }
1196 ct ++ ;
1197 }
1198
1199 for ( i = 0 ; i < 12 ; i ++ )
1200 {
1201 faceflag[ i ] = 1 ;
1202 for ( j = 0 ; j < 4 ; j ++ )
1203 {
1204 nx = ox + sheetNeighbor[i][j][0] ;
1205 ny = oy + sheetNeighbor[i][j][1] ;
1206 nz = oz + sheetNeighbor[i][j][2] ;
1207
1208 if ( getDataAt( nx, ny, nz ) < 0 )
1209 {
1210 faceflag[ i ] = 0 ;
1211 break ;
1212 }
1213 }
1214
1215 if ( faceflag[ i ] )
1216 {
1217 if ( cellflag[ faceCells[i][0] ] ^ cellflag[ faceCells[i][1] ] )
1218 {
1219 int v1 = (int)( fvol->getDataAt(
1220 ox - 1 + (( faceCells[i][0] >> 2 ) & 1 ),
1221 oy - 1 + (( faceCells[i][0] >> 1 ) & 1 ),
1222 oz - 1 + (( faceCells[i][0] ) & 1)) ) ;
1223 int v2 = (int)( fvol->getDataAt(
1224 ox - 1 + (( faceCells[i][1] >> 2 ) & 1 ),
1225 oy - 1 + (( faceCells[i][1] >> 1 ) & 1 ),
1226 oz - 1 + (( faceCells[i][1] ) & 1)) ) ;
1227 if ( ((v1 >> (2 * (2 - i/4))) & 1) ||
1228 ((v2 >> (2 * (2 - i/4) + 1 )) & 1) )
1229 {
1230 faceflag[ i ] = 0 ;
1231 }
1232 }
1233 }
1234
1235 if ( faceflag[ i ] )
1236 {
1237 int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
1238 edge[ e0 ] ++ ;
1239 edge[ e1 ] ++ ;
1240 tot ++ ;
1241 }
1242 }
1243
1244 // Removing 1s
1245 int numones = 0 ;
1246 for ( i = 0 ; i < 6 ; i ++ )
1247 {
1248 if ( edge[ i ] == 1 )
1249 {
1250 numones ++ ;
1251 }
1252 }
1253 while( numones > 0 )
1254 {
1255 int e = 0;
1256 for ( i = 0 ; i < 6 ; i ++ )
1257 {
1258 if ( edge[ i ] == 1 )
1259 {
1260 e = i ;
1261 break ;
1262 }
1263 }
1264 /*
1265 if ( edge[ e ] != 1 )
1266 {
1267 printf("Wrong Again!********\n") ;
1268 }
1269 */
1270
1271 int f, e2 ;
1272 for ( j = 0 ; j < 4 ; j ++ )
1273 {
1274 f = edgeFaces[ e ][ j ] ;
1275 if ( faceflag[ f ] )
1276 {
1277 break ;
1278 }
1279 }
1280
1281 /*
1282 if ( faceflag[ f ] == 0 )
1283 {
1284 printf("Wrong!********\n") ;
1285 }
1286 */
1287
1288 if ( faceEdges[ f ][ 0 ] == e )
1289 {
1290 e2 = faceEdges[ f ][ 1 ] ;
1291 }
1292 else
1293 {
1294 e2 = faceEdges[ f ][ 0 ] ;
1295 }
1296
1297
1298 edge[ e ] -- ;
1299 numones -- ;
1300 edge[ e2 ] -- ;
1301 faceflag[ f ] = 0 ;
1302 tot -- ;
1303
1304 if ( edge[ e2 ] == 1 )
1305 {
1306 numones ++ ;
1307 }
1308 else if ( edge[ e2 ] == 0 )
1309 {
1310 numones -- ;
1311 }
1312 }
1313
1314 if ( tot > 0 )
1315 {
1316 return 1 ;
1317 }
1318
1319 return 0 ;
1320 }
wustl_mm::SkeletonMaker::Volume::hasCell
int hasCell(int ox, int oy, int oz)
Definition: volume.cpp:1016
wustl_mm::SkeletonMaker::edgeFaces
const int edgeFaces[6][4]
Definition: volume.h:63
wustl_mm::SkeletonMaker::faceCells
const int faceCells[12][2]
Definition: volume.h:54
wustl_mm::SkeletonMaker::faceEdges
const int faceEdges[12][2]
Definition: volume.h:59

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

Referenced by erodeSheet(), and isSheetEnd().

◆ isFeatureFace()

int Volume::isFeatureFace ( int  ox,
int  oy,
int  oz 
)

Definition at line 2135 of file volume.cpp.

2136 {
2137 // return 1 ;
2138
2139 int i, j ;
2140 int nx, ny, nz ;
2141
2142 int faces = 12 ;
2143 for ( i = 0 ; i < 12 ; i ++ )
2144 {
2145 int ct = 0 ;
2146 for ( j = 0 ; j < 4 ; j ++ )
2147 {
2148 nx = ox + sheetNeighbor[i][j][0] ;
2149 ny = oy + sheetNeighbor[i][j][1] ;
2150 nz = oz + sheetNeighbor[i][j][2] ;
2151
2152 if ( getDataAt( nx, ny, nz ) < 0 )
2153 {
2154 ct = -1 ;
2155 break ;
2156 }
2157 else if ( getNumNeighbor6( nx, ny, nz ) == 6 )
2158 {
2159 ct = -1 ;
2160 break ;
2161 }
2162 // else if ( getDataAt( nx, ny, nz ) == 0 )
2163 // {
2164 // ct ++ ;
2165 // }
2166
2167
2168 }
2169 if ( ct == -1 || ct >= 1 )
2170 {
2171 faces -- ;
2172 }
2173 }
2174
2175 if ( faces > 0 )
2176 {
2177 return 1 ;
2178 }
2179 return 0 ;
2180
2181 }

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

Referenced by isSheetEnd().

◆ isHelixEnd() [1/2]

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

Definition at line 1789 of file volume.cpp.

1789 {
1790
1791 int i ;
1792 int c1 = 0 , c2 = 0 ;
1793 int nx, ny, nz ;
1794
1795 for ( i = 0 ; i < 6 ; i ++ )
1796 {
1797 nx = ox + neighbor6[i][0] ;
1798 ny = oy + neighbor6[i][1] ;
1799 nz = oz + neighbor6[i][2] ;
1800
1801 double val = getDataAt( nx, ny, nz ) ;
1802
1803 if ( val >= 0 )
1804 {
1805 c1 ++ ;
1806 if ( getNumNeighbor6(nx,ny,nz) < 6 ) // if ( val > 0 && val < MAX_ERODE )
1807 {
1808 c2 ++ ;
1809 }
1810 }
1811
1812 }
1813
1814 if ( c1 == 1 && c2 == 1 )
1815 {
1816 return 1 ;
1817 }
1818
1819 return 0 ;
1820 }

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

◆ isHelixEnd() [2/2]

int Volume::isHelixEnd ( int  ox,
int  oy,
int  oz,
Volume *  nvol 
)

Definition at line 1579 of file volume.cpp.

1579 {
1580 // Returns 1 if it is a curve endpoint
1581 int i ;
1582 int c1 = 0 , c2 = 0 ;
1583 int nx, ny, nz ;
1584
1585 for ( i = 0 ; i < 6 ; i ++ )
1586 {
1587 nx = ox + neighbor6[i][0] ;
1588 ny = oy + neighbor6[i][1] ;
1589 nz = oz + neighbor6[i][2] ;
1590
1591 double val = getDataAt( nx, ny, nz ) ;
1592
1593 if ( val >= 0 )
1594 {
1595 c1 ++ ;
1596 if ( val > 0 && val < MAX_ERODE && nvol->getDataAt( nx, ny, nz ) == 0 )
1597 {
1598 c2 ++ ;
1599 }
1600 }
1601
1602 }
1603
1604 if ( c1 == 1 && c2 == 1 )
1605 {
1606 return 1 ;
1607 }
1608
1609 return 0 ;
1610 }

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

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

◆ isPiercable()

int Volume::isPiercable ( int  ox,
int  oy,
int  oz 
)

Definition at line 2381 of file volume.cpp.

2382 {
2383 /* Test if this is a simple voxel */
2384 // int flag = 0 ;
2385 double vox[3][3][3] ;
2386
2387 int i, j, k ;
2388 for ( i = -1 ; i < 2 ; i ++ )
2389 for ( j = -1 ; j < 2 ; j ++ )
2390 for ( k = -1 ; k < 2 ; k ++ )
2391 {
2392 double tval = getDataAt( ox + i, oy + j, oz + k ) ;
2393
2394 /*
2395 if ( tval == 0 || tval > (va )
2396 {
2397 flag = 1 ;
2398 }
2399 */
2400 /*
2401 if ( tval < 0 && tval == - curwid )
2402 {
2403 printf("Here %d", (int)-tval) ;
2404 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
2405 }
2406 else
2407 */
2408 {
2409 vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
2410 }
2411 }
2412
2413 /* Debugging
2414 printf("{") ;
2415 for ( i = 0 ; i < 3 ; i ++ )
2416 {
2417 if ( i ) printf(",") ;
2418 printf("{") ;
2419 for ( j = 0 ; j < 3 ; j ++ )
2420 {
2421 if ( j ) printf(",") ;
2422 printf("{") ;
2423 for ( k = 0 ; k < 3 ; k ++ )
2424 {
2425 if ( k ) printf(",") ;
2426 printf("%d", (vox[i][j][k] >=0 ? 1: 0));
2427 }
2428 printf("}") ;
2429 }
2430 printf("}") ;
2431 }
2432 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
2433 */
2434
2435 if ( countInt( vox ) == 1 && countExt( vox ) != 1 )
2436 {
2437 return 1 ;
2438 }
2439 else
2440 {
2441 return 0 ;
2442 }
2443 }
wustl_mm::SkeletonMaker::Volume::countInt
int countInt(double vox[3][3][3])
Definition: volume.cpp:2859
wustl_mm::SkeletonMaker::Volume::countExt
int countExt(double vox[3][3][3])
Definition: volume.cpp:2838

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

Referenced by curveSkeleton().

◆ isSheetEnd() [1/2]

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

Definition at line 2238 of file volume.cpp.

2239 {
2240 return ( ( hasCompleteSheet(ox,oy,oz) == 0 ) && isFeatureFace(ox,oy,oz) ) ;
2241
2243 //
2244 //int i, j, k ;
2245 //int nx, ny, nz ;
2246
2247 //double vox[3][3][3] ;
2248 //for ( i = -1 ; i < 2 ; i ++ )
2249 // for ( j = -1 ; j < 2 ; j ++ )
2250 // for ( k = -1 ; k < 2 ; k ++ )
2251 // {
2252 // vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ;
2253 // }
2254
2255 //int edge[6] = { 4,4,4,4,4,4 } ;
2256 //int edge2[6] = { 4,4,4,4,4,4 } ;
2257
2258 //for ( i = 0 ; i < 12 ; i ++ )
2259 //{
2260 // for ( j = 0 ; j < 4 ; j ++ )
2261 // {
2262 // nx = 1 + sheetNeighbor[i][j][0] ;
2263 // ny = 1 + sheetNeighbor[i][j][1] ;
2264 // nz = 1 + sheetNeighbor[i][j][2] ;
2265
2266 // if ( vox[nx][ny][nz] < 0 )
2267 // {
2268 // edge[ faceEdges[ i ][ 0 ] ] -- ;
2269 // edge[ faceEdges[ i ][ 1 ] ] -- ;
2270 // break ;
2271 // }
2272 // }
2273
2274 // for ( j = 0 ; j < 4 ; j ++ )
2275 // {
2276 // nx = 1 + sheetNeighbor[i][j][0] ;
2277 // ny = 1 + sheetNeighbor[i][j][1] ;
2278 // nz = 1 + sheetNeighbor[i][j][2] ;
2279
2280 // if ( vox[nx][ny][nz] <= 0 )
2281 // {
2282 // edge2[ faceEdges[ i ][ 0 ] ] -- ;
2283 // edge2[ faceEdges[ i ][ 1 ] ] -- ;
2284 // break ;
2285 // }
2286 // }
2287 //}
2288
2289 //
2291 //for ( i = 0 ; i < 6 ; i ++ )
2292 //{
2293 // nx = 1 + neighbor6[i][0] ;
2294 // ny = 1 + neighbor6[i][1] ;
2295 // nz = 1 + neighbor6[i][2] ;
2296 // if ( edge[i] == 0 && vox[nx][ny][nz] >= 0 )
2297 // {
2298 // return 0 ;
2299 // }
2300 //}
2301 //*/
2302 //
2303
2304
2305 //for ( i = 0 ; i < 6 ; i ++ )
2306 //{
2307 // if ( edge[ i ] == 1 && edge2[ i ] == 1 )
2308 // {
2309 // return 1 ;
2310 // }
2311 //}
2312
2313 //return 0 ;
2314 }
wustl_mm::SkeletonMaker::Volume::isFeatureFace
int isFeatureFace(int ox, int oy, int oz)
Definition: volume.cpp:2135

References hasCompleteSheet(), and isFeatureFace().

◆ isSheetEnd() [2/2]

int Volume::isSheetEnd ( int  ox,
int  oy,
int  oz,
Volume *  nvol 
)

Definition at line 1821 of file volume.cpp.

1822 {
1823 // Returns 1 if it contains a sheet boundary. Noise-resistant
1824 int i, j ;
1825 int nx, ny, nz ;
1826
1827 int edge[6] = { 0,0,0,0,0,0 } ;
1828 int faceflag[ 12 ] ;
1829 int hasFeatureFace = 0 ;
1830 int tot = 0 ;
1831
1832 for ( i = 0 ; i < 12 ; i ++ )
1833 {
1834 faceflag[ i ] = 1 ;
1835 int hasFeature = 1 ;
1836
1837 for ( j = 0 ; j < 4 ; j ++ )
1838 {
1839 nx = ox + sheetNeighbor[i][j][0] ;
1840 ny = oy + sheetNeighbor[i][j][1] ;
1841 nz = oz + sheetNeighbor[i][j][2] ;
1842
1843 if ( getDataAt( nx, ny, nz ) == 0 || nvol->getDataAt( nx, ny, nz ) == 1 )
1844 {
1845 hasFeature = 0 ;
1846 }
1847 if ( getDataAt( nx, ny, nz ) < 0 )
1848 {
1849 faceflag[ i ] = 0 ;
1850 break ;
1851 }
1852 }
1853 if ( faceflag[ i ] == 1 && hasFeature )
1854 {
1855 hasFeatureFace ++ ;
1856 // return 0 ;
1857 }
1858
1859 if ( faceflag[ i ] )
1860 {
1861 int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
1862 edge[ e0 ] ++ ;
1863 edge[ e1 ] ++ ;
1864 tot ++ ;
1865 }
1866 }
1867
1868 if ( tot == 0 || hasFeatureFace == 0 )
1869 {
1870 return 0 ;
1871 }
1872
1873 // Removing 1s
1874 int numones = 0 ;
1875 for ( i = 0 ; i < 6 ; i ++ )
1876 {
1877 if ( edge[ i ] == 1 )
1878 {
1879 numones ++ ;
1880 }
1881 }
1882 while( numones > 0 )
1883 {
1884 int e = 0;
1885 for ( i = 0 ; i < 6 ; i ++ )
1886 {
1887 if ( edge[ i ] == 1 )
1888 {
1889 e = i ;
1890 break ;
1891 }
1892 }
1893 /*
1894 if ( edge[ e ] != 1 )
1895 {
1896 printf("Wrong Again!********\n") ;
1897 }
1898 */
1899
1900 int f, e2 ;
1901 for ( j = 0 ; j < 4 ; j ++ )
1902 {
1903 f = edgeFaces[ e ][ j ] ;
1904 if ( faceflag[ f ] )
1905 {
1906 break ;
1907 }
1908 }
1909
1910 /*
1911 if ( faceflag[ f ] == 0 )
1912 {
1913 printf("Wrong!********\n") ;
1914 }
1915 */
1916
1917 if ( faceEdges[ f ][ 0 ] == e )
1918 {
1919 e2 = faceEdges[ f ][ 1 ] ;
1920 }
1921 else
1922 {
1923 e2 = faceEdges[ f ][ 0 ] ;
1924 }
1925
1926
1927 edge[ e ] -- ;
1928 numones -- ;
1929 edge[ e2 ] -- ;
1930 faceflag[ f ] = 0 ;
1931 tot -- ;
1932
1933 if ( edge[ e2 ] == 1 )
1934 {
1935 numones ++ ;
1936 }
1937 else if ( edge[ e2 ] == 0 )
1938 {
1939 numones -- ;
1940 }
1941 }
1942
1943 if ( tot > 0 )
1944 {
1945 return 0 ;
1946 }
1947
1948 return 1 ;
1949 }

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

Referenced by surfaceSkeletonPres().

◆ isSimple()

int Volume::isSimple ( int  ox,
int  oy,
int  oz 
)

Definition at line 2316 of file volume.cpp.

2317 {
2318 /* Test if this is a simple voxel */
2319 // int flag = 0 ;
2320 double vox[3][3][3] ;
2321
2322 int i, j, k ;
2323 for ( i = -1 ; i < 2 ; i ++ )
2324 for ( j = -1 ; j < 2 ; j ++ )
2325 for ( k = -1 ; k < 2 ; k ++ )
2326 {
2327 double tval = getDataAt( ox + i, oy + j, oz + k ) ;
2328
2329 /*
2330 if ( tval == 0 || tval > (va )
2331 {
2332 flag = 1 ;
2333 }
2334 */
2335 /*
2336 if ( tval < 0 && tval == - curwid )
2337 {
2338 printf("Here %d", (int)-tval) ;
2339 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
2340 }
2341 else
2342 */
2343 {
2344 vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
2345 }
2346 }
2347
2348 /* Debugging
2349 printf("{") ;
2350 for ( i = 0 ; i < 3 ; i ++ )
2351 {
2352 if ( i ) printf(",") ;
2353 printf("{") ;
2354 for ( j = 0 ; j < 3 ; j ++ )
2355 {
2356 if ( j ) printf(",") ;
2357 printf("{") ;
2358 for ( k = 0 ; k < 3 ; k ++ )
2359 {
2360 if ( k ) printf(",") ;
2361 printf("%d", (vox[i][j][k] >=0 ? 1: 0));
2362 }
2363 printf("}") ;
2364 }
2365 printf("}") ;
2366 }
2367 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
2368 */
2369
2370 if ( countInt( vox ) == 1 && countExt( vox ) == 1 )
2371 {
2372 return 1 ;
2373 }
2374 else
2375 {
2376 return 0 ;
2377 }
2378 }

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

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

◆ markCellFace()

Volume * Volume::markCellFace ( )

Definition at line 1027 of file volume.cpp.

1027 {
1028 int i, j, k ;
1029 Volume* fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
1030
1031 //return fvol ;
1032
1033 for ( i = 0 ; i < getSizeX() ; i ++ )
1034 for ( j = 0 ; j < getSizeY() ; j ++ )
1035 for ( k = 0 ; k < getSizeZ() ; k ++ )
1036 {
1037 if( getDataAt(i,j,k) >= 0 )
1038 {
1039 if ( hasCell(i,j,k) )
1040 {
1041 for ( int m = 0 ; m < 6 ; m ++ )
1042 {
1043 int nx = i + neighbor6[m][0] ;
1044 int ny = j + neighbor6[m][1] ;
1045 int nz = k + neighbor6[m][2] ;
1046 if ( ! hasCell(nx,ny,nz) )
1047 {
1048 fvol->setDataAt(i,j,k,(double)(1<<m)) ;
1049 break ;
1050 }
1051 }
1052 }
1053 }
1054 }
1055
1056
1057 return fvol ;
1058 }

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

Referenced by erodeSheet().

◆ pad()

void Volume::pad ( int  padBy,
double  padValue 
)
  • Next, smooth* Next, smooth *‍/

Definition at line 273 of file volume.cpp.

273 {
274 volData->Pad(padBy, padValue);
275 }
wustl_mm::SkeletonMaker::VolumeData::Pad
void Pad(int padBy, double padValue)
Definition: volume_data.cpp:196

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

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

◆ setDataAt() [1/2]

void Volume::setDataAt ( int  index,
double  d 
)

Definition at line 56 of file volume.cpp.

56 {
57 volData->SetDataAt(index, (float)d);
58 }
wustl_mm::SkeletonMaker::VolumeData::SetDataAt
void SetDataAt(int x, int y, int z, float value)
Definition: volume_data.cpp:188

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

◆ setDataAt() [2/2]

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

Definition at line 52 of file volume.cpp.

52 {
53 volData->SetDataAt(x, y, z, (float)d);
54 }

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

◆ setOrigin()

void Volume::setOrigin ( float  orgX,
float  orgY,
float  orgZ 
)

Definition at line 77 of file volume.cpp.

77 {
78 volData->SetOrigin(orgX, orgY, orgZ);
79 }
wustl_mm::SkeletonMaker::VolumeData::SetOrigin
void SetOrigin(float originX, float originY, float originZ)
Definition: volume_data.cpp:177

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

◆ setSpacing()

void Volume::setSpacing ( float  spx,
float  spy,
float  spz 
)

Definition at line 73 of file volume.cpp.

73 {
74 volData->SetSpacing(spx, spy, spz);
75 }
wustl_mm::SkeletonMaker::VolumeData::SetSpacing
void SetSpacing(float spacingX, float spacingY, float spacingZ)
Definition: volume_data.cpp:171

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

◆ skeleton() [1/2]

void Volume::skeleton ( float  thr,
int  off 
)

Definition at line 5089 of file volume.cpp.

5090 {
5091 int i, j, k ;
5092 // First, threshold the volume
5093 #ifdef VERBOSE
5094 printf("Thresholding the volume to -1/0...\n") ;
5095 #endif
5096 threshold( thr, -1, 0 ) ;
5097
5098 // Next, apply convergent erosion
5099 // by preserving: complex nodes, curve end-points, and sheet points
5100
5101 // Next, initialize the linked queue
5102 #ifdef VERBOSE
5103 printf("Initializing queue...\n") ;
5104 #endif
5105 GridQueue2* queue2 = new GridQueue2( ) ;
5106 GridQueue2* queue = new GridQueue2( ) ;
5107
5108 for ( i = 0 ; i < getSizeX() ; i ++ )
5109 for ( j = 0 ; j < getSizeY() ; j ++ )
5110 for ( k = 0 ; k < getSizeZ() ; k ++ )
5111 {
5112 if ( getDataAt( i, j, k ) >= 0 )
5113 {
5114 {
5115 for ( int m = 0 ; m < 6 ; m ++ )
5116 {
5117 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
5118 {
5119 setDataAt( i, j, k, 1 ) ;
5120 queue2->prepend( i, j, k ) ;
5121 break ;
5122 }
5123 }
5124 }
5125 }
5126 }
5127 int wid = 0 ;
5128 #ifdef VERBOSE
5129 printf("Total %d nodes\n", queue2->getNumElements() ) ;
5130 printf("Start erosion to %d...\n", wid) ;
5131 #endif
5132
5133 // Perform erosion
5134 gridQueueEle* ele ;
5135 int ox, oy, oz ;
5136
5137 while( 1 )
5138 {
5139 wid ++ ;
5140
5141 // At the start of each iteration,
5142 // queue2 holds all the nodes for this layer
5143 // queue is empty
5144
5145 int numComplex = 0, numSimple = 0 ;
5146 #ifdef VERBOSE
5147 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), wid) ;
5148 #endif
5149
5150 // Rename queue2 to be queue,
5151 // Clear queue2
5152 // From now on, queue2 holds nodes of next level
5153 delete queue ;
5154 queue = queue2 ;
5155 queue2 = new GridQueue2( ) ;
5156
5157 // Next, start queue iteration
5158 queue->reset() ;
5159 ele = queue->getNext();
5160 while ( ele != NULL )
5161 {
5162 ox = ele->x ;
5163 oy = ele->y ;
5164 oz = ele->z ;
5165 // delete ele ;
5166
5167 // Check simple
5168 if ( ! isSimple( ox, oy, oz ) )
5169 {
5170 // Complex, set to next layer
5171 queue2->prepend( ox, oy, oz ) ;
5172 numComplex ++ ;
5173 }
5174 /*
5175 else if ( ox == off || oy == off || oz == off ||
5176 ox == getSizeX() - off - 1 || oy == getSizeY() - off - 1 || oz == getSizeZ() - off - 1 )
5177 {
5178 // Wall, don't erode, set to next layer
5179 queue2->prepend( ox, oy, oz ) ;
5180 numComplex ++ ;
5181 }
5182 */
5183 else
5184 {
5185 setDataAt( ox, oy, oz, -1 ) ;
5186 numSimple ++ ;
5187
5188 for ( int m = 0 ; m < 6 ; m ++ )
5189 {
5190 int nx = ox + neighbor6[m][0] ;
5191 int ny = oy + neighbor6[m][1] ;
5192 int nz = oz + neighbor6[m][2] ;
5193 if ( getDataAt( nx, ny, nz ) == 0 )
5194 {
5195 setDataAt( nx, ny, nz, 1 ) ;
5196 queue2->prepend( nx, ny, nz ) ;
5197 }
5198 }
5199
5200 }
5201
5202 ele = queue->remove() ;
5203 }
5204 #ifdef VERBOSE
5205 printf("Level %d: %d complex, %d simple\n", wid, numComplex, numSimple) ;
5206 #endif
5207
5208 if ( numSimple == 0 )
5209 {
5210 break ;
5211 }
5212 }
5213
5214 // Finally, clean up
5215 delete queue;
5216 delete queue2;
5217 #ifdef VERBOSE
5218 printf("Thresholding the volume to 0/1...\n") ;
5219 #endif
5220 threshold( 0, 0, 1 ) ;
5221 }

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

◆ skeleton() [2/2]

void Volume::skeleton ( float  thr,
Volume *  svol,
Volume *  hvol 
)

insert them back into priority queue

  • Added for debugging *‍/ / Check simple* Adding ends *‍/
  • Thin the current volume while preserving voxels with values > highthr or <= lowthr in grayvol
  • Added for debugging *‍/ / Check simple
  • Adding ends *‍/

Definition at line 5692 of file volume.cpp.

5693 {
5694 int i, j, k ;
5695 // First, threshold the volume
5696 #ifdef VERBOSE
5697 printf("Thresholding the volume to -1/0...\n") ;
5698 #endif
5699 threshold( thr, -1, 0 ) ;
5700
5701 // Next, apply convergent erosion
5702 // by preserving: complex nodes, curve end-points, and sheet points
5703
5704 // Next, initialize the linked queue
5705 #ifdef VERBOSE
5706 printf("Initializing queue...\n") ;
5707 #endif
5708 GridQueue2* queue2 = new GridQueue2( ) ;
5709 GridQueue2* queue3 = new GridQueue2( ) ;
5710 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
5711
5712 for ( i = 0 ; i < getSizeX() ; i ++ )
5713 for ( j = 0 ; j < getSizeY() ; j ++ )
5714 for ( k = 0 ; k < getSizeZ() ; k ++ )
5715 {
5716 if ( getDataAt( i, j, k ) >= 0 )
5717 {
5718 if ( svol->getDataAt(i,j,k) > 0 || hvol->getDataAt(i,j,k) > 0 )
5719 {
5720 setDataAt( i, j, k, MAX_ERODE ) ;
5721 }
5722 else
5723 {
5724 for ( int m = 0 ; m < 6 ; m ++ )
5725 {
5726 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
5727 {
5728 setDataAt( i, j, k, 1 ) ;
5729 queue2->prepend( i, j, k ) ;
5730 break ;
5731 }
5732 }
5733 }
5734 }
5735 }
5736 int wid = MAX_ERODE ;
5737 #ifdef VERBOSE
5738 printf("Total %d nodes\n", queue2->getNumElements() ) ;
5739
5740
5741 // Perform erosion
5742 printf("Start erosion to %d...\n", wid) ;
5743 #endif
5744 gridQueueEle* ele ;
5745 gridPoint* gp ;
5746 int ox, oy, oz ;
5747 int score ;
5748 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
5749 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
5750 {
5751 scrvol->setDataAt( i, -1 ) ;
5752 }
5753
5754
5755 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
5756 {
5757 // At the start of each iteration,
5758 // queue2 holds all the nodes for this layer
5759 // queue3 and queue are empty
5760
5761 int numComplex = 0, numSimple = 0 ;
5762 #ifdef VERBOSE
5763 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
5764 #endif
5765
5766
5767 // Next,
5768 // Compute score for each node left in queue2
5769 // move them into priority queue
5770 queue2->reset() ;
5771 ele = queue2->getNext() ;
5772 while ( ele != NULL )
5773 {
5774 ox = ele->x ;
5775 oy = ele->y ;
5776 oz = ele->z ;
5777
5778 // Compute score
5779 score = getNumPotComplex2( ox, oy, oz ) ;
5780 scrvol->setDataAt( ox, oy, oz, score ) ;
5781
5782 // Push to queue
5783 gp = new gridPoint ;
5784 gp->x = ox ;
5785 gp->y = oy ;
5786 gp->z = oz ;
5787 // queue->add( gp, -score ) ;
5788 queue->add( gp, score ) ;
5789
5790 ele = queue2->remove() ;
5791 }
5792
5793 // Rename queue3 to be queue2,
5794 // Clear queue3
5795 // From now on, queue2 holds nodes of next level
5796 delete queue2 ;
5797 queue2 = queue3 ;
5798 queue3 = new GridQueue2( ) ;
5799
5800 // Next, start priority queue iteration
5801 while ( ! queue->isEmpty() )
5802 {
5803 // Retrieve the node with the highest score
5804 queue->remove( gp, score ) ;
5805 ox = gp->x ;
5806 oy = gp->y ;
5807 oz = gp->z ;
5808 delete gp ;
5809 // score = -score ;
5810
5811 // Ignore the node
5812 // if it has been processed before
5813 // or it has an updated score
5814 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
5815 {
5816 continue ;
5817 }
5818
5819 /* Added for debugging */
5820 // Check simple
5821 if ( ! isSimple( ox, oy, oz ) )
5822 {
5823 // Complex, set to next layer
5824 setDataAt( ox, oy, oz, curwid + 1 ) ;
5825 queue2->prepend( ox, oy, oz ) ;
5826 numComplex ++ ;
5827 }
5828 else
5829 {
5830 setDataAt( ox, oy, oz, -1 ) ;
5831 numSimple ++ ;
5832 }
5833 /* Adding ends */
5834
5835 // Move its neighboring unvisited node to queue2
5836 for ( int m = 0 ; m < 6 ; m ++ )
5837 {
5838 int nx = ox + neighbor6[m][0] ;
5839 int ny = oy + neighbor6[m][1] ;
5840 int nz = oz + neighbor6[m][2] ;
5841 if ( getDataAt( nx, ny, nz ) == 0 )
5842 {
5843 setDataAt( nx, ny, nz, curwid + 1 ) ;
5844 queue2->prepend( nx, ny, nz ) ;
5845 }
5846 }
5847
5848 // Update scores for nodes in its 5x5 neighborhood
5849 // insert them back into priority queue
5850
5851 for ( i = -2 ; i < 3 ;i ++ )
5852 for ( j = -2 ; j < 3 ; j ++ )
5853 for ( k = -2 ; k < 3 ; k ++ )
5854 {
5855 int nx = ox + i ;
5856 int ny = oy + j ;
5857 int nz = oz + k ;
5858
5859 if ( getDataAt( nx, ny, nz ) == curwid )
5860 {
5861 // Compute score
5862 score = getNumPotComplex2( nx, ny, nz ) ;
5863
5864 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
5865 {
5866 // printf("Update\n") ;
5867 scrvol->setDataAt( nx, ny, nz, score ) ;
5868 // Push to queue
5869 gp = new gridPoint ;
5870 gp->x = nx ;
5871 gp->y = ny ;
5872 gp->z = nz ;
5873 // queue->add( gp, -score ) ;
5874 queue->add( gp, score ) ;
5875 }
5876 }
5877 }
5878
5879
5880 }
5881
5882 #ifdef VERBOSE
5883 printf("%d complex, %d simple\n", numComplex, numSimple) ;
5884 #endif
5885
5886 if ( numSimple == 0 )
5887 {
5888 delete queue2 ;
5889 break ;
5890 }
5891 }
5892
5893 // Finally, clean up
5894 #ifdef VERBOSE
5895 printf("Thresholding the volume to 0/1...\n") ;
5896 #endif
5897 threshold( 0, 0, 1 ) ;
5898 delete scrvol;
5899 delete queue;
5900 delete queue3;
5901 }

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

◆ surfaceSkeletonPres()

void Volume::surfaceSkeletonPres ( float  thr,
Volume *  preserve 
)

for ( int m = 0 ; m < 6 ; m ++ )

  • Deal with closed rings* Find local minimum: to help determine erosion level *‍/*
  • ‍************************************************************************‍/ / A sequential thinning method for extracting 6-connected skeletons / Properties: medial, topology preserving, shape preserving, noise resistance! /
    Parameters
    thrthreshold /
    type0 for curve preserving, 1 for surface preserving
    noise0 for no-noise-reduction, n for level-n noise reduction ************************************************************************‍/
  • Thin the current volume while preserving voxels with values > highthr or <= lowthr in grayvol
  • Extra step: classify nodes in queue2 into noise and non-noise nodes *‍/
  • Extra step: classify nodes in queue2 into noise and non-noise nodes *‍/
  • Commented for debugging
  • Commented for debugging
  • Added for debugging *‍/ / Check simple
  • Adding ends *‍/
  • Commented for debugging
  • Extra step: classify nodes in queue2 into noise and non-noise nodes *‍/
  • Commented for debugging
  • Commented for debugging
  • Added for debugging *‍/ / Check simple
  • Adding ends *‍/
  • Commented for debugging
  • Extra step: classify nodes in queue2 into noise and non-noise nodes *‍/
  • Commented for debugging
  • Commented for debugging
  • Added for debugging *‍/ / Check simple
  • Adding ends *‍/
  • Commented for debugging

Definition at line 8672 of file volume.cpp.

8673 {
8674 int i, j, k ;
8675 // First, threshold the volume
8676 #ifdef VERBOSE
8677 printf("Thresholding the volume to -MAX_ERODE/0...\n") ;
8678 #endif
8679 threshold( thr, -MAX_ERODE, 0 ) ;
8680
8681 // Next, initialize the linked queue
8682 #ifdef VERBOSE
8683 printf("Initializing queue...\n") ;
8684 #endif
8685 GridQueue2* queue2 = new GridQueue2( ) ;
8686 GridQueue2* queue3 = new GridQueue2( ) ;
8687 GridQueue2* queue4 = new GridQueue2( ) ;
8688
8689 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
8690
8691 for ( i = 0 ; i < getSizeX() ; i ++ )
8692 for ( j = 0 ; j < getSizeY() ; j ++ )
8693 for ( k = 0 ; k < getSizeZ() ; k ++ )
8694 {
8695 if ( getDataAt( i, j, k ) >= 0 ) {
8696 if(preserve->getDataAt(i, j, k) > 0) {
8697 setDataAt(i, j, k, MAX_ERODE);
8698 } else {
8699 for ( int m = 0 ; m < 6 ; m ++ )
8700 {
8701 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
8702 {
8703 // setDataAt( i, j, k, 1 ) ;
8704 queue2->prepend( i, j, k ) ;
8705 break ;
8706 }
8707 }
8708 }
8709 }
8710 }
8711 int wid = MAX_ERODE ;
8712 #ifdef VERBOSE
8713 printf("Total %d nodes\n", queue2->getNumElements() ) ;
8714 printf("Start erosion to %d...\n", wid) ;
8715 #endif
8716
8717
8718 // Perform erosion
8719 gridQueueEle* ele ;
8720 gridPoint* gp ;
8721 int ox, oy, oz ;
8722 int score ;
8723 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
8724 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
8725 {
8726 scrvol->setDataAt( i, -1 ) ;
8727 }
8728
8729 #ifdef NOISE_DIS_SHEET
8730 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
8731 #endif
8732
8733 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
8734 {
8735 // At the start of each iteration,
8736 // queue2 and queue4 holds all the nodes for this layer
8737 // queue3 and queue are empty
8738
8739 int numComplex = 0, numSimple = 0 ;
8740 #ifdef VERBOSE
8741 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
8742 #endif
8743
8744 /*
8745 We first need to assign curwid + 1 to every node in this layer
8746 */
8747 queue2->reset() ;
8748 ele = queue2->getNext() ;
8749 while ( ele != NULL )
8750 {
8751 ox = ele->x ;
8752 oy = ele->y ;
8753 oz = ele->z ;
8754
8755 if ( getDataAt(ox,oy,oz) == curwid )
8756 {
8757 ele = queue2->remove() ;
8758 }
8759 else
8760 {
8761 setDataAt(ox,oy,oz, curwid) ;
8762 ele = queue2->getNext() ;
8763 }
8764 }
8765 queue4->reset() ;
8766 ele = queue4->getNext() ;
8767 while ( ele != NULL )
8768 {
8769 ox = ele->x ;
8770 oy = ele->y ;
8771 oz = ele->z ;
8772
8773 queue2->prepend(ox,oy,oz) ;
8774 ele = queue4->remove() ;
8775 }
8776
8777 // Now queue2 holds all the nodes for this layer
8778
8779 #ifdef NOISE_DIS_SHEET
8780 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
8781 queue2->reset() ;
8782
8783 // First run
8784 int flag = 0 ;
8785 while ( ( ele = queue2->getNext() ) != NULL )
8786 {
8787 ox = ele->x ;
8788 oy = ele->y ;
8789 oz = ele->z ;
8790 if ( NOISE_DIS_SHEET <= 1 )
8791 {
8792 noisevol->setDataAt( ox, oy, oz, 0 ) ;
8793 }
8794 else
8795 {
8796 flag = 0 ;
8797 for ( int m = 0 ; m < 6 ; m ++ )
8798 {
8799 int nx = ox + neighbor6[m][0] ;
8800 int ny = oy + neighbor6[m][1] ;
8801 int nz = oz + neighbor6[m][2] ;
8802 if ( getDataAt( nx, ny, nz ) == 0 )
8803 {
8804 noisevol->setDataAt( ox, oy, oz, 1 ) ;
8805 flag = 1 ;
8806 break ;
8807 }
8808 }
8809 if ( ! flag )
8810 {
8811 noisevol->setDataAt( ox, oy, oz, 0 ) ;
8812 }
8813 }
8814 }
8815
8816 for ( int cur = 1 ; cur < NOISE_DIS_SHEET ; cur ++ )
8817 {
8818 queue2->reset() ;
8819 int count = 0 ;
8820
8821 while ( ( ele = queue2->getNext() ) != NULL )
8822 {
8823 ox = ele->x ;
8824 oy = ele->y ;
8825 oz = ele->z ;
8826
8827 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
8828 {
8829 continue ;
8830 }
8831
8832 flag = 0 ;
8833 for ( int m = 0 ; m < 6 ; m ++ )
8834 {
8835 int nx = ox + neighbor6[m][0] ;
8836 int ny = oy + neighbor6[m][1] ;
8837 int nz = oz + neighbor6[m][2] ;
8838 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
8839 {
8840 noisevol->setDataAt( ox, oy, oz, 1 ) ;
8841 count ++ ;
8842 break ;
8843 }
8844 }
8845 }
8846
8847 if ( count == 0 )
8848 {
8849 break ;
8850 }
8851 }
8852
8853
8854 #endif
8855
8856 /* Commented for debugging
8857
8858 // First,
8859 // check for complex nodes in queue2
8860 // move them from queue2 to queue3
8861 queue2->reset() ;
8862 ele = queue2->getNext() ;
8863 while ( ele != NULL )
8864 {
8865 ox = ele->x ;
8866 oy = ele->y ;
8867 oz = ele->z ;
8868
8869 // Check simple
8870 #ifndef NOISE_DIS_SHEET
8871 if ( isSheetEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
8872 #else
8873 if ( isSheetEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
8874 #endif
8875 {
8876 // Complex, set to next layer
8877 setDataAt( ox, oy, oz, curwid + 1 ) ;
8878 queue3->prepend( ox, oy, oz ) ;
8879 ele = queue2->remove() ;
8880
8881 numComplex ++ ;
8882 }
8883 else
8884 {
8885 ele = queue2->getNext() ;
8886 }
8887 }
8888 */
8889
8890
8891 // Next,
8892 // Compute score for each node left in queue2
8893 // move them into priority queue
8894 queue2->reset() ;
8895 ele = queue2->getNext() ;
8896 while ( ele != NULL )
8897 {
8898 ox = ele->x ;
8899 oy = ele->y ;
8900 oz = ele->z ;
8901
8902 // Compute score
8903 score = getNumPotComplex( ox, oy, oz ) ;
8904 scrvol->setDataAt( ox, oy, oz, score ) ;
8905
8906 // Push to queue
8907 gp = new gridPoint ;
8908 gp->x = ox ;
8909 gp->y = oy ;
8910 gp->z = oz ;
8911 // queue->add( gp, -score ) ;
8912 queue->add( gp, score ) ;
8913
8914 ele = queue2->remove() ;
8915 }
8916
8917 // Rename queue3 to be queue2,
8918 // Clear queue3
8919 // From now on, queue2 holds nodes of next level
8920 delete queue2 ;
8921 queue2 = queue3 ;
8922 queue3 = new GridQueue2( ) ;
8923
8924
8925 // Next, start priority queue iteration
8926 while ( ! queue->isEmpty() )
8927 {
8928 // Retrieve the node with the highest score
8929 queue->remove( gp, score ) ;
8930 ox = gp->x ;
8931 oy = gp->y ;
8932 oz = gp->z ;
8933 delete gp ;
8934 // printf("%d\n", score);
8935 // score = -score ;
8936
8937 // Ignore the node
8938 // if it has been processed before
8939 // or it has an updated score
8940 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
8941 {
8942 continue ;
8943 }
8944
8945 /* Commented for debugging
8946
8947 // Remove this simple node
8948 setDataAt( ox, oy, oz, -1 ) ;
8949 numSimple ++ ;
8950 // printf("Highest score: %d\n", score) ;
8951 */
8952
8953 /* Added for debugging */
8954 // Check simple
8955 #ifndef NOISE_DIS_SHEET
8956 // if ( hasFeatureFace( ox, oy, oz ) )
8957 if ( (! isSimple( ox, oy, oz )) || isSheetEnd( ox, oy, oz ) )
8958 // if ( hasIsolatedFace(ox,oy,oz) && (! isNoiseSheetEnd(ox,oy,oz)))
8959 #else
8960 // if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz, noisevol ) )
8961 if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz, noisevol ) || isHelixEnd( ox, oy, oz, noisevol ))
8962 // if ( isBertrandEndPoint( ox, oy, oz ) )
8963 #endif
8964 {
8965 // Complex, set to next layer
8966 setDataAt( ox, oy, oz, curwid + 1 ) ;
8967 queue4->prepend( ox, oy, oz ) ;
8968 numComplex ++ ;
8969
8970 }
8971 else
8972 {
8973 setDataAt( ox, oy, oz, -1 ) ;
8974 numSimple ++ ;
8975
8976 }
8977 /* Adding ends */
8978
8979 // Move its neighboring unvisited node to queue2
8980 for ( int m = 0 ; m < 6 ; m ++ )
8981 {
8982 int nx = ox + neighbor6[m][0] ;
8983 int ny = oy + neighbor6[m][1] ;
8984 int nz = oz + neighbor6[m][2] ;
8985 if ( getDataAt( nx, ny, nz ) == 0 )
8986 {
8987 // setDataAt( nx, ny, nz, curwid + 1 ) ;
8988 queue2->prepend( nx, ny, nz ) ;
8989 }
8990 }
8991
8992 /* Commented for debugging
8993
8994 // Find complex nodes in its 3x3 neighborhood
8995 // move them to queue2
8996 for ( i = -1 ; i < 2 ; i ++ )
8997 for ( j = -1 ; j < 2 ; j ++ )
8998 for ( k = -1 ; k < 2 ; k ++ )
8999 {
9000 int nx = ox + i ;
9001 int ny = oy + j ;
9002 int nz = oz + k ;
9003
9004 // Check simple
9005 if ( getDataAt( nx, ny, nz ) == curwid &&
9006 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
9007 #ifndef NOISE_DIS_SHEET
9008 ( isSheetEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
9009 #else
9010 ( isSheetEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
9011 #endif
9012
9013 {
9014 // Complex, set to next layer
9015 setDataAt( nx, ny, nz, curwid + 1 ) ;
9016 queue2->prepend( nx, ny, nz ) ;
9017 numComplex ++ ;
9018 }
9019 }
9020 */
9021
9022 // Update scores for nodes in its 5x5 neighborhood
9023 // insert them back into priority queue
9024
9025 for ( i = -2 ; i < 3 ;i ++ )
9026 for ( j = -2 ; j < 3 ; j ++ )
9027 for ( k = -2 ; k < 3 ; k ++ )
9028 {
9029 int nx = ox + i ;
9030 int ny = oy + j ;
9031 int nz = oz + k ;
9032
9033 if ( getDataAt( nx, ny, nz ) == curwid )
9034 {
9035 // Compute score
9036 score = getNumPotComplex( nx, ny, nz ) ;
9037
9038 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
9039 {
9040 // printf("Update\n") ;
9041 scrvol->setDataAt( nx, ny, nz, score ) ;
9042 // Push to queue
9043 gp = new gridPoint ;
9044 gp->x = nx ;
9045 gp->y = ny ;
9046 gp->z = nz ;
9047 // queue->add( gp, -score ) ;
9048 queue->add( gp, score ) ;
9049 }
9050 }
9051 }
9052
9053
9054 }
9055
9056 #ifdef VERBOSE
9057 printf("%d complex, %d simple\n", numComplex, numSimple) ;
9058 #endif
9059
9060 if ( numSimple == 0 )
9061 {
9062 break ;
9063 }
9064 }
9065
9066 // Finally, clean up
9067 #ifdef VERBOSE
9068 printf("Thresholding the volume to 0/1...\n") ;
9069 #endif
9070 threshold( 0, 0, 1 ) ;
9071
9072 delete scrvol;
9073 delete queue;
9074 delete queue2;
9075 delete queue3;
9076 delete queue4;
9077
9078 }
wustl_mm::SkeletonMaker::Volume::isSheetEnd
int isSheetEnd(int ox, int oy, int oz, Volume *nvol)
Definition: volume.cpp:1821

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

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

◆ threshold() [1/4]

void Volume::threshold ( double  thr)

  • Bertrand's parallel 6-connected surface thinning *‍/* Palagyi's parallel surface thinning *‍/
  • Debug *‍/ Normalize to a given range

Definition at line 9515 of file volume.cpp.

9516 {
9517 threshold( thr, 0, 1, 0, true) ;
9518 }

References threshold().

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

◆ threshold() [2/4]

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

Definition at line 9520 of file volume.cpp.

9521 {
9522 threshold( thr, out, in, out, true) ;
9523 }

References threshold().

◆ threshold() [3/4]

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

Definition at line 9525 of file volume.cpp.

9526 {
9527 threshold(thr, out, in, boundary, true);
9528 }

References threshold().

◆ threshold() [4/4]

void Volume::threshold ( double  thr,
int  out,
int  in,
int  boundary,
bool  markBoundary 
)

Definition at line 9530 of file volume.cpp.

9531 {
9532 float val;
9533 for ( int i = 0 ; i < getSizeX() ; i ++ )
9534 for ( int j = 0 ; j < getSizeY() ; j ++ )
9535 for ( int k = 0 ; k < getSizeZ() ; k ++ )
9536 {
9537 val = (float)getDataAt(i, j, k);
9538 if(markBoundary) {
9539 if ( i > 1 && i < getSizeX() - 2 && j > 1 && j < getSizeY() - 2 && k > 1 && k < getSizeZ() - 2 ) {
9540 if(val < thr) {
9541 setDataAt(i, j, k, out);
9542 } else {
9543 setDataAt(i, j, k, in);
9544 }
9545 }
9546 else
9547 {
9548 setDataAt(i, j, k, boundary);
9549 }
9550 } else {
9551 if(val < thr) {
9552 setDataAt(i, j, k, out);
9553 } else {
9554 setDataAt(i, j, k, in);
9555 }
9556 }
9557 }
9558 }

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

Member Data Documentation

◆ volData

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(), getVolumeData(), pad(), setDataAt(), setOrigin(), setSpacing(), Volume(), and ~Volume().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.2