#include <volume.h>
Collaboration diagram for wustl_mm::SkeletonMaker::Volume:

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) |
| * Next, smooth */ | |
| 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) |
| * | |
| int | hasCompleteHelix (int ox, int oy, int oz, Volume *fvol) |
| int | isHelixEnd (int ox, int oy, int oz, Volume *nvol) |
| int | isHelixEnd (int ox, int oy, int oz) |
| int | isSheetEnd (int ox, int oy, int oz, Volume *nvol) |
| int | isFeatureFace (int ox, int oy, int oz) |
| int | isSheetEnd (int ox, int oy, int oz) |
| int | isSimple (int ox, int oy, int oz) |
| int | isPiercable (int ox, int oy, int oz) |
| int | getNumPotComplex (int ox, int oy, int oz) |
| int | getNumPotComplex2 (int ox, int oy, int oz) |
| int | components6 (int vox[3][3][3]) |
| * | |
| int | components26 (int vox[3][3][3]) |
| int | countExt (double vox[3][3][3]) |
| int | countInt (double vox[3][3][3]) |
| int | countIntEuler (int ox, int oy, int oz) |
| void | curveSkeleton (Volume *grayvol, float lowthr, float highthr, Volume *svol) |
| * Thin the current volume while preserving voxels with values > highthr or <= lowthr in grayvol | |
| void | curveSkeleton (float thr, Volume *svol) |
| void | curveSkeleton2D (float thr, Volume *svol) |
| void | skeleton (float thr, int off) |
| void | skeleton (float thr, Volume *svol, Volume *hvol) |
| * Adding ends *// Move its neighboring unvisited node to queue2 | |
| void | erodeHelix () |
| void | erodeHelix (int disthr) |
| int | erodeSheet () |
| int | erodeSheet (int disthr) |
| void | surfaceSkeletonPres (float thr, Volume *preserve) |
| * Commented for debugging/ Find complex nodes in its 3x3 neighborhood/ move them to queue2 | |
| void | threshold (double thr) |
| Normalize to a given range. | |
| void | threshold (double thr, int out, int in) |
| void | threshold (double thr, int out, int in, int boundary) |
| void | threshold (double thr, int out, int in, int boundary, bool markBoundary) |
| VolumeData * | getVolumeData () |
Private Attributes | |
| VolumeData * | volData |
|
|
Definition at line 9 of file volume.cpp. References volData. Referenced by curveSkeleton(), curveSkeleton2D(), erodeHelix(), erodeSheet(), markCellFace(), skeleton(), and surfaceSkeletonPres(). 00010 {
00011 this->volData = new VolumeData(em);
00012 }
|
|
||||||||||||||||
|
Definition at line 13 of file volume.cpp. 00013 {
00014 volData = new VolumeData(x, y, z);
00015 }
|
|
||||||||||||||||||||
|
Definition at line 17 of file volume.cpp. 00017 {
00018 volData = new VolumeData(x, y, z, val);
00019 }
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 21 of file volume.cpp. References getVolumeData(), volData, x, and y. 00021 {
00022 volData = new VolumeData(x, y, z, offx, offy, offz, vol->getVolumeData());
00023 }
|
|
|
Definition at line 26 of file volume.cpp. 00027 {
00028 delete volData;
00029 }
|
|
|
Definition at line 2767 of file volume.cpp. References nx, ny, tot, x, and y. Referenced by countExt(). 02768 {
02769 // Stupid flood fill
02770 int tot = 0 ;
02771 int queue[27][3] ;
02772 int vis[3][3][3] ;
02773 int head = 0, tail = 1 ;
02774 int i, j, k ;
02775 for ( i = 0 ; i < 3 ; i ++ )
02776 for ( j = 0 ; j < 3 ; j ++ )
02777 for ( k = 0 ; k < 3 ; k ++ )
02778 {
02779 vis[i][j][k] = 0 ;
02780 if ( vox[i][j][k] )
02781 {
02782 if ( tot == 0 )
02783 {
02784 queue[0][0] = i ;
02785 queue[0][1] = j ;
02786 queue[0][2] = k ;
02787 vis[i][j][k] = 1 ;
02788 }
02789 tot ++ ;
02790 }
02791 }
02792 if ( tot == 0 )
02793 {
02794 return 0 ;
02795 }
02796
02797 int ct = 1 ;
02798 while ( head != tail )
02799 {
02800 int x = queue[head][0] ;
02801 int y = queue[head][1] ;
02802 int z = queue[head][2] ;
02803 head ++ ;
02804
02805 for ( i = -1 ; i < 2 ; i ++ )
02806 for ( j = -1 ; j < 2 ; j ++ )
02807 for ( k = -1 ; k < 2 ; k ++ )
02808 {
02809 int nx = x + i ;
02810 int ny = y + j ;
02811 int nz = z + k ;
02812 if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
02813 {
02814 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
02815 {
02816 queue[tail][0] = nx ;
02817 queue[tail][1] = ny ;
02818 queue[tail][2] = nz ;
02819 tail ++ ;
02820 vis[nx][ny][nz] = 1 ;
02821 ct ++ ;
02822 }
02823 }
02824 }
02825 }
02826
02827 if ( ct == tot )
02828 {
02829 return 1 ;
02830 }
02831 else
02832 {
02833 return 2 ;
02834 }
02835
02836 }
|
|
|
*
Definition at line 2698 of file volume.cpp. References wustl_mm::SkeletonMaker::neighbor6, nx, ny, tot, x, and y. Referenced by countInt(), and countIntEuler(). 02699 {
02700 // Stupid flood fill
02701 int tot = 0 ;
02702 int queue[27][3] ;
02703 int vis[3][3][3] ;
02704 int head = 0, tail = 1 ;
02705 int i, j, k ;
02706 for ( i = 0 ; i < 3 ; i ++ )
02707 for ( j = 0 ; j < 3 ; j ++ )
02708 for ( k = 0 ; k < 3 ; k ++ )
02709 {
02710 vis[i][j][k] = 0 ;
02711 if ( vox[i][j][k] )
02712 {
02713 if ( tot == 0 )
02714 {
02715 queue[0][0] = i ;
02716 queue[0][1] = j ;
02717 queue[0][2] = k ;
02718 vis[i][j][k] = 1 ;
02719 }
02720 tot ++ ;
02721 }
02722 }
02723 if ( tot == 0 )
02724 {
02725 return 0 ;
02726 }
02727 // printf("total: %d\n", tot) ;
02728
02729 int ct = 1 ;
02730 while ( head != tail )
02731 {
02732 int x = queue[head][0] ;
02733 int y = queue[head][1] ;
02734 int z = queue[head][2] ;
02735 head ++ ;
02736
02737 for ( i = 0 ; i < 6 ; i ++ )
02738 {
02739 int nx = x + neighbor6[i][0] ;
02740 int ny = y + neighbor6[i][1] ;
02741 int nz = z + neighbor6[i][2] ;
02742 if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
02743 {
02744 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
02745 {
02746 queue[tail][0] = nx ;
02747 queue[tail][1] = ny ;
02748 queue[tail][2] = nz ;
02749 tail ++ ;
02750 vis[nx][ny][nz] = 1 ;
02751 ct ++ ;
02752 }
02753 }
02754 }
02755 }
02756
02757 if ( ct == tot )
02758 {
02759 return 1 ;
02760 }
02761 else
02762 {
02763 return 2 ;
02764 }
02765
02766 }
|
|
|
Definition at line 2838 of file volume.cpp. References components26(). Referenced by isPiercable(), and isSimple(). 02839 {
02840 int tvox[3][3][3] ;
02841
02842 for ( int i = 0 ; i < 3 ; i ++ )
02843 for ( int j = 0 ; j < 3 ; j ++ )
02844 for ( int k = 0 ; k < 3 ; k ++ )
02845 {
02846 if ( vox[i][j][k] < 0 )
02847 {
02848 tvox[i][j][k] = 1 ;
02849 }
02850 else
02851 {
02852 tvox[i][j][k] = 0 ;
02853 }
02854 }
02855
02856 return components26( tvox ) ;
02857 }
|
|
|
Definition at line 2859 of file volume.cpp. References components6(), wustl_mm::SkeletonMaker::neighbor6, wustl_mm::SkeletonMaker::neighbor64, nnz, nx, and ny. Referenced by isPiercable(), and isSimple(). 02860 {
02861 int i, j, k ;
02862 int tvox[3][3][3] ;
02863
02864 for ( i = 0 ; i < 3 ; i ++ )
02865 for ( j = 0 ; j < 3 ; j ++ )
02866 for ( k = 0 ; k < 3 ; k ++ )
02867 {
02868 tvox[i][j][k] = 0 ;
02869 }
02870
02871 for ( i = 0 ; i < 6 ; i ++ )
02872 {
02873 int nx = 1 + neighbor6[i][0] ;
02874 int ny = 1 + neighbor6[i][1] ;
02875 int nz = 1 + neighbor6[i][2] ;
02876 if ( vox[ nx ][ ny ][ nz ] >= 0 )
02877 {
02878 tvox[ nx ][ ny ][ nz ] = 1 ;
02879 for ( j = 0 ; j < 4 ; j ++ )
02880 {
02881 int nnx = neighbor64[i][j][0] + nx ;
02882 int nny = neighbor64[i][j][1] + ny ;
02883 int nnz = neighbor64[i][j][2] + nz ;
02884 if ( vox[ nnx ][ nny ][ nnz ] >= 0 )
02885 {
02886 tvox[ nnx ][ nny ][ nnz ] = 1 ;
02887 }
02888 }
02889 }
02890 }
02891
02892 return components6( tvox ) ;
02893 }
|
|
||||||||||||||||
|
Definition at line 2895 of file volume.cpp. References components6(), getDataAt(), wustl_mm::SkeletonMaker::neighbor6, wustl_mm::SkeletonMaker::neighbor64, nnz, nx, and ny. Referenced by hasCompleteSheet(). 02896 {
02897 int nv = 0 , ne = 0, nc = 0 ;
02898
02899 int i, j, k ;
02900 int tvox[3][3][3] ;
02901 double vox[3][3][3] ;
02902
02903 for ( i = 0 ; i < 3 ; i ++ )
02904 for ( j = 0 ; j < 3 ; j ++ )
02905 for ( k = 0 ; k < 3 ; k ++ )
02906 {
02907 vox[i][j][k] = getDataAt( ox - 1 + i, oy - 1 + j, oz - 1 + k ) ;
02908 tvox[i][j][k] = 0 ;
02909 }
02910
02911 for ( i = 0 ; i < 6 ; i ++ )
02912 {
02913 int nx = 1 + neighbor6[i][0] ;
02914 int ny = 1 + neighbor6[i][1] ;
02915 int nz = 1 + neighbor6[i][2] ;
02916 if ( vox[nx][ny][nz] >= 0 )
02917 {
02918 tvox[ nx ][ ny ][ nz ] = 1 ;
02919
02920 nv ++ ;
02921
02922 for ( j = 0 ; j < 4 ; j ++ )
02923 {
02924 int nnx = neighbor64[i][j][0] + nx ;
02925 int nny = neighbor64[i][j][1] + ny ;
02926 int nnz = neighbor64[i][j][2] + nz ;
02927 if ( vox[nnx][nny][nnz] >= 0 )
02928 {
02929 if ( tvox[nnx][nny][nnz] == 0 )
02930 {
02931 tvox[nnx][nny][nnz] = 1 ;
02932 nv ++ ;
02933 }
02934
02935 ne ++ ;
02936 }
02937 }
02938 }
02939 }
02940
02941 nc = components6( tvox ) ;
02942
02943 return ( nc - ( nv - ne ) ) ;
02944 }
|
|
||||||||||||
|
Definition at line 4268 of file volume.cpp. References EMAN::PriorityQueue< ValueT, KeyT >::add(), flag, getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isSimple(), MAX_ERODE, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z. 04269 {
04270 int i, j, k ;
04271 // First, threshold the volume
04272 #ifdef VERBOSE
04273 printf("Thresholding the volume to -1/0...\n") ;
04274 #endif
04275 threshold( thr, -1, 0 ) ;
04276
04277 // Next, apply convergent erosion
04278 // by preserving: complex nodes, curve end-points, and sheet points
04279
04280 // Next, initialize the linked queue
04281 #ifdef VERBOSE
04282 printf("Initializing queue...\n") ;
04283 #endif
04284 GridQueue2* queue2 = new GridQueue2( ) ;
04285 GridQueue2* queue3 = new GridQueue2( ) ;
04286 GridQueue2* queue4 = new GridQueue2( ) ;
04287 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
04288
04289 for ( i = 0 ; i < getSizeX() ; i ++ )
04290 for ( j = 0 ; j < getSizeY() ; j ++ )
04291 for ( k = 0 ; k < getSizeZ() ; k ++ )
04292 {
04293 if ( getDataAt( i, j, k ) >= 0 )
04294 {
04295 if ( svol->getDataAt(i,j,k) > 0 )
04296 {
04297 setDataAt( i, j, k, MAX_ERODE ) ;
04298 }
04299 else
04300 {
04301 for ( int m = 0 ; m < 6 ; m ++ )
04302 {
04303 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
04304 {
04305 // setDataAt( i, j, k, 1 ) ;
04306 queue2->prepend( i, j, k ) ;
04307 break ;
04308 }
04309 }
04310 }
04311 }
04312 }
04313
04314 int wid = MAX_ERODE ;
04315 #ifdef VERBOSE
04316 printf("Total %d nodes\n", queue2->getNumElements() ) ;
04317 printf("Start erosion to %d...\n", wid) ;
04318 #endif
04319
04320
04321 // Perform erosion
04322 gridQueueEle* ele ;
04323 gridPoint* gp ;
04324 int ox, oy, oz ;
04325 int score ;
04326 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
04327 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
04328 {
04329 scrvol->setDataAt( i, -1 ) ;
04330 }
04331
04332 #ifdef NOISE_DIS_HELIX
04333 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
04334 #endif
04335
04336 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
04337 {
04338 // At the start of each iteration,
04339 // queue2 holds all the nodes for this layer
04340 // queue3 and queue are empty
04341
04342 int numComplex = 0, numSimple = 0 ;
04343 #ifdef VERBOSE
04344 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
04345 #endif
04346
04347 /*
04348 We first need to assign curwid + 1 to every node in this layer
04349 */
04350 queue2->reset() ;
04351 ele = queue2->getNext() ;
04352 while ( ele != NULL )
04353 {
04354 ox = ele->x ;
04355 oy = ele->y ;
04356 oz = ele->z ;
04357
04358 if ( getDataAt(ox,oy,oz) == curwid )
04359 {
04360 ele = queue2->remove() ;
04361 }
04362 else
04363 {
04364 setDataAt(ox,oy,oz, curwid) ;
04365 ele = queue2->getNext() ;
04366 }
04367 }
04368 queue4->reset() ;
04369 ele = queue4->getNext() ;
04370 while ( ele != NULL )
04371 {
04372 ox = ele->x ;
04373 oy = ele->y ;
04374 oz = ele->z ;
04375
04376 queue2->prepend(ox,oy,oz) ;
04377 ele = queue4->remove() ;
04378 }
04379
04380 // Now queue2 holds all the nodes for this layer
04381
04382 #ifdef NOISE_DIS_HELIX
04383 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
04384 queue2->reset() ;
04385
04386 // First run
04387 int flag = 0 ;
04388 while ( ( ele = queue2->getNext() ) != NULL )
04389 {
04390 ox = ele->x ;
04391 oy = ele->y ;
04392 oz = ele->z ;
04393 if ( NOISE_DIS_HELIX <= 1 )
04394 {
04395 noisevol->setDataAt( ox, oy, oz, 0 ) ;
04396 }
04397 else
04398 {
04399 flag = 0 ;
04400 for ( int m = 0 ; m < 6 ; m ++ )
04401 {
04402 int nx = ox + neighbor6[m][0] ;
04403 int ny = oy + neighbor6[m][1] ;
04404 int nz = oz + neighbor6[m][2] ;
04405 if ( getDataAt( nx, ny, nz ) == 0 )
04406 {
04407 noisevol->setDataAt( ox, oy, oz, 1 ) ;
04408 flag = 1 ;
04409 break ;
04410 }
04411 }
04412 if ( ! flag )
04413 {
04414 noisevol->setDataAt( ox, oy, oz, 0 ) ;
04415 }
04416 }
04417 }
04418
04419 int cur, visited ;
04420 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
04421 {
04422 queue2->reset() ;
04423 int count = 0 ;
04424 visited = 0 ;
04425
04426 while ( ( ele = queue2->getNext() ) != NULL )
04427 {
04428 ox = ele->x ;
04429 oy = ele->y ;
04430 oz = ele->z ;
04431
04432 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
04433 {
04434 visited ++ ;
04435 continue ;
04436 }
04437
04438 flag = 0 ;
04439 for ( int m = 0 ; m < 6 ; m ++ )
04440 {
04441 int nx = ox + neighbor6[m][0] ;
04442 int ny = oy + neighbor6[m][1] ;
04443 int nz = oz + neighbor6[m][2] ;
04444 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
04445 {
04446 noisevol->setDataAt( ox, oy, oz, 1 ) ;
04447 visited ++ ;
04448 count ++ ;
04449 break ;
04450 }
04451 }
04452 }
04453
04454 if ( count == 0 )
04455 {
04456 break ;
04457 }
04458 }
04459 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
04460
04461
04462 #endif
04463 /* Commented out for debugging
04464
04465 // First,
04466 // check for complex nodes in queue2
04467 // move them from queue2 to queue3
04468 queue2->reset() ;
04469 ele = queue2->getNext() ;
04470 while ( ele != NULL )
04471 {
04472 ox = ele->x ;
04473 oy = ele->y ;
04474 oz = ele->z ;
04475
04476 // Check simple
04477 #ifndef NOISE_DIS_HELIX
04478 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04479 #else
04480 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
04481 #endif
04482 {
04483 // Complex, set to next layer
04484 setDataAt( ox, oy, oz, curwid + 1 ) ;
04485 queue3->prepend( ox, oy, oz ) ;
04486 ele = queue2->remove() ;
04487
04488 numComplex ++ ;
04489 }
04490 else
04491 {
04492 ele = queue2->getNext() ;
04493 }
04494 }
04495 */
04496
04497 // Next,
04498 // Compute score for each node left in queue2
04499 // move them into priority queue
04500 queue2->reset() ;
04501 ele = queue2->getNext() ;
04502 while ( ele != NULL )
04503 {
04504 ox = ele->x ;
04505 oy = ele->y ;
04506 oz = ele->z ;
04507
04508 // Compute score
04509 score = getNumPotComplex2( ox, oy, oz ) ;
04510 scrvol->setDataAt( ox, oy, oz, score ) ;
04511
04512 // Push to queue
04513 gp = new gridPoint ;
04514 gp->x = ox ;
04515 gp->y = oy ;
04516 gp->z = oz ;
04517 // queue->add( gp, -score ) ;
04518 queue->add( gp, score ) ;
04519
04520 ele = queue2->remove() ;
04521 }
04522
04523 // Rename queue3 to be queue2,
04524 // Clear queue3
04525 // From now on, queue2 holds nodes of next level
04526 delete queue2 ;
04527 queue2 = queue3 ;
04528 queue3 = new GridQueue2( ) ;
04529
04530 // Next, start priority queue iteration
04531 while ( ! queue->isEmpty() )
04532 {
04533 // Retrieve the node with the highest score
04534 queue->remove( gp, score ) ;
04535 ox = gp->x ;
04536 oy = gp->y ;
04537 oz = gp->z ;
04538 delete gp ;
04539 // score = -score ;
04540
04541 // Ignore the node
04542 // if it has been processed before
04543 // or it has an updated score
04544 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
04545 {
04546 continue ;
04547 }
04548
04549 /* Commented out for debugging
04550
04551 // Remove this simple node
04552 setDataAt( ox, oy, oz, -1 ) ;
04553 numSimple ++ ;
04554 // printf("Highest score: %d\n", score) ;
04555 */
04556
04557 /* Added for debugging */
04558 // Check simple
04559 #ifndef NOISE_DIS_HELIX
04560 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
04561 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04562 #else
04563 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04564 #endif
04565 {
04566 // Complex, set to next layer
04567 setDataAt( ox, oy, oz, curwid + 1 ) ;
04568 queue4->prepend( ox, oy, oz ) ;
04569 numComplex ++ ;
04570 }
04571 else
04572 {
04573 setDataAt( ox, oy, oz, -1 ) ;
04574 numSimple ++ ;
04575 }
04576 /* Adding ends */
04577
04578 // Move its neighboring unvisited node to queue2
04579 for ( int m = 0 ; m < 6 ; m ++ )
04580 {
04581 int nx = ox + neighbor6[m][0] ;
04582 int ny = oy + neighbor6[m][1] ;
04583 int nz = oz + neighbor6[m][2] ;
04584 if ( getDataAt( nx, ny, nz ) == 0 )
04585 {
04586 // setDataAt( nx, ny, nz, curwid + 1 ) ;
04587 queue2->prepend( nx, ny, nz ) ;
04588 }
04589 }
04590
04591 /* Commented out for debugging
04592
04593 // Find complex nodes in its 3x3 neighborhood
04594 // move them to queue2
04595 for ( i = -1 ; i < 2 ; i ++ )
04596 for ( j = -1 ; j < 2 ; j ++ )
04597 for ( k = -1 ; k < 2 ; k ++ )
04598 {
04599 int nx = ox + i ;
04600 int ny = oy + j ;
04601 int nz = oz + k ;
04602
04603 // Check simple
04604 if ( getDataAt( nx, ny, nz ) == curwid &&
04605 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
04606 #ifndef NOISE_DIS_HELIX
04607 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
04608 #else
04609 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
04610 #endif
04611
04612 {
04613 // Complex, set to next layer
04614 setDataAt( nx, ny, nz, curwid + 1 ) ;
04615 queue2->prepend( nx, ny, nz ) ;
04616 numComplex ++ ;
04617 }
04618 }
04619 */
04620
04621 // Update scores for nodes in its 5x5 neighborhood
04622 // insert them back into priority queue
04623
04624 for ( i = -2 ; i < 3 ;i ++ )
04625 for ( j = -2 ; j < 3 ; j ++ )
04626 for ( k = -2 ; k < 3 ; k ++ )
04627 {
04628 int nx = ox + i ;
04629 int ny = oy + j ;
04630 int nz = oz + k ;
04631
04632 if ( getDataAt( nx, ny, nz ) == curwid )
04633 {
04634 // Compute score
04635 score = getNumPotComplex2( nx, ny, nz ) ;
04636
04637 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
04638 {
04639 // printf("Update\n") ;
04640 scrvol->setDataAt( nx, ny, nz, score ) ;
04641 // Push to queue
04642 gp = new gridPoint ;
04643 gp->x = nx ;
04644 gp->y = ny ;
04645 gp->z = nz ;
04646 // queue->add( gp, -score ) ;
04647 queue->add( gp, score ) ;
04648 }
04649 }
04650 }
04651
04652
04653 }
04654
04655 #ifdef VERBOSE
04656 printf("%d complex, %d simple\n", numComplex, numSimple) ;
04657 #endif
04658
04659 if ( numSimple == 0 )
04660 {
04661 break ;
04662 }
04663 }
04664
04665 // Finally, clean up
04666 delete scrvol;
04667 delete queue;
04668 delete queue2;
04669 delete queue3;
04670 delete queue4;
04671 #ifdef VERBOSE
04672 printf("Thresholding the volume to 0/1...\n") ;
04673 #endif
04674 threshold( 0, 0, 1 ) ;
04675 }
|
|
||||||||||||||||||||
|
* Thin the current volume while preserving voxels with values > highthr or <= lowthr in grayvol
Definition at line 3812 of file volume.cpp. References EMAN::PriorityQueue< ValueT, KeyT >::add(), flag, getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isPiercable(), isSimple(), MAX_ERODE, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), v, Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z. Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning(). 03813 {
03814 int i, j, k ;
03815 // First, threshold the volume
03816 #ifdef VERBOSE
03817 printf("Thresholding the volume to -1/0...\n") ;
03818 #endif
03819 threshold( 0.5f, -1, 0 ) ;
03820
03821 // Next, apply convergent erosion
03822 // by preserving: complex nodes, curve end-points, and sheet points
03823
03824 // Next, initialize the linked queue
03825 #ifdef VERBOSE
03826 printf("Initializing queue...\n") ;
03827 #endif
03828 GridQueue2* queue2 = new GridQueue2( ) ;
03829 GridQueue2* queue3 = new GridQueue2( ) ;
03830 GridQueue2* queue4 = new GridQueue2( ) ;
03831 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
03832
03833 for ( i = 0 ; i < getSizeX() ; i ++ )
03834 for ( j = 0 ; j < getSizeY() ; j ++ )
03835 for ( k = 0 ; k < getSizeZ() ; k ++ )
03836 {
03837 if ( getDataAt( i, j, k ) >= 0 )
03838 {
03839 float v = (float)grayvol->getDataAt(i,j,k) ;
03840 if ( v <= lowthr || v > highthr || svol->getDataAt(i,j,k) > 0 )
03841 {
03842 setDataAt( i, j, k, MAX_ERODE ) ;
03843 }
03844 else
03845 {
03846 for ( int m = 0 ; m < 6 ; m ++ )
03847 {
03848 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
03849 {
03850 // setDataAt( i, j, k, 1 ) ;
03851 queue2->prepend( i, j, k ) ;
03852 break ;
03853 }
03854 }
03855 }
03856 }
03857 }
03858 int wid = MAX_ERODE ;
03859 #ifdef VERBOSE
03860 printf("Total %d nodes\n", queue2->getNumElements() ) ;
03861 printf("Start erosion to %d...\n", wid) ;
03862 #endif
03863
03864
03865 // Perform erosion
03866 gridQueueEle* ele ;
03867 gridPoint* gp ;
03868 int ox, oy, oz ;
03869 int score ;
03870 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
03871 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
03872 {
03873 scrvol->setDataAt( i, -1 ) ;
03874 }
03875
03876 #ifdef NOISE_DIS_HELIX
03877 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
03878 #endif
03879
03880 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
03881 {
03882 // At the start of each iteration,
03883 // queue2 holds all the nodes for this layer
03884 // queue3 and queue are empty
03885
03886 int numComplex = 0, numSimple = 0 ;
03887 #ifdef VERBOSE
03888 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
03889 #endif
03890
03891 /*
03892 We first need to assign curwid + 1 to every node in this layer
03893 */
03894 queue2->reset() ;
03895 ele = queue2->getNext() ;
03896 while ( ele != NULL )
03897 {
03898 ox = ele->x ;
03899 oy = ele->y ;
03900 oz = ele->z ;
03901
03902 if ( getDataAt(ox,oy,oz) == curwid )
03903 {
03904 ele = queue2->remove() ;
03905 }
03906 else
03907 {
03908 setDataAt(ox,oy,oz, curwid) ;
03909 ele = queue2->getNext() ;
03910 }
03911 }
03912 queue4->reset() ;
03913 ele = queue4->getNext() ;
03914 while ( ele != NULL )
03915 {
03916 ox = ele->x ;
03917 oy = ele->y ;
03918 oz = ele->z ;
03919
03920 queue2->prepend(ox,oy,oz) ;
03921 ele = queue4->remove() ;
03922 }
03923
03924 // Now queue2 holds all the nodes for this layer
03925
03926 #ifdef NOISE_DIS_HELIX
03927 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
03928 queue2->reset() ;
03929
03930 // First run
03931 int flag = 0 ;
03932 while ( ( ele = queue2->getNext() ) != NULL )
03933 {
03934 ox = ele->x ;
03935 oy = ele->y ;
03936 oz = ele->z ;
03937 if ( NOISE_DIS_HELIX <= 1 )
03938 {
03939 noisevol->setDataAt( ox, oy, oz, 0 ) ;
03940 }
03941 else
03942 {
03943 flag = 0 ;
03944 for ( int m = 0 ; m < 6 ; m ++ )
03945 {
03946 int nx = ox + neighbor6[m][0] ;
03947 int ny = oy + neighbor6[m][1] ;
03948 int nz = oz + neighbor6[m][2] ;
03949 if ( getDataAt( nx, ny, nz ) == 0 )
03950 {
03951 noisevol->setDataAt( ox, oy, oz, 1 ) ;
03952 flag = 1 ;
03953 break ;
03954 }
03955 }
03956 if ( ! flag )
03957 {
03958 noisevol->setDataAt( ox, oy, oz, 0 ) ;
03959 }
03960 }
03961 }
03962
03963 int cur, visited ;
03964 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
03965 {
03966 queue2->reset() ;
03967 int count = 0 ;
03968 visited = 0 ;
03969
03970 while ( ( ele = queue2->getNext() ) != NULL )
03971 {
03972 ox = ele->x ;
03973 oy = ele->y ;
03974 oz = ele->z ;
03975
03976 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
03977 {
03978 visited ++ ;
03979 continue ;
03980 }
03981
03982 flag = 0 ;
03983 for ( int m = 0 ; m < 6 ; m ++ )
03984 {
03985 int nx = ox + neighbor6[m][0] ;
03986 int ny = oy + neighbor6[m][1] ;
03987 int nz = oz + neighbor6[m][2] ;
03988 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
03989 {
03990 noisevol->setDataAt( ox, oy, oz, 1 ) ;
03991 visited ++ ;
03992 count ++ ;
03993 break ;
03994 }
03995 }
03996 }
03997
03998 if ( count == 0 )
03999 {
04000 break ;
04001 }
04002 }
04003 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
04004
04005
04006 #endif
04007 /* Commented out for debugging
04008
04009 // First,
04010 // check for complex nodes in queue2
04011 // move them from queue2 to queue3
04012 queue2->reset() ;
04013 ele = queue2->getNext() ;
04014 while ( ele != NULL )
04015 {
04016 ox = ele->x ;
04017 oy = ele->y ;
04018 oz = ele->z ;
04019
04020 // Check simple
04021 #ifndef NOISE_DIS_HELIX
04022 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04023 #else
04024 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
04025 #endif
04026 {
04027 // Complex, set to next layer
04028 setDataAt( ox, oy, oz, curwid + 1 ) ;
04029 queue3->prepend( ox, oy, oz ) ;
04030 ele = queue2->remove() ;
04031
04032 numComplex ++ ;
04033 }
04034 else
04035 {
04036 ele = queue2->getNext() ;
04037 }
04038 }
04039 */
04040
04041 // Next,
04042 // Compute score for each node left in queue2
04043 // move them into priority queue
04044 queue2->reset() ;
04045 ele = queue2->getNext() ;
04046 while ( ele != NULL )
04047 {
04048 ox = ele->x ;
04049 oy = ele->y ;
04050 oz = ele->z ;
04051
04052 // Compute score
04053 score = getNumPotComplex2( ox, oy, oz ) ;
04054 scrvol->setDataAt( ox, oy, oz, score ) ;
04055
04056 // Push to queue
04057 gp = new gridPoint ;
04058 gp->x = ox ;
04059 gp->y = oy ;
04060 gp->z = oz ;
04061 // queue->add( gp, -score ) ;
04062 queue->add( gp, score ) ;
04063
04064 ele = queue2->remove() ;
04065 }
04066
04067 // Rename queue3 to be queue2,
04068 // Clear queue3
04069 // From now on, queue2 holds nodes of next level
04070 delete queue2 ;
04071 queue2 = queue3 ;
04072 queue3 = new GridQueue2( ) ;
04073
04074 // Next, start priority queue iteration
04075 while ( ! queue->isEmpty() )
04076 {
04077 // Retrieve the node with the highest score
04078 queue->remove( gp, score ) ;
04079 ox = gp->x ;
04080 oy = gp->y ;
04081 oz = gp->z ;
04082 delete gp ;
04083 // score = -score ;
04084
04085 // Ignore the node
04086 // if it has been processed before
04087 // or it has an updated score
04088 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
04089 {
04090 continue ;
04091 }
04092
04093 /* Commented out for debugging
04094
04095 // Remove this simple node
04096 setDataAt( ox, oy, oz, -1 ) ;
04097 numSimple ++ ;
04098 // printf("Highest score: %d\n", score) ;
04099 */
04100
04101 /* Added for debugging */
04102 // Check simple
04103 #ifndef NOISE_DIS_HELIX
04104 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
04105 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04106 #else
04107 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04108 #endif
04109 {
04110 // Complex, set to next layer
04111 setDataAt( ox, oy, oz, curwid + 1 ) ;
04112 queue4->prepend( ox, oy, oz ) ;
04113 numComplex ++ ;
04114 }
04115 else
04116 {
04117 setDataAt( ox, oy, oz, -1 ) ;
04118 numSimple ++ ;
04119
04120
04121 /* Adding ends */
04122 // Move its neighboring unvisited node to queue2
04123 for ( int m = 0 ; m < 6 ; m ++ )
04124 {
04125 int nx = ox + neighbor6[m][0] ;
04126 int ny = oy + neighbor6[m][1] ;
04127 int nz = oz + neighbor6[m][2] ;
04128 if ( getDataAt( nx, ny, nz ) == 0 )
04129 {
04130 // setDataAt( nx, ny, nz, curwid + 1 ) ;
04131 queue2->prepend( nx, ny, nz ) ;
04132 }
04133 }
04134
04135 }
04136
04137
04138 /* Commented out for debugging
04139
04140 // Find complex nodes in its 3x3 neighborhood
04141 // move them to queue2
04142 for ( i = -1 ; i < 2 ; i ++ )
04143 for ( j = -1 ; j < 2 ; j ++ )
04144 for ( k = -1 ; k < 2 ; k ++ )
04145 {
04146 int nx = ox + i ;
04147 int ny = oy + j ;
04148 int nz = oz + k ;
04149
04150 // Check simple
04151 if ( getDataAt( nx, ny, nz ) == curwid &&
04152 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
04153 #ifndef NOISE_DIS_HELIX
04154 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
04155 #else
04156 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
04157 #endif
04158
04159 {
04160 // Complex, set to next layer
04161 setDataAt( nx, ny, nz, curwid + 1 ) ;
04162 queue2->prepend( nx, ny, nz ) ;
04163 numComplex ++ ;
04164 }
04165 }
04166 */
04167
04168 // Update scores for nodes in its 5x5 neighborhood
04169 // insert them back into priority queue
04170
04171 for ( i = -2 ; i < 3 ;i ++ )
04172 for ( j = -2 ; j < 3 ; j ++ )
04173 for ( k = -2 ; k < 3 ; k ++ )
04174 {
04175 int nx = ox + i ;
04176 int ny = oy + j ;
04177 int nz = oz + k ;
04178
04179 if ( getDataAt( nx, ny, nz ) == curwid )
04180 {
04181 // Compute score
04182 score = getNumPotComplex2( nx, ny, nz ) ;
04183
04184 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
04185 {
04186 // printf("Update\n") ;
04187 scrvol->setDataAt( nx, ny, nz, score ) ;
04188 // Push to queue
04189 gp = new gridPoint ;
04190 gp->x = nx ;
04191 gp->y = ny ;
04192 gp->z = nz ;
04193 // queue->add( gp, -score ) ;
04194 queue->add( gp, score ) ;
04195 }
04196 }
04197 }
04198
04199
04200 }
04201
04202 #ifdef VERBOSE
04203 printf("%d complex, %d simple\n", numComplex, numSimple) ;
04204 #endif
04205
04206 if ( numSimple == 0 )
04207 {
04208 if ( queue2->getNumElements() > 0 )
04209 {
04210 printf("*************************wierd here*************************\n");
04211 }
04212 break ;
04213 }
04214 }
04215
04216 // Remove all internal voxels (contained in manifold surfaces)
04217 queue2->reset() ;
04218 queue4->reset() ;
04219 ele = queue4->getNext() ;
04220 while ( ele != NULL )
04221 {
04222 ox = ele->x ;
04223 oy = ele->y ;
04224 oz = ele->z ;
04225
04226 if ( isPiercable(ox,oy,oz) == 1 ) // hasCompleteSheet( ox, oy, oz ) == 1 ) //
04227 {
04228 queue2->prepend(ox,oy,oz) ;
04229 // setDataAt( ox, oy, oz, -1 ) ;
04230 }
04231 ele = queue4->remove() ;
04232 }
04233
04234 for ( i = 0 ; i < getSizeX() ; i ++ )
04235 for ( j = 0 ; j < getSizeY() ; j ++ )
04236 for ( k = 0 ; k < getSizeZ() ; k ++ )
04237 {
04238 if ( getDataAt( i, j, k ) == 0 && isPiercable(i,j,k) ) //hasCompleteSheet(i,j,k) == 1) //
04239 {
04240 queue2->prepend( i, j, k ) ;
04241 }
04242 }
04243 queue2->reset() ;
04244 ele = queue2->getNext() ;
04245 while ( ele != NULL )
04246 {
04247 ox = ele->x ;
04248 oy = ele->y ;
04249 oz = ele->z ;
04250 setDataAt( ox, oy, oz, -1 ) ;
04251 ele = queue2->remove() ;
04252 }
04253
04254
04255 // Finally, clean up
04256 delete scrvol;
04257 delete queue;
04258 delete queue2;
04259 delete queue3;
04260 delete queue4;
04261 #ifdef VERBOSE
04262 printf("Thresholding the volume to 0/1...\n") ;
04263 #endif
04264 threshold( 0, 0, 1 ) ;
04265 }
|
|
||||||||||||
|
Definition at line 4678 of file volume.cpp. References EMAN::PriorityQueue< ValueT, KeyT >::add(), flag, getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isSimple(), MAX_ERODE, wustl_mm::SkeletonMaker::neighbor4, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z. Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning(). 04679 {
04680 int i, j, k ;
04681 // First, threshold the volume
04682 #ifdef VERBOSE
04683 printf("Thresholding the volume to -1/0...\n") ;
04684 #endif
04685 threshold( thr, -1, 0 ) ;
04686
04687 // Next, apply convergent erosion
04688 // by preserving: complex nodes, curve end-points, and sheet points
04689
04690 // Next, initialize the linked queue
04691 #ifdef VERBOSE
04692 printf("Initializing queue...\n") ;
04693 #endif
04694 GridQueue2* queue2 = new GridQueue2( ) ;
04695 GridQueue2* queue3 = new GridQueue2( ) ;
04696 GridQueue2* queue4 = new GridQueue2( ) ;
04697 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
04698
04699 for ( i = 0 ; i < getSizeX() ; i ++ )
04700 for ( j = 0 ; j < getSizeY() ; j ++ )
04701 for ( k = 0 ; k < getSizeZ() ; k ++ )
04702 {
04703 if ( getDataAt( i, j, k ) >= 0 )
04704 {
04705 if ( svol->getDataAt(i,j,k) > 0 )
04706 {
04707 setDataAt( i, j, k, MAX_ERODE ) ;
04708 }
04709 else
04710 {
04711 for ( int m = 0 ; m < 4 ; m ++ )
04712 {
04713 if ( getDataAt( i + neighbor4[m][0], j + neighbor4[m][1], k ) < 0 )
04714 {
04715 // setDataAt( i, j, k, 1 ) ;
04716 queue2->prepend( i, j, k ) ;
04717 break ;
04718 }
04719 }
04720 }
04721 }
04722 }
04723 int wid = MAX_ERODE ;
04724 #ifdef VERBOSE
04725 printf("Total %d nodes\n", queue2->getNumElements() ) ;
04726 printf("Start erosion to %d...\n", wid) ;
04727 #endif
04728
04729
04730 // Perform erosion
04731 gridQueueEle* ele ;
04732 gridPoint* gp ;
04733 int ox, oy, oz ;
04734 int score ;
04735 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
04736 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
04737 {
04738 scrvol->setDataAt( i, -1 ) ;
04739 }
04740
04741 #ifdef NOISE_DIS_HELIX
04742 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
04743 #endif
04744
04745 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
04746 {
04747 // At the start of each iteration,
04748 // queue2 holds all the nodes for this layer
04749 // queue3 and queue are empty
04750
04751 int numComplex = 0, numSimple = 0 ;
04752 #ifdef VERBOSE
04753 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
04754 #endif
04755
04756 /*
04757 We first need to assign curwid + 1 to every node in this layer
04758 */
04759 queue2->reset() ;
04760 ele = queue2->getNext() ;
04761 while ( ele != NULL )
04762 {
04763 ox = ele->x ;
04764 oy = ele->y ;
04765 oz = ele->z ;
04766
04767 if ( getDataAt(ox,oy,oz) == curwid )
04768 {
04769 ele = queue2->remove() ;
04770 }
04771 else
04772 {
04773 setDataAt(ox,oy,oz, curwid) ;
04774 ele = queue2->getNext() ;
04775 }
04776 }
04777 queue4->reset() ;
04778 ele = queue4->getNext() ;
04779 while ( ele != NULL )
04780 {
04781 ox = ele->x ;
04782 oy = ele->y ;
04783 oz = ele->z ;
04784
04785 queue2->prepend(ox,oy,oz) ;
04786 ele = queue4->remove() ;
04787 }
04788
04789 // Now queue2 holds all the nodes for this layer
04790
04791 #ifdef NOISE_DIS_HELIX
04792 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
04793 queue2->reset() ;
04794
04795 // First run
04796 int flag = 0 ;
04797 while ( ( ele = queue2->getNext() ) != NULL )
04798 {
04799 ox = ele->x ;
04800 oy = ele->y ;
04801 oz = ele->z ;
04802 if ( NOISE_DIS_HELIX <= 1 )
04803 {
04804 noisevol->setDataAt( ox, oy, oz, 0 ) ;
04805 }
04806 else
04807 {
04808 flag = 0 ;
04809 for ( int m = 0 ; m < 6 ; m ++ )
04810 {
04811 int nx = ox + neighbor6[m][0] ;
04812 int ny = oy + neighbor6[m][1] ;
04813 int nz = oz + neighbor6[m][2] ;
04814 if ( getDataAt( nx, ny, nz ) == 0 )
04815 {
04816 noisevol->setDataAt( ox, oy, oz, 1 ) ;
04817 flag = 1 ;
04818 break ;
04819 }
04820 }
04821 if ( ! flag )
04822 {
04823 noisevol->setDataAt( ox, oy, oz, 0 ) ;
04824 }
04825 }
04826 }
04827
04828 int cur, visited ;
04829 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
04830 {
04831 queue2->reset() ;
04832 int count = 0 ;
04833 visited = 0 ;
04834
04835 while ( ( ele = queue2->getNext() ) != NULL )
04836 {
04837 ox = ele->x ;
04838 oy = ele->y ;
04839 oz = ele->z ;
04840
04841 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
04842 {
04843 visited ++ ;
04844 continue ;
04845 }
04846
04847 flag = 0 ;
04848 for ( int m = 0 ; m < 6 ; m ++ )
04849 {
04850 int nx = ox + neighbor6[m][0] ;
04851 int ny = oy + neighbor6[m][1] ;
04852 int nz = oz + neighbor6[m][2] ;
04853 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
04854 {
04855 noisevol->setDataAt( ox, oy, oz, 1 ) ;
04856 visited ++ ;
04857 count ++ ;
04858 break ;
04859 }
04860 }
04861 }
04862
04863 if ( count == 0 )
04864 {
04865 break ;
04866 }
04867 }
04868 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
04869
04870
04871 #endif
04872 /* Commented out for debugging
04873
04874 // First,
04875 // check for complex nodes in queue2
04876 // move them from queue2 to queue3
04877 queue2->reset() ;
04878 ele = queue2->getNext() ;
04879 while ( ele != NULL )
04880 {
04881 ox = ele->x ;
04882 oy = ele->y ;
04883 oz = ele->z ;
04884
04885 // Check simple
04886 #ifndef NOISE_DIS_HELIX
04887 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04888 #else
04889 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
04890 #endif
04891 {
04892 // Complex, set to next layer
04893 setDataAt( ox, oy, oz, curwid + 1 ) ;
04894 queue3->prepend( ox, oy, oz ) ;
04895 ele = queue2->remove() ;
04896
04897 numComplex ++ ;
04898 }
04899 else
04900 {
04901 ele = queue2->getNext() ;
04902 }
04903 }
04904 */
04905
04906 // Next,
04907 // Compute score for each node left in queue2
04908 // move them into priority queue
04909 queue2->reset() ;
04910 ele = queue2->getNext() ;
04911 while ( ele != NULL )
04912 {
04913 ox = ele->x ;
04914 oy = ele->y ;
04915 oz = ele->z ;
04916
04917 // Compute score
04918 score = getNumPotComplex2( ox, oy, oz ) ;
04919 //score = getNumNeighbor6( ox, oy, oz ) ;
04920 scrvol->setDataAt( ox, oy, oz, score ) ;
04921
04922 // Push to queue
04923 gp = new gridPoint ;
04924 gp->x = ox ;
04925 gp->y = oy ;
04926 gp->z = oz ;
04927 // queue->add( gp, -score ) ;
04928 queue->add( gp, score ) ;
04929
04930 ele = queue2->remove() ;
04931 }
04932
04933 // Rename queue3 to be queue2,
04934 // Clear queue3
04935 // From now on, queue2 holds nodes of next level
04936 delete queue2 ;
04937 queue2 = queue3 ;
04938 queue3 = new GridQueue2( ) ;
04939
04940 // Next, start priority queue iteration
04941 while ( ! queue->isEmpty() )
04942 {
04943 // Retrieve the node with the highest score
04944 queue->remove( gp, score ) ;
04945 ox = gp->x ;
04946 oy = gp->y ;
04947 oz = gp->z ;
04948 delete gp ;
04949 // score = -score ;
04950
04951 // Ignore the node
04952 // if it has been processed before
04953 // or it has an updated score
04954 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
04955 {
04956 continue ;
04957 }
04958
04959 /* Commented out for debugging
04960
04961 // Remove this simple node
04962 setDataAt( ox, oy, oz, -1 ) ;
04963 numSimple ++ ;
04964 // printf("Highest score: %d\n", score) ;
04965 */
04966
04967 /* Added for debugging */
04968 // Check simple
04969 #ifndef NOISE_DIS_HELIX
04970 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
04971 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04972 #else
04973 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
04974 #endif
04975 {
04976 // Complex, set to next layer
04977 setDataAt( ox, oy, oz, curwid + 1 ) ;
04978 queue4->prepend( ox, oy, oz ) ;
04979 numComplex ++ ;
04980 }
04981 else
04982 {
04983 setDataAt( ox, oy, oz, -1 ) ;
04984 numSimple ++ ;
04985 }
04986 /* Adding ends */
04987
04988 // Move its neighboring unvisited node to queue2
04989 for ( int m = 0 ; m < 4 ; m ++ )
04990 {
04991 int nx = ox + neighbor4[m][0] ;
04992 int ny = oy + neighbor4[m][1] ;
04993 int nz = oz ;
04994 if ( getDataAt( nx, ny, nz ) == 0 )
04995 {
04996 // setDataAt( nx, ny, nz, curwid + 1 ) ;
04997 queue2->prepend( nx, ny, nz ) ;
04998 }
04999 }
05000
05001 /* Commented out for debugging
05002
05003 // Find complex nodes in its 3x3 neighborhood
05004 // move them to queue2
05005 for ( i = -1 ; i < 2 ; i ++ )
05006 for ( j = -1 ; j < 2 ; j ++ )
05007 for ( k = -1 ; k < 2 ; k ++ )
05008 {
05009 int nx = ox + i ;
05010 int ny = oy + j ;
05011 int nz = oz + k ;
05012
05013 // Check simple
05014 if ( getDataAt( nx, ny, nz ) == curwid &&
05015 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
05016 #ifndef NOISE_DIS_HELIX
05017 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
05018 #else
05019 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
05020 #endif
05021
05022 {
05023 // Complex, set to next layer
05024 setDataAt( nx, ny, nz, curwid + 1 ) ;
05025 queue2->prepend( nx, ny, nz ) ;
05026 numComplex ++ ;
05027 }
05028 }
05029 */
05030
05031 // Update scores for nodes in its 5x5 neighborhood
05032 // insert them back into priority queue
05033
05034 for ( i = -2 ; i < 3 ;i ++ )
05035 for ( j = -2 ; j < 3 ; j ++ )
05036 for ( k = -2 ; k < 3 ; k ++ )
05037 {
05038 int nx = ox + i ;
05039 int ny = oy + j ;
05040 int nz = oz + k ;
05041
05042 if ( getDataAt( nx, ny, nz ) == curwid )
05043 {
05044 // Compute score
05045 score = getNumPotComplex2( nx, ny, nz ) ;
05046 //score = getNumNeighbor6( nx, ny, nz ) ;
05047
05048 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
05049 {
05050 // printf("Update\n") ;
05051 scrvol->setDataAt( nx, ny, nz, score ) ;
05052 // Push to queue
05053 gp = new gridPoint ;
05054 gp->x = nx ;
05055 gp->y = ny ;
05056 gp->z = nz ;
05057 // queue->add( gp, -score ) ;
05058 queue->add( gp, score ) ;
05059 }
05060 }
05061 }
05062
05063
05064 }
05065
05066 #ifdef VERBOSE
05067 printf("%d complex, %d simple\n", numComplex, numSimple) ;
05068 #endif
05069
05070 if ( numSimple == 0 )
05071 {
05072 break ;
05073 }
05074 }
05075
05076 // Finally, clean up
05077 #ifdef VERBOSE
05078 printf("Thresholding the volume to 0/1...\n") ;
05079 #endif
05080 threshold( 0, 0, 1 ) ;
05081 delete scrvol;
05082 delete queue;
05083 delete queue2;
05084 delete queue3;
05085 delete queue4;
05086 }
|
|
|
Definition at line 5911 of file volume.cpp. References getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getSizeX(), getSizeY(), getSizeZ(), hasCompleteHelix(), wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridQueueEle::y, and wustl_mm::SkeletonMaker::gridQueueEle::z. 05912 {
05913 int i, j, k ;
05914 // First, threshold the volume
05915 //printf("Thresholding the volume to -1/0...\n") ;
05916 threshold( 0.1f, -1, 0 ) ;
05917
05918 /* Debug: remove faces */
05919 //Volume* facevol = markFaceEdge() ;
05920 /* End debugging */
05921
05922 // Next, initialize the linked queue
05923 //printf("Initializing queue...\n") ;
05924 Volume * fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
05925 GridQueue2* queue2 = new GridQueue2( ) ;
05926 GridQueue2* queue3 = new GridQueue2( ) ;
05927 GridQueue2** queues = new GridQueue2* [ 10000 ] ;
05928
05929 for ( i = 0 ; i < getSizeX() ; i ++ )
05930 for ( j = 0 ; j < getSizeY() ; j ++ )
05931 for ( k = 0 ; k < getSizeZ() ; k ++ )
05932 {
05933 if ( getDataAt( i, j, k ) >= 0 )
05934 {
05935 if ( ! hasCompleteHelix( i, j, k ) )
05936 // if ( ! hasCompleteHelix( i, j, k, facevol ) )
05937 {
05938 queue2->prepend( i, j, k ) ;
05939 fvol->setDataAt( i, j, k, -1 ) ;
05940 }
05941 }
05942 }
05943 //printf("Total %d nodes\n", queue2->getNumElements() ) ;
05944
05945 // Start erosion
05946 gridQueueEle* ele ;
05947 int dis = -1 ;
05948 while ( queue2->getNumElements() > 0 )
05949 {
05950 // First, set distance
05951 dis -- ;
05952 queues[ - dis ] = new GridQueue2( ) ;
05953 //printf("Distance transform to %d...", dis) ;
05954 queue2->reset() ;
05955 while( ( ele = queue2->getNext() ) != NULL )
05956 {
05957 setDataAt( ele->x, ele->y, ele->z, dis ) ;
05958 queues[ -dis ]->prepend( ele->x, ele->y, ele->z ) ;
05959 }
05960 //printf("%d nodes\n", queues[-dis]->getNumElements()) ;
05961
05962 // Next, find next layer
05963 queue2->reset() ;
05964 ele = queue2->getNext() ;
05965 while ( ele != NULL )
05966 {
05967 for ( int m = 0 ; m < 6 ; m ++ )
05968 {
05969 int nx = ele->x + neighbor6[m][0] ;
05970 int ny = ele->y + neighbor6[m][1] ;
05971 int nz = ele->z + neighbor6[m][2] ;
05972 if ( getDataAt( nx, ny, nz ) == 0 )
05973 {
05974 fvol->setDataAt( nx, ny, nz, dis ) ;
05975
05976 if ( ! hasCompleteHelix( nx, ny, nz ) )
05977 // if ( ! hasCompleteHelix( nx, ny, nz, facevol ) )
05978 {
05979 setDataAt( nx, ny, nz, 1 ) ;
05980 queue3->prepend( nx, ny, nz ) ;
05981 }
05982 }
05983 }
05984
05985 ele = queue2->remove() ;
05986 }
05987
05988 // Next, swap queues
05989 GridQueue2* temp = queue2 ;
05990 queue2 = queue3 ;
05991 queue3 = temp ;
05992 }
05993
05994 // Deal with closed rings
05995 dis -- ;
05996 queues[ - dis ] = new GridQueue2( ) ;
05997 #ifdef VERBOSE
05998 printf("Detecting closed rings %d...", dis) ;
05999 #endif
06000 int ftot = 0 ;
06001 for ( i = 0 ; i < getSizeX() ; i ++ )
06002 for ( j = 0 ; j < getSizeY() ; j ++ )
06003 for ( k = 0 ; k < getSizeZ() ; k ++ )
06004 {
06005 if ( getDataAt( i, j, k ) == 0 )
06006 {
06007 setDataAt( i, j, k, dis ) ;
06008 queues[ -dis ]->prepend( i, j, k ) ;
06009 /*
06010 int fval = (int) fvol->getDataAt( i, j, k ) ;
06011 if ( fval == 0)
06012 {
06013 // queues[ -dis ]->prepend( i, j, k ) ;
06014 }
06015 else
06016 {
06017 setDataAt( i, j, k, fval - 1 ) ;
06018 // queues[ -fval + 1 ]->prepend( i, j, k ) ;
06019 }
06020 */
06021 ftot ++ ;
06022 }
06023 }
06024 #ifdef VERBOSE
06025 printf("%d nodes\n", ftot) ;
06026 #endif
06027
06028
06029 // return ;
06030
06031 /* Find local minimum: to help determine erosion level
06032 int cts[ 64 ] ;
06033 for ( i = 0 ; i <= - dis ; i ++ )
06034 {
06035 cts[ i ] = 0 ;
06036 }
06037 for ( i = 0 ; i < getSizeX() ; i ++ )
06038 for ( j = 0 ; j < getSizeY() ; j ++ )
06039 for ( k = 0 ; k < getSizeZ() ; k ++ )
06040 {
06041 double val = getDataAt( i, j, k ) ;
06042 if ( val < -1 )
06043 {
06044 int min = 1 ;
06045 for ( int m = 0 ; m < 6 ; m ++ )
06046 {
06047 int nx = i + neighbor6[m][0] ;
06048 int ny = j + neighbor6[m][1] ;
06049 int nz = k + neighbor6[m][2] ;
06050 if ( getDataAt( nx, ny, nz ) < val )
06051 {
06052 min = 0 ;
06053 break ;
06054 }
06055 }
06056
06057 if ( min )
06058 {
06059 cts[ (int) - val ] ++ ;
06060 }
06061 }
06062 }
06063
06064 for ( i = 2 ; i <= - dis ; i ++ )
06065 {
06066 printf("Local minima: %d with %d nodes.\n", -i, cts[ i ] ) ;
06067 }
06068 */
06069
06070 // Dilate back
06071 // Starting from nodes with distance - 2 - disthr
06072
06073 if ( disthr + 2 > - dis )
06074 {
06075 disthr = - dis - 2 ;
06076 }
06077 int d;
06078
06079 for ( d = - dis ; d > disthr + 1 ; d -- )
06080 {
06081 queues[ d ]->reset() ;
06082 while ( (ele = queues[ d ]->getNext() ) != NULL )
06083 {
06084 setDataAt( ele->x, ele->y, ele->z, d ) ;
06085 }
06086 }
06087
06088
06089 for ( d = disthr + 1 ; d >= 2 ; d -- )
06090 {
06091 //delete queue3 ;
06092 //queue3 = new GridQueue2( ) ;
06093 queues[ d ]->reset() ;
06094 ele = queues[ d ]->getNext() ;
06095 while ( ele != NULL )
06096 {
06097 int dilatable = 0 ;
06098 for ( int m = 0 ; m < 6 ; m ++ )
06099 {
06100 int nx = ele->x + neighbor6[m][0] ;
06101 int ny = ele->y + neighbor6[m][1] ;
06102 int nz = ele->z + neighbor6[m][2] ;
06103 if ( getDataAt( nx, ny, nz ) == d + 1 )
06104 {
06105 dilatable = 1 ;
06106 break ;
06107 }
06108 }
06109
06110
06111 if ( ! dilatable )
06112 {
06113 /*
06114 setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
06115 */
06116
06117 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
06118 if ( d > 2 )
06119 {
06120 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
06121 }
06122 ele = queues[ d ]->remove() ;
06123 }
06124 else
06125 {
06126 setDataAt( ele->x, ele->y, ele->z, d ) ;
06127 ele = queues[ d ]->getNext() ;
06128 }
06129
06130 }
06131
06132 /* Debug: extract minimal set */
06133 while ( 1 )
06134 {
06135 int num = 0 ;
06136 queues[ d ]->reset() ;
06137 ele = queues[ d ]->getNext() ;
06138 while ( ele != NULL )
06139 {
06140 int critical = 0 ;
06141 setDataAt( ele->x, ele->y, ele->z, -1 ) ;
06142
06143 for ( i = -1 ; i < 2 ; i ++ )
06144 {
06145 for ( j = -1 ; j < 2 ; j ++ )
06146 {
06147 for ( k = -1 ; k < 2 ; k ++ )
06148 {
06149 if ( i != 0 && j != 0 && k != 0 )
06150 {
06151 continue ;
06152 }
06153 int nx = ele->x + i ;
06154 int ny = ele->y + j ;
06155 int nz = ele->z + k ;
06156 if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteHelix( nx,ny,nz ) ) //, facevol ) )
06157 {
06158 critical = 1 ;
06159 break ;
06160 }
06161 }
06162 if ( critical )
06163 {
06164 break ;
06165 }
06166 }
06167 if ( critical )
06168 {
06169 break ;
06170 }
06171 }
06172
06173 if ( critical )
06174 {
06175 setDataAt( ele->x, ele->y, ele->z, d ) ;
06176 ele = queues[ d ]->getNext() ;
06177 }
06178 else
06179 {
06180 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
06181 if ( d > 2 )
06182 {
06183 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
06184 }
06185 ele = queues[ d ]->remove() ;
06186 num ++ ;
06187 }
06188
06189 }
06190
06191 #ifdef VERBOSE
06192 printf("Non-minimal: %d\n", num) ;
06193 #endif
06194
06195 if ( num == 0 )
06196 {
06197 break ;
06198 }
06199 }
06200
06201
06202 /* End of debugging */
06203
06204 /*
06205 queue3->reset() ;
06206 ele = queue3->getNext() ;
06207 while ( ele != NULL )
06208 {
06209 setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
06210 ele = queue3->remove() ;
06211 }
06212 */
06213 }
06214
06215 // Finally, threshold the volume
06216 #ifdef VERBOSE
06217 //printf("Thresholding the volume to 0/1...\n") ;
06218 #endif
06219 //threshold( -1, 1, 0, 0 ) ;
06220 threshold( 0, 0, 1 ) ;
06221 delete fvol ;
06222 delete queue2;
06223 delete queue3;
06224 for ( d = -dis ; d >= 2 ; d -- ) {
06225 delete queues[d];
06226 }
06227 delete [] queues;
06228
06229 }
|
|
|
Definition at line 5905 of file volume.cpp. Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PruneCurves(). 05906 {
05907 erodeHelix( 3 ) ;
05908 }
|
|
|
Definition at line 6240 of file volume.cpp. References flag, getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getSizeX(), getSizeY(), getSizeZ(), hasCompleteSheet(), markCellFace(), nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), wustl_mm::SkeletonMaker::sheetNeighbor, threshold(), Volume(), wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridQueueEle::y, and wustl_mm::SkeletonMaker::gridQueueEle::z. 06241 {
06242 int i, j, k ;
06243 // First, threshold the volume
06244 //printf("Thresholding the volume to -1/0...\n") ;
06245 threshold( 0.1f, -1, 0 ) ;
06246
06247 /* Debug: remove cells */
06248 Volume* facevol = markCellFace() ;
06249 /* End debugging */
06250
06251 // Next, initialize the linked queue
06252 //printf("Initializing queue...\n") ;
06253 Volume * fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
06254 GridQueue2* queue2 = new GridQueue2( ) ;
06255 GridQueue2* queue3 = new GridQueue2( ) ;
06256 GridQueue2** queues = new GridQueue2* [ 10000 ] ;
06257
06258 for ( i = 0 ; i < getSizeX() ; i ++ )
06259 for ( j = 0 ; j < getSizeY() ; j ++ )
06260 for ( k = 0 ; k < getSizeZ() ; k ++ )
06261 {
06262 if ( getDataAt( i, j, k ) >= 0 )
06263 {
06264 if ( ! hasCompleteSheet( i, j, k ) )
06265 //if ( ! hasCompleteSheet( i, j, k, facevol ) )
06266 {
06267 queue2->prepend( i, j, k ) ;
06268 fvol->setDataAt( i, j, k, -1 ) ;
06269 }
06270 }
06271 }
06272 #ifdef VERBOSE
06273 printf("Total %d nodes\n", queue2->getNumElements() ) ;
06274 #endif
06275
06276 // Start erosion
06277 gridQueueEle* ele ;
06278 int dis = -1 ;
06279 while ( queue2->getNumElements() > 0 )
06280 {
06281 // First, set distance
06282 dis -- ;
06283 queues[ - dis ] = new GridQueue2( ) ;
06284 //printf("Distance transform to %d...", dis) ;
06285 queue2->reset() ;
06286 while( ( ele = queue2->getNext() ) != NULL )
06287 {
06288 setDataAt( ele->x, ele->y, ele->z, dis ) ;
06289 queues[ -dis ]->prepend( ele->x, ele->y, ele->z ) ;
06290 }
06291 //printf("%d nodes\n", queues[-dis]->getNumElements()) ;
06292
06293 // Next, find next layer
06294 queue2->reset() ;
06295 ele = queue2->getNext() ;
06296 while ( ele != NULL )
06297 {
06298 // for ( int m = 0 ; m < 6 ; m ++ )
06299 for ( int mx = -1 ; mx < 2 ; mx ++ )
06300 for ( int my = -1 ; my < 2 ; my ++ )
06301 for ( int mz = -1 ; mz < 2 ; mz ++ )
06302 {
06303 if ( mx != 0 && my != 0 && mz != 0 )
06304 {
06305 continue ;
06306 }
06307 //int nx = ele->x + neighbor6[m][0] ;
06308 //int ny = ele->y + neighbor6[m][1] ;
06309 //int nz = ele->z + neighbor6[m][2] ;
06310 int nx = ele->x + mx ;
06311 int ny = ele->y + my ;
06312 int nz = ele->z + mz ;
06313
06314 if ( getDataAt( nx, ny, nz ) == 0 )
06315 {
06316 fvol->setDataAt( nx, ny, nz, dis ) ;
06317
06318 if ( ! hasCompleteSheet( nx, ny, nz ) )
06319 // if ( ! hasCompleteSheet( nx, ny, nz, facevol ) )
06320 {
06321 setDataAt( nx, ny, nz, 1 ) ;
06322 queue3->prepend( nx, ny, nz ) ;
06323 }
06324 }
06325 }
06326
06327 ele = queue2->remove() ;
06328 }
06329
06330 // Next, swap queues
06331 GridQueue2* temp = queue2 ;
06332 queue2 = queue3 ;
06333 queue3 = temp ;
06334 }
06335
06336 /* Deal with closed rings */
06337
06338 dis -- ;
06339 queues[ - dis ] = new GridQueue2( ) ;
06340 #ifdef VERBOSE
06341 printf("Detecting closed rings %d...", dis) ;
06342 #endif
06343 int ftot = 0 ;
06344 for ( i = 0 ; i < getSizeX() ; i ++ )
06345 for ( j = 0 ; j < getSizeY() ; j ++ )
06346 for ( k = 0 ; k < getSizeZ() ; k ++ )
06347 {
06348 if ( getDataAt( i, j, k ) == 0 )
06349 {
06350 /*
06351 int fval = (int) fvol->getDataAt( i, j, k ) ;
06352 if ( fval == 0)
06353 {
06354 setDataAt( i, j, k, dis - 2 ) ;
06355 // queues[ -dis ]->prepend( i, j, k ) ;
06356 }
06357 else
06358 {
06359 setDataAt( i, j, k, fval - 1 ) ;
06360 queues[ -fval + 1 ]->prepend( i, j, k ) ;
06361 }
06362 */
06363 setDataAt( i, j, k, dis ) ;
06364 queues[ -dis ]->prepend( i, j, k ) ;
06365
06366 ftot ++ ;
06367 }
06368 }
06369 #ifdef VERBOSE
06370 printf("%d nodes\n", ftot) ;
06371 #endif
06372
06373
06374 /* Find local minimum: to help determine erosion level
06375 int cts[ 64 ] ;
06376 for ( i = 0 ; i <= - dis ; i ++ )
06377 {
06378 cts[ i ] = 0 ;
06379 }
06380 for ( i = 0 ; i < getSizeX() ; i ++ )
06381 for ( j = 0 ; j < getSizeY() ; j ++ )
06382 for ( k = 0 ; k < getSizeZ() ; k ++ )
06383 {
06384 double val = getDataAt( i, j, k ) ;
06385 if ( val < -1 )
06386 {
06387 int min = 1 ;
06388 for ( int m = 0 ; m < 6 ; m ++ )
06389 {
06390 int nx = i + neighbor6[m][0] ;
06391 int ny = j + neighbor6[m][1] ;
06392 int nz = k + neighbor6[m][2] ;
06393 if ( getDataAt( nx, ny, nz ) < val )
06394 {
06395 min = 0 ;
06396 break ;
06397 }
06398 }
06399
06400 if ( min )
06401 {
06402 cts[ (int) - val ] ++ ;
06403 }
06404 }
06405 }
06406
06407 for ( i = 2 ; i <= - dis ; i ++ )
06408 {
06409 printf("Local minima: %d with %d nodes.\n", -i, cts[ i ] ) ;
06410 }
06411 */
06412
06413 // return ;
06414
06415 // Dilate back
06416 // Starting from nodes with distance - 2 - disthr
06417 int d ;
06418 if ( disthr + 2 > - dis )
06419 {
06420 disthr = - dis - 2 ;
06421
06422 }
06423 for ( d = - dis ; d > disthr + 1 ; d -- )
06424 {
06425 queues[ d ]->reset() ;
06426 while ( (ele = queues[ d ]->getNext() ) != NULL )
06427 {
06428 setDataAt( ele->x, ele->y, ele->z, d ) ;
06429 }
06430 }
06431
06432 for (d = disthr + 1 ; d >= 2 ; d -- )
06433 {
06434
06435 //delete queue3 ;
06436 //queue3 = new GridQueue2( ) ;
06437 queues[ d ]->reset() ;
06438 ele = queues[ d ]->getNext() ;
06439 while ( ele != NULL )
06440 {
06441 int dilatable = 0 ;
06442 // for ( int m = 0 ; m < 6 ; m ++ )
06443 /*
06444 for ( int mx = -1 ; mx < 2 ; mx ++ )
06445 for ( int my = -1 ; my < 2 ; my ++ )
06446 for ( int mz = -1 ; mz < 2 ; mz ++ )
06447 {
06448 if ( mx == 0 || my == 0 || mz == 0 )
06449 {
06450 int nx = ele->x + mx ; // neighbor6[m][0] ;
06451 int ny = ele->y + my ; // neighbor6[m][1] ;
06452 int nz = ele->z + mz ; // neighbor6[m][2] ;
06453 if ( getDataAt( nx, ny, nz ) == - d - 1 )
06454 {
06455 dilatable = 1 ;
06456 break ;
06457 }
06458 }
06459 }
06460 */
06461 for ( i = 0 ; i < 12 ; i ++ )
06462 {
06463 int flag = 1, flag2 = 0 ;
06464 for ( j = 0 ; j < 4 ; j ++ )
06465 {
06466 int nx = ele->x + sheetNeighbor[i][j][0] ;
06467 int ny = ele->y + sheetNeighbor[i][j][1] ;
06468 int nz = ele->z + sheetNeighbor[i][j][2] ;
06469
06470 double val = getDataAt( nx, ny, nz ) ;
06471
06472 if ( val > - d && val < 0)
06473 {
06474 flag = 0 ;
06475 break ;
06476 }
06477 if ( val == d + 1 )
06478 {
06479 flag2 ++ ;
06480 }
06481 }
06482
06483 if ( flag && flag2 )
06484 {
06485 dilatable = 1 ;
06486 break ;
06487 }
06488 }
06489
06490 if ( ! dilatable )
06491 {
06492 // setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
06493 // queue3->prepend( ele->x, ele->y, ele->z ) ;
06494
06495 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
06496 if ( d > 2 )
06497 {
06498 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
06499 }
06500 ele = queues[ d ]->remove() ;
06501 }
06502 else
06503 {
06504 setDataAt( ele->x, ele->y, ele->z, d ) ;
06505 ele = queues[ d ]->getNext() ;
06506 }
06507 }
06508
06509 /* Debug: extract minimal set */
06510 while ( 1 )
06511 {
06512 int num = 0 ;
06513 queues[ d ]->reset() ;
06514 ele = queues[ d ]->getNext() ;
06515 while ( ele != NULL )
06516 {
06517 int critical = 0 ;
06518 setDataAt( ele->x, ele->y, ele->z, -1 ) ;
06519
06520 for ( i = -1 ; i < 2 ; i ++ )
06521 {
06522 for ( j = -1 ; j < 2 ; j ++ )
06523 {
06524 for ( k = -1 ; k < 2 ; k ++ )
06525 {
06526 if ( i != 0 && j != 0 && k != 0 )
06527 {
06528 continue ;
06529 }
06530 int nx = ele->x + i ;
06531 int ny = ele->y + j ;
06532 int nz = ele->z + k ;
06533 // if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteSheet( nx,ny,nz, facevol ) )
06534 if ( getDataAt(nx,ny,nz) == d + 1 && !hasCompleteSheet( nx,ny,nz ) )
06535 {
06536 critical = 1 ;
06537 break ;
06538 }
06539 }
06540 if ( critical )
06541 {
06542 break ;
06543 }
06544 }
06545 if ( critical )
06546 {
06547 break ;
06548 }
06549 }
06550
06551 if ( critical )
06552 {
06553 setDataAt( ele->x, ele->y, ele->z, d ) ;
06554 ele = queues[ d ]->getNext() ;
06555 }
06556 else
06557 {
06558 setDataAt( ele->x, ele->y, ele->z, - d + 1 ) ;
06559 if ( d > 2 )
06560 {
06561 queues[ d - 1 ]->prepend( ele->x, ele->y, ele->z ) ;
06562 }
06563 ele = queues[ d ]->remove() ;
06564 num ++ ;
06565 }
06566
06567 }
06568 #ifdef VERBOSE
06569 printf("Non-minimal: %d\n", num) ;
06570 #endif
06571
06572 if ( num == 0 )
06573 {
06574 break ;
06575 }
06576 }
06577
06578
06579 /* End of debugging */
06580
06581 /*
06582 queue3->reset() ;
06583 ele = queue3->getNext() ;
06584 while ( ele != NULL )
06585 {
06586 setDataAt( ele->x, ele->y, ele->z, - 1 ) ;
06587 ele = queue3->remove() ;
06588 }
06589 */
06590 }
06591
06592
06593 // Finally, threshold the volume
06594 #ifdef VERBOSE
06595 //printf("Thresholding the volume to 0/1...\n") ;
06596 #endif
06597 //threshold( -1, 1, 0, 0 ) ;
06598 threshold( 0, 0, 1 ) ;
06599
06600 delete facevol ;
06601 delete fvol ;
06602 delete queue2;
06603 delete queue3;
06604 for (d = -dis ; d >= 2 ; d -- ) {
06605 delete queues[d];
06606 }
06607 delete [] queues;
06608
06609 return - dis - 1 ;
06610 }
|
|
|
Definition at line 6234 of file volume.cpp. Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PruneSurfaces(). 06235 {
06236 return erodeSheet( 3 ) ;
06237 }
|
|
|
Definition at line 31 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::get_emdata(), and getVolumeData(). Referenced by EMAN::BinarySkeletonizerProcessor::process(). 00032 {
00033 return this->getVolumeData()->get_emdata();
00034 }
|
|
|
Definition at line 65 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::GetDataAt(), and volData.
|
|
||||||||||||||||
|
||||||||||||||||
|
Definition at line 48 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::GetIndex(), volData, x, and y. Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::MarkSurfaces().
|
|
||||||||||||||||
|
Definition at line 661 of file volume.cpp. References getDataAt(), wustl_mm::SkeletonMaker::neighbor6, nx, and ny. Referenced by getNumPotComplex(), isFeatureFace(), and isHelixEnd(). 00661 {
00662 int rvalue = 0 ;
00663 for ( int i = 0 ; i < 6 ; i ++ ) {
00664 int nx = ox + neighbor6[i][0] ;
00665 int ny = oy + neighbor6[i][1] ;
00666 int nz = oz + neighbor6[i][2] ;
00667 if ( getDataAt( nx, ny, nz ) >= 0 ) {
00668 rvalue ++ ;
00669 }
00670 }
00671
00672 return rvalue ;
00673 }
|
|
||||||||||||||||
|
Definition at line 2548 of file volume.cpp. References getDataAt(), getNumNeighbor6(), wustl_mm::SkeletonMaker::neighbor6, nx, ny, and setDataAt(). Referenced by getNumPotComplex2(), and surfaceSkeletonPres(). 02549 {
02550 //return 0 ;
02551
02552 int i;
02553 double val = getDataAt( ox, oy, oz ) ;
02554 if ( val <= 0 )
02555 {
02556 // return 70 ;
02557 }
02558
02559 // return getNumNeighbor6( ox, oy, oz ) ;
02560
02561 // int v = ((getNumNeighbor6( ox, oy, oz ) & 255) << 24) ;
02562 //int v = 0 ;
02563
02564 int rvalue = 0, nx, ny, nz ;
02565 setDataAt( ox, oy, oz, -val ) ;
02566
02567 /*
02568 for ( i = -1 ; i < 2 ; i ++ )
02569 for ( j = -1 ; j < 2 ; j ++ )
02570 for ( k = -1 ; k < 2 ; k ++ )
02571 {
02572 nx = ox + i ;
02573 ny = oy + j ;
02574 nz = oz + k ;
02575 if ( getDataAt( nx, ny, nz ) == val )
02576 {
02577 if ( isSheetEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
02578 {
02579 rvalue ++ ;
02580 }
02581 }
02582 }
02583 */
02584
02585 for ( i = 0 ; i < 6 ; i ++ )
02586 {
02587 nx = ox + neighbor6[i][0] ;
02588 ny = oy + neighbor6[i][1] ;
02589 nz = oz + neighbor6[i][2] ;
02590
02591 if ( getDataAt(nx,ny,nz) >= 0 )
02592 {
02593 int num = getNumNeighbor6( nx, ny, nz ) ;
02594 if ( num > rvalue )
02595 {
02596 rvalue = num ;
02597 }
02598 }
02599 }
02600
02601
02602 setDataAt( ox, oy, oz, val ) ;
02603
02604 return rvalue + getNumNeighbor6( ox, oy, oz ) * 10 ;
02605 /*
02606 int v = (((rvalue + getNumNeighbor6( ox, oy, oz ) * 10) & 255) << 24) ;
02607 v |= ( ( ox & 255 ) << 16 ) ;
02608 v |= ( ( oy & 255 ) << 8 ) ;
02609 v |= ( ( oz & 255 ) ) ;
02610 return v ;
02611 */
02612
02613 }
|
|
||||||||||||||||
|
Definition at line 2615 of file volume.cpp. References getNumPotComplex(). Referenced by curveSkeleton(), curveSkeleton2D(), and skeleton(). 02616 {
02617 return getNumPotComplex( ox, oy, oz ) ;
02618
02619 //int i, j, k ;
02620 //double val = getDataAt( ox, oy, oz ) ;
02621 //if ( val <= 0 )
02622 //{
02623 // return 0 ;
02624 //}
02625
02626 //int rvalue = 0, nx, ny, nz ;
02627 //setDataAt( ox, oy, oz, -val ) ;
02628
02629 //for ( i = -1 ; i < 2 ; i ++ )
02630 // for ( j = -1 ; j < 2 ; j ++ )
02631 // for ( k = -1 ; k < 2 ; k ++ )
02632 // {
02633 // nx = ox + i ;
02634 // ny = oy + j ;
02635 // nz = oz + k ;
02636 // if ( getDataAt( nx, ny, nz ) == val )
02637 // {
02638 // if ( isHelixEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
02639 // {
02640 // rvalue ++ ;
02641 // }
02642 // }
02643 // }
02644
02645 //setDataAt( ox, oy, oz, val ) ;
02646
02647 //return rvalue + getNumNeighbor6( ox, oy, oz ) * 30 ;
02648 }
|
|
|
Definition at line 93 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::GetOriginX(), and volData. 00093 {
00094 return volData->GetOriginX();
00095 }
|
|
|
Definition at line 97 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::GetOriginY(), and volData. 00097 {
00098 return volData->GetOriginY();
00099 }
|
|
|
Definition at line 101 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::GetOriginZ(), and volData. 00101 {
00102 return volData->GetOriginZ();
00103 }
|
|
|
|
|
|
Definition at line 81 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::GetSpacingX(), and volData. 00081 {
00082 return volData->GetSpacingX();
00083 }
|
|
|
Definition at line 85 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::GetSpacingY(), and volData. 00085 {
00086 return volData->GetSpacingY();
00087 }
|
|
|
Definition at line 89 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::GetSpacingZ(), and volData. 00089 {
00090 return volData->GetSpacingZ();
00091 }
|
|
|
Definition at line 69 of file volume.cpp. Referenced by get_emdata(), EMAN::BinarySkeletonizerProcessor::process(), and Volume(). 00069 {
00070 return volData;
00071 }
|
|
||||||||||||||||
|
Definition at line 1016 of file volume.cpp. References getDataAt(). Referenced by hasCompleteSheet(), and markCellFace(). 01016 {
01017 for ( int i = 0 ; i < 2 ; i ++ )
01018 for ( int j = 0 ; j < 2 ; j ++ )
01019 for ( int k = 0 ; k < 2 ; k ++ ) {
01020 if ( getDataAt( ox + i, oy + j, oz + k ) < 0 ) {
01021 return 0 ;
01022 }
01023 }
01024 return 1 ;
01025 }
|
|
||||||||||||||||||||
|
Definition at line 1512 of file volume.cpp. References getDataAt(), wustl_mm::SkeletonMaker::neighbor6, nx, and ny. 01513 {
01514
01515 int i ;
01516 int c1 = 0;
01517 int nx, ny, nz ;
01518 int j ;
01519
01520 for ( i = 0 ; i < 6 ; i ++ )
01521 {
01522 nx = ox + neighbor6[i][0] ;
01523 ny = oy + neighbor6[i][1] ;
01524 nz = oz + neighbor6[i][2] ;
01525 if ( getDataAt( nx, ny, nz ) >= 0 )
01526 {
01527 if ( i % 2 == 0 )
01528 {
01529 nx = ox ;
01530 ny = oy ;
01531 nz = oz ;
01532 }
01533
01534 int val = (int)fvol->getDataAt( nx, ny, nz) ;
01535 if ( (( val >> ( 2 * ( i / 2 ) ) ) & 1) == 0 )
01536 {
01537 c1 ++ ;
01538 j = i ;
01539 }
01540 }
01541
01542 }
01543
01544 if ( c1 > 1 )
01545 {
01546 return 1 ;
01547 }
01548
01549 return 0 ;
01550
01551 /*
01552 ox = ox + neighbor6[j][0] ;
01553 oy = oy + neighbor6[j][1] ;
01554 oz = oz + neighbor6[j][2] ;
01555 c1 = 0 ;
01556 for ( i = 0 ; i < 6 ; i ++ )
01557 {
01558 nx = ox + neighbor6[i][0] ;
01559 ny = oy + neighbor6[i][1] ;
01560 nz = oz + neighbor6[i][2] ;
01561 if ( getDataAt( nx, ny, nz ) >= 0 )
01562 {
01563 c1 ++ ;
01564 }
01565
01566 }
01567
01568 if ( c1 > 1 )
01569 {
01570 return 0 ;
01571 }
01572 else
01573 {
01574 return 1 ;
01575 }
01576 */
01577 }
|
|
||||||||||||||||
|
*
Definition at line 1456 of file volume.cpp. References getDataAt(), wustl_mm::SkeletonMaker::neighbor6, nx, and ny. Referenced by erodeHelix(). 01457 {
01458 // Returns 1 if it has a complete helix
01459 int i ;
01460 int c1 = 0;
01461 int nx, ny, nz ;
01462 int j ;
01463
01464 for ( i = 0 ; i < 6 ; i ++ )
01465 {
01466 nx = ox + neighbor6[i][0] ;
01467 ny = oy + neighbor6[i][1] ;
01468 nz = oz + neighbor6[i][2] ;
01469 if ( getDataAt( nx, ny, nz ) >= 0 )
01470 {
01471 c1 ++ ;
01472 j = i ;
01473 }
01474
01475 }
01476
01477 if ( c1 > 1 ) // || c1 == 0 )
01478 {
01479 return 1 ;
01480 }
01481
01482 return 0 ;
01483
01484 /*
01485 ox = ox + neighbor6[j][0] ;
01486 oy = oy + neighbor6[j][1] ;
01487 oz = oz + neighbor6[j][2] ;
01488 c1 = 0 ;
01489 for ( i = 0 ; i < 6 ; i ++ )
01490 {
01491 nx = ox + neighbor6[i][0] ;
01492 ny = oy + neighbor6[i][1] ;
01493 nz = oz + neighbor6[i][2] ;
01494 if ( getDataAt( nx, ny, nz ) >= 0 )
01495 {
01496 c1 ++ ;
01497 }
01498
01499 }
01500
01501 if ( c1 > 1 )
01502 {
01503 return 0 ;
01504 }
01505 else
01506 {
01507 return 1 ;
01508 }
01509 */
01510 }
|
|
||||||||||||||||
|
Definition at line 1322 of file volume.cpp. References countIntEuler(). 01322 {
01323 // Returns 1 if it lies in the middle of a sheet
01324 int temp = countIntEuler( ox, oy, oz ) ;
01325 if ( temp > 0 )
01326 {
01327 return 1 ;
01328 }
01329 else
01330 {
01331 return 0 ;
01332 }
01333 }
|
|
||||||||||||||||||||
|
Definition at line 1174 of file volume.cpp. References wustl_mm::SkeletonMaker::edgeFaces, wustl_mm::SkeletonMaker::faceCells, wustl_mm::SkeletonMaker::faceEdges, getDataAt(), hasCell(), nx, ny, wustl_mm::SkeletonMaker::sheetNeighbor, and tot. Referenced by erodeSheet(), and isSheetEnd(). 01174 {
01175 int i, j, k ;
01176 int nx, ny, nz ;
01177
01178 int edge[6] = { 0,0,0,0,0,0 } ;
01179 int faceflag[ 12 ] ;
01180 int tot = 0 ;
01181 int cellflag[ 8 ] ;
01182
01183 int ct = 0 ;
01184 for ( i = -1 ; i < 1 ; i ++ )
01185 for ( j = -1 ; j < 1 ; j ++ )
01186 for ( k = -1 ; k < 1 ; k ++ )
01187 {
01188 if ( hasCell( ox + i, oy + j, oz + k ) )
01189 {
01190 cellflag[ ct ] = 1 ;
01191 }
01192 else
01193 {
01194 cellflag[ ct ] = 0 ;
01195 }
01196 ct ++ ;
01197 }
01198
01199 for ( i = 0 ; i < 12 ; i ++ )
01200 {
01201 faceflag[ i ] = 1 ;
01202 for ( j = 0 ; j < 4 ; j ++ )
01203 {
01204 nx = ox + sheetNeighbor[i][j][0] ;
01205 ny = oy + sheetNeighbor[i][j][1] ;
01206 nz = oz + sheetNeighbor[i][j][2] ;
01207
01208 if ( getDataAt( nx, ny, nz ) < 0 )
01209 {
01210 faceflag[ i ] = 0 ;
01211 break ;
01212 }
01213 }
01214
01215 if ( faceflag[ i ] )
01216 {
01217 if ( cellflag[ faceCells[i][0] ] ^ cellflag[ faceCells[i][1] ] )
01218 {
01219 int v1 = (int)( fvol->getDataAt(
01220 ox - 1 + (( faceCells[i][0] >> 2 ) & 1 ),
01221 oy - 1 + (( faceCells[i][0] >> 1 ) & 1 ),
01222 oz - 1 + (( faceCells[i][0] ) & 1)) ) ;
01223 int v2 = (int)( fvol->getDataAt(
01224 ox - 1 + (( faceCells[i][1] >> 2 ) & 1 ),
01225 oy - 1 + (( faceCells[i][1] >> 1 ) & 1 ),
01226 oz - 1 + (( faceCells[i][1] ) & 1)) ) ;
01227 if ( ((v1 >> (2 * (2 - i/4))) & 1) ||
01228 ((v2 >> (2 * (2 - i/4) + 1 )) & 1) )
01229 {
01230 faceflag[ i ] = 0 ;
01231 }
01232 }
01233 }
01234
01235 if ( faceflag[ i ] )
01236 {
01237 int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
01238 edge[ e0 ] ++ ;
01239 edge[ e1 ] ++ ;
01240 tot ++ ;
01241 }
01242 }
01243
01244 // Removing 1s
01245 int numones = 0 ;
01246 for ( i = 0 ; i < 6 ; i ++ )
01247 {
01248 if ( edge[ i ] == 1 )
01249 {
01250 numones ++ ;
01251 }
01252 }
01253 while( numones > 0 )
01254 {
01255 int e = 0;
01256 for ( i = 0 ; i < 6 ; i ++ )
01257 {
01258 if ( edge[ i ] == 1 )
01259 {
01260 e = i ;
01261 break ;
01262 }
01263 }
01264 /*
01265 if ( edge[ e ] != 1 )
01266 {
01267 printf("Wrong Again!********\n") ;
01268 }
01269 */
01270
01271 int f, e2 ;
01272 for ( j = 0 ; j < 4 ; j ++ )
01273 {
01274 f = edgeFaces[ e ][ j ] ;
01275 if ( faceflag[ f ] )
01276 {
01277 break ;
01278 }
01279 }
01280
01281 /*
01282 if ( faceflag[ f ] == 0 )
01283 {
01284 printf("Wrong!********\n") ;
01285 }
01286 */
01287
01288 if ( faceEdges[ f ][ 0 ] == e )
01289 {
01290 e2 = faceEdges[ f ][ 1 ] ;
01291 }
01292 else
01293 {
01294 e2 = faceEdges[ f ][ 0 ] ;
01295 }
01296
01297
01298 edge[ e ] -- ;
01299 numones -- ;
01300 edge[ e2 ] -- ;
01301 faceflag[ f ] = 0 ;
01302 tot -- ;
01303
01304 if ( edge[ e2 ] == 1 )
01305 {
01306 numones ++ ;
01307 }
01308 else if ( edge[ e2 ] == 0 )
01309 {
01310 numones -- ;
01311 }
01312 }
01313
01314 if ( tot > 0 )
01315 {
01316 return 1 ;
01317 }
01318
01319 return 0 ;
01320 }
|
|
||||||||||||||||
|
Definition at line 2135 of file volume.cpp. References getDataAt(), getNumNeighbor6(), nx, ny, and wustl_mm::SkeletonMaker::sheetNeighbor. Referenced by isSheetEnd(). 02136 {
02137 // return 1 ;
02138
02139 int i, j ;
02140 int nx, ny, nz ;
02141
02142 int faces = 12 ;
02143 for ( i = 0 ; i < 12 ; i ++ )
02144 {
02145 int ct = 0 ;
02146 for ( j = 0 ; j < 4 ; j ++ )
02147 {
02148 nx = ox + sheetNeighbor[i][j][0] ;
02149 ny = oy + sheetNeighbor[i][j][1] ;
02150 nz = oz + sheetNeighbor[i][j][2] ;
02151
02152 if ( getDataAt( nx, ny, nz ) < 0 )
02153 {
02154 ct = -1 ;
02155 break ;
02156 }
02157 else if ( getNumNeighbor6( nx, ny, nz ) == 6 )
02158 {
02159 ct = -1 ;
02160 break ;
02161 }
02162 // else if ( getDataAt( nx, ny, nz ) == 0 )
02163 // {
02164 // ct ++ ;
02165 // }
02166
02167
02168 }
02169 if ( ct == -1 || ct >= 1 )
02170 {
02171 faces -- ;
02172 }
02173 }
02174
02175 if ( faces > 0 )
02176 {
02177 return 1 ;
02178 }
02179 return 0 ;
02180
02181 }
|
|
||||||||||||||||
|
Definition at line 1789 of file volume.cpp. References getDataAt(), getNumNeighbor6(), wustl_mm::SkeletonMaker::neighbor6, nx, and ny. 01789 {
01790
01791 int i ;
01792 int c1 = 0 , c2 = 0 ;
01793 int nx, ny, nz ;
01794
01795 for ( i = 0 ; i < 6 ; i ++ )
01796 {
01797 nx = ox + neighbor6[i][0] ;
01798 ny = oy + neighbor6[i][1] ;
01799 nz = oz + neighbor6[i][2] ;
01800
01801 double val = getDataAt( nx, ny, nz ) ;
01802
01803 if ( val >= 0 )
01804 {
01805 c1 ++ ;
01806 if ( getNumNeighbor6(nx,ny,nz) < 6 ) // if ( val > 0 && val < MAX_ERODE )
01807 {
01808 c2 ++ ;
01809 }
01810 }
01811
01812 }
01813
01814 if ( c1 == 1 && c2 == 1 )
01815 {
01816 return 1 ;
01817 }
01818
01819 return 0 ;
01820 }
|
|
||||||||||||||||||||
|
Definition at line 1579 of file volume.cpp. References getDataAt(), MAX_ERODE, wustl_mm::SkeletonMaker::neighbor6, nx, and ny. Referenced by curveSkeleton(), curveSkeleton2D(), and surfaceSkeletonPres(). 01579 {
01580 // Returns 1 if it is a curve endpoint
01581 int i ;
01582 int c1 = 0 , c2 = 0 ;
01583 int nx, ny, nz ;
01584
01585 for ( i = 0 ; i < 6 ; i ++ )
01586 {
01587 nx = ox + neighbor6[i][0] ;
01588 ny = oy + neighbor6[i][1] ;
01589 nz = oz + neighbor6[i][2] ;
01590
01591 double val = getDataAt( nx, ny, nz ) ;
01592
01593 if ( val >= 0 )
01594 {
01595 c1 ++ ;
01596 if ( val > 0 && val < MAX_ERODE && nvol->getDataAt( nx, ny, nz ) == 0 )
01597 {
01598 c2 ++ ;
01599 }
01600 }
01601
01602 }
01603
01604 if ( c1 == 1 && c2 == 1 )
01605 {
01606 return 1 ;
01607 }
01608
01609 return 0 ;
01610 }
|
|
||||||||||||||||
|
Definition at line 2381 of file volume.cpp. References countExt(), countInt(), and getDataAt(). Referenced by curveSkeleton(). 02382 {
02383 /* Test if this is a simple voxel */
02384 // int flag = 0 ;
02385 double vox[3][3][3] ;
02386
02387 int i, j, k ;
02388 for ( i = -1 ; i < 2 ; i ++ )
02389 for ( j = -1 ; j < 2 ; j ++ )
02390 for ( k = -1 ; k < 2 ; k ++ )
02391 {
02392 double tval = getDataAt( ox + i, oy + j, oz + k ) ;
02393
02394 /*
02395 if ( tval == 0 || tval > (va )
02396 {
02397 flag = 1 ;
02398 }
02399 */
02400 /*
02401 if ( tval < 0 && tval == - curwid )
02402 {
02403 printf("Here %d", (int)-tval) ;
02404 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
02405 }
02406 else
02407 */
02408 {
02409 vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
02410 }
02411 }
02412
02413 /* Debugging
02414 printf("{") ;
02415 for ( i = 0 ; i < 3 ; i ++ )
02416 {
02417 if ( i ) printf(",") ;
02418 printf("{") ;
02419 for ( j = 0 ; j < 3 ; j ++ )
02420 {
02421 if ( j ) printf(",") ;
02422 printf("{") ;
02423 for ( k = 0 ; k < 3 ; k ++ )
02424 {
02425 if ( k ) printf(",") ;
02426 printf("%d", (vox[i][j][k] >=0 ? 1: 0));
02427 }
02428 printf("}") ;
02429 }
02430 printf("}") ;
02431 }
02432 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
02433 */
02434
02435 if ( countInt( vox ) == 1 && countExt( vox ) != 1 )
02436 {
02437 return 1 ;
02438 }
02439 else
02440 {
02441 return 0 ;
02442 }
02443 }
|
|
||||||||||||||||
|
* Definition at line 2238 of file volume.cpp. References hasCompleteSheet(), and isFeatureFace(). 02239 {
02240 return ( ( hasCompleteSheet(ox,oy,oz) == 0 ) && isFeatureFace(ox,oy,oz) ) ;
02241
02243 //
02244 //int i, j, k ;
02245 //int nx, ny, nz ;
02246
02247 //double vox[3][3][3] ;
02248 //for ( i = -1 ; i < 2 ; i ++ )
02249 // for ( j = -1 ; j < 2 ; j ++ )
02250 // for ( k = -1 ; k < 2 ; k ++ )
02251 // {
02252 // vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ;
02253 // }
02254
02255 //int edge[6] = { 4,4,4,4,4,4 } ;
02256 //int edge2[6] = { 4,4,4,4,4,4 } ;
02257
02258 //for ( i = 0 ; i < 12 ; i ++ )
02259 //{
02260 // for ( j = 0 ; j < 4 ; j ++ )
02261 // {
02262 // nx = 1 + sheetNeighbor[i][j][0] ;
02263 // ny = 1 + sheetNeighbor[i][j][1] ;
02264 // nz = 1 + sheetNeighbor[i][j][2] ;
02265
02266 // if ( vox[nx][ny][nz] < 0 )
02267 // {
02268 // edge[ faceEdges[ i ][ 0 ] ] -- ;
02269 // edge[ faceEdges[ i ][ 1 ] ] -- ;
02270 // break ;
02271 // }
02272 // }
02273
02274 // for ( j = 0 ; j < 4 ; j ++ )
02275 // {
02276 // nx = 1 + sheetNeighbor[i][j][0] ;
02277 // ny = 1 + sheetNeighbor[i][j][1] ;
02278 // nz = 1 + sheetNeighbor[i][j][2] ;
02279
02280 // if ( vox[nx][ny][nz] <= 0 )
02281 // {
02282 // edge2[ faceEdges[ i ][ 0 ] ] -- ;
02283 // edge2[ faceEdges[ i ][ 1 ] ] -- ;
02284 // break ;
02285 // }
02286 // }
02287 //}
02288
02289 //
02291 //for ( i = 0 ; i < 6 ; i ++ )
02292 //{
02293 // nx = 1 + neighbor6[i][0] ;
02294 // ny = 1 + neighbor6[i][1] ;
02295 // nz = 1 + neighbor6[i][2] ;
02296 // if ( edge[i] == 0 && vox[nx][ny][nz] >= 0 )
02297 // {
02298 // return 0 ;
02299 // }
02300 //}
02301 //*/
02302 //
02303
02304
02305 //for ( i = 0 ; i < 6 ; i ++ )
02306 //{
02307 // if ( edge[ i ] == 1 && edge2[ i ] == 1 )
02308 // {
02309 // return 1 ;
02310 // }
02311 //}
02312
02313 //return 0 ;
02314 }
|
|
||||||||||||||||||||
|
Definition at line 1821 of file volume.cpp. References wustl_mm::SkeletonMaker::edgeFaces, wustl_mm::SkeletonMaker::faceEdges, getDataAt(), nx, ny, wustl_mm::SkeletonMaker::sheetNeighbor, and tot. Referenced by surfaceSkeletonPres(). 01822 {
01823 // Returns 1 if it contains a sheet boundary. Noise-resistant
01824 int i, j ;
01825 int nx, ny, nz ;
01826
01827 int edge[6] = { 0,0,0,0,0,0 } ;
01828 int faceflag[ 12 ] ;
01829 int hasFeatureFace = 0 ;
01830 int tot = 0 ;
01831
01832 for ( i = 0 ; i < 12 ; i ++ )
01833 {
01834 faceflag[ i ] = 1 ;
01835 int hasFeature = 1 ;
01836
01837 for ( j = 0 ; j < 4 ; j ++ )
01838 {
01839 nx = ox + sheetNeighbor[i][j][0] ;
01840 ny = oy + sheetNeighbor[i][j][1] ;
01841 nz = oz + sheetNeighbor[i][j][2] ;
01842
01843 if ( getDataAt( nx, ny, nz ) == 0 || nvol->getDataAt( nx, ny, nz ) == 1 )
01844 {
01845 hasFeature = 0 ;
01846 }
01847 if ( getDataAt( nx, ny, nz ) < 0 )
01848 {
01849 faceflag[ i ] = 0 ;
01850 break ;
01851 }
01852 }
01853 if ( faceflag[ i ] == 1 && hasFeature )
01854 {
01855 hasFeatureFace ++ ;
01856 // return 0 ;
01857 }
01858
01859 if ( faceflag[ i ] )
01860 {
01861 int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
01862 edge[ e0 ] ++ ;
01863 edge[ e1 ] ++ ;
01864 tot ++ ;
01865 }
01866 }
01867
01868 if ( tot == 0 || hasFeatureFace == 0 )
01869 {
01870 return 0 ;
01871 }
01872
01873 // Removing 1s
01874 int numones = 0 ;
01875 for ( i = 0 ; i < 6 ; i ++ )
01876 {
01877 if ( edge[ i ] == 1 )
01878 {
01879 numones ++ ;
01880 }
01881 }
01882 while( numones > 0 )
01883 {
01884 int e = 0;
01885 for ( i = 0 ; i < 6 ; i ++ )
01886 {
01887 if ( edge[ i ] == 1 )
01888 {
01889 e = i ;
01890 break ;
01891 }
01892 }
01893 /*
01894 if ( edge[ e ] != 1 )
01895 {
01896 printf("Wrong Again!********\n") ;
01897 }
01898 */
01899
01900 int f, e2 ;
01901 for ( j = 0 ; j < 4 ; j ++ )
01902 {
01903 f = edgeFaces[ e ][ j ] ;
01904 if ( faceflag[ f ] )
01905 {
01906 break ;
01907 }
01908 }
01909
01910 /*
01911 if ( faceflag[ f ] == 0 )
01912 {
01913 printf("Wrong!********\n") ;
01914 }
01915 */
01916
01917 if ( faceEdges[ f ][ 0 ] == e )
01918 {
01919 e2 = faceEdges[ f ][ 1 ] ;
01920 }
01921 else
01922 {
01923 e2 = faceEdges[ f ][ 0 ] ;
01924 }
01925
01926
01927 edge[ e ] -- ;
01928 numones -- ;
01929 edge[ e2 ] -- ;
01930 faceflag[ f ] = 0 ;
01931 tot -- ;
01932
01933 if ( edge[ e2 ] == 1 )
01934 {
01935 numones ++ ;
01936 }
01937 else if ( edge[ e2 ] == 0 )
01938 {
01939 numones -- ;
01940 }
01941 }
01942
01943 if ( tot > 0 )
01944 {
01945 return 0 ;
01946 }
01947
01948 return 1 ;
01949 }
|
|
||||||||||||||||
|
Definition at line 2316 of file volume.cpp. References countExt(), countInt(), and getDataAt(). Referenced by curveSkeleton(), curveSkeleton2D(), skeleton(), and surfaceSkeletonPres(). 02317 {
02318 /* Test if this is a simple voxel */
02319 // int flag = 0 ;
02320 double vox[3][3][3] ;
02321
02322 int i, j, k ;
02323 for ( i = -1 ; i < 2 ; i ++ )
02324 for ( j = -1 ; j < 2 ; j ++ )
02325 for ( k = -1 ; k < 2 ; k ++ )
02326 {
02327 double tval = getDataAt( ox + i, oy + j, oz + k ) ;
02328
02329 /*
02330 if ( tval == 0 || tval > (va )
02331 {
02332 flag = 1 ;
02333 }
02334 */
02335 /*
02336 if ( tval < 0 && tval == - curwid )
02337 {
02338 printf("Here %d", (int)-tval) ;
02339 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
02340 }
02341 else
02342 */
02343 {
02344 vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
02345 }
02346 }
02347
02348 /* Debugging
02349 printf("{") ;
02350 for ( i = 0 ; i < 3 ; i ++ )
02351 {
02352 if ( i ) printf(",") ;
02353 printf("{") ;
02354 for ( j = 0 ; j < 3 ; j ++ )
02355 {
02356 if ( j ) printf(",") ;
02357 printf("{") ;
02358 for ( k = 0 ; k < 3 ; k ++ )
02359 {
02360 if ( k ) printf(",") ;
02361 printf("%d", (vox[i][j][k] >=0 ? 1: 0));
02362 }
02363 printf("}") ;
02364 }
02365 printf("}") ;
02366 }
02367 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
02368 */
02369
02370 if ( countInt( vox ) == 1 && countExt( vox ) == 1 )
02371 {
02372 return 1 ;
02373 }
02374 else
02375 {
02376 return 0 ;
02377 }
02378 }
|
|
|
Definition at line 1027 of file volume.cpp. References getDataAt(), getSizeX(), getSizeY(), getSizeZ(), hasCell(), wustl_mm::SkeletonMaker::neighbor6, nx, ny, setDataAt(), and Volume(). Referenced by erodeSheet(). 01027 {
01028 int i, j, k ;
01029 Volume* fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
01030
01031 //return fvol ;
01032
01033 for ( i = 0 ; i < getSizeX() ; i ++ )
01034 for ( j = 0 ; j < getSizeY() ; j ++ )
01035 for ( k = 0 ; k < getSizeZ() ; k ++ )
01036 {
01037 if( getDataAt(i,j,k) >= 0 )
01038 {
01039 if ( hasCell(i,j,k) )
01040 {
01041 for ( int m = 0 ; m < 6 ; m ++ )
01042 {
01043 int nx = i + neighbor6[m][0] ;
01044 int ny = j + neighbor6[m][1] ;
01045 int nz = k + neighbor6[m][2] ;
01046 if ( ! hasCell(nx,ny,nz) )
01047 {
01048 fvol->setDataAt(i,j,k,(double)(1<<m)) ;
01049 break ;
01050 }
01051 }
01052 }
01053 }
01054 }
01055
01056
01057 return fvol ;
01058 }
|
|
||||||||||||
|
* Next, smooth */
Definition at line 273 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::Pad(), and volData. Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::PerformPureJuSkeletonization().
|
|
||||||||||||
|
Definition at line 56 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::SetDataAt(), and volData.
|
|
||||||||||||||||||||
|
||||||||||||||||
|
Definition at line 77 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::SetOrigin(), and volData.
|
|
||||||||||||||||
|
Definition at line 73 of file volume.cpp. References wustl_mm::SkeletonMaker::VolumeData::SetSpacing(), and volData. 00073 {
00074 volData->SetSpacing(spx, spy, spz);
00075 }
|
|
||||||||||||||||
|
* Adding ends *// Move its neighboring unvisited node to queue2
Definition at line 5692 of file volume.cpp. References EMAN::PriorityQueue< ValueT, KeyT >::add(), getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex2(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isSimple(), MAX_ERODE, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z. 05693 {
05694 int i, j, k ;
05695 // First, threshold the volume
05696 #ifdef VERBOSE
05697 printf("Thresholding the volume to -1/0...\n") ;
05698 #endif
05699 threshold( thr, -1, 0 ) ;
05700
05701 // Next, apply convergent erosion
05702 // by preserving: complex nodes, curve end-points, and sheet points
05703
05704 // Next, initialize the linked queue
05705 #ifdef VERBOSE
05706 printf("Initializing queue...\n") ;
05707 #endif
05708 GridQueue2* queue2 = new GridQueue2( ) ;
05709 GridQueue2* queue3 = new GridQueue2( ) ;
05710 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
05711
05712 for ( i = 0 ; i < getSizeX() ; i ++ )
05713 for ( j = 0 ; j < getSizeY() ; j ++ )
05714 for ( k = 0 ; k < getSizeZ() ; k ++ )
05715 {
05716 if ( getDataAt( i, j, k ) >= 0 )
05717 {
05718 if ( svol->getDataAt(i,j,k) > 0 || hvol->getDataAt(i,j,k) > 0 )
05719 {
05720 setDataAt( i, j, k, MAX_ERODE ) ;
05721 }
05722 else
05723 {
05724 for ( int m = 0 ; m < 6 ; m ++ )
05725 {
05726 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
05727 {
05728 setDataAt( i, j, k, 1 ) ;
05729 queue2->prepend( i, j, k ) ;
05730 break ;
05731 }
05732 }
05733 }
05734 }
05735 }
05736 int wid = MAX_ERODE ;
05737 #ifdef VERBOSE
05738 printf("Total %d nodes\n", queue2->getNumElements() ) ;
05739
05740
05741 // Perform erosion
05742 printf("Start erosion to %d...\n", wid) ;
05743 #endif
05744 gridQueueEle* ele ;
05745 gridPoint* gp ;
05746 int ox, oy, oz ;
05747 int score ;
05748 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
05749 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
05750 {
05751 scrvol->setDataAt( i, -1 ) ;
05752 }
05753
05754
05755 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
05756 {
05757 // At the start of each iteration,
05758 // queue2 holds all the nodes for this layer
05759 // queue3 and queue are empty
05760
05761 int numComplex = 0, numSimple = 0 ;
05762 #ifdef VERBOSE
05763 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
05764 #endif
05765
05766
05767 // Next,
05768 // Compute score for each node left in queue2
05769 // move them into priority queue
05770 queue2->reset() ;
05771 ele = queue2->getNext() ;
05772 while ( ele != NULL )
05773 {
05774 ox = ele->x ;
05775 oy = ele->y ;
05776 oz = ele->z ;
05777
05778 // Compute score
05779 score = getNumPotComplex2( ox, oy, oz ) ;
05780 scrvol->setDataAt( ox, oy, oz, score ) ;
05781
05782 // Push to queue
05783 gp = new gridPoint ;
05784 gp->x = ox ;
05785 gp->y = oy ;
05786 gp->z = oz ;
05787 // queue->add( gp, -score ) ;
05788 queue->add( gp, score ) ;
05789
05790 ele = queue2->remove() ;
05791 }
05792
05793 // Rename queue3 to be queue2,
05794 // Clear queue3
05795 // From now on, queue2 holds nodes of next level
05796 delete queue2 ;
05797 queue2 = queue3 ;
05798 queue3 = new GridQueue2( ) ;
05799
05800 // Next, start priority queue iteration
05801 while ( ! queue->isEmpty() )
05802 {
05803 // Retrieve the node with the highest score
05804 queue->remove( gp, score ) ;
05805 ox = gp->x ;
05806 oy = gp->y ;
05807 oz = gp->z ;
05808 delete gp ;
05809 // score = -score ;
05810
05811 // Ignore the node
05812 // if it has been processed before
05813 // or it has an updated score
05814 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
05815 {
05816 continue ;
05817 }
05818
05819 /* Added for debugging */
05820 // Check simple
05821 if ( ! isSimple( ox, oy, oz ) )
05822 {
05823 // Complex, set to next layer
05824 setDataAt( ox, oy, oz, curwid + 1 ) ;
05825 queue2->prepend( ox, oy, oz ) ;
05826 numComplex ++ ;
05827 }
05828 else
05829 {
05830 setDataAt( ox, oy, oz, -1 ) ;
05831 numSimple ++ ;
05832 }
05833 /* Adding ends */
05834
05835 // Move its neighboring unvisited node to queue2
05836 for ( int m = 0 ; m < 6 ; m ++ )
05837 {
05838 int nx = ox + neighbor6[m][0] ;
05839 int ny = oy + neighbor6[m][1] ;
05840 int nz = oz + neighbor6[m][2] ;
05841 if ( getDataAt( nx, ny, nz ) == 0 )
05842 {
05843 setDataAt( nx, ny, nz, curwid + 1 ) ;
05844 queue2->prepend( nx, ny, nz ) ;
05845 }
05846 }
05847
05848 // Update scores for nodes in its 5x5 neighborhood
05849 // insert them back into priority queue
05850
05851 for ( i = -2 ; i < 3 ;i ++ )
05852 for ( j = -2 ; j < 3 ; j ++ )
05853 for ( k = -2 ; k < 3 ; k ++ )
05854 {
05855 int nx = ox + i ;
05856 int ny = oy + j ;
05857 int nz = oz + k ;
05858
05859 if ( getDataAt( nx, ny, nz ) == curwid )
05860 {
05861 // Compute score
05862 score = getNumPotComplex2( nx, ny, nz ) ;
05863
05864 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
05865 {
05866 // printf("Update\n") ;
05867 scrvol->setDataAt( nx, ny, nz, score ) ;
05868 // Push to queue
05869 gp = new gridPoint ;
05870 gp->x = nx ;
05871 gp->y = ny ;
05872 gp->z = nz ;
05873 // queue->add( gp, -score ) ;
05874 queue->add( gp, score ) ;
05875 }
05876 }
05877 }
05878
05879
05880 }
05881
05882 #ifdef VERBOSE
05883 printf("%d complex, %d simple\n", numComplex, numSimple) ;
05884 #endif
05885
05886 if ( numSimple == 0 )
05887 {
05888 delete queue2 ;
05889 break ;
05890 }
05891 }
05892
05893 // Finally, clean up
05894 #ifdef VERBOSE
05895 printf("Thresholding the volume to 0/1...\n") ;
05896 #endif
05897 threshold( 0, 0, 1 ) ;
05898 delete scrvol;
05899 delete queue;
05900 delete queue3;
05901 }
|
|
||||||||||||
|
Definition at line 5089 of file volume.cpp. References getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getSizeX(), getSizeY(), getSizeZ(), isSimple(), wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridQueueEle::y, and wustl_mm::SkeletonMaker::gridQueueEle::z. Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning(). 05090 {
05091 int i, j, k ;
05092 // First, threshold the volume
05093 #ifdef VERBOSE
05094 printf("Thresholding the volume to -1/0...\n") ;
05095 #endif
05096 threshold( thr, -1, 0 ) ;
05097
05098 // Next, apply convergent erosion
05099 // by preserving: complex nodes, curve end-points, and sheet points
05100
05101 // Next, initialize the linked queue
05102 #ifdef VERBOSE
05103 printf("Initializing queue...\n") ;
05104 #endif
05105 GridQueue2* queue2 = new GridQueue2( ) ;
05106 GridQueue2* queue = new GridQueue2( ) ;
05107
05108 for ( i = 0 ; i < getSizeX() ; i ++ )
05109 for ( j = 0 ; j < getSizeY() ; j ++ )
05110 for ( k = 0 ; k < getSizeZ() ; k ++ )
05111 {
05112 if ( getDataAt( i, j, k ) >= 0 )
05113 {
05114 {
05115 for ( int m = 0 ; m < 6 ; m ++ )
05116 {
05117 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
05118 {
05119 setDataAt( i, j, k, 1 ) ;
05120 queue2->prepend( i, j, k ) ;
05121 break ;
05122 }
05123 }
05124 }
05125 }
05126 }
05127 int wid = 0 ;
05128 #ifdef VERBOSE
05129 printf("Total %d nodes\n", queue2->getNumElements() ) ;
05130 printf("Start erosion to %d...\n", wid) ;
05131 #endif
05132
05133 // Perform erosion
05134 gridQueueEle* ele ;
05135 int ox, oy, oz ;
05136
05137 while( 1 )
05138 {
05139 wid ++ ;
05140
05141 // At the start of each iteration,
05142 // queue2 holds all the nodes for this layer
05143 // queue is empty
05144
05145 int numComplex = 0, numSimple = 0 ;
05146 #ifdef VERBOSE
05147 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), wid) ;
05148 #endif
05149
05150 // Rename queue2 to be queue,
05151 // Clear queue2
05152 // From now on, queue2 holds nodes of next level
05153 delete queue ;
05154 queue = queue2 ;
05155 queue2 = new GridQueue2( ) ;
05156
05157 // Next, start queue iteration
05158 queue->reset() ;
05159 ele = queue->getNext();
05160 while ( ele != NULL )
05161 {
05162 ox = ele->x ;
05163 oy = ele->y ;
05164 oz = ele->z ;
05165 // delete ele ;
05166
05167 // Check simple
05168 if ( ! isSimple( ox, oy, oz ) )
05169 {
05170 // Complex, set to next layer
05171 queue2->prepend( ox, oy, oz ) ;
05172 numComplex ++ ;
05173 }
05174 /*
05175 else if ( ox == off || oy == off || oz == off ||
05176 ox == getSizeX() - off - 1 || oy == getSizeY() - off - 1 || oz == getSizeZ() - off - 1 )
05177 {
05178 // Wall, don't erode, set to next layer
05179 queue2->prepend( ox, oy, oz ) ;
05180 numComplex ++ ;
05181 }
05182 */
05183 else
05184 {
05185 setDataAt( ox, oy, oz, -1 ) ;
05186 numSimple ++ ;
05187
05188 for ( int m = 0 ; m < 6 ; m ++ )
05189 {
05190 int nx = ox + neighbor6[m][0] ;
05191 int ny = oy + neighbor6[m][1] ;
05192 int nz = oz + neighbor6[m][2] ;
05193 if ( getDataAt( nx, ny, nz ) == 0 )
05194 {
05195 setDataAt( nx, ny, nz, 1 ) ;
05196 queue2->prepend( nx, ny, nz ) ;
05197 }
05198 }
05199
05200 }
05201
05202 ele = queue->remove() ;
05203 }
05204 #ifdef VERBOSE
05205 printf("Level %d: %d complex, %d simple\n", wid, numComplex, numSimple) ;
05206 #endif
05207
05208 if ( numSimple == 0 )
05209 {
05210 break ;
05211 }
05212 }
05213
05214 // Finally, clean up
05215 delete queue;
05216 delete queue2;
05217 #ifdef VERBOSE
05218 printf("Thresholding the volume to 0/1...\n") ;
05219 #endif
05220 threshold( 0, 0, 1 ) ;
05221 }
|
|
||||||||||||
|
* Commented for debugging/ Find complex nodes in its 3x3 neighborhood/ move them to queue2
Definition at line 8672 of file volume.cpp. References EMAN::PriorityQueue< ValueT, KeyT >::add(), flag, getDataAt(), wustl_mm::SkeletonMaker::GridQueue2::getNext(), wustl_mm::SkeletonMaker::GridQueue2::getNumElements(), getNumPotComplex(), getSizeX(), getSizeY(), getSizeZ(), EMAN::PriorityQueue< ValueT, KeyT >::isEmpty(), isHelixEnd(), isSheetEnd(), isSimple(), MAX_ERODE, wustl_mm::SkeletonMaker::neighbor6, nx, ny, wustl_mm::SkeletonMaker::GridQueue2::prepend(), EMAN::PriorityQueue< ValueT, KeyT >::remove(), wustl_mm::SkeletonMaker::GridQueue2::remove(), wustl_mm::SkeletonMaker::GridQueue2::reset(), setDataAt(), threshold(), Volume(), wustl_mm::SkeletonMaker::gridPoint::x, wustl_mm::SkeletonMaker::gridQueueEle::x, wustl_mm::SkeletonMaker::gridPoint::y, wustl_mm::SkeletonMaker::gridQueueEle::y, wustl_mm::SkeletonMaker::gridPoint::z, and wustl_mm::SkeletonMaker::gridQueueEle::z. Referenced by wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::GetJuThinning(). 08673 {
08674 int i, j, k ;
08675 // First, threshold the volume
08676 #ifdef VERBOSE
08677 printf("Thresholding the volume to -MAX_ERODE/0...\n") ;
08678 #endif
08679 threshold( thr, -MAX_ERODE, 0 ) ;
08680
08681 // Next, initialize the linked queue
08682 #ifdef VERBOSE
08683 printf("Initializing queue...\n") ;
08684 #endif
08685 GridQueue2* queue2 = new GridQueue2( ) ;
08686 GridQueue2* queue3 = new GridQueue2( ) ;
08687 GridQueue2* queue4 = new GridQueue2( ) ;
08688
08689 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
08690
08691 for ( i = 0 ; i < getSizeX() ; i ++ )
08692 for ( j = 0 ; j < getSizeY() ; j ++ )
08693 for ( k = 0 ; k < getSizeZ() ; k ++ )
08694 {
08695 if ( getDataAt( i, j, k ) >= 0 ) {
08696 if(preserve->getDataAt(i, j, k) > 0) {
08697 setDataAt(i, j, k, MAX_ERODE);
08698 } else {
08699 for ( int m = 0 ; m < 6 ; m ++ )
08700 {
08701 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
08702 {
08703 // setDataAt( i, j, k, 1 ) ;
08704 queue2->prepend( i, j, k ) ;
08705 break ;
08706 }
08707 }
08708 }
08709 }
08710 }
08711 int wid = MAX_ERODE ;
08712 #ifdef VERBOSE
08713 printf("Total %d nodes\n", queue2->getNumElements() ) ;
08714 printf("Start erosion to %d...\n", wid) ;
08715 #endif
08716
08717
08718 // Perform erosion
08719 gridQueueEle* ele ;
08720 gridPoint* gp ;
08721 int ox, oy, oz ;
08722 int score ;
08723 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
08724 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
08725 {
08726 scrvol->setDataAt( i, -1 ) ;
08727 }
08728
08729 #ifdef NOISE_DIS_SHEET
08730 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
08731 #endif
08732
08733 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
08734 {
08735 // At the start of each iteration,
08736 // queue2 and queue4 holds all the nodes for this layer
08737 // queue3 and queue are empty
08738
08739 int numComplex = 0, numSimple = 0 ;
08740 #ifdef VERBOSE
08741 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
08742 #endif
08743
08744 /*
08745 We first need to assign curwid + 1 to every node in this layer
08746 */
08747 queue2->reset() ;
08748 ele = queue2->getNext() ;
08749 while ( ele != NULL )
08750 {
08751 ox = ele->x ;
08752 oy = ele->y ;
08753 oz = ele->z ;
08754
08755 if ( getDataAt(ox,oy,oz) == curwid )
08756 {
08757 ele = queue2->remove() ;
08758 }
08759 else
08760 {
08761 setDataAt(ox,oy,oz, curwid) ;
08762 ele = queue2->getNext() ;
08763 }
08764 }
08765 queue4->reset() ;
08766 ele = queue4->getNext() ;
08767 while ( ele != NULL )
08768 {
08769 ox = ele->x ;
08770 oy = ele->y ;
08771 oz = ele->z ;
08772
08773 queue2->prepend(ox,oy,oz) ;
08774 ele = queue4->remove() ;
08775 }
08776
08777 // Now queue2 holds all the nodes for this layer
08778
08779 #ifdef NOISE_DIS_SHEET
08780 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
08781 queue2->reset() ;
08782
08783 // First run
08784 int flag = 0 ;
08785 while ( ( ele = queue2->getNext() ) != NULL )
08786 {
08787 ox = ele->x ;
08788 oy = ele->y ;
08789 oz = ele->z ;
08790 if ( NOISE_DIS_SHEET <= 1 )
08791 {
08792 noisevol->setDataAt( ox, oy, oz, 0 ) ;
08793 }
08794 else
08795 {
08796 flag = 0 ;
08797 for ( int m = 0 ; m < 6 ; m ++ )
08798 {
08799 int nx = ox + neighbor6[m][0] ;
08800 int ny = oy + neighbor6[m][1] ;
08801 int nz = oz + neighbor6[m][2] ;
08802 if ( getDataAt( nx, ny, nz ) == 0 )
08803 {
08804 noisevol->setDataAt( ox, oy, oz, 1 ) ;
08805 flag = 1 ;
08806 break ;
08807 }
08808 }
08809 if ( ! flag )
08810 {
08811 noisevol->setDataAt( ox, oy, oz, 0 ) ;
08812 }
08813 }
08814 }
08815
08816 for ( int cur = 1 ; cur < NOISE_DIS_SHEET ; cur ++ )
08817 {
08818 queue2->reset() ;
08819 int count = 0 ;
08820
08821 while ( ( ele = queue2->getNext() ) != NULL )
08822 {
08823 ox = ele->x ;
08824 oy = ele->y ;
08825 oz = ele->z ;
08826
08827 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
08828 {
08829 continue ;
08830 }
08831
08832 flag = 0 ;
08833 for ( int m = 0 ; m < 6 ; m ++ )
08834 {
08835 int nx = ox + neighbor6[m][0] ;
08836 int ny = oy + neighbor6[m][1] ;
08837 int nz = oz + neighbor6[m][2] ;
08838 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
08839 {
08840 noisevol->setDataAt( ox, oy, oz, 1 ) ;
08841 count ++ ;
08842 break ;
08843 }
08844 }
08845 }
08846
08847 if ( count == 0 )
08848 {
08849 break ;
08850 }
08851 }
08852
08853
08854 #endif
08855
08856 /* Commented for debugging
08857
08858 // First,
08859 // check for complex nodes in queue2
08860 // move them from queue2 to queue3
08861 queue2->reset() ;
08862 ele = queue2->getNext() ;
08863 while ( ele != NULL )
08864 {
08865 ox = ele->x ;
08866 oy = ele->y ;
08867 oz = ele->z ;
08868
08869 // Check simple
08870 #ifndef NOISE_DIS_SHEET
08871 if ( isSheetEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
08872 #else
08873 if ( isSheetEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
08874 #endif
08875 {
08876 // Complex, set to next layer
08877 setDataAt( ox, oy, oz, curwid + 1 ) ;
08878 queue3->prepend( ox, oy, oz ) ;
08879 ele = queue2->remove() ;
08880
08881 numComplex ++ ;
08882 }
08883 else
08884 {
08885 ele = queue2->getNext() ;
08886 }
08887 }
08888 */
08889
08890
08891 // Next,
08892 // Compute score for each node left in queue2
08893 // move them into priority queue
08894 queue2->reset() ;
08895 ele = queue2->getNext() ;
08896 while ( ele != NULL )
08897 {
08898 ox = ele->x ;
08899 oy = ele->y ;
08900 oz = ele->z ;
08901
08902 // Compute score
08903 score = getNumPotComplex( ox, oy, oz ) ;
08904 scrvol->setDataAt( ox, oy, oz, score ) ;
08905
08906 // Push to queue
08907 gp = new gridPoint ;
08908 gp->x = ox ;
08909 gp->y = oy ;
08910 gp->z = oz ;
08911 // queue->add( gp, -score ) ;
08912 queue->add( gp, score ) ;
08913
08914 ele = queue2->remove() ;
08915 }
08916
08917 // Rename queue3 to be queue2,
08918 // Clear queue3
08919 // From now on, queue2 holds nodes of next level
08920 delete queue2 ;
08921 queue2 = queue3 ;
08922 queue3 = new GridQueue2( ) ;
08923
08924
08925 // Next, start priority queue iteration
08926 while ( ! queue->isEmpty() )
08927 {
08928 // Retrieve the node with the highest score
08929 queue->remove( gp, score ) ;
08930 ox = gp->x ;
08931 oy = gp->y ;
08932 oz = gp->z ;
08933 delete gp ;
08934 // printf("%d\n", score);
08935 // score = -score ;
08936
08937 // Ignore the node
08938 // if it has been processed before
08939 // or it has an updated score
08940 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
08941 {
08942 continue ;
08943 }
08944
08945 /* Commented for debugging
08946
08947 // Remove this simple node
08948 setDataAt( ox, oy, oz, -1 ) ;
08949 numSimple ++ ;
08950 // printf("Highest score: %d\n", score) ;
08951 */
08952
08953 /* Added for debugging */
08954 // Check simple
08955 #ifndef NOISE_DIS_SHEET
08956 // if ( hasFeatureFace( ox, oy, oz ) )
08957 if ( (! isSimple( ox, oy, oz )) || isSheetEnd( ox, oy, oz ) )
08958 // if ( hasIsolatedFace(ox,oy,oz) && (! isNoiseSheetEnd(ox,oy,oz)))
08959 #else
08960 // if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz, noisevol ) )
08961 if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz, noisevol ) || isHelixEnd( ox, oy, oz, noisevol ))
08962 // if ( isBertrandEndPoint( ox, oy, oz ) )
08963 #endif
08964 {
08965 // Complex, set to next layer
08966 setDataAt( ox, oy, oz, curwid + 1 ) ;
08967 queue4->prepend( ox, oy, oz ) ;
08968 numComplex ++ ;
08969
08970 }
08971 else
08972 {
08973 setDataAt( ox, oy, oz, -1 ) ;
08974 numSimple ++ ;
08975
08976 }
08977 /* Adding ends */
08978
08979 // Move its neighboring unvisited node to queue2
08980 for ( int m = 0 ; m < 6 ; m ++ )
08981 {
08982 int nx = ox + neighbor6[m][0] ;
08983 int ny = oy + neighbor6[m][1] ;
08984 int nz = oz + neighbor6[m][2] ;
08985 if ( getDataAt( nx, ny, nz ) == 0 )
08986 {
08987 // setDataAt( nx, ny, nz, curwid + 1 ) ;
08988 queue2->prepend( nx, ny, nz ) ;
08989 }
08990 }
08991
08992 /* Commented for debugging
08993
08994 // Find complex nodes in its 3x3 neighborhood
08995 // move them to queue2
08996 for ( i = -1 ; i < 2 ; i ++ )
08997 for ( j = -1 ; j < 2 ; j ++ )
08998 for ( k = -1 ; k < 2 ; k ++ )
08999 {
09000 int nx = ox + i ;
09001 int ny = oy + j ;
09002 int nz = oz + k ;
09003
09004 // Check simple
09005 if ( getDataAt( nx, ny, nz ) == curwid &&
09006 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
09007 #ifndef NOISE_DIS_SHEET
09008 ( isSheetEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
09009 #else
09010 ( isSheetEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
09011 #endif
09012
09013 {
09014 // Complex, set to next layer
09015 setDataAt( nx, ny, nz, curwid + 1 ) ;
09016 queue2->prepend( nx, ny, nz ) ;
09017 numComplex ++ ;
09018 }
09019 }
09020 */
09021
09022 // Update scores for nodes in its 5x5 neighborhood
09023 // insert them back into priority queue
09024
09025 for ( i = -2 ; i < 3 ;i ++ )
09026 for ( j = -2 ; j < 3 ; j ++ )
09027 for ( k = -2 ; k < 3 ; k ++ )
09028 {
09029 int nx = ox + i ;
09030 int ny = oy + j ;
09031 int nz = oz + k ;
09032
09033 if ( getDataAt( nx, ny, nz ) == curwid )
09034 {
09035 // Compute score
09036 score = getNumPotComplex( nx, ny, nz ) ;
09037
09038 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
09039 {
09040 // printf("Update\n") ;
09041 scrvol->setDataAt( nx, ny, nz, score ) ;
09042 // Push to queue
09043 gp = new gridPoint ;
09044 gp->x = nx ;
09045 gp->y = ny ;
09046 gp->z = nz ;
09047 // queue->add( gp, -score ) ;
09048 queue->add( gp, score ) ;
09049 }
09050 }
09051 }
09052
09053
09054 }
09055
09056 #ifdef VERBOSE
09057 printf("%d complex, %d simple\n", numComplex, numSimple) ;
09058 #endif
09059
09060 if ( numSimple == 0 )
09061 {
09062 break ;
09063 }
09064 }
09065
09066 // Finally, clean up
09067 #ifdef VERBOSE
09068 printf("Thresholding the volume to 0/1...\n") ;
09069 #endif
09070 threshold( 0, 0, 1 ) ;
09071
09072 delete scrvol;
09073 delete queue;
09074 delete queue2;
09075 delete queue3;
09076 delete queue4;
09077
09078 }
|
|
||||||||||||||||||||||||
|
Definition at line 9530 of file volume.cpp. References getDataAt(), getSizeX(), getSizeY(), getSizeZ(), in, and setDataAt(). 09531 {
09532 float val;
09533 for ( int i = 0 ; i < getSizeX() ; i ++ )
09534 for ( int j = 0 ; j < getSizeY() ; j ++ )
09535 for ( int k = 0 ; k < getSizeZ() ; k ++ )
09536 {
09537 val = (float)getDataAt(i, j, k);
09538 if(markBoundary) {
09539 if ( i > 1 && i < getSizeX() - 2 && j > 1 && j < getSizeY() - 2 && k > 1 && k < getSizeZ() - 2 ) {
09540 if(val < thr) {
09541 setDataAt(i, j, k, out);
09542 } else {
09543 setDataAt(i, j, k, in);
09544 }
09545 }
09546 else
09547 {
09548 setDataAt(i, j, k, boundary);
09549 }
09550 } else {
09551 if(val < thr) {
09552 setDataAt(i, j, k, out);
09553 } else {
09554 setDataAt(i, j, k, in);
09555 }
09556 }
09557 }
09558 }
|
|
||||||||||||||||||||
|
Definition at line 9525 of file volume.cpp. References in, and threshold(). 09526 {
09527 threshold(thr, out, in, boundary, true);
09528 }
|
|
||||||||||||||||
|
Definition at line 9520 of file volume.cpp. References in, and threshold(). 09521 {
09522 threshold( thr, out, in, out, true) ;
09523 }
|
|
|
Normalize to a given range.
Definition at line 9515 of file volume.cpp. Referenced by curveSkeleton(), curveSkeleton2D(), erodeHelix(), erodeSheet(), skeleton(), surfaceSkeletonPres(), and threshold(). 09516 {
09517 threshold( thr, 0, 1, 0, true) ;
09518 }
|
|
|
Definition at line 222 of file volume.h. Referenced by getDataAt(), getIndex(), getOriginX(), getOriginY(), getOriginZ(), getSizeX(), getSizeY(), getSizeZ(), getSpacingX(), getSpacingY(), getSpacingZ(), pad(), setDataAt(), setOrigin(), setSpacing(), and Volume(). |
1.3.9.1