TEST: vector image normalisation

Code/BasicFilters/otbHorizontalSobelVectorImageFilter.h

@@ -75,6 +75,9 @@ public:
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
 
+  /** Creation through object factory macro */
+  itkTypeMacro(HorizontalSobelVectorImageFilter, ImageToImageFilter);
+
 protected:
   HorizontalSobelVectorImageFilter()
   {

Code/BasicFilters/otbLocalActivityVectorImageFilter.h

@@ -97,6 +97,9 @@ public:
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
 
+  /** Creation through object factory macro */
+  itkTypeMacro(LocalActivityVectorImageFilter, ImageToImageFilter);
+
 protected:
   LocalActivityVectorImageFilter() { }
   virtual ~LocalActivityVectorImageFilter() { }

Code/BasicFilters/otbNormalizeVectorImageFilter.h

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbNormalizeVectorImageFilter_h
+#define __otbNormalizeVectorImageFilter_h
+
+#include <vnl/vnl_math.h>
+
+#include <itkNumericTraits.h>
+#include <itkExceptionObject.h>
+#include <itkVariableLengthVector.h>
+#include <otbVectorImage.h>
+#include <otbUnaryFunctorVectorImageFilter.h>
+
+namespace otb 
+{
+namespace Functor
+{
+/** \brief NormalizeVectorImageFunctor
+ * \brief this functor performs affine transformation on VariableLengthVectors
+ */
+template < class TInput, class TOutput >
+class NormalizeVectorImageFunctor
+{
+public:
+  NormalizeVectorImageFunctor() { }
+  virtual ~NormalizeVectorImageFunctor () { }
+
+  typedef typename itk::NumericTraits< TInput >::RealType RealVectorType;
+  typedef typename  itk::NumericTraits< typename RealVectorType::ValueType >::RealType RealType;
+
+  TOutput operator() ( const TInput & input )
+  {
+    unsigned int length = input.Size();
+    TOutput output ( length );
+    for ( unsigned int i = 0; i < length; i++ )
+    {
+      output[i] = static_cast<typename TOutput::ValueType>(
+                    ( static_cast< RealType >( input[i] ) - m_Mean[i] ) 
+                      / m_StdDev[i] );
+    }
+    return output;
+  }
+
+  template < class T > 
+  void SetMean ( const itk::VariableLengthVector<T> & m )
+  {
+    m_Mean.SetSize( m.Size() );
+    for ( unsigned int i = 0; i < m_Mean.Size(); i++ )
+      m_Mean[i] = static_cast< RealType >( m[i] );
+  }
+
+  template < class T>
+  void SetStdDev ( const itk::VariableLengthVector<T> & sigma )
+  {
+    m_StdDev.SetSize( sigma.Size() );
+    for ( unsigned int i = 0; i < m_StdDev.Size(); i++ )
+    {
+      m_StdDev[i] = static_cast< RealType >( sigma[i] );
+      if ( m_StdDev[i] == itk::NumericTraits< RealType >::Zero )
+      {
+        throw itk::ExceptionObject(__FILE__, __LINE__, 
+                "Cannot divide by zero !", ITK_LOCATION);
+      }
+    }
+  }
+
+  template < class T >
+  void SetVariance ( const itk::VariableLengthVector<T> & var )
+  {
+    m_StdDev.SetSize( var.Size() );
+    for ( unsigned int i = 0; i < m_StdDev.Size(); i++ )
+    {
+      m_StdDev[i] = vcl_sqrt( static_cast< RealType >( var[i] ) );
+      if ( m_StdDev[i] == itk::NumericTraits< RealType >::Zero )
+      {
+        throw itk::ExceptionObject(__FILE__, __LINE__, 
+                "Cannot divide by zero !", ITK_LOCATION);
+      }
+    }
+  }
+
+protected:
+
+  RealVectorType m_Mean;
+  RealVectorType m_StdDev;
+}; // end of class NormalizeVectorImageFunctor
+
+} // end of namespace Functor
+
+
+/** \class NormalizeVectorImageFilter
+ * \brief Normalize an VectorImage by setting its mean to zero and possibily variance to one (band by band).
+ * 
+ * NormalizeVectorImageFilter shifts and scales an image so that the pixels in the image 
+ * have a zero mean and unit variance. 
+ *
+ * This filter uses StreamingStatisticsVectorImageFilter2 to compute the mean and variance of
+ * the input and then applies the dedicated functor.
+ *
+ * \sa StreamingStatisticsVectorImageFilter2
+ */
+template < class TInputImage, class TOutputImage >
+class ITK_EXPORT NormalizeVectorImageFilter 
+  : public UnaryFunctorVectorImageFilter< TInputImage, TOutputImage,
+            Functor::NormalizeVectorImageFunctor< 
+              typename TInputImage::PixelType,
+              typename TOutputImage::PixelType > >
+{
+public:
+  /** Standart class typedefs */
+  typedef NormalizeVectorImageFilter Self;
+  typedef UnaryFunctorVectorImageFilter< TInputImage, TOutputImage,
+            Functor::NormalizeVectorImageFunctor< 
+              typename TInputImage::PixelType,
+              typename TOutputImage::PixelType > > Superclass;
+  typedef itk::SmartPointer<Self> Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Creation through object factory macro */
+  itkTypeMacro(NormalizeVectorImageFilter, ImageToImageFilter);
+
+  typedef TInputImage InputImageType;
+  typedef TOutputImage OutputImageType;
+
+  typedef typename itk::NumericTraits< typename TInputImage::PixelType >::RealType RealVectorType;
+  typedef typename  itk::NumericTraits< typename RealVectorType::ValueType >::RealType RealType;
+
+  template < class T >
+  void SetMean ( const itk::VariableLengthVector<T> & m )
+  {
+    this->m_Functor->SetMean( m );
+    m_IsGivenMean = true;
+    m_UseMean = true;
+    this->Modified();
+  }
+
+  template < class T >
+  void SetStdDev ( const itk::VariableLengthVector<T> & sigma )
+  {
+    this->m_Functor->SetStdDev( sigma );
+    m_IsGivenStdDev = true;
+    m_UseStdDev = true;
+    this->Modified();
+  }
+
+  template < class T >
+  void SetVariance ( const itk::VariableLengthVector<T> & var )
+  {
+    this->m_Functor->SetVariance( var );
+    m_IsGivenStdDev = true;
+    m_UseStdDev = true;
+    this->Modified();
+  }
+
+
+  itkSetMacro(UseMean,bool);
+  itkSetMacro(UseStdDev,bool);
+
+protected:
+  NormalizeVectorImageFilter ();
+  virtual ~NormalizeVectorImageFilter() { }
+
+  void BeforeThreadedGenerateData();
+
+private:
+  NormalizeVectorImageFilter ( const Self & );
+  void operator=( const Self & );
+
+  bool m_IsGivenMean;
+  bool m_IsGivenStdDev;
+
+  bool m_UseMean;
+  bool m_UseStdDev;
+}; // end of class NormalizeVectorImageFilter
+
+} // end of namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbNormalizeVectorImageFilter.txx"
+#endif
+
+#endif // __otbNormalizeVectorImageFilter_h
+

Code/BasicFilters/otbNormalizeVectorImageFilter.txx

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbNormalizeVectorImageFilter_txx
+#define __otbNormalizeVectorImageFilter_txx
+#include "otbNormalizeVectorImageFilter.h"
+
+#include "otbStreamingStatisticsVectorImageFilter2.h"
+
+namespace otb 
+{
+
+template < class TInputImage, class TOutputImage >
+NormalizeVectorImageFilter< TInputImage, TOutputImage >
+::NormalizeVectorImageFilter ()
+{
+  m_IsGivenMean = false;
+  m_IsGivenStdDev = false;
+
+  m_UseMean = true;
+  m_UseStdDev = true;
+}
+
+template < class TInputImage, class TOutputImage >
+void
+NormalizeVectorImageFilter< TInputImage, TOutputImage >
+::BeforeThreadedGenerateData ()
+{
+  if ( !m_UseMean )
+  {
+    typename TInputImage::PixelType vector ( this->GetInput()->GetNumberOfComponentsPerPixel() );
+    vector.Fill( itk::NumericTraits< typename TInputImage::PixelType::ValueType >::Zero );
+    this->GetFunctor().SetMean( vector );
+  }
+
+  if ( !m_UseStdDev )
+  {
+    typename TInputImage::PixelType vector ( this->GetInput()->GetNumberOfComponentsPerPixel() );
+    vector.Fill( itk::NumericTraits< typename TInputImage::PixelType::ValueType >::One );
+    this->GetFunctor().SetStdDev( vector );
+  }
+
+  if ( !m_IsGivenMean )
+  {
+    typedef StreamingStatisticsVectorImageFilter2< InputImageType > CovarianceEstimatorFilterType;
+    typename CovarianceEstimatorFilterType::Pointer estimator = CovarianceEstimatorFilterType::New();
+    estimator->SetInput( const_cast<InputImageType*>( this->GetInput() ) );
+    estimator->Update();
+
+    this->GetFunctor().SetMean( estimator->GetMean() );
+
+    if ( !m_IsGivenStdDev && m_UseStdDev )
+    {
+      typename StreamingStatisticsVectorImageFilter2< InputImageType >::RealPixelType sigma
+        ( this->GetInput()->GetNumberOfComponentsPerPixel() );
+      for ( unsigned int i = 0; i < sigma.Size(); i++ )
+        sigma[i] = vcl_sqrt( estimator->GetCovariance()(i,i) );
+
+      this->GetFunctor().SetStdDev( sigma );
+    }
+  }
+}
+
+} // end of namespace otb
+
+#endif // __otbNormalizeVectorImageFilter_txx

Code/BasicFilters/otbSobelVectorImageFilter.h

@@ -84,6 +84,9 @@ public:
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
 
+  /** Creation through object factory macro */
+  itkTypeMacro(SobelVectorImageFilter, ImageToImageFilter);
+
 protected:
   SobelVectorImageFilter()
   {

Code/BasicFilters/otbUnaryFunctorVectorImageFilter.txx

@@ -75,7 +75,7 @@ UnaryFunctorVectorImageFilter<TInputImage, TOutputImage, TFunction>
 
   while ( !outputIt.IsAtEnd() && !inputIt.IsAtEnd() )
   {
-    outputIt.Set( m_Functor( neighInputIt ) );
+    outputIt.Set( m_Functor( inputIt.Get() ) );
 
     ++inputIt;
     ++outputIt;

Code/BasicFilters/otbVerticalSobelVectorImageFilter.h

@@ -75,6 +75,9 @@ public:
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
 
+  /** Creation through object factory macro */
+  itkTypeMacro(VerticalSobelVectorImageFilter, ImageToImageFilter);
+
 protected:
   VerticalSobelVectorImageFilter()
   {

Testing/CMakeLists.txt

@@ -20,6 +20,7 @@ SET( otbHyperTests1_SRC
   otbVerticalSobelVectorImageFilter.cxx
   otbSobelVectorImageFilter.cxx
   otbLocalActivityVectorImageFilter.cxx
+  otbNormalizeVectorImageFilter.cxx
   otbPCAImageFilter.cxx
   otbMNFImageFilter.cxx
 )
@@ -200,6 +201,16 @@ ADD_TEST(bfTvLocalActivityVectorImageFilter
 		 -in ${DATA}/CupriteSubHsi/cupriteSubHsi.tif 
 		 -out ${TEMP}/cupriteLocalActivityFlt.hdr)
 
+ADD_TEST(bfTuNormalizeVectorImageFilter
+         ${TESTEXE_DIR}/otbHyperTests1
+		 otbNormalizeVectorImageFilterNewTest)
+
+ADD_TEST(bfTvNormalizeVectorImageFilter
+         ${TESTEXE_DIR}/otbHyperTests1
+		 otbNormalizeVectorImageFilterTest
+         -in ${DATA}/CupriteSubHsi/cupriteSubHsi.tif 
+		 -out ${TEMP}/cupriteNormalizedFlt.hdr)
+
 ADD_TEST(hyTuPCAImageFilterNew
          ${TESTEXE_DIR}/otbHyperTests1
 		 otbPCAImageFilterNewTest)

Testing/otbHyperTests1.cxx

@@ -51,6 +51,8 @@ void RegisterTests()
   REGISTER_TEST(otbSobelVectorImageFilterTest);
   REGISTER_TEST(otbLocalActivityVectorImageFilterNewTest);
   REGISTER_TEST(otbLocalActivityVectorImageFilterTest);
+  REGISTER_TEST(otbNormalizeVectorImageFilterNewTest);
+  REGISTER_TEST(otbNormalizeVectorImageFilterTest);
   REGISTER_TEST(otbPCAImageFilterNewTest);
   REGISTER_TEST(otbPCAImageFilterTest);
   REGISTER_TEST(otbMNFImageFilterNewTest);

Testing/otbNormalizeVectorImageFilter.cxx

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbCommandProgressUpdate.h"
+#include "otbCommandLineArgumentParser.h"
+
+#include "otbNormalizeVectorImageFilter.h"
+
+int otbNormalizeVectorImageFilterNewTest ( int argc, char* argv[] )
+{
+  const unsigned int Dimension = 2;
+  typedef double PixelType;
+  typedef otb::VectorImage< PixelType, Dimension > ImageType;
+  typedef otb::NormalizeVectorImageFilter< ImageType, ImageType >
+    FilterType;
+  FilterType::Pointer filter = FilterType::New();
+
+  return EXIT_SUCCESS;
+}
+
+int otbNormalizeVectorImageFilterTest ( int argc, char* argv[] )
+{
+  typedef otb::CommandLineArgumentParser ParserType;
+  ParserType::Pointer parser = ParserType::New();
+
+  parser->AddInputImage();
+  parser->AddOutputImage();
+
+  typedef otb::CommandLineArgumentParseResult ParserResultType;  
+  ParserResultType::Pointer  parseResult = ParserResultType::New();
+    
+  try
+  {  
+    parser->ParseCommandLine( argc, argv, parseResult );
+  }
+  catch( itk::ExceptionObject & err )
+  {
+    std::cerr << argv[0] << " performs horizonal sobel on a vector image\n";
+    std::string descriptionException = err.GetDescription();
+    if ( descriptionException.find("ParseCommandLine(): Help Parser")
+        != std::string::npos )
+      return EXIT_SUCCESS;
+    if(descriptionException.find("ParseCommandLine(): Version Parser")
+        != std::string::npos )
+      return EXIT_SUCCESS;
+    return EXIT_FAILURE;
+  }
+
+  const char * inputImageName = parseResult->GetInputImage().c_str();
+  const char * outputImageName = parseResult->GetOutputImage().c_str();
+
+  // Main type definition
+  const unsigned int Dimension = 2;
+  typedef double PixelType;
+  typedef otb::VectorImage< PixelType, Dimension > ImageType;
+
+  // Reading input images
+  typedef otb::ImageFileReader<ImageType> ReaderType;
+  ReaderType::Pointer reader = ReaderType::New();
+  reader->SetFileName(inputImageName);
+
+  // Image filtering
+  typedef otb::NormalizeVectorImageFilter< ImageType, ImageType >
+    FilterType;
+  FilterType::Pointer filter = FilterType::New();
+  filter->SetInput( reader->GetOutput() );
+
+  typedef otb::CommandProgressUpdate< FilterType > CommandType;
+  CommandType::Pointer observer = CommandType::New();
+  filter->AddObserver( itk::ProgressEvent(), observer );
+
+  typedef otb::ImageFileWriter< ImageType > ImageWriterType;
+  ImageWriterType::Pointer writer = ImageWriterType::New();
+  writer->SetFileName( outputImageName );
+  writer->SetInput( filter->GetOutput() );
+  writer->Update();
+
+  return EXIT_SUCCESS;
+}
+
+
+