From 9b4daebb476dd805386a2dccb4b02588419bab6c Mon Sep 17 00:00:00 2001
From: Sebastien Harasse <sebastien.harasse@c-s.fr>
Date: Tue, 17 Apr 2012 17:43:47 +0200
Subject: [PATCH] ENH: Added iteration map output to the mean shift filter. Use
 GetIterationOutput(). Updated tests accordingly.

---
 Code/BasicFilters/otbMeanShiftImageFilter2.h  |  17 ++-
 .../BasicFilters/otbMeanShiftImageFilter2.txx | 130 +++++++++++-------
 Code/IO/otbSpotImageMetadataInterface.cxx     |   2 +-
 Testing/Code/BasicFilters/CMakeLists.txt      |   5 +-
 .../BasicFilters/otbMeanShiftImageFilter2.cxx |  23 +++-
 5 files changed, 109 insertions(+), 68 deletions(-)

diff --git a/Code/BasicFilters/otbMeanShiftImageFilter2.h b/Code/BasicFilters/otbMeanShiftImageFilter2.h
index fb06eef495..ded4608810 100644
--- a/Code/BasicFilters/otbMeanShiftImageFilter2.h
+++ b/Code/BasicFilters/otbMeanShiftImageFilter2.h
@@ -18,12 +18,11 @@
 #ifndef __otbMeanShiftImageFilter2_h
 #define __otbMeanShiftImageFilter2_h
 
-
+#include "otbImage.h"
 #include "itkImageToImageFilter.h"
 #include <vcl_algorithm.h>
 /*
 #include "itkVariableLengthVector.h"
-#include "otbImage.h"
 #include "otbObjectList.h"
 #include "otbPolygon.h"
 */
@@ -96,10 +95,13 @@ public:
   typedef typename OutputImageType::PixelType   OutputPixelType;
   typedef typename OutputImageType::RegionType  OutputRegionType;
 
-   typedef TOutputMetricImage                          OutputMetricImageType;
-   typedef typename OutputMetricImageType::Pointer     OutputMetricImagePointerType;
-   typedef typename OutputMetricImageType::PixelType   OutputMetricPixelType;
-   typedef typename OutputMetricImageType::RegionType  MetricRegionType;
+  typedef TOutputMetricImage                          OutputMetricImageType;
+  typedef typename OutputMetricImageType::Pointer     OutputMetricImagePointerType;
+  typedef typename OutputMetricImageType::PixelType   OutputMetricPixelType;
+  typedef typename OutputMetricImageType::RegionType  MetricRegionType;
+  
+  typedef unsigned int                                OutputIterationPixelType;
+  typedef otb::Image<OutputIterationPixelType, 2>                OutputIterationImageType; 
 
   typedef TKernel      KernelType;
 
@@ -139,6 +141,9 @@ public:
   /** Return the mean shift vector image output */
   OutputMetricImageType * GetMetricOutput();
 
+  /** Returns the number of iterations done at each pixel */
+  OutputIterationImageType * GetIterationOutput();
+
   bool IsImageLatticeInitialized()
     {return m_ImageLatticeInitialized; };
 
diff --git a/Code/BasicFilters/otbMeanShiftImageFilter2.txx b/Code/BasicFilters/otbMeanShiftImageFilter2.txx
index f0b2c00d6c..16ce46f90b 100644
--- a/Code/BasicFilters/otbMeanShiftImageFilter2.txx
+++ b/Code/BasicFilters/otbMeanShiftImageFilter2.txx
@@ -62,10 +62,11 @@ MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>
 
 
   m_NeighborhoodHasTobeUpdated = true;
-  this->SetNumberOfOutputs(3);
+  this->SetNumberOfOutputs(4);
   this->SetNthOutput(0, OutputImageType::New());
   this->SetNthOutput(1, OutputImageType::New());
   this->SetNthOutput(2, OutputMetricImageType::New());
+  this->SetNthOutput(3, OutputIterationImageType::New());
 
 }
 
@@ -249,6 +250,7 @@ MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>
   return static_cast<const OutputMetricImageType *>(this->itk::ProcessObject::GetOutput(2));
 }
 
+
 template <class TInputImage, class TOutputMetricImage, class TOutputImage, class TKernel>
 typename MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>::OutputMetricImageType *
 MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>
@@ -261,6 +263,17 @@ MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>
   return static_cast<OutputMetricImageType *>(this->itk::ProcessObject::GetOutput(2));
 }
 
+template <class TInputImage,  class TOutputMetricImage,class TOutputImage, class TKernel>
+typename MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>::OutputIterationImageType *
+MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>
+::GetIterationOutput()
+{
+  if (this->GetNumberOfOutputs() < 4)
+    {
+      return 0;
+    }
+  return static_cast<OutputIterationImageType *>(this->itk::ProcessObject::GetOutput(3));
+}
 
 template <class TInputImage,class TOutputMetricImage, class TOutputImage, class TKernel>
 void
@@ -271,7 +284,7 @@ MeanShiftImageFilter2<TInputImage, TOutputMetricImage, TOutputImage, TKernel>
   typename OutputImageType::Pointer   spatialOutputPtr = this->GetSpatialOutput();
   typename OutputImageType::Pointer   rangeOutputPtr = this->GetRangeOutput();
   typename OutputImageType::Pointer   metricOutputPtr = this->GetMetricOutput();
-
+  typename OutputIterationImageType::Pointer iterationOutputPtr = this->GetIterationOutput();
 
   metricOutputPtr->SetBufferedRegion(metricOutputPtr->GetRequestedRegion());
   metricOutputPtr->Allocate();
@@ -282,6 +295,9 @@ MeanShiftImageFilter2<TInputImage, TOutputMetricImage, TOutputImage, TKernel>
   rangeOutputPtr->SetBufferedRegion(rangeOutputPtr->GetRequestedRegion());
   rangeOutputPtr->Allocate();
 
+  iterationOutputPtr->SetBufferedRegion(iterationOutputPtr->GetRequestedRegion());
+  iterationOutputPtr->Allocate();
+
  }
 
 
@@ -338,9 +354,10 @@ MeanShiftImageFilter2<TInputImage, TOutputMetricImage, TOutputImage, TKernel>
   TOutputMetricImage    * outMetricPtr = this->GetMetricOutput();
   TOutputImage * outSpatialPtr = this->GetSpatialOutput();
   TOutputImage * outRangePtr = this->GetRangeOutput();
+  OutputIterationImageType * outIterationPtr = this->GetIterationOutput();
 
   // Check pointers before using them
-  if(!inPtr || !outMetricPtr || !outSpatialPtr || !outRangePtr)
+  if(!inPtr || !outMetricPtr || !outSpatialPtr || !outRangePtr || !outIterationPtr)
     {
     return;
     }
@@ -440,61 +457,61 @@ typename MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKe
        // std::cout<<"it.Size "<<it->Size()<<std::endl;
         isInside=true;
         for(unsigned int comp=0; comp<rangeNumberOfComponents; comp++)
-                   {
-
-                     neighborhoodValue=it->GetElement(comp+spatialNumberOfComponents);
-                     if(vcl_abs(neighborhoodValue-rangePixel[comp])>m_SpectralBandwidth)
-                       isInside=false;
-                   }
-
-          if(it->GetElement(boundaryWeightIndex) && isInside)
-           {
-
+	  {
+	    
+	    neighborhoodValue=it->GetElement(comp+spatialNumberOfComponents);
+	    if(vcl_abs(neighborhoodValue-rangePixel[comp])>m_SpectralBandwidth)
+	      isInside=false;
+	  }
+ 
+	if(it->GetElement(boundaryWeightIndex) && isInside)
+	  {
+	    
             for(unsigned int comp=0; comp<spatialNumberOfComponents; comp++)
-            {
-             neighborhoodValue=it->GetElement(comp);
-             //value=spatialIt->GetElement(comp);
-             value=1;
-             meanShiftVector[comp]+=(neighborhoodValue);
-             weightingMeanShiftVector[comp]+=value;
-
-            }
-           for(unsigned int comp=0; comp<rangeNumberOfComponents; comp++)
-           {
-             neighborhoodValue=it->GetElement(comp+spatialNumberOfComponents);
-             //value=rangeIt->GetElement(comp);
-             value=1;
-            //  meanShiftVector[spatialNumberOfComponents+comp]+=(neighborhoodValue-rangePixel[comp])*(neighborhoodValue-rangePixel[comp])*neighborhoodValue*value;
-            //  weightingMeanShiftVector[spatialNumberOfComponents+comp]+=(neighborhoodValue-rangePixel[comp])*(neighborhoodValue-rangePixel[comp])*value;
-             meanShiftVector[spatialNumberOfComponents+comp]+=(neighborhoodValue);
-             weightingMeanShiftVector[spatialNumberOfComponents+comp]+=value;
-           }
-
-           }
-           ++it;
-           ++rangeIt;
-           ++spatialIt;
+	      {
+		neighborhoodValue=it->GetElement(comp);
+		//value=spatialIt->GetElement(comp);
+		value=1;
+		meanShiftVector[comp]+=(neighborhoodValue);
+		weightingMeanShiftVector[comp]+=value;
+
+	      }
+	    for(unsigned int comp=0; comp<rangeNumberOfComponents; comp++)
+	      {
+		neighborhoodValue=it->GetElement(comp+spatialNumberOfComponents);
+		//value=rangeIt->GetElement(comp);
+		value=1;
+		//  meanShiftVector[spatialNumberOfComponents+comp]+=(neighborhoodValue-rangePixel[comp])*(neighborhoodValue-rangePixel[comp])*neighborhoodValue*value;
+		//  weightingMeanShiftVector[spatialNumberOfComponents+comp]+=(neighborhoodValue-rangePixel[comp])*(neighborhoodValue-rangePixel[comp])*value;
+		meanShiftVector[spatialNumberOfComponents+comp]+=(neighborhoodValue);
+		weightingMeanShiftVector[spatialNumberOfComponents+comp]+=value;
+	      }
+	    
+	  }
+	++it;
+	++rangeIt;
+	++spatialIt;
         }
     }
-
+  
   for(unsigned int comp=0; comp<spatialNumberOfComponents; comp++)
-              {
-               if( weightingMeanShiftVector[comp]>0)
-                meanShiftVector[comp]=meanShiftVector[comp]/weightingMeanShiftVector[comp]-spatialPixel[comp];
-               else
-                meanShiftVector[comp]=0;
-              }
+    {
+      if( weightingMeanShiftVector[comp]>0)
+	meanShiftVector[comp]=meanShiftVector[comp]/weightingMeanShiftVector[comp]-spatialPixel[comp];
+      else
+	meanShiftVector[comp]=0;
+    }
   for(unsigned int comp=0; comp<rangeNumberOfComponents; comp++)
     {
-    if( weightingMeanShiftVector[spatialNumberOfComponents+comp]>0)
-                   meanShiftVector[spatialNumberOfComponents+comp]=meanShiftVector[spatialNumberOfComponents+comp]/weightingMeanShiftVector[spatialNumberOfComponents+comp]-rangePixel[comp];
-                  else
-                   meanShiftVector[spatialNumberOfComponents+comp]=0;
-   }
-
+      if( weightingMeanShiftVector[spatialNumberOfComponents+comp]>0)
+	meanShiftVector[spatialNumberOfComponents+comp]=meanShiftVector[spatialNumberOfComponents+comp]/weightingMeanShiftVector[spatialNumberOfComponents+comp]-rangePixel[comp];
+      else
+	meanShiftVector[spatialNumberOfComponents+comp]=0;
+    }
+  
   return meanShiftVector;
  }
-
+  
 
 // returns input spatial neighborhood, range, and binarry map for boundaries
 template <class TInputImage,class TOutputMetricImage, class TOutputImage, class TKernel>
@@ -606,6 +623,7 @@ MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>
   typename OutputImageType::Pointer spatialOutput = this->GetSpatialOutput();
   typename OutputImageType::Pointer rangeOutput = this->GetRangeOutput();
   typename OutputMetricImageType::Pointer metricOutput = this->GetMetricOutput();
+  typename OutputIterationImageType::Pointer iterationOutput = this->GetIterationOutput();
 
   typename InputImageType::ConstPointer input = this->GetInput();
 
@@ -614,10 +632,12 @@ MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>
   typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
   typedef itk::ImageRegionIterator<OutputMetricImageType> OutputMetricIteratorType;
   typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputIteratorWithIndexType;
+  typedef itk::ImageRegionIterator<OutputIterationImageType> OutputIterationIteratorType;
 
   typename OutputImageType::PixelType rangePixel;
   typename OutputImageType::PixelType spatialPixel;
   typename OutputMetricImageType::PixelType metricPixel;
+  typename OutputIterationImageType::PixelType iterationPixel;
 
   InputIteratorWithIndexType inputIt(input, inputRegionForThread);
 
@@ -628,12 +648,14 @@ MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>
   OutputIteratorType rangeIt(rangeOutput, outputRegionForThread);
   OutputIteratorType spatialIt(spatialOutput, outputRegionForThread);
   OutputMetricIteratorType metricIt(metricOutput, outputRegionForThread);
-
+  OutputIterationIteratorType iterationIt(iterationOutput, outputRegionForThread);
+  
   // fill pixel
   rangeIt.GoToBegin();
   spatialIt.GoToBegin();
   metricIt.GoToBegin();
-
+  iterationIt.GoToBegin();
+  
   unsigned int spatialNumberOfComponents = spatialOutput->GetNumberOfComponentsPerPixel();
   unsigned int rangeNumberOfComponents = rangeOutput->GetNumberOfComponentsPerPixel();
   unsigned int numberOfComponents = spatialNumberOfComponents + rangeNumberOfComponents;
@@ -711,11 +733,15 @@ MeanShiftImageFilter2<TInputImage,TOutputMetricImage, TOutputImage, TKernel>
     spatialIt.Set(spatialPixel);
     metricIt.Set(metricPixel);
 
+    iterationPixel = iteration;
+    iterationIt.Set(iterationPixel);
+
     ++inputIt;
 
     ++rangeIt;
     ++spatialIt;
     ++metricIt;
+    ++iterationIt;
 
     progress.CompletedPixel();
 
diff --git a/Code/IO/otbSpotImageMetadataInterface.cxx b/Code/IO/otbSpotImageMetadataInterface.cxx
index 5c00ea4b02..cdde8a29ea 100644
--- a/Code/IO/otbSpotImageMetadataInterface.cxx
+++ b/Code/IO/otbSpotImageMetadataInterface.cxx
@@ -173,7 +173,7 @@ SpotImageMetadataInterface::GetDay() const
   std::string valueString = imageKeywordlist.GetMetadataByKey("support_data.image_date");
   std::vector<std::string> outputValues;
 
-  boost::split(outputValues, valueString, boost::is_any_of(" T:-"));
+  boost::split(outputValues, valueString , boost::is_any_of(" T:-"));
 
   if (outputValues.size() <= 2) itkExceptionMacro(<< "Invalid Day");
 
diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt
index 456e9447a4..7457e9e445 100644
--- a/Testing/Code/BasicFilters/CMakeLists.txt
+++ b/Testing/Code/BasicFilters/CMakeLists.txt
@@ -1133,14 +1133,14 @@ ADD_TEST(bfTvMeanShiftImageFilterValid ${BASICFILTERS_TESTS9}
         ${TEMP}/bfMeanShiftImageFilterClusterBoundariesOutputValid.tif
         9 50 10 1.0
         )
-        
-       
+              
 ADD_TEST(bfTvMeanShiftImageFilter2 ${BASICFILTERS_TESTS9}
         otbMeanShiftImageFilter2
         ${INPUTDATA}/QB_Suburb.png
         ${TEMP}/bfMeanShiftImageSpatialOutput.tif
         ${TEMP}/bfMeanShiftImageFilterSpectralOutput.tif
         ${TEMP}/bfMeanShiftImageFilterMetricOutput.tif
+        ${TEMP}/bfMeanShiftImageFilterIterationOutput.tif
         9 9 50 10 0.00001 10
         )
         
@@ -1161,6 +1161,7 @@ ADD_TEST(bfTvMeanShiftImageFilter2Mul ${BASICFILTERS_TESTS9}
         ${TEMP}/bfMeanShiftImageSpatialOutputMul.tif
         ${TEMP}/bfMeanShiftImageFilterSpectralOutputMul.tif
         ${TEMP}/bfMeanShiftImageFilterMetricOutputMul.tif
+        ${TEMP}/bfMeanShiftImageFilterIterationOutputMul.tif
         9 9 50 10 0.01 100
         )
         
diff --git a/Testing/Code/BasicFilters/otbMeanShiftImageFilter2.cxx b/Testing/Code/BasicFilters/otbMeanShiftImageFilter2.cxx
index 7e4b0aecd4..77c3456bf9 100644
--- a/Testing/Code/BasicFilters/otbMeanShiftImageFilter2.cxx
+++ b/Testing/Code/BasicFilters/otbMeanShiftImageFilter2.cxx
@@ -24,7 +24,8 @@
 
 int otbMeanShiftImageFilter2(int argc, char * argv[])
 {
-  if (argc != 11)
+  std::cout << "coucou" << std::endl;
+  if (argc != 12)
     {
     std::cerr << "Usage: " << argv[0] <<
     " infname spatialfname spectralfname metricfname spatialRadius spectralRadius spectralbandwidth spatialBandwidth threshold maxiterationnumber"
@@ -36,12 +37,13 @@ int otbMeanShiftImageFilter2(int argc, char * argv[])
   const char *       spatialfname              = argv[2];
   const char *       spectralfname             = argv[3];
   const char *       metricfname               = argv[4];
-  const unsigned int spatialRadius             = atoi(argv[5]);
-  const unsigned int spectralRadius            = atoi(argv[6]);
-  const double       spectralbandwidth         = atof(argv[7]);
-  const double       spatialbandwidth          = atof(argv[8]);
-  const double       threshold                 = atof(argv[9]);
-  const unsigned int maxiterationnumber        = atoi(argv[10]);
+  const char *       iterationfname            = argv[5];
+  const unsigned int spatialRadius             = atoi(argv[6]);
+  const unsigned int spectralRadius            = atoi(argv[7]);
+  const double       spectralbandwidth         = atof(argv[8]);
+  const double       spatialbandwidth          = atof(argv[9]);
+  const double       threshold                 = atof(argv[10]);
+  const unsigned int maxiterationnumber        = atoi(argv[11]);
   /* maxit - threshold */
 
   const unsigned int Dimension = 2;
@@ -51,6 +53,8 @@ int otbMeanShiftImageFilter2(int argc, char * argv[])
   typedef otb::ImageFileReader<ImageType>                 ReaderType;
   typedef otb::ImageFileWriter<ImageType>                 WriterType;
   typedef otb::MeanShiftImageFilter2<ImageType, ImageType,ImageType,KernelType> FilterType;
+  typedef FilterType::OutputIterationImageType IterationImageType;
+  typedef otb::ImageFileWriter<IterationImageType> IterationWriterType;
 
   // Instantiating object
   FilterType::Pointer filter = FilterType::New();
@@ -77,18 +81,23 @@ int otbMeanShiftImageFilter2(int argc, char * argv[])
   WriterType::Pointer writer1 = WriterType::New();
   WriterType::Pointer writer2 = WriterType::New();
   WriterType::Pointer writer3 = WriterType::New();
+  IterationWriterType::Pointer writer4 = IterationWriterType::New();
 
 
   writer1->SetFileName(spatialfname);
   writer2->SetFileName(spectralfname);
   writer3->SetFileName(metricfname);
+  writer4->SetFileName(iterationfname);
 
   writer1->SetInput(filter->GetSpatialOutput());
   writer2->SetInput(filter->GetRangeOutput());
   writer3->SetInput(filter->GetMetricOutput());
+  writer4->SetInput(filter->GetIterationOutput());
+
   writer1->Update();
   writer2->Update();
   writer3->Update();
+  writer4->Update();
 
   return EXIT_SUCCESS;
 }
-- 
GitLab