Add simple utility tools

CMakeLists.txt

+project(SimpleExtractionTools)
+otb_module_impl()

app/CMakeLists.txt

+cmake_minimum_required (VERSION 2.8)
+
+OTB_CREATE_APPLICATION(NAME           ExtractBand
+                       SOURCES        otbExtractBand.cxx
+                       LINK_LIBRARIES OTBCommon)
+
+OTB_CREATE_APPLICATION(NAME           ExtractGeom
+                       SOURCES        otbExtractGeom.cxx
+                       LINK_LIBRARIES OTBCommon)
+

app/otbExtractBand.cxx

+/*=========================================================================
+
+  Copyright (c) IRSTEA. All rights reserved.
+
+
+     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 "itkFixedArray.h"
+#include "itkObjectFactory.h"
+
+// Elevation handler
+#include "otbWrapperElevationParametersHandler.h"
+#include "otbWrapperApplicationFactory.h"
+
+// Application engine
+#include "otbStandardFilterWatcher.h"
+#include "itkFixedArray.h"
+
+// Filter
+#include "otbMultiChannelExtractROI.h"
+
+using namespace std;
+
+namespace otb
+{
+
+namespace Wrapper
+{
+
+class ExtractBand : public Application
+{
+public:
+  /** Standard class typedefs. */
+  typedef ExtractBand                        Self;
+  typedef Application                   Superclass;
+  typedef itk::SmartPointer<Self>       Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  typedef otb::MultiChannelExtractROI<FloatVectorImageType::InternalPixelType,
+      FloatVectorImageType::InternalPixelType> ExtractFilterType;
+
+  /** Standard macro */
+  itkNewMacro(Self);
+  itkTypeMacro(ExtractBand, Application);
+
+  void DoInit()
+  {
+
+    SetName("ExtractBand");
+    SetDescription("Perform single band extraction");
+
+    // Documentation
+    SetDocName("ExtractBand");
+    SetDocLongDescription("This application performs single band extraction");
+    SetDocLimitations("None");
+    SetDocAuthors("Remi Cresson");
+    SetDocSeeAlso(" ");
+
+    AddDocTag(Tags::Manip);
+
+    AddParameter(ParameterType_InputImage,   "inxs",  "Input XS Image");
+    SetParameterDescription("inxs"," Input XS image.");
+
+    AddParameter(ParameterType_Int, "band", "band index" );
+    SetParameterDescription("band","Set the band index to extract" );
+    SetMinimumParameterIntValue("band", 1);
+
+    AddParameter(ParameterType_OutputImage,  "out",   "Output image");
+    SetParameterDescription("out"," Output image.");
+
+    AddRAMParameter();
+
+    // Doc example parameter settings
+    SetDocExampleParameterValue("band", "2");
+    SetDocExampleParameterValue("inxs", "QB_Toulouse_Ortho_XS.tif");
+    SetDocExampleParameterValue("out", "Pansharpening.tif uint16");
+
+
+  }
+
+  void DoUpdateParameters()
+  {
+    // Nothing to do here : all parameters are independent
+  }
+
+  void DoExecute()
+  {
+
+    FloatVectorImageType* xs = GetParameterImage("inxs");
+
+    unsigned int band = GetParameterInt("band");
+    m_Filter = ExtractFilterType::New();
+    m_Filter->SetFirstChannel(band);
+    m_Filter->SetLastChannel(band);
+    m_Filter->SetInput(xs);
+
+    SetParameterOutputImage("out", m_Filter->GetOutput());
+
+  }
+
+  ExtractFilterType::Pointer m_Filter;
+
+};
+}
+}
+
+OTB_APPLICATION_EXPORT( otb::Wrapper::ExtractBand )

app/otbExtractGeom.cxx

+/*=========================================================================
+
+  Copyright (c) IRSTEA. All rights reserved.
+
+
+     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 "itkFixedArray.h"
+#include "itkObjectFactory.h"
+
+// Elevation handler
+#include "otbWrapperElevationParametersHandler.h"
+#include "otbWrapperApplicationFactory.h"
+
+// Application engine
+#include "otbStandardFilterWatcher.h"
+#include "itkFixedArray.h"
+
+// Filter
+#include "otbMultiChannelExtractROI.h"
+#include "otbVectorDataToLabelImageCustomFilter.h"
+#include "itkMaskImageFilter.h"
+#include "otbVectorDataIntoImageProjectionFilter.h"
+#include "otbVectorDataProperties.h"
+#include "otbRegionComparator.h"
+
+using namespace std;
+
+namespace otb
+{
+
+namespace Wrapper
+{
+
+class ExtractGeom : public Application
+{
+public:
+  /** Standard class typedefs. */
+  typedef ExtractGeom                        Self;
+  typedef Application                   Superclass;
+  typedef itk::SmartPointer<Self>       Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  /** Standard macro */
+  itkNewMacro(Self);
+  itkTypeMacro(ExtractGeom, Application);
+
+  /** Filters */
+  typedef otb::MultiChannelExtractROI<FloatVectorImageType::InternalPixelType,
+      FloatVectorImageType::InternalPixelType> ExtractFilterType;
+
+  /** mask */
+  typedef bool TMaskPixelValueType;
+  typedef otb::Image<TMaskPixelValueType, 2> MaskImageType;
+  typedef itk::MaskImageFilter<FloatVectorImageType, MaskImageType,
+      FloatVectorImageType> MaskFilterType;
+
+  /* vector data filters */
+  typedef otb::VectorDataIntoImageProjectionFilter<VectorDataType,
+      FloatVectorImageType> VectorDataReprojFilterType;
+  typedef otb::VectorDataToLabelImageCustomFilter<VectorDataType,
+      MaskImageType> RasteriseFilterType;
+  typedef otb::VectorDataProperties<VectorDataType> VectorDataPropertiesType;
+
+  void DoInit()
+  {
+    SetName("ExtractGeom");
+    SetDescription("Perform geometric extraction");
+
+    // Documentation
+    SetDocName("ExtractGeom");
+    SetDocLongDescription("This application performs geometric extraction");
+    SetDocLimitations("None");
+    SetDocAuthors("RemiCresson");
+    SetDocSeeAlso(" ");
+
+    AddDocTag(Tags::Manip);
+
+    AddParameter(ParameterType_InputImage,   "in",  "Input Image");
+    SetParameterDescription("in"," Input image.");
+
+    AddParameter(ParameterType_InputVectorData, "shp", "Input Shapefile" );
+    SetParameterDescription("shp","Input Shapefile" );
+
+    AddParameter(ParameterType_OutputImage,  "out",   "Output image");
+    SetParameterDescription("out"," Output image.");
+
+    AddRAMParameter();
+
+    // Doc example parameter settings
+    SetDocExampleParameterValue("shp", "vecteur.shp");
+    SetDocExampleParameterValue("in", "QB_Toulouse_Ortho_XS.tif");
+    SetDocExampleParameterValue("out", "QB_Toulouse_Ortho_XS_decoup.tif uint16");
+
+
+
+
+  }
+
+  void DoUpdateParameters()
+  {
+    // Nothing to do here : all parameters are independent
+  }
+
+  void DoExecute()
+  {
+
+    // Get inputs
+    FloatVectorImageType::Pointer xs = GetParameterImage("in");
+    VectorDataType* shp = GetParameterVectorData("shp");
+
+    /* Reproject vector data */
+    vdReproj = VectorDataReprojFilterType::New();
+    vdReproj->SetInputVectorData(shp);
+    vdReproj->SetInputImage(xs);
+    vdReproj->Update();
+
+    /* Get VectorData bounding box */
+    vdProperties = VectorDataPropertiesType::New();
+    vdProperties->SetVectorDataObject(vdReproj->GetOutput());
+    vdProperties->ComputeBoundingRegion();
+
+    /* Compute intersecting region */
+    otb::RegionComparator<FloatVectorImageType, FloatVectorImageType> comparator;
+    comparator.SetImage1(xs);
+    FloatVectorImageType::RegionType roi =
+        comparator.RSRegionToImageRegion(vdProperties->GetBoundingRegion());
+    roi.PadByRadius(1); // avoid extrapolation
+    if (!roi.Crop(xs->GetLargestPossibleRegion()))
+      {
+      otbAppLogFATAL( << " Input VectorData is outside image !" );
+      return;
+      }
+
+    /* Rasterize vector data */
+    rasterizer = RasteriseFilterType::New();
+    rasterizer->AddVectorData(vdReproj->GetOutput());
+    rasterizer->SetOutputOrigin(xs->GetOrigin());
+    rasterizer->SetOutputSpacing(xs->GetSpacing());
+    rasterizer->SetOutputSize(
+        xs->GetLargestPossibleRegion().GetSize());
+    rasterizer->SetBurnMaxValueMode(true);
+    rasterizer->SetOutputProjectionRef(
+        xs->GetProjectionRef());
+    rasterizer->Update();
+
+    /* Mask input image */
+    maskFilter = MaskFilterType::New();
+    maskFilter->SetInput(xs);
+    maskFilter->SetMaskImage(rasterizer->GetOutput());
+
+    /* Extract ROI */
+    m_Filter = ExtractFilterType::New();
+    m_Filter->SetInput(maskFilter->GetOutput());
+    m_Filter->SetExtractionRegion(roi);
+
+    SetParameterOutputImage("out", m_Filter->GetOutput());
+
+
+  }
+
+  ExtractFilterType::Pointer m_Filter;
+  VectorDataReprojFilterType::Pointer vdReproj;
+  RasteriseFilterType::Pointer rasterizer;
+  VectorDataPropertiesType::Pointer vdProperties;
+  MaskFilterType::Pointer maskFilter;
+
+};
+}
+}
+
+OTB_APPLICATION_EXPORT( otb::Wrapper::ExtractGeom )

include/otbRegionComparator.h

+/*=========================================================================
+
+  Copyright (c) IRSTEA. All rights reserved.
+
+
+     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 otbRegionComparator_H_
+#define otbRegionComparator_H_
+
+#include "otbObjectList.h"
+#include "itkContinuousIndex.h"
+#include "itkImageLinearConstIteratorWithIndex.h"
+#include "itkMetaDataDictionary.h"
+#include "ogr_api.h"
+#include "ogr_feature.h"
+#include "otbMetaDataKey.h"
+#include "itkMetaDataObject.h"
+#include "itkNumericTraits.h"
+
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <string>
+#include <algorithm>
+
+#include "otbRemoteSensingRegion.h"
+
+#define msg(x) do { std::cout << x << std::endl;} while(0)
+
+using namespace std;
+
+namespace otb {
+
+template<class TInputImage1, class TInputImage2=TInputImage1> class RegionComparator{
+public:
+
+  typedef TInputImage1    InputImage1Type;
+  typedef typename InputImage1Type::Pointer   InputImage1PointerType;
+  typedef typename InputImage1Type::RegionType    InputImage1RegionType;
+
+  typedef TInputImage2    InputImage2Type;
+  typedef typename InputImage2Type::Pointer   InputImage2PointerType;
+  typedef typename InputImage2Type::RegionType    InputImage2RegionType;
+
+  typedef typename otb::RemoteSensingRegion<double> RSRegion;
+
+  void SetImage1(InputImage1PointerType im1) {m_InputImage1 = im1;}
+  void SetImage2(InputImage2PointerType im2) {m_InputImage2 = im2;}
+
+  InputImage2RegionType GetOverlapInImage1Indices()
+  {
+    // Largest region of im2 --> region in im1
+    InputImage1RegionType region2 = m_InputImage2->GetLargestPossibleRegion();
+    InputImage1RegionType region2_in_1;
+    OutputRegionToInputRegion(region2, region2_in_1, m_InputImage2, m_InputImage1);
+
+    // Collision test
+    InputImage2RegionType region1 = m_InputImage1->GetLargestPossibleRegion();
+    region2_in_1.Crop(region1);
+
+    // Return overlap
+    return region2_in_1;
+  }
+
+  InputImage2RegionType GetOverlapInImage2Indices()
+  {
+    // Largest region of im1 --> region in im2
+    InputImage1RegionType region1 = m_InputImage1->GetLargestPossibleRegion();
+    InputImage1RegionType region1_in_2;
+    OutputRegionToInputRegion(region1, region1_in_2, m_InputImage1, m_InputImage2);
+
+    // Collision test
+    InputImage2RegionType region2 = m_InputImage2->GetLargestPossibleRegion();
+    region1_in_2.Crop(region2);
+
+    // Return overlap
+    return region1_in_2;
+  }
+
+  bool DoesOverlap()
+  {
+    // Largest region of im1 --> region in im2
+    InputImage1RegionType region1 = m_InputImage1->GetLargestPossibleRegion();
+    InputImage1RegionType region1_in_2;
+    OutputRegionToInputRegion(region1, region1_in_2, m_InputImage1, m_InputImage2);
+
+    // Collision test
+    InputImage2RegionType region2 = m_InputImage2->GetLargestPossibleRegion();
+    return region1_in_2.Crop(region2);
+  }
+
+  /*
+   * Converts sourceRegion of im1 into targetRegion of im2
+   */
+  void OutputRegionToInputRegion(const InputImage1RegionType &sourceRegion, InputImage2RegionType &targetRegion,
+      const InputImage1PointerType &sourceImage, const InputImage2PointerType &targetImage)
+  {
+
+    // Source Region (source image indices)
+    typename InputImage1RegionType::IndexType sourceRegionIndexStart = sourceRegion.GetIndex();
+    typename InputImage1RegionType::IndexType sourceRegionIndexEnd;
+    for(unsigned int dim = 0; dim < InputImage1Type::ImageDimension; ++dim)
+      sourceRegionIndexEnd[dim]= sourceRegionIndexStart[dim] + sourceRegion.GetSize()[dim]-1;
+
+    // Source Region (Geo)
+    typename InputImage1Type::PointType sourceRegionIndexStartGeo, sourceRegionIndexEndGeo;
+    sourceImage->TransformIndexToPhysicalPoint(sourceRegionIndexStart, sourceRegionIndexStartGeo);
+    sourceImage->TransformIndexToPhysicalPoint(sourceRegionIndexEnd  , sourceRegionIndexEndGeo  );
+
+    // Source Region (target image indices)
+    typename InputImage2RegionType::IndexType targetIndexStart, targetIndexEnd;
+    targetImage->TransformPhysicalPointToIndex(sourceRegionIndexStartGeo, targetIndexStart);
+    targetImage->TransformPhysicalPointToIndex(sourceRegionIndexEndGeo  , targetIndexEnd);
+
+    // Target Region (target image indices)
+    typename InputImage2RegionType::IndexType targetRegionStart;
+    typename InputImage2RegionType::SizeType targetRegionSize;
+    for(unsigned int dim = 0; dim<InputImage2Type::ImageDimension; ++dim)
+      {
+      targetRegionStart[dim] = std::min(targetIndexStart[dim], targetIndexEnd[dim]);
+      targetRegionSize[dim]  = std::max(targetIndexStart[dim], targetIndexEnd[dim]) - targetRegionStart[dim] + 1;
+      }
+    InputImage2RegionType computedInputRegion(targetRegionStart, targetRegionSize);
+
+    // Avoid extrapolation
+    computedInputRegion.PadByRadius(1);
+
+    // Target Region
+    targetRegion = computedInputRegion;
+
+  }
+
+  /*
+   * Converts a RemoteSensingREgion into a ImageRegion (image1)
+   */
+  InputImage1RegionType RSRegionToImageRegion(const RSRegion rsRegion)
+  {
+    typename itk::ContinuousIndex<double> rsRegionOrigin = rsRegion.GetOrigin();
+    typename itk::ContinuousIndex<double> rsRegionSize = rsRegion.GetSize();
+    typename itk::ContinuousIndex<double> rsRegionEnd;
+    for(unsigned int dim = 0; dim<InputImage2Type::ImageDimension; ++dim)
+      {
+      rsRegionEnd[dim] = rsRegionOrigin[dim] + rsRegionSize[dim];
+      }
+    typename InputImage1RegionType::IndexType imageRegionStart, imageRegionEnd;
+    m_InputImage1->TransformPhysicalPointToIndex(rsRegionOrigin, imageRegionStart);
+    m_InputImage1->TransformPhysicalPointToIndex(rsRegionEnd, imageRegionEnd);
+
+    // Target Region (target image indices)
+    typename InputImage1RegionType::IndexType targetRegionStart;
+    typename InputImage1RegionType::SizeType targetRegionSize;
+    for(unsigned int dim = 0; dim<InputImage2Type::ImageDimension; ++dim)
+      {
+      targetRegionStart[dim] = std::min(imageRegionStart[dim], imageRegionEnd[dim]);
+      targetRegionSize[dim]  = std::max(imageRegionStart[dim], imageRegionEnd[dim]) - targetRegionStart[dim] + 1;
+      }
+
+    InputImage1RegionType region(targetRegionStart, targetRegionSize);
+    return region;
+
+  }
+
+  /*
+   * Converts a given region of an image into a remoste sensing region
+   */
+  RSRegion ImageRegionToRSRegion(const InputImage1RegionType &sourceRegion, InputImage1PointerType sourceImage)
+  {
+
+    // Source Region (source image indices)
+    typename InputImage1RegionType::IndexType sourceRegionIndexStart = sourceRegion.GetIndex();
+    typename InputImage1RegionType::IndexType sourceRegionIndexEnd;
+    for(unsigned int dim = 0; dim < InputImage1Type::ImageDimension; ++dim)
+      sourceRegionIndexEnd[dim]= sourceRegionIndexStart[dim] + sourceRegion.GetSize()[dim]-1;
+
+    // Source Region (Geo)
+    typename InputImage1Type::PointType sourceRegionIndexStartGeo, sourceRegionIndexEndGeo;
+    sourceImage->TransformIndexToPhysicalPoint(sourceRegionIndexStart, sourceRegionIndexStartGeo);
+    sourceImage->TransformIndexToPhysicalPoint(sourceRegionIndexEnd  , sourceRegionIndexEndGeo  );
+
+    // Source Region (Geo min & max)
+    typename itk::ContinuousIndex<double> sourceRegionMin, sourceRegionMax, sourceRegionSize;
+    for(unsigned int dim = 0; dim<InputImage2Type::ImageDimension; ++dim)
+      {
+      sourceRegionMin[dim] = std::min(sourceRegionIndexStartGeo[dim], sourceRegionIndexEndGeo[dim]);
+      sourceRegionMax[dim] = std::max(sourceRegionIndexStartGeo[dim], sourceRegionIndexEndGeo[dim]);
+      }
+    RSRegion region;
+    sourceRegionSize[0] = sourceRegionMax[0] - sourceRegionMin[0];
+    sourceRegionSize[1] = sourceRegionMax[1] - sourceRegionMin[1];
+    region.SetOrigin(sourceRegionMin);
+    region.SetSize(sourceRegionSize);
+    return region;
+
+  }
+  bool HaveSameProjection()
+  {
+    itk::MetaDataDictionary metaData1 = m_InputImage1 -> GetMetaDataDictionary();
+    itk::MetaDataDictionary metaData2 = m_InputImage2 -> GetMetaDataDictionary();
+
+    if (metaData1.HasKey(otb::MetaDataKey::ProjectionRefKey))
+      {
+      if (metaData2.HasKey(otb::MetaDataKey::ProjectionRefKey))
+        {
+        string proj1, proj2;
+        itk::ExposeMetaData<string>(metaData1, static_cast<string>(otb::MetaDataKey::ProjectionRefKey), proj1);
+        itk::ExposeMetaData<string>(metaData2, static_cast<string>(otb::MetaDataKey::ProjectionRefKey), proj2);
+        if (proj1.compare(proj2)==0)
+          return true;
+        }
+      }
+    return false;
+
+  }
+
+  bool HaveSameNbOfBands()
+  {
+    return ( (m_InputImage1->GetNumberOfComponentsPerPixel()) == (m_InputImage2->GetNumberOfComponentsPerPixel()) );
+  }
+
+  RSRegion ComputeIntersectingRemoteSensingRegion()
+  {
+    InputImage2RegionType region1 = m_InputImage1->GetLargestPossibleRegion();
+    InputImage1RegionType region2 = m_InputImage2->GetLargestPossibleRegion();
+    RSRegion rsr1 = ImageRegionToRSRegion(region1, m_InputImage1);
+    RSRegion rsr2 = ImageRegionToRSRegion(region2, m_InputImage2);
+    rsr1.Crop(rsr2);
+    return rsr1;
+  }
+private:
+
+  InputImage1PointerType m_InputImage1;
+  InputImage2PointerType m_InputImage2;
+
+};
+
+}
+
+#endif /* GTBRegionComparator_H_ */

include/otbStatisticsImageCustomFilter.h

+/*=========================================================================
+
+  Copyright (c) IRSTEA. All rights reserved.
+
+Some parts of this code are derived from ITK. See ITKCopyright.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 __otbStatisticsImageCustomFilter_h
+#define __otbStatisticsImageCustomFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkNumericTraits.h"
+#include "itkArray.h"
+#include "itkSimpleDataObjectDecorator.h"
+
+namespace otb
+{
+/** \class StatisticsImageCustomFilter
+ * \brief Compute min. max, variance and mean of an Image.
+ *
+ * StatisticsImageCustomFilter computes the minimum, maximum, sum, mean, variance
+ * sigma of an image.  The filter needs all of its input image.  It
+ * behaves as a filter with an input and output. Thus it can be inserted
+ * in a pipline with other filters and the statistics will only be
+ * recomputed if a downstream filter changes.
+ *
+ * The filter passes its input through unmodified.  The filter is
+ * threaded. It computes statistics in each thread then combines them in
+ * its AfterThreadedGenerate method.
+ *
+ * A value to ignore (in the stats) can be chosen
+ *
+ * \ingroup MathematicalStatisticsImageCustomFilters
+ * \ingroup ITKImageStatistics
+ *
+ * \wiki
+ * \wikiexample{Statistics/StatisticsImageCustomFilter,Compute min\, max\, variance and mean of an Image.}
+ * \endwiki
+ */
+template< typename TInputImage >
+class StatisticsImageCustomFilter:
+  public itk::ImageToImageFilter< TInputImage, TInputImage >
+{
+public:
+  /** Standard Self typedef */
+  typedef StatisticsImageCustomFilter                          Self;
+  typedef itk::ImageToImageFilter< TInputImage, TInputImage > Superclass;
+  typedef itk::SmartPointer< Self >                           Pointer;
+  typedef itk::SmartPointer< const Self >                     ConstPointer;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Runtime information support. */
+  itkTypeMacro(StatisticsImageCustomFilter,itk::ImageToImageFilter);
+
+  /** Image related typedefs. */
+  typedef typename TInputImage::Pointer InputImagePointer;
+
+  typedef typename TInputImage::RegionType RegionType;
+  typedef typename TInputImage::SizeType   SizeType;
+  typedef typename TInputImage::IndexType  IndexType;
+  typedef typename TInputImage::PixelType  PixelType;
+
+  /** Image related typedefs. */
+  itkStaticConstMacro(ImageDimension, unsigned int,
+                      TInputImage::ImageDimension);
+
+  /** Type to use for computations. */
+  typedef typename itk::NumericTraits< PixelType >::RealType RealType;
+  typedef unsigned long LongType;
+
+  /** Smart Pointer type to a DataObject. */
+  typedef typename itk::DataObject::Pointer DataObjectPointer;
+
+  /** Type of DataObjects used for scalar outputs */
+  typedef itk::SimpleDataObjectDecorator< RealType >  RealObjectType;
+  typedef itk::SimpleDataObjectDecorator< LongType >  LongObjectType;
+  typedef itk::SimpleDataObjectDecorator< PixelType > PixelObjectType;
+
+  /** Return the computed Minimum. */
+  PixelType GetMinimum() const
+  { return this->GetMinimumOutput()->Get(); }
+  PixelObjectType * GetMinimumOutput();
+
+  const PixelObjectType * GetMinimumOutput() const;
+
+  /** Return the computed Maximum. */
+  PixelType GetMaximum() const
+  { return this->GetMaximumOutput()->Get(); }
+  PixelObjectType * GetMaximumOutput();
+
+  const PixelObjectType * GetMaximumOutput() const;
+
+  /** Return the computed Mean. */
+  RealType GetMean() const
+  { return this->GetMeanOutput()->Get(); }
+  RealObjectType * GetMeanOutput();
+
+  const RealObjectType * GetMeanOutput() const;
+
+  /** Return the computed Standard Deviation. */
+  RealType GetSigma() const
+  { return this->GetSigmaOutput()->Get(); }
+  RealObjectType * GetSigmaOutput();
+
+  const RealObjectType * GetSigmaOutput() const;
+
+  /** Return the computed Variance. */
+  RealType GetVariance() const
+  { return this->GetVarianceOutput()->Get(); }
+  RealObjectType * GetVarianceOutput();
+
+  const RealObjectType * GetVarianceOutput() const;
+
+  /** Return the compute Sum. */
+  RealType GetSum() const
+  { return this->GetSumOutput()->Get(); }
+  RealObjectType * GetSumOutput();
+
+  const RealObjectType * GetSumOutput() const;
+
+  /** Return the compute Sum of squares. */
+  RealType GetSumOfSquares() const
+  { return this->GetSumOfSquaresOutput()->Get(); }
+  RealObjectType * GetSumOfSquaresOutput();
+
+  const RealObjectType * GetSumOfSquaresOutput() const;
+
+  /** Return the compute count. */
+  RealType GetCount() const
+  { return this->GetCountOutput()->Get(); }
+  LongObjectType * GetCountOutput();
+
+  const LongObjectType * GetCountOutput() const;
+
+  /** Make a DataObject of the correct type to be used as the specified
+   * output. */
+  typedef itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType;
+  using Superclass::MakeOutput;
+  virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx);
+
+#ifdef ITK_USE_CONCEPT_CHECKING
+  // Begin concept checking
+  itkConceptMacro( InputHasNumericTraitsCheck,
+                   ( itk::Concept::HasNumericTraits< PixelType > ) );
+  // End concept checking
+#endif
+
+  /** Accessors for ignored value */
+  itkGetMacro(IgnoredInputValue, PixelType);
+  itkSetMacro(IgnoredInputValue, PixelType);
+
+protected:
+  StatisticsImageCustomFilter();
+  ~StatisticsImageCustomFilter(){}
+  void PrintSelf(std::ostream & os, itk::Indent indent) const;
+
+  /** Pass the input through unmodified. Do this by Grafting in the
+   *  AllocateOutputs method.
+   */
+  void AllocateOutputs();
+
+  /** Initialize some accumulators before the threads run. */
+  void BeforeThreadedGenerateData();
+
+  /** Do final mean and variance computation from data accumulated in threads.
+   */
+  void AfterThreadedGenerateData();
+
+  /** Multi-thread version GenerateData. */
+  void  ThreadedGenerateData(const RegionType &
+                             outputRegionForThread,
+                             itk::ThreadIdType threadId);
+
+private:
+  StatisticsImageCustomFilter(const Self &); //purposely not implemented
+  void operator=(const Self &);        //purposely not implemented
+
+  itk::Array< RealType >       m_ThreadSum;
+  itk::Array< RealType >       m_SumOfSquares;
+  itk::Array< itk::SizeValueType >  m_Count;
+  itk::Array< PixelType >      m_ThreadMin;
+  itk::Array< PixelType >      m_ThreadMax;
+
+  PixelType m_IgnoredInputValue;
+}; // end of class
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "otbStatisticsImageCustomFilter.hxx"
+#endif
+
+#endif

include/otbStatisticsImageCustomFilter.hxx

+/*=========================================================================
+
+  Copyright (c) IRSTEA. All rights reserved.
+
+Some parts of this code are derived from ITK. See ITKCopyright.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 __otbStatisticsImageCustomFilter_hxx
+#define __otbStatisticsImageCustomFilter_hxx
+#include "otbStatisticsImageCustomFilter.h"
+
+
+#include "itkImageScanlineIterator.h"
+#include "itkProgressReporter.h"
+
+namespace otb
+{
+template< typename TInputImage >
+StatisticsImageCustomFilter< TInputImage >
+::StatisticsImageCustomFilter():m_ThreadSum(1), m_SumOfSquares(1),
+ m_Count(1), m_ThreadMin(1), m_ThreadMax(1)
+{
+  // first output is a copy of the image, DataObject created by
+  // superclass
+  //
+  // allocate the data objects for the outputs which are
+  // just decorators around pixel types
+  for ( int i = 1; i < 3; ++i )
+    {
+    typename PixelObjectType::Pointer output =
+      static_cast< PixelObjectType * >( this->MakeOutput(i).GetPointer() );
+    this->itk::ProcessObject::SetNthOutput( i, output.GetPointer() );
+    }
+  // allocate the data objects for the outputs which are
+  // just decorators around real types
+  for ( int i = 3; i < 7; ++i )
+    {
+    typename RealObjectType::Pointer output =
+      static_cast< RealObjectType * >( this->MakeOutput(i).GetPointer() );
+    this->itk::ProcessObject::SetNthOutput( i, output.GetPointer() );
+    }
+
+  // allocate the data objects for the count output (long type)
+  typename LongObjectType::Pointer output =
+        static_cast< LongObjectType * >( this->MakeOutput(7).GetPointer() );
+  this->itk::ProcessObject::SetNthOutput( 7, output.GetPointer() );
+
+  this->GetMinimumOutput()->Set( itk::NumericTraits< PixelType >::max() );
+  this->GetMaximumOutput()->Set( itk::NumericTraits< PixelType >::NonpositiveMin() );
+  this->GetMeanOutput()->Set( itk::NumericTraits< RealType >::max() );
+  this->GetSigmaOutput()->Set( itk::NumericTraits< RealType >::max() );
+  this->GetVarianceOutput()->Set( itk::NumericTraits< RealType >::max() );
+  this->GetSumOutput()->Set(itk::NumericTraits< RealType >::Zero);
+  this->GetCountOutput()->Set(itk::NumericTraits< LongType >::Zero);
+
+  m_IgnoredInputValue = itk::NumericTraits<PixelType>::max();
+}
+
+template< typename TInputImage >
+itk::DataObject::Pointer
+StatisticsImageCustomFilter< TInputImage >
+::MakeOutput(DataObjectPointerArraySizeType output)
+{
+  switch ( output )
+    {
+    case 0:
+      return TInputImage::New().GetPointer();
+      break;
+    case 1:
+      return PixelObjectType::New().GetPointer();
+      break;
+    case 2:
+      return PixelObjectType::New().GetPointer();
+      break;
+    case 3:
+    case 4:
+    case 5:
+    case 6:
+      return RealObjectType::New().GetPointer();
+      break;
+    case 7:
+      return LongObjectType::New().GetPointer();
+      break;
+    default:
+      // might as well make an image
+      return TInputImage::New().GetPointer();
+      break;
+    }
+}
+
+template< typename TInputImage >
+typename StatisticsImageCustomFilter< TInputImage >::PixelObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetMinimumOutput()
+{
+  return static_cast< PixelObjectType * >( this->itk::ProcessObject::GetOutput(1) );
+}
+
+template< typename TInputImage >
+const typename StatisticsImageCustomFilter< TInputImage >::PixelObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetMinimumOutput() const
+{
+  return static_cast< const PixelObjectType * >( this->itk::ProcessObject::GetOutput(1) );
+}
+
+template< typename TInputImage >
+typename StatisticsImageCustomFilter< TInputImage >::PixelObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetMaximumOutput()
+{
+  return static_cast< PixelObjectType * >( this->itk::ProcessObject::GetOutput(2) );
+}
+
+template< typename TInputImage >
+const typename StatisticsImageCustomFilter< TInputImage >::PixelObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetMaximumOutput() const
+{
+  return static_cast< const PixelObjectType * >( this->itk::ProcessObject::GetOutput(2) );
+}
+
+template< typename TInputImage >
+typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetMeanOutput()
+{
+  return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(3) );
+}
+
+template< typename TInputImage >
+const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetMeanOutput() const
+{
+  return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(3) );
+}
+
+template< typename TInputImage >
+typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetSigmaOutput()
+{
+  return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(4) );
+}
+
+template< typename TInputImage >
+const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetSigmaOutput() const
+{
+  return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(4) );
+}
+
+template< typename TInputImage >
+typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetVarianceOutput()
+{
+  return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(5) );
+}
+
+template< typename TInputImage >
+const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetVarianceOutput() const
+{
+  return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(5) );
+}
+
+template< typename TInputImage >
+typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetSumOutput()
+{
+  return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(6) );
+}
+
+template< typename TInputImage >
+const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetSumOutput() const
+{
+  return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(6) );
+}
+
+template< typename TInputImage >
+typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetSumOfSquaresOutput()
+{
+  return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(6) );
+}
+
+template< typename TInputImage >
+const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetSumOfSquaresOutput() const
+{
+  return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(6) );
+}
+
+template< typename TInputImage >
+typename StatisticsImageCustomFilter< TInputImage >::LongObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetCountOutput()
+{
+  return static_cast< LongObjectType * >( this->itk::ProcessObject::GetOutput(7) );
+}
+
+template< typename TInputImage >
+const typename StatisticsImageCustomFilter< TInputImage >::LongObjectType *
+StatisticsImageCustomFilter< TInputImage >
+::GetCountOutput() const
+{
+  return static_cast< const LongObjectType * >( this->itk::ProcessObject::GetOutput(7) );
+}
+
+template< typename TInputImage >
+void
+StatisticsImageCustomFilter< TInputImage >
+::AllocateOutputs()
+{
+  // Pass the input through as the output
+  InputImagePointer image =
+    const_cast< TInputImage * >( this->GetInput() );
+
+  this->GraftOutput(image);
+
+  // Nothing that needs to be allocated for the remaining outputs
+}
+
+template< typename TInputImage >
+void
+StatisticsImageCustomFilter< TInputImage >
+::BeforeThreadedGenerateData()
+{
+  itk::ThreadIdType numberOfThreads = this->GetNumberOfThreads();
+
+  // Resize the thread temporaries
+  m_Count.SetSize(numberOfThreads);
+  m_SumOfSquares.SetSize(numberOfThreads);
+  m_ThreadSum.SetSize(numberOfThreads);
+  m_ThreadMin.SetSize(numberOfThreads);
+  m_ThreadMax.SetSize(numberOfThreads);
+
+  // Initialize the temporaries
+  m_Count.Fill(itk::NumericTraits< itk::SizeValueType >::Zero);
+  m_ThreadSum.Fill(itk::NumericTraits< RealType >::Zero);
+  m_SumOfSquares.Fill(itk::NumericTraits< RealType >::Zero);
+  m_ThreadMin.Fill( itk::NumericTraits< PixelType >::max() );
+  m_ThreadMax.Fill( itk::NumericTraits< PixelType >::NonpositiveMin() );
+}
+
+template< typename TInputImage >
+void
+StatisticsImageCustomFilter< TInputImage >
+::AfterThreadedGenerateData()
+{
+  itk::ThreadIdType    i;
+
+  itk::ThreadIdType numberOfThreads = this->GetNumberOfThreads();
+
+  PixelType minimum = this->GetMinimumOutput()->Get();
+  PixelType maximum = this->GetMaximumOutput()->Get();
+  RealType  sum = this ->GetSumOutput()->Get();
+  RealType  sumOfSquares = this->GetSumOfSquaresOutput()->Get();
+  LongType  count = this->GetCountOutput()->Get();
+
+  RealType  mean, sigma, variance;
+  mean = sigma = variance = itk::NumericTraits< RealType >::ZeroValue();
+
+  // Find the min/max over all threads and accumulate count, sum and
+  // sum of squares
+  for ( i = 0; i < numberOfThreads; i++ )
+    {
+    count += m_Count[i];
+    sum += m_ThreadSum[i];
+    sumOfSquares += m_SumOfSquares[i];
+
+    if ( m_ThreadMin[i] < minimum )
+      {
+      minimum = m_ThreadMin[i];
+      }
+    if ( m_ThreadMax[i] > maximum )
+      {
+      maximum = m_ThreadMax[i];
+      }
+    }
+
+  // compute statistics
+  if (count > 1)
+    {
+    mean = sum / static_cast< RealType >( count );
+
+    // unbiased estimate
+    variance = ( sumOfSquares - ( sum * sum / static_cast< RealType >( count ) ) )
+                 / ( static_cast< RealType >( count ) - 1 );
+    sigma = std::sqrt(variance);
+    }
+
+  // Update the outputs
+  this->GetMinimumOutput()->Set(minimum);
+  this->GetMaximumOutput()->Set(maximum);
+  this->GetMeanOutput()->Set(mean);
+  this->GetSigmaOutput()->Set(sigma);
+  this->GetVarianceOutput()->Set(variance);
+  this->GetSumOutput()->Set(sum);
+  this->GetSumOfSquaresOutput()->Set(sumOfSquares);
+  this->GetCountOutput()->Set(count);
+}
+
+template< typename TInputImage >
+void
+StatisticsImageCustomFilter< TInputImage >
+::ThreadedGenerateData(const RegionType & outputRegionForThread,
+    itk::ThreadIdType threadId)
+{
+  const itk::SizeValueType size0 = outputRegionForThread.GetSize(0);
+  if( size0 == 0)
+    {
+    return;
+    }
+  RealType  realValue;
+  PixelType value;
+
+  RealType sum = itk::NumericTraits< RealType >::ZeroValue();
+  RealType sumOfSquares = itk::NumericTraits< RealType >::ZeroValue();
+  itk::SizeValueType count = itk::NumericTraits< itk::SizeValueType >::ZeroValue();
+  PixelType min = itk::NumericTraits< PixelType >::max();
+  PixelType max = itk::NumericTraits< PixelType >::NonpositiveMin();
+
+  itk::ImageScanlineConstIterator< TInputImage > it (this->GetInput(),  outputRegionForThread);
+
+  // support progress methods/callbacks
+  const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0;
+  itk::ProgressReporter progress( this, threadId, numberOfLinesToProcess );
+
+  // do the work
+  while ( !it.IsAtEnd() )
+    {
+    while ( !it.IsAtEndOfLine() )
+      {
+      value = it.Get();
+      if (value != m_IgnoredInputValue)
+        {
+        realValue = static_cast< RealType >( value );
+        if ( value < min )
+          {
+          min = value;
+          }
+        if ( value > max )
+          {
+          max  = value;
+          }
+        sum += realValue;
+        sumOfSquares += ( realValue * realValue );
+        ++count;
+        }
+      ++it;
+      }
+    it.NextLine();
+    progress.CompletedPixel();
+    }
+
+  m_ThreadSum[threadId] = sum;
+  m_SumOfSquares[threadId] = sumOfSquares;
+  m_Count[threadId] = count;
+  m_ThreadMin[threadId] = min;
+  m_ThreadMax[threadId] = max;
+}
+
+template< typename TImage >
+void
+StatisticsImageCustomFilter< TImage >
+::PrintSelf(std::ostream & os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+
+  os << indent << "Minimum: "
+     << static_cast< typename itk::NumericTraits< PixelType >::PrintType >( this->GetMinimum() ) << std::endl;
+  os << indent << "Maximum: "
+     << static_cast< typename itk::NumericTraits< PixelType >::PrintType >( this->GetMaximum() ) << std::endl;
+  os << indent << "Sum: "      << this->GetSum() << std::endl;
+  os << indent << "SumOfSquares: "      << this->GetSumOfSquares() << std::endl;
+  os << indent << "Count: "      << this->GetCount() << std::endl;
+  os << indent << "Mean: "     << this->GetMean() << std::endl;
+  os << indent << "Sigma: "    << this->GetSigma() << std::endl;
+  os << indent << "Variance: " << this->GetVariance() << std::endl;
+}
+} // end namespace itk
+#endif

include/otbVectorDataToLabelImageCustomFilter.h

+/*=========================================================================
+
+  Copyright (c) IRSTEA. All rights reserved.
+
+
+     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 __gtbVectorDataToLabelImageCustomFilter_h
+#define __gtbVectorDataToLabelImageCustomFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkImageSource.h"
+#include "otbMacro.h"
+#include "otbImageMetadataInterfaceFactory.h"
+#include "otbVectorData.h"
+
+#include "gdal.h"
+#include "ogr_api.h"
+#include <ogrsf_frmts.h>
+
+namespace otb {
+
+/** \class VectorDataToLabelImageCustomFilter
+ *  \brief Burn geometries from the specified VectorData into raster
+ *
+ *  This class handles burning several input VectorDatas into the
+ *  output raster.  The burn values are extracted from a field set by
+ *  the user.If no burning field is set, the "FID" is used by default.
+ *
+ *  A "Burn max value" mode overrides this and makes only the max value
+ *  burnt whenever the burning field.
+ *
+ *  Setting the output raster informations can be done in two ways by:
+ *    - Setting the Origin/Size/Spacing of the output image
+ *    - Using an existing image as support via SetOutputParametersFromImage(ImageBase)
+ *
+ *  OGRRegisterAll() method must have been called before applying filter.
+ *
+ */
+template <class TVectorData, class TOutputImage  >
+class  ITK_EXPORT VectorDataToLabelImageCustomFilter :
+    public itk::ImageSource<TOutputImage>
+{
+public:
+  /** Standard class typedefs */
+  typedef VectorDataToLabelImageCustomFilter                           Self;
+  typedef itk::ImageSource<TOutputImage>     Superclass;
+  typedef itk::SmartPointer< Self >                           Pointer;
+  typedef itk::SmartPointer<const Self>                       ConstPointer;
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(VectorDataToLabelImageCustomFilter, itk::ImageSource);
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  typedef TOutputImage                                OutputImageType;
+  typedef typename OutputImageType::Pointer           OutputImagePointer;
+  typedef typename OutputImageType::SizeType          OutputSizeType;
+  typedef typename OutputImageType::IndexType         OutputIndexType;
+  typedef typename OutputImageType::SpacingType       OutputSpacingType;
+  typedef typename OutputImageType::PointType         OutputOriginType;
+  typedef typename OutputImageType::RegionType        OutputImageRegionType;
+  typedef typename OutputImageType::PixelType         OutputImagePixelType;
+  typedef typename OutputImageType::InternalPixelType OutputImageInternalPixelType;
+
+  /** VectorData typedefs*/
+  typedef TVectorData                            VectorDataType;
+  typedef typename VectorDataType::DataTreeType  DataTreeType;
+  typedef typename DataTreeType::TreeNodeType    InternalTreeNodeType;
+  typedef typename DataTreeType::Pointer         DataTreePointerType;
+  typedef typename DataTreeType::ConstPointer    DataTreeConstPointerType;
+
+  typedef itk::ImageBase<OutputImageType::ImageDimension>      ImageBaseType;
+
+  /** Get Nth input VectorData */
+  const VectorDataType*  GetInput(unsigned int idx);
+
+  /** Method for adding VectorData to rasterize  */
+  virtual void AddVectorData(const VectorDataType* vd);
+
+  /** Set the size of the output image. */
+  itkSetMacro(OutputSize, OutputSizeType);
+
+  /** Get the size of the output image. */
+  itkGetConstReferenceMacro(OutputSize, OutputSizeType);
+
+  /** Burn mode */
+  itkSetMacro(BurnMaxValueMode, bool);
+  itkGetMacro(BurnMaxValueMode, bool);
+
+
+  /** Set the origin of the output image.
+   * \sa GetOrigin()
+   */
+  itkSetMacro(OutputOrigin, OutputOriginType);
+  virtual void SetOutputOrigin(const double origin[2]);
+  virtual void SetOutputOrigin(const float origin[2]);
+
+  itkGetConstReferenceMacro(OutputOrigin, OutputOriginType);
+
+  /** Set the spacing (size of a pixel) of the output image.
+  * \sa GetSpacing()
+  */
+  virtual void SetOutputSpacing(const OutputSpacingType& spacing);
+  virtual void SetOutputSpacing(const double spacing[2]);
+  virtual void SetOutputSpacing(const float spacing[2]);
+
+  /** Set/Get Output Projection Ref */
+  itkSetStringMacro(OutputProjectionRef);
+  itkGetStringMacro(OutputProjectionRef);
+
+  itkSetStringMacro(BurnAttribute);
+  itkGetStringMacro(BurnAttribute);
+
+  /** Useful to set the output parameters from an existing image*/
+  void SetOutputParametersFromImage(const ImageBaseType * image);
+
+protected:
+  virtual void GenerateData();
+
+  VectorDataToLabelImageCustomFilter();
+  virtual ~VectorDataToLabelImageCustomFilter()
+  {
+    // Destroy the geometries stored
+    for (unsigned int idx = 0; idx < m_SrcDataSetGeometries.size(); ++idx)
+      {
+      OGR_G_DestroyGeometry(m_SrcDataSetGeometries[idx]);
+      }
+
+    if (m_OGRDataSourcePointer != NULL)
+      {
+      OGRDataSource::DestroyDataSource(m_OGRDataSourcePointer);
+      }
+  }
+
+  virtual void GenerateOutputInformation();
+
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+private:
+  VectorDataToLabelImageCustomFilter(const Self&); //purposely not implemented
+  void operator=(const Self&); //purposely not implemented
+
+  OGRDataSource*                m_OGRDataSourcePointer;
+
+  // Vector Of OGRGeometyH
+  std::vector< OGRGeometryH >   m_SrcDataSetGeometries;
+
+  std::vector<double>           m_BurnValues;
+  std::vector<double>           m_FullBurnValues;
+  std::vector<int>              m_BandsToBurn;
+
+  // Field used to extract the burn value
+  std::string                   m_BurnAttribute;
+
+  // Burn mode
+  bool							m_BurnMaxValueMode;
+
+  // Default burn value
+  double                        m_DefaultBurnValue;
+
+  // Output params
+  std::string                   m_OutputProjectionRef;
+  OutputSpacingType             m_OutputSpacing;
+  OutputOriginType              m_OutputOrigin;
+  OutputSizeType                m_OutputSize;
+  OutputIndexType               m_OutputStartIndex;
+}; // end of class VectorDataToLabelImageCustomFilter
+
+} // end of namespace otb
+
+
+#ifndef  OTB_MANUAL_INSTANTIATION
+#include "otbVectorDataToLabelImageCustomFilter.hxx"
+#endif
+
+#endif
+

include/otbVectorDataToLabelImageCustomFilter.hxx

+/*=========================================================================
+
+  Copyright (c) IRSTEA. All rights reserved.
+
+
+     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 __gtbVectorDataToLabelImageCustomFilter_txx
+#define __gtbVectorDataToLabelImageCustomFilter_txx
+
+#include "otbVectorDataToLabelImageCustomFilter.h"
+#include "otbOGRIOHelper.h"
+#include "otbGdalDataTypeBridge.h"
+#include "otbMacro.h"
+
+#include "gdal_alg.h"
+#include "ogr_srs_api.h"
+
+#include "itkImageRegionIterator.h"
+
+namespace otb
+{
+template<class TVectorData, class TOutputImage>
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::VectorDataToLabelImageCustomFilter()
+ : m_OGRDataSourcePointer(0),
+   m_BurnAttribute("FID"),
+   m_BurnMaxValueMode(false)
+   {
+  this->SetNumberOfRequiredInputs(1);
+
+  // Output parameters initialization
+  m_OutputSpacing.Fill(1.0);
+  m_OutputSize.Fill(0);
+  m_OutputStartIndex.Fill(0);
+
+  // This filter is intended to work with otb::Image
+  m_BandsToBurn.push_back(1);
+
+  // Default burn value if no burnAttribute available
+  m_DefaultBurnValue = 1.;
+   }
+
+template<class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::AddVectorData(const VectorDataType* vd)
+ {
+  this->itk::ProcessObject::PushBackInput( vd );
+ }
+
+template <class TVectorData, class TOutputImage>
+const typename VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>::VectorDataType *
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::GetInput(unsigned int idx)
+ {
+  return static_cast<const TVectorData *>
+  (this->itk::ProcessObject::GetInput(idx));
+ }
+
+template <class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::SetOutputSpacing(const OutputSpacingType& spacing)
+ {
+  if (this->m_OutputSpacing != spacing)
+    {
+    this->m_OutputSpacing = spacing;
+    this->Modified();
+    }
+ }
+
+template <class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::SetOutputSpacing(const double spacing[2])
+ {
+  OutputSpacingType s(spacing);
+  this->SetOutputSpacing(s);
+ }
+
+template <class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::SetOutputSpacing(const float spacing[2])
+ {
+  itk::Vector<float, 2> sf(spacing);
+  OutputSpacingType s;
+  s.CastFrom(sf);
+  this->SetOutputSpacing(s);
+ }
+
+template <class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::SetOutputOrigin(const double origin[2])
+ {
+  OutputOriginType p(origin);
+  this->SetOutputOrigin(p);
+ }
+
+template <class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::SetOutputOrigin(const float origin[2])
+ {
+  itk::Point<float, 2> of(origin);
+  OutputOriginType p;
+  p.CastFrom(of);
+  this->SetOutputOrigin(p);
+ }
+
+template <class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::SetOutputParametersFromImage(const ImageBaseType * src)
+ {
+  this->SetOutputOrigin ( src->GetOrigin() );
+  this->SetOutputSpacing ( src->GetSpacing() );
+  this->SetOutputSize ( src->GetLargestPossibleRegion().GetSize() );
+  otb::ImageMetadataInterfaceBase::Pointer imi = otb::ImageMetadataInterfaceFactory::CreateIMI(src->GetMetaDataDictionary());
+  this->SetOutputProjectionRef(imi->GetProjectionRef());
+ }
+
+template<class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::GenerateOutputInformation()
+ {
+  // get pointer to the output
+  OutputImagePointer outputPtr = this->GetOutput();
+  if (!outputPtr)
+    {
+    return;
+    }
+
+  // Set the size of the output region
+  typename TOutputImage::RegionType outputLargestPossibleRegion;
+  outputLargestPossibleRegion.SetSize(m_OutputSize);
+  //outputLargestPossibleRegion.SetIndex(m_OutputStartIndex);
+  outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
+
+  // Set spacing and origin
+  outputPtr->SetSpacing(m_OutputSpacing);
+  outputPtr->SetOrigin(m_OutputOrigin);
+
+  itk::MetaDataDictionary& dict = outputPtr->GetMetaDataDictionary();
+  itk::EncapsulateMetaData<std::string> (dict, otb::MetaDataKey::ProjectionRefKey,
+      static_cast<std::string>(this->GetOutputProjectionRef()));
+
+  // Generate the OGRLayers from the input VectorDatas
+  // iteration begin from 1 cause the 0th input is a image
+  for (unsigned int inputIdx = 0; inputIdx < this->GetNumberOfInputs(); ++inputIdx)
+    {
+    const VectorDataType* vd = dynamic_cast< const VectorDataType*>(this->itk::ProcessObject::GetInput(inputIdx));
+
+    // Get the projection ref of the current VectorData
+    std::string projectionRefWkt = vd->GetProjectionRef();
+    bool        projectionInformationAvailable = !projectionRefWkt.empty();
+    OGRSpatialReference * oSRS = NULL;
+
+    if (projectionInformationAvailable)
+      {
+      oSRS = static_cast<OGRSpatialReference *>(OSRNewSpatialReference(projectionRefWkt.c_str()));
+      }
+    else
+      {
+      otbMsgDevMacro(<< "Projection information unavailable");
+      }
+
+    // Retrieving root node
+    DataTreeConstPointerType tree = vd->GetDataTree();
+
+    // Get the input tree root
+    InternalTreeNodeType * inputRoot = const_cast<InternalTreeNodeType *>(tree->GetRoot());
+
+    // Iterative method to build the layers from a VectorData
+    OGRLayer *   ogrCurrentLayer = NULL;
+    std::vector<OGRLayer *> ogrLayerVector;
+    otb::OGRIOHelper::Pointer IOConversion = otb::OGRIOHelper::New();
+
+    // The method ConvertDataTreeNodeToOGRLayers create the
+    // OGRDataSource but don t release it. Destruction is done in the
+    // desctructor
+    m_OGRDataSourcePointer = NULL;
+    ogrLayerVector = IOConversion->ConvertDataTreeNodeToOGRLayers(inputRoot,
+        m_OGRDataSourcePointer,
+        ogrCurrentLayer,
+        oSRS);
+
+    // From OGRLayer* to OGRGeometryH vector
+    for (unsigned int idx = 0; idx < ogrLayerVector.size(); ++idx)
+      {
+      // test if the layers contain a field m_BurnField;
+      int burnField = -1;
+
+      if( !m_BurnAttribute.empty() )
+        {
+        burnField = OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( (OGRLayerH)(ogrLayerVector[idx]) ),
+            m_BurnAttribute.c_str() );
+
+        // Get the geometries of the layer
+        OGRFeatureH hFeat;
+        OGR_L_ResetReading( (OGRLayerH)(ogrLayerVector[idx]) );
+        while( ( hFeat = OGR_L_GetNextFeature( (OGRLayerH)(ogrLayerVector[idx]) )) != NULL )
+          {
+          OGRGeometryH hGeom;
+          if( OGR_F_GetGeometryRef( hFeat ) == NULL )
+            {
+            OGR_F_Destroy( hFeat );
+            continue;
+            }
+
+          hGeom = OGR_G_Clone( OGR_F_GetGeometryRef( hFeat ) );
+          m_SrcDataSetGeometries.push_back( hGeom );
+
+          if (m_BurnMaxValueMode)
+            {
+            m_FullBurnValues.push_back(static_cast<double>(itk::NumericTraits<OutputImageInternalPixelType>::max()));
+            }
+          else if (burnField == -1 )
+            {
+            // TODO : if no burnAttribute available, warning or raise an exception??
+            m_FullBurnValues.push_back(m_DefaultBurnValue++);
+            itkWarningMacro(<<"Failed to find attribute "<<m_BurnAttribute << " in layer "
+                << OGR_FD_GetName( OGR_L_GetLayerDefn( (OGRLayerH)(ogrLayerVector[idx]) ))
+                <<" .Setting burn value to default =  "
+                << m_DefaultBurnValue);
+            }
+          else
+            {
+            m_FullBurnValues.push_back( OGR_F_GetFieldAsDouble( hFeat, burnField ) );
+            }
+
+          OGR_F_Destroy( hFeat );
+          }
+        }
+
+      // Destroy the oSRS
+      if (oSRS != NULL)
+        {
+        OSRRelease(oSRS);
+        }
+      }
+    }
+ }
+
+template<class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>::GenerateData()
+{
+  // Call Superclass GenerateData
+  this->AllocateOutputs();
+  itk::ImageRegionIterator<TOutputImage> it (this->GetOutput(), this->GetOutput()->GetBufferedRegion());
+  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
+    {
+    it.Set(0);
+    }
+
+  // Get the buffered region
+  OutputImageRegionType bufferedRegion = this->GetOutput()->GetBufferedRegion();
+
+  // nb bands
+  unsigned int nbBands =  this->GetOutput()->GetNumberOfComponentsPerPixel();
+
+  // register drivers
+  GDALAllRegister();
+
+  std::ostringstream stream;
+  stream << "MEM:::"
+      <<  "DATAPOINTER=" << (unsigned long)(this->GetOutput()->GetBufferPointer()) << ","
+      <<  "PIXELS=" << bufferedRegion.GetSize()[0] << ","
+      <<  "LINES=" << bufferedRegion.GetSize()[1]<< ","
+      <<  "BANDS=" << nbBands << ","
+      <<  "DATATYPE=" << GDALGetDataTypeName(otb::GdalDataTypeBridge::GetGDALDataType<OutputImageInternalPixelType>()) << ","
+      <<  "PIXELOFFSET=" << sizeof(OutputImageInternalPixelType) *  nbBands << ","
+      <<  "LINEOFFSET=" << sizeof(OutputImageInternalPixelType)*nbBands*bufferedRegion.GetSize()[0] << ","
+      <<  "BANDOFFSET=" << sizeof(OutputImageInternalPixelType);
+
+  GDALDatasetH dataset = GDALOpen(stream.str().c_str(), GA_Update);
+
+  // Add the projection ref to the dataset
+  GDALSetProjection (dataset, this->GetOutput()->GetProjectionRef().c_str());
+
+  // add the geoTransform to the dataset
+  itk::VariableLengthVector<double> geoTransform(6);
+
+  // Reporting origin and spacing of the buffered region
+  // the spacing is unchanged, the origin is relative to the buffered region
+  OutputIndexType  bufferIndexOrigin = bufferedRegion.GetIndex();
+  OutputOriginType bufferOrigin;
+  this->GetOutput()->TransformIndexToPhysicalPoint(bufferIndexOrigin, bufferOrigin);
+  geoTransform[0] = bufferOrigin[0];
+  geoTransform[3] = bufferOrigin[1];
+  geoTransform[1] = this->GetOutput()->GetSpacing()[0];
+  geoTransform[5] = this->GetOutput()->GetSpacing()[1];
+
+  // FIXME: Here component 1 and 4 should be replaced by the orientation parameters
+  geoTransform[2] = 0.;
+  geoTransform[4] = 0.;
+  GDALSetGeoTransform(dataset,const_cast<double*>(geoTransform.GetDataPointer()));
+
+  // Burn the geometries into the dataset
+  if (dataset != NULL)
+    {
+    GDALRasterizeGeometries( dataset, m_BandsToBurn.size(),
+        &(m_BandsToBurn[0]),
+        m_SrcDataSetGeometries.size(),
+        &(m_SrcDataSetGeometries[0]),
+        NULL, NULL, &(m_FullBurnValues[0]),
+        NULL,
+        GDALDummyProgress, NULL );
+
+    // release the dataset
+    GDALClose( dataset );
+    }
+}
+
+template<class TVectorData, class TOutputImage>
+void
+VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+ {
+  Superclass::PrintSelf(os, indent);
+ }
+
+} // end namespace otb
+
+#endif
+
+

otb-module.cmake

+set(DOCUMENTATION "Extraction of subset of remote sensing images")
+
+otb_module(SimpleExtractionTools
+  DEPENDS
+    OTBCommon
+	OTBApplicationEngine
+	OTBConversion
+	OTBImageBase
+  TEST_DEPENDS
+    OTBTestKernel
+    OTBCommandLine
+  DESCRIPTION
+    $DOCUMENTATION
+)