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