]> Creatis software - clitk.git/blobdiff - itk/clitkSegmentationUtils.txx
Merge branch 'master' of /home/dsarrut/clitk3.server
[clitk.git] / itk / clitkSegmentationUtils.txx
index b432188aec93aa8e34a21f1047ea3d22991a2c25..edb208a7064770fb3cb737107948ff06983395de 100644 (file)
@@ -353,8 +353,9 @@ namespace clitk {
     sliceRelPosFilter->AddOrientationTypeString(orientation);
     sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1));
     sliceRelPosFilter->SetIntermediateSpacing(spacing);
-    sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent);
-    sliceRelPosFilter->SetUseASingleObjectConnectedComponentBySliceFlag(singleObjectCCL);
+    sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent);
+    sliceRelPosFilter->ObjectCCLSelectionFlagOff();
+    sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL);
     //    sliceRelPosFilter->SetInverseOrientationFlag(inverseflag); 
     sliceRelPosFilter->SetAutoCropFlag(autocropFlag); 
     sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
@@ -363,6 +364,7 @@ namespace clitk {
   }
   //--------------------------------------------------------------------
 
+
   //--------------------------------------------------------------------
   template<class ImageType>
   bool
@@ -377,6 +379,7 @@ namespace clitk {
   }
   //--------------------------------------------------------------------
 
+
   //--------------------------------------------------------------------
   template<class ImageType>
   bool
@@ -729,7 +732,8 @@ namespace clitk {
     dilateFilter->SetForegroundValue(FG);
     dilateFilter->SetBoundaryToForeground(false);
     dilateFilter->SetKernel(structuringElement);
-    dilateFilter->SetInput(output);
+    if (extendSupport) dilateFilter->SetInput(output);
+    else dilateFilter->SetInput(image);
     dilateFilter->Update();
     return dilateFilter->GetOutput();
   }
@@ -816,6 +820,7 @@ namespace clitk {
                                               int mainDirection, 
                                               double offsetToKeep)
   {
+    assert((mainDirection==0) || (mainDirection==1));
     typedef itk::ImageSliceIteratorWithIndex<ImageType> SliceIteratorType;
     SliceIteratorType siter = SliceIteratorType(input, 
                                                 input->GetLargestPossibleRegion());
@@ -830,9 +835,6 @@ namespace clitk {
     while ((i<lA.size()) && (!siter.IsAtEnd())) {
       // Check that the current slice correspond to the current point
       input->TransformIndexToPhysicalPoint(siter.GetIndex(), C);
-      // DD(C);
-      // DD(i);
-      // DD(lA[i]);
       if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm
       }
       else {
@@ -840,9 +842,6 @@ namespace clitk {
         A = lA[i];
         B = lB[i];
         C = A;
-        // DD(A);
-        // DD(B);
-        // DD(C);
       
         // Check that the line is not a point (A=B)
         bool p = (A[0] == B[0]) && (A[1] == B[1]);
@@ -876,6 +875,7 @@ namespace clitk {
   }                                                   
   //--------------------------------------------------------------------
 
+
   //--------------------------------------------------------------------
   template<class ImageType>
   void 
@@ -904,6 +904,62 @@ namespace clitk {
   //--------------------------------------------------------------------
 
 
+  //--------------------------------------------------------------------
+  template<class ImageType>
+  void 
+  And(ImageType * input, 
+      const ImageType * object, 
+      typename ImageType::PixelType BG)
+  {
+    typename ImageType::Pointer o;
+    bool resized=false;
+    if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
+      o = clitk::ResizeImageLike<ImageType>(object, input, BG);
+      resized = true;
+    }
+
+    typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
+    typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
+    boolFilter->InPlaceOn();
+    boolFilter->SetInput1(input);
+    if (resized) boolFilter->SetInput2(o);  
+    else boolFilter->SetInput2(object);
+    boolFilter->SetBackgroundValue1(BG);
+    boolFilter->SetBackgroundValue2(BG);
+    boolFilter->SetOperationType(BoolFilterType::And);
+    boolFilter->Update();
+  }
+  //--------------------------------------------------------------------
+
+
+  //--------------------------------------------------------------------
+  template<class ImageType>
+  void 
+  Or(ImageType * input, 
+     const ImageType * object, 
+     typename ImageType::PixelType BG)
+  {
+    typename ImageType::Pointer o;
+    bool resized=false;
+    if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, object)) {
+      o = clitk::ResizeImageLike<ImageType>(object, input, BG);
+      resized = true;
+    }
+
+    typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
+    typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); 
+    boolFilter->InPlaceOn();
+    boolFilter->SetInput1(input);
+    if (resized) boolFilter->SetInput2(o);  
+    else boolFilter->SetInput2(object);
+    boolFilter->SetBackgroundValue1(BG);
+    boolFilter->SetBackgroundValue2(BG);
+    boolFilter->SetOperationType(BoolFilterType::Or);
+    boolFilter->Update();
+  }
+  //--------------------------------------------------------------------
+
+
   //--------------------------------------------------------------------
   template<class ImageType>
   typename ImageType::Pointer
@@ -1033,7 +1089,7 @@ namespace clitk {
                                                             extremaDirection, extremaOppositeFlag, p);
       if (found) {
         position2D[i] = p;
-      }    
+      }
     }
     
     // Convert 2D points in slice into 3D points
@@ -1046,8 +1102,8 @@ namespace clitk {
       p[lineDirection] += 10;
       B.push_back(p);
       // Margins ?
-      A[i][1] += margin;
-      B[i][1] += margin;
+      A[i][extremaDirection] += margin;
+      B[i][extremaDirection] += margin;
     }
 
   }
@@ -1095,5 +1151,181 @@ namespace clitk {
   //--------------------------------------------------------------------
 
 
+  //--------------------------------------------------------------------
+  /* Consider an input object, start at A, for each slice (dim1): 
+     - compute the intersection between the AB line and the current slice
+     - remove what is at lower or greater according to dim2 of this point
+     - stop at B
+  */
+  template<class ImageType>
+  typename ImageType::Pointer
+  SliceBySliceSetBackgroundFromSingleLine(const ImageType * input, 
+                                          typename ImageType::PixelType BG, 
+                                          typename ImageType::PointType & A, 
+                                          typename ImageType::PointType & B, 
+                                          int dim1, int dim2, bool removeLowerPartFlag)
+    
+  {
+    // Extract slices
+    typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
+    typedef typename SliceType::Pointer SlicePointer;
+    std::vector<SlicePointer> slices;
+    clitk::ExtractSlices<ImageType>(input, dim1, slices);
+
+    // Start at slice that contains A, and stop at B
+    typename ImageType::IndexType Ap;
+    typename ImageType::IndexType Bp;
+    input->TransformPhysicalPointToIndex(A, Ap);
+    input->TransformPhysicalPointToIndex(B, Bp);
+    
+    // Determine slice largest region
+    typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
+    typename SliceType::SizeType size = region.GetSize();
+    typename SliceType::IndexType index = region.GetIndex();
+
+    // Line slope
+    double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]);
+    double b = Ap[dim2];
+
+    // Loop from slice A to slice B
+    for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) {
+      // Compute intersection between line AB and current slice for the dim2
+      double p = a*i+b;
+      // Change region (lower than dim2)
+      if (removeLowerPartFlag) {
+        size[dim2] = p-Ap[dim2];
+      }
+      else {
+        size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p;
+        index[dim2] = p;
+      }
+      region.SetSize(size);
+      region.SetIndex(index);
+      // Fill region with BG (simple region iterator)
+      FillRegionWithValue<SliceType>(slices[i+Ap[dim1]], BG, region);
+      /*
+      typedef itk::ImageRegionIterator<SliceType> IteratorType;
+      IteratorType iter(slices[i+Ap[dim1]], region);
+      iter.GoToBegin();
+      while (!iter.IsAtEnd()) {
+        iter.Set(BG);
+        ++iter;
+      }
+      */
+      // Loop
+    }
+    
+    // Merge slices
+    typename ImageType::Pointer output;
+    output = clitk::JoinSlices<ImageType>(slices, input, dim1);
+    return output;
+  }
+  //--------------------------------------------------------------------
+
+  //--------------------------------------------------------------------
+  /* Consider an input object, slice by slice, use the point A and set
+     pixel to BG according to their position relatively to A
+  */
+  template<class ImageType>
+  typename ImageType::Pointer
+  SliceBySliceSetBackgroundFromPoints(const ImageType * input, 
+                                      typename ImageType::PixelType BG, 
+                                      int sliceDim,
+                                      std::vector<typename ImageType::PointType> & A, 
+                                      bool removeGreaterThanXFlag,
+                                      bool removeGreaterThanYFlag)
+    
+  {
+    // Extract slices
+    typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
+    typedef typename SliceType::Pointer SlicePointer;
+    std::vector<SlicePointer> slices;
+    clitk::ExtractSlices<ImageType>(input, sliceDim, slices);
+
+    // Start at slice that contains A
+    typename ImageType::IndexType Ap;
+    
+    // Determine slice largest region
+    typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion();
+    typename SliceType::SizeType size = region.GetSize();
+    typename SliceType::IndexType index = region.GetIndex();
+
+    // Loop from slice A to slice B
+    for(uint i=0; i<A.size(); i++) {
+      input->TransformPhysicalPointToIndex(A[i], Ap);
+      uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2];
+      if ((sliceIndex < 0) || (sliceIndex >= slices.size())) {
+        continue; // do not consider this slice
+      }
+      
+      // Compute region for BG
+      if (removeGreaterThanXFlag) {
+        index[0] = Ap[0];
+        size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]);
+      }
+      else {
+        index[0] = region.GetIndex()[0];
+        size[0] = Ap[0] - index[0];
+      }
+
+      if (removeGreaterThanYFlag) {
+        index[1] = Ap[1];
+        size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]);
+      }
+      else {
+        index[1] = region.GetIndex()[1];
+        size[1] = Ap[1] - index[1];
+      }
+
+      // Set region
+      region.SetSize(size);
+      region.SetIndex(index);
+
+      // Fill region with BG (simple region iterator)
+      FillRegionWithValue<SliceType>(slices[sliceIndex], BG, region);
+      // Loop
+    }
+    
+    // Merge slices
+    typename ImageType::Pointer output;
+    output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
+    return output;
+  }
+  //--------------------------------------------------------------------
+
+
+  //--------------------------------------------------------------------
+  template<class ImageType>
+  void
+  FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
+  {
+    typedef itk::ImageRegionIterator<ImageType> IteratorType;
+    IteratorType iter(input, region);
+    iter.GoToBegin();
+    while (!iter.IsAtEnd()) {
+      iter.Set(value);
+      ++iter;
+    }    
+  }
+  //--------------------------------------------------------------------
+
+
+  //--------------------------------------------------------------------
+  template<class ImageType>
+  void
+  GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min, 
+                    typename ImageType::PointType & max)
+  {
+    typedef typename ImageType::PointType PointType;
+    typedef typename ImageType::IndexType IndexType;
+    IndexType min_i, max_i;
+    min_i = input->GetLargestPossibleRegion().GetIndex();
+    for(uint i=0; i<ImageType::ImageDimension; i++)
+      max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
+    input->TransformIndexToPhysicalPoint(min_i, min);
+    input->TransformIndexToPhysicalPoint(max_i, max);  
+  }
+  //--------------------------------------------------------------------
+
 } // end of namespace