From c1387ff3ede56f5a38afd62358e27eda91ad973a Mon Sep 17 00:00:00 2001
From: Julien Michel <julien.michel@c-s.fr>
Date: Thu, 30 Nov 2006 15:18:10 +0000
Subject: [PATCH] Ajout d'une optimisation pour le calcul du RCC8.

---
 ...naryImageMinimalBoundingRegionCalculator.h |   9 +-
 ...ryImageMinimalBoundingRegionCalculator.txx |  30 +++--
 .../otbImageToImageRCC8Calculator.txx         | 112 +++++-------------
 ...ryImageMinimalBoundingRegionCalculator.cxx |   2 +-
 4 files changed, 52 insertions(+), 101 deletions(-)

diff --git a/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.h b/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.h
index 7d7a268be4..154c1ba7d6 100644
--- a/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.h
+++ b/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.h
@@ -30,7 +30,7 @@ namespace otb
  *
  * This class is used for instance in the RCC8 calculator filter,
  * where the input region used for computation has to be the smallest possible
- * for costs reasons. The Pad flag allows the user to get a region 1 pixel larger
+ * for costs reasons. The Pad arg allows the user to get a region of pad  pixel larger
  * at each bound in case a security margin has to be kept. 
  *
  * \sa ImageToImageRCC8Calculator
@@ -55,12 +55,11 @@ public:
 	typedef typename InputImageType::RegionType RegionType;
 	typedef typename InputImageType::Pointer InputImagePointerType;
 	/** Toogle the pad option */
-	itkBooleanMacro(Pad);
 	itkGetMacro(Region,RegionType);
 	itkSetMacro(InsideValue,PixelType);
 	itkGetMacro(InsideValue,PixelType);
-	itkSetMacro(Pad,bool);
-	itkGetMacro(Pad,bool);
+	itkSetMacro(Pad,unsigned int);
+	itkGetMacro(Pad,unsigned int);
 	
 protected:
 	/** Constructor */
@@ -77,7 +76,7 @@ private:
 	/** The computed region */
 	RegionType m_Region;
 	/** Toogle if pad wanted */
-	bool m_Pad;
+	unsigned int  m_Pad;
 	/** Inside value */
 	PixelType m_InsideValue;
 };
diff --git a/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.txx b/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.txx
index bd1a826d35..c4b239b7b0 100644
--- a/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.txx
+++ b/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.txx
@@ -31,13 +31,15 @@ template <class TInputImage>
 BinaryImageMinimalBoundingRegionCalculator<TInputImage>
 ::BinaryImageMinimalBoundingRegionCalculator() 
 {
-  m_Pad=false;
+  // The pad option is desactivated by default
+  m_Pad=0;
+  // Set the default region
   typename InputImageType::SizeType size;
   typename InputImageType::IndexType index;
   size[0]=0;
   size[1]=0;
   index[0]=0;
-  index[1]=1;
+  index[1]=0;
   m_Region.SetSize(size);
   m_Region.SetIndex(index);
   m_InsideValue=static_cast<PixelType>(255);
@@ -118,28 +120,32 @@ BinaryImageMinimalBoundingRegionCalculator<TInputImage>
 	    rit.PreviousSlice();
 	  }
       }
-    otbMsgDebugMacro(<<"BinaryImageMinimalBoundingBoxCalculator: min "<<min);
-    otbMsgDebugMacro(<<"BinaryImageMinimalBoundingBoxCalculator: max "<<max);
+    
+    // Compute size and index of the region
     typename InputImageType::SizeType size;
     typename InputImageType::IndexType index;
     RegionType maxRegion = image->GetLargestPossibleRegion();
-
-    if(m_Pad)
+    // If the pad option is activated
+    if(m_Pad>0)
       {
 	for(int i=0;i<InputImageType::ImageDimension;i++)
 	  {
+	    // If we are not on boundary case, we can do what we want
 	    if(min[i]> maxRegion.GetIndex()[i])
 	      {
-		index[i]= min[i]-1;
+		index[i]= min[i]-m_Pad;
 	      }
+	    // else we are at beginning of the image, so don't pad
 	    else
 	      {
 		index[i]= maxRegion.GetIndex()[i]; 
 	      }
-	    if (index[i]+max[i]-min[i]+2<=maxRegion.GetIndex()[i]+maxRegion.GetSize()[i])
+	    // If are not on boundary case, we can pad the size
+	    if (index[i]+max[i]-min[i]+2*m_Pad+1<=maxRegion.GetIndex()[i]+maxRegion.GetSize()[i])
 	      {
-		size[i]=max[i]-min[i]+2;
+		size[i]=max[i]-min[i]+2*m_Pad+1;
 	      }
+	    // Else we must restrain ourselves to the image boundaries
 	    else 
 	      {
 		size[i]=maxRegion.GetSize()[i]+maxRegion.GetIndex()[i]
@@ -148,6 +154,7 @@ BinaryImageMinimalBoundingRegionCalculator<TInputImage>
 	  }
       }
     else
+      // If the pad option is not activated, the result is simple
       {
 	for(int i=0;i<InputImageType::ImageDimension;i++)
 	  {
@@ -155,6 +162,9 @@ BinaryImageMinimalBoundingRegionCalculator<TInputImage>
 	    index[i]=min[i];
 	  }
       }
+    otbMsgDebugMacro(<<"BinaryImageMinimalBoundingBoxCalculator: index "<<index);
+    otbMsgDebugMacro(<<"BinaryImageMinimalBoundingBoxCalculator: size "<<size);
+    // Set the size and index of the output region
     m_Region.SetIndex(index);
     m_Region.SetSize(size);
 }
@@ -165,7 +175,7 @@ template <class TInputImage>
 void
 BinaryImageMinimalBoundingRegionCalculator<TInputImage>
 ::PrintSelf( std::ostream& os,itk::Indent indent ) const
-  {
+  { 
     Superclass::PrintSelf(os,indent);
   }
 
diff --git a/Code/SpatialReasoning/otbImageToImageRCC8Calculator.txx b/Code/SpatialReasoning/otbImageToImageRCC8Calculator.txx
index 08e52a6e35..e62f40cfd8 100644
--- a/Code/SpatialReasoning/otbImageToImageRCC8Calculator.txx
+++ b/Code/SpatialReasoning/otbImageToImageRCC8Calculator.txx
@@ -26,7 +26,7 @@
 #include "itkAndImageFilter.h"
 #include "itkImageRegionIterator.h"
 #include "itkImageRegionConstIterator.h"
-#include "itkImageSliceConstIteratorWithIndex.h"
+#include "otbBinaryImageMinimalBoundingRegionCalculator.h"
 #include "otbMacro.h"
 
 namespace otb
@@ -116,91 +116,33 @@ namespace otb
     // Input images pointers
     typename ImageType::Pointer image1 = this->GetInput1();
     typename ImageType::Pointer image2 = this->GetInput2();
-    // Iterator definition
-    typedef itk::ImageSliceConstIteratorWithIndex<ImageType> SliceIteratorType;
-    // Indexes containing upper-left and lower-right corner
-    typename ImageType::IndexType min;
-    typename ImageType::IndexType max;
-    min[0]=0;
-    min[1]=0;
-    max[1]=0;
-    max[1]=0;
-    for ( unsigned int axis = 0; axis < ImageType::ImageDimension; axis++ )
-      { // Create the forward iterator to find lower bound
-	SliceIteratorType fit1(image1,image1->GetLargestPossibleRegion());
-	SliceIteratorType fit2(image2,image2->GetLargestPossibleRegion());
-	fit1.SetFirstDirection( !axis );
-	fit1.SetSecondDirection( axis );
-	fit1.GoToBegin();
-	fit2.SetFirstDirection( !axis );
-	fit2.SetSecondDirection( axis );
-	fit2.GoToBegin();
-	// Walk through the two images line by line
-	while (!fit1.IsAtEnd()&&!fit2.IsAtEnd())
-	  {
-	    while (!fit1.IsAtEndOfSlice()&&!fit2.IsAtEndOfSlice())
-	      {
-		while(!fit1.IsAtEndOfLine()&&!fit2.IsAtEndOfLine())
-		  {
-		    // If a common intersection is found
-		    if ((fit1.Get()==m_InsideValue)||(fit2.Get()==m_InsideValue))
-		      {
-			// then the lower bound is found
-			min[axis]=fit1.GetIndex()[axis];
-			fit1.GoToReverseBegin(); // skip to the end
-			fit2.GoToReverseBegin();
-			break;
-		      }
-		    ++fit1;
-		    ++fit2;
-		  }
-		fit1.NextLine();
-		fit2.NextLine();
-	      }
-	    fit1.NextSlice();
-	    fit2.NextSlice();
-	  }
-	// Create the reverse iterator to find upper bound
-	SliceIteratorType rit1(image1,image1->GetLargestPossibleRegion());
-	SliceIteratorType rit2(image2,image2->GetLargestPossibleRegion());
-	rit1.SetFirstDirection(!axis);
-	rit1.SetSecondDirection(axis);
-	rit1.GoToReverseBegin();
-	rit2.SetFirstDirection(!axis);
-	rit2.SetSecondDirection(axis);
-	rit2.GoToReverseBegin();
-	// Walk through the two images line by line
-	while (!rit1.IsAtReverseEnd()&&!rit2.IsAtReverseEnd())
-	  {
-	    while (!rit1.IsAtReverseEndOfSlice()&&!rit2.IsAtReverseEndOfSlice())
-	      {
-		while (!rit1.IsAtReverseEndOfLine()&&!rit2.IsAtReverseEndOfLine())
-		  {
-		    // If a common intersection is found
-		    if ((rit1.Get()==m_InsideValue)||(rit2.Get()==m_InsideValue))
-		      {
-			max[axis]=rit1.GetIndex()[axis];
-			rit1.GoToBegin();
-			rit2.GoToBegin(); //Skip to reverse end
-			break;
-		      }
-		    --rit1;
-		    --rit2;
-		  }
-		rit1.PreviousLine();
-		rit2.PreviousLine();
-	      }
-	    rit1.PreviousSlice();
-	    rit2.PreviousSlice();
-	  }
-      }
-    typename ImageType::RegionType region;
+    typename ImageType::RegionType region1, region2, region;
+    typedef otb::BinaryImageMinimalBoundingRegionCalculator<ImageType> RegionCalculator;
+    typename RegionCalculator::Pointer rc =RegionCalculator::New();
+    rc->SetInput(image1);
+    rc->SetPad(2);
+    rc->SetInsideValue(this->GetInsideValue());
+    rc->Update();
+    region1=rc->GetRegion();
+    rc=RegionCalculator::New();
+    rc->SetInput(image2);
+    rc->SetPad(2);
+    rc->SetInsideValue(this->GetInsideValue());
+    rc->Update();
+    region2=rc->GetRegion();
     typename ImageType::SizeType size;
-    size[0]=max[0]-min[0];
-    size[1]=max[1]-min[1];
-    region.SetIndex(min);
+    typename ImageType::IndexType index;
+    
+    for(int i=0;i<ImageType::ImageDimension;i++)
+      {
+	index[i]=std::max(region1.GetIndex()[i],region2.GetIndex()[i]);
+	int potSize = std::min(region1.GetIndex()[i]+region1.GetSize()[i],
+			 region2.GetIndex()[i]+region2.GetSize()[i]);
+	size[i]=(potSize-index[i]<0 ? 0 : potSize-index[i]);
+      }
+    region.SetIndex(index);
     region.SetSize(size);
-    otbMsgDebugMacro(<<"RCC8Calculator->ComputeMinimalRegion(): index: "<<min<<" size: "<<size);
+    otbMsgDebugMacro(<<"RCC8Calculator->ComputeMinimalRegion(): index: "<<index<<" size: "<<size);
     return region;
   }
 /**
@@ -447,7 +389,7 @@ ImageToImageRCC8Calculator<TInputImage>
     /// First we compute the minimal region of interest we will use for the relation computation
     m_MinimalROI=this->ComputeMinimalRegion();
     /// If they are disjoint, the answer is trivial
-    if((m_MinimalROI.GetSize()[0]==0)||(m_MinimalROI.GetSize()[1]==0))
+    if((m_MinimalROI.GetSize()[0]<=1)||(m_MinimalROI.GetSize()[1]<=1))
       {
 	/// The relation is DC
 	m_Value=OTB_RCC8_DC;
diff --git a/Testing/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.cxx b/Testing/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.cxx
index 6eb794571e..b2930628f9 100644
--- a/Testing/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.cxx
+++ b/Testing/Code/BasicFilters/otbBinaryImageMinimalBoundingRegionCalculator.cxx
@@ -70,7 +70,7 @@ try
 for(IteratorType it=images->Begin();it!=images->End();++it)
       {
 	brct = BoundingRegionCalculatorType::New();
-	brct->SetPad(true);
+	brct->SetPad(1);
 	brct->SetInput(it.Get());
 	brct->Update();
 	RegionType region = brct->GetRegion();
-- 
GitLab