EMAN2
skeletonizer.cpp
Go to the documentation of this file.
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                 //}