Change preset 3 to 400/20 (like in [Chapet et al])

[clitk.git] / itk / clitkSegmentationUtils.txx

diff --git a/itk/clitkSegmentationUtils.txx b/itk/clitkSegmentationUtils.txx
index 94e915431c72a78b082515d1002f802398ed4397..c8f976e69449f2d95456e1a9c1a3a1c6d4e77064 100644 (file)
--- a/itk/clitkSegmentationUtils.txx
+++ b/itk/clitkSegmentationUtils.txx
@@ -31,6 +31,7 @@
 #include <itkBinaryBallStructuringElement.h>
 #include <itkBinaryDilateImageFilter.h>
 #include <itkConstantPadImageFilter.h>
+#include <itkImageSliceIteratorWithIndex.h>
 
 //--------------------------------------------------------------------
 template<class ImageType>
@@ -191,6 +192,53 @@ clitk::Labelize(const ImageType * input,
   relabelFilter->SetMinimumObjectSize(minimalComponentSize);
   relabelFilter->Update();
 
+  // DD(relabelFilter->GetNumberOfObjects());
+  // DD(relabelFilter->GetOriginalNumberOfObjects());
+  // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
+
+  // Return result
+  typename ImageType::Pointer output = relabelFilter->GetOutput();
+  return output;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+/*
+  Warning : in this cas, we consider outputType like inputType, not
+  InternalImageType. Be sure it fits.
+ */
+template<class ImageType>
+typename ImageType::Pointer
+clitk::LabelizeAndCountNumberOfObjects(const ImageType * input, 
+                                       typename ImageType::PixelType BG, 
+                                       bool isFullyConnected, 
+                                       int minimalComponentSize, 
+                                       int & nb) {
+  // InternalImageType for storing large number of component
+  typedef itk::Image<int, ImageType::ImageDimension> InternalImageType;
+  
+  // Connected Component label 
+  typedef itk::ConnectedComponentImageFilter<ImageType, InternalImageType> ConnectFilterType;
+  typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New();
+  //  connectFilter->ReleaseDataFlagOn(); 
+  connectFilter->SetInput(input);
+  connectFilter->SetBackgroundValue(BG);
+  connectFilter->SetFullyConnected(isFullyConnected);
+  
+  // Sort by size and remove too small area.
+  typedef itk::RelabelComponentImageFilter<InternalImageType, ImageType> RelabelFilterType;
+  typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
+  //  relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ???
+  relabelFilter->SetInput(connectFilter->GetOutput());
+  relabelFilter->SetMinimumObjectSize(minimalComponentSize);
+  relabelFilter->Update();
+
+  nb = relabelFilter->GetNumberOfObjects();
+  // DD(relabelFilter->GetNumberOfObjects());
+  // DD(relabelFilter->GetOriginalNumberOfObjects());
+  // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]);
+
   // Return result
   typename ImageType::Pointer output = relabelFilter->GetOutput();
   return output;
@@ -198,6 +246,7 @@ clitk::Labelize(const ImageType * input,
 //--------------------------------------------------------------------
 
 
+
 //--------------------------------------------------------------------
 template<class ImageType>
 typename ImageType::Pointer
@@ -669,3 +718,73 @@ void clitk::ConvertOption(std::string optionName, uint given,
 //--------------------------------------------------------------------
 
 
+//--------------------------------------------------------------------
+/*
+  http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
+  Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
+  (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
+  This will equal zero if the point C is on the line formed by
+  points A and B, and will have a different sign depending on the
+  side. Which side this is depends on the orientation of your (x,y)
+  coordinates, but you can plug test values for A,B and C into this
+  formula to determine whether negative values are to the left or to
+  the right.
+  => to accelerate, start with formula, when change sign -> stop and fill
+*/
+template<class ImageType>
+void 
+clitk::SliceBySliceSetBackgroundFromLineSeparation(typename ImageType::Pointer input, 
+                                                   std::vector<typename ImageType::PointType> & lA, 
+                                                   std::vector<typename ImageType::PointType> & lB, 
+                                                   typename ImageType::PixelType BG, 
+                                                   int mainDirection, 
+                                                   double offsetToKeep)
+{
+  
+  typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
+  SliceIteratorType siter = SliceIteratorType(input, 
+                                              input->GetLargestPossibleRegion());
+  siter.SetFirstDirection(0);
+  siter.SetSecondDirection(1);
+  siter.GoToBegin();
+  int i=0;
+  typename ImageType::PointType A;
+  typename ImageType::PointType B;
+  typename ImageType::PointType C;
+  while (!siter.IsAtEnd()) {
+    // Check that the current slice correspond to the current point
+    input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
+    if (C[2] != lA[i][2]) {
+      // DD(C);
+      // DD(lA[i]);
+    }
+    else {
+      // Define A,B,C points
+      A = lA[i];
+      B = lB[i];
+      C = A;
+      C[mainDirection] += offsetToKeep; // I know I must keep this point
+      double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
+      bool isPositive = s<0;
+      while (!siter.IsAtEndOfSlice()) {
+        while (!siter.IsAtEndOfLine()) {
+          // Very slow, I know ... but image should be very small
+          input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
+          double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
+          if (s == 0) siter.Set(BG); // on the line, we decide to remove
+          if (isPositive) {
+            if (s > 0) siter.Set(BG);
+          }
+          else {
+            if (s < 0) siter.Set(BG); 
+          }
+          ++siter;
+        }
+        siter.NextLine();
+      }
+      ++i;
+    }
+    siter.NextSlice();
+  }
+}                                                   
+//--------------------------------------------------------------------