EMAN2
volume.cpp
Go to the documentation of this file.
1// Copyright (C) 2005-2008 Washington University in St Louis, Baylor College of Medicine. All rights reserved
2// Author: Tao Ju (taoju@cse.wustl.edu), Refactored by Sasakthi Abeysinghe (sasakthi.abeysinghe@wustl.edu)
3// Description: Volumetric data definition
4
5#include "volume.h"
6
7using namespace wustl_mm::SkeletonMaker;
8
9 Volume::Volume(EMData* em) //eman2
10 {
11 this->volData = new VolumeData(em);
12 }
13 Volume::Volume(int x, int y, int z) {
14 volData = new VolumeData(x, y, z);
15 }
16
17 Volume::Volume(int x, int y, int z, float val) {
18 volData = new VolumeData(x, y, z, val);
19 }
20
21 Volume::Volume(int x, int y, int z, int offx, int offy, int offz, Volume * vol) {
22 volData = new VolumeData(x, y, z, offx, offy, offz, vol->getVolumeData());
23 }
24
25
26 Volume::~Volume( )
27 {
28 delete volData;
29 }
30
31 EMData* Volume::get_emdata() //eman2
32 {
33 return this->getVolumeData()->get_emdata();
34 }
35
36 int Volume::getSizeX() {
37 return volData->GetSizeX();
38 }
39
40 int Volume::getSizeY() {
41 return volData->GetSizeY();
42 }
43
44 int Volume::getSizeZ() {
45 return volData->GetSizeZ();
46 }
47
48 int Volume::getIndex(int x, int y, int z) {
49 return volData->GetIndex(x, y, z);
50 }
51
52 void Volume::setDataAt( int x, int y, int z, double d ) {
53 volData->SetDataAt(x, y, z, (float)d);
54 }
55
56 void Volume::setDataAt( int index, double d ) {
57 volData->SetDataAt(index, (float)d);
58 }
59
60 double Volume::getDataAt( int x, int y, int z )
61 {
62 return volData->GetDataAt(x, y, z);
63 }
64
65 double Volume::getDataAt( int index ) {
66 return volData->GetDataAt(index);
67 }
68
69 VolumeData * Volume::getVolumeData() {
70 return volData;
71 }
72
73 void Volume::setSpacing(float spx, float spy, float spz ) {
74 volData->SetSpacing(spx, spy, spz);
75 }
76
77 void Volume::setOrigin(float orgX, float orgY, float orgZ) {
78 volData->SetOrigin(orgX, orgY, orgZ);
79 }
80
81 float Volume::getSpacingX() {
82 return volData->GetSpacingX();
83 }
84
85 float Volume::getSpacingY() {
86 return volData->GetSpacingY();
87 }
88
89 float Volume::getSpacingZ() {
90 return volData->GetSpacingZ();
91 }
92
93 float Volume::getOriginX() {
94 return volData->GetOriginX();
95 }
96
97 float Volume::getOriginY() {
98 return volData->GetOriginY();
99 }
100
101 float Volume::getOriginZ() {
102 return volData->GetOriginZ();
103 }
104
105 //Volume * Volume::getPseudoDensity( ) {
108 //int i, j, k ;
109 //Volume * res = new Volume(getSizeX(), getSizeY(), getSizeZ(), 0, 0, 0, this);
110 //int size = getSizeX() * getSizeY() * getSizeZ() ;
111 //srand(123) ;
112
113 //for ( i = 0 ; i < getSizeX() ; i ++ )
114 //for ( j = 0 ; j < getSizeY() ; j ++ )
115 //for ( k = 0 ; k < getSizeZ() ; k ++ ) {
116 //if ( res->getDataAt( i, j, k ) > 0 ) {
117 //int ct = 0 ;
118 //for ( int m = 0 ; m < 6 ; m ++ ) {
119 //if ( res->getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) > 0 ) {
120 //ct ++ ;
121 //}
122 //}
123 //res->setDataAt( i,j,k, (k/(float)getSizeZ())*(k/(float)getSizeZ()) ) ;
124 //if ( ct > 2 ) {
126 //} else {
128 //}
129 //}
130 //}
131
133 //for ( i = 0 ; i < 20 ; i ++ )
134 //{
135 //printf("Smoothing round %d\n", i) ;
136 //res->smooth( 0.5f ) ;
137 //}
138 //*/
139
140 //Volume * tvol = new Volume( getSizeX(), getSizeY(), getSizeZ(), 0, 0, 0, res);
141 //float d, ad, ct, temp;
142 //for ( int it = 0 ; it < 3 ; it ++ )
143 //for ( i = 0 ; i < getSizeX() ; i ++ )
144 //for ( j = 0 ; j < getSizeY() ; j ++ )
145 //for ( k = 0 ; k < getSizeZ() ; k ++ ) {
146 //if ( (d = (float)tvol->getDataAt( i, j, k )) > 0 ) {
147 //ad = 0 ; ct = 0 ;
148 //for ( int m = 0 ; m < 6 ; m ++ ) {
149 //if ( (temp = (float)tvol->getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] )) > 0 ) {
150 //ad += temp;
151 //ct ++ ;
152 //}
153 //}
154 //if ( ct > 0 ) {
155 //res->setDataAt( i, j, k, ( d + ad/ct ) / 2 ) ;
156 //}
157 //}
158 //}
159
160 //delete tvol;
161 //tvol = new Volume( getSizeX(), getSizeY(), getSizeZ(), 0, 0, 0, res ) ;
162 //for ( i = 0 ; i < 40 ; i ++ )
163 //{
164 //printf("Smoothing round %d\n", i) ;
165 //res->smooth( 0.5f ) ;
166 //continue ;
167
176
177 //}
178
179
180 //return res ;
181 //}
182
183
184 //Volume * Volume::getDistanceField(int rad, float randf) {
188
190 //int i, j, k ;
191 //Volume * res = new Volume(getSizeX(), getSizeY(), getSizeZ(), 0, 0, 0, this);
192 //srand( 123 ) ;
193
194 //for ( i = 0 ; i < getSizeX() ; i ++ )
195 //for ( j = 0 ; j < getSizeY() ; j ++ )
196 //for ( k = 0 ; k < getSizeZ() ; k ++ ) {
197 //if ( getDataAt(i, j, k) > 0 ) {
198 //float mag = 1 + randf * (float) rand() / (float) RAND_MAX ;
199 //int lx = max(0,i-rad) ;
200 //int ly = max(0,j-rad) ;
201 //int lz = max(0,k-rad) ;
202 //int hx = min(getSizeX()-1,i+rad) ;
203 //int hy = min(getSizeY()-1,j+rad) ;
204 //int hz = min(getSizeZ()-1,k+rad) ;
205 //int x,y,z;
206 //for ( x = lx ; x <= hx ; x ++ )
207 //for ( y = ly ; y <= hy ; y ++ )
208 //for ( z = lz ; z <= hz ; z ++ ) {
209 //float val = 1 - (float) sqrt((double)((x-i)*(x-i) + (y-j)*(y-j) + (z-k)*(z-k))) / (float) rad ;
210 //val *= mag ;
211 //if ( res->getDataAt( x, y, z ) < val ) {
212 //res->setDataAt( x, y, z, val ) ;
213 //}
214 //}
215 //}
216 //}
217
219 //for ( i = 0 ; i < 2 ; i ++ )
220 //{
221 //printf("Smoothing round %d\n", i) ;
222 //res->smooth( 0.5f ) ;
223 //}
224
225
226 //return res ;
227 //}
228
229
230 //int Volume::getNonZeroVoxelCount() {
231 //int count = 0;
232 //for(int x = 0; x < getSizeX(); x++){
233 //for(int y = 0; y < getSizeY(); y++){
234 //for(int z = 0; z < getSizeZ(); z++){
235 //if(this->getDataAt(x, y, z) > 0.0) {
236 //count++;
237 //}
238 //}
239 //}
240 //}
241 //return count;
242 //}
243 //void Volume::print() {
244 //for(int x = 0; x < getSizeX(); x++) {
245 //printf("{ ");
246 //for(int y = 0; y < getSizeY(); y++) {
247 //printf("{ ");
248 //for(int z = 0; z < getSizeZ(); z++) {
249 //printf("%f, ", getDataAt(x, y, z));
250 //}
251 //printf("} ");
252 //}
253 //printf("} ");
254 //}
255 //printf("\n");
256 //}
257
258
259 //void Volume::subtract( Volume* vol ) {
260 //int i, j, k ;
261 //for ( i = 0 ; i < getSizeX() ; i ++ )
262 //for ( j = 0 ; j < getSizeY() ; j ++ )
263 //for ( k = 0 ; k < getSizeZ() ; k ++ ) {
264 //if ( getDataAt( i, j, k ) > 0 ) {
265 //if ( vol->getDataAt(i,j,k) > 0 ) {
266 //setDataAt( i, j, k, 0 ) ;
267 //}
268 //}
269 //}
270
271 //}
272
273 void Volume::pad (int padBy, double padValue) {
274 volData->Pad(padBy, padValue);
275 }
276
277 //void Volume::applyMask(Volume * maskVol, double maskValue, bool keepMaskValue) {
278 //for(int x = 0; x < maskVol->getSizeX(); x++) {
279 //for(int y = 0; y < maskVol->getSizeY(); y++) {
280 //for(int z = 0; z < maskVol->getSizeZ(); z++) {
281 //if(((maskVol->getDataAt(x, y, z) == maskValue) && !keepMaskValue) ||
282 //((maskVol->getDataAt(x, y, z) != maskValue) && keepMaskValue)) {
283 //setDataAt(x, y, z, 0);
284 //}
285 //}
286 //}
287 //}
288 //}
289
290 //double Volume::getMin() {
291 //int size = volData->GetMaxIndex();
292 //double rvalue = volData->GetDataAt(0);
293 //for (int i=1; i < size; i++) {
294 //float val = volData->GetDataAt(i);
295 //if ( rvalue > val) {
296 //rvalue = val;
297 //}
298 //}
299 //return rvalue;
300 //}
301
302 //double Volume::getMax() {
303 //int size = volData->GetMaxIndex();
304 //double rvalue = volData->GetDataAt(0);
305 //for (int i=1; i < size; i++) {
306 //float val = volData->GetDataAt(i);
307 //if ( rvalue < val) {
308 //rvalue = val;
309 //}
310 //}
311 //return rvalue ;
312 //}
313
314 //double Volume::getMaxValuePosition(int& maxX, int& maxY, int& maxZ) {
315 //double maxVal = getDataAt(0,0,0);
316 //maxX = 0; maxY = 0; maxZ = 0;
317 //double data;
318
319 //for(int x = 0; x < getSizeX(); x++) {
320 //for(int y = 0; y < getSizeY(); y++) {
321 //for(int z = 0; z < getSizeZ(); z++) {
322 //data = getDataAt(x, y, z);
323 //if(data > maxVal) {
324 //maxVal = data;
325 //maxX = x; maxY = y; maxZ = z;
326 //}
327 //}
328 //}
329 //}
330 //return maxVal;
331 //}
332
333 //double Volume::getLocalMax(int x, int y, int z, int radius) {
334 //double mx = getDataAt(x, y, z);
335 //for(int xx = x - radius; xx <= x + radius; xx++) {
336 //for(int yy = y - radius; yy <= y + radius; yy++) {
337 //for(int zz = z - radius; zz <= z + radius; zz++) {
338 //mx = max(mx, getDataAt(xx, yy, zz));
339 //}
340 //}
341 //}
342 //return mx;
343 //}
344
345 //double Volume::getLocalMin(int x, int y, int z, int radius) {
346 //double mn = getDataAt(x, y, z);
347 //for(int xx = x - radius; xx <= x + radius; xx++) {
348 //for(int yy = y - radius; yy <= y + radius; yy++) {
349 //for(int zz = z - radius; zz <= z + radius; zz++) {
350 //mn = min(mn, getDataAt(xx, yy, zz));
351 //}
352 //}
353 //}
354 //return mn;
355 //}
356
357 //void Volume::fill( double val )
358 //{
359 //for(int x = 0; x < getSizeX(); x++) {
360 //for(int y = 0; y < getSizeY(); y++) {
361 //for(int z = 0; z < getSizeZ(); z++) {
362 //setDataAt(x, y, z, val);
363 //}
364 //}
365 //}
366 //}
367
368 //int Volume::isBertrandBorder( int ox, int oy, int oz, int dir ) {
369 //int nx = ox + neighbor6[dir][0] ;
370 //int ny = oy + neighbor6[dir][1] ;
371 //int nz = oz + neighbor6[dir][2] ;
372 //if ( getDataAt( nx, ny, nz) < 0 ) {
373 //return 1 ;
374 //}
375 //return 0 ;
376 //}
377
378 //int Volume::isBertrandEndPoint( int ox, int oy, int oz ) {
379 //double vox[3][3][3] ;
380
381 //int i, j, k ;
382 //for ( i = -1 ; i < 2 ; i ++ )
383 //for ( j = -1 ; j < 2 ; j ++ )
384 //for ( k = -1 ; k < 2 ; k ++ ) {
385 //double tval = getDataAt( ox + i, oy + j, oz + k ) ;
386 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
387 //}
388
390 //int xx6 = 0 ;
391 //for ( i = 0 ; i < 6 ; i ++ ) {
392 //if ( vox[ neighbor6[i][0] + 1 ][ neighbor6[i][1] + 1 ][ neighbor6[i][2] + 1 ] >= 0 ) {
393 //xx6 ++ ;
394 //}
395 //}
396
398 //int a26 = 0, na26 = 0 ;
399 //for ( i = 0 ; i < 3 ; i += 2 )
400 //for ( j = 0 ; j < 3 ; j += 2 )
401 //for ( k = 0 ; k < 3 ; k += 2 ) {
402 //if ( vox[1][j][k] < 0 ) {
403 //continue ;
404 //}
405 //if ( vox[i][1][k] < 0 ) {
406 //continue ;
407 //}
408 //if ( vox[i][j][1] < 0 ) {
409 //continue ;
410 //}
411 //if ( vox[i][1][1] < 0 ) {
412 //continue ;
413 //}
414 //if ( vox[1][j][1] < 0 ) {
415 //continue ;
416 //}
417 //if ( vox[1][1][k] < 0 ) {
418 //continue ;
419 //}
420 //if ( vox[i][j][k] >= 0 ) {
421 //a26 ++ ;
422 //} else {
423 //na26 ++ ;
424 //}
425 //}
426
427 //if (( na26 == 0 ) && ( a26 != 0 ) && ( xx6 <= a26 + 2 )) {
428 //return 0 ;
429 //}
430 //if ( isFeatureFace(ox,oy,oz) ) {
432 //}
433 //return 1 ;
434 //}
435
436 //int Volume::isHelix( int ox, int oy, int oz ) {
437 //int cn = 12 ;
438 //int nx, ny, nz ;
439 //int i, j, k ;
440
441 //double vox[3][3][3] ;
442 //for ( i = -1 ; i < 2 ; i ++ )
443 //for ( j = -1 ; j < 2; j ++ )
444 //for ( k = -1 ; k < 2 ; k ++ ) {
445 //vox[i+1][j+1][k+1] = getDataAt( ox + i, oy + j, oz + k ) ;
446 //}
447
448 //for ( i = 0 ; i < 12 ; i ++ ) {
449 //for ( j = 0 ; j < 4 ; j ++ ) {
450 //nx = sheetNeighbor[i][j][0] + 1;
451 //ny = sheetNeighbor[i][j][1] + 1;
452 //nz = sheetNeighbor[i][j][2] + 1;
453
454 //if ( vox[nx][ny][nz] <= 0 ) {
455 //cn -- ;
456 //break ;
457 //}
458 //}
459 //}
460
461 //if ( cn >= 1 ) {
462 //return 0 ;
463 //} else {
464 //return 1 ;
465 //}
466 //}
467
468 //int Volume::isSheet( int ox, int oy, int oz ) {
469 //int cn = 12 ;
470 //int nx, ny, nz ;
471
472 //for ( int i = 0 ; i < 12 ; i ++ ) {
473 //for ( int j = 0 ; j < 4 ; j ++ ) {
474 //nx = ox + sheetNeighbor[i][j][0] ;
475 //ny = oy + sheetNeighbor[i][j][1] ;
476 //nz = oz + sheetNeighbor[i][j][2] ;
477
478 //if ( getDataAt( nx, ny, nz ) <= 0 ) {
479 //cn -- ;
480 //break ;
481 //}
482 //}
483 //}
484 //return ( cn >= 3 ) ;
485 //}
486
487 //Volume * Volume::getSheets( int minSize ) {
488 //int i, j, k ;
489
491 //printf("Initialize volume at %d %d %d\n", getSizeX(), getSizeY(), getSizeZ() ) ;
492 //Volume* svol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
493
495 //int sheets[MAX_SHEETS] ;
496 //for ( i = 0 ; i < MAX_SHEETS ; i ++ ) {
497 //sheets[ i ] = 0 ;
498 //}
499 //int totSheets = 1 ;
500
502 //printf("Start clustering...\n" ) ;
503 //int ox, oy, oz ;
504 //for ( i = 0 ; i < getSizeX() ; i ++ )
505 //for ( j = 0 ; j < getSizeY() ; j ++ )
506 //for ( k = 0 ; k < getSizeZ() ; k ++ ) {
507 //if ( getDataAt(i,j,k) <= 0 || svol->getDataAt(i,j,k) != 0 ) {
509 //continue ;
510 //}
511 //if ( ! isSheet( i, j, k ) ) {
513 //continue ;
514 //}
515
517 //int numNodes = 1 ;
518 //svol->setDataAt( i, j, k, totSheets ) ;
519 //GridQueue* queue = new GridQueue() ;
520 //queue->pushQueue( i, j, k ) ;
521 //while ( queue->popQueue(ox, oy, oz) ) {
523 //if ( isSheet( ox, oy, oz ) ) {
524 //for ( int m = 0 ; m < 6 ; m ++ ) {
525 //int nx = ox + neighbor6[m][0] ;
526 //int ny = oy + neighbor6[m][1] ;
527 //int nz = oz + neighbor6[m][2] ;
528
529 //if ( getDataAt(nx,ny,nz) > 0 && svol->getDataAt(nx,ny,nz) == 0 ) {
530 //svol->setDataAt(nx,ny,nz,totSheets);
531 //queue->pushQueue(nx,ny,nz) ;
532 //numNodes ++ ;
533 //}
534 //}
535 //}
536 //}
537
538 //delete queue ;
539 //if ( numNodes > 0 ) {
541 //sheets[ totSheets ] = numNodes ;
542 //totSheets ++ ;
543 //}
544 //}
545
547 //printf("Removing small clusters.\n") ;
548 //for ( i = 0 ; i < getSizeX() ; i ++ )
549 //for ( j = 0 ; j < getSizeY() ; j ++ )
550 //for ( k = 0 ; k < getSizeZ() ; k ++ ) {
551 //int cnt = (int) svol->getDataAt(i,j,k) ;
552 //if ( cnt > 0 && sheets[ cnt ] < minSize ) {
553 //svol->setDataAt(i,j,k,-1) ;
554 //}
555 //}
556
558 //#ifdef VERBOSE
559 //printf("Thresholding the volume to 0/1...\n") ;
560 //#endif
561 //svol->threshold( 0.1, 0, 1 ) ;
562
563 //return svol ;
564 //}
565
566 //Volume * Volume::getHelices( int minSize ) {
567 //printf("Segmenting helices from eroded volume.\n") ;
568 //int i, j, k ;
569
571 //printf("Initialize volume at %d %d %d\n", getSizeX(), getSizeY(), getSizeZ() ) ;
572 //Volume* svol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
573
575 //int helices[MAX_SHEETS] ;
576 //for ( i = 0 ; i < MAX_SHEETS ; i ++ ) {
577 //helices[ i ] = 0 ;
578 //}
579 //int totHelices = 1 ;
580
582 //printf("Start clustering...\n" ) ;
583 //int ox, oy, oz ;
584 //for ( i = 0 ; i < getSizeX() ; i ++ )
585 //for ( j = 0 ; j < getSizeY() ; j ++ )
586 //for ( k = 0 ; k < getSizeZ() ; k ++ ) {
587 //if ( getDataAt(i,j,k) <= 0 || svol->getDataAt(i,j,k) != 0 ) {
589 //continue ;
590 //}
591 //if ( ! isHelix( i, j, k ) ) {
593 //continue ;
594 //}
595
597 //int numNodes = 1 ;
598 //svol->setDataAt( i, j, k, totHelices ) ;
599 //GridQueue* queue = new GridQueue() ;
600 //queue->pushQueue( i, j, k ) ;
601 //while ( queue->popQueue(ox, oy, oz) )
602 //{
604 //if ( isHelix( ox, oy, oz ) ) {
605 //for ( int m = 0 ; m < 6 ; m ++ ) {
606 //int nx = ox + neighbor6[m][0] ;
607 //int ny = oy + neighbor6[m][1] ;
608 //int nz = oz + neighbor6[m][2] ;
609
610 //if ( getDataAt(nx,ny,nz) > 0 && svol->getDataAt(nx,ny,nz) == 0 ) {
611 //svol->setDataAt(nx,ny,nz,totHelices);
612 //queue->pushQueue(nx,ny,nz) ;
613 //numNodes ++ ;
614 //}
615 //}
616 //}
617 //}
618
619 //delete queue ;
620 //if ( numNodes > 0 ) {
622 //helices[ totHelices ] = numNodes ;
623 //totHelices ++ ;
624 //}
625 //}
626
628 //printf("Removing small clusters.\n") ;
629 //for ( i = 0 ; i < getSizeX() ; i ++ )
630 //for ( j = 0 ; j < getSizeY() ; j ++ )
631 //for ( k = 0 ; k < getSizeZ() ; k ++ ) {
632 //int cnt = (int) svol->getDataAt(i,j,k) ;
633 //if ( cnt > 0 && helices[ cnt ] < minSize ) {
634 //svol->setDataAt(i,j,k,-1) ;
635 //}
636 //}
637
639 //#ifdef VERBOSE
640 //printf("Thresholding the volume to 0/1...\n") ;
641 //#endif
642 //svol->threshold( 0.1, 0, 1 ) ;
643
644 //return svol ;
645
646 //}
647
648 //int Volume::isEndPoint( int ox, int oy, int oz ) {
649 //if ( getDataAt( ox - 1, oy, oz ) < 0 && getDataAt( ox + 1, oy, oz ) < 0 ) {
650 //return 1 ;
651 //}
652 //if ( getDataAt( ox, oy - 1, oz ) < 0 && getDataAt( ox, oy + 1, oz ) < 0 ) {
653 //return 1 ;
654 //}
655 //if ( getDataAt( ox, oy, oz - 1 ) < 0 && getDataAt( ox, oy, oz + 1 ) < 0 ) {
656 //return 1 ;
657 //}
658 //return 0 ;
659 //}
660
661 int Volume::getNumNeighbor6( int ox, int oy, int oz ) {
662 int rvalue = 0 ;
663 for ( int i = 0 ; i < 6 ; i ++ ) {
664 int nx = ox + neighbor6[i][0] ;
665 int ny = oy + neighbor6[i][1] ;
666 int nz = oz + neighbor6[i][2] ;
667 if ( getDataAt( nx, ny, nz ) >= 0 ) {
668 rvalue ++ ;
669 }
670 }
671
672 return rvalue ;
673 }
674 //int Volume::testIsSheetEnd( int ox, int oy, int oz ) {
676 //int i, j ;
677 //int nx, ny, nz ;
678
679 //int edge[6] = { 0,0,0,0,0,0 } ;
680 //int faceflag[ 12 ] ;
681 //int hasNoiseFace = 0 ;
682 //int tot = 0 ;
683
684 //for ( i = 0 ; i < 12 ; i ++ ) {
685 //faceflag[ i ] = 1 ;
686 //int hasNoise = 0 ;
687
688 //for ( j = 0 ; j < 4 ; j ++ ) {
689 //nx = ox + sheetNeighbor[i][j][0] ;
690 //ny = oy + sheetNeighbor[i][j][1] ;
691 //nz = oz + sheetNeighbor[i][j][2] ;
692
693 //if ( getDataAt( nx, ny, nz ) == 0 ) {
694 //hasNoise = 1 ;
695 //}
696 //if ( getDataAt( nx, ny, nz ) < 0 ) {
697 //faceflag[ i ] = 0 ;
698 //break ;
699 //}
700 //}
701 //if ( faceflag[ i ] == 1 && hasNoise ) {
702 //hasNoiseFace ++ ;
704 //}
705
706 //if ( faceflag[ i ] ) {
707 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
708 //edge[ e0 ] ++ ;
709 //edge[ e1 ] ++ ;
710 //tot ++ ;
711 //}
712 //}
713
714 //if ( hasNoiseFace == tot ) {
715 //return 0 ;
716 //}
717
718 //if ( tot == 0 ) {
719 //return 0 ;
720 //}
721
723 //int numones = 0 ;
724 //for ( i = 0 ; i < 6 ; i ++ ) {
725 //if ( edge[ i ] == 1 ) {
726 //numones ++ ;
727 //}
728 //}
729 //while( numones > 0 ) {
730 //int e ;
731 //for ( i = 0 ; i < 6 ; i ++ ) {
732 //if ( edge[ i ] == 1 ) {
733 //e = i ;
734 //break ;
735 //}
736 //}
737
738 //int f, e2 ;
739 //for ( j = 0 ; j < 4 ; j ++ ) {
740 //f = edgeFaces[ e ][ j ] ;
741 //if ( faceflag[ f ] ) {
742 //break ;
743 //}
744 //}
745
746 //if ( faceEdges[ f ][ 0 ] == e ) {
747 //e2 = faceEdges[ f ][ 1 ] ;
748 //} else {
749 //e2 = faceEdges[ f ][ 0 ] ;
750 //}
751
752 //edge[ e ] -- ;
753 //numones -- ;
754 //edge[ e2 ] -- ;
755 //faceflag[ f ] = 0 ;
756 //tot -- ;
757
758 //if ( edge[ e2 ] == 1 ) {
759 //numones ++ ;
760 //} else if ( edge[ e2 ] == 0 ) {
761 //numones -- ;
762 //}
763 //}
764
765 //if ( tot > 0 ) {
766 //return 0 ;
767 //}
768 //return 1 ;
769 //}
770
771 //int Volume::isNoiseSheetEnd( int ox, int oy, int oz ) {
772 //int i, j ;
773 //int nx, ny, nz ;
774
775 //int edge[6] = { 0,0,0,0,0,0 } ;
776 //int faceflag[ 12 ] ;
777 //int hasNoiseFace = 0 ;
778 //int tot = 0 ;
779
780 //for ( i = 0 ; i < 12 ; i ++ ) {
781 //faceflag[ i ] = 1 ;
782 //int hasNoise = 0 ;
783
784 //for ( j = 0 ; j < 4 ; j ++ ) {
785 //nx = ox + sheetNeighbor[i][j][0] ;
786 //ny = oy + sheetNeighbor[i][j][1] ;
787 //nz = oz + sheetNeighbor[i][j][2] ;
788
789 //if ( getDataAt( nx, ny, nz ) == 0 ) {
790 //hasNoise = 1 ;
791 //}
792 //if ( getDataAt( nx, ny, nz ) < 0 ) {
793 //faceflag[ i ] = 0 ;
794 //break ;
795 //}
796 //}
797 //if ( faceflag[ i ] == 1 && hasNoise ) {
798 //hasNoiseFace ++ ;
800 //}
801
802 //if ( faceflag[ i ] ) {
803 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
804 //edge[ e0 ] ++ ;
805 //edge[ e1 ] ++ ;
806 //tot ++ ;
807 //}
808 //}
809
810 //if ( hasNoiseFace < tot ) {
811 //return 0 ;
812 //}
813
814 //if ( tot == 0 ) {
815 //return 0 ;
816 //}
817
819 //int numones = 0 ;
820 //for ( i = 0 ; i < 6 ; i ++ ) {
821 //if ( edge[ i ] == 1 )
822 //{
823 //numones ++ ;
824 //}
825 //}
826 //while( numones > 0 ) {
827 //int e ;
828 //for ( i = 0 ; i < 6 ; i ++ ) {
829 //if ( edge[ i ] == 1 ) {
830 //e = i ;
831 //break ;
832 //}
833 //}
834
835 //int f, e2 ;
836 //for ( j = 0 ; j < 4 ; j ++ ) {
837 //f = edgeFaces[ e ][ j ] ;
838 //if ( faceflag[ f ] ) {
839 //break ;
840 //}
841 //}
842
843 //if ( faceEdges[ f ][ 0 ] == e ) {
844 //e2 = faceEdges[ f ][ 1 ] ;
845 //} else {
846 //e2 = faceEdges[ f ][ 0 ] ;
847 //}
848
849 //edge[ e ] -- ;
850 //numones -- ;
851 //edge[ e2 ] -- ;
852 //faceflag[ f ] = 0 ;
853 //tot -- ;
854
855 //if ( edge[ e2 ] == 1 ) {
856 //numones ++ ;
857 //} else if ( edge[ e2 ] == 0 ) {
858 //numones -- ;
859 //}
860 //}
861
862 //if ( tot > 0 ) {
863 //return 0 ;
864 //}
865 //return 1 ;
866 //}
867
868 //int Volume::isInternal( int ox, int oy, int oz ) {
870 //int i, j, k ;
871
872 //for ( i = -1 ; i < 2 ; i ++ )
873 //for ( j = -1 ; j < 2 ; j ++ )
874 //for ( k = -1 ; k < 2 ; k ++ ) {
875 //if ( getDataAt( ox + i, oy + j, oz + k ) <= 0 ) {
876 //return 0 ;
877 //}
878 //}
879 //return 1 ;
880 //}
881
882 //int Volume::isInternal2( int ox, int oy, int oz ) {
884 //int i, j, k ;
885
886 //for ( i = -1 ; i < 2 ; i ++ )
887 //for ( j = -1 ; j < 2 ; j ++ )
888 //for ( k = -1 ; k < 2 ; k ++ ) {
889 //if ( getDataAt( ox + i, oy + j, oz + k ) < 0 ) {
890 //return 0 ;
891 //}
892 //}
893
894 //return 1 ;
895 //}
896
897 //int Volume::hasIsolatedFace( int ox, int oy, int oz ) {
898 //int i, j, k ;
899 //int nx, ny, nz ;
900
901 //double vox[3][3][3] ;
902 //for ( i = -1 ; i < 2 ; i ++ )
903 //for ( j = -1 ; j < 2 ; j ++ )
904 //for ( k = -1 ; k < 2 ; k ++ ) {
905 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ;
906 //}
907
908 //int cells[8] = { 1, 1, 1, 1, 1, 1, 1, 1 } ;
909 //for ( i = 0 ; i < 8 ; i ++ ) {
910 //int x = ( ( i >> 2 ) & 1 ) ;
911 //int y = ( ( i >> 1 ) & 1 ) ;
912 //int z = ( i & 1 ) ;
913 //for ( j = 0 ; j < 8 ; j ++ ) {
914 //nx = x + ( ( j >> 2 ) & 1 ) ;
915 //ny = y + ( ( j >> 1 ) & 1 ) ;
916 //nz = z + ( j & 1 ) ;
917
918 //if ( vox[nx][ny][nz] < 0 ) {
919 //cells[i] = 0 ;
920 //break;
921 //}
922 //}
923 //}
924
925 //for ( i = 0 ; i < 12 ; i ++ ) {
926 //if ( cells[ faceCells[i][0] ] == 1 || cells[ faceCells[i][1] ] == 1 ) {
927 //continue ;
928 //}
929 //int flag = 1 ;
930 //for ( j = 0 ; j < 4 ; j ++ ) {
931 //nx = 1 + sheetNeighbor[i][j][0] ;
932 //ny = 1 + sheetNeighbor[i][j][1] ;
933 //nz = 1 + sheetNeighbor[i][j][2] ;
934
935 //if ( vox[nx][ny][nz] < 0 ) {
936 //flag = 0 ;
937 //break;
938 //}
939 //}
940 //if ( flag ) {
941 //return 1 ;
942 //}
943 //}
944
945 //return 0 ;
946 //}
947
948 //int Volume::hasIsolatedEdge( int ox, int oy, int oz ) {
949 //int i, j, k ;
950 //int nx, ny, nz ;
951
952 //double vox[3][3][3] ;
953 //for ( i = -1 ; i < 2 ; i ++ )
954 //for ( j = -1 ; j < 2 ; j ++ )
955 //for ( k = -1 ; k < 2 ; k ++ ) {
956 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ;
957 //}
958
959 //int edge[6] = { 0,0,0,0,0,0 } ;
960
961 //for ( i = 0 ; i < 12 ; i ++ ) {
962 //int flag = 1 ;
963 //for ( j = 0 ; j < 4 ; j ++ ) {
964 //nx = 1 + sheetNeighbor[i][j][0] ;
965 //ny = 1 + sheetNeighbor[i][j][1] ;
966 //nz = 1 + sheetNeighbor[i][j][2] ;
967
968 //if ( vox[nx][ny][nz] < 0 ) {
969 //flag = 0 ;
970 //break ;
971 //}
972 //}
973
974 //if ( flag ) {
975 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
976 //edge[ e0 ] ++ ;
977 //edge[ e1 ] ++ ;
978 //}
979 //}
980
981 //for ( i = 0 ; i < 6 ; i ++ ) {
982 //if ( edge[i] ) {
983 //continue ;
984 //}
985
986 //nx = 1 + neighbor6[i][0] ;
987 //ny = 1 + neighbor6[i][1] ;
988 //nz = 1 + neighbor6[i][2] ;
989
990 //if ( vox[nx][ny][nz] >= 0 ) {
991 //return 1 ;
992 //}
993
994 //}
995
996 //return 0 ;
997 //}
998
999 //int Volume::countFace( int ox, int oy, int oz, int m ) {
1000 //int facenum = 4 ;
1001 //for ( int i = 0 ; i < 4 ; i ++ ) {
1002 //for ( int j = 0 ; j < 4 ; j ++ ) {
1003 //int nx = ox + sheetNeighbor[edgeFaces[m][i]][j][0] ;
1004 //int ny = oy + sheetNeighbor[edgeFaces[m][i]][j][1] ;
1005 //int nz = oz + sheetNeighbor[edgeFaces[m][i]][j][2] ;
1006
1007 //if ( getDataAt( nx, ny, nz ) < 0 ) {
1008 //facenum -- ;
1009 //break ;
1010 //}
1011 //}
1012 //}
1013
1014 //return facenum;
1015 //}
1016 int Volume::hasCell( int ox, int oy, int oz ) {
1017 for ( int i = 0 ; i < 2 ; i ++ )
1018 for ( int j = 0 ; j < 2 ; j ++ )
1019 for ( int k = 0 ; k < 2 ; k ++ ) {
1020 if ( getDataAt( ox + i, oy + j, oz + k ) < 0 ) {
1021 return 0 ;
1022 }
1023 }
1024 return 1 ;
1025 }
1026
1027 Volume * Volume::markCellFace( ) {
1028 int i, j, k ;
1029 Volume* fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
1030
1031 //return fvol ;
1032
1033 for ( i = 0 ; i < getSizeX() ; i ++ )
1034 for ( j = 0 ; j < getSizeY() ; j ++ )
1035 for ( k = 0 ; k < getSizeZ() ; k ++ )
1036 {
1037 if( getDataAt(i,j,k) >= 0 )
1038 {
1039 if ( hasCell(i,j,k) )
1040 {
1041 for ( int m = 0 ; m < 6 ; m ++ )
1042 {
1043 int nx = i + neighbor6[m][0] ;
1044 int ny = j + neighbor6[m][1] ;
1045 int nz = k + neighbor6[m][2] ;
1046 if ( ! hasCell(nx,ny,nz) )
1047 {
1048 fvol->setDataAt(i,j,k,(double)(1<<m)) ;
1049 break ;
1050 }
1051 }
1052 }
1053 }
1054 }
1055
1056
1057 return fvol ;
1058 }
1059
1060 //Volume * Volume::markFaceEdge( ) {
1061 //int x,y,z, i,j ;
1062 //Volume* fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
1063
1065
1066 //for ( x = 0 ; x < getSizeX() ; x ++ )
1067 //for ( y = 0 ; y < getSizeY() ; y ++ )
1068 //for ( z = 0 ; z < getSizeZ() ; z ++ )
1069 //{
1070 //if( getDataAt(x,y,z) >= 0 )
1071 //{
1072
1073 //for ( i = 0 ; i < 3 ; i ++ )
1074 //{
1075 //int hasFace = 1 ;
1076 //for ( j = 0 ; j < 4 ; j ++ )
1077 //{
1078 //int nx = x + sheetNeighbor[4 * i + 3][j][0] ;
1079 //int ny = y + sheetNeighbor[4 * i + 3][j][1] ;
1080 //int nz = z + sheetNeighbor[4 * i + 3][j][2] ;
1081
1082 //if ( getDataAt( nx, ny, nz ) < 0 )
1083 //{
1084 //hasFace = 0 ;
1085 //break ;
1086 //}
1087 //}
1088
1089 //if ( hasFace )
1090 //{
1092 //switch( i )
1093 //{
1094 //case 0:
1095 //if ( countFace( x, y, z, 0 ) == 1 )
1096 //{
1097 //fvol->setDataAt(x, y, z, (double)(1<<0)) ;
1098 //break ;
1099 //}
1100 //if ( countFace( x, y, z, 2 ) == 1 )
1101 //{
1102 //fvol->setDataAt(x, y, z, (double)(1<<2)) ;
1103 //break ;
1104 //}
1105 //if ( countFace( x, y + 1, z, 0 ) == 1 )
1106 //{
1107 //fvol->setDataAt(x, y + 1, z, (double)(1<<0)) ;
1108 //break ;
1109 //}
1110 //if ( countFace( x, y, z + 1, 2 ) == 1 )
1111 //{
1112 //fvol->setDataAt(x, y, z + 1, (double)(1<<2)) ;
1113 //break ;
1114 //}
1115 //printf("Hmmm... a face with no open edges.\n");
1116 //break ;
1117 //case 1:
1118 //if ( countFace( x, y, z, 0 ) == 1 )
1119 //{
1120 //fvol->setDataAt(x, y, z, (double)(1<<0)) ;
1121 //break ;
1122 //}
1123 //if ( countFace( x, y, z, 4 ) == 1 )
1124 //{
1125 //fvol->setDataAt(x, y, z, (double)(1<<4)) ;
1126 //break ;
1127 //}
1128 //if ( countFace( x + 1, y, z, 0 ) == 1 )
1129 //{
1130 //fvol->setDataAt(x + 1, y, z, (double)(1<<0)) ;
1131 //break ;
1132 //}
1133 //if ( countFace( x, y, z + 1, 4 ) == 1 )
1134 //{
1135 //fvol->setDataAt(x, y, z + 1, (double)(1<<4)) ;
1136 //break ;
1137 //}
1138 //printf("Hmmm... a face with no open edges.\n");
1139 //break ;
1140 //case 2:
1141 //if ( countFace( x, y, z, 2 ) == 1 )
1142 //{
1143 //fvol->setDataAt(x, y, z, (double)(1<<2)) ;
1144 //break ;
1145 //}
1146 //if ( countFace( x, y, z, 4 ) == 1 )
1147 //{
1148 //fvol->setDataAt(x, y, z, (double)(1<<4)) ;
1149 //break ;
1150 //}
1151 //if ( countFace( x + 1, y, z, 2 ) == 1 )
1152 //{
1153 //fvol->setDataAt(x + 1, y, z, (double)(1<<2)) ;
1154 //break ;
1155 //}
1156 //if ( countFace( x, y + 1, z, 4 ) == 1 )
1157 //{
1158 //fvol->setDataAt(x, y + 1, z, (double)(1<<4)) ;
1159 //break ;
1160 //}
1161 //printf("Hmmm... a face with no open edges.\n");
1162 //break ;
1163 //}
1164 //}
1165 //}
1166
1167 //}
1168 //}
1169
1170
1171 //return fvol ;
1172 //}
1173
1174 int Volume::hasCompleteSheet( int ox, int oy, int oz, Volume* fvol ) {
1175 int i, j, k ;
1176 int nx, ny, nz ;
1177
1178 int edge[6] = { 0,0,0,0,0,0 } ;
1179 int faceflag[ 12 ] ;
1180 int tot = 0 ;
1181 int cellflag[ 8 ] ;
1182
1183 int ct = 0 ;
1184 for ( i = -1 ; i < 1 ; i ++ )
1185 for ( j = -1 ; j < 1 ; j ++ )
1186 for ( k = -1 ; k < 1 ; k ++ )
1187 {
1188 if ( hasCell( ox + i, oy + j, oz + k ) )
1189 {
1190 cellflag[ ct ] = 1 ;
1191 }
1192 else
1193 {
1194 cellflag[ ct ] = 0 ;
1195 }
1196 ct ++ ;
1197 }
1198
1199 for ( i = 0 ; i < 12 ; i ++ )
1200 {
1201 faceflag[ i ] = 1 ;
1202 for ( j = 0 ; j < 4 ; j ++ )
1203 {
1204 nx = ox + sheetNeighbor[i][j][0] ;
1205 ny = oy + sheetNeighbor[i][j][1] ;
1206 nz = oz + sheetNeighbor[i][j][2] ;
1207
1208 if ( getDataAt( nx, ny, nz ) < 0 )
1209 {
1210 faceflag[ i ] = 0 ;
1211 break ;
1212 }
1213 }
1214
1215 if ( faceflag[ i ] )
1216 {
1217 if ( cellflag[ faceCells[i][0] ] ^ cellflag[ faceCells[i][1] ] )
1218 {
1219 int v1 = (int)( fvol->getDataAt(
1220 ox - 1 + (( faceCells[i][0] >> 2 ) & 1 ),
1221 oy - 1 + (( faceCells[i][0] >> 1 ) & 1 ),
1222 oz - 1 + (( faceCells[i][0] ) & 1)) ) ;
1223 int v2 = (int)( fvol->getDataAt(
1224 ox - 1 + (( faceCells[i][1] >> 2 ) & 1 ),
1225 oy - 1 + (( faceCells[i][1] >> 1 ) & 1 ),
1226 oz - 1 + (( faceCells[i][1] ) & 1)) ) ;
1227 if ( ((v1 >> (2 * (2 - i/4))) & 1) ||
1228 ((v2 >> (2 * (2 - i/4) + 1 )) & 1) )
1229 {
1230 faceflag[ i ] = 0 ;
1231 }
1232 }
1233 }
1234
1235 if ( faceflag[ i ] )
1236 {
1237 int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
1238 edge[ e0 ] ++ ;
1239 edge[ e1 ] ++ ;
1240 tot ++ ;
1241 }
1242 }
1243
1244 // Removing 1s
1245 int numones = 0 ;
1246 for ( i = 0 ; i < 6 ; i ++ )
1247 {
1248 if ( edge[ i ] == 1 )
1249 {
1250 numones ++ ;
1251 }
1252 }
1253 while( numones > 0 )
1254 {
1255 int e = 0;
1256 for ( i = 0 ; i < 6 ; i ++ )
1257 {
1258 if ( edge[ i ] == 1 )
1259 {
1260 e = i ;
1261 break ;
1262 }
1263 }
1264 /*
1265 if ( edge[ e ] != 1 )
1266 {
1267 printf("Wrong Again!********\n") ;
1268 }
1269 */
1270
1271 int f, e2 ;
1272 for ( j = 0 ; j < 4 ; j ++ )
1273 {
1274 f = edgeFaces[ e ][ j ] ;
1275 if ( faceflag[ f ] )
1276 {
1277 break ;
1278 }
1279 }
1280
1281 /*
1282 if ( faceflag[ f ] == 0 )
1283 {
1284 printf("Wrong!********\n") ;
1285 }
1286 */
1287
1288 if ( faceEdges[ f ][ 0 ] == e )
1289 {
1290 e2 = faceEdges[ f ][ 1 ] ;
1291 }
1292 else
1293 {
1294 e2 = faceEdges[ f ][ 0 ] ;
1295 }
1296
1297
1298 edge[ e ] -- ;
1299 numones -- ;
1300 edge[ e2 ] -- ;
1301 faceflag[ f ] = 0 ;
1302 tot -- ;
1303
1304 if ( edge[ e2 ] == 1 )
1305 {
1306 numones ++ ;
1307 }
1308 else if ( edge[ e2 ] == 0 )
1309 {
1310 numones -- ;
1311 }
1312 }
1313
1314 if ( tot > 0 )
1315 {
1316 return 1 ;
1317 }
1318
1319 return 0 ;
1320 }
1321
1322 int Volume::hasCompleteSheet( int ox, int oy, int oz ) {
1323 // Returns 1 if it lies in the middle of a sheet
1324 int temp = countIntEuler( ox, oy, oz ) ;
1325 if ( temp > 0 )
1326 {
1327 return 1 ;
1328 }
1329 else
1330 {
1331 return 0 ;
1332 }
1333 }
1334
1335 //int Volume::hasCompleteSheetSlow( int ox, int oy, int oz ) {
1336 //int i, j ;
1337 //int nx, ny, nz ;
1338
1339 //int edge[6] = { 0,0,0,0,0,0 } ;
1340 //int faceflag[ 12 ] ;
1341 //int tot = 0 ;
1342
1343 //for ( i = 0 ; i < 12 ; i ++ )
1344 //{
1345 //faceflag[ i ] = 1 ;
1346 //for ( j = 0 ; j < 4 ; j ++ )
1347 //{
1348 //nx = ox + sheetNeighbor[i][j][0] ;
1349 //ny = oy + sheetNeighbor[i][j][1] ;
1350 //nz = oz + sheetNeighbor[i][j][2] ;
1351
1352 //if ( getDataAt( nx, ny, nz ) < 0 )
1353 //{
1354 //faceflag[ i ] = 0 ;
1355 //break ;
1356 //}
1357 //}
1358
1359 //if ( faceflag[ i ] )
1360 //{
1361 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
1362 //edge[ e0 ] ++ ;
1363 //edge[ e1 ] ++ ;
1364 //tot ++ ;
1365 //}
1366 //}
1367
1369 //int numones = 0 ;
1370 //for ( i = 0 ; i < 6 ; i ++ )
1371 //{
1372 //if ( edge[ i ] == 1 )
1373 //{
1374 //numones ++ ;
1375 //}
1376 //}
1377 //while( numones > 0 )
1378 //{
1379 //int e ;
1380 //for ( i = 0 ; i < 6 ; i ++ )
1381 //{
1382 //if ( edge[ i ] == 1 )
1383 //{
1384 //e = i ;
1385 //break ;
1386 //}
1387 //}
1389 //if ( edge[ e ] != 1 )
1390 //{
1391 //printf("Wrong Again!********\n") ;
1392 //}
1393 //*/
1394
1395 //int f, e2 ;
1396 //for ( j = 0 ; j < 4 ; j ++ )
1397 //{
1398 //f = edgeFaces[ e ][ j ] ;
1399 //if ( faceflag[ f ] )
1400 //{
1401 //break ;
1402 //}
1403 //}
1404
1406 //if ( faceflag[ f ] == 0 )
1407 //{
1408 //printf("Wrong!********\n") ;
1409 //}
1410 //*/
1411
1412 //if ( faceEdges[ f ][ 0 ] == e )
1413 //{
1414 //e2 = faceEdges[ f ][ 1 ] ;
1415 //}
1416 //else
1417 //{
1418 //e2 = faceEdges[ f ][ 0 ] ;
1419 //}
1420
1421
1422 //edge[ e ] -- ;
1423 //numones -- ;
1424 //edge[ e2 ] -- ;
1425 //faceflag[ f ] = 0 ;
1426 //tot -- ;
1427
1428 //if ( edge[ e2 ] == 1 )
1429 //{
1430 //numones ++ ;
1431 //}
1432 //else if ( edge[ e2 ] == 0 )
1433 //{
1434 //numones -- ;
1435 //}
1436 //}
1437
1439 //if ( tot > 0 )
1440 //{
1442 //{
1444 //}
1445 //return 1 ;
1446 //}
1447 //else
1448 //{
1450 //{
1452 //}
1453 //return 0 ;
1454 //}
1455 //}
1456 int Volume::hasCompleteHelix( int ox, int oy, int oz )
1457 {
1458 // Returns 1 if it has a complete helix
1459 int i ;
1460 int c1 = 0;
1461 int nx, ny, nz ;
1462 int j ;
1463
1464 for ( i = 0 ; i < 6 ; i ++ )
1465 {
1466 nx = ox + neighbor6[i][0] ;
1467 ny = oy + neighbor6[i][1] ;
1468 nz = oz + neighbor6[i][2] ;
1469 if ( getDataAt( nx, ny, nz ) >= 0 )
1470 {
1471 c1 ++ ;
1472 j = i ;
1473 }
1474
1475 }
1476
1477 if ( c1 > 1 ) // || c1 == 0 )
1478 {
1479 return 1 ;
1480 }
1481
1482 return 0 ;
1483
1484 /*
1485 ox = ox + neighbor6[j][0] ;
1486 oy = oy + neighbor6[j][1] ;
1487 oz = oz + neighbor6[j][2] ;
1488 c1 = 0 ;
1489 for ( i = 0 ; i < 6 ; i ++ )
1490 {
1491 nx = ox + neighbor6[i][0] ;
1492 ny = oy + neighbor6[i][1] ;
1493 nz = oz + neighbor6[i][2] ;
1494 if ( getDataAt( nx, ny, nz ) >= 0 )
1495 {
1496 c1 ++ ;
1497 }
1498
1499 }
1500
1501 if ( c1 > 1 )
1502 {
1503 return 0 ;
1504 }
1505 else
1506 {
1507 return 1 ;
1508 }
1509 */
1510 }
1511
1512 int Volume::hasCompleteHelix( int ox, int oy, int oz, Volume* fvol )
1513 {
1514
1515 int i ;
1516 int c1 = 0;
1517 int nx, ny, nz ;
1518 int j ;
1519
1520 for ( i = 0 ; i < 6 ; i ++ )
1521 {
1522 nx = ox + neighbor6[i][0] ;
1523 ny = oy + neighbor6[i][1] ;
1524 nz = oz + neighbor6[i][2] ;
1525 if ( getDataAt( nx, ny, nz ) >= 0 )
1526 {
1527 if ( i % 2 == 0 )
1528 {
1529 nx = ox ;
1530 ny = oy ;
1531 nz = oz ;
1532 }
1533
1534 int val = (int)fvol->getDataAt( nx, ny, nz) ;
1535 if ( (( val >> ( 2 * ( i / 2 ) ) ) & 1) == 0 )
1536 {
1537 c1 ++ ;
1538 j = i ;
1539 }
1540 }
1541
1542 }
1543
1544 if ( c1 > 1 )
1545 {
1546 return 1 ;
1547 }
1548
1549 return 0 ;
1550
1551 /*
1552 ox = ox + neighbor6[j][0] ;
1553 oy = oy + neighbor6[j][1] ;
1554 oz = oz + neighbor6[j][2] ;
1555 c1 = 0 ;
1556 for ( i = 0 ; i < 6 ; i ++ )
1557 {
1558 nx = ox + neighbor6[i][0] ;
1559 ny = oy + neighbor6[i][1] ;
1560 nz = oz + neighbor6[i][2] ;
1561 if ( getDataAt( nx, ny, nz ) >= 0 )
1562 {
1563 c1 ++ ;
1564 }
1565
1566 }
1567
1568 if ( c1 > 1 )
1569 {
1570 return 0 ;
1571 }
1572 else
1573 {
1574 return 1 ;
1575 }
1576 */
1577 }
1578
1579 int Volume::isHelixEnd( int ox, int oy, int oz, Volume* nvol ) {
1580 // Returns 1 if it is a curve endpoint
1581 int i ;
1582 int c1 = 0 , c2 = 0 ;
1583 int nx, ny, nz ;
1584
1585 for ( i = 0 ; i < 6 ; i ++ )
1586 {
1587 nx = ox + neighbor6[i][0] ;
1588 ny = oy + neighbor6[i][1] ;
1589 nz = oz + neighbor6[i][2] ;
1590
1591 double val = getDataAt( nx, ny, nz ) ;
1592
1593 if ( val >= 0 )
1594 {
1595 c1 ++ ;
1596 if ( val > 0 && val < MAX_ERODE && nvol->getDataAt( nx, ny, nz ) == 0 )
1597 {
1598 c2 ++ ;
1599 }
1600 }
1601
1602 }
1603
1604 if ( c1 == 1 && c2 == 1 )
1605 {
1606 return 1 ;
1607 }
1608
1609 return 0 ;
1610 }
1611
1612 //int Volume::isFeature18( int ox, int oy, int oz )
1613 //{
1617 //if ( getDataAt( ox, oy, oz ) <= 0 )
1618 //{
1619 //return 0 ;
1620 //}
1621
1622 //for ( int i = 0 ; i < 3 ; i ++ )
1623 //for ( int j = 0 ; j < 3 ; j ++ )
1624 //for ( int k = 0 ; k < 3 ; k ++ )
1625 //{
1626 //if ( i == 1 || j == 1 || k == 1 )
1627 //{
1628 //if( getDataAt( ox - 1 + i, oy - 1 + j, oz - 1 + k ) == 0 )
1629 //{
1630 //return 0 ;
1631 //}
1632 //}
1633 //}
1634 //return 1 ;
1635 //}
1636
1637 //int Volume::isEdgeEnd( int ox, int oy, int oz )
1638 //{
1639 //int i ;
1640 //int c1 = 0 ;
1641 //int nx, ny, nz ;
1642
1643 //for ( i = 0 ; i < 6 ; i ++ )
1644 //{
1645 //nx = ox + neighbor6[i][0] ;
1646 //ny = oy + neighbor6[i][1] ;
1647 //nz = oz + neighbor6[i][2] ;
1648
1649 //double val = getDataAt( nx, ny, nz ) ;
1650
1651 //if ( val >= 0 )
1652 //{
1653 //c1 ++ ;
1654 //}
1655
1656 //}
1657
1658 //if ( c1 == 1 )
1659 //{
1660 //return 1 ;
1661 //}
1662
1663 //return 0 ;
1664 //}
1665
1666 //int Volume::isFaceEnd( int ox, int oy, int oz )
1667 //{
1669
1670 //int i, j, k ;
1671 //int nx, ny, nz ;
1672
1673 //double vox[3][3][3] ;
1674 //for ( i = -1 ; i < 2 ; i ++ )
1675 //for ( j = -1 ; j < 2 ; j ++ )
1676 //for ( k = -1 ; k < 2 ; k ++ )
1677 //{
1678 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ;
1679 //}
1680
1681 //int edge[6] = { 4,4,4,4,4,4 } ;
1682 //int edge2[6] = { 4,4,4,4,4,4 } ;
1683
1684 //for ( i = 0 ; i < 12 ; i ++ )
1685 //{
1686 //for ( j = 0 ; j < 4 ; j ++ )
1687 //{
1688 //nx = 1 + sheetNeighbor[i][j][0] ;
1689 //ny = 1 + sheetNeighbor[i][j][1] ;
1690 //nz = 1 + sheetNeighbor[i][j][2] ;
1691
1692 //if ( vox[nx][ny][nz] < 0 )
1693 //{
1694 //edge[ faceEdges[ i ][ 0 ] ] -- ;
1695 //edge[ faceEdges[ i ][ 1 ] ] -- ;
1696 //break ;
1697 //}
1698 //}
1699 //for ( j = 0 ; j < 4 ; j ++ )
1700 //{
1701 //nx = 1 + sheetNeighbor[i][j][0] ;
1702 //ny = 1 + sheetNeighbor[i][j][1] ;
1703 //nz = 1 + sheetNeighbor[i][j][2] ;
1704
1705 //if ( vox[nx][ny][nz] < 2 )
1706 //{
1707 //edge2[ faceEdges[ i ][ 0 ] ] -- ;
1708 //edge2[ faceEdges[ i ][ 1 ] ] -- ;
1709 //break ;
1710 //}
1711 //}
1712 //}
1713
1714 //for ( i = 0 ; i < 6 ; i ++ )
1715 //{
1716 //if ( edge[ i ] == 1 && edge2[ i ] == 1 )
1717 //{
1718 //return 1 ;
1719 //}
1720 //}
1721
1722 //return 0 ;
1723 //}
1724
1725 //int Volume::isNoise( int ox, int oy, int oz, int noise )
1726 //{
1727 //if ( getDataAt(ox,oy,oz) == 1 )
1728 //{
1729 //return 1 ;
1730 //}
1731
1732 //if ( noise > 0 )
1733 //{
1734 //for ( int i = 0 ; i < 6 ; i ++ )
1735 //{
1736 //int nx = ox + neighbor6[i][0] ;
1737 //int ny = oy + neighbor6[i][1] ;
1738 //int nz = oz + neighbor6[i][2] ;
1739
1740 //if ( getDataAt( nx, ny, nz ) > 0 )
1741 //{
1742 //if ( isNoise( nx, ny, nz, noise - 1 ) )
1743 //{
1744 //return 1 ;
1745 //}
1746 //}
1747 //}
1748 //}
1749
1750 //return 0 ;
1751
1752 //}
1753
1754 //int Volume::isNoiseHelixEnd( int ox, int oy, int oz )
1755 //{
1756
1757 //int i ;
1758 //int c1 = 0 , c2 = 0 ;
1759 //int nx, ny, nz ;
1760
1761 //for ( i = 0 ; i < 6 ; i ++ )
1762 //{
1763 //nx = ox + neighbor6[i][0] ;
1764 //ny = oy + neighbor6[i][1] ;
1765 //nz = oz + neighbor6[i][2] ;
1766
1767 //double val = getDataAt( nx, ny, nz ) ;
1768
1769 //if ( val >= 0 )
1770 //{
1771 //c1 ++ ;
1772 //if ( val > 0 && val < MAX_ERODE )
1773 //{
1774 //c2 ++ ;
1775 //}
1776 //}
1777
1778 //}
1779
1780 //if ( c1 == 1 && c2 == 0 )
1781 //{
1782 //return 1 ;
1783 //}
1784
1785 //return 0 ;
1786 //}
1787
1788
1789 int Volume::isHelixEnd( int ox, int oy, int oz ) {
1790
1791 int i ;
1792 int c1 = 0 , c2 = 0 ;
1793 int nx, ny, nz ;
1794
1795 for ( i = 0 ; i < 6 ; i ++ )
1796 {
1797 nx = ox + neighbor6[i][0] ;
1798 ny = oy + neighbor6[i][1] ;
1799 nz = oz + neighbor6[i][2] ;
1800
1801 double val = getDataAt( nx, ny, nz ) ;
1802
1803 if ( val >= 0 )
1804 {
1805 c1 ++ ;
1806 if ( getNumNeighbor6(nx,ny,nz) < 6 ) // if ( val > 0 && val < MAX_ERODE )
1807 {
1808 c2 ++ ;
1809 }
1810 }
1811
1812 }
1813
1814 if ( c1 == 1 && c2 == 1 )
1815 {
1816 return 1 ;
1817 }
1818
1819 return 0 ;
1820 }
1821 int Volume::isSheetEnd( int ox, int oy, int oz, Volume* nvol )
1822 {
1823 // Returns 1 if it contains a sheet boundary. Noise-resistant
1824 int i, j ;
1825 int nx, ny, nz ;
1826
1827 int edge[6] = { 0,0,0,0,0,0 } ;
1828 int faceflag[ 12 ] ;
1829 int hasFeatureFace = 0 ;
1830 int tot = 0 ;
1831
1832 for ( i = 0 ; i < 12 ; i ++ )
1833 {
1834 faceflag[ i ] = 1 ;
1835 int hasFeature = 1 ;
1836
1837 for ( j = 0 ; j < 4 ; j ++ )
1838 {
1839 nx = ox + sheetNeighbor[i][j][0] ;
1840 ny = oy + sheetNeighbor[i][j][1] ;
1841 nz = oz + sheetNeighbor[i][j][2] ;
1842
1843 if ( getDataAt( nx, ny, nz ) == 0 || nvol->getDataAt( nx, ny, nz ) == 1 )
1844 {
1845 hasFeature = 0 ;
1846 }
1847 if ( getDataAt( nx, ny, nz ) < 0 )
1848 {
1849 faceflag[ i ] = 0 ;
1850 break ;
1851 }
1852 }
1853 if ( faceflag[ i ] == 1 && hasFeature )
1854 {
1855 hasFeatureFace ++ ;
1856 // return 0 ;
1857 }
1858
1859 if ( faceflag[ i ] )
1860 {
1861 int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
1862 edge[ e0 ] ++ ;
1863 edge[ e1 ] ++ ;
1864 tot ++ ;
1865 }
1866 }
1867
1868 if ( tot == 0 || hasFeatureFace == 0 )
1869 {
1870 return 0 ;
1871 }
1872
1873 // Removing 1s
1874 int numones = 0 ;
1875 for ( i = 0 ; i < 6 ; i ++ )
1876 {
1877 if ( edge[ i ] == 1 )
1878 {
1879 numones ++ ;
1880 }
1881 }
1882 while( numones > 0 )
1883 {
1884 int e = 0;
1885 for ( i = 0 ; i < 6 ; i ++ )
1886 {
1887 if ( edge[ i ] == 1 )
1888 {
1889 e = i ;
1890 break ;
1891 }
1892 }
1893 /*
1894 if ( edge[ e ] != 1 )
1895 {
1896 printf("Wrong Again!********\n") ;
1897 }
1898 */
1899
1900 int f, e2 ;
1901 for ( j = 0 ; j < 4 ; j ++ )
1902 {
1903 f = edgeFaces[ e ][ j ] ;
1904 if ( faceflag[ f ] )
1905 {
1906 break ;
1907 }
1908 }
1909
1910 /*
1911 if ( faceflag[ f ] == 0 )
1912 {
1913 printf("Wrong!********\n") ;
1914 }
1915 */
1916
1917 if ( faceEdges[ f ][ 0 ] == e )
1918 {
1919 e2 = faceEdges[ f ][ 1 ] ;
1920 }
1921 else
1922 {
1923 e2 = faceEdges[ f ][ 0 ] ;
1924 }
1925
1926
1927 edge[ e ] -- ;
1928 numones -- ;
1929 edge[ e2 ] -- ;
1930 faceflag[ f ] = 0 ;
1931 tot -- ;
1932
1933 if ( edge[ e2 ] == 1 )
1934 {
1935 numones ++ ;
1936 }
1937 else if ( edge[ e2 ] == 0 )
1938 {
1939 numones -- ;
1940 }
1941 }
1942
1943 if ( tot > 0 )
1944 {
1945 return 0 ;
1946 }
1947
1948 return 1 ;
1949 }
1950
1951 //int Volume::getNumFaces( int ox, int oy, int oz )
1952 //{
1953 //int i, j ;
1954 //int nx, ny, nz ;
1955
1956 //int faces = 12 ;
1957 //for ( i = 0 ; i < 12 ; i ++ )
1958 //{
1959 //for ( j = 0 ; j < 4 ; j ++ )
1960 //{
1961 //nx = ox + sheetNeighbor[i][j][0] ;
1962 //ny = oy + sheetNeighbor[i][j][1] ;
1963 //nz = oz + sheetNeighbor[i][j][2] ;
1964
1965 //if ( getDataAt( nx, ny, nz ) < 0 )
1966 //{
1967 //faces -- ;
1968 //break ;
1969 //}
1970 //}
1971 //}
1972
1973 //return faces ;
1974 //}
1975
1976 //int Volume::getNumCells( int ox, int oy, int oz )
1977 //{
1978 //int i, j, k ;
1979
1980 //int cells = 0 ;
1981
1982 //for ( i = -1 ; i < 1 ; i ++ )
1983 //for ( j = -1 ; j < 1 ; j ++ )
1984 //for ( k = -1 ; k < 1 ; k ++ )
1985 //{
1986 //if ( hasCell( ox + i, oy + j, oz + k ) )
1987 //{
1988 //cells ++ ;
1989 //}
1990 //}
1991
1993 //return cells ;
1994 //}
1995
1996
1997 //int Volume::getNumIsolatedEdges( int ox, int oy, int oz )
1998 //{
1999 //int i, j, k ;
2000 //int nx, ny, nz ;
2001
2002 //double vox[3][3][3] ;
2003 //for ( i = -1 ; i < 2 ; i ++ )
2004 //for ( j = -1 ; j < 2 ; j ++ )
2005 //for ( k = -1 ; k < 2 ; k ++ )
2006 //{
2007 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ;
2008 //}
2009
2010
2011 //int edge[6] = { 0,0,0,0,0,0 } ;
2012
2013 //for ( i = 0 ; i < 12 ; i ++ )
2014 //{
2015 //int flag = 1 ;
2016 //for ( j = 0 ; j < 4 ; j ++ )
2017 //{
2018 //nx = 1 + sheetNeighbor[i][j][0] ;
2019 //ny = 1 + sheetNeighbor[i][j][1] ;
2020 //nz = 1 + sheetNeighbor[i][j][2] ;
2021
2022 //if ( vox[nx][ny][nz] < 0 )
2023 //{
2024 //flag = 0 ;
2025 //break ;
2026 //}
2027 //}
2028
2029 //if ( flag )
2030 //{
2031 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ;
2032 //edge[ e0 ] ++ ;
2033 //edge[ e1 ] ++ ;
2034 //}
2035 //}
2036
2037 //int edges = 0 ;
2038 //for ( i = 0 ; i < 6 ; i ++ )
2039 //{
2040 //if ( edge[i] )
2041 //{
2042 //continue ;
2043 //}
2044
2045 //nx = 1 + neighbor6[i][0] ;
2046 //ny = 1 + neighbor6[i][1] ;
2047 //nz = 1 + neighbor6[i][2] ;
2048
2049 //if ( vox[nx][ny][nz] >= 0 )
2050 //{
2051 //edges ++ ;
2052 //}
2053
2054 //}
2055
2056 //return edges ;
2057 //}
2058
2059 //int Volume::getNumIsolatedFaces( int ox, int oy, int oz )
2060 //{
2061 //int i, j, k ;
2062 //int nx, ny, nz ;
2063
2064 //int faces = 0 ;
2065 //int cellflag[ 8 ] ;
2066
2067 //int ct = 0 ;
2068 //for ( i = -1 ; i < 1 ; i ++ )
2069 //for ( j = -1 ; j < 1 ; j ++ )
2070 //for ( k = -1 ; k < 1 ; k ++ )
2071 //{
2072 //if ( hasCell( ox + i, oy + j, oz + k ) )
2073 //{
2074 //cellflag[ ct ] = 1 ;
2075 //}
2076 //else
2077 //{
2078 //cellflag[ ct ] = 0 ;
2079 //}
2080 //ct ++ ;
2081 //}
2082
2083 //for ( i = 0 ; i < 12 ; i ++ )
2084 //{
2085 //int hasFace = 1 ;
2086 //for ( j = 0 ; j < 4 ; j ++ )
2087 //{
2088 //nx = ox + sheetNeighbor[i][j][0] ;
2089 //ny = oy + sheetNeighbor[i][j][1] ;
2090 //nz = oz + sheetNeighbor[i][j][2] ;
2091
2092 //if ( getDataAt( nx, ny, nz ) < 0 )
2093 //{
2094 //hasFace = 0 ;
2095 //break ;
2096 //}
2097 //}
2098
2099 //if ( hasFace )
2100 //{
2101 //if (cellflag[ faceCells[i][0] ] == 0 && cellflag[ faceCells[i][1] ] == 0 )
2102 //{
2103 //faces ++ ;
2104 //}
2105 //}
2106 //}
2107
2109 //return faces ;
2110 //}
2111
2112 //int Volume::isFeatureFace2( int ox, int oy, int oz )
2113 //{
2114 //int i ;
2115 //int nx, ny, nz ;
2116
2117 //for ( i = 0 ; i < 6 ; i ++ )
2118 //{
2119 //nx = ox + neighbor6[i][0] ;
2120 //ny = oy + neighbor6[i][1] ;
2121 //nz = oz + neighbor6[i][2] ;
2122
2123 //double val = getDataAt( nx, ny, nz ) ;
2124
2125 //if ( val == 0 )
2126 //{
2127 //return 0 ;
2128 //}
2129
2130 //}
2131
2132 //return 1 ;
2133 //}
2134
2135 int Volume::isFeatureFace( int ox, int oy, int oz )
2136 {
2137 // return 1 ;
2138
2139 int i, j ;
2140 int nx, ny, nz ;
2141
2142 int faces = 12 ;
2143 for ( i = 0 ; i < 12 ; i ++ )
2144 {
2145 int ct = 0 ;
2146 for ( j = 0 ; j < 4 ; j ++ )
2147 {
2148 nx = ox + sheetNeighbor[i][j][0] ;
2149 ny = oy + sheetNeighbor[i][j][1] ;
2150 nz = oz + sheetNeighbor[i][j][2] ;
2151
2152 if ( getDataAt( nx, ny, nz ) < 0 )
2153 {
2154 ct = -1 ;
2155 break ;
2156 }
2157 else if ( getNumNeighbor6( nx, ny, nz ) == 6 )
2158 {
2159 ct = -1 ;
2160 break ;
2161 }
2162 // else if ( getDataAt( nx, ny, nz ) == 0 )
2163 // {
2164 // ct ++ ;
2165 // }
2166
2167
2168 }
2169 if ( ct == -1 || ct >= 1 )
2170 {
2171 faces -- ;
2172 }
2173 }
2174
2175 if ( faces > 0 )
2176 {
2177 return 1 ;
2178 }
2179 return 0 ;
2180
2181 }
2182
2183 //int Volume::hasFeatureFace( int ox, int oy, int oz )
2184 //{
2185 //int i, j, k ;
2186 //int nx, ny, nz ;
2187
2188 //int faceflag ;
2189 //int cellflag[ 8 ] ;
2190
2192 //int ct = 0 ;
2193 //for ( i = -1 ; i < 1 ; i ++ )
2194 //for ( j = -1 ; j < 1 ; j ++ )
2195 //for ( k = -1 ; k < 1 ; k ++ )
2196 //{
2197 //if ( hasCell( ox + i, oy + j, oz + k ) )
2198 //{
2199 //cellflag[ ct ] = 1 ;
2200 //}
2201 //else
2202 //{
2203 //cellflag[ ct ] = 0 ;
2204 //}
2205 //ct ++ ;
2206 //}
2207
2209 //for ( i = 0 ; i < 12 ; i ++ )
2210 //{
2211 //faceflag = 1 ;
2212 //for ( j = 0 ; j < 4 ; j ++ )
2213 //{
2214 //nx = ox + sheetNeighbor[i][j][0] ;
2215 //ny = oy + sheetNeighbor[i][j][1] ;
2216 //nz = oz + sheetNeighbor[i][j][2] ;
2217
2218 //if ( getDataAt( nx, ny, nz ) < 0 )
2219 //{
2220 //faceflag = 0 ;
2221 //break ;
2222 //}
2223 //}
2224
2225 //if ( faceflag )
2226 //{
2227 //if ( cellflag[ faceCells[i][0] ] == 0 && cellflag[ faceCells[i][1] ] == 0 )
2228 //{
2229 //return 1 ;
2230 //}
2231 //}
2232 //}
2233
2234 //return 0 ;
2235
2236 //}
2237
2238 int Volume::isSheetEnd( int ox, int oy, int oz )
2239 {
2240 return ( ( hasCompleteSheet(ox,oy,oz) == 0 ) && isFeatureFace(ox,oy,oz) ) ;
2241
2243 //
2244 //int i, j, k ;
2245 //int nx, ny, nz ;
2246
2247 //double vox[3][3][3] ;
2248 //for ( i = -1 ; i < 2 ; i ++ )
2249 // for ( j = -1 ; j < 2 ; j ++ )
2250 // for ( k = -1 ; k < 2 ; k ++ )
2251 // {
2252 // vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ;
2253 // }
2254
2255 //int edge[6] = { 4,4,4,4,4,4 } ;
2256 //int edge2[6] = { 4,4,4,4,4,4 } ;
2257
2258 //for ( i = 0 ; i < 12 ; i ++ )
2259 //{
2260 // for ( j = 0 ; j < 4 ; j ++ )
2261 // {
2262 // nx = 1 + sheetNeighbor[i][j][0] ;
2263 // ny = 1 + sheetNeighbor[i][j][1] ;
2264 // nz = 1 + sheetNeighbor[i][j][2] ;
2265
2266 // if ( vox[nx][ny][nz] < 0 )
2267 // {
2268 // edge[ faceEdges[ i ][ 0 ] ] -- ;
2269 // edge[ faceEdges[ i ][ 1 ] ] -- ;
2270 // break ;
2271 // }
2272 // }
2273
2274 // for ( j = 0 ; j < 4 ; j ++ )
2275 // {
2276 // nx = 1 + sheetNeighbor[i][j][0] ;
2277 // ny = 1 + sheetNeighbor[i][j][1] ;
2278 // nz = 1 + sheetNeighbor[i][j][2] ;
2279
2280 // if ( vox[nx][ny][nz] <= 0 )
2281 // {
2282 // edge2[ faceEdges[ i ][ 0 ] ] -- ;
2283 // edge2[ faceEdges[ i ][ 1 ] ] -- ;
2284 // break ;
2285 // }
2286 // }
2287 //}
2288
2289 //
2291 //for ( i = 0 ; i < 6 ; i ++ )
2292 //{
2293 // nx = 1 + neighbor6[i][0] ;
2294 // ny = 1 + neighbor6[i][1] ;
2295 // nz = 1 + neighbor6[i][2] ;
2296 // if ( edge[i] == 0 && vox[nx][ny][nz] >= 0 )
2297 // {
2298 // return 0 ;
2299 // }
2300 //}
2301 //*/
2302 //
2303
2304
2305 //for ( i = 0 ; i < 6 ; i ++ )
2306 //{
2307 // if ( edge[ i ] == 1 && edge2[ i ] == 1 )
2308 // {
2309 // return 1 ;
2310 // }
2311 //}
2312
2313 //return 0 ;
2314 }
2315
2316 int Volume::isSimple( int ox, int oy, int oz )
2317 {
2318 /* Test if this is a simple voxel */
2319 // int flag = 0 ;
2320 double vox[3][3][3] ;
2321
2322 int i, j, k ;
2323 for ( i = -1 ; i < 2 ; i ++ )
2324 for ( j = -1 ; j < 2 ; j ++ )
2325 for ( k = -1 ; k < 2 ; k ++ )
2326 {
2327 double tval = getDataAt( ox + i, oy + j, oz + k ) ;
2328
2329 /*
2330 if ( tval == 0 || tval > (va )
2331 {
2332 flag = 1 ;
2333 }
2334 */
2335 /*
2336 if ( tval < 0 && tval == - curwid )
2337 {
2338 printf("Here %d", (int)-tval) ;
2339 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
2340 }
2341 else
2342 */
2343 {
2344 vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
2345 }
2346 }
2347
2348 /* Debugging
2349 printf("{") ;
2350 for ( i = 0 ; i < 3 ; i ++ )
2351 {
2352 if ( i ) printf(",") ;
2353 printf("{") ;
2354 for ( j = 0 ; j < 3 ; j ++ )
2355 {
2356 if ( j ) printf(",") ;
2357 printf("{") ;
2358 for ( k = 0 ; k < 3 ; k ++ )
2359 {
2360 if ( k ) printf(",") ;
2361 printf("%d", (vox[i][j][k] >=0 ? 1: 0));
2362 }
2363 printf("}") ;
2364 }
2365 printf("}") ;
2366 }
2367 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
2368 */
2369
2370 if ( countInt( vox ) == 1 && countExt( vox ) == 1 )
2371 {
2372 return 1 ;
2373 }
2374 else
2375 {
2376 return 0 ;
2377 }
2378 }
2379
2380
2381 int Volume::isPiercable( int ox, int oy, int oz )
2382 {
2383 /* Test if this is a simple voxel */
2384 // int flag = 0 ;
2385 double vox[3][3][3] ;
2386
2387 int i, j, k ;
2388 for ( i = -1 ; i < 2 ; i ++ )
2389 for ( j = -1 ; j < 2 ; j ++ )
2390 for ( k = -1 ; k < 2 ; k ++ )
2391 {
2392 double tval = getDataAt( ox + i, oy + j, oz + k ) ;
2393
2394 /*
2395 if ( tval == 0 || tval > (va )
2396 {
2397 flag = 1 ;
2398 }
2399 */
2400 /*
2401 if ( tval < 0 && tval == - curwid )
2402 {
2403 printf("Here %d", (int)-tval) ;
2404 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ;
2405 }
2406 else
2407 */
2408 {
2409 vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ;
2410 }
2411 }
2412
2413 /* Debugging
2414 printf("{") ;
2415 for ( i = 0 ; i < 3 ; i ++ )
2416 {
2417 if ( i ) printf(",") ;
2418 printf("{") ;
2419 for ( j = 0 ; j < 3 ; j ++ )
2420 {
2421 if ( j ) printf(",") ;
2422 printf("{") ;
2423 for ( k = 0 ; k < 3 ; k ++ )
2424 {
2425 if ( k ) printf(",") ;
2426 printf("%d", (vox[i][j][k] >=0 ? 1: 0));
2427 }
2428 printf("}") ;
2429 }
2430 printf("}") ;
2431 }
2432 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ;
2433 */
2434
2435 if ( countInt( vox ) == 1 && countExt( vox ) != 1 )
2436 {
2437 return 1 ;
2438 }
2439 else
2440 {
2441 return 0 ;
2442 }
2443 }
2444
2445
2446 //int Volume::isSimple2( int v[3][3][3] )
2447 //{
2449 //double vox[3][3][3] ;
2450
2451 //int i, j, k ;
2452 //for ( i = 0 ; i < 3 ; i ++ )
2453 //for ( j = 0 ; j < 3 ; j ++ )
2454 //for ( k = 0 ; k < 3 ; k ++ )
2455 //{
2456 //if ( v[i][j][k] == 0 )
2457 //{
2458 //vox[ i ][ j ][ k ] = 1 ;
2459 //}
2460 //else
2461 //{
2462 //vox[i][j][k] = -1 ;
2463 //}
2464 //}
2465 //if ( countInt( vox ) == 1 && countExt( vox ) == 1 )
2466 //{
2467 //return 1 ;
2468 //}
2469 //else
2470 //{
2471 //printf("Int: %d Ext: %d\n",countInt( vox ),countExt( vox ) );
2472 //return 0 ;
2473 //}
2474 //}
2475
2476 //int Volume::getNumPotComplex3( int ox, int oy, int oz )
2477 //{
2479
2480
2481 //int i, j, k ;
2482 //if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz ) )
2483 //{
2484 //return 60 ;
2485 //}
2486 //double val = getDataAt( ox, oy, oz ) ;
2487
2488 //int nx, ny, nz ;
2489
2490
2491 //int numSimple = 0 ;
2492 //for ( i = -1 ; i < 2 ; i ++ )
2493 //for ( j = -1 ; j < 2 ; j ++ )
2494 //for ( k = -1 ; k < 2 ; k ++ )
2495 //{
2496 //nx = ox + i ;
2497 //ny = oy + j ;
2498 //nz = oz + k ;
2499 //if ( getDataAt( nx, ny, nz ) > 0 )
2500 //{
2501 //if ( isSimple( nx, ny, nz ) && ! isSheetEnd( nx, ny, nz ) )
2502 //{
2503 //numSimple ++ ;
2504 //}
2505 //}
2506 //}
2507 //numSimple -- ;
2508
2509 //setDataAt( ox, oy, oz, -val ) ;
2510
2511 //int numPotSimple = 0 ;
2512 //for ( i = -1 ; i < 2 ; i ++ )
2513 //for ( j = -1 ; j < 2 ; j ++ )
2514 //for ( k = -1 ; k < 2 ; k ++ )
2515 //{
2516 //nx = ox + i ;
2517 //ny = oy + j ;
2518 //nz = oz + k ;
2519 //if ( getDataAt( nx, ny, nz ) > 0 )
2520 //{
2521 //if ( isSimple( nx, ny, nz ) && ! isSheetEnd( nx, ny, nz ) )
2522 //{
2523 //numPotSimple ++ ;
2524 //}
2525 //}
2526 //}
2527
2528
2529 //setDataAt( ox, oy, oz, val ) ;
2530
2531
2532
2533
2534 //return 30 - ( numPotSimple - numSimple ) ;
2535 //}
2536
2537 //int Volume::getNumPotComplex4( int ox, int oy, int oz )
2538 //{
2539
2540 //int cells = getNumCells(ox,oy,oz) ;
2542 //int faces = getNumFaces(ox,oy,oz) ;
2544
2545 //return ( cells * 6 - 2 * faces ) ; // + 2 * ( faces * 4 - 4 * iedges ) ;
2546 //}
2547
2548 int Volume::getNumPotComplex( int ox, int oy, int oz )
2549 {
2550 //return 0 ;
2551
2552 int i;
2553 double val = getDataAt( ox, oy, oz ) ;
2554 if ( val <= 0 )
2555 {
2556 // return 70 ;
2557 }
2558
2559 // return getNumNeighbor6( ox, oy, oz ) ;
2560
2561 // int v = ((getNumNeighbor6( ox, oy, oz ) & 255) << 24) ;
2562 //int v = 0 ;
2563
2564 int rvalue = 0, nx, ny, nz ;
2565 setDataAt( ox, oy, oz, -val ) ;
2566
2567 /*
2568 for ( i = -1 ; i < 2 ; i ++ )
2569 for ( j = -1 ; j < 2 ; j ++ )
2570 for ( k = -1 ; k < 2 ; k ++ )
2571 {
2572 nx = ox + i ;
2573 ny = oy + j ;
2574 nz = oz + k ;
2575 if ( getDataAt( nx, ny, nz ) == val )
2576 {
2577 if ( isSheetEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
2578 {
2579 rvalue ++ ;
2580 }
2581 }
2582 }
2583 */
2584
2585 for ( i = 0 ; i < 6 ; i ++ )
2586 {
2587 nx = ox + neighbor6[i][0] ;
2588 ny = oy + neighbor6[i][1] ;
2589 nz = oz + neighbor6[i][2] ;
2590
2591 if ( getDataAt(nx,ny,nz) >= 0 )
2592 {
2593 int num = getNumNeighbor6( nx, ny, nz ) ;
2594 if ( num > rvalue )
2595 {
2596 rvalue = num ;
2597 }
2598 }
2599 }
2600
2601
2602 setDataAt( ox, oy, oz, val ) ;
2603
2604 return rvalue + getNumNeighbor6( ox, oy, oz ) * 10 ;
2605 /*
2606 int v = (((rvalue + getNumNeighbor6( ox, oy, oz ) * 10) & 255) << 24) ;
2607 v |= ( ( ox & 255 ) << 16 ) ;
2608 v |= ( ( oy & 255 ) << 8 ) ;
2609 v |= ( ( oz & 255 ) ) ;
2610 return v ;
2611 */
2612
2613 }
2614
2615 int Volume::getNumPotComplex2( int ox, int oy, int oz )
2616 {
2617 return getNumPotComplex( ox, oy, oz ) ;
2618
2619 //int i, j, k ;
2620 //double val = getDataAt( ox, oy, oz ) ;
2621 //if ( val <= 0 )
2622 //{
2623 // return 0 ;
2624 //}
2625
2626 //int rvalue = 0, nx, ny, nz ;
2627 //setDataAt( ox, oy, oz, -val ) ;
2628
2629 //for ( i = -1 ; i < 2 ; i ++ )
2630 // for ( j = -1 ; j < 2 ; j ++ )
2631 // for ( k = -1 ; k < 2 ; k ++ )
2632 // {
2633 // nx = ox + i ;
2634 // ny = oy + j ;
2635 // nz = oz + k ;
2636 // if ( getDataAt( nx, ny, nz ) == val )
2637 // {
2638 // if ( isHelixEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) )
2639 // {
2640 // rvalue ++ ;
2641 // }
2642 // }
2643 // }
2644
2645 //setDataAt( ox, oy, oz, val ) ;
2646
2647 //return rvalue + getNumNeighbor6( ox, oy, oz ) * 30 ;
2648 }
2649
2650 //int Volume::getNumNeighbor( int ox, int oy, int oz )
2651 //{
2652 //int rvalue = 0 ;
2653 //double val = getDataAt( ox, oy, oz ) ;
2654 //for ( int i = 0 ; i < 6 ; i ++ )
2655 //{
2656 //int nx = ox + neighbor6[i][0] ;
2657 //int ny = oy + neighbor6[i][1] ;
2658 //int nz = oz + neighbor6[i][2] ;
2659
2660 //if ( getDataAt( nx, ny, nz ) == val )
2661 //{
2662 //rvalue ++ ;
2663 //}
2664 //}
2666 //for ( int i = -1 ; i < 2 ; i ++ )
2667 //for ( int j = -1 ; j < 2 ; j ++ )
2668 //for ( int k = -1 ; k < 2 ; k ++ )
2669 //{
2670 //int nx = ox + i ;
2671 //int ny = oy + j ;
2672 //int nz = oz + k ;
2673
2674 //if ( getDataAt( nx, ny, nz ) == val )
2675 //{
2676 //rvalue ++ ;
2677 //}
2678 //}
2679 //*/
2680 //return rvalue ;
2681 //}
2682
2683
2684 //void Volume::setScoreNeighbor( GridQueue* queue )
2685 //{
2687 //gridQueueEle* ele = queue->getHead() ;
2688 //while ( ele != NULL )
2689 //{
2690 //ele->score = getNumNeighbor( ele->x, ele->y, ele->z ) ;
2691 //ele = ele->next ;
2692 //}
2693
2694 //queue->sort( queue->getNumElements() ) ;
2695 //}
2696
2697
2698 int Volume::components6( int vox[3][3][3] )
2699 {
2700 // Stupid flood fill
2701 int tot = 0 ;
2702 int queue[27][3] ;
2703 int vis[3][3][3] ;
2704 int head = 0, tail = 1 ;
2705 int i, j, k ;
2706 for ( i = 0 ; i < 3 ; i ++ )
2707 for ( j = 0 ; j < 3 ; j ++ )
2708 for ( k = 0 ; k < 3 ; k ++ )
2709 {
2710 vis[i][j][k] = 0 ;
2711 if ( vox[i][j][k] )
2712 {
2713 if ( tot == 0 )
2714 {
2715 queue[0][0] = i ;
2716 queue[0][1] = j ;
2717 queue[0][2] = k ;
2718 vis[i][j][k] = 1 ;
2719 }
2720 tot ++ ;
2721 }
2722 }
2723 if ( tot == 0 )
2724 {
2725 return 0 ;
2726 }
2727 // printf("total: %d\n", tot) ;
2728
2729 int ct = 1 ;
2730 while ( head != tail )
2731 {
2732 int x = queue[head][0] ;
2733 int y = queue[head][1] ;
2734 int z = queue[head][2] ;
2735 head ++ ;
2736
2737 for ( i = 0 ; i < 6 ; i ++ )
2738 {
2739 int nx = x + neighbor6[i][0] ;
2740 int ny = y + neighbor6[i][1] ;
2741 int nz = z + neighbor6[i][2] ;
2742 if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
2743 {
2744 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
2745 {
2746 queue[tail][0] = nx ;
2747 queue[tail][1] = ny ;
2748 queue[tail][2] = nz ;
2749 tail ++ ;
2750 vis[nx][ny][nz] = 1 ;
2751 ct ++ ;
2752 }
2753 }
2754 }
2755 }
2756
2757 if ( ct == tot )
2758 {
2759 return 1 ;
2760 }
2761 else
2762 {
2763 return 2 ;
2764 }
2765
2766 }
2767 int Volume::components26( int vox[3][3][3] )
2768 {
2769 // Stupid flood fill
2770 int tot = 0 ;
2771 int queue[27][3] ;
2772 int vis[3][3][3] ;
2773 int head = 0, tail = 1 ;
2774 int i, j, k ;
2775 for ( i = 0 ; i < 3 ; i ++ )
2776 for ( j = 0 ; j < 3 ; j ++ )
2777 for ( k = 0 ; k < 3 ; k ++ )
2778 {
2779 vis[i][j][k] = 0 ;
2780 if ( vox[i][j][k] )
2781 {
2782 if ( tot == 0 )
2783 {
2784 queue[0][0] = i ;
2785 queue[0][1] = j ;
2786 queue[0][2] = k ;
2787 vis[i][j][k] = 1 ;
2788 }
2789 tot ++ ;
2790 }
2791 }
2792 if ( tot == 0 )
2793 {
2794 return 0 ;
2795 }
2796
2797 int ct = 1 ;
2798 while ( head != tail )
2799 {
2800 int x = queue[head][0] ;
2801 int y = queue[head][1] ;
2802 int z = queue[head][2] ;
2803 head ++ ;
2804
2805 for ( i = -1 ; i < 2 ; i ++ )
2806 for ( j = -1 ; j < 2 ; j ++ )
2807 for ( k = -1 ; k < 2 ; k ++ )
2808 {
2809 int nx = x + i ;
2810 int ny = y + j ;
2811 int nz = z + k ;
2812 if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 )
2813 {
2814 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] )
2815 {
2816 queue[tail][0] = nx ;
2817 queue[tail][1] = ny ;
2818 queue[tail][2] = nz ;
2819 tail ++ ;
2820 vis[nx][ny][nz] = 1 ;
2821 ct ++ ;
2822 }
2823 }
2824 }
2825 }
2826
2827 if ( ct == tot )
2828 {
2829 return 1 ;
2830 }
2831 else
2832 {
2833 return 2 ;
2834 }
2835
2836 }
2837
2838 int Volume::countExt( double vox[3][3][3] )
2839 {
2840 int tvox[3][3][3] ;
2841
2842 for ( int i = 0 ; i < 3 ; i ++ )
2843 for ( int j = 0 ; j < 3 ; j ++ )
2844 for ( int k = 0 ; k < 3 ; k ++ )
2845 {
2846 if ( vox[i][j][k] < 0 )
2847 {
2848 tvox[i][j][k] = 1 ;
2849 }
2850 else
2851 {
2852 tvox[i][j][k] = 0 ;
2853 }
2854 }
2855
2856 return components26( tvox ) ;
2857 }
2858
2859 int Volume::countInt( double vox[3][3][3] )
2860 {
2861 int i, j, k ;
2862 int tvox[3][3][3] ;
2863
2864 for ( i = 0 ; i < 3 ; i ++ )
2865 for ( j = 0 ; j < 3 ; j ++ )
2866 for ( k = 0 ; k < 3 ; k ++ )
2867 {
2868 tvox[i][j][k] = 0 ;
2869 }
2870
2871 for ( i = 0 ; i < 6 ; i ++ )
2872 {
2873 int nx = 1 + neighbor6[i][0] ;
2874 int ny = 1 + neighbor6[i][1] ;
2875 int nz = 1 + neighbor6[i][2] ;
2876 if ( vox[ nx ][ ny ][ nz ] >= 0 )
2877 {
2878 tvox[ nx ][ ny ][ nz ] = 1 ;
2879 for ( j = 0 ; j < 4 ; j ++ )
2880 {
2881 int nnx = neighbor64[i][j][0] + nx ;
2882 int nny = neighbor64[i][j][1] + ny ;
2883 int nnz = neighbor64[i][j][2] + nz ;
2884 if ( vox[ nnx ][ nny ][ nnz ] >= 0 )
2885 {
2886 tvox[ nnx ][ nny ][ nnz ] = 1 ;
2887 }
2888 }
2889 }
2890 }
2891
2892 return components6( tvox ) ;
2893 }
2894
2895 int Volume::countIntEuler( int ox, int oy, int oz )
2896 {
2897 int nv = 0 , ne = 0, nc = 0 ;
2898
2899 int i, j, k ;
2900 int tvox[3][3][3] ;
2901 double vox[3][3][3] ;
2902
2903 for ( i = 0 ; i < 3 ; i ++ )
2904 for ( j = 0 ; j < 3 ; j ++ )
2905 for ( k = 0 ; k < 3 ; k ++ )
2906 {
2907 vox[i][j][k] = getDataAt( ox - 1 + i, oy - 1 + j, oz - 1 + k ) ;
2908 tvox[i][j][k] = 0 ;
2909 }
2910
2911 for ( i = 0 ; i < 6 ; i ++ )
2912 {
2913 int nx = 1 + neighbor6[i][0] ;
2914 int ny = 1 + neighbor6[i][1] ;
2915 int nz = 1 + neighbor6[i][2] ;
2916 if ( vox[nx][ny][nz] >= 0 )
2917 {
2918 tvox[ nx ][ ny ][ nz ] = 1 ;
2919
2920 nv ++ ;
2921
2922 for ( j = 0 ; j < 4 ; j ++ )
2923 {
2924 int nnx = neighbor64[i][j][0] + nx ;
2925 int nny = neighbor64[i][j][1] + ny ;
2926 int nnz = neighbor64[i][j][2] + nz ;
2927 if ( vox[nnx][nny][nnz] >= 0 )
2928 {
2929 if ( tvox[nnx][nny][nnz] == 0 )
2930 {
2931 tvox[nnx][nny][nnz] = 1 ;
2932 nv ++ ;
2933 }
2934
2935 ne ++ ;
2936 }
2937 }
2938 }
2939 }
2940
2941 nc = components6( tvox ) ;
2942
2943 return ( nc - ( nv - ne ) ) ;
2944 }
2945
2946 //void Volume::erodeNoTopo( float thr, int wid )
2947 //{
2949 //int i, j, k ;
2951 //#ifdef VERBOSE
2952 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ;
2953 //#endif
2954 //threshold( thr, -MAX_ERODE, 0 ) ;
2955
2957 //#ifdef VERBOSE
2958 //printf("Initializing queue...\n") ;
2959 //#endif
2960 //GridQueue* queue = new GridQueue( ) ;
2961
2962 //for ( i = 0 ; i < getSizeX() ; i ++ )
2963 //for ( j = 0 ; j < getSizeY() ; j ++ )
2964 //for ( k = 0 ; k < getSizeZ() ; k ++ )
2965 //{
2966 //if ( getDataAt( i, j, k ) >= 0 )
2967 //{
2968 //for ( int m = 0 ; m < 6 ; m ++ )
2969 //{
2970 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
2971 //{
2972 //setDataAt( i, j, k, 1 ) ;
2973 //queue->pushQueue( i, j, k ) ;
2974 //break ;
2975 //}
2976 //}
2977 //}
2978 //}
2979 //#ifdef VERBOSE
2980 //printf("Total %d nodes\n", queue->getNumElements() ) ;
2981
2983 //printf("Start erosion to %d...\n", wid) ;
2984 //#endif
2985 //double val = 0;
2986 //int ox, oy, oz ;
2987 //int curwid = 0 ;
2988 //int total = 0, ct = 0 ;
2989 //while ( 1 )
2990 //{
2991 //if ( ct == total )
2992 //{
2993 //#ifdef VERBOSE
2994 //printf("Layer %d has %d nodes.\n", (int) curwid, total ) ;
2995 //#endif
2996 //curwid ++ ;
2997 //ct = 0 ;
2998 //total = queue->getNumElements() ;
2999 //if ( total == 0 )
3000 //{
3001 //break ;
3002 //}
3003 //}
3004
3005 //queue->popQueue(ox, oy, oz) ;
3006 //val = getDataAt( ox, oy, oz ) ;
3007 //if ( val > wid )
3008 //{
3009 //break ;
3010 //}
3011 //ct ++ ;
3012
3013 //setDataAt( ox, oy, oz, -val ) ;
3014
3015
3016 //for ( int m = 0 ; m < 6 ; m ++ )
3017 //{
3018 //int nx = ox + neighbor6[m][0] ;
3019 //int ny = oy + neighbor6[m][1] ;
3020 //int nz = oz + neighbor6[m][2] ;
3021 //if ( getDataAt( nx, ny, nz ) == 0 )
3022 //{
3023 //setDataAt( nx, ny, nz, val + 1 ) ;
3024 //queue->pushQueue( nx, ny, nz ) ;
3025 //}
3026 //}
3027 //}
3028
3030 //#ifdef VERBOSE
3031 //printf("Thresholding the volume to 0/1...\n") ;
3032 //#endif
3033 //threshold( 0, 0, 1 ) ;
3034
3035 //}
3036
3037 //void Volume::erodeTopo( float thr, int wid )
3038 //{
3040 //int i, j, k ;
3042 //#ifdef VERBOSE
3043 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ;
3044 //#endif
3045 //threshold( thr, -MAX_ERODE, 0 ) ;
3046
3048 //#ifdef VERBOSE
3049 //printf("Initializing queue...\n") ;
3050 //#endif
3051 //GridQueue* queue = new GridQueue( ) ;
3052
3053 //for ( i = 0 ; i < getSizeX() ; i ++ )
3054 //for ( j = 0 ; j < getSizeY() ; j ++ )
3055 //for ( k = 0 ; k < getSizeZ() ; k ++ )
3056 //{
3057 //if ( getDataAt( i, j, k ) >= 0 )
3058 //{
3059 //for ( int m = 0 ; m < 6 ; m ++ )
3060 //{
3061 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
3062 //{
3063 //setDataAt( i, j, k, 1 ) ;
3064 //queue->pushQueue( i, j, k ) ;
3065 //break ;
3066 //}
3067 //}
3068 //}
3069 //}
3070 //#ifdef VERBOSE
3071 //printf("Total %d nodes\n", queue->getNumElements() ) ;
3072
3074 //printf("Start erosion to %d...\n", wid) ;
3075 //#endif
3076
3077 //double val = 0;
3078 //int ox, oy, oz ;
3079 //int curwid = 0 ;
3080 //int total = 0, ct = 0 ;
3081 //while ( 1 )
3082 //{
3083 //if ( ct == total )
3084 //{
3085 //#ifdef VERBOSE
3086 //printf("Layer %d has %d nodes.\n", (int) curwid, total ) ;
3087 //#endif
3088 //curwid ++ ;
3089 //ct = 0 ;
3090 //total = queue->getNumElements() ;
3091 //if ( total == 0 )
3092 //{
3093 //break ;
3094 //}
3095 //}
3096
3097 //queue->popQueue(ox, oy, oz) ;
3098 //val = getDataAt( ox, oy, oz ) ;
3099 //if ( val > wid )
3100 //{
3101 //break ;
3102 //}
3103 //ct ++ ;
3104
3105 //if ( isSimple( ox, oy, oz ) )
3106 //{
3108 //setDataAt( ox, oy, oz, -val ) ;
3109 //}
3110 //else
3111 //{
3113 //setDataAt( ox, oy, oz, val + 1 ) ;
3114 //queue->pushQueue( ox, oy, oz ) ;
3115 //}
3116
3117
3118 //for ( int m = 0 ; m < 6 ; m ++ )
3119 //{
3120 //int nx = ox + neighbor6[m][0] ;
3121 //int ny = oy + neighbor6[m][1] ;
3122 //int nz = oz + neighbor6[m][2] ;
3123 //if ( getDataAt( nx, ny, nz ) == 0 )
3124 //{
3125 //setDataAt( nx, ny, nz, val + 1 ) ;
3126 //queue->pushQueue( nx, ny, nz ) ;
3127 //}
3128 //}
3129 //}
3130
3132 //#ifdef VERBOSE
3133 //printf("Thresholding the volume to 0/1...\n") ;
3134 //#endif
3135 //threshold( 0, 0, 1 ) ;
3136
3137 //}
3138
3139 //void Volume::erode2( float thr, int wid )
3140 //{
3141 //int i, j, k ;
3143 //#ifdef VERBOSE
3144 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ;
3145 //#endif
3146 //threshold( thr, -MAX_ERODE, 0 ) ;
3147
3149 //#ifdef VERBOSE
3150 //printf("Initializing queue...\n") ;
3151 //#endif
3152 //GridQueue2* queue = new GridQueue2( ) ;
3153 //GridQueue2* queue2 = new GridQueue2( ) ;
3154
3155 //for ( i = 0 ; i < getSizeX() ; i ++ )
3156 //for ( j = 0 ; j < getSizeY() ; j ++ )
3157 //for ( k = 0 ; k < getSizeZ() ; k ++ )
3158 //{
3159 //if ( getDataAt( i, j, k ) >= 0 )
3160 //{
3161 //for ( int m = 0 ; m < 6 ; m ++ )
3162 //{
3163 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
3164 //{
3165 //setDataAt( i, j, k, 1 ) ;
3166 //queue->prepend( i, j, k ) ;
3167 //break ;
3168 //}
3169 //}
3170 //}
3171 //}
3172 //#ifdef VERBOSE
3173 //printf("Total %d nodes\n", queue->getNumElements() ) ;
3174
3177 //printf("Start erosion to %d...\n", wid) ;
3178 //#endif
3179 //gridQueueEle* ele ;
3180 //int ox, oy, oz ;
3181
3182 //for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
3183 //{
3184 //int numComplex = 0, numSimple = 0 ;
3185 //#ifdef VERBOSE
3186 //printf("Processing %d nodes in layer %d\n", queue->getNumElements(), curwid) ;
3187 //#endif
3188
3190 //while ( ( ele = queue->getNext() ) != NULL )
3191 //{
3192 //for ( int m = 0 ; m < 6 ; m ++ )
3193 //{
3194 //int nx = ele->x + neighbor6[m][0] ;
3195 //int ny = ele->y + neighbor6[m][1] ;
3196 //int nz = ele->z + neighbor6[m][2] ;
3197 //if ( getDataAt( nx, ny, nz ) == 0 )
3198 //{
3199 //setDataAt( nx, ny, nz, curwid + 1 ) ;
3200 //queue2->prepend( nx, ny, nz ) ;
3201 //}
3202 //}
3203
3204 //}
3205 //*/
3206
3208 //int seed[3] = {-1,-1,-1};
3209 //queue->reset() ;
3210 //while ( queue->getNumElements() > 0 )
3211 //{
3212 //if ( seed[0] < 0 ) printf("After initial scoring...\n");
3213 //queue->reset() ;
3214 //ele = queue->getNext() ;
3215
3217 //while ( ele != NULL )
3218 //{
3219 //ox = ele->x ;
3220 //oy = ele->y ;
3221 //oz = ele->z ;
3222
3224 //if ( seed[0] < 0 ||
3225 //( ox < seed[0] + 2 && ox > seed[0] - 2 &&
3226 //oy < seed[1] + 2 && oy > seed[1] - 2 &&
3227 //oz < seed[2] + 2 && oz > seed[2] - 2 ) )
3228 //{
3229 //if ( ! isSimple( ox, oy, oz ) )
3230 //{
3232 //setDataAt( ox, oy, oz, curwid + 1 ) ;
3233 //queue2->prepend( ox, oy, oz ) ;
3234 //ele = queue->remove() ;
3235
3236 //numComplex ++ ;
3237 //}
3238 //else
3239 //{
3240 //ele = queue->getNext() ;
3241 //}
3242 //}
3243 //else
3244 //{
3245 //ele = queue->getNext() ;
3246 //}
3247 //}
3248
3250 //queue->reset() ;
3251 //ele = queue->getNext() ;
3252 //int preScore = -1;
3253 //while ( ele != NULL )
3254 //{
3255 //ox = ele->x ;
3256 //oy = ele->y ;
3257 //oz = ele->z ;
3258
3260 //if ( seed[0] < 0 ||
3261 //( ox < seed[0] + 3 && ox > seed[0] - 3 &&
3262 //oy < seed[1] + 3 && oy > seed[1] - 3 &&
3263 //oz < seed[2] + 3 && oz > seed[2] - 3 ) )
3264 //{
3265 //ele->score = getNumPotComplex( ox, oy, oz ) ;
3266 //}
3267
3268
3269 //if ( ele->score < preScore )
3270 //{
3272 //ele = queue->swap() ;
3273 //}
3274 //else
3275 //{
3276 //preScore = ele->score ;
3277 //}
3278
3280 //if ( ele->next == NULL )
3281 //{
3282 //ox = ele->x ;
3283 //oy = ele->y ;
3284 //oz = ele->z ;
3285 //setDataAt( ox, oy, oz, -1 ) ;
3287 //seed[0] = ox ;
3288 //seed[1] = oy ;
3289 //seed[2] = oz ;
3290 //queue->remove() ;
3292
3293
3294 //for ( int m = 0 ; m < 6 ; m ++ )
3295 //{
3296 //int nx = ox + neighbor6[m][0] ;
3297 //int ny = oy + neighbor6[m][1] ;
3298 //int nz = oz + neighbor6[m][2] ;
3299 //if ( getDataAt( nx, ny, nz ) == 0 )
3300 //{
3301 //setDataAt( nx, ny, nz, curwid + 1 ) ;
3302 //queue2->prepend( nx, ny, nz ) ;
3303 //}
3304 //}
3305
3306
3307 //numSimple ++ ;
3308 //ele = NULL ;
3309 //}
3310 //else
3311 //{
3312 //ele = queue->getNext() ;
3313 //}
3314 //}
3315
3316 //}
3317
3318 //delete queue ;
3319 //queue = queue2 ;
3320 //queue2 = new GridQueue2() ;
3321 //#ifdef VERBOSE
3322 //printf("%d complex, %d simple\n", numComplex, numSimple) ;
3323 //#endif
3324
3325 //if ( numSimple == 0 )
3326 //{
3327 //break ;
3328 //}
3329 //}
3330
3332 //#ifdef VERBOSE
3333 //printf("Thresholding the volume to 0/1...\n") ;
3334 //#endif
3335 //threshold( 0, 0, 1 ) ;
3336
3337 //}
3338
3339 //void Volume::erodeShapeTopo( float thr, int wid )
3340 //{
3342 //int i, j, k ;
3344 //#ifdef VERBOSE
3345 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ;
3346 //#endif
3347 //threshold( thr, -MAX_ERODE, 0 ) ;
3348
3350 //#ifdef VERBOSE
3351 //printf("Initializing queue...\n") ;
3352 //#endif
3353 //GridQueue2* queue2 = new GridQueue2( ) ;
3354 //GridQueue2* queue3 = new GridQueue2( ) ;
3355 //PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
3356
3357 //for ( i = 0 ; i < getSizeX() ; i ++ )
3358 //for ( j = 0 ; j < getSizeY() ; j ++ )
3359 //for ( k = 0 ; k < getSizeZ() ; k ++ )
3360 //{
3361 //if ( getDataAt( i, j, k ) >= 0 )
3362 //{
3363 //for ( int m = 0 ; m < 6 ; m ++ )
3364 //{
3365 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
3366 //{
3367 //setDataAt( i, j, k, 1 ) ;
3368 //queue2->prepend( i, j, k ) ;
3369 //break ;
3370 //}
3371 //}
3372 //}
3373 //}
3374 //#ifdef VERBOSE
3375 //printf("Total %d nodes\n", queue2->getNumElements() ) ;
3376
3377
3380 //printf("Start erosion to %d...\n", wid) ;
3381 //#endif
3382 //gridQueueEle* ele ;
3383 //gridPoint* gp ;
3384 //int ox, oy, oz ;
3385 //int score ;
3386 //Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
3387 //for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
3388 //{
3389 //scrvol->setDataAt( i, -1 ) ;
3390 //}
3391
3392 //for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
3393 //{
3397
3398 //int numComplex = 0, numSimple = 0 ;
3399 //#ifdef VERBOSE
3400 //printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
3401 //#endif
3402
3406 //queue2->reset() ;
3407 //ele = queue2->getNext() ;
3408 //while ( ele != NULL )
3409 //{
3410 //ox = ele->x ;
3411 //oy = ele->y ;
3412 //oz = ele->z ;
3413
3415 //if ( ! isSimple( ox, oy, oz ) )
3416 //{
3418 //setDataAt( ox, oy, oz, curwid + 1 ) ;
3419 //queue3->prepend( ox, oy, oz ) ;
3420 //ele = queue2->remove() ;
3421
3422 //numComplex ++ ;
3423 //}
3424 //else
3425 //{
3426 //ele = queue2->getNext() ;
3427 //}
3428 //}
3429
3433 //queue2->reset() ;
3434 //ele = queue2->getNext() ;
3435 //while ( ele != NULL )
3436 //{
3437 //ox = ele->x ;
3438 //oy = ele->y ;
3439 //oz = ele->z ;
3440
3442 //score = getNumPotComplex( ox, oy, oz ) ;
3443 //scrvol->setDataAt( ox, oy, oz, score ) ;
3444
3446 //gp = new gridPoint ;
3447 //gp->x = ox ;
3448 //gp->y = oy ;
3449 //gp->z = oz ;
3450 //queue->add( gp, -score ) ;
3451
3452 //ele = queue2->remove() ;
3453 //}
3454
3457 //delete queue2 ;
3458 //queue2 = queue3 ;
3459 //queue3 = new GridQueue2( ) ;
3460
3462 //while ( ! queue->isEmpty() )
3463 //{
3465 //queue->remove( gp, score ) ;
3466 //ox = gp->x ;
3467 //oy = gp->y ;
3468 //oz = gp->z ;
3469 //delete gp ;
3470 //score = -score ;
3471
3475 //if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
3476 //{
3477 //continue ;
3478 //}
3479
3481 //setDataAt( ox, oy, oz, -1 ) ;
3482 //numSimple ++ ;
3484
3486 //for ( int m = 0 ; m < 6 ; m ++ )
3487 //{
3488 //int nx = ox + neighbor6[m][0] ;
3489 //int ny = oy + neighbor6[m][1] ;
3490 //int nz = oz + neighbor6[m][2] ;
3491 //if ( getDataAt( nx, ny, nz ) == 0 )
3492 //{
3493 //setDataAt( nx, ny, nz, curwid + 1 ) ;
3494 //queue2->prepend( nx, ny, nz ) ;
3495 //}
3496 //}
3497
3500 //for ( i = -1 ; i < 2 ; i ++ )
3501 //for ( j = -1 ; j < 2 ; j ++ )
3502 //for ( k = -1 ; k < 2 ; k ++ )
3503 //{
3504 //int nx = ox + i ;
3505 //int ny = oy + j ;
3506 //int nz = oz + k ;
3507
3509 //if ( getDataAt( nx, ny, nz ) == curwid && ! isSimple( nx, ny, nz ) )
3510 //{
3512 //setDataAt( nx, ny, nz, curwid + 1 ) ;
3513 //queue2->prepend( nx, ny, nz ) ;
3514 //numComplex ++ ;
3515 //}
3516 //}
3517
3521 //for ( i = -2 ; i < 3 ;i ++ )
3522 //for ( j = -2 ; j < 3 ; j ++ )
3523 //for ( k = -2 ; k < 3 ; k ++ )
3524 //{
3525 //int nx = ox + i ;
3526 //int ny = oy + j ;
3527 //int nz = oz + k ;
3528
3529 //if ( getDataAt( nx, ny, nz ) == curwid )
3530 //{
3532 //score = getNumPotComplex( nx, ny, nz ) ;
3533
3534 //if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
3535 //{
3537 //scrvol->setDataAt( nx, ny, nz, score ) ;
3539 //gp = new gridPoint ;
3540 //gp->x = nx ;
3541 //gp->y = ny ;
3542 //gp->z = nz ;
3543 //queue->add( gp, -score ) ;
3544 //}
3545 //}
3546 //}
3547 //*/
3548
3549 //}
3550
3551 //#ifdef VERBOSE
3552 //printf("%d complex, %d simple\n", numComplex, numSimple) ;
3553 //#endif
3554 //if ( numSimple == 0 )
3555 //{
3556 //break ;
3557 //}
3558 //}
3559
3561 //#ifdef VERBOSE
3562 //printf("Thresholding the volume to 0/1...\n") ;
3563 //#endif
3564 //threshold( 0, 0, 1 ) ;
3565 //delete queue;
3566
3567 //}
3568
3569
3570 //void Volume::erodeAtom( float thr, int wid, Volume* avol )
3571 //{
3573 //int i, j, k ;
3575 //#ifdef VERBOSE
3576 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ;
3577 //#endif
3578 //threshold( thr, -MAX_ERODE, 0 ) ;
3579
3581 //#ifdef VERBOSE
3582 //printf("Initializing queue...\n") ;
3583 //#endif
3584 //GridQueue2* queue2 = new GridQueue2( ) ;
3585 //GridQueue2* queue3 = new GridQueue2( ) ;
3586 //PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
3587
3588 //for ( i = 0 ; i < getSizeX() ; i ++ )
3589 //for ( j = 0 ; j < getSizeY() ; j ++ )
3590 //for ( k = 0 ; k < getSizeZ() ; k ++ )
3591 //{
3592 //if ( getDataAt( i, j, k ) >= 0 )
3593 //{
3594 //if ( avol->getDataAt(i,j,k) > 0 )
3595 //{
3596 //setDataAt( i, j, k, MAX_ERODE ) ;
3597 //}
3598 //else
3599 //{
3600 //for ( int m = 0 ; m < 6 ; m ++ )
3601 //{
3602 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
3603 //{
3604 //setDataAt( i, j, k, 1 ) ;
3605 //queue2->prepend( i, j, k ) ;
3606 //break ;
3607 //}
3608 //}
3609 //}
3610 //}
3611 //}
3612 //#ifdef VERBOSE
3613 //printf("Total %d nodes\n", queue2->getNumElements() ) ;
3614
3615
3618 //printf("Start erosion to %d...\n", wid) ;
3619 //#endif
3620
3621 //gridQueueEle* ele ;
3622 //gridPoint* gp ;
3623 //int ox, oy, oz ;
3624 //int score ;
3625 //Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
3626 //for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
3627 //{
3628 //scrvol->setDataAt( i, -1 ) ;
3629 //}
3630
3631 //for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
3632 //{
3636
3637 //int numComplex = 0, numSimple = 0 ;
3638 //#ifdef VERBOSE
3639 //printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
3640 //#endif
3641
3645 //queue2->reset() ;
3646 //ele = queue2->getNext() ;
3647 //while ( ele != NULL )
3648 //{
3649 //ox = ele->x ;
3650 //oy = ele->y ;
3651 //oz = ele->z ;
3652
3654 //if ( ! isSimple( ox, oy, oz ) )
3655 //{
3657 //setDataAt( ox, oy, oz, curwid + 1 ) ;
3658 //queue3->prepend( ox, oy, oz ) ;
3659 //ele = queue2->remove() ;
3660
3661 //numComplex ++ ;
3662 //}
3663 //else
3664 //{
3665 //ele = queue2->getNext() ;
3666 //}
3667 //}
3668
3672 //queue2->reset() ;
3673 //ele = queue2->getNext() ;
3674 //while ( ele != NULL )
3675 //{
3676 //ox = ele->x ;
3677 //oy = ele->y ;
3678 //oz = ele->z ;
3679
3681 //score = getNumPotComplex( ox, oy, oz ) ;
3682 //scrvol->setDataAt( ox, oy, oz, score ) ;
3683
3685 //gp = new gridPoint ;
3686 //gp->x = ox ;
3687 //gp->y = oy ;
3688 //gp->z = oz ;
3689 //queue->add( gp, -score ) ;
3690
3691 //ele = queue2->remove() ;
3692 //}
3693
3696 //delete queue2 ;
3697 //queue2 = queue3 ;
3698 //queue3 = new GridQueue2( ) ;
3699
3701 //while ( ! queue->isEmpty() )
3702 //{
3704 //queue->remove( gp, score ) ;
3705 //ox = gp->x ;
3706 //oy = gp->y ;
3707 //oz = gp->z ;
3708 //delete gp ;
3709 //score = -score ;
3710
3714 //if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
3715 //{
3716 //continue ;
3717 //}
3718
3720 //setDataAt( ox, oy, oz, -1 ) ;
3721 //numSimple ++ ;
3723
3725 //for ( int m = 0 ; m < 6 ; m ++ )
3726 //{
3727 //int nx = ox + neighbor6[m][0] ;
3728 //int ny = oy + neighbor6[m][1] ;
3729 //int nz = oz + neighbor6[m][2] ;
3730 //if ( getDataAt( nx, ny, nz ) == 0 )
3731 //{
3732 //setDataAt( nx, ny, nz, curwid + 1 ) ;
3733 //queue2->prepend( nx, ny, nz ) ;
3734 //}
3735 //}
3736
3739 //for ( i = -1 ; i < 2 ; i ++ )
3740 //for ( j = -1 ; j < 2 ; j ++ )
3741 //for ( k = -1 ; k < 2 ; k ++ )
3742 //{
3743 //int nx = ox + i ;
3744 //int ny = oy + j ;
3745 //int nz = oz + k ;
3746
3748 //if ( getDataAt( nx, ny, nz ) == curwid && ! isSimple( nx, ny, nz ) )
3749
3750 //{
3752 //setDataAt( nx, ny, nz, curwid + 1 ) ;
3753 //queue2->prepend( nx, ny, nz ) ;
3754 //numComplex ++ ;
3755 //}
3756 //}
3757
3761 //for ( i = -2 ; i < 3 ;i ++ )
3762 //for ( j = -2 ; j < 3 ; j ++ )
3763 //for ( k = -2 ; k < 3 ; k ++ )
3764 //{
3765 //int nx = ox + i ;
3766 //int ny = oy + j ;
3767 //int nz = oz + k ;
3768
3769 //if ( getDataAt( nx, ny, nz ) == curwid )
3770 //{
3772 //score = getNumPotComplex( nx, ny, nz ) ;
3773
3774 //if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
3775 //{
3777 //scrvol->setDataAt( nx, ny, nz, score ) ;
3779 //gp = new gridPoint ;
3780 //gp->x = nx ;
3781 //gp->y = ny ;
3782 //gp->z = nz ;
3783 //queue->add( gp, -score ) ;
3784 //}
3785 //}
3786 //}
3787 //*/
3788
3789 //}
3790 //#ifdef VERBOSE
3791 //printf("%d complex, %d simple\n", numComplex, numSimple) ;
3792 //#endif
3793
3794 //if ( numSimple == 0 )
3795 //{
3796 //break ;
3797 //}
3798 //}
3799
3801 //#ifdef VERBOSE
3802 //printf("Thresholding the volume to 0/1...\n") ;
3803 //#endif
3804 //threshold( 0, 0, 1 ) ;
3805 //delete queue;
3806
3807 //}
3808
3810 //* Assuming the current volume has already been thresholded to 0/1
3811 //*/
3812 void Volume::curveSkeleton( Volume* grayvol, float lowthr, float highthr, Volume* svol )
3813 {
3814 int i, j, k ;
3815 // First, threshold the volume
3816 #ifdef VERBOSE
3817 printf("Thresholding the volume to -1/0...\n") ;
3818 #endif
3819 threshold( 0.5f, -1, 0 ) ;
3820
3821 // Next, apply convergent erosion
3822 // by preserving: complex nodes, curve end-points, and sheet points
3823
3824 // Next, initialize the linked queue
3825 #ifdef VERBOSE
3826 printf("Initializing queue...\n") ;
3827 #endif
3828 GridQueue2* queue2 = new GridQueue2( ) ;
3829 GridQueue2* queue3 = new GridQueue2( ) ;
3830 GridQueue2* queue4 = new GridQueue2( ) ;
3831 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
3832
3833 for ( i = 0 ; i < getSizeX() ; i ++ )
3834 for ( j = 0 ; j < getSizeY() ; j ++ )
3835 for ( k = 0 ; k < getSizeZ() ; k ++ )
3836 {
3837 if ( getDataAt( i, j, k ) >= 0 )
3838 {
3839 float v = (float)grayvol->getDataAt(i,j,k) ;
3840 if ( v <= lowthr || v > highthr || svol->getDataAt(i,j,k) > 0 )
3841 {
3842 setDataAt( i, j, k, MAX_ERODE ) ;
3843 }
3844 else
3845 {
3846 for ( int m = 0 ; m < 6 ; m ++ )
3847 {
3848 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
3849 {
3850 // setDataAt( i, j, k, 1 ) ;
3851 queue2->prepend( i, j, k ) ;
3852 break ;
3853 }
3854 }
3855 }
3856 }
3857 }
3858 int wid = MAX_ERODE ;
3859 #ifdef VERBOSE
3860 printf("Total %d nodes\n", queue2->getNumElements() ) ;
3861 printf("Start erosion to %d...\n", wid) ;
3862 #endif
3863
3864
3865 // Perform erosion
3866 gridQueueEle* ele ;
3867 gridPoint* gp ;
3868 int ox, oy, oz ;
3869 int score ;
3870 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
3871 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
3872 {
3873 scrvol->setDataAt( i, -1 ) ;
3874 }
3875
3876 #ifdef NOISE_DIS_HELIX
3877 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
3878 #endif
3879
3880 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
3881 {
3882 // At the start of each iteration,
3883 // queue2 holds all the nodes for this layer
3884 // queue3 and queue are empty
3885
3886 int numComplex = 0, numSimple = 0 ;
3887 #ifdef VERBOSE
3888 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
3889 #endif
3890
3891 /*
3892 We first need to assign curwid + 1 to every node in this layer
3893 */
3894 queue2->reset() ;
3895 ele = queue2->getNext() ;
3896 while ( ele != NULL )
3897 {
3898 ox = ele->x ;
3899 oy = ele->y ;
3900 oz = ele->z ;
3901
3902 if ( getDataAt(ox,oy,oz) == curwid )
3903 {
3904 ele = queue2->remove() ;
3905 }
3906 else
3907 {
3908 setDataAt(ox,oy,oz, curwid) ;
3909 ele = queue2->getNext() ;
3910 }
3911 }
3912 queue4->reset() ;
3913 ele = queue4->getNext() ;
3914 while ( ele != NULL )
3915 {
3916 ox = ele->x ;
3917 oy = ele->y ;
3918 oz = ele->z ;
3919
3920 queue2->prepend(ox,oy,oz) ;
3921 ele = queue4->remove() ;
3922 }
3923
3924 // Now queue2 holds all the nodes for this layer
3925
3926 #ifdef NOISE_DIS_HELIX
3927 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
3928 queue2->reset() ;
3929
3930 // First run
3931 int flag = 0 ;
3932 while ( ( ele = queue2->getNext() ) != NULL )
3933 {
3934 ox = ele->x ;
3935 oy = ele->y ;
3936 oz = ele->z ;
3937 if ( NOISE_DIS_HELIX <= 1 )
3938 {
3939 noisevol->setDataAt( ox, oy, oz, 0 ) ;
3940 }
3941 else
3942 {
3943 flag = 0 ;
3944 for ( int m = 0 ; m < 6 ; m ++ )
3945 {
3946 int nx = ox + neighbor6[m][0] ;
3947 int ny = oy + neighbor6[m][1] ;
3948 int nz = oz + neighbor6[m][2] ;
3949 if ( getDataAt( nx, ny, nz ) == 0 )
3950 {
3951 noisevol->setDataAt( ox, oy, oz, 1 ) ;
3952 flag = 1 ;
3953 break ;
3954 }
3955 }
3956 if ( ! flag )
3957 {
3958 noisevol->setDataAt( ox, oy, oz, 0 ) ;
3959 }
3960 }
3961 }
3962
3963 int cur, visited ;
3964 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
3965 {
3966 queue2->reset() ;
3967 int count = 0 ;
3968 visited = 0 ;
3969
3970 while ( ( ele = queue2->getNext() ) != NULL )
3971 {
3972 ox = ele->x ;
3973 oy = ele->y ;
3974 oz = ele->z ;
3975
3976 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
3977 {
3978 visited ++ ;
3979 continue ;
3980 }
3981
3982 flag = 0 ;
3983 for ( int m = 0 ; m < 6 ; m ++ )
3984 {
3985 int nx = ox + neighbor6[m][0] ;
3986 int ny = oy + neighbor6[m][1] ;
3987 int nz = oz + neighbor6[m][2] ;
3988 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
3989 {
3990 noisevol->setDataAt( ox, oy, oz, 1 ) ;
3991 visited ++ ;
3992 count ++ ;
3993 break ;
3994 }
3995 }
3996 }
3997
3998 if ( count == 0 )
3999 {
4000 break ;
4001 }
4002 }
4003 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
4004
4005
4006 #endif
4007 /* Commented out for debugging
4008
4009 // First,
4010 // check for complex nodes in queue2
4011 // move them from queue2 to queue3
4012 queue2->reset() ;
4013 ele = queue2->getNext() ;
4014 while ( ele != NULL )
4015 {
4016 ox = ele->x ;
4017 oy = ele->y ;
4018 oz = ele->z ;
4019
4020 // Check simple
4021 #ifndef NOISE_DIS_HELIX
4022 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4023 #else
4024 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
4025 #endif
4026 {
4027 // Complex, set to next layer
4028 setDataAt( ox, oy, oz, curwid + 1 ) ;
4029 queue3->prepend( ox, oy, oz ) ;
4030 ele = queue2->remove() ;
4031
4032 numComplex ++ ;
4033 }
4034 else
4035 {
4036 ele = queue2->getNext() ;
4037 }
4038 }
4039 */
4040
4041 // Next,
4042 // Compute score for each node left in queue2
4043 // move them into priority queue
4044 queue2->reset() ;
4045 ele = queue2->getNext() ;
4046 while ( ele != NULL )
4047 {
4048 ox = ele->x ;
4049 oy = ele->y ;
4050 oz = ele->z ;
4051
4052 // Compute score
4053 score = getNumPotComplex2( ox, oy, oz ) ;
4054 scrvol->setDataAt( ox, oy, oz, score ) ;
4055
4056 // Push to queue
4057 gp = new gridPoint ;
4058 gp->x = ox ;
4059 gp->y = oy ;
4060 gp->z = oz ;
4061 // queue->add( gp, -score ) ;
4062 queue->add( gp, score ) ;
4063
4064 ele = queue2->remove() ;
4065 }
4066
4067 // Rename queue3 to be queue2,
4068 // Clear queue3
4069 // From now on, queue2 holds nodes of next level
4070 delete queue2 ;
4071 queue2 = queue3 ;
4072 queue3 = new GridQueue2( ) ;
4073
4074 // Next, start priority queue iteration
4075 while ( ! queue->isEmpty() )
4076 {
4077 // Retrieve the node with the highest score
4078 queue->remove( gp, score ) ;
4079 ox = gp->x ;
4080 oy = gp->y ;
4081 oz = gp->z ;
4082 delete gp ;
4083 // score = -score ;
4084
4085 // Ignore the node
4086 // if it has been processed before
4087 // or it has an updated score
4088 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
4089 {
4090 continue ;
4091 }
4092
4093 /* Commented out for debugging
4094
4095 // Remove this simple node
4096 setDataAt( ox, oy, oz, -1 ) ;
4097 numSimple ++ ;
4098 // printf("Highest score: %d\n", score) ;
4099 */
4100
4101 /* Added for debugging */
4102 // Check simple
4103 #ifndef NOISE_DIS_HELIX
4104 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
4105 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4106 #else
4107 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4108 #endif
4109 {
4110 // Complex, set to next layer
4111 setDataAt( ox, oy, oz, curwid + 1 ) ;
4112 queue4->prepend( ox, oy, oz ) ;
4113 numComplex ++ ;
4114 }
4115 else
4116 {
4117 setDataAt( ox, oy, oz, -1 ) ;
4118 numSimple ++ ;
4119
4120
4121 /* Adding ends */
4122 // Move its neighboring unvisited node to queue2
4123 for ( int m = 0 ; m < 6 ; m ++ )
4124 {
4125 int nx = ox + neighbor6[m][0] ;
4126 int ny = oy + neighbor6[m][1] ;
4127 int nz = oz + neighbor6[m][2] ;
4128 if ( getDataAt( nx, ny, nz ) == 0 )
4129 {
4130 // setDataAt( nx, ny, nz, curwid + 1 ) ;
4131 queue2->prepend( nx, ny, nz ) ;
4132 }
4133 }
4134
4135 }
4136
4137
4138 /* Commented out for debugging
4139
4140 // Find complex nodes in its 3x3 neighborhood
4141 // move them to queue2
4142 for ( i = -1 ; i < 2 ; i ++ )
4143 for ( j = -1 ; j < 2 ; j ++ )
4144 for ( k = -1 ; k < 2 ; k ++ )
4145 {
4146 int nx = ox + i ;
4147 int ny = oy + j ;
4148 int nz = oz + k ;
4149
4150 // Check simple
4151 if ( getDataAt( nx, ny, nz ) == curwid &&
4152 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
4153 #ifndef NOISE_DIS_HELIX
4154 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
4155 #else
4156 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
4157 #endif
4158
4159 {
4160 // Complex, set to next layer
4161 setDataAt( nx, ny, nz, curwid + 1 ) ;
4162 queue2->prepend( nx, ny, nz ) ;
4163 numComplex ++ ;
4164 }
4165 }
4166 */
4167
4168 // Update scores for nodes in its 5x5 neighborhood
4169 // insert them back into priority queue
4170
4171 for ( i = -2 ; i < 3 ;i ++ )
4172 for ( j = -2 ; j < 3 ; j ++ )
4173 for ( k = -2 ; k < 3 ; k ++ )
4174 {
4175 int nx = ox + i ;
4176 int ny = oy + j ;
4177 int nz = oz + k ;
4178
4179 if ( getDataAt( nx, ny, nz ) == curwid )
4180 {
4181 // Compute score
4182 score = getNumPotComplex2( nx, ny, nz ) ;
4183
4184 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
4185 {
4186 // printf("Update\n") ;
4187 scrvol->setDataAt( nx, ny, nz, score ) ;
4188 // Push to queue
4189 gp = new gridPoint ;
4190 gp->x = nx ;
4191 gp->y = ny ;
4192 gp->z = nz ;
4193 // queue->add( gp, -score ) ;
4194 queue->add( gp, score ) ;
4195 }
4196 }
4197 }
4198
4199
4200 }
4201
4202 #ifdef VERBOSE
4203 printf("%d complex, %d simple\n", numComplex, numSimple) ;
4204 #endif
4205
4206 if ( numSimple == 0 )
4207 {
4208 if ( queue2->getNumElements() > 0 )
4209 {
4210 printf("*************************wierd here*************************\n");
4211 }
4212 break ;
4213 }
4214 }
4215
4216 // Remove all internal voxels (contained in manifold surfaces)
4217 queue2->reset() ;
4218 queue4->reset() ;
4219 ele = queue4->getNext() ;
4220 while ( ele != NULL )
4221 {
4222 ox = ele->x ;
4223 oy = ele->y ;
4224 oz = ele->z ;
4225
4226 if ( isPiercable(ox,oy,oz) == 1 ) // hasCompleteSheet( ox, oy, oz ) == 1 ) //
4227 {
4228 queue2->prepend(ox,oy,oz) ;
4229 // setDataAt( ox, oy, oz, -1 ) ;
4230 }
4231 ele = queue4->remove() ;
4232 }
4233
4234 for ( i = 0 ; i < getSizeX() ; i ++ )
4235 for ( j = 0 ; j < getSizeY() ; j ++ )
4236 for ( k = 0 ; k < getSizeZ() ; k ++ )
4237 {
4238 if ( getDataAt( i, j, k ) == 0 && isPiercable(i,j,k) ) //hasCompleteSheet(i,j,k) == 1) //
4239 {
4240 queue2->prepend( i, j, k ) ;
4241 }
4242 }
4243 queue2->reset() ;
4244 ele = queue2->getNext() ;
4245 while ( ele != NULL )
4246 {
4247 ox = ele->x ;
4248 oy = ele->y ;
4249 oz = ele->z ;
4250 setDataAt( ox, oy, oz, -1 ) ;
4251 ele = queue2->remove() ;
4252 }
4253
4254
4255 // Finally, clean up
4256 delete scrvol;
4257 delete queue;
4258 delete queue2;
4259 delete queue3;
4260 delete queue4;
4261 #ifdef VERBOSE
4262 printf("Thresholding the volume to 0/1...\n") ;
4263 #endif
4264 threshold( 0, 0, 1 ) ;
4265 }
4266
4267 // Compute curve skeleton
4268 void Volume::curveSkeleton( float thr, Volume* svol )
4269 {
4270 int i, j, k ;
4271 // First, threshold the volume
4272 #ifdef VERBOSE
4273 printf("Thresholding the volume to -1/0...\n") ;
4274 #endif
4275 threshold( thr, -1, 0 ) ;
4276
4277 // Next, apply convergent erosion
4278 // by preserving: complex nodes, curve end-points, and sheet points
4279
4280 // Next, initialize the linked queue
4281 #ifdef VERBOSE
4282 printf("Initializing queue...\n") ;
4283 #endif
4284 GridQueue2* queue2 = new GridQueue2( ) ;
4285 GridQueue2* queue3 = new GridQueue2( ) ;
4286 GridQueue2* queue4 = new GridQueue2( ) ;
4287 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
4288
4289 for ( i = 0 ; i < getSizeX() ; i ++ )
4290 for ( j = 0 ; j < getSizeY() ; j ++ )
4291 for ( k = 0 ; k < getSizeZ() ; k ++ )
4292 {
4293 if ( getDataAt( i, j, k ) >= 0 )
4294 {
4295 if ( svol->getDataAt(i,j,k) > 0 )
4296 {
4297 setDataAt( i, j, k, MAX_ERODE ) ;
4298 }
4299 else
4300 {
4301 for ( int m = 0 ; m < 6 ; m ++ )
4302 {
4303 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
4304 {
4305 // setDataAt( i, j, k, 1 ) ;
4306 queue2->prepend( i, j, k ) ;
4307 break ;
4308 }
4309 }
4310 }
4311 }
4312 }
4313
4314 int wid = MAX_ERODE ;
4315 #ifdef VERBOSE
4316 printf("Total %d nodes\n", queue2->getNumElements() ) ;
4317 printf("Start erosion to %d...\n", wid) ;
4318 #endif
4319
4320
4321 // Perform erosion
4322 gridQueueEle* ele ;
4323 gridPoint* gp ;
4324 int ox, oy, oz ;
4325 int score ;
4326 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
4327 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
4328 {
4329 scrvol->setDataAt( i, -1 ) ;
4330 }
4331
4332 #ifdef NOISE_DIS_HELIX
4333 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
4334 #endif
4335
4336 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
4337 {
4338 // At the start of each iteration,
4339 // queue2 holds all the nodes for this layer
4340 // queue3 and queue are empty
4341
4342 int numComplex = 0, numSimple = 0 ;
4343 #ifdef VERBOSE
4344 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
4345 #endif
4346
4347 /*
4348 We first need to assign curwid + 1 to every node in this layer
4349 */
4350 queue2->reset() ;
4351 ele = queue2->getNext() ;
4352 while ( ele != NULL )
4353 {
4354 ox = ele->x ;
4355 oy = ele->y ;
4356 oz = ele->z ;
4357
4358 if ( getDataAt(ox,oy,oz) == curwid )
4359 {
4360 ele = queue2->remove() ;
4361 }
4362 else
4363 {
4364 setDataAt(ox,oy,oz, curwid) ;
4365 ele = queue2->getNext() ;
4366 }
4367 }
4368 queue4->reset() ;
4369 ele = queue4->getNext() ;
4370 while ( ele != NULL )
4371 {
4372 ox = ele->x ;
4373 oy = ele->y ;
4374 oz = ele->z ;
4375
4376 queue2->prepend(ox,oy,oz) ;
4377 ele = queue4->remove() ;
4378 }
4379
4380 // Now queue2 holds all the nodes for this layer
4381
4382 #ifdef NOISE_DIS_HELIX
4383 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
4384 queue2->reset() ;
4385
4386 // First run
4387 int flag = 0 ;
4388 while ( ( ele = queue2->getNext() ) != NULL )
4389 {
4390 ox = ele->x ;
4391 oy = ele->y ;
4392 oz = ele->z ;
4393 if ( NOISE_DIS_HELIX <= 1 )
4394 {
4395 noisevol->setDataAt( ox, oy, oz, 0 ) ;
4396 }
4397 else
4398 {
4399 flag = 0 ;
4400 for ( int m = 0 ; m < 6 ; m ++ )
4401 {
4402 int nx = ox + neighbor6[m][0] ;
4403 int ny = oy + neighbor6[m][1] ;
4404 int nz = oz + neighbor6[m][2] ;
4405 if ( getDataAt( nx, ny, nz ) == 0 )
4406 {
4407 noisevol->setDataAt( ox, oy, oz, 1 ) ;
4408 flag = 1 ;
4409 break ;
4410 }
4411 }
4412 if ( ! flag )
4413 {
4414 noisevol->setDataAt( ox, oy, oz, 0 ) ;
4415 }
4416 }
4417 }
4418
4419 int cur, visited ;
4420 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
4421 {
4422 queue2->reset() ;
4423 int count = 0 ;
4424 visited = 0 ;
4425
4426 while ( ( ele = queue2->getNext() ) != NULL )
4427 {
4428 ox = ele->x ;
4429 oy = ele->y ;
4430 oz = ele->z ;
4431
4432 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
4433 {
4434 visited ++ ;
4435 continue ;
4436 }
4437
4438 flag = 0 ;
4439 for ( int m = 0 ; m < 6 ; m ++ )
4440 {
4441 int nx = ox + neighbor6[m][0] ;
4442 int ny = oy + neighbor6[m][1] ;
4443 int nz = oz + neighbor6[m][2] ;
4444 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
4445 {
4446 noisevol->setDataAt( ox, oy, oz, 1 ) ;
4447 visited ++ ;
4448 count ++ ;
4449 break ;
4450 }
4451 }
4452 }
4453
4454 if ( count == 0 )
4455 {
4456 break ;
4457 }
4458 }
4459 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
4460
4461
4462 #endif
4463 /* Commented out for debugging
4464
4465 // First,
4466 // check for complex nodes in queue2
4467 // move them from queue2 to queue3
4468 queue2->reset() ;
4469 ele = queue2->getNext() ;
4470 while ( ele != NULL )
4471 {
4472 ox = ele->x ;
4473 oy = ele->y ;
4474 oz = ele->z ;
4475
4476 // Check simple
4477 #ifndef NOISE_DIS_HELIX
4478 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4479 #else
4480 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
4481 #endif
4482 {
4483 // Complex, set to next layer
4484 setDataAt( ox, oy, oz, curwid + 1 ) ;
4485 queue3->prepend( ox, oy, oz ) ;
4486 ele = queue2->remove() ;
4487
4488 numComplex ++ ;
4489 }
4490 else
4491 {
4492 ele = queue2->getNext() ;
4493 }
4494 }
4495 */
4496
4497 // Next,
4498 // Compute score for each node left in queue2
4499 // move them into priority queue
4500 queue2->reset() ;
4501 ele = queue2->getNext() ;
4502 while ( ele != NULL )
4503 {
4504 ox = ele->x ;
4505 oy = ele->y ;
4506 oz = ele->z ;
4507
4508 // Compute score
4509 score = getNumPotComplex2( ox, oy, oz ) ;
4510 scrvol->setDataAt( ox, oy, oz, score ) ;
4511
4512 // Push to queue
4513 gp = new gridPoint ;
4514 gp->x = ox ;
4515 gp->y = oy ;
4516 gp->z = oz ;
4517 // queue->add( gp, -score ) ;
4518 queue->add( gp, score ) ;
4519
4520 ele = queue2->remove() ;
4521 }
4522
4523 // Rename queue3 to be queue2,
4524 // Clear queue3
4525 // From now on, queue2 holds nodes of next level
4526 delete queue2 ;
4527 queue2 = queue3 ;
4528 queue3 = new GridQueue2( ) ;
4529
4530 // Next, start priority queue iteration
4531 while ( ! queue->isEmpty() )
4532 {
4533 // Retrieve the node with the highest score
4534 queue->remove( gp, score ) ;
4535 ox = gp->x ;
4536 oy = gp->y ;
4537 oz = gp->z ;
4538 delete gp ;
4539 // score = -score ;
4540
4541 // Ignore the node
4542 // if it has been processed before
4543 // or it has an updated score
4544 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
4545 {
4546 continue ;
4547 }
4548
4549 /* Commented out for debugging
4550
4551 // Remove this simple node
4552 setDataAt( ox, oy, oz, -1 ) ;
4553 numSimple ++ ;
4554 // printf("Highest score: %d\n", score) ;
4555 */
4556
4557 /* Added for debugging */
4558 // Check simple
4559 #ifndef NOISE_DIS_HELIX
4560 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
4561 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4562 #else
4563 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4564 #endif
4565 {
4566 // Complex, set to next layer
4567 setDataAt( ox, oy, oz, curwid + 1 ) ;
4568 queue4->prepend( ox, oy, oz ) ;
4569 numComplex ++ ;
4570 }
4571 else
4572 {
4573 setDataAt( ox, oy, oz, -1 ) ;
4574 numSimple ++ ;
4575 }
4576 /* Adding ends */
4577
4578 // Move its neighboring unvisited node to queue2
4579 for ( int m = 0 ; m < 6 ; m ++ )
4580 {
4581 int nx = ox + neighbor6[m][0] ;
4582 int ny = oy + neighbor6[m][1] ;
4583 int nz = oz + neighbor6[m][2] ;
4584 if ( getDataAt( nx, ny, nz ) == 0 )
4585 {
4586 // setDataAt( nx, ny, nz, curwid + 1 ) ;
4587 queue2->prepend( nx, ny, nz ) ;
4588 }
4589 }
4590
4591 /* Commented out for debugging
4592
4593 // Find complex nodes in its 3x3 neighborhood
4594 // move them to queue2
4595 for ( i = -1 ; i < 2 ; i ++ )
4596 for ( j = -1 ; j < 2 ; j ++ )
4597 for ( k = -1 ; k < 2 ; k ++ )
4598 {
4599 int nx = ox + i ;
4600 int ny = oy + j ;
4601 int nz = oz + k ;
4602
4603 // Check simple
4604 if ( getDataAt( nx, ny, nz ) == curwid &&
4605 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
4606 #ifndef NOISE_DIS_HELIX
4607 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
4608 #else
4609 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
4610 #endif
4611
4612 {
4613 // Complex, set to next layer
4614 setDataAt( nx, ny, nz, curwid + 1 ) ;
4615 queue2->prepend( nx, ny, nz ) ;
4616 numComplex ++ ;
4617 }
4618 }
4619 */
4620
4621 // Update scores for nodes in its 5x5 neighborhood
4622 // insert them back into priority queue
4623
4624 for ( i = -2 ; i < 3 ;i ++ )
4625 for ( j = -2 ; j < 3 ; j ++ )
4626 for ( k = -2 ; k < 3 ; k ++ )
4627 {
4628 int nx = ox + i ;
4629 int ny = oy + j ;
4630 int nz = oz + k ;
4631
4632 if ( getDataAt( nx, ny, nz ) == curwid )
4633 {
4634 // Compute score
4635 score = getNumPotComplex2( nx, ny, nz ) ;
4636
4637 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
4638 {
4639 // printf("Update\n") ;
4640 scrvol->setDataAt( nx, ny, nz, score ) ;
4641 // Push to queue
4642 gp = new gridPoint ;
4643 gp->x = nx ;
4644 gp->y = ny ;
4645 gp->z = nz ;
4646 // queue->add( gp, -score ) ;
4647 queue->add( gp, score ) ;
4648 }
4649 }
4650 }
4651
4652
4653 }
4654
4655 #ifdef VERBOSE
4656 printf("%d complex, %d simple\n", numComplex, numSimple) ;
4657 #endif
4658
4659 if ( numSimple == 0 )
4660 {
4661 break ;
4662 }
4663 }
4664
4665 // Finally, clean up
4666 delete scrvol;
4667 delete queue;
4668 delete queue2;
4669 delete queue3;
4670 delete queue4;
4671 #ifdef VERBOSE
4672 printf("Thresholding the volume to 0/1...\n") ;
4673 #endif
4674 threshold( 0, 0, 1 ) ;
4675 }
4676
4677 // Compute curve skeleton in 2D
4678 void Volume::curveSkeleton2D( float thr, Volume* svol )
4679 {
4680 int i, j, k ;
4681 // First, threshold the volume
4682 #ifdef VERBOSE
4683 printf("Thresholding the volume to -1/0...\n") ;
4684 #endif
4685 threshold( thr, -1, 0 ) ;
4686
4687 // Next, apply convergent erosion
4688 // by preserving: complex nodes, curve end-points, and sheet points
4689
4690 // Next, initialize the linked queue
4691 #ifdef VERBOSE
4692 printf("Initializing queue...\n") ;
4693 #endif
4694 GridQueue2* queue2 = new GridQueue2( ) ;
4695 GridQueue2* queue3 = new GridQueue2( ) ;
4696 GridQueue2* queue4 = new GridQueue2( ) ;
4697 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
4698
4699 for ( i = 0 ; i < getSizeX() ; i ++ )
4700 for ( j = 0 ; j < getSizeY() ; j ++ )
4701 for ( k = 0 ; k < getSizeZ() ; k ++ )
4702 {
4703 if ( getDataAt( i, j, k ) >= 0 )
4704 {
4705 if ( svol->getDataAt(i,j,k) > 0 )
4706 {
4707 setDataAt( i, j, k, MAX_ERODE ) ;
4708 }
4709 else
4710 {
4711 for ( int m = 0 ; m < 4 ; m ++ )
4712 {
4713 if ( getDataAt( i + neighbor4[m][0], j + neighbor4[m][1], k ) < 0 )
4714 {
4715 // setDataAt( i, j, k, 1 ) ;
4716 queue2->prepend( i, j, k ) ;
4717 break ;
4718 }
4719 }
4720 }
4721 }
4722 }
4723 int wid = MAX_ERODE ;
4724 #ifdef VERBOSE
4725 printf("Total %d nodes\n", queue2->getNumElements() ) ;
4726 printf("Start erosion to %d...\n", wid) ;
4727 #endif
4728
4729
4730 // Perform erosion
4731 gridQueueEle* ele ;
4732 gridPoint* gp ;
4733 int ox, oy, oz ;
4734 int score ;
4735 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
4736 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
4737 {
4738 scrvol->setDataAt( i, -1 ) ;
4739 }
4740
4741 #ifdef NOISE_DIS_HELIX
4742 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
4743 #endif
4744
4745 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
4746 {
4747 // At the start of each iteration,
4748 // queue2 holds all the nodes for this layer
4749 // queue3 and queue are empty
4750
4751 int numComplex = 0, numSimple = 0 ;
4752 #ifdef VERBOSE
4753 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
4754 #endif
4755
4756 /*
4757 We first need to assign curwid + 1 to every node in this layer
4758 */
4759 queue2->reset() ;
4760 ele = queue2->getNext() ;
4761 while ( ele != NULL )
4762 {
4763 ox = ele->x ;
4764 oy = ele->y ;
4765 oz = ele->z ;
4766
4767 if ( getDataAt(ox,oy,oz) == curwid )
4768 {
4769 ele = queue2->remove() ;
4770 }
4771 else
4772 {
4773 setDataAt(ox,oy,oz, curwid) ;
4774 ele = queue2->getNext() ;
4775 }
4776 }
4777 queue4->reset() ;
4778 ele = queue4->getNext() ;
4779 while ( ele != NULL )
4780 {
4781 ox = ele->x ;
4782 oy = ele->y ;
4783 oz = ele->z ;
4784
4785 queue2->prepend(ox,oy,oz) ;
4786 ele = queue4->remove() ;
4787 }
4788
4789 // Now queue2 holds all the nodes for this layer
4790
4791 #ifdef NOISE_DIS_HELIX
4792 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */
4793 queue2->reset() ;
4794
4795 // First run
4796 int flag = 0 ;
4797 while ( ( ele = queue2->getNext() ) != NULL )
4798 {
4799 ox = ele->x ;
4800 oy = ele->y ;
4801 oz = ele->z ;
4802 if ( NOISE_DIS_HELIX <= 1 )
4803 {
4804 noisevol->setDataAt( ox, oy, oz, 0 ) ;
4805 }
4806 else
4807 {
4808 flag = 0 ;
4809 for ( int m = 0 ; m < 6 ; m ++ )
4810 {
4811 int nx = ox + neighbor6[m][0] ;
4812 int ny = oy + neighbor6[m][1] ;
4813 int nz = oz + neighbor6[m][2] ;
4814 if ( getDataAt( nx, ny, nz ) == 0 )
4815 {
4816 noisevol->setDataAt( ox, oy, oz, 1 ) ;
4817 flag = 1 ;
4818 break ;
4819 }
4820 }
4821 if ( ! flag )
4822 {
4823 noisevol->setDataAt( ox, oy, oz, 0 ) ;
4824 }
4825 }
4826 }
4827
4828 int cur, visited ;
4829 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ )
4830 {
4831 queue2->reset() ;
4832 int count = 0 ;
4833 visited = 0 ;
4834
4835 while ( ( ele = queue2->getNext() ) != NULL )
4836 {
4837 ox = ele->x ;
4838 oy = ele->y ;
4839 oz = ele->z ;
4840
4841 if ( noisevol->getDataAt( ox, oy, oz ) == 1 )
4842 {
4843 visited ++ ;
4844 continue ;
4845 }
4846
4847 flag = 0 ;
4848 for ( int m = 0 ; m < 6 ; m ++ )
4849 {
4850 int nx = ox + neighbor6[m][0] ;
4851 int ny = oy + neighbor6[m][1] ;
4852 int nz = oz + neighbor6[m][2] ;
4853 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 )
4854 {
4855 noisevol->setDataAt( ox, oy, oz, 1 ) ;
4856 visited ++ ;
4857 count ++ ;
4858 break ;
4859 }
4860 }
4861 }
4862
4863 if ( count == 0 )
4864 {
4865 break ;
4866 }
4867 }
4868 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ;
4869
4870
4871 #endif
4872 /* Commented out for debugging
4873
4874 // First,
4875 // check for complex nodes in queue2
4876 // move them from queue2 to queue3
4877 queue2->reset() ;
4878 ele = queue2->getNext() ;
4879 while ( ele != NULL )
4880 {
4881 ox = ele->x ;
4882 oy = ele->y ;
4883 oz = ele->z ;
4884
4885 // Check simple
4886 #ifndef NOISE_DIS_HELIX
4887 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4888 #else
4889 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) )
4890 #endif
4891 {
4892 // Complex, set to next layer
4893 setDataAt( ox, oy, oz, curwid + 1 ) ;
4894 queue3->prepend( ox, oy, oz ) ;
4895 ele = queue2->remove() ;
4896
4897 numComplex ++ ;
4898 }
4899 else
4900 {
4901 ele = queue2->getNext() ;
4902 }
4903 }
4904 */
4905
4906 // Next,
4907 // Compute score for each node left in queue2
4908 // move them into priority queue
4909 queue2->reset() ;
4910 ele = queue2->getNext() ;
4911 while ( ele != NULL )
4912 {
4913 ox = ele->x ;
4914 oy = ele->y ;
4915 oz = ele->z ;
4916
4917 // Compute score
4918 score = getNumPotComplex2( ox, oy, oz ) ;
4919 //score = getNumNeighbor6( ox, oy, oz ) ;
4920 scrvol->setDataAt( ox, oy, oz, score ) ;
4921
4922 // Push to queue
4923 gp = new gridPoint ;
4924 gp->x = ox ;
4925 gp->y = oy ;
4926 gp->z = oz ;
4927 // queue->add( gp, -score ) ;
4928 queue->add( gp, score ) ;
4929
4930 ele = queue2->remove() ;
4931 }
4932
4933 // Rename queue3 to be queue2,
4934 // Clear queue3
4935 // From now on, queue2 holds nodes of next level
4936 delete queue2 ;
4937 queue2 = queue3 ;
4938 queue3 = new GridQueue2( ) ;
4939
4940 // Next, start priority queue iteration
4941 while ( ! queue->isEmpty() )
4942 {
4943 // Retrieve the node with the highest score
4944 queue->remove( gp, score ) ;
4945 ox = gp->x ;
4946 oy = gp->y ;
4947 oz = gp->z ;
4948 delete gp ;
4949 // score = -score ;
4950
4951 // Ignore the node
4952 // if it has been processed before
4953 // or it has an updated score
4954 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
4955 {
4956 continue ;
4957 }
4958
4959 /* Commented out for debugging
4960
4961 // Remove this simple node
4962 setDataAt( ox, oy, oz, -1 ) ;
4963 numSimple ++ ;
4964 // printf("Highest score: %d\n", score) ;
4965 */
4966
4967 /* Added for debugging */
4968 // Check simple
4969 #ifndef NOISE_DIS_HELIX
4970 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) )
4971 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4972 #else
4973 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) )
4974 #endif
4975 {
4976 // Complex, set to next layer
4977 setDataAt( ox, oy, oz, curwid + 1 ) ;
4978 queue4->prepend( ox, oy, oz ) ;
4979 numComplex ++ ;
4980 }
4981 else
4982 {
4983 setDataAt( ox, oy, oz, -1 ) ;
4984 numSimple ++ ;
4985 }
4986 /* Adding ends */
4987
4988 // Move its neighboring unvisited node to queue2
4989 for ( int m = 0 ; m < 4 ; m ++ )
4990 {
4991 int nx = ox + neighbor4[m][0] ;
4992 int ny = oy + neighbor4[m][1] ;
4993 int nz = oz ;
4994 if ( getDataAt( nx, ny, nz ) == 0 )
4995 {
4996 // setDataAt( nx, ny, nz, curwid + 1 ) ;
4997 queue2->prepend( nx, ny, nz ) ;
4998 }
4999 }
5000
5001 /* Commented out for debugging
5002
5003 // Find complex nodes in its 3x3 neighborhood
5004 // move them to queue2
5005 for ( i = -1 ; i < 2 ; i ++ )
5006 for ( j = -1 ; j < 2 ; j ++ )
5007 for ( k = -1 ; k < 2 ; k ++ )
5008 {
5009 int nx = ox + i ;
5010 int ny = oy + j ;
5011 int nz = oz + k ;
5012
5013 // Check simple
5014 if ( getDataAt( nx, ny, nz ) == curwid &&
5015 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) )
5016 #ifndef NOISE_DIS_HELIX
5017 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) )
5018 #else
5019 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) )
5020 #endif
5021
5022 {
5023 // Complex, set to next layer
5024 setDataAt( nx, ny, nz, curwid + 1 ) ;
5025 queue2->prepend( nx, ny, nz ) ;
5026 numComplex ++ ;
5027 }
5028 }
5029 */
5030
5031 // Update scores for nodes in its 5x5 neighborhood
5032 // insert them back into priority queue
5033
5034 for ( i = -2 ; i < 3 ;i ++ )
5035 for ( j = -2 ; j < 3 ; j ++ )
5036 for ( k = -2 ; k < 3 ; k ++ )
5037 {
5038 int nx = ox + i ;
5039 int ny = oy + j ;
5040 int nz = oz + k ;
5041
5042 if ( getDataAt( nx, ny, nz ) == curwid )
5043 {
5044 // Compute score
5045 score = getNumPotComplex2( nx, ny, nz ) ;
5046 //score = getNumNeighbor6( nx, ny, nz ) ;
5047
5048 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
5049 {
5050 // printf("Update\n") ;
5051 scrvol->setDataAt( nx, ny, nz, score ) ;
5052 // Push to queue
5053 gp = new gridPoint ;
5054 gp->x = nx ;
5055 gp->y = ny ;
5056 gp->z = nz ;
5057 // queue->add( gp, -score ) ;
5058 queue->add( gp, score ) ;
5059 }
5060 }
5061 }
5062
5063
5064 }
5065
5066 #ifdef VERBOSE
5067 printf("%d complex, %d simple\n", numComplex, numSimple) ;
5068 #endif
5069
5070 if ( numSimple == 0 )
5071 {
5072 break ;
5073 }
5074 }
5075
5076 // Finally, clean up
5077 #ifdef VERBOSE
5078 printf("Thresholding the volume to 0/1...\n") ;
5079 #endif
5080 threshold( 0, 0, 1 ) ;
5081 delete scrvol;
5082 delete queue;
5083 delete queue2;
5084 delete queue3;
5085 delete queue4;
5086 }
5087
5088 // Compute minimal skeleton
5089 void Volume::skeleton( float thr, int )
5090 {
5091 int i, j, k ;
5092 // First, threshold the volume
5093 #ifdef VERBOSE
5094 printf("Thresholding the volume to -1/0...\n") ;
5095 #endif
5096 threshold( thr, -1, 0 ) ;
5097
5098 // Next, apply convergent erosion
5099 // by preserving: complex nodes, curve end-points, and sheet points
5100
5101 // Next, initialize the linked queue
5102 #ifdef VERBOSE
5103 printf("Initializing queue...\n") ;
5104 #endif
5105 GridQueue2* queue2 = new GridQueue2( ) ;
5106 GridQueue2* queue = new GridQueue2( ) ;
5107
5108 for ( i = 0 ; i < getSizeX() ; i ++ )
5109 for ( j = 0 ; j < getSizeY() ; j ++ )
5110 for ( k = 0 ; k < getSizeZ() ; k ++ )
5111 {
5112 if ( getDataAt( i, j, k ) >= 0 )
5113 {
5114 {
5115 for ( int m = 0 ; m < 6 ; m ++ )
5116 {
5117 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
5118 {
5119 setDataAt( i, j, k, 1 ) ;
5120 queue2->prepend( i, j, k ) ;
5121 break ;
5122 }
5123 }
5124 }
5125 }
5126 }
5127 int wid = 0 ;
5128 #ifdef VERBOSE
5129 printf("Total %d nodes\n", queue2->getNumElements() ) ;
5130 printf("Start erosion to %d...\n", wid) ;
5131 #endif
5132
5133 // Perform erosion
5134 gridQueueEle* ele ;
5135 int ox, oy, oz ;
5136
5137 while( 1 )
5138 {
5139 wid ++ ;
5140
5141 // At the start of each iteration,
5142 // queue2 holds all the nodes for this layer
5143 // queue is empty
5144
5145 int numComplex = 0, numSimple = 0 ;
5146 #ifdef VERBOSE
5147 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), wid) ;
5148 #endif
5149
5150 // Rename queue2 to be queue,
5151 // Clear queue2
5152 // From now on, queue2 holds nodes of next level
5153 delete queue ;
5154 queue = queue2 ;
5155 queue2 = new GridQueue2( ) ;
5156
5157 // Next, start queue iteration
5158 queue->reset() ;
5159 ele = queue->getNext();
5160 while ( ele != NULL )
5161 {
5162 ox = ele->x ;
5163 oy = ele->y ;
5164 oz = ele->z ;
5165 // delete ele ;
5166
5167 // Check simple
5168 if ( ! isSimple( ox, oy, oz ) )
5169 {
5170 // Complex, set to next layer
5171 queue2->prepend( ox, oy, oz ) ;
5172 numComplex ++ ;
5173 }
5174 /*
5175 else if ( ox == off || oy == off || oz == off ||
5176 ox == getSizeX() - off - 1 || oy == getSizeY() - off - 1 || oz == getSizeZ() - off - 1 )
5177 {
5178 // Wall, don't erode, set to next layer
5179 queue2->prepend( ox, oy, oz ) ;
5180 numComplex ++ ;
5181 }
5182 */
5183 else
5184 {
5185 setDataAt( ox, oy, oz, -1 ) ;
5186 numSimple ++ ;
5187
5188 for ( int m = 0 ; m < 6 ; m ++ )
5189 {
5190 int nx = ox + neighbor6[m][0] ;
5191 int ny = oy + neighbor6[m][1] ;
5192 int nz = oz + neighbor6[m][2] ;
5193 if ( getDataAt( nx, ny, nz ) == 0 )
5194 {
5195 setDataAt( nx, ny, nz, 1 ) ;
5196 queue2->prepend( nx, ny, nz ) ;
5197 }
5198 }
5199
5200 }
5201
5202 ele = queue->remove() ;
5203 }
5204 #ifdef VERBOSE
5205 printf("Level %d: %d complex, %d simple\n", wid, numComplex, numSimple) ;
5206 #endif
5207
5208 if ( numSimple == 0 )
5209 {
5210 break ;
5211 }
5212 }
5213
5214 // Finally, clean up
5215 delete queue;
5216 delete queue2;
5217 #ifdef VERBOSE
5218 printf("Thresholding the volume to 0/1...\n") ;
5219 #endif
5220 threshold( 0, 0, 1 ) ;
5221 }
5222
5224 //void Volume::skeleton2( float thr, int off )
5225 //{
5226 //int i, j, k ;
5228 //#ifdef VERBOSE
5229 //printf("Thresholding the volume to -1/0...\n") ;
5230 //#endif
5231 //threshold( thr, -1, 0 ) ;
5232
5235
5237 //#ifdef VERBOSE
5238 //printf("Initializing queue...\n") ;
5239 //#endif
5240 //GridQueue2* queue2 = new GridQueue2( ) ;
5241 //GridQueue2* queue3 = new GridQueue2( ) ;
5242 //PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
5243
5244 //for ( i = 0 ; i < getSizeX() ; i ++ )
5245 //for ( j = 0 ; j < getSizeY() ; j ++ )
5246 //for ( k = 0 ; k < getSizeZ() ; k ++ )
5247 //{
5248 //if ( getDataAt( i, j, k ) >= 0 )
5249 //{
5250 //if ( i == off || j == off || k == off ||
5251 //i == getSizeX() - off - 1 || j == getSizeY() - off - 1 || k == getSizeZ() - off - 1 )
5252 //{
5253 //setDataAt( i, j, k, MAX_ERODE ) ;
5254 //}
5255 //else
5256 //{
5257 //for ( int m = 0 ; m < 6 ; m ++ )
5258 //{
5259 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
5260 //{
5261 //setDataAt( i, j, k, 1 ) ;
5262 //queue2->prepend( i, j, k ) ;
5263 //break ;
5264 //}
5265 //}
5266 //}
5267 //}
5268 //}
5269 //int wid = MAX_ERODE ;
5270 //#ifdef VERBOSE
5271 //printf("Total %d nodes\n", queue2->getNumElements() ) ;
5272
5273
5275 //printf("Start erosion to %d...\n", wid) ;
5276 //#endif
5277
5278 //gridQueueEle* ele ;
5279 //gridPoint* gp ;
5280 //int ox, oy, oz ;
5281 //int score ;
5282
5283 //Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
5284 //for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
5285 //{
5286 //scrvol->setDataAt( i, -1 ) ;
5287 //}
5288
5289
5290 //for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
5291 //{
5295
5296 //int numComplex = 0, numSimple = 0 ;
5297 //#ifdef VERBOSE
5298 //printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
5299 //#endif
5300
5304 //queue2->reset() ;
5305 //ele = queue2->getNext() ;
5306 //while ( ele != NULL )
5307 //{
5308 //ox = ele->x ;
5309 //oy = ele->y ;
5310 //oz = ele->z ;
5311
5313 //score = getNumPotComplex2( ox, oy, oz ) ;
5314 //scrvol->setDataAt( ox, oy, oz, score ) ;
5315
5317 //gp = new gridPoint ;
5318 //gp->x = ox ;
5319 //gp->y = oy ;
5320 //gp->z = oz ;
5322 //queue->add( gp, score ) ;
5323
5324 //ele = queue2->remove() ;
5325 //}
5326
5330 //delete queue2 ;
5331 //queue2 = queue3 ;
5332 //queue3 = new GridQueue2( ) ;
5333
5335 //while ( ! queue->isEmpty() )
5336 //{
5338 //queue->remove( gp, score ) ;
5339 //ox = gp->x ;
5340 //oy = gp->y ;
5341 //oz = gp->z ;
5342 //delete gp ;
5344
5348 //if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
5349 //{
5350 //continue ;
5351 //}
5352
5355 //if ( ! isSimple( ox, oy, oz ) )
5356 //{
5358 //setDataAt( ox, oy, oz, curwid + 1 ) ;
5359 //queue2->prepend( ox, oy, oz ) ;
5360 //numComplex ++ ;
5361 //}
5362 //else
5363 //{
5364 //setDataAt( ox, oy, oz, -1 ) ;
5365 //numSimple ++ ;
5366 //}
5368
5370 //for ( int m = 0 ; m < 6 ; m ++ )
5371 //{
5372 //int nx = ox + neighbor6[m][0] ;
5373 //int ny = oy + neighbor6[m][1] ;
5374 //int nz = oz + neighbor6[m][2] ;
5375 //if ( getDataAt( nx, ny, nz ) == 0 )
5376 //{
5377 //setDataAt( nx, ny, nz, curwid + 1 ) ;
5378 //queue2->prepend( nx, ny, nz ) ;
5379 //}
5380 //}
5381
5385 //for ( i = -2 ; i < 3 ;i ++ )
5386 //for ( j = -2 ; j < 3 ; j ++ )
5387 //for ( k = -2 ; k < 3 ; k ++ )
5388 //{
5389 //int nx = ox + i ;
5390 //int ny = oy + j ;
5391 //int nz = oz + k ;
5392
5393 //if ( getDataAt( nx, ny, nz ) == curwid )
5394 //{
5396 //score = getNumPotComplex2( nx, ny, nz ) ;
5397
5398 //if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
5399 //{
5401 //scrvol->setDataAt( nx, ny, nz, score ) ;
5403 //gp = new gridPoint ;
5404 //gp->x = nx ;
5405 //gp->y = ny ;
5406 //gp->z = nz ;
5408 //queue->add( gp, score ) ;
5409 //}
5410 //}
5411 //}
5412
5413 //*/
5414 //}
5415 //#ifdef VERBOSE
5416 //printf("%d complex, %d simple\n", numComplex, numSimple) ;
5417 //#endif
5418
5419 //if ( numSimple == 0 )
5420 //{
5421 //break ;
5422 //}
5423 //}
5424
5426 //delete scrvol ;
5427 //delete queue;
5428 //#ifdef VERBOSE
5429 //printf("Thresholding the volume to 0/1...\n") ;
5430 //#endif
5431 //threshold( 0, 0, 1 ) ;
5432 //}
5433
5435 //* Assuming the current volume has already been thresholded to 0/1
5436 //*/
5437 //void Volume::pointSkeleton( Volume* grayvol, float lowthr, float highthr, Volume* svol, Volume* hvol )
5438 //{
5439 //int i, j, k ;
5441 //#ifdef VERBOSE
5442 //printf("Thresholding the volume to -1/0...\n") ;
5443 //#endif
5444 //threshold( 0.5f, -1, 0 ) ;
5445
5448
5450 //#ifdef VERBOSE
5451 //printf("Initializing queue...\n") ;
5452 //#endif
5453 //GridQueue2* queue2 = new GridQueue2( ) ;
5454 //GridQueue2* queue3 = new GridQueue2( ) ;
5455 //PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
5456
5457 //for ( i = 0 ; i < getSizeX() ; i ++ )
5458 //for ( j = 0 ; j < getSizeY() ; j ++ )
5459 //for ( k = 0 ; k < getSizeZ() ; k ++ )
5460 //{
5461 //if ( getDataAt( i, j, k ) >= 0 )
5462 //{
5463 //float v = (float)grayvol->getDataAt( i, j, k ) ;
5464 //if ( v <= lowthr || v > highthr || svol->getDataAt(i,j,k) > 0 || hvol->getDataAt(i,j,k) > 0 )
5465 //{
5466 //setDataAt( i, j, k, MAX_ERODE ) ;
5467 //}
5468 //else
5469 //{
5470 //for ( int m = 0 ; m < 6 ; m ++ )
5471 //{
5472 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
5473 //{
5474 //setDataAt( i, j, k, 1 ) ;
5475 //queue2->prepend( i, j, k ) ;
5476 //break ;
5477 //}
5478 //}
5479 //}
5480 //}
5481 //}
5482 //#ifdef VERBOSE
5483 //printf("Total %d nodes\n", queue2->getNumElements() ) ;
5484 //#endif
5485
5486
5488 //int wid = MAX_ERODE ;
5489 //#ifdef VERBOSE
5490 //printf("Start erosion to %d...\n", wid) ;
5491 //#endif
5492 //gridQueueEle* ele ;
5493 //gridPoint* gp ;
5494 //int ox, oy, oz ;
5495 //int score ;
5496 //Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
5497 //for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
5498 //{
5499 //scrvol->setDataAt( i, -1 ) ;
5500 //}
5501
5502
5503 //for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
5504 //{
5508
5509 //int numComplex = 0, numSimple = 0 ;
5510 //#ifdef VERBOSE
5511 //printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
5512 //#endif
5513
5514
5518 //queue2->reset() ;
5519 //ele = queue2->getNext() ;
5520 //while ( ele != NULL )
5521 //{
5522 //ox = ele->x ;
5523 //oy = ele->y ;
5524 //oz = ele->z ;
5525
5527 //score = getNumPotComplex2( ox, oy, oz ) ;
5528 //scrvol->setDataAt( ox, oy, oz, score ) ;
5529
5531 //gp = new gridPoint ;
5532 //gp->x = ox ;
5533 //gp->y = oy ;
5534 //gp->z = oz ;
5536 //queue->add( gp, score ) ;
5537
5538 //ele = queue2->remove() ;
5539 //}
5540
5544 //delete queue2 ;
5545 //queue2 = queue3 ;
5546 //queue3 = new GridQueue2( ) ;
5547
5549 //while ( ! queue->isEmpty() )
5550 //{
5552 //queue->remove( gp, score ) ;
5553 //ox = gp->x ;
5554 //oy = gp->y ;
5555 //oz = gp->z ;
5556 //delete gp ;
5558
5562 //if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
5563 //{
5564 //continue ;
5565 //}
5566
5569 //if ( ! isSimple( ox, oy, oz ) )
5570 //{
5572 //setDataAt( ox, oy, oz, curwid + 1 ) ;
5573 //queue2->prepend( ox, oy, oz ) ;
5574 //numComplex ++ ;
5575 //}
5576 //else
5577 //{
5578 //setDataAt( ox, oy, oz, -1 ) ;
5579 //numSimple ++ ;
5581
5583 //for ( int m = 0 ; m < 6 ; m ++ )
5584 //{
5585 //int nx = ox + neighbor6[m][0] ;
5586 //int ny = oy + neighbor6[m][1] ;
5587 //int nz = oz + neighbor6[m][2] ;
5588 //if ( getDataAt( nx, ny, nz ) == 0 )
5589 //{
5590 //setDataAt( nx, ny, nz, curwid + 1 ) ;
5591 //queue2->prepend( nx, ny, nz ) ;
5592 //}
5593 //}
5594
5595 //}
5596
5599
5600 //for ( i = -2 ; i < 3 ;i ++ )
5601 //for ( j = -2 ; j < 3 ; j ++ )
5602 //for ( k = -2 ; k < 3 ; k ++ )
5603 //{
5604 //int nx = ox + i ;
5605 //int ny = oy + j ;
5606 //int nz = oz + k ;
5607
5608 //if ( getDataAt( nx, ny, nz ) == curwid )
5609 //{
5611 //score = getNumPotComplex2( nx, ny, nz ) ;
5612
5613 //if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
5614 //{
5616 //scrvol->setDataAt( nx, ny, nz, score ) ;
5618 //gp = new gridPoint ;
5619 //gp->x = nx ;
5620 //gp->y = ny ;
5621 //gp->z = nz ;
5623 //queue->add( gp, score ) ;
5624 //}
5625 //}
5626 //}
5627
5628
5629 //}
5630 //#ifdef VERBOSE
5631 //printf("%d complex, %d simple\n", numComplex, numSimple) ;
5632 //#endif
5633
5634 //if ( numSimple == 0 )
5635 //{
5636 //break ;
5637 //}
5638 //}
5639
5641 //queue2->reset() ;
5642 //ele = queue2->getNext() ;
5643 //while ( ele != NULL )
5644 //{
5645 //ox = ele->x ;
5646 //oy = ele->y ;
5647 //oz = ele->z ;
5648
5649 //if ( hasCompleteHelix( ox, oy, oz ) == 1 )
5650 //{
5651 //ele = queue2->getNext() ;
5652 //}
5653 //else
5654 //{
5655 //ele = queue2->remove() ;
5656 //}
5657 //}
5658
5659 //for ( i = 0 ; i < getSizeX() ; i ++ )
5660 //for ( j = 0 ; j < getSizeY() ; j ++ )
5661 //for ( k = 0 ; k < getSizeZ() ; k ++ )
5662 //{
5663 //if ( getDataAt( i, j, k ) == 0 && hasCompleteHelix(i,j,k) == 1)
5664 //{
5665 //queue2->prepend( i, j, k ) ;
5666 //}
5667 //}
5668 //queue2->reset() ;
5669 //ele = queue2->getNext() ;
5670 //while ( ele != NULL )
5671 //{
5672 //ox = ele->x ;
5673 //oy = ele->y ;
5674 //oz = ele->z ;
5675 //setDataAt( ox, oy, oz, -1 ) ;
5676 //ele = queue2->remove() ;
5677 //}
5678
5680 //delete scrvol;
5681 //delete queue;
5682 //delete queue2;
5683 //delete queue3;
5684 //#ifdef VERBOSE
5685 //printf("Thresholding the volume to 0/1...\n") ;
5686 //#endif
5687 //threshold( 0, 0, 1 ) ;
5688 //}
5689
5690
5691 // Compute minimal skeleton
5692 void Volume::skeleton( float thr, Volume* svol, Volume* hvol )
5693 {
5694 int i, j, k ;
5695 // First, threshold the volume
5696 #ifdef VERBOSE
5697 printf("Thresholding the volume to -1/0...\n") ;
5698 #endif
5699 threshold( thr, -1, 0 ) ;
5700
5701 // Next, apply convergent erosion
5702 // by preserving: complex nodes, curve end-points, and sheet points
5703
5704 // Next, initialize the linked queue
5705 #ifdef VERBOSE
5706 printf("Initializing queue...\n") ;
5707 #endif
5708 GridQueue2* queue2 = new GridQueue2( ) ;
5709 GridQueue2* queue3 = new GridQueue2( ) ;
5710 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN );
5711
5712 for ( i = 0 ; i < getSizeX() ; i ++ )
5713 for ( j = 0 ; j < getSizeY() ; j ++ )
5714 for ( k = 0 ; k < getSizeZ() ; k ++ )
5715 {
5716 if ( getDataAt( i, j, k ) >= 0 )
5717 {
5718 if ( svol->getDataAt(i,j,k) > 0 || hvol->getDataAt(i,j,k) > 0 )
5719 {
5720 setDataAt( i, j, k, MAX_ERODE ) ;
5721 }
5722 else
5723 {
5724 for ( int m = 0 ; m < 6 ; m ++ )
5725 {
5726 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 )
5727 {
5728 setDataAt( i, j, k, 1 ) ;
5729 queue2->prepend( i, j, k ) ;
5730 break ;
5731 }
5732 }
5733 }
5734 }
5735 }
5736 int wid = MAX_ERODE ;
5737 #ifdef VERBOSE
5738 printf("Total %d nodes\n", queue2->getNumElements() ) ;
5739
5740
5741 // Perform erosion
5742 printf("Start erosion to %d...\n", wid) ;
5743 #endif
5744 gridQueueEle* ele ;
5745 gridPoint* gp ;
5746 int ox, oy, oz ;
5747 int score ;
5748 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ;
5749 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ )
5750 {
5751 scrvol->setDataAt( i, -1 ) ;
5752 }
5753
5754
5755 for ( int curwid = 1 ; curwid <= wid ; curwid ++ )
5756 {
5757 // At the start of each iteration,
5758 // queue2 holds all the nodes for this layer
5759 // queue3 and queue are empty
5760
5761 int numComplex = 0, numSimple = 0 ;
5762 #ifdef VERBOSE
5763 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ;
5764 #endif
5765
5766
5767 // Next,
5768 // Compute score for each node left in queue2
5769 // move them into priority queue
5770 queue2->reset() ;
5771 ele = queue2->getNext() ;
5772 while ( ele != NULL )
5773 {
5774 ox = ele->x ;
5775 oy = ele->y ;
5776 oz = ele->z ;
5777
5778 // Compute score
5779 score = getNumPotComplex2( ox, oy, oz ) ;
5780 scrvol->setDataAt( ox, oy, oz, score ) ;
5781
5782 // Push to queue
5783 gp = new gridPoint ;
5784 gp->x = ox ;
5785 gp->y = oy ;
5786 gp->z = oz ;
5787 // queue->add( gp, -score ) ;
5788 queue->add( gp, score ) ;
5789
5790 ele = queue2->remove() ;
5791 }
5792
5793 // Rename queue3 to be queue2,
5794 // Clear queue3
5795 // From now on, queue2 holds nodes of next level
5796 delete queue2 ;
5797 queue2 = queue3 ;
5798 queue3 = new GridQueue2( ) ;
5799
5800 // Next, start priority queue iteration
5801 while ( ! queue->isEmpty() )
5802 {
5803 // Retrieve the node with the highest score
5804 queue->remove( gp, score ) ;
5805 ox = gp->x ;
5806 oy = gp->y ;
5807 oz = gp->z ;
5808 delete gp ;
5809 // score = -score ;
5810
5811 // Ignore the node
5812 // if it has been processed before
5813 // or it has an updated score
5814 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score )
5815 {
5816 continue ;
5817 }
5818
5819 /* Added for debugging */
5820 // Check simple
5821 if ( ! isSimple( ox, oy, oz ) )
5822 {
5823 // Complex, set to next layer
5824 setDataAt( ox, oy, oz, curwid + 1 ) ;
5825 queue2->prepend( ox, oy, oz ) ;
5826 numComplex ++ ;
5827 }
5828 else
5829 {
5830 setDataAt( ox, oy, oz, -1 ) ;
5831 numSimple ++ ;
5832 }
5833 /* Adding ends */
5834
5835 // Move its neighboring unvisited node to queue2
5836 for ( int m = 0 ; m < 6 ; m ++ )
5837 {
5838 int nx = ox + neighbor6[m][0] ;
5839 int ny = oy + neighbor6[m][1] ;
5840 int nz = oz + neighbor6[m][2] ;
5841 if ( getDataAt( nx, ny, nz ) == 0 )
5842 {
5843 setDataAt( nx, ny, nz, curwid + 1 ) ;
5844 queue2->prepend( nx, ny, nz ) ;
5845 }
5846 }
5847
5848 // Update scores for nodes in its 5x5 neighborhood
5849 // insert them back into priority queue
5850
5851 for ( i = -2 ; i < 3 ;i ++ )
5852 for ( j = -2 ; j < 3 ; j ++ )
5853 for ( k = -2 ; k < 3 ; k ++ )
5854 {
5855 int nx = ox + i ;
5856 int ny = oy + j ;
5857 int nz = oz + k ;
5858
5859 if ( getDataAt( nx, ny, nz ) == curwid )
5860 {
5861 // Compute score
5862 score = getNumPotComplex2( nx, ny, nz ) ;
5863
5864 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) )
5865 {
5866 // printf("Update\n") ;
5867 scrvol->setDataAt( nx, ny, nz, score ) ;
5868 // Push to queue
5869 gp = new gridPoint ;
5870 gp->x = nx ;
5871 gp->y = ny ;
5872 gp->z = nz ;
5873 // queue->add( gp, -score ) ;
5874 queue->add( gp, score ) ;
5875 }
5876 }
5877 }
5878
5879
5880 }
5881
5882 #ifdef VERBOSE
5883 printf("%d complex, %d simple\n", numComplex, numSimple) ;
5884 #endif
5885
5886 if ( numSimple == 0 )
5887 {
5888 delete queue2 ;
5889 break ;
5890 }
5891 }
5892
5893 // Finally, clean up
5894 #ifdef VERBOSE
5895 printf("Thresholding the volume to 0/1...\n") ;
5896 #endif
5897 threshold( 0, 0, 1 ) ;
5898 delete scrvol;
5899 delete queue;
5900 delete queue3;
5901 }
5902
5903
5904 // Apply helix erosion
5905 void Volume::erodeHelix( )
5906 {
5907 erodeHelix( 3 ) ;
5908 }
5909
5910
5911 void Volume::erodeHelix( int disthr )
5912 {
5913 int i, j, k ;
5914 // First, threshold the volume
5915 //printf("Thresholding the volume to -1/0...\n") ;
5916 threshold( 0.1f, -1, 0 ) ;
5917
5918 /* Debug: remove faces */
5919 //Volume* facevol = markFaceEdge() ;
5920 /* End debugging */
5921
5922 // Next, initialize the linked queue
5923 //printf("Initializing queue...\n") ;
5924 Volume * fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ;
5925 GridQueue2* queue2 = new GridQueue2( ) ;
5926 GridQueue2* queue3 = new GridQueue2( ) ;
5927 GridQueue2** queues = new GridQueue2* [ 10000 ] ;
5928
5929 for ( i = 0 ; i < getSizeX() ; i ++ )
5930 for ( j = 0 ; j < getSizeY() ; j ++ )
5931 for ( k = 0 ; k < getSizeZ() ; k ++ )
5932 {
5933 if ( getDataAt( i, j, k ) >= 0 )
5934 {
5935 if ( ! hasCompleteHelix( i, j, k ) )
5936 // if ( ! hasCompleteHelix( i, j, k, facevol ) )
5937 {
5938 queue2->prepend( i, j, k ) ;
5939 fvol->setDataAt( i, j, k, -1 ) ;
5940 }
5941 }
5942 }
5943 //printf("Total %d nodes\n", queue2->getNumElements() ) ;
5944
5945 // Start erosion
5946 gridQueueEle* ele ;
5947 int dis = -1 ;
5948 while ( queue2->getNumElements() > 0 )
5949 {
5950 // First, set distance
5951 dis -- ;
5952 queues[ - dis ] = new GridQueue2( ) ;
5953 //printf("Distance transform to %d...", dis) ;
5954 queue2->reset() ;
5955 while( ( ele = queue2->getNext() ) != NULL )
5956 {
5957 setDataAt( ele->x, ele->y, ele->z, dis ) ;
5958 queues[ -dis ]->prepend( ele->x, ele->y, ele->z ) ;
5959 }
5960 //printf("%d nodes\n", queues[-dis]->getNumElements()) ;
5961
5962 // Next, find next layer
5963 queue2->reset() ;
5964 ele = queue2->getNext() ;
5965 while ( ele != NULL )
5966 {
5967 for ( int m = 0 ; m < 6 ; m ++ )
5968 {
5969 int nx = ele->x + neighbor6[m][0] ;
5970 int ny = ele->y + neighbor6[m][1] ;
5971 int nz = ele->z + neighbor6[m][2] ;
5972 if ( getDataAt( nx, ny, nz ) == 0 )
5973 {
5974 fvol->setDataAt( nx, ny, nz, dis ) ;
5975
5976 if ( ! hasCompleteHelix( nx, ny, nz ) )
5977 // if ( ! hasCompleteHelix( nx, ny, nz, facevol ) )
5978 {
5979 setDataAt( nx, ny, nz, 1 ) ;
5980 queue3->prepend( nx, ny, nz ) ;
5981 }
5982 }
5983 }
5984
5985 ele = queue2->remove() ;
5986 }
5987
5988 // Next, swap queues
5989 GridQueue2* temp = queue2 ;
5990 queue2 = queue3 ;
5991 queue3 = temp ;
5992 }
5993
5994 // Deal with closed rings
5995 dis -- ;
5996 queues[ - dis ] = new GridQueue2( ) ;
5997 #ifdef VERBOSE
5998 printf("Detecting closed rings %d...", dis) ;
5999 #endif
6000 int ftot = 0 ;
6001 for ( i = 0 ; i < getSizeX() ; i ++ )
6002 for ( j = 0 ; j < getSizeY() ; j ++ )
6003 for ( k = 0 ; k < getSizeZ() ; k ++ )
6004 {
6005 if ( getDataAt( i, j, k ) == 0 )
6006 {
6007 setDataAt( i, j, k, dis ) ;
6008 queues[ -dis ]->prepend( i, j, k ) ;
6009 /*
6010 int fval = (int) fvol->getDataAt( i, j, k ) ;
6011 if ( fval == 0)
6012 {
6013 // queues[ -dis ]->prepend( i, j, k ) ;
6014 }
6015 else
6016 {
6017 setDataAt( i, j, k, fval - 1 ) ;
6018 // queues[ -fval + 1 ]->prepend( i, j, k ) ;
6019 }
6020 */
6021 ftot ++ ;
6022