00001 // Copyright (C) 2005-2008 Washington University in St Louis, Baylor College of Medicine. All rights reserved 00002 // Author: Sasakthi S. Abeysinghe (sasakthi@gmail.com) 00003 // Description: Performs skeletonization on a grayscale volume 00004 00005 #include "skeletonizer.h" 00006 #include <list> 00007 using std::list; 00008 00009 using namespace wustl_mm::GraySkeletonCPP; 00010 00011 const char VolumeSkeletonizer::THINNING_CLASS_SURFACE_PRESERVATION = 4; 00012 const char VolumeSkeletonizer::THINNING_CLASS_CURVE_PRESERVATION_2D = 3; 00013 const char VolumeSkeletonizer::THINNING_CLASS_CURVE_PRESERVATION = 2; 00014 const char VolumeSkeletonizer::THINNING_CLASS_POINT_PRESERVATION = 1; 00015 const char VolumeSkeletonizer::THINNING_CLASS_TOPOLOGY_PRESERVATION = 0; 00016 const char VolumeSkeletonizer::PRUNING_CLASS_PRUNE_SURFACES = 5; 00017 const char VolumeSkeletonizer::PRUNING_CLASS_PRUNE_CURVES = 6; 00018 const char VolumeSkeletonizer::PRUNING_CLASS_PRUNE_POINTS = 7; 00019 00020 VolumeSkeletonizer::VolumeSkeletonizer(int pointRadius, int curveRadius, int surfaceRadius, int skeletonDirectionRadius) { 00021 // math = new MathLib(); 00022 // surfaceNormalFinder = new NormalFinder(); 00023 this->pointRadius = pointRadius; 00024 this->curveRadius = curveRadius; 00025 this->surfaceRadius = surfaceRadius; 00026 this->skeletonDirectionRadius = skeletonDirectionRadius; 00027 00028 // gaussianFilterPointRadius.radius = pointRadius; 00029 // math->GetBinomialDistribution(gaussianFilterPointRadius); 00030 // 00031 // gaussianFilterCurveRadius.radius = curveRadius; 00032 // math->GetBinomialDistribution(gaussianFilterCurveRadius); 00033 // 00034 // gaussianFilterSurfaceRadius.radius = surfaceRadius; 00035 // math->GetBinomialDistribution(gaussianFilterSurfaceRadius); 00036 // 00037 // gaussianFilterMaxRadius.radius = MAX_GAUSSIAN_FILTER_RADIUS; 00038 // math->GetBinomialDistribution(gaussianFilterMaxRadius); 00039 // 00040 // uniformFilterSkeletonDirectionRadius.radius = skeletonDirectionRadius; 00041 // math->GetUniformDistribution(uniformFilterSkeletonDirectionRadius); 00042 } 00043 00044 VolumeSkeletonizer::~VolumeSkeletonizer() { 00045 //delete math; 00046 //delete surfaceNormalFinder; 00047 } 00048 // **************************** Added by Ross **************************************** 00049 00050 const float CURVE_VAL = 1.0f; 00051 const float SURFACE_VAL = 2.0f; 00052 const float MAP_ERR_VAL = 100.0f; 00053 00054 bool VolumeSkeletonizer::Are26Neighbors(Vec3<int> u, Vec3<int> v) { 00055 if ( u==v || abs(u[0]-v[0])>1 || abs(u[1]-v[1])>1 || abs(u[2]-v[2])>1) { 00056 return false; 00057 } else { 00058 return true; 00059 } 00060 } 00061 00062 // Based on NonManifoldMesh<TVertex, TEdge, TFace>::NonManifoldMesh(Volume * sourceVol) 00063 void VolumeSkeletonizer::MarkSurfaces(Volume* skeleton) { 00064 00065 int faceNeighbors[3][3][3] = { {{1,0,0}, {1,0,1}, {0,0,1}}, 00066 {{1,0,0}, {1,1,0}, {0,1,0}}, 00067 {{0,1,0}, {0,1,1}, {0,0,1}} }; 00068 int indices[4]; 00069 bool faceFound; 00070 00071 for (int z = 0; z < skeleton->getSizeZ(); z++) { 00072 for (int y = 0; y < skeleton->getSizeY(); y++) { 00073 for (int x = 0; x < skeleton->getSizeX(); x++) { 00074 00075 indices[0] = skeleton->getIndex(x,y,z); 00076 for (int n = 0; n < 3; n++) { 00077 faceFound = true; 00078 for (int m = 0; m < 3; m++) { 00079 indices[m+1] = skeleton->getIndex(x+faceNeighbors[n][m][0], y+faceNeighbors[n][m][1], z+faceNeighbors[n][m][2]); 00080 faceFound = faceFound && (skeleton->getDataAt(indices[m+1]) > 0); 00081 } 00082 if (faceFound) { 00083 for (int m = 0; m < 4; m++) { 00084 skeleton->setDataAt(indices[m], SURFACE_VAL); 00085 } 00086 } 00087 } 00088 00089 } 00090 } 00091 } 00092 } 00093 00094 void VolumeSkeletonizer::CleanUpSkeleton(Volume * skeleton, int minNumVoxels, float valueThreshold) { 00095 00096 //Get the indices of voxels that are above the threshold, i.e. part of the skeleton 00097 list< Vec3<int> > skel_indices; 00098 Vec3<int> voxel_indices; 00099 for (int k=0; k < skeleton->getSizeZ(); k++) { 00100 for (int j=0; j < skeleton->getSizeY(); j++) { 00101 for (int i=0; i < skeleton->getSizeX(); i++) { 00102 if (skeleton->getDataAt(i,j,k) > valueThreshold) { 00103 voxel_indices.set_value(i,j,k); 00104 skel_indices.push_front(voxel_indices); 00105 } 00106 } 00107 } 00108 } 00109 00110 vector< Vec3<int> > segment; //the coordinates for a set of connected skeleton voxels 00111 list< Vec3<int> >::iterator itr; 00112 00113 /* 00114 1. Group the connected voxels together 00115 2. Check the number of voxels in that group 00116 3. If below the minimum number of voxels, remove them from the skeleton 00117 4. Repeat until all the voxels in the skeleton have been grouped (possibly alone if not connected) 00118 */ 00119 while (skel_indices.size() > 0) { 00120 segment.clear(); 00121 segment.push_back(skel_indices.front()); 00122 skel_indices.pop_front(); 00123 // group connected voxels -- each member of segment neighbors at least one other member of segment 00124 //For each voxel in segment, we test if each voxel in skel_indices is a neighbor 00125 for (unsigned int seg_ix=0; seg_ix < segment.size(); seg_ix++) { 00126 for (itr = skel_indices.begin(); itr != skel_indices.end(); itr++) { 00127 00128 if (Are26Neighbors(segment[seg_ix], *itr)) { 00129 segment.push_back(*itr); 00130 skel_indices.erase(itr); 00131 } 00132 } 00133 } //Now, a segment is complete 00134 00135 00136 //If the region of connected voxels is too small, remove them from the map 00137 if (segment.size() < static_cast<unsigned int>(minNumVoxels)) { 00138 for (unsigned int ix=0; ix<segment.size(); ix++) { 00139 skeleton->setDataAt(segment[ix][0], segment[ix][1], segment[ix][2],0.0f); 00140 } 00141 } 00142 00143 } 00144 } 00145 // **************************** End: added by Ross **************************************** 00146 00147 00148 Volume * VolumeSkeletonizer::PerformPureJuSkeletonization(Volume * imageVol, string, double threshold, int minCurveWidth, int minSurfaceWidth) { 00149 imageVol->pad(MAX_GAUSSIAN_FILTER_RADIUS, 0); 00150 Volume * preservedVol = new Volume(imageVol->getSizeX(), imageVol->getSizeY(), imageVol->getSizeZ()); 00151 Volume * surfaceVol; 00152 Volume * curveVol; 00153 Volume * topologyVol; 00154 //printf("\t\t\tUSING THRESHOLD : %f\n", threshold); 00155 // Skeletonizing while preserving surface features curve features and topology 00156 surfaceVol = GetJuSurfaceSkeleton(imageVol, preservedVol, threshold); 00157 PruneSurfaces(surfaceVol, minSurfaceWidth); 00158 VoxelOr(preservedVol, surfaceVol); 00159 curveVol = VolumeSkeletonizer::GetJuCurveSkeleton(imageVol, preservedVol, threshold, true); 00160 VolumeSkeletonizer::PruneCurves(curveVol, minCurveWidth); 00161 VoxelOr(preservedVol, curveVol); 00162 00163 topologyVol = VolumeSkeletonizer::GetJuTopologySkeleton(imageVol, preservedVol, threshold); 00164 00165 //Code below by Ross as a test -- to replace GetJuTopologySkeleton return value 00166 // int curveVolMax = curveVol->getVolumeData()->GetMaxIndex(); 00167 // int surfaceVolMax = curveVol->getVolumeData()->GetMaxIndex(); 00168 // int maximum = curveVolMax <= surfaceVolMax? curveVolMax : surfaceVolMax; 00169 // if (curveVolMax != surfaceVolMax) 00170 // cout << "Curve Skeleton: " << curveVolMax << '\n' << "Surface Skeleton" << surfaceVolMax << endl; 00171 // topologyVol = new Volume(curveVol->getSizeX(), curveVol->getSizeY(), curveVol->getSizeZ()); 00172 // float val, cval, sval; 00173 // for (int i = 0; i < maximum; i++) 00174 // { 00175 // cval = float(curveVol->getDataAt(i)); 00176 // sval = float(surfaceVol->getDataAt(i)); 00177 // if (cval && sval) 00178 // val = 100; //Something went wrong! 00179 // else if (cval) 00180 // val = 1; 00181 // else if (sval) 00182 // val = -1; 00183 // else 00184 // val = 0; 00185 // topologyVol->setDataAt(i, val); 00186 // } 00187 00188 00189 00190 00191 00192 imageVol->pad(-MAX_GAUSSIAN_FILTER_RADIUS, 0); 00193 topologyVol->pad(-MAX_GAUSSIAN_FILTER_RADIUS, 0); 00194 delete preservedVol; 00195 delete surfaceVol; 00196 delete curveVol; 00197 return topologyVol; 00198 } 00199 00200 Volume * VolumeSkeletonizer::GetJuCurveSkeleton(Volume * sourceVolume, Volume * preserve, double threshold, bool is3D){ 00201 char thinningClass = is3D ? THINNING_CLASS_CURVE_PRESERVATION : THINNING_CLASS_CURVE_PRESERVATION_2D; 00202 return GetJuThinning(sourceVolume, preserve, threshold, thinningClass); 00203 } 00204 00205 Volume * VolumeSkeletonizer::GetJuSurfaceSkeleton(Volume * sourceVolume, Volume * preserve, double threshold){ 00206 return GetJuThinning(sourceVolume, preserve, threshold, THINNING_CLASS_SURFACE_PRESERVATION); 00207 } 00208 00209 Volume * VolumeSkeletonizer::GetJuTopologySkeleton(Volume * sourceVolume, Volume * preserve, double threshold){ 00210 return GetJuThinning(sourceVolume, preserve, threshold, THINNING_CLASS_TOPOLOGY_PRESERVATION); 00211 } 00212 00213 Volume * VolumeSkeletonizer::GetJuThinning(Volume * sourceVolume, Volume * preserve, double threshold, char thinningClass) { 00214 Volume * thinnedVolume = new Volume(sourceVolume->getSizeX(), sourceVolume->getSizeY(), sourceVolume->getSizeZ(), 0, 0, 0, sourceVolume); 00215 switch(thinningClass) { 00216 case THINNING_CLASS_SURFACE_PRESERVATION : 00217 thinnedVolume->surfaceSkeletonPres((float)threshold, preserve); 00218 break; 00219 case THINNING_CLASS_CURVE_PRESERVATION : 00220 thinnedVolume->curveSkeleton((float)threshold, preserve); 00221 break; 00222 case THINNING_CLASS_CURVE_PRESERVATION_2D : 00223 thinnedVolume->curveSkeleton2D((float)threshold, preserve); 00224 break; 00225 case THINNING_CLASS_TOPOLOGY_PRESERVATION : 00226 thinnedVolume->skeleton((float)threshold, preserve, preserve); 00227 } 00228 return thinnedVolume; 00229 } 00230 00231 void VolumeSkeletonizer::PruneCurves(Volume * sourceVolume, int pruneLength) { 00232 sourceVolume->erodeHelix(pruneLength); 00233 } 00234 void VolumeSkeletonizer::PruneSurfaces(Volume * sourceVolume, int pruneLength) { 00235 sourceVolume->erodeSheet(pruneLength); 00236 } 00237 00238 void VolumeSkeletonizer::VoxelOr(Volume * sourceAndDestVolume1, Volume * sourceVolume2){ 00239 if(sourceVolume2 != NULL) { 00240 for(int x = 0; x < sourceAndDestVolume1->getSizeX(); x++) { 00241 for(int y = 0; y < sourceAndDestVolume1->getSizeY(); y++) { 00242 for(int z = 0; z < sourceAndDestVolume1->getSizeZ(); z++) { 00243 sourceAndDestVolume1->setDataAt(x, y, z, max(sourceAndDestVolume1->getDataAt(x, y, z), sourceVolume2->getDataAt(x, y, z))); 00244 } 00245 } 00246 } 00247 } 00248 } 00249 00250 00251 00252 00253 00254 00255 00256 00257 00258 00259 00260 00261 00262 //Volume * VolumeSkeletonizer::PerformImmersionSkeletonizationAndPruning(Volume * sourceVol, Volume * preserveVol, double startGray, double endGray, double stepSize, int smoothingIterations, int smoothingRadius, int minCurveSize, int minSurfaceSize, int maxCurveHole, int maxSurfaceHole, string outputPath, bool doPruning, double pointThreshold, double curveThreshold, double surfaceThreshold) { 00263 //appTimeManager.PushCurrentTime(); 00264 //for(int i = 0; i < smoothingIterations; i++) { 00265 //SmoothenVolume(sourceVol, startGray, endGray, smoothingRadius); 00266 //} 00267 //appTimeManager.PopAndDisplayTime("Smoothing : %f seconds!\n"); 00268 //Vector3DFloat * volumeGradient = NULL; 00269 //EigenResults3D * volumeEigens; 00270 //sourceVol->pad(MAX_GAUSSIAN_FILTER_RADIUS, 0); 00271 //if(preserveVol != NULL) { 00272 //preserveVol->pad(MAX_GAUSSIAN_FILTER_RADIUS, 0); 00273 //} 00274 00275 //if(doPruning) { 00276 //volumeGradient = GetVolumeGradient(sourceVol); 00277 //} 00278 00279 //Volume * nullVol = new Volume(sourceVol->getSizeX(), sourceVol->getSizeY(), sourceVol->getSizeZ()); 00280 //appTimeManager.PushCurrentTime(); 00281 //Volume * surfaceVol = GetImmersionThinning(sourceVol, preserveVol, startGray, endGray, stepSize, THINNING_CLASS_SURFACE_PRESERVATION); 00282 //appTimeManager.PopAndDisplayTime("Surface Thinning : %f seconds!\n"); 00283 00284 //#ifdef SAVE_INTERMEDIATE_RESULTS 00285 //surfaceVol->toMRCFile((char *)(outputPath + "-S-Pre-Prune-Pre-Erode.mrc").c_str()); 00286 //#endif 00287 00288 //PruneSurfaces(surfaceVol, minSurfaceSize); 00289 00290 //appTimeManager.PushCurrentTime(); 00291 //if(doPruning) { 00292 //#ifdef SAVE_INTERMEDIATE_RESULTS 00293 //surfaceVol->toMRCFile((char *)(outputPath + "-S-Pre-Prune.mrc").c_str()); 00294 //WriteVolumeToVRMLFile(surfaceVol, outputPath + "-S-Pre-Prune.wrl"); 00295 //#endif 00296 //appTimeManager.PushCurrentTime(); 00297 //volumeEigens = GetEigenResults(surfaceVol, volumeGradient, gaussianFilterSurfaceRadius, surfaceRadius, true); 00298 //appTimeManager.PopAndDisplayTime(" Getting Eigens : %f seconds!\n"); 00299 00300 //appTimeManager.PushCurrentTime(); 00301 //Volume * prunedSurfaceVol = new Volume(surfaceVol->getSizeX(), surfaceVol->getSizeY(), surfaceVol->getSizeZ(), 0, 0, 0, surfaceVol); 00302 //appTimeManager.PopAndDisplayTime(" Getting Copy of surface : %f seconds!\n"); 00303 00304 00305 //appTimeManager.PushCurrentTime(); 00306 //PruneUsingStructureTensor(prunedSurfaceVol, sourceVol, preserveVol, volumeGradient, volumeEigens, gaussianFilterSurfaceRadius, surfaceThreshold, PRUNING_CLASS_PRUNE_SURFACES, outputPath + "-S"); 00307 //appTimeManager.PopAndDisplayTime(" Pruning: %f seconds!\n"); 00308 00309 //appTimeManager.PushCurrentTime(); 00310 //delete [] volumeEigens; 00311 //#ifdef SAVE_INTERMEDIATE_RESULTS 00312 //prunedSurfaceVol->toMRCFile((char *)(outputPath + "-S-Post-Prune.mrc").c_str()); 00313 //#endif 00314 00315 //delete surfaceVol; 00316 //surfaceVol = prunedSurfaceVol; 00317 //appTimeManager.PopAndDisplayTime(" Memory Cleanup: %f seconds!\n"); 00318 00319 //} 00320 00321 //PruneSurfaces(surfaceVol, minSurfaceSize); 00322 //appTimeManager.PopAndDisplayTime("Surface Pruning : %f seconds!\n"); 00323 00324 //#ifdef SAVE_INTERMEDIATE_RESULTS 00325 //surfaceVol->toMRCFile((char *)(outputPath + "-S-Post-Erosion.mrc").c_str()); 00326 //#endif 00327 00328 //Volume * cleanedSurfaceVol = GetJuSurfaceSkeleton(surfaceVol, nullVol, 0.5); 00329 //PruneSurfaces(cleanedSurfaceVol, minSurfaceSize); 00330 //#ifdef SAVE_INTERMEDIATE_RESULTS 00331 //cleanedSurfaceVol->toMRCFile((char *)(outputPath + "-S-Cleaned.mrc").c_str()); 00332 //#endif 00333 00334 //delete surfaceVol; 00335 //surfaceVol = cleanedSurfaceVol; 00336 //VoxelOr(surfaceVol, preserveVol); 00337 00338 //appTimeManager.PushCurrentTime(); 00339 00340 //Volume * curveVol = GetImmersionThinning(sourceVol, surfaceVol, startGray, endGray, stepSize, THINNING_CLASS_CURVE_PRESERVATION); 00341 //appTimeManager.PopAndDisplayTime("Curve Thinning : %f seconds!\n"); 00342 00343 //#ifdef SAVE_INTERMEDIATE_RESULTS 00344 //curveVol->toMRCFile((char *)(outputPath + "-C-Pre-Prune_Pre-Erode.mrc").c_str()); 00345 //#endif 00346 00347 //PruneCurves(curveVol, minCurveSize); 00348 //VoxelBinarySubtract(curveVol, surfaceVol); 00349 00350 //appTimeManager.PushCurrentTime(); 00351 //if(doPruning) { 00352 //#ifdef SAVE_INTERMEDIATE_RESULTS 00353 //curveVol->toMRCFile((char *)(outputPath + "-C-Pre-Prune.mrc").c_str()); 00354 //#endif 00355 00356 //volumeEigens = GetEigenResults(curveVol, volumeGradient, gaussianFilterCurveRadius, curveRadius, true); 00357 //Volume * prunedCurveVol = new Volume(curveVol->getSizeX(), curveVol->getSizeY(), curveVol->getSizeZ(), 0, 0, 0, curveVol); 00358 //PruneUsingStructureTensor(prunedCurveVol, sourceVol, preserveVol, volumeGradient, volumeEigens, gaussianFilterCurveRadius, curveThreshold, PRUNING_CLASS_PRUNE_CURVES, outputPath + "-C"); 00359 //delete [] volumeEigens; 00360 //#ifdef SAVE_INTERMEDIATE_RESULTS 00361 //prunedCurveVol->toMRCFile((char *)(outputPath + "-C-Post-Prune.mrc").c_str()); 00362 //#endif 00363 00364 //Volume * filledCurveVol = FillCurveHoles(prunedCurveVol, curveVol, maxCurveHole); 00365 //#ifdef SAVE_INTERMEDIATE_RESULTS 00366 //filledCurveVol->toMRCFile((char *)(outputPath + "-C-Post-Fill.mrc").c_str()); 00367 //#endif 00368 //delete curveVol; 00369 //delete prunedCurveVol; 00370 //curveVol = filledCurveVol; 00371 00372 00373 //} 00374 00375 //VoxelOr(curveVol, surfaceVol); 00376 //PruneCurves(curveVol, minCurveSize); 00377 //appTimeManager.PopAndDisplayTime("Curve Pruning : %f seconds!\n"); 00378 //#ifdef SAVE_INTERMEDIATE_RESULTS 00379 //curveVol->toMRCFile((char *)(outputPath + "-C-Post-Erosion.mrc").c_str()); 00380 //#endif 00381 00382 //Volume * cleanedCurveVol = GetJuCurveSkeleton(curveVol, surfaceVol, 0.5, true); 00383 //PruneCurves(cleanedCurveVol, minCurveSize); 00384 //#ifdef SAVE_INTERMEDIATE_RESULTS 00385 //cleanedCurveVol->toMRCFile((char *)(outputPath + "-C-Cleaned.mrc").c_str()); 00386 //#endif 00387 00388 //delete curveVol; 00389 //curveVol = cleanedCurveVol; 00390 00391 //VoxelOr(curveVol, surfaceVol); 00392 //#ifdef SAVE_INTERMEDIATE_RESULTS 00393 //curveVol->toMRCFile((char *)(outputPath + "-SC.mrc").c_str()); 00394 //#endif 00395 00405 00407 00408 //delete surfaceVol; 00410 //delete nullVol; 00411 //delete [] volumeGradient; 00412 00413 //#ifdef SAVE_INTERMEDIATE_RESULTS 00414 //curveVol->toOFFCells2((char *)(outputPath + "-SC.off").c_str()); 00415 //#endif 00416 00417 //sourceVol->pad(-MAX_GAUSSIAN_FILTER_RADIUS, 0); 00418 //curveVol->pad(-MAX_GAUSSIAN_FILTER_RADIUS, 0); 00419 //return curveVol; 00420 //}