00001 // Copyright (C) 2005-2008 Washington University in St Louis, Baylor College of Medicine. All rights reserved 00002 // Author: Tao Ju (taoju@cse.wustl.edu), Refactored by Sasakthi Abeysinghe (sasakthi.abeysinghe@wustl.edu) 00003 // Description: Volumetric data definition 00004 00005 #include "volume.h" 00006 00007 using namespace wustl_mm::SkeletonMaker; 00008 00009 Volume::Volume(EMData* em) //eman2 00010 { 00011 this->volData = new VolumeData(em); 00012 } 00013 Volume::Volume(int x, int y, int z) { 00014 volData = new VolumeData(x, y, z); 00015 } 00016 00017 Volume::Volume(int x, int y, int z, float val) { 00018 volData = new VolumeData(x, y, z, val); 00019 } 00020 00021 Volume::Volume(int x, int y, int z, int offx, int offy, int offz, Volume * vol) { 00022 volData = new VolumeData(x, y, z, offx, offy, offz, vol->getVolumeData()); 00023 } 00024 00025 00026 Volume::~Volume( ) 00027 { 00028 delete volData; 00029 } 00030 00031 EMData* Volume::get_emdata() //eman2 00032 { 00033 return this->getVolumeData()->get_emdata(); 00034 } 00035 00036 int Volume::getSizeX() { 00037 return volData->GetSizeX(); 00038 } 00039 00040 int Volume::getSizeY() { 00041 return volData->GetSizeY(); 00042 } 00043 00044 int Volume::getSizeZ() { 00045 return volData->GetSizeZ(); 00046 } 00047 00048 int Volume::getIndex(int x, int y, int z) { 00049 return volData->GetIndex(x, y, z); 00050 } 00051 00052 void Volume::setDataAt( int x, int y, int z, double d ) { 00053 volData->SetDataAt(x, y, z, (float)d); 00054 } 00055 00056 void Volume::setDataAt( int index, double d ) { 00057 volData->SetDataAt(index, (float)d); 00058 } 00059 00060 double Volume::getDataAt( int x, int y, int z ) 00061 { 00062 return volData->GetDataAt(x, y, z); 00063 } 00064 00065 double Volume::getDataAt( int index ) { 00066 return volData->GetDataAt(index); 00067 } 00068 00069 VolumeData * Volume::getVolumeData() { 00070 return volData; 00071 } 00072 00073 void Volume::setSpacing(float spx, float spy, float spz ) { 00074 volData->SetSpacing(spx, spy, spz); 00075 } 00076 00077 void Volume::setOrigin(float orgX, float orgY, float orgZ) { 00078 volData->SetOrigin(orgX, orgY, orgZ); 00079 } 00080 00081 float Volume::getSpacingX() { 00082 return volData->GetSpacingX(); 00083 } 00084 00085 float Volume::getSpacingY() { 00086 return volData->GetSpacingY(); 00087 } 00088 00089 float Volume::getSpacingZ() { 00090 return volData->GetSpacingZ(); 00091 } 00092 00093 float Volume::getOriginX() { 00094 return volData->GetOriginX(); 00095 } 00096 00097 float Volume::getOriginY() { 00098 return volData->GetOriginY(); 00099 } 00100 00101 float Volume::getOriginZ() { 00102 return volData->GetOriginZ(); 00103 } 00104 00105 //Volume * Volume::getPseudoDensity( ) { 00108 //int i, j, k ; 00109 //Volume * res = new Volume(getSizeX(), getSizeY(), getSizeZ(), 0, 0, 0, this); 00110 //int size = getSizeX() * getSizeY() * getSizeZ() ; 00111 //srand(123) ; 00112 00113 //for ( i = 0 ; i < getSizeX() ; i ++ ) 00114 //for ( j = 0 ; j < getSizeY() ; j ++ ) 00115 //for ( k = 0 ; k < getSizeZ() ; k ++ ) { 00116 //if ( res->getDataAt( i, j, k ) > 0 ) { 00117 //int ct = 0 ; 00118 //for ( int m = 0 ; m < 6 ; m ++ ) { 00119 //if ( res->getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) > 0 ) { 00120 //ct ++ ; 00121 //} 00122 //} 00123 //res->setDataAt( i,j,k, (k/(float)getSizeZ())*(k/(float)getSizeZ()) ) ; 00124 //if ( ct > 2 ) { 00126 //} else { 00128 //} 00129 //} 00130 //} 00131 00133 //for ( i = 0 ; i < 20 ; i ++ ) 00134 //{ 00135 //printf("Smoothing round %d\n", i) ; 00136 //res->smooth( 0.5f ) ; 00137 //} 00138 //*/ 00139 00140 //Volume * tvol = new Volume( getSizeX(), getSizeY(), getSizeZ(), 0, 0, 0, res); 00141 //float d, ad, ct, temp; 00142 //for ( int it = 0 ; it < 3 ; it ++ ) 00143 //for ( i = 0 ; i < getSizeX() ; i ++ ) 00144 //for ( j = 0 ; j < getSizeY() ; j ++ ) 00145 //for ( k = 0 ; k < getSizeZ() ; k ++ ) { 00146 //if ( (d = (float)tvol->getDataAt( i, j, k )) > 0 ) { 00147 //ad = 0 ; ct = 0 ; 00148 //for ( int m = 0 ; m < 6 ; m ++ ) { 00149 //if ( (temp = (float)tvol->getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] )) > 0 ) { 00150 //ad += temp; 00151 //ct ++ ; 00152 //} 00153 //} 00154 //if ( ct > 0 ) { 00155 //res->setDataAt( i, j, k, ( d + ad/ct ) / 2 ) ; 00156 //} 00157 //} 00158 //} 00159 00160 //delete tvol; 00161 //tvol = new Volume( getSizeX(), getSizeY(), getSizeZ(), 0, 0, 0, res ) ; 00162 //for ( i = 0 ; i < 40 ; i ++ ) 00163 //{ 00164 //printf("Smoothing round %d\n", i) ; 00165 //res->smooth( 0.5f ) ; 00166 //continue ; 00167 00176 00177 //} 00178 00179 00180 //return res ; 00181 //} 00182 00183 00184 //Volume * Volume::getDistanceField(int rad, float randf) { 00188 00190 //int i, j, k ; 00191 //Volume * res = new Volume(getSizeX(), getSizeY(), getSizeZ(), 0, 0, 0, this); 00192 //srand( 123 ) ; 00193 00194 //for ( i = 0 ; i < getSizeX() ; i ++ ) 00195 //for ( j = 0 ; j < getSizeY() ; j ++ ) 00196 //for ( k = 0 ; k < getSizeZ() ; k ++ ) { 00197 //if ( getDataAt(i, j, k) > 0 ) { 00198 //float mag = 1 + randf * (float) rand() / (float) RAND_MAX ; 00199 //int lx = max(0,i-rad) ; 00200 //int ly = max(0,j-rad) ; 00201 //int lz = max(0,k-rad) ; 00202 //int hx = min(getSizeX()-1,i+rad) ; 00203 //int hy = min(getSizeY()-1,j+rad) ; 00204 //int hz = min(getSizeZ()-1,k+rad) ; 00205 //int x,y,z; 00206 //for ( x = lx ; x <= hx ; x ++ ) 00207 //for ( y = ly ; y <= hy ; y ++ ) 00208 //for ( z = lz ; z <= hz ; z ++ ) { 00209 //float val = 1 - (float) sqrt((double)((x-i)*(x-i) + (y-j)*(y-j) + (z-k)*(z-k))) / (float) rad ; 00210 //val *= mag ; 00211 //if ( res->getDataAt( x, y, z ) < val ) { 00212 //res->setDataAt( x, y, z, val ) ; 00213 //} 00214 //} 00215 //} 00216 //} 00217 00219 //for ( i = 0 ; i < 2 ; i ++ ) 00220 //{ 00221 //printf("Smoothing round %d\n", i) ; 00222 //res->smooth( 0.5f ) ; 00223 //} 00224 00225 00226 //return res ; 00227 //} 00228 00229 00230 //int Volume::getNonZeroVoxelCount() { 00231 //int count = 0; 00232 //for(int x = 0; x < getSizeX(); x++){ 00233 //for(int y = 0; y < getSizeY(); y++){ 00234 //for(int z = 0; z < getSizeZ(); z++){ 00235 //if(this->getDataAt(x, y, z) > 0.0) { 00236 //count++; 00237 //} 00238 //} 00239 //} 00240 //} 00241 //return count; 00242 //} 00243 //void Volume::print() { 00244 //for(int x = 0; x < getSizeX(); x++) { 00245 //printf("{ "); 00246 //for(int y = 0; y < getSizeY(); y++) { 00247 //printf("{ "); 00248 //for(int z = 0; z < getSizeZ(); z++) { 00249 //printf("%f, ", getDataAt(x, y, z)); 00250 //} 00251 //printf("} "); 00252 //} 00253 //printf("} "); 00254 //} 00255 //printf("\n"); 00256 //} 00257 00258 00259 //void Volume::subtract( Volume* vol ) { 00260 //int i, j, k ; 00261 //for ( i = 0 ; i < getSizeX() ; i ++ ) 00262 //for ( j = 0 ; j < getSizeY() ; j ++ ) 00263 //for ( k = 0 ; k < getSizeZ() ; k ++ ) { 00264 //if ( getDataAt( i, j, k ) > 0 ) { 00265 //if ( vol->getDataAt(i,j,k) > 0 ) { 00266 //setDataAt( i, j, k, 0 ) ; 00267 //} 00268 //} 00269 //} 00270 00271 //} 00272 00273 void Volume::pad (int padBy, double padValue) { 00274 volData->Pad(padBy, padValue); 00275 } 00276 00277 //void Volume::applyMask(Volume * maskVol, double maskValue, bool keepMaskValue) { 00278 //for(int x = 0; x < maskVol->getSizeX(); x++) { 00279 //for(int y = 0; y < maskVol->getSizeY(); y++) { 00280 //for(int z = 0; z < maskVol->getSizeZ(); z++) { 00281 //if(((maskVol->getDataAt(x, y, z) == maskValue) && !keepMaskValue) || 00282 //((maskVol->getDataAt(x, y, z) != maskValue) && keepMaskValue)) { 00283 //setDataAt(x, y, z, 0); 00284 //} 00285 //} 00286 //} 00287 //} 00288 //} 00289 00290 //double Volume::getMin() { 00291 //int size = volData->GetMaxIndex(); 00292 //double rvalue = volData->GetDataAt(0); 00293 //for (int i=1; i < size; i++) { 00294 //float val = volData->GetDataAt(i); 00295 //if ( rvalue > val) { 00296 //rvalue = val; 00297 //} 00298 //} 00299 //return rvalue; 00300 //} 00301 00302 //double Volume::getMax() { 00303 //int size = volData->GetMaxIndex(); 00304 //double rvalue = volData->GetDataAt(0); 00305 //for (int i=1; i < size; i++) { 00306 //float val = volData->GetDataAt(i); 00307 //if ( rvalue < val) { 00308 //rvalue = val; 00309 //} 00310 //} 00311 //return rvalue ; 00312 //} 00313 00314 //double Volume::getMaxValuePosition(int& maxX, int& maxY, int& maxZ) { 00315 //double maxVal = getDataAt(0,0,0); 00316 //maxX = 0; maxY = 0; maxZ = 0; 00317 //double data; 00318 00319 //for(int x = 0; x < getSizeX(); x++) { 00320 //for(int y = 0; y < getSizeY(); y++) { 00321 //for(int z = 0; z < getSizeZ(); z++) { 00322 //data = getDataAt(x, y, z); 00323 //if(data > maxVal) { 00324 //maxVal = data; 00325 //maxX = x; maxY = y; maxZ = z; 00326 //} 00327 //} 00328 //} 00329 //} 00330 //return maxVal; 00331 //} 00332 00333 //double Volume::getLocalMax(int x, int y, int z, int radius) { 00334 //double mx = getDataAt(x, y, z); 00335 //for(int xx = x - radius; xx <= x + radius; xx++) { 00336 //for(int yy = y - radius; yy <= y + radius; yy++) { 00337 //for(int zz = z - radius; zz <= z + radius; zz++) { 00338 //mx = max(mx, getDataAt(xx, yy, zz)); 00339 //} 00340 //} 00341 //} 00342 //return mx; 00343 //} 00344 00345 //double Volume::getLocalMin(int x, int y, int z, int radius) { 00346 //double mn = getDataAt(x, y, z); 00347 //for(int xx = x - radius; xx <= x + radius; xx++) { 00348 //for(int yy = y - radius; yy <= y + radius; yy++) { 00349 //for(int zz = z - radius; zz <= z + radius; zz++) { 00350 //mn = min(mn, getDataAt(xx, yy, zz)); 00351 //} 00352 //} 00353 //} 00354 //return mn; 00355 //} 00356 00357 //void Volume::fill( double val ) 00358 //{ 00359 //for(int x = 0; x < getSizeX(); x++) { 00360 //for(int y = 0; y < getSizeY(); y++) { 00361 //for(int z = 0; z < getSizeZ(); z++) { 00362 //setDataAt(x, y, z, val); 00363 //} 00364 //} 00365 //} 00366 //} 00367 00368 //int Volume::isBertrandBorder( int ox, int oy, int oz, int dir ) { 00369 //int nx = ox + neighbor6[dir][0] ; 00370 //int ny = oy + neighbor6[dir][1] ; 00371 //int nz = oz + neighbor6[dir][2] ; 00372 //if ( getDataAt( nx, ny, nz) < 0 ) { 00373 //return 1 ; 00374 //} 00375 //return 0 ; 00376 //} 00377 00378 //int Volume::isBertrandEndPoint( int ox, int oy, int oz ) { 00379 //double vox[3][3][3] ; 00380 00381 //int i, j, k ; 00382 //for ( i = -1 ; i < 2 ; i ++ ) 00383 //for ( j = -1 ; j < 2 ; j ++ ) 00384 //for ( k = -1 ; k < 2 ; k ++ ) { 00385 //double tval = getDataAt( ox + i, oy + j, oz + k ) ; 00386 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ; 00387 //} 00388 00390 //int xx6 = 0 ; 00391 //for ( i = 0 ; i < 6 ; i ++ ) { 00392 //if ( vox[ neighbor6[i][0] + 1 ][ neighbor6[i][1] + 1 ][ neighbor6[i][2] + 1 ] >= 0 ) { 00393 //xx6 ++ ; 00394 //} 00395 //} 00396 00398 //int a26 = 0, na26 = 0 ; 00399 //for ( i = 0 ; i < 3 ; i += 2 ) 00400 //for ( j = 0 ; j < 3 ; j += 2 ) 00401 //for ( k = 0 ; k < 3 ; k += 2 ) { 00402 //if ( vox[1][j][k] < 0 ) { 00403 //continue ; 00404 //} 00405 //if ( vox[i][1][k] < 0 ) { 00406 //continue ; 00407 //} 00408 //if ( vox[i][j][1] < 0 ) { 00409 //continue ; 00410 //} 00411 //if ( vox[i][1][1] < 0 ) { 00412 //continue ; 00413 //} 00414 //if ( vox[1][j][1] < 0 ) { 00415 //continue ; 00416 //} 00417 //if ( vox[1][1][k] < 0 ) { 00418 //continue ; 00419 //} 00420 //if ( vox[i][j][k] >= 0 ) { 00421 //a26 ++ ; 00422 //} else { 00423 //na26 ++ ; 00424 //} 00425 //} 00426 00427 //if (( na26 == 0 ) && ( a26 != 0 ) && ( xx6 <= a26 + 2 )) { 00428 //return 0 ; 00429 //} 00430 //if ( isFeatureFace(ox,oy,oz) ) { 00432 //} 00433 //return 1 ; 00434 //} 00435 00436 //int Volume::isHelix( int ox, int oy, int oz ) { 00437 //int cn = 12 ; 00438 //int nx, ny, nz ; 00439 //int i, j, k ; 00440 00441 //double vox[3][3][3] ; 00442 //for ( i = -1 ; i < 2 ; i ++ ) 00443 //for ( j = -1 ; j < 2; j ++ ) 00444 //for ( k = -1 ; k < 2 ; k ++ ) { 00445 //vox[i+1][j+1][k+1] = getDataAt( ox + i, oy + j, oz + k ) ; 00446 //} 00447 00448 //for ( i = 0 ; i < 12 ; i ++ ) { 00449 //for ( j = 0 ; j < 4 ; j ++ ) { 00450 //nx = sheetNeighbor[i][j][0] + 1; 00451 //ny = sheetNeighbor[i][j][1] + 1; 00452 //nz = sheetNeighbor[i][j][2] + 1; 00453 00454 //if ( vox[nx][ny][nz] <= 0 ) { 00455 //cn -- ; 00456 //break ; 00457 //} 00458 //} 00459 //} 00460 00461 //if ( cn >= 1 ) { 00462 //return 0 ; 00463 //} else { 00464 //return 1 ; 00465 //} 00466 //} 00467 00468 //int Volume::isSheet( int ox, int oy, int oz ) { 00469 //int cn = 12 ; 00470 //int nx, ny, nz ; 00471 00472 //for ( int i = 0 ; i < 12 ; i ++ ) { 00473 //for ( int j = 0 ; j < 4 ; j ++ ) { 00474 //nx = ox + sheetNeighbor[i][j][0] ; 00475 //ny = oy + sheetNeighbor[i][j][1] ; 00476 //nz = oz + sheetNeighbor[i][j][2] ; 00477 00478 //if ( getDataAt( nx, ny, nz ) <= 0 ) { 00479 //cn -- ; 00480 //break ; 00481 //} 00482 //} 00483 //} 00484 //return ( cn >= 3 ) ; 00485 //} 00486 00487 //Volume * Volume::getSheets( int minSize ) { 00488 //int i, j, k ; 00489 00491 //printf("Initialize volume at %d %d %d\n", getSizeX(), getSizeY(), getSizeZ() ) ; 00492 //Volume* svol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ; 00493 00495 //int sheets[MAX_SHEETS] ; 00496 //for ( i = 0 ; i < MAX_SHEETS ; i ++ ) { 00497 //sheets[ i ] = 0 ; 00498 //} 00499 //int totSheets = 1 ; 00500 00502 //printf("Start clustering...\n" ) ; 00503 //int ox, oy, oz ; 00504 //for ( i = 0 ; i < getSizeX() ; i ++ ) 00505 //for ( j = 0 ; j < getSizeY() ; j ++ ) 00506 //for ( k = 0 ; k < getSizeZ() ; k ++ ) { 00507 //if ( getDataAt(i,j,k) <= 0 || svol->getDataAt(i,j,k) != 0 ) { 00509 //continue ; 00510 //} 00511 //if ( ! isSheet( i, j, k ) ) { 00513 //continue ; 00514 //} 00515 00517 //int numNodes = 1 ; 00518 //svol->setDataAt( i, j, k, totSheets ) ; 00519 //GridQueue* queue = new GridQueue() ; 00520 //queue->pushQueue( i, j, k ) ; 00521 //while ( queue->popQueue(ox, oy, oz) ) { 00523 //if ( isSheet( ox, oy, oz ) ) { 00524 //for ( int m = 0 ; m < 6 ; m ++ ) { 00525 //int nx = ox + neighbor6[m][0] ; 00526 //int ny = oy + neighbor6[m][1] ; 00527 //int nz = oz + neighbor6[m][2] ; 00528 00529 //if ( getDataAt(nx,ny,nz) > 0 && svol->getDataAt(nx,ny,nz) == 0 ) { 00530 //svol->setDataAt(nx,ny,nz,totSheets); 00531 //queue->pushQueue(nx,ny,nz) ; 00532 //numNodes ++ ; 00533 //} 00534 //} 00535 //} 00536 //} 00537 00538 //delete queue ; 00539 //if ( numNodes > 0 ) { 00541 //sheets[ totSheets ] = numNodes ; 00542 //totSheets ++ ; 00543 //} 00544 //} 00545 00547 //printf("Removing small clusters.\n") ; 00548 //for ( i = 0 ; i < getSizeX() ; i ++ ) 00549 //for ( j = 0 ; j < getSizeY() ; j ++ ) 00550 //for ( k = 0 ; k < getSizeZ() ; k ++ ) { 00551 //int cnt = (int) svol->getDataAt(i,j,k) ; 00552 //if ( cnt > 0 && sheets[ cnt ] < minSize ) { 00553 //svol->setDataAt(i,j,k,-1) ; 00554 //} 00555 //} 00556 00558 //#ifdef VERBOSE 00559 //printf("Thresholding the volume to 0/1...\n") ; 00560 //#endif 00561 //svol->threshold( 0.1, 0, 1 ) ; 00562 00563 //return svol ; 00564 //} 00565 00566 //Volume * Volume::getHelices( int minSize ) { 00567 //printf("Segmenting helices from eroded volume.\n") ; 00568 //int i, j, k ; 00569 00571 //printf("Initialize volume at %d %d %d\n", getSizeX(), getSizeY(), getSizeZ() ) ; 00572 //Volume* svol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ; 00573 00575 //int helices[MAX_SHEETS] ; 00576 //for ( i = 0 ; i < MAX_SHEETS ; i ++ ) { 00577 //helices[ i ] = 0 ; 00578 //} 00579 //int totHelices = 1 ; 00580 00582 //printf("Start clustering...\n" ) ; 00583 //int ox, oy, oz ; 00584 //for ( i = 0 ; i < getSizeX() ; i ++ ) 00585 //for ( j = 0 ; j < getSizeY() ; j ++ ) 00586 //for ( k = 0 ; k < getSizeZ() ; k ++ ) { 00587 //if ( getDataAt(i,j,k) <= 0 || svol->getDataAt(i,j,k) != 0 ) { 00589 //continue ; 00590 //} 00591 //if ( ! isHelix( i, j, k ) ) { 00593 //continue ; 00594 //} 00595 00597 //int numNodes = 1 ; 00598 //svol->setDataAt( i, j, k, totHelices ) ; 00599 //GridQueue* queue = new GridQueue() ; 00600 //queue->pushQueue( i, j, k ) ; 00601 //while ( queue->popQueue(ox, oy, oz) ) 00602 //{ 00604 //if ( isHelix( ox, oy, oz ) ) { 00605 //for ( int m = 0 ; m < 6 ; m ++ ) { 00606 //int nx = ox + neighbor6[m][0] ; 00607 //int ny = oy + neighbor6[m][1] ; 00608 //int nz = oz + neighbor6[m][2] ; 00609 00610 //if ( getDataAt(nx,ny,nz) > 0 && svol->getDataAt(nx,ny,nz) == 0 ) { 00611 //svol->setDataAt(nx,ny,nz,totHelices); 00612 //queue->pushQueue(nx,ny,nz) ; 00613 //numNodes ++ ; 00614 //} 00615 //} 00616 //} 00617 //} 00618 00619 //delete queue ; 00620 //if ( numNodes > 0 ) { 00622 //helices[ totHelices ] = numNodes ; 00623 //totHelices ++ ; 00624 //} 00625 //} 00626 00628 //printf("Removing small clusters.\n") ; 00629 //for ( i = 0 ; i < getSizeX() ; i ++ ) 00630 //for ( j = 0 ; j < getSizeY() ; j ++ ) 00631 //for ( k = 0 ; k < getSizeZ() ; k ++ ) { 00632 //int cnt = (int) svol->getDataAt(i,j,k) ; 00633 //if ( cnt > 0 && helices[ cnt ] < minSize ) { 00634 //svol->setDataAt(i,j,k,-1) ; 00635 //} 00636 //} 00637 00639 //#ifdef VERBOSE 00640 //printf("Thresholding the volume to 0/1...\n") ; 00641 //#endif 00642 //svol->threshold( 0.1, 0, 1 ) ; 00643 00644 //return svol ; 00645 00646 //} 00647 00648 //int Volume::isEndPoint( int ox, int oy, int oz ) { 00649 //if ( getDataAt( ox - 1, oy, oz ) < 0 && getDataAt( ox + 1, oy, oz ) < 0 ) { 00650 //return 1 ; 00651 //} 00652 //if ( getDataAt( ox, oy - 1, oz ) < 0 && getDataAt( ox, oy + 1, oz ) < 0 ) { 00653 //return 1 ; 00654 //} 00655 //if ( getDataAt( ox, oy, oz - 1 ) < 0 && getDataAt( ox, oy, oz + 1 ) < 0 ) { 00656 //return 1 ; 00657 //} 00658 //return 0 ; 00659 //} 00660 00661 int Volume::getNumNeighbor6( int ox, int oy, int oz ) { 00662 int rvalue = 0 ; 00663 for ( int i = 0 ; i < 6 ; i ++ ) { 00664 int nx = ox + neighbor6[i][0] ; 00665 int ny = oy + neighbor6[i][1] ; 00666 int nz = oz + neighbor6[i][2] ; 00667 if ( getDataAt( nx, ny, nz ) >= 0 ) { 00668 rvalue ++ ; 00669 } 00670 } 00671 00672 return rvalue ; 00673 } 00674 //int Volume::testIsSheetEnd( int ox, int oy, int oz ) { 00676 //int i, j ; 00677 //int nx, ny, nz ; 00678 00679 //int edge[6] = { 0,0,0,0,0,0 } ; 00680 //int faceflag[ 12 ] ; 00681 //int hasNoiseFace = 0 ; 00682 //int tot = 0 ; 00683 00684 //for ( i = 0 ; i < 12 ; i ++ ) { 00685 //faceflag[ i ] = 1 ; 00686 //int hasNoise = 0 ; 00687 00688 //for ( j = 0 ; j < 4 ; j ++ ) { 00689 //nx = ox + sheetNeighbor[i][j][0] ; 00690 //ny = oy + sheetNeighbor[i][j][1] ; 00691 //nz = oz + sheetNeighbor[i][j][2] ; 00692 00693 //if ( getDataAt( nx, ny, nz ) == 0 ) { 00694 //hasNoise = 1 ; 00695 //} 00696 //if ( getDataAt( nx, ny, nz ) < 0 ) { 00697 //faceflag[ i ] = 0 ; 00698 //break ; 00699 //} 00700 //} 00701 //if ( faceflag[ i ] == 1 && hasNoise ) { 00702 //hasNoiseFace ++ ; 00704 //} 00705 00706 //if ( faceflag[ i ] ) { 00707 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ; 00708 //edge[ e0 ] ++ ; 00709 //edge[ e1 ] ++ ; 00710 //tot ++ ; 00711 //} 00712 //} 00713 00714 //if ( hasNoiseFace == tot ) { 00715 //return 0 ; 00716 //} 00717 00718 //if ( tot == 0 ) { 00719 //return 0 ; 00720 //} 00721 00723 //int numones = 0 ; 00724 //for ( i = 0 ; i < 6 ; i ++ ) { 00725 //if ( edge[ i ] == 1 ) { 00726 //numones ++ ; 00727 //} 00728 //} 00729 //while( numones > 0 ) { 00730 //int e ; 00731 //for ( i = 0 ; i < 6 ; i ++ ) { 00732 //if ( edge[ i ] == 1 ) { 00733 //e = i ; 00734 //break ; 00735 //} 00736 //} 00737 00738 //int f, e2 ; 00739 //for ( j = 0 ; j < 4 ; j ++ ) { 00740 //f = edgeFaces[ e ][ j ] ; 00741 //if ( faceflag[ f ] ) { 00742 //break ; 00743 //} 00744 //} 00745 00746 //if ( faceEdges[ f ][ 0 ] == e ) { 00747 //e2 = faceEdges[ f ][ 1 ] ; 00748 //} else { 00749 //e2 = faceEdges[ f ][ 0 ] ; 00750 //} 00751 00752 //edge[ e ] -- ; 00753 //numones -- ; 00754 //edge[ e2 ] -- ; 00755 //faceflag[ f ] = 0 ; 00756 //tot -- ; 00757 00758 //if ( edge[ e2 ] == 1 ) { 00759 //numones ++ ; 00760 //} else if ( edge[ e2 ] == 0 ) { 00761 //numones -- ; 00762 //} 00763 //} 00764 00765 //if ( tot > 0 ) { 00766 //return 0 ; 00767 //} 00768 //return 1 ; 00769 //} 00770 00771 //int Volume::isNoiseSheetEnd( int ox, int oy, int oz ) { 00772 //int i, j ; 00773 //int nx, ny, nz ; 00774 00775 //int edge[6] = { 0,0,0,0,0,0 } ; 00776 //int faceflag[ 12 ] ; 00777 //int hasNoiseFace = 0 ; 00778 //int tot = 0 ; 00779 00780 //for ( i = 0 ; i < 12 ; i ++ ) { 00781 //faceflag[ i ] = 1 ; 00782 //int hasNoise = 0 ; 00783 00784 //for ( j = 0 ; j < 4 ; j ++ ) { 00785 //nx = ox + sheetNeighbor[i][j][0] ; 00786 //ny = oy + sheetNeighbor[i][j][1] ; 00787 //nz = oz + sheetNeighbor[i][j][2] ; 00788 00789 //if ( getDataAt( nx, ny, nz ) == 0 ) { 00790 //hasNoise = 1 ; 00791 //} 00792 //if ( getDataAt( nx, ny, nz ) < 0 ) { 00793 //faceflag[ i ] = 0 ; 00794 //break ; 00795 //} 00796 //} 00797 //if ( faceflag[ i ] == 1 && hasNoise ) { 00798 //hasNoiseFace ++ ; 00800 //} 00801 00802 //if ( faceflag[ i ] ) { 00803 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ; 00804 //edge[ e0 ] ++ ; 00805 //edge[ e1 ] ++ ; 00806 //tot ++ ; 00807 //} 00808 //} 00809 00810 //if ( hasNoiseFace < tot ) { 00811 //return 0 ; 00812 //} 00813 00814 //if ( tot == 0 ) { 00815 //return 0 ; 00816 //} 00817 00819 //int numones = 0 ; 00820 //for ( i = 0 ; i < 6 ; i ++ ) { 00821 //if ( edge[ i ] == 1 ) 00822 //{ 00823 //numones ++ ; 00824 //} 00825 //} 00826 //while( numones > 0 ) { 00827 //int e ; 00828 //for ( i = 0 ; i < 6 ; i ++ ) { 00829 //if ( edge[ i ] == 1 ) { 00830 //e = i ; 00831 //break ; 00832 //} 00833 //} 00834 00835 //int f, e2 ; 00836 //for ( j = 0 ; j < 4 ; j ++ ) { 00837 //f = edgeFaces[ e ][ j ] ; 00838 //if ( faceflag[ f ] ) { 00839 //break ; 00840 //} 00841 //} 00842 00843 //if ( faceEdges[ f ][ 0 ] == e ) { 00844 //e2 = faceEdges[ f ][ 1 ] ; 00845 //} else { 00846 //e2 = faceEdges[ f ][ 0 ] ; 00847 //} 00848 00849 //edge[ e ] -- ; 00850 //numones -- ; 00851 //edge[ e2 ] -- ; 00852 //faceflag[ f ] = 0 ; 00853 //tot -- ; 00854 00855 //if ( edge[ e2 ] == 1 ) { 00856 //numones ++ ; 00857 //} else if ( edge[ e2 ] == 0 ) { 00858 //numones -- ; 00859 //} 00860 //} 00861 00862 //if ( tot > 0 ) { 00863 //return 0 ; 00864 //} 00865 //return 1 ; 00866 //} 00867 00868 //int Volume::isInternal( int ox, int oy, int oz ) { 00870 //int i, j, k ; 00871 00872 //for ( i = -1 ; i < 2 ; i ++ ) 00873 //for ( j = -1 ; j < 2 ; j ++ ) 00874 //for ( k = -1 ; k < 2 ; k ++ ) { 00875 //if ( getDataAt( ox + i, oy + j, oz + k ) <= 0 ) { 00876 //return 0 ; 00877 //} 00878 //} 00879 //return 1 ; 00880 //} 00881 00882 //int Volume::isInternal2( int ox, int oy, int oz ) { 00884 //int i, j, k ; 00885 00886 //for ( i = -1 ; i < 2 ; i ++ ) 00887 //for ( j = -1 ; j < 2 ; j ++ ) 00888 //for ( k = -1 ; k < 2 ; k ++ ) { 00889 //if ( getDataAt( ox + i, oy + j, oz + k ) < 0 ) { 00890 //return 0 ; 00891 //} 00892 //} 00893 00894 //return 1 ; 00895 //} 00896 00897 //int Volume::hasIsolatedFace( int ox, int oy, int oz ) { 00898 //int i, j, k ; 00899 //int nx, ny, nz ; 00900 00901 //double vox[3][3][3] ; 00902 //for ( i = -1 ; i < 2 ; i ++ ) 00903 //for ( j = -1 ; j < 2 ; j ++ ) 00904 //for ( k = -1 ; k < 2 ; k ++ ) { 00905 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ; 00906 //} 00907 00908 //int cells[8] = { 1, 1, 1, 1, 1, 1, 1, 1 } ; 00909 //for ( i = 0 ; i < 8 ; i ++ ) { 00910 //int x = ( ( i >> 2 ) & 1 ) ; 00911 //int y = ( ( i >> 1 ) & 1 ) ; 00912 //int z = ( i & 1 ) ; 00913 //for ( j = 0 ; j < 8 ; j ++ ) { 00914 //nx = x + ( ( j >> 2 ) & 1 ) ; 00915 //ny = y + ( ( j >> 1 ) & 1 ) ; 00916 //nz = z + ( j & 1 ) ; 00917 00918 //if ( vox[nx][ny][nz] < 0 ) { 00919 //cells[i] = 0 ; 00920 //break; 00921 //} 00922 //} 00923 //} 00924 00925 //for ( i = 0 ; i < 12 ; i ++ ) { 00926 //if ( cells[ faceCells[i][0] ] == 1 || cells[ faceCells[i][1] ] == 1 ) { 00927 //continue ; 00928 //} 00929 //int flag = 1 ; 00930 //for ( j = 0 ; j < 4 ; j ++ ) { 00931 //nx = 1 + sheetNeighbor[i][j][0] ; 00932 //ny = 1 + sheetNeighbor[i][j][1] ; 00933 //nz = 1 + sheetNeighbor[i][j][2] ; 00934 00935 //if ( vox[nx][ny][nz] < 0 ) { 00936 //flag = 0 ; 00937 //break; 00938 //} 00939 //} 00940 //if ( flag ) { 00941 //return 1 ; 00942 //} 00943 //} 00944 00945 //return 0 ; 00946 //} 00947 00948 //int Volume::hasIsolatedEdge( int ox, int oy, int oz ) { 00949 //int i, j, k ; 00950 //int nx, ny, nz ; 00951 00952 //double vox[3][3][3] ; 00953 //for ( i = -1 ; i < 2 ; i ++ ) 00954 //for ( j = -1 ; j < 2 ; j ++ ) 00955 //for ( k = -1 ; k < 2 ; k ++ ) { 00956 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ; 00957 //} 00958 00959 //int edge[6] = { 0,0,0,0,0,0 } ; 00960 00961 //for ( i = 0 ; i < 12 ; i ++ ) { 00962 //int flag = 1 ; 00963 //for ( j = 0 ; j < 4 ; j ++ ) { 00964 //nx = 1 + sheetNeighbor[i][j][0] ; 00965 //ny = 1 + sheetNeighbor[i][j][1] ; 00966 //nz = 1 + sheetNeighbor[i][j][2] ; 00967 00968 //if ( vox[nx][ny][nz] < 0 ) { 00969 //flag = 0 ; 00970 //break ; 00971 //} 00972 //} 00973 00974 //if ( flag ) { 00975 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ; 00976 //edge[ e0 ] ++ ; 00977 //edge[ e1 ] ++ ; 00978 //} 00979 //} 00980 00981 //for ( i = 0 ; i < 6 ; i ++ ) { 00982 //if ( edge[i] ) { 00983 //continue ; 00984 //} 00985 00986 //nx = 1 + neighbor6[i][0] ; 00987 //ny = 1 + neighbor6[i][1] ; 00988 //nz = 1 + neighbor6[i][2] ; 00989 00990 //if ( vox[nx][ny][nz] >= 0 ) { 00991 //return 1 ; 00992 //} 00993 00994 //} 00995 00996 //return 0 ; 00997 //} 00998 00999 //int Volume::countFace( int ox, int oy, int oz, int m ) { 01000 //int facenum = 4 ; 01001 //for ( int i = 0 ; i < 4 ; i ++ ) { 01002 //for ( int j = 0 ; j < 4 ; j ++ ) { 01003 //int nx = ox + sheetNeighbor[edgeFaces[m][i]][j][0] ; 01004 //int ny = oy + sheetNeighbor[edgeFaces[m][i]][j][1] ; 01005 //int nz = oz + sheetNeighbor[edgeFaces[m][i]][j][2] ; 01006 01007 //if ( getDataAt( nx, ny, nz ) < 0 ) { 01008 //facenum -- ; 01009 //break ; 01010 //} 01011 //} 01012 //} 01013 01014 //return facenum; 01015 //} 01016 int Volume::hasCell( int ox, int oy, int oz ) { 01017 for ( int i = 0 ; i < 2 ; i ++ ) 01018 for ( int j = 0 ; j < 2 ; j ++ ) 01019 for ( int k = 0 ; k < 2 ; k ++ ) { 01020 if ( getDataAt( ox + i, oy + j, oz + k ) < 0 ) { 01021 return 0 ; 01022 } 01023 } 01024 return 1 ; 01025 } 01026 01027 Volume * Volume::markCellFace( ) { 01028 int i, j, k ; 01029 Volume* fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ; 01030 01031 //return fvol ; 01032 01033 for ( i = 0 ; i < getSizeX() ; i ++ ) 01034 for ( j = 0 ; j < getSizeY() ; j ++ ) 01035 for ( k = 0 ; k < getSizeZ() ; k ++ ) 01036 { 01037 if( getDataAt(i,j,k) >= 0 ) 01038 { 01039 if ( hasCell(i,j,k) ) 01040 { 01041 for ( int m = 0 ; m < 6 ; m ++ ) 01042 { 01043 int nx = i + neighbor6[m][0] ; 01044 int ny = j + neighbor6[m][1] ; 01045 int nz = k + neighbor6[m][2] ; 01046 if ( ! hasCell(nx,ny,nz) ) 01047 { 01048 fvol->setDataAt(i,j,k,(double)(1<<m)) ; 01049 break ; 01050 } 01051 } 01052 } 01053 } 01054 } 01055 01056 01057 return fvol ; 01058 } 01059 01060 //Volume * Volume::markFaceEdge( ) { 01061 //int x,y,z, i,j ; 01062 //Volume* fvol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ; 01063 01065 01066 //for ( x = 0 ; x < getSizeX() ; x ++ ) 01067 //for ( y = 0 ; y < getSizeY() ; y ++ ) 01068 //for ( z = 0 ; z < getSizeZ() ; z ++ ) 01069 //{ 01070 //if( getDataAt(x,y,z) >= 0 ) 01071 //{ 01072 01073 //for ( i = 0 ; i < 3 ; i ++ ) 01074 //{ 01075 //int hasFace = 1 ; 01076 //for ( j = 0 ; j < 4 ; j ++ ) 01077 //{ 01078 //int nx = x + sheetNeighbor[4 * i + 3][j][0] ; 01079 //int ny = y + sheetNeighbor[4 * i + 3][j][1] ; 01080 //int nz = z + sheetNeighbor[4 * i + 3][j][2] ; 01081 01082 //if ( getDataAt( nx, ny, nz ) < 0 ) 01083 //{ 01084 //hasFace = 0 ; 01085 //break ; 01086 //} 01087 //} 01088 01089 //if ( hasFace ) 01090 //{ 01092 //switch( i ) 01093 //{ 01094 //case 0: 01095 //if ( countFace( x, y, z, 0 ) == 1 ) 01096 //{ 01097 //fvol->setDataAt(x, y, z, (double)(1<<0)) ; 01098 //break ; 01099 //} 01100 //if ( countFace( x, y, z, 2 ) == 1 ) 01101 //{ 01102 //fvol->setDataAt(x, y, z, (double)(1<<2)) ; 01103 //break ; 01104 //} 01105 //if ( countFace( x, y + 1, z, 0 ) == 1 ) 01106 //{ 01107 //fvol->setDataAt(x, y + 1, z, (double)(1<<0)) ; 01108 //break ; 01109 //} 01110 //if ( countFace( x, y, z + 1, 2 ) == 1 ) 01111 //{ 01112 //fvol->setDataAt(x, y, z + 1, (double)(1<<2)) ; 01113 //break ; 01114 //} 01115 //printf("Hmmm... a face with no open edges.\n"); 01116 //break ; 01117 //case 1: 01118 //if ( countFace( x, y, z, 0 ) == 1 ) 01119 //{ 01120 //fvol->setDataAt(x, y, z, (double)(1<<0)) ; 01121 //break ; 01122 //} 01123 //if ( countFace( x, y, z, 4 ) == 1 ) 01124 //{ 01125 //fvol->setDataAt(x, y, z, (double)(1<<4)) ; 01126 //break ; 01127 //} 01128 //if ( countFace( x + 1, y, z, 0 ) == 1 ) 01129 //{ 01130 //fvol->setDataAt(x + 1, y, z, (double)(1<<0)) ; 01131 //break ; 01132 //} 01133 //if ( countFace( x, y, z + 1, 4 ) == 1 ) 01134 //{ 01135 //fvol->setDataAt(x, y, z + 1, (double)(1<<4)) ; 01136 //break ; 01137 //} 01138 //printf("Hmmm... a face with no open edges.\n"); 01139 //break ; 01140 //case 2: 01141 //if ( countFace( x, y, z, 2 ) == 1 ) 01142 //{ 01143 //fvol->setDataAt(x, y, z, (double)(1<<2)) ; 01144 //break ; 01145 //} 01146 //if ( countFace( x, y, z, 4 ) == 1 ) 01147 //{ 01148 //fvol->setDataAt(x, y, z, (double)(1<<4)) ; 01149 //break ; 01150 //} 01151 //if ( countFace( x + 1, y, z, 2 ) == 1 ) 01152 //{ 01153 //fvol->setDataAt(x + 1, y, z, (double)(1<<2)) ; 01154 //break ; 01155 //} 01156 //if ( countFace( x, y + 1, z, 4 ) == 1 ) 01157 //{ 01158 //fvol->setDataAt(x, y + 1, z, (double)(1<<4)) ; 01159 //break ; 01160 //} 01161 //printf("Hmmm... a face with no open edges.\n"); 01162 //break ; 01163 //} 01164 //} 01165 //} 01166 01167 //} 01168 //} 01169 01170 01171 //return fvol ; 01172 //} 01173 01174 int Volume::hasCompleteSheet( int ox, int oy, int oz, Volume* fvol ) { 01175 int i, j, k ; 01176 int nx, ny, nz ; 01177 01178 int edge[6] = { 0,0,0,0,0,0 } ; 01179 int faceflag[ 12 ] ; 01180 int tot = 0 ; 01181 int cellflag[ 8 ] ; 01182 01183 int ct = 0 ; 01184 for ( i = -1 ; i < 1 ; i ++ ) 01185 for ( j = -1 ; j < 1 ; j ++ ) 01186 for ( k = -1 ; k < 1 ; k ++ ) 01187 { 01188 if ( hasCell( ox + i, oy + j, oz + k ) ) 01189 { 01190 cellflag[ ct ] = 1 ; 01191 } 01192 else 01193 { 01194 cellflag[ ct ] = 0 ; 01195 } 01196 ct ++ ; 01197 } 01198 01199 for ( i = 0 ; i < 12 ; i ++ ) 01200 { 01201 faceflag[ i ] = 1 ; 01202 for ( j = 0 ; j < 4 ; j ++ ) 01203 { 01204 nx = ox + sheetNeighbor[i][j][0] ; 01205 ny = oy + sheetNeighbor[i][j][1] ; 01206 nz = oz + sheetNeighbor[i][j][2] ; 01207 01208 if ( getDataAt( nx, ny, nz ) < 0 ) 01209 { 01210 faceflag[ i ] = 0 ; 01211 break ; 01212 } 01213 } 01214 01215 if ( faceflag[ i ] ) 01216 { 01217 if ( cellflag[ faceCells[i][0] ] ^ cellflag[ faceCells[i][1] ] ) 01218 { 01219 int v1 = (int)( fvol->getDataAt( 01220 ox - 1 + (( faceCells[i][0] >> 2 ) & 1 ), 01221 oy - 1 + (( faceCells[i][0] >> 1 ) & 1 ), 01222 oz - 1 + (( faceCells[i][0] ) & 1)) ) ; 01223 int v2 = (int)( fvol->getDataAt( 01224 ox - 1 + (( faceCells[i][1] >> 2 ) & 1 ), 01225 oy - 1 + (( faceCells[i][1] >> 1 ) & 1 ), 01226 oz - 1 + (( faceCells[i][1] ) & 1)) ) ; 01227 if ( ((v1 >> (2 * (2 - i/4))) & 1) || 01228 ((v2 >> (2 * (2 - i/4) + 1 )) & 1) ) 01229 { 01230 faceflag[ i ] = 0 ; 01231 } 01232 } 01233 } 01234 01235 if ( faceflag[ i ] ) 01236 { 01237 int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ; 01238 edge[ e0 ] ++ ; 01239 edge[ e1 ] ++ ; 01240 tot ++ ; 01241 } 01242 } 01243 01244 // Removing 1s 01245 int numones = 0 ; 01246 for ( i = 0 ; i < 6 ; i ++ ) 01247 { 01248 if ( edge[ i ] == 1 ) 01249 { 01250 numones ++ ; 01251 } 01252 } 01253 while( numones > 0 ) 01254 { 01255 int e ; 01256 for ( i = 0 ; i < 6 ; i ++ ) 01257 { 01258 if ( edge[ i ] == 1 ) 01259 { 01260 e = i ; 01261 break ; 01262 } 01263 } 01264 /* 01265 if ( edge[ e ] != 1 ) 01266 { 01267 printf("Wrong Again!********\n") ; 01268 } 01269 */ 01270 01271 int f, e2 ; 01272 for ( j = 0 ; j < 4 ; j ++ ) 01273 { 01274 f = edgeFaces[ e ][ j ] ; 01275 if ( faceflag[ f ] ) 01276 { 01277 break ; 01278 } 01279 } 01280 01281 /* 01282 if ( faceflag[ f ] == 0 ) 01283 { 01284 printf("Wrong!********\n") ; 01285 } 01286 */ 01287 01288 if ( faceEdges[ f ][ 0 ] == e ) 01289 { 01290 e2 = faceEdges[ f ][ 1 ] ; 01291 } 01292 else 01293 { 01294 e2 = faceEdges[ f ][ 0 ] ; 01295 } 01296 01297 01298 edge[ e ] -- ; 01299 numones -- ; 01300 edge[ e2 ] -- ; 01301 faceflag[ f ] = 0 ; 01302 tot -- ; 01303 01304 if ( edge[ e2 ] == 1 ) 01305 { 01306 numones ++ ; 01307 } 01308 else if ( edge[ e2 ] == 0 ) 01309 { 01310 numones -- ; 01311 } 01312 } 01313 01314 if ( tot > 0 ) 01315 { 01316 return 1 ; 01317 } 01318 01319 return 0 ; 01320 } 01321 01322 int Volume::hasCompleteSheet( int ox, int oy, int oz ) { 01323 // Returns 1 if it lies in the middle of a sheet 01324 int temp = countIntEuler( ox, oy, oz ) ; 01325 if ( temp > 0 ) 01326 { 01327 return 1 ; 01328 } 01329 else 01330 { 01331 return 0 ; 01332 } 01333 } 01334 01335 //int Volume::hasCompleteSheetSlow( int ox, int oy, int oz ) { 01336 //int i, j ; 01337 //int nx, ny, nz ; 01338 01339 //int edge[6] = { 0,0,0,0,0,0 } ; 01340 //int faceflag[ 12 ] ; 01341 //int tot = 0 ; 01342 01343 //for ( i = 0 ; i < 12 ; i ++ ) 01344 //{ 01345 //faceflag[ i ] = 1 ; 01346 //for ( j = 0 ; j < 4 ; j ++ ) 01347 //{ 01348 //nx = ox + sheetNeighbor[i][j][0] ; 01349 //ny = oy + sheetNeighbor[i][j][1] ; 01350 //nz = oz + sheetNeighbor[i][j][2] ; 01351 01352 //if ( getDataAt( nx, ny, nz ) < 0 ) 01353 //{ 01354 //faceflag[ i ] = 0 ; 01355 //break ; 01356 //} 01357 //} 01358 01359 //if ( faceflag[ i ] ) 01360 //{ 01361 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ; 01362 //edge[ e0 ] ++ ; 01363 //edge[ e1 ] ++ ; 01364 //tot ++ ; 01365 //} 01366 //} 01367 01369 //int numones = 0 ; 01370 //for ( i = 0 ; i < 6 ; i ++ ) 01371 //{ 01372 //if ( edge[ i ] == 1 ) 01373 //{ 01374 //numones ++ ; 01375 //} 01376 //} 01377 //while( numones > 0 ) 01378 //{ 01379 //int e ; 01380 //for ( i = 0 ; i < 6 ; i ++ ) 01381 //{ 01382 //if ( edge[ i ] == 1 ) 01383 //{ 01384 //e = i ; 01385 //break ; 01386 //} 01387 //} 01389 //if ( edge[ e ] != 1 ) 01390 //{ 01391 //printf("Wrong Again!********\n") ; 01392 //} 01393 //*/ 01394 01395 //int f, e2 ; 01396 //for ( j = 0 ; j < 4 ; j ++ ) 01397 //{ 01398 //f = edgeFaces[ e ][ j ] ; 01399 //if ( faceflag[ f ] ) 01400 //{ 01401 //break ; 01402 //} 01403 //} 01404 01406 //if ( faceflag[ f ] == 0 ) 01407 //{ 01408 //printf("Wrong!********\n") ; 01409 //} 01410 //*/ 01411 01412 //if ( faceEdges[ f ][ 0 ] == e ) 01413 //{ 01414 //e2 = faceEdges[ f ][ 1 ] ; 01415 //} 01416 //else 01417 //{ 01418 //e2 = faceEdges[ f ][ 0 ] ; 01419 //} 01420 01421 01422 //edge[ e ] -- ; 01423 //numones -- ; 01424 //edge[ e2 ] -- ; 01425 //faceflag[ f ] = 0 ; 01426 //tot -- ; 01427 01428 //if ( edge[ e2 ] == 1 ) 01429 //{ 01430 //numones ++ ; 01431 //} 01432 //else if ( edge[ e2 ] == 0 ) 01433 //{ 01434 //numones -- ; 01435 //} 01436 //} 01437 01439 //if ( tot > 0 ) 01440 //{ 01442 //{ 01444 //} 01445 //return 1 ; 01446 //} 01447 //else 01448 //{ 01450 //{ 01452 //} 01453 //return 0 ; 01454 //} 01455 //} 01456 int Volume::hasCompleteHelix( int ox, int oy, int oz ) 01457 { 01458 // Returns 1 if it has a complete helix 01459 int i ; 01460 int c1 = 0; 01461 int nx, ny, nz ; 01462 int j ; 01463 01464 for ( i = 0 ; i < 6 ; i ++ ) 01465 { 01466 nx = ox + neighbor6[i][0] ; 01467 ny = oy + neighbor6[i][1] ; 01468 nz = oz + neighbor6[i][2] ; 01469 if ( getDataAt( nx, ny, nz ) >= 0 ) 01470 { 01471 c1 ++ ; 01472 j = i ; 01473 } 01474 01475 } 01476 01477 if ( c1 > 1 ) // || c1 == 0 ) 01478 { 01479 return 1 ; 01480 } 01481 01482 return 0 ; 01483 01484 /* 01485 ox = ox + neighbor6[j][0] ; 01486 oy = oy + neighbor6[j][1] ; 01487 oz = oz + neighbor6[j][2] ; 01488 c1 = 0 ; 01489 for ( i = 0 ; i < 6 ; i ++ ) 01490 { 01491 nx = ox + neighbor6[i][0] ; 01492 ny = oy + neighbor6[i][1] ; 01493 nz = oz + neighbor6[i][2] ; 01494 if ( getDataAt( nx, ny, nz ) >= 0 ) 01495 { 01496 c1 ++ ; 01497 } 01498 01499 } 01500 01501 if ( c1 > 1 ) 01502 { 01503 return 0 ; 01504 } 01505 else 01506 { 01507 return 1 ; 01508 } 01509 */ 01510 } 01511 01512 int Volume::hasCompleteHelix( int ox, int oy, int oz, Volume* fvol ) 01513 { 01514 01515 int i ; 01516 int c1 = 0; 01517 int nx, ny, nz ; 01518 int j ; 01519 01520 for ( i = 0 ; i < 6 ; i ++ ) 01521 { 01522 nx = ox + neighbor6[i][0] ; 01523 ny = oy + neighbor6[i][1] ; 01524 nz = oz + neighbor6[i][2] ; 01525 if ( getDataAt( nx, ny, nz ) >= 0 ) 01526 { 01527 if ( i % 2 == 0 ) 01528 { 01529 nx = ox ; 01530 ny = oy ; 01531 nz = oz ; 01532 } 01533 01534 int val = (int)fvol->getDataAt( nx, ny, nz) ; 01535 if ( (( val >> ( 2 * ( i / 2 ) ) ) & 1) == 0 ) 01536 { 01537 c1 ++ ; 01538 j = i ; 01539 } 01540 } 01541 01542 } 01543 01544 if ( c1 > 1 ) 01545 { 01546 return 1 ; 01547 } 01548 01549 return 0 ; 01550 01551 /* 01552 ox = ox + neighbor6[j][0] ; 01553 oy = oy + neighbor6[j][1] ; 01554 oz = oz + neighbor6[j][2] ; 01555 c1 = 0 ; 01556 for ( i = 0 ; i < 6 ; i ++ ) 01557 { 01558 nx = ox + neighbor6[i][0] ; 01559 ny = oy + neighbor6[i][1] ; 01560 nz = oz + neighbor6[i][2] ; 01561 if ( getDataAt( nx, ny, nz ) >= 0 ) 01562 { 01563 c1 ++ ; 01564 } 01565 01566 } 01567 01568 if ( c1 > 1 ) 01569 { 01570 return 0 ; 01571 } 01572 else 01573 { 01574 return 1 ; 01575 } 01576 */ 01577 } 01578 01579 int Volume::isHelixEnd( int ox, int oy, int oz, Volume* nvol ) { 01580 // Returns 1 if it is a curve endpoint 01581 int i ; 01582 int c1 = 0 , c2 = 0 ; 01583 int nx, ny, nz ; 01584 01585 for ( i = 0 ; i < 6 ; i ++ ) 01586 { 01587 nx = ox + neighbor6[i][0] ; 01588 ny = oy + neighbor6[i][1] ; 01589 nz = oz + neighbor6[i][2] ; 01590 01591 double val = getDataAt( nx, ny, nz ) ; 01592 01593 if ( val >= 0 ) 01594 { 01595 c1 ++ ; 01596 if ( val > 0 && val < MAX_ERODE && nvol->getDataAt( nx, ny, nz ) == 0 ) 01597 { 01598 c2 ++ ; 01599 } 01600 } 01601 01602 } 01603 01604 if ( c1 == 1 && c2 == 1 ) 01605 { 01606 return 1 ; 01607 } 01608 01609 return 0 ; 01610 } 01611 01612 //int Volume::isFeature18( int ox, int oy, int oz ) 01613 //{ 01617 //if ( getDataAt( ox, oy, oz ) <= 0 ) 01618 //{ 01619 //return 0 ; 01620 //} 01621 01622 //for ( int i = 0 ; i < 3 ; i ++ ) 01623 //for ( int j = 0 ; j < 3 ; j ++ ) 01624 //for ( int k = 0 ; k < 3 ; k ++ ) 01625 //{ 01626 //if ( i == 1 || j == 1 || k == 1 ) 01627 //{ 01628 //if( getDataAt( ox - 1 + i, oy - 1 + j, oz - 1 + k ) == 0 ) 01629 //{ 01630 //return 0 ; 01631 //} 01632 //} 01633 //} 01634 //return 1 ; 01635 //} 01636 01637 //int Volume::isEdgeEnd( int ox, int oy, int oz ) 01638 //{ 01639 //int i ; 01640 //int c1 = 0 ; 01641 //int nx, ny, nz ; 01642 01643 //for ( i = 0 ; i < 6 ; i ++ ) 01644 //{ 01645 //nx = ox + neighbor6[i][0] ; 01646 //ny = oy + neighbor6[i][1] ; 01647 //nz = oz + neighbor6[i][2] ; 01648 01649 //double val = getDataAt( nx, ny, nz ) ; 01650 01651 //if ( val >= 0 ) 01652 //{ 01653 //c1 ++ ; 01654 //} 01655 01656 //} 01657 01658 //if ( c1 == 1 ) 01659 //{ 01660 //return 1 ; 01661 //} 01662 01663 //return 0 ; 01664 //} 01665 01666 //int Volume::isFaceEnd( int ox, int oy, int oz ) 01667 //{ 01669 01670 //int i, j, k ; 01671 //int nx, ny, nz ; 01672 01673 //double vox[3][3][3] ; 01674 //for ( i = -1 ; i < 2 ; i ++ ) 01675 //for ( j = -1 ; j < 2 ; j ++ ) 01676 //for ( k = -1 ; k < 2 ; k ++ ) 01677 //{ 01678 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ; 01679 //} 01680 01681 //int edge[6] = { 4,4,4,4,4,4 } ; 01682 //int edge2[6] = { 4,4,4,4,4,4 } ; 01683 01684 //for ( i = 0 ; i < 12 ; i ++ ) 01685 //{ 01686 //for ( j = 0 ; j < 4 ; j ++ ) 01687 //{ 01688 //nx = 1 + sheetNeighbor[i][j][0] ; 01689 //ny = 1 + sheetNeighbor[i][j][1] ; 01690 //nz = 1 + sheetNeighbor[i][j][2] ; 01691 01692 //if ( vox[nx][ny][nz] < 0 ) 01693 //{ 01694 //edge[ faceEdges[ i ][ 0 ] ] -- ; 01695 //edge[ faceEdges[ i ][ 1 ] ] -- ; 01696 //break ; 01697 //} 01698 //} 01699 //for ( j = 0 ; j < 4 ; j ++ ) 01700 //{ 01701 //nx = 1 + sheetNeighbor[i][j][0] ; 01702 //ny = 1 + sheetNeighbor[i][j][1] ; 01703 //nz = 1 + sheetNeighbor[i][j][2] ; 01704 01705 //if ( vox[nx][ny][nz] < 2 ) 01706 //{ 01707 //edge2[ faceEdges[ i ][ 0 ] ] -- ; 01708 //edge2[ faceEdges[ i ][ 1 ] ] -- ; 01709 //break ; 01710 //} 01711 //} 01712 //} 01713 01714 //for ( i = 0 ; i < 6 ; i ++ ) 01715 //{ 01716 //if ( edge[ i ] == 1 && edge2[ i ] == 1 ) 01717 //{ 01718 //return 1 ; 01719 //} 01720 //} 01721 01722 //return 0 ; 01723 //} 01724 01725 //int Volume::isNoise( int ox, int oy, int oz, int noise ) 01726 //{ 01727 //if ( getDataAt(ox,oy,oz) == 1 ) 01728 //{ 01729 //return 1 ; 01730 //} 01731 01732 //if ( noise > 0 ) 01733 //{ 01734 //for ( int i = 0 ; i < 6 ; i ++ ) 01735 //{ 01736 //int nx = ox + neighbor6[i][0] ; 01737 //int ny = oy + neighbor6[i][1] ; 01738 //int nz = oz + neighbor6[i][2] ; 01739 01740 //if ( getDataAt( nx, ny, nz ) > 0 ) 01741 //{ 01742 //if ( isNoise( nx, ny, nz, noise - 1 ) ) 01743 //{ 01744 //return 1 ; 01745 //} 01746 //} 01747 //} 01748 //} 01749 01750 //return 0 ; 01751 01752 //} 01753 01754 //int Volume::isNoiseHelixEnd( int ox, int oy, int oz ) 01755 //{ 01756 01757 //int i ; 01758 //int c1 = 0 , c2 = 0 ; 01759 //int nx, ny, nz ; 01760 01761 //for ( i = 0 ; i < 6 ; i ++ ) 01762 //{ 01763 //nx = ox + neighbor6[i][0] ; 01764 //ny = oy + neighbor6[i][1] ; 01765 //nz = oz + neighbor6[i][2] ; 01766 01767 //double val = getDataAt( nx, ny, nz ) ; 01768 01769 //if ( val >= 0 ) 01770 //{ 01771 //c1 ++ ; 01772 //if ( val > 0 && val < MAX_ERODE ) 01773 //{ 01774 //c2 ++ ; 01775 //} 01776 //} 01777 01778 //} 01779 01780 //if ( c1 == 1 && c2 == 0 ) 01781 //{ 01782 //return 1 ; 01783 //} 01784 01785 //return 0 ; 01786 //} 01787 01788 01789 int Volume::isHelixEnd( int ox, int oy, int oz ) { 01790 01791 int i ; 01792 int c1 = 0 , c2 = 0 ; 01793 int nx, ny, nz ; 01794 01795 for ( i = 0 ; i < 6 ; i ++ ) 01796 { 01797 nx = ox + neighbor6[i][0] ; 01798 ny = oy + neighbor6[i][1] ; 01799 nz = oz + neighbor6[i][2] ; 01800 01801 double val = getDataAt( nx, ny, nz ) ; 01802 01803 if ( val >= 0 ) 01804 { 01805 c1 ++ ; 01806 if ( getNumNeighbor6(nx,ny,nz) < 6 ) // if ( val > 0 && val < MAX_ERODE ) 01807 { 01808 c2 ++ ; 01809 } 01810 } 01811 01812 } 01813 01814 if ( c1 == 1 && c2 == 1 ) 01815 { 01816 return 1 ; 01817 } 01818 01819 return 0 ; 01820 } 01821 int Volume::isSheetEnd( int ox, int oy, int oz, Volume* nvol ) 01822 { 01823 // Returns 1 if it contains a sheet boundary. Noise-resistant 01824 int i, j ; 01825 int nx, ny, nz ; 01826 01827 int edge[6] = { 0,0,0,0,0,0 } ; 01828 int faceflag[ 12 ] ; 01829 int hasFeatureFace = 0 ; 01830 int tot = 0 ; 01831 01832 for ( i = 0 ; i < 12 ; i ++ ) 01833 { 01834 faceflag[ i ] = 1 ; 01835 int hasFeature = 1 ; 01836 01837 for ( j = 0 ; j < 4 ; j ++ ) 01838 { 01839 nx = ox + sheetNeighbor[i][j][0] ; 01840 ny = oy + sheetNeighbor[i][j][1] ; 01841 nz = oz + sheetNeighbor[i][j][2] ; 01842 01843 if ( getDataAt( nx, ny, nz ) == 0 || nvol->getDataAt( nx, ny, nz ) == 1 ) 01844 { 01845 hasFeature = 0 ; 01846 } 01847 if ( getDataAt( nx, ny, nz ) < 0 ) 01848 { 01849 faceflag[ i ] = 0 ; 01850 break ; 01851 } 01852 } 01853 if ( faceflag[ i ] == 1 && hasFeature ) 01854 { 01855 hasFeatureFace ++ ; 01856 // return 0 ; 01857 } 01858 01859 if ( faceflag[ i ] ) 01860 { 01861 int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ; 01862 edge[ e0 ] ++ ; 01863 edge[ e1 ] ++ ; 01864 tot ++ ; 01865 } 01866 } 01867 01868 if ( tot == 0 || hasFeatureFace == 0 ) 01869 { 01870 return 0 ; 01871 } 01872 01873 // Removing 1s 01874 int numones = 0 ; 01875 for ( i = 0 ; i < 6 ; i ++ ) 01876 { 01877 if ( edge[ i ] == 1 ) 01878 { 01879 numones ++ ; 01880 } 01881 } 01882 while( numones > 0 ) 01883 { 01884 int e ; 01885 for ( i = 0 ; i < 6 ; i ++ ) 01886 { 01887 if ( edge[ i ] == 1 ) 01888 { 01889 e = i ; 01890 break ; 01891 } 01892 } 01893 /* 01894 if ( edge[ e ] != 1 ) 01895 { 01896 printf("Wrong Again!********\n") ; 01897 } 01898 */ 01899 01900 int f, e2 ; 01901 for ( j = 0 ; j < 4 ; j ++ ) 01902 { 01903 f = edgeFaces[ e ][ j ] ; 01904 if ( faceflag[ f ] ) 01905 { 01906 break ; 01907 } 01908 } 01909 01910 /* 01911 if ( faceflag[ f ] == 0 ) 01912 { 01913 printf("Wrong!********\n") ; 01914 } 01915 */ 01916 01917 if ( faceEdges[ f ][ 0 ] == e ) 01918 { 01919 e2 = faceEdges[ f ][ 1 ] ; 01920 } 01921 else 01922 { 01923 e2 = faceEdges[ f ][ 0 ] ; 01924 } 01925 01926 01927 edge[ e ] -- ; 01928 numones -- ; 01929 edge[ e2 ] -- ; 01930 faceflag[ f ] = 0 ; 01931 tot -- ; 01932 01933 if ( edge[ e2 ] == 1 ) 01934 { 01935 numones ++ ; 01936 } 01937 else if ( edge[ e2 ] == 0 ) 01938 { 01939 numones -- ; 01940 } 01941 } 01942 01943 if ( tot > 0 ) 01944 { 01945 return 0 ; 01946 } 01947 01948 return 1 ; 01949 } 01950 01951 //int Volume::getNumFaces( int ox, int oy, int oz ) 01952 //{ 01953 //int i, j ; 01954 //int nx, ny, nz ; 01955 01956 //int faces = 12 ; 01957 //for ( i = 0 ; i < 12 ; i ++ ) 01958 //{ 01959 //for ( j = 0 ; j < 4 ; j ++ ) 01960 //{ 01961 //nx = ox + sheetNeighbor[i][j][0] ; 01962 //ny = oy + sheetNeighbor[i][j][1] ; 01963 //nz = oz + sheetNeighbor[i][j][2] ; 01964 01965 //if ( getDataAt( nx, ny, nz ) < 0 ) 01966 //{ 01967 //faces -- ; 01968 //break ; 01969 //} 01970 //} 01971 //} 01972 01973 //return faces ; 01974 //} 01975 01976 //int Volume::getNumCells( int ox, int oy, int oz ) 01977 //{ 01978 //int i, j, k ; 01979 01980 //int cells = 0 ; 01981 01982 //for ( i = -1 ; i < 1 ; i ++ ) 01983 //for ( j = -1 ; j < 1 ; j ++ ) 01984 //for ( k = -1 ; k < 1 ; k ++ ) 01985 //{ 01986 //if ( hasCell( ox + i, oy + j, oz + k ) ) 01987 //{ 01988 //cells ++ ; 01989 //} 01990 //} 01991 01993 //return cells ; 01994 //} 01995 01996 01997 //int Volume::getNumIsolatedEdges( int ox, int oy, int oz ) 01998 //{ 01999 //int i, j, k ; 02000 //int nx, ny, nz ; 02001 02002 //double vox[3][3][3] ; 02003 //for ( i = -1 ; i < 2 ; i ++ ) 02004 //for ( j = -1 ; j < 2 ; j ++ ) 02005 //for ( k = -1 ; k < 2 ; k ++ ) 02006 //{ 02007 //vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ; 02008 //} 02009 02010 02011 //int edge[6] = { 0,0,0,0,0,0 } ; 02012 02013 //for ( i = 0 ; i < 12 ; i ++ ) 02014 //{ 02015 //int flag = 1 ; 02016 //for ( j = 0 ; j < 4 ; j ++ ) 02017 //{ 02018 //nx = 1 + sheetNeighbor[i][j][0] ; 02019 //ny = 1 + sheetNeighbor[i][j][1] ; 02020 //nz = 1 + sheetNeighbor[i][j][2] ; 02021 02022 //if ( vox[nx][ny][nz] < 0 ) 02023 //{ 02024 //flag = 0 ; 02025 //break ; 02026 //} 02027 //} 02028 02029 //if ( flag ) 02030 //{ 02031 //int e0 = faceEdges[ i ][ 0 ], e1 = faceEdges[ i ][ 1 ] ; 02032 //edge[ e0 ] ++ ; 02033 //edge[ e1 ] ++ ; 02034 //} 02035 //} 02036 02037 //int edges = 0 ; 02038 //for ( i = 0 ; i < 6 ; i ++ ) 02039 //{ 02040 //if ( edge[i] ) 02041 //{ 02042 //continue ; 02043 //} 02044 02045 //nx = 1 + neighbor6[i][0] ; 02046 //ny = 1 + neighbor6[i][1] ; 02047 //nz = 1 + neighbor6[i][2] ; 02048 02049 //if ( vox[nx][ny][nz] >= 0 ) 02050 //{ 02051 //edges ++ ; 02052 //} 02053 02054 //} 02055 02056 //return edges ; 02057 //} 02058 02059 //int Volume::getNumIsolatedFaces( int ox, int oy, int oz ) 02060 //{ 02061 //int i, j, k ; 02062 //int nx, ny, nz ; 02063 02064 //int faces = 0 ; 02065 //int cellflag[ 8 ] ; 02066 02067 //int ct = 0 ; 02068 //for ( i = -1 ; i < 1 ; i ++ ) 02069 //for ( j = -1 ; j < 1 ; j ++ ) 02070 //for ( k = -1 ; k < 1 ; k ++ ) 02071 //{ 02072 //if ( hasCell( ox + i, oy + j, oz + k ) ) 02073 //{ 02074 //cellflag[ ct ] = 1 ; 02075 //} 02076 //else 02077 //{ 02078 //cellflag[ ct ] = 0 ; 02079 //} 02080 //ct ++ ; 02081 //} 02082 02083 //for ( i = 0 ; i < 12 ; i ++ ) 02084 //{ 02085 //int hasFace = 1 ; 02086 //for ( j = 0 ; j < 4 ; j ++ ) 02087 //{ 02088 //nx = ox + sheetNeighbor[i][j][0] ; 02089 //ny = oy + sheetNeighbor[i][j][1] ; 02090 //nz = oz + sheetNeighbor[i][j][2] ; 02091 02092 //if ( getDataAt( nx, ny, nz ) < 0 ) 02093 //{ 02094 //hasFace = 0 ; 02095 //break ; 02096 //} 02097 //} 02098 02099 //if ( hasFace ) 02100 //{ 02101 //if (cellflag[ faceCells[i][0] ] == 0 && cellflag[ faceCells[i][1] ] == 0 ) 02102 //{ 02103 //faces ++ ; 02104 //} 02105 //} 02106 //} 02107 02109 //return faces ; 02110 //} 02111 02112 //int Volume::isFeatureFace2( int ox, int oy, int oz ) 02113 //{ 02114 //int i ; 02115 //int nx, ny, nz ; 02116 02117 //for ( i = 0 ; i < 6 ; i ++ ) 02118 //{ 02119 //nx = ox + neighbor6[i][0] ; 02120 //ny = oy + neighbor6[i][1] ; 02121 //nz = oz + neighbor6[i][2] ; 02122 02123 //double val = getDataAt( nx, ny, nz ) ; 02124 02125 //if ( val == 0 ) 02126 //{ 02127 //return 0 ; 02128 //} 02129 02130 //} 02131 02132 //return 1 ; 02133 //} 02134 02135 int Volume::isFeatureFace( int ox, int oy, int oz ) 02136 { 02137 // return 1 ; 02138 02139 int i, j ; 02140 int nx, ny, nz ; 02141 02142 int faces = 12 ; 02143 for ( i = 0 ; i < 12 ; i ++ ) 02144 { 02145 int ct = 0 ; 02146 for ( j = 0 ; j < 4 ; j ++ ) 02147 { 02148 nx = ox + sheetNeighbor[i][j][0] ; 02149 ny = oy + sheetNeighbor[i][j][1] ; 02150 nz = oz + sheetNeighbor[i][j][2] ; 02151 02152 if ( getDataAt( nx, ny, nz ) < 0 ) 02153 { 02154 ct = -1 ; 02155 break ; 02156 } 02157 else if ( getNumNeighbor6( nx, ny, nz ) == 6 ) 02158 { 02159 ct = -1 ; 02160 break ; 02161 } 02162 // else if ( getDataAt( nx, ny, nz ) == 0 ) 02163 // { 02164 // ct ++ ; 02165 // } 02166 02167 02168 } 02169 if ( ct == -1 || ct >= 1 ) 02170 { 02171 faces -- ; 02172 } 02173 } 02174 02175 if ( faces > 0 ) 02176 { 02177 return 1 ; 02178 } 02179 return 0 ; 02180 02181 } 02182 02183 //int Volume::hasFeatureFace( int ox, int oy, int oz ) 02184 //{ 02185 //int i, j, k ; 02186 //int nx, ny, nz ; 02187 02188 //int faceflag ; 02189 //int cellflag[ 8 ] ; 02190 02192 //int ct = 0 ; 02193 //for ( i = -1 ; i < 1 ; i ++ ) 02194 //for ( j = -1 ; j < 1 ; j ++ ) 02195 //for ( k = -1 ; k < 1 ; k ++ ) 02196 //{ 02197 //if ( hasCell( ox + i, oy + j, oz + k ) ) 02198 //{ 02199 //cellflag[ ct ] = 1 ; 02200 //} 02201 //else 02202 //{ 02203 //cellflag[ ct ] = 0 ; 02204 //} 02205 //ct ++ ; 02206 //} 02207 02209 //for ( i = 0 ; i < 12 ; i ++ ) 02210 //{ 02211 //faceflag = 1 ; 02212 //for ( j = 0 ; j < 4 ; j ++ ) 02213 //{ 02214 //nx = ox + sheetNeighbor[i][j][0] ; 02215 //ny = oy + sheetNeighbor[i][j][1] ; 02216 //nz = oz + sheetNeighbor[i][j][2] ; 02217 02218 //if ( getDataAt( nx, ny, nz ) < 0 ) 02219 //{ 02220 //faceflag = 0 ; 02221 //break ; 02222 //} 02223 //} 02224 02225 //if ( faceflag ) 02226 //{ 02227 //if ( cellflag[ faceCells[i][0] ] == 0 && cellflag[ faceCells[i][1] ] == 0 ) 02228 //{ 02229 //return 1 ; 02230 //} 02231 //} 02232 //} 02233 02234 //return 0 ; 02235 02236 //} 02237 02238 int Volume::isSheetEnd( int ox, int oy, int oz ) 02239 { 02240 return ( ( hasCompleteSheet(ox,oy,oz) == 0 ) && isFeatureFace(ox,oy,oz) ) ; 02241 02243 // 02244 //int i, j, k ; 02245 //int nx, ny, nz ; 02246 02247 //double vox[3][3][3] ; 02248 //for ( i = -1 ; i < 2 ; i ++ ) 02249 // for ( j = -1 ; j < 2 ; j ++ ) 02250 // for ( k = -1 ; k < 2 ; k ++ ) 02251 // { 02252 // vox[ i + 1 ][ j + 1 ][ k + 1 ] = getDataAt( ox + i, oy + j, oz + k ) ; 02253 // } 02254 02255 //int edge[6] = { 4,4,4,4,4,4 } ; 02256 //int edge2[6] = { 4,4,4,4,4,4 } ; 02257 02258 //for ( i = 0 ; i < 12 ; i ++ ) 02259 //{ 02260 // for ( j = 0 ; j < 4 ; j ++ ) 02261 // { 02262 // nx = 1 + sheetNeighbor[i][j][0] ; 02263 // ny = 1 + sheetNeighbor[i][j][1] ; 02264 // nz = 1 + sheetNeighbor[i][j][2] ; 02265 02266 // if ( vox[nx][ny][nz] < 0 ) 02267 // { 02268 // edge[ faceEdges[ i ][ 0 ] ] -- ; 02269 // edge[ faceEdges[ i ][ 1 ] ] -- ; 02270 // break ; 02271 // } 02272 // } 02273 02274 // for ( j = 0 ; j < 4 ; j ++ ) 02275 // { 02276 // nx = 1 + sheetNeighbor[i][j][0] ; 02277 // ny = 1 + sheetNeighbor[i][j][1] ; 02278 // nz = 1 + sheetNeighbor[i][j][2] ; 02279 02280 // if ( vox[nx][ny][nz] <= 0 ) 02281 // { 02282 // edge2[ faceEdges[ i ][ 0 ] ] -- ; 02283 // edge2[ faceEdges[ i ][ 1 ] ] -- ; 02284 // break ; 02285 // } 02286 // } 02287 //} 02288 02289 // 02291 //for ( i = 0 ; i < 6 ; i ++ ) 02292 //{ 02293 // nx = 1 + neighbor6[i][0] ; 02294 // ny = 1 + neighbor6[i][1] ; 02295 // nz = 1 + neighbor6[i][2] ; 02296 // if ( edge[i] == 0 && vox[nx][ny][nz] >= 0 ) 02297 // { 02298 // return 0 ; 02299 // } 02300 //} 02301 //*/ 02302 // 02303 02304 02305 //for ( i = 0 ; i < 6 ; i ++ ) 02306 //{ 02307 // if ( edge[ i ] == 1 && edge2[ i ] == 1 ) 02308 // { 02309 // return 1 ; 02310 // } 02311 //} 02312 02313 //return 0 ; 02314 } 02315 02316 int Volume::isSimple( int ox, int oy, int oz ) 02317 { 02318 /* Test if this is a simple voxel */ 02319 // int flag = 0 ; 02320 double vox[3][3][3] ; 02321 02322 int i, j, k ; 02323 for ( i = -1 ; i < 2 ; i ++ ) 02324 for ( j = -1 ; j < 2 ; j ++ ) 02325 for ( k = -1 ; k < 2 ; k ++ ) 02326 { 02327 double tval = getDataAt( ox + i, oy + j, oz + k ) ; 02328 02329 /* 02330 if ( tval == 0 || tval > (va ) 02331 { 02332 flag = 1 ; 02333 } 02334 */ 02335 /* 02336 if ( tval < 0 && tval == - curwid ) 02337 { 02338 printf("Here %d", (int)-tval) ; 02339 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ; 02340 } 02341 else 02342 */ 02343 { 02344 vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ; 02345 } 02346 } 02347 02348 /* Debugging 02349 printf("{") ; 02350 for ( i = 0 ; i < 3 ; i ++ ) 02351 { 02352 if ( i ) printf(",") ; 02353 printf("{") ; 02354 for ( j = 0 ; j < 3 ; j ++ ) 02355 { 02356 if ( j ) printf(",") ; 02357 printf("{") ; 02358 for ( k = 0 ; k < 3 ; k ++ ) 02359 { 02360 if ( k ) printf(",") ; 02361 printf("%d", (vox[i][j][k] >=0 ? 1: 0)); 02362 } 02363 printf("}") ; 02364 } 02365 printf("}") ; 02366 } 02367 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ; 02368 */ 02369 02370 if ( countInt( vox ) == 1 && countExt( vox ) == 1 ) 02371 { 02372 return 1 ; 02373 } 02374 else 02375 { 02376 return 0 ; 02377 } 02378 } 02379 02380 02381 int Volume::isPiercable( int ox, int oy, int oz ) 02382 { 02383 /* Test if this is a simple voxel */ 02384 // int flag = 0 ; 02385 double vox[3][3][3] ; 02386 02387 int i, j, k ; 02388 for ( i = -1 ; i < 2 ; i ++ ) 02389 for ( j = -1 ; j < 2 ; j ++ ) 02390 for ( k = -1 ; k < 2 ; k ++ ) 02391 { 02392 double tval = getDataAt( ox + i, oy + j, oz + k ) ; 02393 02394 /* 02395 if ( tval == 0 || tval > (va ) 02396 { 02397 flag = 1 ; 02398 } 02399 */ 02400 /* 02401 if ( tval < 0 && tval == - curwid ) 02402 { 02403 printf("Here %d", (int)-tval) ; 02404 vox[ i + 1 ][ j + 1 ][ k + 1 ] = - tval ; 02405 } 02406 else 02407 */ 02408 { 02409 vox[ i + 1 ][ j + 1 ][ k + 1 ] = tval ; 02410 } 02411 } 02412 02413 /* Debugging 02414 printf("{") ; 02415 for ( i = 0 ; i < 3 ; i ++ ) 02416 { 02417 if ( i ) printf(",") ; 02418 printf("{") ; 02419 for ( j = 0 ; j < 3 ; j ++ ) 02420 { 02421 if ( j ) printf(",") ; 02422 printf("{") ; 02423 for ( k = 0 ; k < 3 ; k ++ ) 02424 { 02425 if ( k ) printf(",") ; 02426 printf("%d", (vox[i][j][k] >=0 ? 1: 0)); 02427 } 02428 printf("}") ; 02429 } 02430 printf("}") ; 02431 } 02432 printf("} Int: %d, Ext: %d\n", countInt( vox ), countExt( vox )) ; 02433 */ 02434 02435 if ( countInt( vox ) == 1 && countExt( vox ) != 1 ) 02436 { 02437 return 1 ; 02438 } 02439 else 02440 { 02441 return 0 ; 02442 } 02443 } 02444 02445 02446 //int Volume::isSimple2( int v[3][3][3] ) 02447 //{ 02449 //double vox[3][3][3] ; 02450 02451 //int i, j, k ; 02452 //for ( i = 0 ; i < 3 ; i ++ ) 02453 //for ( j = 0 ; j < 3 ; j ++ ) 02454 //for ( k = 0 ; k < 3 ; k ++ ) 02455 //{ 02456 //if ( v[i][j][k] == 0 ) 02457 //{ 02458 //vox[ i ][ j ][ k ] = 1 ; 02459 //} 02460 //else 02461 //{ 02462 //vox[i][j][k] = -1 ; 02463 //} 02464 //} 02465 //if ( countInt( vox ) == 1 && countExt( vox ) == 1 ) 02466 //{ 02467 //return 1 ; 02468 //} 02469 //else 02470 //{ 02471 //printf("Int: %d Ext: %d\n",countInt( vox ),countExt( vox ) ); 02472 //return 0 ; 02473 //} 02474 //} 02475 02476 //int Volume::getNumPotComplex3( int ox, int oy, int oz ) 02477 //{ 02479 02480 02481 //int i, j, k ; 02482 //if ( ! isSimple( ox, oy, oz ) || isSheetEnd( ox, oy, oz ) ) 02483 //{ 02484 //return 60 ; 02485 //} 02486 //double val = getDataAt( ox, oy, oz ) ; 02487 02488 //int nx, ny, nz ; 02489 02490 02491 //int numSimple = 0 ; 02492 //for ( i = -1 ; i < 2 ; i ++ ) 02493 //for ( j = -1 ; j < 2 ; j ++ ) 02494 //for ( k = -1 ; k < 2 ; k ++ ) 02495 //{ 02496 //nx = ox + i ; 02497 //ny = oy + j ; 02498 //nz = oz + k ; 02499 //if ( getDataAt( nx, ny, nz ) > 0 ) 02500 //{ 02501 //if ( isSimple( nx, ny, nz ) && ! isSheetEnd( nx, ny, nz ) ) 02502 //{ 02503 //numSimple ++ ; 02504 //} 02505 //} 02506 //} 02507 //numSimple -- ; 02508 02509 //setDataAt( ox, oy, oz, -val ) ; 02510 02511 //int numPotSimple = 0 ; 02512 //for ( i = -1 ; i < 2 ; i ++ ) 02513 //for ( j = -1 ; j < 2 ; j ++ ) 02514 //for ( k = -1 ; k < 2 ; k ++ ) 02515 //{ 02516 //nx = ox + i ; 02517 //ny = oy + j ; 02518 //nz = oz + k ; 02519 //if ( getDataAt( nx, ny, nz ) > 0 ) 02520 //{ 02521 //if ( isSimple( nx, ny, nz ) && ! isSheetEnd( nx, ny, nz ) ) 02522 //{ 02523 //numPotSimple ++ ; 02524 //} 02525 //} 02526 //} 02527 02528 02529 //setDataAt( ox, oy, oz, val ) ; 02530 02531 02532 02533 02534 //return 30 - ( numPotSimple - numSimple ) ; 02535 //} 02536 02537 //int Volume::getNumPotComplex4( int ox, int oy, int oz ) 02538 //{ 02539 02540 //int cells = getNumCells(ox,oy,oz) ; 02542 //int faces = getNumFaces(ox,oy,oz) ; 02544 02545 //return ( cells * 6 - 2 * faces ) ; // + 2 * ( faces * 4 - 4 * iedges ) ; 02546 //} 02547 02548 int Volume::getNumPotComplex( int ox, int oy, int oz ) 02549 { 02550 //return 0 ; 02551 02552 int i; 02553 double val = getDataAt( ox, oy, oz ) ; 02554 if ( val <= 0 ) 02555 { 02556 // return 70 ; 02557 } 02558 02559 // return getNumNeighbor6( ox, oy, oz ) ; 02560 02561 // int v = ((getNumNeighbor6( ox, oy, oz ) & 255) << 24) ; 02562 //int v = 0 ; 02563 02564 int rvalue = 0, nx, ny, nz ; 02565 setDataAt( ox, oy, oz, -val ) ; 02566 02567 /* 02568 for ( i = -1 ; i < 2 ; i ++ ) 02569 for ( j = -1 ; j < 2 ; j ++ ) 02570 for ( k = -1 ; k < 2 ; k ++ ) 02571 { 02572 nx = ox + i ; 02573 ny = oy + j ; 02574 nz = oz + k ; 02575 if ( getDataAt( nx, ny, nz ) == val ) 02576 { 02577 if ( isSheetEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) ) 02578 { 02579 rvalue ++ ; 02580 } 02581 } 02582 } 02583 */ 02584 02585 for ( i = 0 ; i < 6 ; i ++ ) 02586 { 02587 nx = ox + neighbor6[i][0] ; 02588 ny = oy + neighbor6[i][1] ; 02589 nz = oz + neighbor6[i][2] ; 02590 02591 if ( getDataAt(nx,ny,nz) >= 0 ) 02592 { 02593 int num = getNumNeighbor6( nx, ny, nz ) ; 02594 if ( num > rvalue ) 02595 { 02596 rvalue = num ; 02597 } 02598 } 02599 } 02600 02601 02602 setDataAt( ox, oy, oz, val ) ; 02603 02604 return rvalue + getNumNeighbor6( ox, oy, oz ) * 10 ; 02605 /* 02606 int v = (((rvalue + getNumNeighbor6( ox, oy, oz ) * 10) & 255) << 24) ; 02607 v |= ( ( ox & 255 ) << 16 ) ; 02608 v |= ( ( oy & 255 ) << 8 ) ; 02609 v |= ( ( oz & 255 ) ) ; 02610 return v ; 02611 */ 02612 02613 } 02614 02615 int Volume::getNumPotComplex2( int ox, int oy, int oz ) 02616 { 02617 return getNumPotComplex( ox, oy, oz ) ; 02618 02619 //int i, j, k ; 02620 //double val = getDataAt( ox, oy, oz ) ; 02621 //if ( val <= 0 ) 02622 //{ 02623 // return 0 ; 02624 //} 02625 02626 //int rvalue = 0, nx, ny, nz ; 02627 //setDataAt( ox, oy, oz, -val ) ; 02628 02629 //for ( i = -1 ; i < 2 ; i ++ ) 02630 // for ( j = -1 ; j < 2 ; j ++ ) 02631 // for ( k = -1 ; k < 2 ; k ++ ) 02632 // { 02633 // nx = ox + i ; 02634 // ny = oy + j ; 02635 // nz = oz + k ; 02636 // if ( getDataAt( nx, ny, nz ) == val ) 02637 // { 02638 // if ( isHelixEnd( nx, ny, nz) || ! isSimple ( nx, ny, nz ) ) 02639 // { 02640 // rvalue ++ ; 02641 // } 02642 // } 02643 // } 02644 02645 //setDataAt( ox, oy, oz, val ) ; 02646 02647 //return rvalue + getNumNeighbor6( ox, oy, oz ) * 30 ; 02648 } 02649 02650 //int Volume::getNumNeighbor( int ox, int oy, int oz ) 02651 //{ 02652 //int rvalue = 0 ; 02653 //double val = getDataAt( ox, oy, oz ) ; 02654 //for ( int i = 0 ; i < 6 ; i ++ ) 02655 //{ 02656 //int nx = ox + neighbor6[i][0] ; 02657 //int ny = oy + neighbor6[i][1] ; 02658 //int nz = oz + neighbor6[i][2] ; 02659 02660 //if ( getDataAt( nx, ny, nz ) == val ) 02661 //{ 02662 //rvalue ++ ; 02663 //} 02664 //} 02666 //for ( int i = -1 ; i < 2 ; i ++ ) 02667 //for ( int j = -1 ; j < 2 ; j ++ ) 02668 //for ( int k = -1 ; k < 2 ; k ++ ) 02669 //{ 02670 //int nx = ox + i ; 02671 //int ny = oy + j ; 02672 //int nz = oz + k ; 02673 02674 //if ( getDataAt( nx, ny, nz ) == val ) 02675 //{ 02676 //rvalue ++ ; 02677 //} 02678 //} 02679 //*/ 02680 //return rvalue ; 02681 //} 02682 02683 02684 //void Volume::setScoreNeighbor( GridQueue* queue ) 02685 //{ 02687 //gridQueueEle* ele = queue->getHead() ; 02688 //while ( ele != NULL ) 02689 //{ 02690 //ele->score = getNumNeighbor( ele->x, ele->y, ele->z ) ; 02691 //ele = ele->next ; 02692 //} 02693 02694 //queue->sort( queue->getNumElements() ) ; 02695 //} 02696 02697 02698 int Volume::components6( int vox[3][3][3] ) 02699 { 02700 // Stupid flood fill 02701 int tot = 0 ; 02702 int queue[27][3] ; 02703 int vis[3][3][3] ; 02704 int head = 0, tail = 1 ; 02705 int i, j, k ; 02706 for ( i = 0 ; i < 3 ; i ++ ) 02707 for ( j = 0 ; j < 3 ; j ++ ) 02708 for ( k = 0 ; k < 3 ; k ++ ) 02709 { 02710 vis[i][j][k] = 0 ; 02711 if ( vox[i][j][k] ) 02712 { 02713 if ( tot == 0 ) 02714 { 02715 queue[0][0] = i ; 02716 queue[0][1] = j ; 02717 queue[0][2] = k ; 02718 vis[i][j][k] = 1 ; 02719 } 02720 tot ++ ; 02721 } 02722 } 02723 if ( tot == 0 ) 02724 { 02725 return 0 ; 02726 } 02727 // printf("total: %d\n", tot) ; 02728 02729 int ct = 1 ; 02730 while ( head != tail ) 02731 { 02732 int x = queue[head][0] ; 02733 int y = queue[head][1] ; 02734 int z = queue[head][2] ; 02735 head ++ ; 02736 02737 for ( i = 0 ; i < 6 ; i ++ ) 02738 { 02739 int nx = x + neighbor6[i][0] ; 02740 int ny = y + neighbor6[i][1] ; 02741 int nz = z + neighbor6[i][2] ; 02742 if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 ) 02743 { 02744 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] ) 02745 { 02746 queue[tail][0] = nx ; 02747 queue[tail][1] = ny ; 02748 queue[tail][2] = nz ; 02749 tail ++ ; 02750 vis[nx][ny][nz] = 1 ; 02751 ct ++ ; 02752 } 02753 } 02754 } 02755 } 02756 02757 if ( ct == tot ) 02758 { 02759 return 1 ; 02760 } 02761 else 02762 { 02763 return 2 ; 02764 } 02765 02766 } 02767 int Volume::components26( int vox[3][3][3] ) 02768 { 02769 // Stupid flood fill 02770 int tot = 0 ; 02771 int queue[27][3] ; 02772 int vis[3][3][3] ; 02773 int head = 0, tail = 1 ; 02774 int i, j, k ; 02775 for ( i = 0 ; i < 3 ; i ++ ) 02776 for ( j = 0 ; j < 3 ; j ++ ) 02777 for ( k = 0 ; k < 3 ; k ++ ) 02778 { 02779 vis[i][j][k] = 0 ; 02780 if ( vox[i][j][k] ) 02781 { 02782 if ( tot == 0 ) 02783 { 02784 queue[0][0] = i ; 02785 queue[0][1] = j ; 02786 queue[0][2] = k ; 02787 vis[i][j][k] = 1 ; 02788 } 02789 tot ++ ; 02790 } 02791 } 02792 if ( tot == 0 ) 02793 { 02794 return 0 ; 02795 } 02796 02797 int ct = 1 ; 02798 while ( head != tail ) 02799 { 02800 int x = queue[head][0] ; 02801 int y = queue[head][1] ; 02802 int z = queue[head][2] ; 02803 head ++ ; 02804 02805 for ( i = -1 ; i < 2 ; i ++ ) 02806 for ( j = -1 ; j < 2 ; j ++ ) 02807 for ( k = -1 ; k < 2 ; k ++ ) 02808 { 02809 int nx = x + i ; 02810 int ny = y + j ; 02811 int nz = z + k ; 02812 if ( nx >=0 && nx < 3 && ny >=0 && ny < 3 && nz >=0 && nz < 3 ) 02813 { 02814 if ( vox[nx][ny][nz] && ! vis[nx][ny][nz] ) 02815 { 02816 queue[tail][0] = nx ; 02817 queue[tail][1] = ny ; 02818 queue[tail][2] = nz ; 02819 tail ++ ; 02820 vis[nx][ny][nz] = 1 ; 02821 ct ++ ; 02822 } 02823 } 02824 } 02825 } 02826 02827 if ( ct == tot ) 02828 { 02829 return 1 ; 02830 } 02831 else 02832 { 02833 return 2 ; 02834 } 02835 02836 } 02837 02838 int Volume::countExt( double vox[3][3][3] ) 02839 { 02840 int tvox[3][3][3] ; 02841 02842 for ( int i = 0 ; i < 3 ; i ++ ) 02843 for ( int j = 0 ; j < 3 ; j ++ ) 02844 for ( int k = 0 ; k < 3 ; k ++ ) 02845 { 02846 if ( vox[i][j][k] < 0 ) 02847 { 02848 tvox[i][j][k] = 1 ; 02849 } 02850 else 02851 { 02852 tvox[i][j][k] = 0 ; 02853 } 02854 } 02855 02856 return components26( tvox ) ; 02857 } 02858 02859 int Volume::countInt( double vox[3][3][3] ) 02860 { 02861 int i, j, k ; 02862 int tvox[3][3][3] ; 02863 02864 for ( i = 0 ; i < 3 ; i ++ ) 02865 for ( j = 0 ; j < 3 ; j ++ ) 02866 for ( k = 0 ; k < 3 ; k ++ ) 02867 { 02868 tvox[i][j][k] = 0 ; 02869 } 02870 02871 for ( i = 0 ; i < 6 ; i ++ ) 02872 { 02873 int nx = 1 + neighbor6[i][0] ; 02874 int ny = 1 + neighbor6[i][1] ; 02875 int nz = 1 + neighbor6[i][2] ; 02876 if ( vox[ nx ][ ny ][ nz ] >= 0 ) 02877 { 02878 tvox[ nx ][ ny ][ nz ] = 1 ; 02879 for ( j = 0 ; j < 4 ; j ++ ) 02880 { 02881 int nnx = neighbor64[i][j][0] + nx ; 02882 int nny = neighbor64[i][j][1] + ny ; 02883 int nnz = neighbor64[i][j][2] + nz ; 02884 if ( vox[ nnx ][ nny ][ nnz ] >= 0 ) 02885 { 02886 tvox[ nnx ][ nny ][ nnz ] = 1 ; 02887 } 02888 } 02889 } 02890 } 02891 02892 return components6( tvox ) ; 02893 } 02894 02895 int Volume::countIntEuler( int ox, int oy, int oz ) 02896 { 02897 int nv = 0 , ne = 0, nc = 0 ; 02898 02899 int i, j, k ; 02900 int tvox[3][3][3] ; 02901 double vox[3][3][3] ; 02902 02903 for ( i = 0 ; i < 3 ; i ++ ) 02904 for ( j = 0 ; j < 3 ; j ++ ) 02905 for ( k = 0 ; k < 3 ; k ++ ) 02906 { 02907 vox[i][j][k] = getDataAt( ox - 1 + i, oy - 1 + j, oz - 1 + k ) ; 02908 tvox[i][j][k] = 0 ; 02909 } 02910 02911 for ( i = 0 ; i < 6 ; i ++ ) 02912 { 02913 int nx = 1 + neighbor6[i][0] ; 02914 int ny = 1 + neighbor6[i][1] ; 02915 int nz = 1 + neighbor6[i][2] ; 02916 if ( vox[nx][ny][nz] >= 0 ) 02917 { 02918 tvox[ nx ][ ny ][ nz ] = 1 ; 02919 02920 nv ++ ; 02921 02922 for ( j = 0 ; j < 4 ; j ++ ) 02923 { 02924 int nnx = neighbor64[i][j][0] + nx ; 02925 int nny = neighbor64[i][j][1] + ny ; 02926 int nnz = neighbor64[i][j][2] + nz ; 02927 if ( vox[nnx][nny][nnz] >= 0 ) 02928 { 02929 if ( tvox[nnx][nny][nnz] == 0 ) 02930 { 02931 tvox[nnx][nny][nnz] = 1 ; 02932 nv ++ ; 02933 } 02934 02935 ne ++ ; 02936 } 02937 } 02938 } 02939 } 02940 02941 nc = components6( tvox ) ; 02942 02943 return ( nc - ( nv - ne ) ) ; 02944 } 02945 02946 //void Volume::erodeNoTopo( float thr, int wid ) 02947 //{ 02949 //int i, j, k ; 02951 //#ifdef VERBOSE 02952 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ; 02953 //#endif 02954 //threshold( thr, -MAX_ERODE, 0 ) ; 02955 02957 //#ifdef VERBOSE 02958 //printf("Initializing queue...\n") ; 02959 //#endif 02960 //GridQueue* queue = new GridQueue( ) ; 02961 02962 //for ( i = 0 ; i < getSizeX() ; i ++ ) 02963 //for ( j = 0 ; j < getSizeY() ; j ++ ) 02964 //for ( k = 0 ; k < getSizeZ() ; k ++ ) 02965 //{ 02966 //if ( getDataAt( i, j, k ) >= 0 ) 02967 //{ 02968 //for ( int m = 0 ; m < 6 ; m ++ ) 02969 //{ 02970 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 ) 02971 //{ 02972 //setDataAt( i, j, k, 1 ) ; 02973 //queue->pushQueue( i, j, k ) ; 02974 //break ; 02975 //} 02976 //} 02977 //} 02978 //} 02979 //#ifdef VERBOSE 02980 //printf("Total %d nodes\n", queue->getNumElements() ) ; 02981 02983 //printf("Start erosion to %d...\n", wid) ; 02984 //#endif 02985 //double val = 0; 02986 //int ox, oy, oz ; 02987 //int curwid = 0 ; 02988 //int total = 0, ct = 0 ; 02989 //while ( 1 ) 02990 //{ 02991 //if ( ct == total ) 02992 //{ 02993 //#ifdef VERBOSE 02994 //printf("Layer %d has %d nodes.\n", (int) curwid, total ) ; 02995 //#endif 02996 //curwid ++ ; 02997 //ct = 0 ; 02998 //total = queue->getNumElements() ; 02999 //if ( total == 0 ) 03000 //{ 03001 //break ; 03002 //} 03003 //} 03004 03005 //queue->popQueue(ox, oy, oz) ; 03006 //val = getDataAt( ox, oy, oz ) ; 03007 //if ( val > wid ) 03008 //{ 03009 //break ; 03010 //} 03011 //ct ++ ; 03012 03013 //setDataAt( ox, oy, oz, -val ) ; 03014 03015 03016 //for ( int m = 0 ; m < 6 ; m ++ ) 03017 //{ 03018 //int nx = ox + neighbor6[m][0] ; 03019 //int ny = oy + neighbor6[m][1] ; 03020 //int nz = oz + neighbor6[m][2] ; 03021 //if ( getDataAt( nx, ny, nz ) == 0 ) 03022 //{ 03023 //setDataAt( nx, ny, nz, val + 1 ) ; 03024 //queue->pushQueue( nx, ny, nz ) ; 03025 //} 03026 //} 03027 //} 03028 03030 //#ifdef VERBOSE 03031 //printf("Thresholding the volume to 0/1...\n") ; 03032 //#endif 03033 //threshold( 0, 0, 1 ) ; 03034 03035 //} 03036 03037 //void Volume::erodeTopo( float thr, int wid ) 03038 //{ 03040 //int i, j, k ; 03042 //#ifdef VERBOSE 03043 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ; 03044 //#endif 03045 //threshold( thr, -MAX_ERODE, 0 ) ; 03046 03048 //#ifdef VERBOSE 03049 //printf("Initializing queue...\n") ; 03050 //#endif 03051 //GridQueue* queue = new GridQueue( ) ; 03052 03053 //for ( i = 0 ; i < getSizeX() ; i ++ ) 03054 //for ( j = 0 ; j < getSizeY() ; j ++ ) 03055 //for ( k = 0 ; k < getSizeZ() ; k ++ ) 03056 //{ 03057 //if ( getDataAt( i, j, k ) >= 0 ) 03058 //{ 03059 //for ( int m = 0 ; m < 6 ; m ++ ) 03060 //{ 03061 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 ) 03062 //{ 03063 //setDataAt( i, j, k, 1 ) ; 03064 //queue->pushQueue( i, j, k ) ; 03065 //break ; 03066 //} 03067 //} 03068 //} 03069 //} 03070 //#ifdef VERBOSE 03071 //printf("Total %d nodes\n", queue->getNumElements() ) ; 03072 03074 //printf("Start erosion to %d...\n", wid) ; 03075 //#endif 03076 03077 //double val = 0; 03078 //int ox, oy, oz ; 03079 //int curwid = 0 ; 03080 //int total = 0, ct = 0 ; 03081 //while ( 1 ) 03082 //{ 03083 //if ( ct == total ) 03084 //{ 03085 //#ifdef VERBOSE 03086 //printf("Layer %d has %d nodes.\n", (int) curwid, total ) ; 03087 //#endif 03088 //curwid ++ ; 03089 //ct = 0 ; 03090 //total = queue->getNumElements() ; 03091 //if ( total == 0 ) 03092 //{ 03093 //break ; 03094 //} 03095 //} 03096 03097 //queue->popQueue(ox, oy, oz) ; 03098 //val = getDataAt( ox, oy, oz ) ; 03099 //if ( val > wid ) 03100 //{ 03101 //break ; 03102 //} 03103 //ct ++ ; 03104 03105 //if ( isSimple( ox, oy, oz ) ) 03106 //{ 03108 //setDataAt( ox, oy, oz, -val ) ; 03109 //} 03110 //else 03111 //{ 03113 //setDataAt( ox, oy, oz, val + 1 ) ; 03114 //queue->pushQueue( ox, oy, oz ) ; 03115 //} 03116 03117 03118 //for ( int m = 0 ; m < 6 ; m ++ ) 03119 //{ 03120 //int nx = ox + neighbor6[m][0] ; 03121 //int ny = oy + neighbor6[m][1] ; 03122 //int nz = oz + neighbor6[m][2] ; 03123 //if ( getDataAt( nx, ny, nz ) == 0 ) 03124 //{ 03125 //setDataAt( nx, ny, nz, val + 1 ) ; 03126 //queue->pushQueue( nx, ny, nz ) ; 03127 //} 03128 //} 03129 //} 03130 03132 //#ifdef VERBOSE 03133 //printf("Thresholding the volume to 0/1...\n") ; 03134 //#endif 03135 //threshold( 0, 0, 1 ) ; 03136 03137 //} 03138 03139 //void Volume::erode2( float thr, int wid ) 03140 //{ 03141 //int i, j, k ; 03143 //#ifdef VERBOSE 03144 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ; 03145 //#endif 03146 //threshold( thr, -MAX_ERODE, 0 ) ; 03147 03149 //#ifdef VERBOSE 03150 //printf("Initializing queue...\n") ; 03151 //#endif 03152 //GridQueue2* queue = new GridQueue2( ) ; 03153 //GridQueue2* queue2 = new GridQueue2( ) ; 03154 03155 //for ( i = 0 ; i < getSizeX() ; i ++ ) 03156 //for ( j = 0 ; j < getSizeY() ; j ++ ) 03157 //for ( k = 0 ; k < getSizeZ() ; k ++ ) 03158 //{ 03159 //if ( getDataAt( i, j, k ) >= 0 ) 03160 //{ 03161 //for ( int m = 0 ; m < 6 ; m ++ ) 03162 //{ 03163 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 ) 03164 //{ 03165 //setDataAt( i, j, k, 1 ) ; 03166 //queue->prepend( i, j, k ) ; 03167 //break ; 03168 //} 03169 //} 03170 //} 03171 //} 03172 //#ifdef VERBOSE 03173 //printf("Total %d nodes\n", queue->getNumElements() ) ; 03174 03177 //printf("Start erosion to %d...\n", wid) ; 03178 //#endif 03179 //gridQueueEle* ele ; 03180 //int ox, oy, oz ; 03181 03182 //for ( int curwid = 1 ; curwid <= wid ; curwid ++ ) 03183 //{ 03184 //int numComplex = 0, numSimple = 0 ; 03185 //#ifdef VERBOSE 03186 //printf("Processing %d nodes in layer %d\n", queue->getNumElements(), curwid) ; 03187 //#endif 03188 03190 //while ( ( ele = queue->getNext() ) != NULL ) 03191 //{ 03192 //for ( int m = 0 ; m < 6 ; m ++ ) 03193 //{ 03194 //int nx = ele->x + neighbor6[m][0] ; 03195 //int ny = ele->y + neighbor6[m][1] ; 03196 //int nz = ele->z + neighbor6[m][2] ; 03197 //if ( getDataAt( nx, ny, nz ) == 0 ) 03198 //{ 03199 //setDataAt( nx, ny, nz, curwid + 1 ) ; 03200 //queue2->prepend( nx, ny, nz ) ; 03201 //} 03202 //} 03203 03204 //} 03205 //*/ 03206 03208 //int seed[3] = {-1,-1,-1}; 03209 //queue->reset() ; 03210 //while ( queue->getNumElements() > 0 ) 03211 //{ 03212 //if ( seed[0] < 0 ) printf("After initial scoring...\n"); 03213 //queue->reset() ; 03214 //ele = queue->getNext() ; 03215 03217 //while ( ele != NULL ) 03218 //{ 03219 //ox = ele->x ; 03220 //oy = ele->y ; 03221 //oz = ele->z ; 03222 03224 //if ( seed[0] < 0 || 03225 //( ox < seed[0] + 2 && ox > seed[0] - 2 && 03226 //oy < seed[1] + 2 && oy > seed[1] - 2 && 03227 //oz < seed[2] + 2 && oz > seed[2] - 2 ) ) 03228 //{ 03229 //if ( ! isSimple( ox, oy, oz ) ) 03230 //{ 03232 //setDataAt( ox, oy, oz, curwid + 1 ) ; 03233 //queue2->prepend( ox, oy, oz ) ; 03234 //ele = queue->remove() ; 03235 03236 //numComplex ++ ; 03237 //} 03238 //else 03239 //{ 03240 //ele = queue->getNext() ; 03241 //} 03242 //} 03243 //else 03244 //{ 03245 //ele = queue->getNext() ; 03246 //} 03247 //} 03248 03250 //queue->reset() ; 03251 //ele = queue->getNext() ; 03252 //int preScore = -1; 03253 //while ( ele != NULL ) 03254 //{ 03255 //ox = ele->x ; 03256 //oy = ele->y ; 03257 //oz = ele->z ; 03258 03260 //if ( seed[0] < 0 || 03261 //( ox < seed[0] + 3 && ox > seed[0] - 3 && 03262 //oy < seed[1] + 3 && oy > seed[1] - 3 && 03263 //oz < seed[2] + 3 && oz > seed[2] - 3 ) ) 03264 //{ 03265 //ele->score = getNumPotComplex( ox, oy, oz ) ; 03266 //} 03267 03268 03269 //if ( ele->score < preScore ) 03270 //{ 03272 //ele = queue->swap() ; 03273 //} 03274 //else 03275 //{ 03276 //preScore = ele->score ; 03277 //} 03278 03280 //if ( ele->next == NULL ) 03281 //{ 03282 //ox = ele->x ; 03283 //oy = ele->y ; 03284 //oz = ele->z ; 03285 //setDataAt( ox, oy, oz, -1 ) ; 03287 //seed[0] = ox ; 03288 //seed[1] = oy ; 03289 //seed[2] = oz ; 03290 //queue->remove() ; 03292 03293 03294 //for ( int m = 0 ; m < 6 ; m ++ ) 03295 //{ 03296 //int nx = ox + neighbor6[m][0] ; 03297 //int ny = oy + neighbor6[m][1] ; 03298 //int nz = oz + neighbor6[m][2] ; 03299 //if ( getDataAt( nx, ny, nz ) == 0 ) 03300 //{ 03301 //setDataAt( nx, ny, nz, curwid + 1 ) ; 03302 //queue2->prepend( nx, ny, nz ) ; 03303 //} 03304 //} 03305 03306 03307 //numSimple ++ ; 03308 //ele = NULL ; 03309 //} 03310 //else 03311 //{ 03312 //ele = queue->getNext() ; 03313 //} 03314 //} 03315 03316 //} 03317 03318 //delete queue ; 03319 //queue = queue2 ; 03320 //queue2 = new GridQueue2() ; 03321 //#ifdef VERBOSE 03322 //printf("%d complex, %d simple\n", numComplex, numSimple) ; 03323 //#endif 03324 03325 //if ( numSimple == 0 ) 03326 //{ 03327 //break ; 03328 //} 03329 //} 03330 03332 //#ifdef VERBOSE 03333 //printf("Thresholding the volume to 0/1...\n") ; 03334 //#endif 03335 //threshold( 0, 0, 1 ) ; 03336 03337 //} 03338 03339 //void Volume::erodeShapeTopo( float thr, int wid ) 03340 //{ 03342 //int i, j, k ; 03344 //#ifdef VERBOSE 03345 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ; 03346 //#endif 03347 //threshold( thr, -MAX_ERODE, 0 ) ; 03348 03350 //#ifdef VERBOSE 03351 //printf("Initializing queue...\n") ; 03352 //#endif 03353 //GridQueue2* queue2 = new GridQueue2( ) ; 03354 //GridQueue2* queue3 = new GridQueue2( ) ; 03355 //PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN ); 03356 03357 //for ( i = 0 ; i < getSizeX() ; i ++ ) 03358 //for ( j = 0 ; j < getSizeY() ; j ++ ) 03359 //for ( k = 0 ; k < getSizeZ() ; k ++ ) 03360 //{ 03361 //if ( getDataAt( i, j, k ) >= 0 ) 03362 //{ 03363 //for ( int m = 0 ; m < 6 ; m ++ ) 03364 //{ 03365 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 ) 03366 //{ 03367 //setDataAt( i, j, k, 1 ) ; 03368 //queue2->prepend( i, j, k ) ; 03369 //break ; 03370 //} 03371 //} 03372 //} 03373 //} 03374 //#ifdef VERBOSE 03375 //printf("Total %d nodes\n", queue2->getNumElements() ) ; 03376 03377 03380 //printf("Start erosion to %d...\n", wid) ; 03381 //#endif 03382 //gridQueueEle* ele ; 03383 //gridPoint* gp ; 03384 //int ox, oy, oz ; 03385 //int score ; 03386 //Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ; 03387 //for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ ) 03388 //{ 03389 //scrvol->setDataAt( i, -1 ) ; 03390 //} 03391 03392 //for ( int curwid = 1 ; curwid <= wid ; curwid ++ ) 03393 //{ 03397 03398 //int numComplex = 0, numSimple = 0 ; 03399 //#ifdef VERBOSE 03400 //printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ; 03401 //#endif 03402 03406 //queue2->reset() ; 03407 //ele = queue2->getNext() ; 03408 //while ( ele != NULL ) 03409 //{ 03410 //ox = ele->x ; 03411 //oy = ele->y ; 03412 //oz = ele->z ; 03413 03415 //if ( ! isSimple( ox, oy, oz ) ) 03416 //{ 03418 //setDataAt( ox, oy, oz, curwid + 1 ) ; 03419 //queue3->prepend( ox, oy, oz ) ; 03420 //ele = queue2->remove() ; 03421 03422 //numComplex ++ ; 03423 //} 03424 //else 03425 //{ 03426 //ele = queue2->getNext() ; 03427 //} 03428 //} 03429 03433 //queue2->reset() ; 03434 //ele = queue2->getNext() ; 03435 //while ( ele != NULL ) 03436 //{ 03437 //ox = ele->x ; 03438 //oy = ele->y ; 03439 //oz = ele->z ; 03440 03442 //score = getNumPotComplex( ox, oy, oz ) ; 03443 //scrvol->setDataAt( ox, oy, oz, score ) ; 03444 03446 //gp = new gridPoint ; 03447 //gp->x = ox ; 03448 //gp->y = oy ; 03449 //gp->z = oz ; 03450 //queue->add( gp, -score ) ; 03451 03452 //ele = queue2->remove() ; 03453 //} 03454 03457 //delete queue2 ; 03458 //queue2 = queue3 ; 03459 //queue3 = new GridQueue2( ) ; 03460 03462 //while ( ! queue->isEmpty() ) 03463 //{ 03465 //queue->remove( gp, score ) ; 03466 //ox = gp->x ; 03467 //oy = gp->y ; 03468 //oz = gp->z ; 03469 //delete gp ; 03470 //score = -score ; 03471 03475 //if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score ) 03476 //{ 03477 //continue ; 03478 //} 03479 03481 //setDataAt( ox, oy, oz, -1 ) ; 03482 //numSimple ++ ; 03484 03486 //for ( int m = 0 ; m < 6 ; m ++ ) 03487 //{ 03488 //int nx = ox + neighbor6[m][0] ; 03489 //int ny = oy + neighbor6[m][1] ; 03490 //int nz = oz + neighbor6[m][2] ; 03491 //if ( getDataAt( nx, ny, nz ) == 0 ) 03492 //{ 03493 //setDataAt( nx, ny, nz, curwid + 1 ) ; 03494 //queue2->prepend( nx, ny, nz ) ; 03495 //} 03496 //} 03497 03500 //for ( i = -1 ; i < 2 ; i ++ ) 03501 //for ( j = -1 ; j < 2 ; j ++ ) 03502 //for ( k = -1 ; k < 2 ; k ++ ) 03503 //{ 03504 //int nx = ox + i ; 03505 //int ny = oy + j ; 03506 //int nz = oz + k ; 03507 03509 //if ( getDataAt( nx, ny, nz ) == curwid && ! isSimple( nx, ny, nz ) ) 03510 //{ 03512 //setDataAt( nx, ny, nz, curwid + 1 ) ; 03513 //queue2->prepend( nx, ny, nz ) ; 03514 //numComplex ++ ; 03515 //} 03516 //} 03517 03521 //for ( i = -2 ; i < 3 ;i ++ ) 03522 //for ( j = -2 ; j < 3 ; j ++ ) 03523 //for ( k = -2 ; k < 3 ; k ++ ) 03524 //{ 03525 //int nx = ox + i ; 03526 //int ny = oy + j ; 03527 //int nz = oz + k ; 03528 03529 //if ( getDataAt( nx, ny, nz ) == curwid ) 03530 //{ 03532 //score = getNumPotComplex( nx, ny, nz ) ; 03533 03534 //if ( score != (int) scrvol->getDataAt( nx, ny, nz ) ) 03535 //{ 03537 //scrvol->setDataAt( nx, ny, nz, score ) ; 03539 //gp = new gridPoint ; 03540 //gp->x = nx ; 03541 //gp->y = ny ; 03542 //gp->z = nz ; 03543 //queue->add( gp, -score ) ; 03544 //} 03545 //} 03546 //} 03547 //*/ 03548 03549 //} 03550 03551 //#ifdef VERBOSE 03552 //printf("%d complex, %d simple\n", numComplex, numSimple) ; 03553 //#endif 03554 //if ( numSimple == 0 ) 03555 //{ 03556 //break ; 03557 //} 03558 //} 03559 03561 //#ifdef VERBOSE 03562 //printf("Thresholding the volume to 0/1...\n") ; 03563 //#endif 03564 //threshold( 0, 0, 1 ) ; 03565 //delete queue; 03566 03567 //} 03568 03569 03570 //void Volume::erodeAtom( float thr, int wid, Volume* avol ) 03571 //{ 03573 //int i, j, k ; 03575 //#ifdef VERBOSE 03576 //printf("Thresholding the volume to -MAX_ERODE/0...\n") ; 03577 //#endif 03578 //threshold( thr, -MAX_ERODE, 0 ) ; 03579 03581 //#ifdef VERBOSE 03582 //printf("Initializing queue...\n") ; 03583 //#endif 03584 //GridQueue2* queue2 = new GridQueue2( ) ; 03585 //GridQueue2* queue3 = new GridQueue2( ) ; 03586 //PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN ); 03587 03588 //for ( i = 0 ; i < getSizeX() ; i ++ ) 03589 //for ( j = 0 ; j < getSizeY() ; j ++ ) 03590 //for ( k = 0 ; k < getSizeZ() ; k ++ ) 03591 //{ 03592 //if ( getDataAt( i, j, k ) >= 0 ) 03593 //{ 03594 //if ( avol->getDataAt(i,j,k) > 0 ) 03595 //{ 03596 //setDataAt( i, j, k, MAX_ERODE ) ; 03597 //} 03598 //else 03599 //{ 03600 //for ( int m = 0 ; m < 6 ; m ++ ) 03601 //{ 03602 //if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 ) 03603 //{ 03604 //setDataAt( i, j, k, 1 ) ; 03605 //queue2->prepend( i, j, k ) ; 03606 //break ; 03607 //} 03608 //} 03609 //} 03610 //} 03611 //} 03612 //#ifdef VERBOSE 03613 //printf("Total %d nodes\n", queue2->getNumElements() ) ; 03614 03615 03618 //printf("Start erosion to %d...\n", wid) ; 03619 //#endif 03620 03621 //gridQueueEle* ele ; 03622 //gridPoint* gp ; 03623 //int ox, oy, oz ; 03624 //int score ; 03625 //Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ; 03626 //for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ ) 03627 //{ 03628 //scrvol->setDataAt( i, -1 ) ; 03629 //} 03630 03631 //for ( int curwid = 1 ; curwid <= wid ; curwid ++ ) 03632 //{ 03636 03637 //int numComplex = 0, numSimple = 0 ; 03638 //#ifdef VERBOSE 03639 //printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ; 03640 //#endif 03641 03645 //queue2->reset() ; 03646 //ele = queue2->getNext() ; 03647 //while ( ele != NULL ) 03648 //{ 03649 //ox = ele->x ; 03650 //oy = ele->y ; 03651 //oz = ele->z ; 03652 03654 //if ( ! isSimple( ox, oy, oz ) ) 03655 //{ 03657 //setDataAt( ox, oy, oz, curwid + 1 ) ; 03658 //queue3->prepend( ox, oy, oz ) ; 03659 //ele = queue2->remove() ; 03660 03661 //numComplex ++ ; 03662 //} 03663 //else 03664 //{ 03665 //ele = queue2->getNext() ; 03666 //} 03667 //} 03668 03672 //queue2->reset() ; 03673 //ele = queue2->getNext() ; 03674 //while ( ele != NULL ) 03675 //{ 03676 //ox = ele->x ; 03677 //oy = ele->y ; 03678 //oz = ele->z ; 03679 03681 //score = getNumPotComplex( ox, oy, oz ) ; 03682 //scrvol->setDataAt( ox, oy, oz, score ) ; 03683 03685 //gp = new gridPoint ; 03686 //gp->x = ox ; 03687 //gp->y = oy ; 03688 //gp->z = oz ; 03689 //queue->add( gp, -score ) ; 03690 03691 //ele = queue2->remove() ; 03692 //} 03693 03696 //delete queue2 ; 03697 //queue2 = queue3 ; 03698 //queue3 = new GridQueue2( ) ; 03699 03701 //while ( ! queue->isEmpty() ) 03702 //{ 03704 //queue->remove( gp, score ) ; 03705 //ox = gp->x ; 03706 //oy = gp->y ; 03707 //oz = gp->z ; 03708 //delete gp ; 03709 //score = -score ; 03710 03714 //if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score ) 03715 //{ 03716 //continue ; 03717 //} 03718 03720 //setDataAt( ox, oy, oz, -1 ) ; 03721 //numSimple ++ ; 03723 03725 //for ( int m = 0 ; m < 6 ; m ++ ) 03726 //{ 03727 //int nx = ox + neighbor6[m][0] ; 03728 //int ny = oy + neighbor6[m][1] ; 03729 //int nz = oz + neighbor6[m][2] ; 03730 //if ( getDataAt( nx, ny, nz ) == 0 ) 03731 //{ 03732 //setDataAt( nx, ny, nz, curwid + 1 ) ; 03733 //queue2->prepend( nx, ny, nz ) ; 03734 //} 03735 //} 03736 03739 //for ( i = -1 ; i < 2 ; i ++ ) 03740 //for ( j = -1 ; j < 2 ; j ++ ) 03741 //for ( k = -1 ; k < 2 ; k ++ ) 03742 //{ 03743 //int nx = ox + i ; 03744 //int ny = oy + j ; 03745 //int nz = oz + k ; 03746 03748 //if ( getDataAt( nx, ny, nz ) == curwid && ! isSimple( nx, ny, nz ) ) 03749 03750 //{ 03752 //setDataAt( nx, ny, nz, curwid + 1 ) ; 03753 //queue2->prepend( nx, ny, nz ) ; 03754 //numComplex ++ ; 03755 //} 03756 //} 03757 03761 //for ( i = -2 ; i < 3 ;i ++ ) 03762 //for ( j = -2 ; j < 3 ; j ++ ) 03763 //for ( k = -2 ; k < 3 ; k ++ ) 03764 //{ 03765 //int nx = ox + i ; 03766 //int ny = oy + j ; 03767 //int nz = oz + k ; 03768 03769 //if ( getDataAt( nx, ny, nz ) == curwid ) 03770 //{ 03772 //score = getNumPotComplex( nx, ny, nz ) ; 03773 03774 //if ( score != (int) scrvol->getDataAt( nx, ny, nz ) ) 03775 //{ 03777 //scrvol->setDataAt( nx, ny, nz, score ) ; 03779 //gp = new gridPoint ; 03780 //gp->x = nx ; 03781 //gp->y = ny ; 03782 //gp->z = nz ; 03783 //queue->add( gp, -score ) ; 03784 //} 03785 //} 03786 //} 03787 //*/ 03788 03789 //} 03790 //#ifdef VERBOSE 03791 //printf("%d complex, %d simple\n", numComplex, numSimple) ; 03792 //#endif 03793 03794 //if ( numSimple == 0 ) 03795 //{ 03796 //break ; 03797 //} 03798 //} 03799 03801 //#ifdef VERBOSE 03802 //printf("Thresholding the volume to 0/1...\n") ; 03803 //#endif 03804 //threshold( 0, 0, 1 ) ; 03805 //delete queue; 03806 03807 //} 03808 03810 //* Assuming the current volume has already been thresholded to 0/1 03811 //*/ 03812 void Volume::curveSkeleton( Volume* grayvol, float lowthr, float highthr, Volume* svol ) 03813 { 03814 int i, j, k ; 03815 // First, threshold the volume 03816 #ifdef VERBOSE 03817 printf("Thresholding the volume to -1/0...\n") ; 03818 #endif 03819 threshold( 0.5f, -1, 0 ) ; 03820 03821 // Next, apply convergent erosion 03822 // by preserving: complex nodes, curve end-points, and sheet points 03823 03824 // Next, initialize the linked queue 03825 #ifdef VERBOSE 03826 printf("Initializing queue...\n") ; 03827 #endif 03828 GridQueue2* queue2 = new GridQueue2( ) ; 03829 GridQueue2* queue3 = new GridQueue2( ) ; 03830 GridQueue2* queue4 = new GridQueue2( ) ; 03831 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN ); 03832 03833 for ( i = 0 ; i < getSizeX() ; i ++ ) 03834 for ( j = 0 ; j < getSizeY() ; j ++ ) 03835 for ( k = 0 ; k < getSizeZ() ; k ++ ) 03836 { 03837 if ( getDataAt( i, j, k ) >= 0 ) 03838 { 03839 float v = (float)grayvol->getDataAt(i,j,k) ; 03840 if ( v <= lowthr || v > highthr || svol->getDataAt(i,j,k) > 0 ) 03841 { 03842 setDataAt( i, j, k, MAX_ERODE ) ; 03843 } 03844 else 03845 { 03846 for ( int m = 0 ; m < 6 ; m ++ ) 03847 { 03848 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 ) 03849 { 03850 // setDataAt( i, j, k, 1 ) ; 03851 queue2->prepend( i, j, k ) ; 03852 break ; 03853 } 03854 } 03855 } 03856 } 03857 } 03858 int wid = MAX_ERODE ; 03859 #ifdef VERBOSE 03860 printf("Total %d nodes\n", queue2->getNumElements() ) ; 03861 printf("Start erosion to %d...\n", wid) ; 03862 #endif 03863 03864 03865 // Perform erosion 03866 gridQueueEle* ele ; 03867 gridPoint* gp ; 03868 int ox, oy, oz ; 03869 int score ; 03870 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ; 03871 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ ) 03872 { 03873 scrvol->setDataAt( i, -1 ) ; 03874 } 03875 03876 #ifdef NOISE_DIS_HELIX 03877 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ; 03878 #endif 03879 03880 for ( int curwid = 1 ; curwid <= wid ; curwid ++ ) 03881 { 03882 // At the start of each iteration, 03883 // queue2 holds all the nodes for this layer 03884 // queue3 and queue are empty 03885 03886 int numComplex = 0, numSimple = 0 ; 03887 #ifdef VERBOSE 03888 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ; 03889 #endif 03890 03891 /* 03892 We first need to assign curwid + 1 to every node in this layer 03893 */ 03894 queue2->reset() ; 03895 ele = queue2->getNext() ; 03896 while ( ele != NULL ) 03897 { 03898 ox = ele->x ; 03899 oy = ele->y ; 03900 oz = ele->z ; 03901 03902 if ( getDataAt(ox,oy,oz) == curwid ) 03903 { 03904 ele = queue2->remove() ; 03905 } 03906 else 03907 { 03908 setDataAt(ox,oy,oz, curwid) ; 03909 ele = queue2->getNext() ; 03910 } 03911 } 03912 queue4->reset() ; 03913 ele = queue4->getNext() ; 03914 while ( ele != NULL ) 03915 { 03916 ox = ele->x ; 03917 oy = ele->y ; 03918 oz = ele->z ; 03919 03920 queue2->prepend(ox,oy,oz) ; 03921 ele = queue4->remove() ; 03922 } 03923 03924 // Now queue2 holds all the nodes for this layer 03925 03926 #ifdef NOISE_DIS_HELIX 03927 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */ 03928 queue2->reset() ; 03929 03930 // First run 03931 int flag = 0 ; 03932 while ( ( ele = queue2->getNext() ) != NULL ) 03933 { 03934 ox = ele->x ; 03935 oy = ele->y ; 03936 oz = ele->z ; 03937 if ( NOISE_DIS_HELIX <= 1 ) 03938 { 03939 noisevol->setDataAt( ox, oy, oz, 0 ) ; 03940 } 03941 else 03942 { 03943 flag = 0 ; 03944 for ( int m = 0 ; m < 6 ; m ++ ) 03945 { 03946 int nx = ox + neighbor6[m][0] ; 03947 int ny = oy + neighbor6[m][1] ; 03948 int nz = oz + neighbor6[m][2] ; 03949 if ( getDataAt( nx, ny, nz ) == 0 ) 03950 { 03951 noisevol->setDataAt( ox, oy, oz, 1 ) ; 03952 flag = 1 ; 03953 break ; 03954 } 03955 } 03956 if ( ! flag ) 03957 { 03958 noisevol->setDataAt( ox, oy, oz, 0 ) ; 03959 } 03960 } 03961 } 03962 03963 int cur, visited ; 03964 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ ) 03965 { 03966 queue2->reset() ; 03967 int count = 0 ; 03968 visited = 0 ; 03969 03970 while ( ( ele = queue2->getNext() ) != NULL ) 03971 { 03972 ox = ele->x ; 03973 oy = ele->y ; 03974 oz = ele->z ; 03975 03976 if ( noisevol->getDataAt( ox, oy, oz ) == 1 ) 03977 { 03978 visited ++ ; 03979 continue ; 03980 } 03981 03982 flag = 0 ; 03983 for ( int m = 0 ; m < 6 ; m ++ ) 03984 { 03985 int nx = ox + neighbor6[m][0] ; 03986 int ny = oy + neighbor6[m][1] ; 03987 int nz = oz + neighbor6[m][2] ; 03988 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 ) 03989 { 03990 noisevol->setDataAt( ox, oy, oz, 1 ) ; 03991 visited ++ ; 03992 count ++ ; 03993 break ; 03994 } 03995 } 03996 } 03997 03998 if ( count == 0 ) 03999 { 04000 break ; 04001 } 04002 } 04003 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ; 04004 04005 04006 #endif 04007 /* Commented out for debugging 04008 04009 // First, 04010 // check for complex nodes in queue2 04011 // move them from queue2 to queue3 04012 queue2->reset() ; 04013 ele = queue2->getNext() ; 04014 while ( ele != NULL ) 04015 { 04016 ox = ele->x ; 04017 oy = ele->y ; 04018 oz = ele->z ; 04019 04020 // Check simple 04021 #ifndef NOISE_DIS_HELIX 04022 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) ) 04023 #else 04024 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) ) 04025 #endif 04026 { 04027 // Complex, set to next layer 04028 setDataAt( ox, oy, oz, curwid + 1 ) ; 04029 queue3->prepend( ox, oy, oz ) ; 04030 ele = queue2->remove() ; 04031 04032 numComplex ++ ; 04033 } 04034 else 04035 { 04036 ele = queue2->getNext() ; 04037 } 04038 } 04039 */ 04040 04041 // Next, 04042 // Compute score for each node left in queue2 04043 // move them into priority queue 04044 queue2->reset() ; 04045 ele = queue2->getNext() ; 04046 while ( ele != NULL ) 04047 { 04048 ox = ele->x ; 04049 oy = ele->y ; 04050 oz = ele->z ; 04051 04052 // Compute score 04053 score = getNumPotComplex2( ox, oy, oz ) ; 04054 scrvol->setDataAt( ox, oy, oz, score ) ; 04055 04056 // Push to queue 04057 gp = new gridPoint ; 04058 gp->x = ox ; 04059 gp->y = oy ; 04060 gp->z = oz ; 04061 // queue->add( gp, -score ) ; 04062 queue->add( gp, score ) ; 04063 04064 ele = queue2->remove() ; 04065 } 04066 04067 // Rename queue3 to be queue2, 04068 // Clear queue3 04069 // From now on, queue2 holds nodes of next level 04070 delete queue2 ; 04071 queue2 = queue3 ; 04072 queue3 = new GridQueue2( ) ; 04073 04074 // Next, start priority queue iteration 04075 while ( ! queue->isEmpty() ) 04076 { 04077 // Retrieve the node with the highest score 04078 queue->remove( gp, score ) ; 04079 ox = gp->x ; 04080 oy = gp->y ; 04081 oz = gp->z ; 04082 delete gp ; 04083 // score = -score ; 04084 04085 // Ignore the node 04086 // if it has been processed before 04087 // or it has an updated score 04088 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score ) 04089 { 04090 continue ; 04091 } 04092 04093 /* Commented out for debugging 04094 04095 // Remove this simple node 04096 setDataAt( ox, oy, oz, -1 ) ; 04097 numSimple ++ ; 04098 // printf("Highest score: %d\n", score) ; 04099 */ 04100 04101 /* Added for debugging */ 04102 // Check simple 04103 #ifndef NOISE_DIS_HELIX 04104 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) ) 04105 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) ) 04106 #else 04107 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) ) 04108 #endif 04109 { 04110 // Complex, set to next layer 04111 setDataAt( ox, oy, oz, curwid + 1 ) ; 04112 queue4->prepend( ox, oy, oz ) ; 04113 numComplex ++ ; 04114 } 04115 else 04116 { 04117 setDataAt( ox, oy, oz, -1 ) ; 04118 numSimple ++ ; 04119 04120 04121 /* Adding ends */ 04122 // Move its neighboring unvisited node to queue2 04123 for ( int m = 0 ; m < 6 ; m ++ ) 04124 { 04125 int nx = ox + neighbor6[m][0] ; 04126 int ny = oy + neighbor6[m][1] ; 04127 int nz = oz + neighbor6[m][2] ; 04128 if ( getDataAt( nx, ny, nz ) == 0 ) 04129 { 04130 // setDataAt( nx, ny, nz, curwid + 1 ) ; 04131 queue2->prepend( nx, ny, nz ) ; 04132 } 04133 } 04134 04135 } 04136 04137 04138 /* Commented out for debugging 04139 04140 // Find complex nodes in its 3x3 neighborhood 04141 // move them to queue2 04142 for ( i = -1 ; i < 2 ; i ++ ) 04143 for ( j = -1 ; j < 2 ; j ++ ) 04144 for ( k = -1 ; k < 2 ; k ++ ) 04145 { 04146 int nx = ox + i ; 04147 int ny = oy + j ; 04148 int nz = oz + k ; 04149 04150 // Check simple 04151 if ( getDataAt( nx, ny, nz ) == curwid && 04152 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) ) 04153 #ifndef NOISE_DIS_HELIX 04154 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) ) 04155 #else 04156 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) ) 04157 #endif 04158 04159 { 04160 // Complex, set to next layer 04161 setDataAt( nx, ny, nz, curwid + 1 ) ; 04162 queue2->prepend( nx, ny, nz ) ; 04163 numComplex ++ ; 04164 } 04165 } 04166 */ 04167 04168 // Update scores for nodes in its 5x5 neighborhood 04169 // insert them back into priority queue 04170 04171 for ( i = -2 ; i < 3 ;i ++ ) 04172 for ( j = -2 ; j < 3 ; j ++ ) 04173 for ( k = -2 ; k < 3 ; k ++ ) 04174 { 04175 int nx = ox + i ; 04176 int ny = oy + j ; 04177 int nz = oz + k ; 04178 04179 if ( getDataAt( nx, ny, nz ) == curwid ) 04180 { 04181 // Compute score 04182 score = getNumPotComplex2( nx, ny, nz ) ; 04183 04184 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) ) 04185 { 04186 // printf("Update\n") ; 04187 scrvol->setDataAt( nx, ny, nz, score ) ; 04188 // Push to queue 04189 gp = new gridPoint ; 04190 gp->x = nx ; 04191 gp->y = ny ; 04192 gp->z = nz ; 04193 // queue->add( gp, -score ) ; 04194 queue->add( gp, score ) ; 04195 } 04196 } 04197 } 04198 04199 04200 } 04201 04202 #ifdef VERBOSE 04203 printf("%d complex, %d simple\n", numComplex, numSimple) ; 04204 #endif 04205 04206 if ( numSimple == 0 ) 04207 { 04208 if ( queue2->getNumElements() > 0 ) 04209 { 04210 printf("*************************wierd here*************************\n"); 04211 } 04212 break ; 04213 } 04214 } 04215 04216 // Remove all internal voxels (contained in manifold surfaces) 04217 queue2->reset() ; 04218 queue4->reset() ; 04219 ele = queue4->getNext() ; 04220 while ( ele != NULL ) 04221 { 04222 ox = ele->x ; 04223 oy = ele->y ; 04224 oz = ele->z ; 04225 04226 if ( isPiercable(ox,oy,oz) == 1 ) // hasCompleteSheet( ox, oy, oz ) == 1 ) // 04227 { 04228 queue2->prepend(ox,oy,oz) ; 04229 // setDataAt( ox, oy, oz, -1 ) ; 04230 } 04231 ele = queue4->remove() ; 04232 } 04233 04234 for ( i = 0 ; i < getSizeX() ; i ++ ) 04235 for ( j = 0 ; j < getSizeY() ; j ++ ) 04236 for ( k = 0 ; k < getSizeZ() ; k ++ ) 04237 { 04238 if ( getDataAt( i, j, k ) == 0 && isPiercable(i,j,k) ) //hasCompleteSheet(i,j,k) == 1) // 04239 { 04240 queue2->prepend( i, j, k ) ; 04241 } 04242 } 04243 queue2->reset() ; 04244 ele = queue2->getNext() ; 04245 while ( ele != NULL ) 04246 { 04247 ox = ele->x ; 04248 oy = ele->y ; 04249 oz = ele->z ; 04250 setDataAt( ox, oy, oz, -1 ) ; 04251 ele = queue2->remove() ; 04252 } 04253 04254 04255 // Finally, clean up 04256 delete scrvol; 04257 delete queue; 04258 delete queue2; 04259 delete queue3; 04260 delete queue4; 04261 #ifdef VERBOSE 04262 printf("Thresholding the volume to 0/1...\n") ; 04263 #endif 04264 threshold( 0, 0, 1 ) ; 04265 } 04266 04267 // Compute curve skeleton 04268 void Volume::curveSkeleton( float thr, Volume* svol ) 04269 { 04270 int i, j, k ; 04271 // First, threshold the volume 04272 #ifdef VERBOSE 04273 printf("Thresholding the volume to -1/0...\n") ; 04274 #endif 04275 threshold( thr, -1, 0 ) ; 04276 04277 // Next, apply convergent erosion 04278 // by preserving: complex nodes, curve end-points, and sheet points 04279 04280 // Next, initialize the linked queue 04281 #ifdef VERBOSE 04282 printf("Initializing queue...\n") ; 04283 #endif 04284 GridQueue2* queue2 = new GridQueue2( ) ; 04285 GridQueue2* queue3 = new GridQueue2( ) ; 04286 GridQueue2* queue4 = new GridQueue2( ) ; 04287 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN ); 04288 04289 for ( i = 0 ; i < getSizeX() ; i ++ ) 04290 for ( j = 0 ; j < getSizeY() ; j ++ ) 04291 for ( k = 0 ; k < getSizeZ() ; k ++ ) 04292 { 04293 if ( getDataAt( i, j, k ) >= 0 ) 04294 { 04295 if ( svol->getDataAt(i,j,k) > 0 ) 04296 { 04297 setDataAt( i, j, k, MAX_ERODE ) ; 04298 } 04299 else 04300 { 04301 for ( int m = 0 ; m < 6 ; m ++ ) 04302 { 04303 if ( getDataAt( i + neighbor6[m][0], j + neighbor6[m][1], k + neighbor6[m][2] ) < 0 ) 04304 { 04305 // setDataAt( i, j, k, 1 ) ; 04306 queue2->prepend( i, j, k ) ; 04307 break ; 04308 } 04309 } 04310 } 04311 } 04312 } 04313 04314 int wid = MAX_ERODE ; 04315 #ifdef VERBOSE 04316 printf("Total %d nodes\n", queue2->getNumElements() ) ; 04317 printf("Start erosion to %d...\n", wid) ; 04318 #endif 04319 04320 04321 // Perform erosion 04322 gridQueueEle* ele ; 04323 gridPoint* gp ; 04324 int ox, oy, oz ; 04325 int score ; 04326 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ; 04327 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ ) 04328 { 04329 scrvol->setDataAt( i, -1 ) ; 04330 } 04331 04332 #ifdef NOISE_DIS_HELIX 04333 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ; 04334 #endif 04335 04336 for ( int curwid = 1 ; curwid <= wid ; curwid ++ ) 04337 { 04338 // At the start of each iteration, 04339 // queue2 holds all the nodes for this layer 04340 // queue3 and queue are empty 04341 04342 int numComplex = 0, numSimple = 0 ; 04343 #ifdef VERBOSE 04344 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ; 04345 #endif 04346 04347 /* 04348 We first need to assign curwid + 1 to every node in this layer 04349 */ 04350 queue2->reset() ; 04351 ele = queue2->getNext() ; 04352 while ( ele != NULL ) 04353 { 04354 ox = ele->x ; 04355 oy = ele->y ; 04356 oz = ele->z ; 04357 04358 if ( getDataAt(ox,oy,oz) == curwid ) 04359 { 04360 ele = queue2->remove() ; 04361 } 04362 else 04363 { 04364 setDataAt(ox,oy,oz, curwid) ; 04365 ele = queue2->getNext() ; 04366 } 04367 } 04368 queue4->reset() ; 04369 ele = queue4->getNext() ; 04370 while ( ele != NULL ) 04371 { 04372 ox = ele->x ; 04373 oy = ele->y ; 04374 oz = ele->z ; 04375 04376 queue2->prepend(ox,oy,oz) ; 04377 ele = queue4->remove() ; 04378 } 04379 04380 // Now queue2 holds all the nodes for this layer 04381 04382 #ifdef NOISE_DIS_HELIX 04383 /* Extra step: classify nodes in queue2 into noise and non-noise nodes */ 04384 queue2->reset() ; 04385 04386 // First run 04387 int flag = 0 ; 04388 while ( ( ele = queue2->getNext() ) != NULL ) 04389 { 04390 ox = ele->x ; 04391 oy = ele->y ; 04392 oz = ele->z ; 04393 if ( NOISE_DIS_HELIX <= 1 ) 04394 { 04395 noisevol->setDataAt( ox, oy, oz, 0 ) ; 04396 } 04397 else 04398 { 04399 flag = 0 ; 04400 for ( int m = 0 ; m < 6 ; m ++ ) 04401 { 04402 int nx = ox + neighbor6[m][0] ; 04403 int ny = oy + neighbor6[m][1] ; 04404 int nz = oz + neighbor6[m][2] ; 04405 if ( getDataAt( nx, ny, nz ) == 0 ) 04406 { 04407 noisevol->setDataAt( ox, oy, oz, 1 ) ; 04408 flag = 1 ; 04409 break ; 04410 } 04411 } 04412 if ( ! flag ) 04413 { 04414 noisevol->setDataAt( ox, oy, oz, 0 ) ; 04415 } 04416 } 04417 } 04418 04419 int cur, visited ; 04420 for ( cur = 1 ; cur < NOISE_DIS_HELIX ; cur ++ ) 04421 { 04422 queue2->reset() ; 04423 int count = 0 ; 04424 visited = 0 ; 04425 04426 while ( ( ele = queue2->getNext() ) != NULL ) 04427 { 04428 ox = ele->x ; 04429 oy = ele->y ; 04430 oz = ele->z ; 04431 04432 if ( noisevol->getDataAt( ox, oy, oz ) == 1 ) 04433 { 04434 visited ++ ; 04435 continue ; 04436 } 04437 04438 flag = 0 ; 04439 for ( int m = 0 ; m < 6 ; m ++ ) 04440 { 04441 int nx = ox + neighbor6[m][0] ; 04442 int ny = oy + neighbor6[m][1] ; 04443 int nz = oz + neighbor6[m][2] ; 04444 if ( getDataAt( nx, ny, nz ) > 0 && noisevol->getDataAt( nx, ny, nz ) == 1 ) 04445 { 04446 noisevol->setDataAt( ox, oy, oz, 1 ) ; 04447 visited ++ ; 04448 count ++ ; 04449 break ; 04450 } 04451 } 04452 } 04453 04454 if ( count == 0 ) 04455 { 04456 break ; 04457 } 04458 } 04459 printf("Maximum feature distance: %d Un-touched: %d\n", cur, queue2->getNumElements() - visited ) ; 04460 04461 04462 #endif 04463 /* Commented out for debugging 04464 04465 // First, 04466 // check for complex nodes in queue2 04467 // move them from queue2 to queue3 04468 queue2->reset() ; 04469 ele = queue2->getNext() ; 04470 while ( ele != NULL ) 04471 { 04472 ox = ele->x ; 04473 oy = ele->y ; 04474 oz = ele->z ; 04475 04476 // Check simple 04477 #ifndef NOISE_DIS_HELIX 04478 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) ) 04479 #else 04480 if ( isHelixEnd( ox, oy, oz, noisevol ) || ! isSimple( ox, oy, oz ) ) 04481 #endif 04482 { 04483 // Complex, set to next layer 04484 setDataAt( ox, oy, oz, curwid + 1 ) ; 04485 queue3->prepend( ox, oy, oz ) ; 04486 ele = queue2->remove() ; 04487 04488 numComplex ++ ; 04489 } 04490 else 04491 { 04492 ele = queue2->getNext() ; 04493 } 04494 } 04495 */ 04496 04497 // Next, 04498 // Compute score for each node left in queue2 04499 // move them into priority queue 04500 queue2->reset() ; 04501 ele = queue2->getNext() ; 04502 while ( ele != NULL ) 04503 { 04504 ox = ele->x ; 04505 oy = ele->y ; 04506 oz = ele->z ; 04507 04508 // Compute score 04509 score = getNumPotComplex2( ox, oy, oz ) ; 04510 scrvol->setDataAt( ox, oy, oz, score ) ; 04511 04512 // Push to queue 04513 gp = new gridPoint ; 04514 gp->x = ox ; 04515 gp->y = oy ; 04516 gp->z = oz ; 04517 // queue->add( gp, -score ) ; 04518 queue->add( gp, score ) ; 04519 04520 ele = queue2->remove() ; 04521 } 04522 04523 // Rename queue3 to be queue2, 04524 // Clear queue3 04525 // From now on, queue2 holds nodes of next level 04526 delete queue2 ; 04527 queue2 = queue3 ; 04528 queue3 = new GridQueue2( ) ; 04529 04530 // Next, start priority queue iteration 04531 while ( ! queue->isEmpty() ) 04532 { 04533 // Retrieve the node with the highest score 04534 queue->remove( gp, score ) ; 04535 ox = gp->x ; 04536 oy = gp->y ; 04537 oz = gp->z ; 04538 delete gp ; 04539 // score = -score ; 04540 04541 // Ignore the node 04542 // if it has been processed before 04543 // or it has an updated score 04544 if ( getDataAt( ox, oy, oz ) != curwid || (int) scrvol->getDataAt( ox, oy, oz ) != score ) 04545 { 04546 continue ; 04547 } 04548 04549 /* Commented out for debugging 04550 04551 // Remove this simple node 04552 setDataAt( ox, oy, oz, -1 ) ; 04553 numSimple ++ ; 04554 // printf("Highest score: %d\n", score) ; 04555 */ 04556 04557 /* Added for debugging */ 04558 // Check simple 04559 #ifndef NOISE_DIS_HELIX 04560 // if ( hasIsolatedEdge( ox, oy, oz ) && ! isNoiseHelixEnd( ox, oy, oz ) ) 04561 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) ) 04562 #else 04563 if ( isHelixEnd( ox, oy, oz ) || ! isSimple( ox, oy, oz ) ) 04564 #endif 04565 { 04566 // Complex, set to next layer 04567 setDataAt( ox, oy, oz, curwid + 1 ) ; 04568 queue4->prepend( ox, oy, oz ) ; 04569 numComplex ++ ; 04570 } 04571 else 04572 { 04573 setDataAt( ox, oy, oz, -1 ) ; 04574 numSimple ++ ; 04575 } 04576 /* Adding ends */ 04577 04578 // Move its neighboring unvisited node to queue2 04579 for ( int m = 0 ; m < 6 ; m ++ ) 04580 { 04581 int nx = ox + neighbor6[m][0] ; 04582 int ny = oy + neighbor6[m][1] ; 04583 int nz = oz + neighbor6[m][2] ; 04584 if ( getDataAt( nx, ny, nz ) == 0 ) 04585 { 04586 // setDataAt( nx, ny, nz, curwid + 1 ) ; 04587 queue2->prepend( nx, ny, nz ) ; 04588 } 04589 } 04590 04591 /* Commented out for debugging 04592 04593 // Find complex nodes in its 3x3 neighborhood 04594 // move them to queue2 04595 for ( i = -1 ; i < 2 ; i ++ ) 04596 for ( j = -1 ; j < 2 ; j ++ ) 04597 for ( k = -1 ; k < 2 ; k ++ ) 04598 { 04599 int nx = ox + i ; 04600 int ny = oy + j ; 04601 int nz = oz + k ; 04602 04603 // Check simple 04604 if ( getDataAt( nx, ny, nz ) == curwid && 04605 // ( isSheetEnd( ox, oy, oz ) || ! isSimple( nx, ny, nz )) ) 04606 #ifndef NOISE_DIS_HELIX 04607 ( isHelixEnd( nx, ny, nz ) || ! isSimple( nx, ny, nz ) ) ) 04608 #else 04609 ( isHelixEnd( nx, ny, nz, noisevol ) || ! isSimple( nx, ny, nz ) ) ) 04610 #endif 04611 04612 { 04613 // Complex, set to next layer 04614 setDataAt( nx, ny, nz, curwid + 1 ) ; 04615 queue2->prepend( nx, ny, nz ) ; 04616 numComplex ++ ; 04617 } 04618 } 04619 */ 04620 04621 // Update scores for nodes in its 5x5 neighborhood 04622 // insert them back into priority queue 04623 04624 for ( i = -2 ; i < 3 ;i ++ ) 04625 for ( j = -2 ; j < 3 ; j ++ ) 04626 for ( k = -2 ; k < 3 ; k ++ ) 04627 { 04628 int nx = ox + i ; 04629 int ny = oy + j ; 04630 int nz = oz + k ; 04631 04632 if ( getDataAt( nx, ny, nz ) == curwid ) 04633 { 04634 // Compute score 04635 score = getNumPotComplex2( nx, ny, nz ) ; 04636 04637 if ( score != (int) scrvol->getDataAt( nx, ny, nz ) ) 04638 { 04639 // printf("Update\n") ; 04640 scrvol->setDataAt( nx, ny, nz, score ) ; 04641 // Push to queue 04642 gp = new gridPoint ; 04643 gp->x = nx ; 04644 gp->y = ny ; 04645 gp->z = nz ; 04646 // queue->add( gp, -score ) ; 04647 queue->add( gp, score ) ; 04648 } 04649 } 04650 } 04651 04652 04653 } 04654 04655 #ifdef VERBOSE 04656 printf("%d complex, %d simple\n", numComplex, numSimple) ; 04657 #endif 04658 04659 if ( numSimple == 0 ) 04660 { 04661 break ; 04662 } 04663 } 04664 04665 // Finally, clean up 04666 delete scrvol; 04667 delete queue; 04668 delete queue2; 04669 delete queue3; 04670 delete queue4; 04671 #ifdef VERBOSE 04672 printf("Thresholding the volume to 0/1...\n") ; 04673 #endif 04674 threshold( 0, 0, 1 ) ; 04675 } 04676 04677 // Compute curve skeleton in 2D 04678 void Volume::curveSkeleton2D( float thr, Volume* svol ) 04679 { 04680 int i, j, k ; 04681 // First, threshold the volume 04682 #ifdef VERBOSE 04683 printf("Thresholding the volume to -1/0...\n") ; 04684 #endif 04685 threshold( thr, -1, 0 ) ; 04686 04687 // Next, apply convergent erosion 04688 // by preserving: complex nodes, curve end-points, and sheet points 04689 04690 // Next, initialize the linked queue 04691 #ifdef VERBOSE 04692 printf("Initializing queue...\n") ; 04693 #endif 04694 GridQueue2* queue2 = new GridQueue2( ) ; 04695 GridQueue2* queue3 = new GridQueue2( ) ; 04696 GridQueue2* queue4 = new GridQueue2( ) ; 04697 PriorityQueue <gridPoint,int> * queue = new PriorityQueue <gridPoint,int> ( MAX_QUEUELEN ); 04698 04699 for ( i = 0 ; i < getSizeX() ; i ++ ) 04700 for ( j = 0 ; j < getSizeY() ; j ++ ) 04701 for ( k = 0 ; k < getSizeZ() ; k ++ ) 04702 { 04703 if ( getDataAt( i, j, k ) >= 0 ) 04704 { 04705 if ( svol->getDataAt(i,j,k) > 0 ) 04706 { 04707 setDataAt( i, j, k, MAX_ERODE ) ; 04708 } 04709 else 04710 { 04711 for ( int m = 0 ; m < 4 ; m ++ ) 04712 { 04713 if ( getDataAt( i + neighbor4[m][0], j + neighbor4[m][1], k ) < 0 ) 04714 { 04715 // setDataAt( i, j, k, 1 ) ; 04716 queue2->prepend( i, j, k ) ; 04717 break ; 04718 } 04719 } 04720 } 04721 } 04722 } 04723 int wid = MAX_ERODE ; 04724 #ifdef VERBOSE 04725 printf("Total %d nodes\n", queue2->getNumElements() ) ; 04726 printf("Start erosion to %d...\n", wid) ; 04727 #endif 04728 04729 04730 // Perform erosion 04731 gridQueueEle* ele ; 04732 gridPoint* gp ; 04733 int ox, oy, oz ; 04734 int score ; 04735 Volume* scrvol = new Volume( this->getSizeX() , this->getSizeY(), this->getSizeZ() ) ; 04736 for ( i = 0 ; i < getSizeX() * getSizeY() * getSizeZ() ; i ++ ) 04737 { 04738 scrvol->setDataAt( i, -1 ) ; 04739 } 04740 04741 #ifdef NOISE_DIS_HELIX 04742 Volume* noisevol = new Volume( getSizeX(), getSizeY(), getSizeZ() ) ; 04743 #endif 04744 04745 for ( int curwid = 1 ; curwid <= wid ; curwid ++ ) 04746 { 04747 // At the start of each iteration, 04748 // queue2 holds all the nodes for this layer 04749 // queue3 and queue are empty 04750 04751 int numComplex = 0, numSimple = 0 ; 04752 #ifdef VERBOSE 04753 printf("Processing %d nodes in layer %d\n", queue2->getNumElements(), curwid) ; 04754 #endif 04755 04756 /* 04757 We first need to assign curwid + 1 to every node in this layer 04758