X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkAddRelativePositionConstraintToLabelImageFilter.txx;h=8e12c225e59852544637c5db492d6b26afb532df;hb=72375037f90c596a034b2ebe5e54e209e7b45511;hp=5e7011ad37f861ef900a7ecf6d083378f657261f;hpb=38786c4da19b87319bbe3cecc145e3d1771d10da;p=clitk.git

diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx
index 5e7011a..8e12c22 100644
--- a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx
+++ b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx
@@ -25,12 +25,14 @@
 
 // itk
 #include <deque>
-#include "itkStatisticsLabelMapFilter.h"
-#include "itkLabelImageToStatisticsLabelMapFilter.h"
-#include "itkRegionOfInterestImageFilter.h"
-#include "itkBinaryThresholdImageFilter.h"
-#include "itkBinaryErodeImageFilter.h"
-#include "itkBinaryBallStructuringElement.h"
+#include <itkStatisticsLabelMapFilter.h>
+#include <itkLabelImageToStatisticsLabelMapFilter.h>
+#include <itkRegionOfInterestImageFilter.h>
+#include <itkBinaryThresholdImageFilter.h>
+#include <itkBinaryErodeImageFilter.h>
+#include <itkBinaryBallStructuringElement.h>
+#include <itkAddImageFilter.h>
+#include <itkDivideByConstantImageFilter.h>
 
 // itk [Bloch et al] 
 #include "RelativePositionPropImageFilter.h"
@@ -46,10 +48,13 @@ AddRelativePositionConstraintToLabelImageFilter():
   SetFuzzyThreshold(0.6);
   SetBackgroundValue(0);
   SetObjectBackgroundValue(0);
-  SetOrientationType(LeftTo);
+  ClearOrientationType();
   ResampleBeforeRelativePositionFilterOn();
   SetIntermediateSpacing(10);
-  AutoCropOn();
+  AutoCropFlagOn();
+  InverseOrientationFlagOff();
+  RemoveObjectFlagOn();
+  CombineWithOrFlagOff();
 }
 //--------------------------------------------------------------------
 
@@ -78,6 +83,51 @@ SetInputObject(const ImageType * image)
 //--------------------------------------------------------------------
   
 
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+ClearOrientationType() 
+{
+  m_OrientationTypeString.clear();
+  m_OrientationType.clear();
+  m_Angle1.clear();
+  m_Angle2.clear();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+int
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+GetNumberOfAngles()
+{
+  return m_OrientationType.size();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void 
+clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
+AddOrientationTypeString(std::string t) 
+{
+  m_OrientationTypeString.push_back(t);
+  switch (t[0]) {
+  case 'L' : AddOrientationType(LeftTo); break;
+  case 'R' : AddOrientationType(RightTo);break;
+  case 'A' : AddOrientationType(AntTo);break;
+  case 'P' : AddOrientationType(PostTo);break;
+  case 'S' : AddOrientationType(SupTo);break;
+  case 'I' : AddOrientationType(InfTo);break;
+  default: clitkExceptionMacro("Error, you must provide L,R or A,P or S,I");
+  }
+}
+//--------------------------------------------------------------------
+  
+
 //--------------------------------------------------------------------
 template <class ImageType>
 void 
@@ -112,10 +162,11 @@ GenerateInputRequestedRegion()
 template <class ImageType>
 void 
 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
-SetAngle1(double a) 
+AddAngles(double a, double b) 
 {
-  SetOrientationType(Angle);
-  m_Angle1 = a;
+  AddOrientationTypeString("Angle");
+  m_Angle1.push_back(a);
+  m_Angle2.push_back(b);
 }
 //--------------------------------------------------------------------
 
@@ -124,29 +175,35 @@ SetAngle1(double a)
 template <class ImageType>
 void 
 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
-SetAngle2(double a) 
+AddOrientationType(OrientationTypeEnumeration orientation) 
 {
-  SetOrientationType(Angle);
-  m_Angle2 = a;
-}
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template <class ImageType>
-void 
-clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
-SetOrientationType(OrientationTypeEnumeration orientation) 
-{
-  m_OrientationType = orientation;
-  switch (m_OrientationType) {
-  case LeftTo:   m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(0);   break;
-  case RightTo:  m_Angle1 = clitk::deg2rad(180); m_Angle2 = clitk::deg2rad(0);   break;
-  case AntTo:    m_Angle1 = clitk::deg2rad(90);  m_Angle2 = clitk::deg2rad(0);   break;
-  case PostTo:   m_Angle1 = clitk::deg2rad(-90); m_Angle2 = clitk::deg2rad(0);   break;
-  case InfTo:    m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(90);  break;
-  case SupTo:    m_Angle1 = clitk::deg2rad(0);   m_Angle2 = clitk::deg2rad(-90); break;
-  case Angle:  break;      
+  m_OrientationType.push_back(orientation);
+  switch (orientation) {
+  case LeftTo:   
+    m_Angle1.push_back(clitk::deg2rad(0));   
+    m_Angle2.push_back(clitk::deg2rad(0));
+    break;
+  case RightTo:  
+    m_Angle1.push_back(clitk::deg2rad(180)); 
+    m_Angle2.push_back(clitk::deg2rad(0));
+    break;
+  case AntTo:
+    m_Angle1.push_back(clitk::deg2rad(90));
+    m_Angle2.push_back(clitk::deg2rad(0));
+    break;
+  case PostTo:
+    m_Angle1.push_back(clitk::deg2rad(-90)); 
+    m_Angle2.push_back(clitk::deg2rad(0));
+    break;
+  case InfTo:    
+    m_Angle1.push_back(clitk::deg2rad(0));   
+    m_Angle2.push_back(clitk::deg2rad(90));
+    break;
+  case SupTo:    
+    m_Angle1.push_back(clitk::deg2rad(0));   
+    m_Angle2.push_back(clitk::deg2rad(-90));
+    break;
+  case Angle:  break;
   }
   /*         A1   A2
              Left      0    0
@@ -166,6 +223,10 @@ void
 clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>::
 GenerateData() 
 {
+  if (GetNumberOfAngles() <1) {
+    clitkExceptionMacro("Add at least one orientation type");
+  }  
+
   // Get input pointer
   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
@@ -173,7 +234,7 @@ GenerateData()
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   static const unsigned int dim = ImageType::ImageDimension;
-  StartNewStep("Initial resample and pad");  
+  StartNewStep("Initial resample");  
   // Step 1 : resample
   if (m_ResampleBeforeRelativePositionFilter) {
     typedef clitk::ResampleImageWithOptionsFilter<ImageType> ResampleFilterType;
@@ -181,36 +242,56 @@ GenerateData()
     resampleFilter->SetInput(object);
     resampleFilter->SetOutputIsoSpacing(m_IntermediateSpacing);
     resampleFilter->SetGaussianFilteringEnabled(false);
-    // resampleFilter->SetVerboseOptions(true);
+    //    resampleFilter->SetVerboseOptions(true);
     resampleFilter->Update();
     working_image = resampleFilter->GetOutput();
   }
   else {
     working_image = object;
   }
+  StopCurrentStep<ImageType>(working_image);
 
   // Step 2: object pad to input image -> we want to compute the
   // relative position for each point belonging to the input image
   // domain, so we have to extend (pad) the object image to fit the
   // domain size
   if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input, working_image)) {
+    StartNewStep("Pad object to image size");  
     typename ImageType::Pointer output = ImageType::New();
     SizeType size;
     for(unsigned int i=0; i<dim; i++) {
-      size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
+      size[i] = lrint((input->GetLargestPossibleRegion().GetSize()[i]*
+                       input->GetSpacing()[i])/(double)working_image->GetSpacing()[i]);
     }
+
+    // The index of the input is not necessarily zero, so we have to
+    // take it into account (not done)
     RegionType region;
+    IndexType index = input->GetLargestPossibleRegion().GetIndex();
     region.SetSize(size);
+    for(unsigned int i=0; i<dim; i++) {
+      if (index[i] != 0) {
+        std::cerr << "Index diff from zero : " << index << ". not done yet !" << std::endl;
+        exit(0);
+      }
+    }
     // output->SetLargestPossibleRegion(region);
     output->SetRegions(region);
-    output->SetSpacing(working_image->GetSpacing());
-    output->SetOrigin(input->GetOrigin());
+    output->SetSpacing(working_image->GetSpacing());    
+    PointType origin = input->GetOrigin();
+    for(unsigned int i=0; i<dim; i++) {
+      origin[i] = index[i]*input->GetSpacing()[i] + input->GetOrigin()[i];
+    }
+    output->SetOrigin(origin);
+    //    output->SetOrigin(input->GetOrigin());
+
     output->Allocate();
     output->FillBuffer(m_BackgroundValue);
     typename PadFilterType::Pointer padFilter = PadFilterType::New();
-    typename PadFilterType::InputImageIndexType index;
+    // typename PadFilterType::InputImageIndexType index;
     for(unsigned int i=0; i<dim; i++) {
-      index[i] = lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
+      index[i] = -index[i]*input->GetSpacing()[i]/(double)working_image->GetSpacing()[i]
+        + lrint((working_image->GetOrigin()[i] - input->GetOrigin()[i])/working_image->GetSpacing()[i]);
     }
     padFilter->SetSourceImage(working_image);
     padFilter->SetDestinationImage(output);
@@ -218,27 +299,63 @@ GenerateData()
     padFilter->SetSourceRegion(working_image->GetLargestPossibleRegion());
     padFilter->Update();
     working_image = padFilter->GetOutput();
+    StopCurrentStep<ImageType>(working_image);
   }
   else {
     // DD("[debug] RelPos : same size and spacing : no padding");
   }
   // Keep object image (with resampline and pad)
   object_resampled = working_image;
-  StopCurrentStep<ImageType>(working_image);
+  //  StopCurrentStep<ImageType>(working_image);
 
   // Step 3: compute rel pos in object
   StartNewStep("Relative Position Map");  
   typedef itk::RelativePositionPropImageFilter<ImageType, FloatImageType> RelPosFilterType;
-  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
-  relPosFilter->SetInput(working_image);
-  relPosFilter->SetAlpha1(m_Angle1); // xy plane
-  relPosFilter->SetAlpha2(m_Angle2);
-  relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
-  relPosFilter->SetFast(true);
-  relPosFilter->SetRadius(1); // seems sufficient in this cas
-  // relPosFilter->SetVerboseProgress(true);
-  relPosFilter->Update();
-  relPos = relPosFilter->GetOutput();
+  typename RelPosFilterType::Pointer relPosFilter;
+
+  typename FloatImageType::Pointer m_FuzzyMap;
+  for(int i=0; i<GetNumberOfAngles(); i++) {
+    // Compute fuzzy map
+    relPosFilter = RelPosFilterType::New();
+    relPosFilter->SetInput(working_image);
+    relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
+    relPosFilter->SetAlpha2(m_Angle2[i]);
+    relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
+    relPosFilter->SetFast(true);
+    relPosFilter->SetRadius(1); // seems sufficient in this case
+    // relPosFilter->SetVerboseProgress(true);
+    relPosFilter->Update();
+    relPos = relPosFilter->GetOutput();
+
+    if (GetNumberOfAngles() != 1) {
+      // Creation of the first m_FuzzyMap
+      if (i==0) {
+        m_FuzzyMap = clitk::NewImageLike<FloatImageType>(relPos, true);
+        m_FuzzyMap->FillBuffer(0.0);
+      }
+      
+      // Add to current fuzzy map
+      typedef itk::AddImageFilter<FloatImageType, FloatImageType, FloatImageType> AddImageFilter;
+      typename AddImageFilter::Pointer addFilter = AddImageFilter::New();
+      addFilter->SetInput1(m_FuzzyMap);
+      addFilter->SetInput2(relPos);
+      addFilter->Update();
+      m_FuzzyMap = addFilter->GetOutput();
+    }
+    else m_FuzzyMap = relPos;
+  }
+
+  // Divide by the number of relpos
+  if (GetNumberOfAngles() != 1) {
+    typedef itk::DivideByConstantImageFilter<FloatImageType, float, FloatImageType> DivideFilter;
+    typename DivideFilter::Pointer divideFilter = DivideFilter::New();
+    divideFilter->SetInput(m_FuzzyMap);
+    divideFilter->SetConstant(GetNumberOfAngles());
+    divideFilter->Update();
+    m_FuzzyMap = divideFilter->GetOutput();
+  }
+
+  relPos = m_FuzzyMap;
   StopCurrentStep<FloatImageType>(relPos);
                
   //--------------------------------------------------------------------
@@ -288,22 +405,24 @@ GenerateData()
     // resampleFilter->SetVerboseOptions(true);
     resampleFilter->Update();
     working_image = resampleFilter->GetOutput();
-  StopCurrentStep<ImageType>(working_image);
+    StopCurrentStep<ImageType>(working_image);
   }
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   // Pre Step 6: pad if not the same size : it can occur when downsample and upsample
+  //if (!HaveSameSizeAndSpacing(working_image, input)) {
   if (working_image->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
     StartNewStep("Pad to get the same size than input");
     typename ImageType::Pointer temp = ImageType::New();
     temp->CopyInformation(input);
     temp->SetRegions(input->GetLargestPossibleRegion()); // Do not forget !!
     temp->Allocate();
-    temp->FillBuffer(m_BackgroundValue);
-    typename PadFilterType::Pointer padFilter2 = PadFilterType::New(); // if yes : redo relpos
+    temp->FillBuffer(m_BackgroundValue); 
+    typename PadFilterType::Pointer padFilter2 = PadFilterType::New();
     padFilter2->SetSourceImage(working_image);
     padFilter2->SetDestinationImage(temp);
+    // DD(input->GetLargestPossibleRegion().GetIndex());
     padFilter2->SetDestinationIndex(input->GetLargestPossibleRegion().GetIndex());
     padFilter2->SetSourceRegion(working_image->GetLargestPossibleRegion());
     padFilter2->Update();
@@ -320,34 +439,59 @@ GenerateData()
   StartNewStep("Combine with initial input (boolean And)");
   typedef clitk::BooleanOperatorLabelImageFilter<ImageType> BoolFilterType;
   typename BoolFilterType::Pointer combineFilter = BoolFilterType::New();
-  writeImage<ImageType>(input, "i.mhd");
-  writeImage<ImageType>(working_image, "w.mhd");
   combineFilter->SetBackgroundValue(m_BackgroundValue);
   combineFilter->SetBackgroundValue1(m_BackgroundValue);
   combineFilter->SetBackgroundValue2(m_BackgroundValue);
   combineFilter->SetForegroundValue(m_BackgroundValue+1);
   combineFilter->SetInput1(input);
   combineFilter->SetInput2(working_image);
-  combineFilter->SetOperationType(BoolFilterType::And);
-  combineFilter->InPlaceOn();
+  if (GetInverseOrientationFlag())
+    combineFilter->SetOperationType(BoolFilterType::AndNot);
+  else {
+    if (GetCombineWithOrFlag())
+      combineFilter->SetOperationType(BoolFilterType::Or);
+    else
+      combineFilter->SetOperationType(BoolFilterType::And);
+  }
+  combineFilter->InPlaceOff(); // Do not modify initial input (!)
   combineFilter->Update(); 
   working_image = combineFilter->GetOutput();
-  // writeImage<ImageType>(working_image, "res.mhd");
- 
-  combineFilter = BoolFilterType::New();
-  combineFilter->SetInput1(working_image);
-  combineFilter->SetInput2(object);
-  combineFilter->SetOperationType(BoolFilterType::AndNot);
-  combineFilter->InPlaceOn();
-  combineFilter->Update(); 
 
-  working_image = combineFilter->GetOutput();
+  // Remove (if needed the object from the support)
+  if (GetRemoveObjectFlag()) {
+    combineFilter = BoolFilterType::New();
+    combineFilter->SetInput1(working_image);
+    combineFilter->SetInput2(object);
+    combineFilter->SetOperationType(BoolFilterType::AndNot);
+    combineFilter->InPlaceOn();
+    combineFilter->Update(); 
+    working_image = combineFilter->GetOutput();
+  }
+  // In the other case, we must *add* the initial object to keep it
+  // but not more than the initial support
+  else { 
+    combineFilter = BoolFilterType::New();
+    combineFilter->SetInput1(working_image);
+    combineFilter->SetInput2(object);
+    combineFilter->SetOperationType(BoolFilterType::Or);
+    combineFilter->InPlaceOn();
+    combineFilter->Update(); 
+    working_image = combineFilter->GetOutput(); // not needed because InPlaceOn ?
+    combineFilter = BoolFilterType::New();
+    combineFilter->SetInput1(working_image);
+    combineFilter->SetInput2(input);
+    combineFilter->SetOperationType(BoolFilterType::And);
+    combineFilter->InPlaceOn();
+    combineFilter->Update(); 
+    working_image = combineFilter->GetOutput();
+  }
+
   StopCurrentStep<ImageType>(working_image);
 
   //--------------------------------------------------------------------
   //--------------------------------------------------------------------
   // Step 7: autocrop
-  if (GetAutoCrop()) {
+  if (GetAutoCropFlag()) {
     StartNewStep("Final AutoCrop");
     typedef clitk::AutoCropFilter<ImageType> CropFilterType;
     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
@@ -363,7 +507,7 @@ GenerateData()
   
   // Final Step -> set output
   this->SetNthOutput(0, working_image);
-  return;
+  //  this->GraftOutput(working_image);
 }
 //--------------------------------------------------------------------