initial commit

CMakeLists.txt

@@ -0,0 +1,14 @@
+cmake_minimum_required(VERSION 2.8)
+ 
+project(FollicleSeparator)
+ 
+find_package(ITK REQUIRED)
+include(${ITK_USE_FILE})
+
+file(GLOB_RECURSE HEADERS  "src/*.h")
+file(GLOB_RECURSE INLINES  "src/*.inl")
+file(GLOB_RECURSE SOURCES  "src/*.cpp")
+
+add_executable(FollicleSeparator ${HEADERS} ${INLINES} ${SOURCES})
+ 
+target_link_libraries(FollicleSeparator ${ITK_LIBRARIES})

src/FollicleSeparator.cpp

@@ -0,0 +1,157 @@
+#pragma warning(push)
+#pragma warning(disable:4996)
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <itkBinaryThresholdImageFilter.h>
+#include <itkMaskNegatedImageFilter.h>
+#include <itkMinimumMaximumImageCalculator.h>
+#include <itkIntensityWindowingImageFilter.h>
+#include <itkSignedMaurerDistanceMapImageFilter.h>
+#include <itkMorphologicalWatershedImageFilter.h>
+#include <itkTimeProbe.h>
+#pragma warning(pop)
+
+#include <iostream>
+using namespace std;
+
+///
+/// Binary object splitting using
+/// watershed segmentation applied to the distance transform
+///
+int main(int argc, char* argv[])
+{
+  try
+  {
+    const int DIM = 3;
+    typedef itk::Image<unsigned char, DIM> UInt8Volume;   // Binary Mask
+    typedef itk::Image<short, DIM> Int16Volume;           // Input Volume
+    typedef itk::Image<unsigned short, DIM> UInt16Volume; // Label Volume
+    typedef itk::Image<float, DIM> Float32Volume;         // Distance
+
+    // Read labeled volume
+
+    itk::ImageFileReader<Int16Volume>::Pointer pReader = itk::ImageFileReader<Int16Volume>::New();
+    {
+      itk::TimeProbe clock;
+      clock.Start();
+
+      cout << "Read input from labeledVolume.vtk" << endl;
+      pReader->SetFileName("labeledVolume.vtk");
+      pReader->Update();
+
+      clock.Stop();
+      cout << "ImageFileReader: " << clock.GetTotal() << "s" << endl;
+    }
+
+    // Binarize input volume: Follicles = 0, Background = 255
+
+    itk::BinaryThresholdImageFilter<Int16Volume,UInt8Volume>::Pointer pBinarizationFilter =
+      itk::BinaryThresholdImageFilter<Int16Volume,UInt8Volume>::New();
+    {
+      itk::TimeProbe clock;
+      clock.Start();
+
+      pBinarizationFilter->SetInput(pReader->GetOutput());
+      pBinarizationFilter->SetLowerThreshold(0);
+      pBinarizationFilter->SetUpperThreshold(0);
+      pBinarizationFilter->SetInsideValue(255);
+      pBinarizationFilter->SetOutsideValue(0);
+      pBinarizationFilter->Update();
+
+      clock.Stop();
+      cout << "BinaryThresholdImageFilter: " << clock.GetTotal() << "s" << endl;
+    }
+
+    // Compute signed distance from background pixels
+    // Inside the follicle: positive
+    // Outside of the follicle: negative
+
+    itk::SignedMaurerDistanceMapImageFilter<UInt8Volume,Float32Volume>::Pointer pSignedMaurer =
+      itk::SignedMaurerDistanceMapImageFilter<UInt8Volume,Float32Volume>::New();
+    {
+      itk::TimeProbe clock;
+      clock.Start();
+
+      pSignedMaurer->SetInput(pBinarizationFilter->GetOutput());
+
+      pSignedMaurer->Update();
+
+      clock.Stop();
+      cout << "SignedMaurerDistanceMapImageFilter: " << clock.GetTotal() << "s" << endl;
+    }
+
+    // Invert distance map: [0 ... max distance] -> [255 ... 0]
+    // background pixels will have intensity 255
+    // pixels inside the follicle far away from the border will have low intensities
+    itk::IntensityWindowingImageFilter<Float32Volume,UInt8Volume>::Pointer pConvertToUInt8 =
+      itk::IntensityWindowingImageFilter<Float32Volume,UInt8Volume>::New();
+    {
+      itk::TimeProbe clock;
+      clock.Start();
+
+      itk::MinimumMaximumImageCalculator<Float32Volume>::Pointer pMinMaxDistance =
+        itk::MinimumMaximumImageCalculator<Float32Volume>::New();
+      pMinMaxDistance->SetImage(pSignedMaurer->GetOutput());
+      pMinMaxDistance->ComputeMaximum();
+      float maxDistance = pMinMaxDistance->GetMaximum();
+
+      pConvertToUInt8->SetInput(pSignedMaurer->GetOutput());
+      pConvertToUInt8->SetWindowMinimum(0.0);
+      pConvertToUInt8->SetWindowMaximum(maxDistance);
+      pConvertToUInt8->SetOutputMinimum(255);
+      pConvertToUInt8->SetOutputMaximum(0);
+      pConvertToUInt8->Update();
+
+      clock.Stop();
+      cout << "ConvertToUInt8: " << clock.GetTotal() << "s" << endl;
+
+      itk::ImageFileWriter<UInt8Volume>::Pointer pWriter = itk::ImageFileWriter<UInt8Volume>::New();
+      pWriter->SetInput(pConvertToUInt8->GetOutput());
+      pWriter->SetFileName("labeledVolume.distanceMap.vtk");
+      pWriter->Update();
+    }
+
+    // Watershed transformation of the distance map
+    itk::MorphologicalWatershedImageFilter<UInt8Volume,UInt16Volume>::Pointer pWatershed = itk::MorphologicalWatershedImageFilter<UInt8Volume,UInt16Volume>::New();
+    {
+      itk::TimeProbe clock;
+      clock.Start();
+
+      pWatershed->SetLevel(25); // high level -> low Number of regions
+      pWatershed->SetMarkWatershedLine(false);
+      pWatershed->SetInput(pConvertToUInt8->GetOutput());
+      pWatershed->Update();
+
+      clock.Stop();
+      cout << "WatershedImageFilter: " << clock.GetTotal() << "s" << endl;
+    }
+
+    // Mask out background voxels from the watershed result
+    itk::MaskNegatedImageFilter<UInt16Volume, UInt8Volume, UInt16Volume>::Pointer pMaskImageFilter =
+      itk::MaskNegatedImageFilter<UInt16Volume, UInt8Volume, UInt16Volume>::New();
+    {
+      pMaskImageFilter->SetInput(pWatershed->GetOutput());
+      pMaskImageFilter->SetMaskImage(pBinarizationFilter->GetOutput());
+      pMaskImageFilter->Update();
+    }
+
+    itk::MinimumMaximumImageCalculator<UInt16Volume>::Pointer pMaxLabelCalculator =
+      itk::MinimumMaximumImageCalculator<UInt16Volume>::New();
+    pMaxLabelCalculator->SetImage(pMaskImageFilter->GetOutput());
+    pMaxLabelCalculator->ComputeMaximum();
+    cout << "Number of labels: " << pMaxLabelCalculator->GetMaximum() << endl;
+
+    cout << "Write result to labeledVolume.watershed.vtk" << endl;
+    itk::ImageFileWriter<UInt16Volume>::Pointer pWriter = itk::ImageFileWriter<UInt16Volume>::New();
+    pWriter->SetInput(pMaskImageFilter->GetOutput());
+    pWriter->SetFileName("labeledVolume.watershed.vtk");
+    pWriter->Update();
+  }
+  catch (const itk::ExceptionObject& ex)
+  {
+    std::cerr << ex.GetDescription() << std::endl;
+    return 1;
+  }
+
+  return 0;
+}