diff --git a/CMake/OTBModuleMacros.cmake b/CMake/OTBModuleMacros.cmake
index 0938a0f2fb9ed50fe44edf7017d9a4ac325e41b7..afb9a29f5c45b2509d0cd0bfe350df5a59ea846f 100644
--- a/CMake/OTBModuleMacros.cmake
+++ b/CMake/OTBModuleMacros.cmake
@@ -277,6 +277,12 @@ macro(otb_module_test)
   foreach(dep IN LISTS OTB_MODULE_${otb-module-test}_DEPENDS)
     list(APPEND ${otb-module-test}_LIBRARIES "${${dep}_LIBRARIES}")
   endforeach()
+  # make sure the test can link with optional libs
+  foreach(dep IN LISTS OTB_MODULE_${otb-module}_OPTIONAL_DEPENDS)
+    if (${dep}_ENABLED)
+      list(APPEND ${otb-module-test}_LIBRARIES "${${dep}_LIBRARIES}")
+    endif()
+  endforeach()
 endmacro()
 
 macro(otb_module_warnings_disable)
diff --git a/Documentation/Cookbook/rst/recipes/python.rst b/Documentation/Cookbook/rst/recipes/python.rst
index 41306ac1ead88fee83e005cc9000ba3c25498b4b..f4d5ff597e3a8cc81b46dc70fa1e760c4eeeccec 100644
--- a/Documentation/Cookbook/rst/recipes/python.rst
+++ b/Documentation/Cookbook/rst/recipes/python.rst
@@ -29,14 +29,14 @@ application, changing the algorithm at each iteration.
     import otbApplication
 
     # otbApplication.Registry can tell you what application are available
-    print "Available applications: "
-    print str( otbApplication.Registry.GetAvailableApplications() )
+    print('Available applications: ')
+    print (str( otbApplication.Registry.GetAvailableApplications()))
 
     # Let's create the application  "Smoothing"
     app = otbApplication.Registry.CreateApplication("Smoothing")
 
     # We print the keys of all its parameters
-    print app.GetParametersKeys()
+    print (app.GetParametersKeys())
 
     # First, we set the input image filename
     app.SetParameterString("in", argv[1])
@@ -45,7 +45,7 @@ application, changing the algorithm at each iteration.
     # and can take 3 values: 'mean', 'gaussian', 'anidif'
     for type in ['mean', 'gaussian', 'anidif']:
 
-      print 'Running with ' + type + ' smoothing type'
+      print('Running with ' + type + ' smoothing type')
 
       # Now we configure the smoothing algorithm
       app.SetParameterString("type", type)
@@ -303,7 +303,7 @@ Here is a small example of what can be done:
   # Check the result of ReadImageInfo
   someKeys = ['sizex', 'sizey', 'spacingx', 'spacingy', 'sensor', 'projectionref']
   for key in someKeys:
-    print(key + ' : ' + str(app2.GetParameterValue(key)) )
+    print(key + ' : ' + str(app2.GetParameterValue(key)))
   
   # Only a portion of "out" was exported but ReadImageInfo is still able to detect the 
   # correct full size of the image
diff --git a/Documentation/SoftwareGuide/Latex/CompilingMonteverdi.tex b/Documentation/SoftwareGuide/Latex/CompilingMonteverdi.tex
deleted file mode 100644
index 7c38627bd2838e95a5936225d7b4dbb664cc2b08..0000000000000000000000000000000000000000
--- a/Documentation/SoftwareGuide/Latex/CompilingMonteverdi.tex
+++ /dev/null
@@ -1,59 +0,0 @@
-\setcounter{secnumdepth}{3}
-
-\chapter{Compiling Monteverdi from source}
-\label{chapter:CompilingMonteverdi}
-
-\section{Linux and Mac OS X}
-
-Compiling Monteverdi from source follows the same procedure as a regular CMake project: setup a directory for an
-out-of-source build, configure and compile.
-
-\subsection{Monteverdi}
-
-Make sure OTB is compiled with at least \texttt{OTB\_USE\_QT4}, \texttt{OTB\_USE\_GLUT}, \texttt{OTB\_USE\_GLFW} and \texttt{OTB\_USE\_GLEW} set to ON.
-Setup another out-of-source build environment for Monteverdi:
-\begin{verbatim}
-$ mkdir ~/monteverdi
-$ cd ~/monteverdi
-$ git clone https://git@git.orfeo-toolbox.org/git/monteverdi2.git
-$ mkdir build
-$ mkdir install
-\end{verbatim}
-
-Remember to checkout the develop branch if you want the current development version:
-\begin{verbatim}
-$ cd ~/monteverdi/monteverdi2
-$ git checkout develop
-\end{verbatim}
-
-CMake needs to find both OTB and QWT installation locations.
-For example, set an CMake cache pre-population script with the following content:
-\begin{verbatim}
-# monteverdi-configuration.cmake
-set(CMAKE_INSTALL_PREFIX "~/monteverdi/install" CACHE STRING "" FORCE)
-set(OTB_DIR "~/OTB/install/lib/cmake/OTB-5.0" CACHE STRING "" FORCE)
-set(QWT_INCLUDE_DIR "/usr/include/qwt5-qt4" CACHE STRING "" FORCE)
-set(QWT_LIBRARY "/usr/lib64/libqwt.so.5" CACHE STRING "" FORCE)
-\end{verbatim}
-
-Configure and compile monteverdi:
-\begin{verbatim}
-$ cd ~/monteverdi/build
-$ cmake -C monteverdi-configuration.cmake ../monteverdi2
-$ make
-$ make install
-\end{verbatim}
-
-Remember to set the \texttt{OTB\_APPLICATION\_PATH} environment variable to
-allow monteverdi to find the OTB applications:
-\begin{verbatim}
-export OTB_APPLICATION_PATH=~/OTB/build/lib/otb/applications
-./bin/monteverdi
-\end{verbatim}
-
-\section{Windows}
-
-Everything that is needed for Monteverdi development on Windows, including compiling from source, is covered in details on the OTB wiki at:
-\begin{center}
-\url{http://wiki.orfeo-toolbox.org/index.php/OTB_development_on_Windows}
-\end{center}
diff --git a/Documentation/SoftwareGuide/Latex/SoftwareGuide.tex b/Documentation/SoftwareGuide/Latex/SoftwareGuide.tex
index c5bd1dff8218ca02c31995c44a83913874e97c27..8cf39aa7a0f2cfff48375656602a0a2404cf053c 100644
--- a/Documentation/SoftwareGuide/Latex/SoftwareGuide.tex
+++ b/Documentation/SoftwareGuide/Latex/SoftwareGuide.tex
@@ -217,12 +217,8 @@ colorlinks,linkcolor={blue},citecolor={blue},urlcolor={blue},
 
 \input{Introduction.tex}
 \input{Installation.tex}
-%Monteverdi is integrated in OTB as a module since release 5.8 so there is no
-%need to maintain a specific section regarding monteverdi compilation.
-%\input{CompilingMonteverdi.tex}
 \input{SystemOverview.tex}
 
-
 \part{Tutorials}\label{part:tutorials}
 
 \input{Tutorial.tex}
diff --git a/Modules/Applications/AppClassification/app/otbKMeansClassification.cxx b/Modules/Applications/AppClassification/app/otbKMeansClassification.cxx
index 52e0d6b426d7defb6343113029c17abe6dcacd15..ee04d9cd0fb2ff3035ca2208a8cd64c735cd93d8 100644
--- a/Modules/Applications/AppClassification/app/otbKMeansClassification.cxx
+++ b/Modules/Applications/AppClassification/app/otbKMeansClassification.cxx
@@ -291,12 +291,15 @@ protected:
         itkExceptionMacro(<< "File : " << modelFileName << " couldn't be opened");
       }
 
-      // get the end line with the centroids
+      // get the line with the centroids (starts with "2 ")
       std::string line, centroidLine;
       while(std::getline(infile,line))
       {
-        if (!line.empty())
+        if (line.size() > 2 && line[0] == '2' && line[1] == ' ')
+          {
           centroidLine = line;
+          break;
+          }
       }
 
       std::vector<std::string> centroidElm;
diff --git a/Modules/Applications/AppDimensionalityReduction/app/otbImageDimensionalityReduction.cxx b/Modules/Applications/AppDimensionalityReduction/app/otbImageDimensionalityReduction.cxx
index 1df2c463fc26c531095154515039a9245d0ba5c1..d9dd8e8816fec2cce2502a81b624a59ce33817ab 100644
--- a/Modules/Applications/AppDimensionalityReduction/app/otbImageDimensionalityReduction.cxx
+++ b/Modules/Applications/AppDimensionalityReduction/app/otbImageDimensionalityReduction.cxx
@@ -124,12 +124,12 @@ protected:
 private:
   void DoInit() override
   {
-    SetName("DimensionalityReduction");
+    SetName("ImageDimensionalityReduction");
     SetDescription("Performs dimensionality reduction of the input image "
       "according to a dimensionality reduction model file.");
 
     // Documentation
-    SetDocName("DimensionalityReduction");
+    SetDocName("Image Dimensionality Reduction");
     SetDocLongDescription("This application reduces the dimension of an input"
                           " image, based on a machine learning model file produced by"
                           " the TrainDimensionalityReduction application. Pixels of the "
diff --git a/Modules/Core/Streaming/include/otbStreamingManager.h b/Modules/Core/Streaming/include/otbStreamingManager.h
index cd57e83f9e7ddeda81282d35af4d048a4d837faa..0bffd7e95779c755c41e2a94a276a72360845e1d 100644
--- a/Modules/Core/Streaming/include/otbStreamingManager.h
+++ b/Modules/Core/Streaming/include/otbStreamingManager.h
@@ -67,6 +67,7 @@ public:
   typedef typename ImageType::InternalPixelType PixelType;
 
   typedef otb::PipelineMemoryPrintCalculator::MemoryPrintType MemoryPrintType;
+  typedef itk::ImageRegionSplitterBase          AbstractSplitterType;
 
   /** Type macro */
   itkTypeMacro(StreamingManager, itk::LightObject);
@@ -74,6 +75,8 @@ public:
   /** Dimension of input image. */
   itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension);
 
+  const AbstractSplitterType * GetSplitter() const;
+
   /** Actually computes the stream divisions, according to the specified streaming mode,
    * eventually using the input parameter to estimate memory consumption */
   virtual void PrepareStreaming(itk::DataObject * input, const RegionType &region) = 0;
@@ -106,7 +109,6 @@ protected:
   RegionType m_Region;
 
   /** The splitter used to compute the different strips */
-  typedef itk::ImageRegionSplitterBase           AbstractSplitterType;
   typedef typename AbstractSplitterType::Pointer AbstractSplitterPointerType;
   AbstractSplitterPointerType m_Splitter;
 
diff --git a/Modules/Core/Streaming/include/otbStreamingManager.txx b/Modules/Core/Streaming/include/otbStreamingManager.txx
index 1d76ff1fe498007c5b9bb60fcc7fa80e040733b2..24e487a4e55998ee7d28b1be1a58cbda21ee8768 100644
--- a/Modules/Core/Streaming/include/otbStreamingManager.txx
+++ b/Modules/Core/Streaming/include/otbStreamingManager.txx
@@ -40,6 +40,13 @@ StreamingManager<TImage>::~StreamingManager()
 {
 }
 
+template <class TImage>
+const typename StreamingManager<TImage>::AbstractSplitterType *
+StreamingManager<TImage>::GetSplitter() const
+{
+  return m_Splitter;
+}
+
 template <class TImage>
 typename StreamingManager<TImage>::MemoryPrintType
 StreamingManager<TImage>::GetActualAvailableRAMInBytes(MemoryPrintType availableRAMInMB)
diff --git a/Modules/Feature/Edge/include/otbLineSegmentDetector.txx b/Modules/Feature/Edge/include/otbLineSegmentDetector.txx
index c06344e1831c17e0f9c07448a40f7ceb1b154155..8ff969c22c23d5558295af53589c62a5d8168e38 100644
--- a/Modules/Feature/Edge/include/otbLineSegmentDetector.txx
+++ b/Modules/Feature/Edge/include/otbLineSegmentDetector.txx
@@ -792,8 +792,14 @@ LineSegmentDetector<TInputImage, TPrecision>
     l =  (static_cast<double>((*it)[0]) - x) * dx + (static_cast<double>((*it)[1]) - y) * dy;
     w = -(static_cast<double>((*it)[0]) - x) * dy + (static_cast<double>((*it)[1]) - y) * dx;
 
-    if (l < l_min) l_min = l; if (l > l_max) l_max = l;
-    if (w < w_min) w_min = w; if (w > w_max) w_max = w;
+    if (l < l_min)
+      l_min = l;
+    if (l > l_max)
+      l_max = l;
+    if (w < w_min)
+      w_min = w;
+    if (w > w_max)
+      w_max = w;
 
     sum_l[static_cast < int > (vcl_floor(l) + 0.5) + Diagonal] +=  static_cast<MagnitudePixelType>(weight);
     sum_w[static_cast < int > (vcl_floor(w) + 0.5) + Diagonal] +=  static_cast<MagnitudePixelType>(weight);
diff --git a/Modules/Filtering/Path/include/otbClosePathFunctor.h b/Modules/Filtering/Path/include/otbClosePathFunctor.h
index 90924a2b427f09eedc3289678f17115078c5cd1f..d880196c9f4b9d4992fd538fc8c988f74369fcdc 100644
--- a/Modules/Filtering/Path/include/otbClosePathFunctor.h
+++ b/Modules/Filtering/Path/include/otbClosePathFunctor.h
@@ -59,6 +59,9 @@ public:
 
     if(input->GetVertexList()->Size()>0)
       {
+      //Initialization of lastVertex to GetVertexList
+      lastVertex = input->GetVertexList()->Begin().Value();
+      
       for (VertexListConstIteratorType vertexIt = input->GetVertexList()->Begin();
           vertexIt != input->GetVertexList()->End();
           ++vertexIt)
diff --git a/Modules/IO/ImageIO/include/otbImageFileWriter.h b/Modules/IO/ImageIO/include/otbImageFileWriter.h
index cbc4c6c4e9a7fb49c3cdfe8821fec27ca84ea189..169e946c7579a241c7bed42424c561101d608912 100644
--- a/Modules/IO/ImageIO/include/otbImageFileWriter.h
+++ b/Modules/IO/ImageIO/include/otbImageFileWriter.h
@@ -213,6 +213,9 @@ protected:
   /** Does the real work. */
   void GenerateData(void) override;
 
+  /** Prepare the streaming and write the output information on disk */
+  void GenerateOutputInformation(void) override;
+
 private:
   ImageFileWriter(const ImageFileWriter &); //purposely not implemented
   void operator =(const ImageFileWriter&); //purposely not implemented
diff --git a/Modules/IO/ImageIO/include/otbImageFileWriter.txx b/Modules/IO/ImageIO/include/otbImageFileWriter.txx
index 2ddaf018bf89a48109c86a902d2c6fa57c88089b..6cfff64f8ae1df38f26344110537bb1c47d0336c 100644
--- a/Modules/IO/ImageIO/include/otbImageFileWriter.txx
+++ b/Modules/IO/ImageIO/include/otbImageFileWriter.txx
@@ -263,13 +263,11 @@ ImageFileWriter<TInputImage>
   return static_cast<const InputImageType*>(this->ProcessObject::GetInput(0));
 }
 
-/**
- * Update method : update output information of input and write to file
- */
+/** Prepare everything and call m_ImageIO.WriteInformation() */
 template<class TInputImage>
 void
 ImageFileWriter<TInputImage>
-::Update()
+::GenerateOutputInformation(void)
 {
   // Update output information on input image
   InputImagePointer inputPtr =
@@ -398,14 +396,6 @@ ImageFileWriter<TInputImage>
       }
     }
 
-  this->SetAbortGenerateData(0);
-  this->SetProgress(0.0);
-
-  /**
-   * Tell all Observers that the filter is starting
-   */
-  this->InvokeEvent(itk::StartEvent());
-
   /** Prepare ImageIO  : create ImageFactory */
 
   if (m_FileName == "")
@@ -477,7 +467,6 @@ ImageFileWriter<TInputImage>
   /**
    * Grab the input
    */
-  inputPtr->UpdateOutputInformation();
   InputImageRegionType inputRegion = inputPtr->GetLargestPossibleRegion();
 
   /** Parse region size modes */
@@ -541,12 +530,6 @@ ImageFileWriter<TInputImage>
   const auto firstSplitSize = m_StreamingManager->GetSplit(0).GetSize();
   otbLogMacro(Info,<<"File "<<m_FileName<<" will be written in "<<m_NumberOfDivisions<<" blocks of "<<firstSplitSize[0]<<"x"<<firstSplitSize[1]<<" pixels");
 
-  /**
-   * Loop over the number of pieces, execute the upstream pipeline on each
-   * piece, and copy the results into the output image.
-   */
-  InputImageRegionType streamRegion;
-
   //
   // Setup the ImageIO with information from inputPtr
   //
@@ -586,12 +569,33 @@ ImageFileWriter<TInputImage>
   m_ImageIO->SetFileName(m_FileName.c_str());
 
   m_ImageIO->WriteImageInformation();
+}
 
+/**
+ * Update method : update output information of input and write to file
+ */
+template<class TInputImage>
+void
+ImageFileWriter<TInputImage>
+::Update()
+{
+  this->UpdateOutputInformation();
+
+  this->SetAbortGenerateData(0);
+  this->SetProgress(0.0);
+
+  /**
+   * Tell all Observers that the filter is starting
+   */
+  this->InvokeEvent(itk::StartEvent());
+  
   this->UpdateProgress(0);
   m_CurrentDivision = 0;
   m_DivisionProgress = 0;
 
   // Get the source process object
+  InputImagePointer inputPtr =
+    const_cast<InputImageType *>(this->GetInput());
   itk::ProcessObject* source = inputPtr->GetSource();
   m_IsObserving = false;
   m_ObserverID = 0;
@@ -613,6 +617,12 @@ ImageFileWriter<TInputImage>
     otbLogMacro(Warning,<< "Could not get the source process object. Progress report might be buggy");
     }
 
+  /**
+   * Loop over the number of pieces, execute the upstream pipeline on each
+   * piece, and copy the results into the output image.
+   */
+  InputImageRegionType streamRegion;
+
   for (m_CurrentDivision = 0;
        m_CurrentDivision < m_NumberOfDivisions && !this->GetAbortGenerateData();
        m_CurrentDivision++, m_DivisionProgress = 0, this->UpdateFilterProgress())
diff --git a/Modules/IO/ImageIO/include/otbImageIOFactory.h b/Modules/IO/ImageIO/include/otbImageIOFactory.h
index 287c149792eada50c4605aaa1213d6e621c4807a..591999ac6ef4aa91d5017f171fad14df516bb1d7 100644
--- a/Modules/IO/ImageIO/include/otbImageIOFactory.h
+++ b/Modules/IO/ImageIO/include/otbImageIOFactory.h
@@ -23,6 +23,7 @@
 
 #include "itkObject.h"
 #include "otbImageIOBase.h"
+#include "OTBImageIOExport.h"
 
 namespace otb
 {
@@ -31,7 +32,7 @@ namespace otb
  *
  * \ingroup OTBImageIO
  */
-class ITK_EXPORT ImageIOFactory : public itk::Object
+class OTBImageIO_EXPORT ImageIOFactory : public itk::Object
 {
 public:
   /** Standard class typedefs. */
diff --git a/Modules/IO/ImageIO/include/otbMultiImageFileWriter.h b/Modules/IO/ImageIO/include/otbMultiImageFileWriter.h
new file mode 100644
index 0000000000000000000000000000000000000000..d8b05f5c4bdaa69a60606dc5a5ed4db1ef47b416
--- /dev/null
+++ b/Modules/IO/ImageIO/include/otbMultiImageFileWriter.h
@@ -0,0 +1,296 @@
+/*
+ * Copyright (C) CS SI
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbMultiImageFileWriter_h
+#define otbMultiImageFileWriter_h
+
+#include "otbImageFileWriter.h"
+#include "otbImage.h"
+#include "itkImageBase.h"
+#include "itkProcessObject.h"
+#include "itkImageIOBase.h"
+#include "OTBImageIOExport.h"
+
+#include <boost/shared_ptr.hpp>
+
+namespace otb
+{
+
+/** \class MultiImageFileWriter
+ *  \brief Streams a pipeline with multiple outputs.
+ *  It writes each output to a file. Inputs
+ *  are connected to the writer using the AddInputImage method.
+ *  The type of streaming can be chosen. Each output may have a different size.
+ *  When the user gives a number of lines per strip or a tile size, the value
+ *  is interpreted on the first input to deduce the number of streams. This
+ *  number of streams is then used to split the other inputs.
+ *
+ * \ingroup OTBImageIO
+ */
+class OTBImageIO_EXPORT MultiImageFileWriter: public itk::ProcessObject
+{
+public:
+  /** Standard class typedefs. */
+  typedef MultiImageFileWriter Self;
+  typedef itk::ProcessObject Superclass;
+  typedef itk::SmartPointer<Self> Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  itkNewMacro(Self);
+
+  itkTypeMacro(MultiImageFileWriter, itk::ProcessObject);
+
+  /** Public typedefs */
+  typedef itk::ImageBase<2> ImageBaseType;
+  typedef ImageBaseType::RegionType RegionType;
+  typedef ImageBaseType::IndexType  IndexType;
+  typedef ImageBaseType::IndexValueType  IndexValueType;
+  typedef ImageBaseType::SizeType   SizeType;
+  typedef ImageBaseType::SizeValueType   SizeValueType;
+
+  /**  Set the streaming mode to 'stripped' and configure the number of strips
+   *   which will be used to stream the image */
+  void SetNumberOfDivisionsStrippedStreaming(unsigned int nbDivisions);
+
+  /**  Set the streaming mode to 'tiled' and configure the number of tiles
+   *   which will be used to stream the image */
+  void SetNumberOfDivisionsTiledStreaming(unsigned int nbDivisions);
+
+  /**  Set the streaming mode to 'stripped' and configure the number of strips
+   *   which will be used to stream the image with respect to a number of line
+   *   per strip */
+  void SetNumberOfLinesStrippedStreaming(unsigned int nbLinesPerStrip);
+
+  /**  Set the streaming mode to 'stripped' and configure the number of MB
+   *   available. The actual number of divisions is computed automatically
+   *   by estimating the memory consumption of the pipeline.
+   *   Setting the availableRAM parameter to 0 means that the available RAM
+   *   is set from the CMake configuration option.
+   *   The bias parameter is a multiplier applied on the estimated memory size
+   *   of the pipeline and can be used to fine tune the potential gap between
+   *   estimated memory and actual memory used, which can happen because of
+   *   composite filters for example */
+  void SetAutomaticStrippedStreaming(unsigned int availableRAM = 0, double bias = 1.0);
+
+  /**  Set the streaming mode to 'tiled' and configure the dimension of the tiles
+   *   in pixels for each dimension (square tiles will be generated) */
+  void SetTileDimensionTiledStreaming(unsigned int tileDimension);
+
+  /**  Set the streaming mode to 'tiled' and configure the number of MB
+   *   available. The actual number of divisions is computed automatically
+   *   by estimating the memory consumption of the pipeline.
+   *   Tiles will be square.
+   *   Setting the availableRAM parameter to 0 means that the available RAM
+   *   is set from the CMake configuration option
+   *   The bias parameter is a multiplier applied on the estimated memory size
+   *   of the pipeline and can be used to fine tune the potential gap between
+   *   estimated memory and actual memory used, which can happen because of
+   *   composite filters for example */
+  void SetAutomaticTiledStreaming(unsigned int availableRAM = 0, double bias = 1.0);
+
+  /**  Set the streaming mode to 'adaptative' and configure the number of MB
+   *   available. The actual number of divisions is computed automatically
+   *   by estimating the memory consumption of the pipeline.
+   *   Tiles will try to match the input file tile scheme.
+   *   Setting the availableRAM parameter to 0 means that the available RAM
+   *   is set from the CMake configuration option */
+  void SetAutomaticAdaptativeStreaming(unsigned int availableRAM = 0, double bias = 1.0);
+
+  virtual void UpdateOutputData(itk::DataObject * itkNotUsed(output));
+
+  /** Connect a new input to the multi-writer. Only the input pointer is
+   *  required. If the filename list is empty,
+   *  streaming will occur without writing. It the filename list contains more
+   *  than one element, then the output will be divided into this number of
+   *  granule files. The resolution factor specifies the ratio between the height of this image and the
+   *  height of a reference image. The number of lines per strip class parameter will be modified according to this factor, so
+   *  that images with different resolutions can be streamed together. For example, use 1.0 for a 10m pixels image, 0.5 for a 20m
+   *  pixels image, and specify the number of lines per strip according to the 10m pixels image.
+   *  You may specify top and bottom margins that will be removed from the input image, according to its largest possible region.
+   */
+  template <class TImage>
+  void AddInputImage(const TImage* inputPtr, const std::string & fileName)
+  {
+    Sink<TImage> * sink = new Sink<TImage>(inputPtr, fileName);
+    m_SinkList.push_back(SinkBase::Pointer(sink));
+    unsigned int size = m_SinkList.size();
+    this->SetNthInput(size - 1, const_cast<itk::DataObject*>(dynamic_cast<const itk::DataObject*>(inputPtr)));
+  }
+
+  /** Add a new ImageFileWriter to the multi-writer. This is an alternative method
+   *  when you already have an instanciated writer.
+   */
+  template <class TWriter>
+  void AddInputWriter(const TWriter* writer)
+  {
+    Sink<typename TWriter::InputImageType > * sink =
+      new Sink<typename TWriter::InputImageType >(writer);
+    m_SinkList.push_back(SinkBase::Pointer(sink));
+    unsigned int size = m_SinkList.size();
+    this->SetNthInput(size - 1, const_cast<itk::DataObject*>(dynamic_cast<const itk::DataObject*>(writer->GetInput())));
+  }
+
+  virtual void UpdateOutputInformation();
+
+  virtual void Update()
+  {
+    this->UpdateOutputInformation();
+    this->UpdateOutputData(NULL);
+  }
+
+protected:
+  /** SetInput is changed to protected. Use AddInputImage to connect the
+   *  pipeline to the writer
+   */
+  virtual void SetInput(const itk::ProcessObject::DataObjectIdentifierType & key, itk::DataObject* image)
+    { this->Superclass::SetInput(key, image); }
+
+  /** SetNthInput is changed to protected. Use AddInputImage to connect the
+   *  pipeline to the writer
+   */
+  virtual void SetNthInput(itk::ProcessObject::DataObjectPointerArraySizeType i, itk::DataObject* image)
+    { this->Superclass::SetNthInput(i, image); }
+
+  MultiImageFileWriter();
+  virtual ~MultiImageFileWriter() {}
+
+  /** GenerateData calls the Write method for each connected input */
+  virtual void GenerateData(void);
+
+  /** GenerateInputRequestedRegion can predict approximate input regions
+   *  based on the requested region of the fake output. Only usefull for
+   *  pipeline memory print estimation
+   */
+  virtual void GenerateInputRequestedRegion();
+
+  /** Computes the number of divisions */
+  virtual void InitializeStreaming();
+
+  /** Goes up the pipeline starting at imagePtr, resetting all requested regions
+   *  to a null region. This may be usefull when mixing inputs of different
+   *  resolutions. */
+  void ResetAllRequestedRegions(ImageBaseType* imagePtr);
+
+  /** Returns the current stream region of the given input */
+  virtual RegionType GetStreamRegion(int inputIndex);
+
+  void operator =(const MultiImageFileWriter&); //purposely not implemented
+
+  void ObserveSourceFilterProgress(itk::Object* object, const itk::EventObject & event)
+  {
+    if (typeid(event) != typeid(itk::ProgressEvent))
+      {
+      return;
+      }
+
+    itk::ProcessObject* processObject = dynamic_cast<itk::ProcessObject*>(object);
+    if (processObject)
+      {
+      m_DivisionProgress = processObject->GetProgress();
+      }
+
+    this->UpdateFilterProgress();
+  }
+
+  void UpdateFilterProgress()
+  {
+    this->UpdateProgress((m_DivisionProgress + m_CurrentDivision) / m_NumberOfDivisions);
+  }
+
+private:
+  typedef otb::Image<unsigned char, 2> FakeOutputType;
+  typedef StreamingManager<FakeOutputType>        StreamingManagerType;
+  /** Streaming manager used for the N inputs */
+  StreamingManagerType::Pointer m_StreamingManager;
+
+  //Division parameters
+  unsigned int m_NumberOfDivisions;
+  unsigned int m_CurrentDivision;
+  float m_DivisionProgress;
+
+  bool m_IsObserving;
+  unsigned long m_ObserverID;
+
+  /** \class SinkBase
+   * Internal base wrapper class to handle each ImageFileWriter
+   *
+   * \ingroup OTBImageIO
+   */
+  class SinkBase
+  {
+  public:
+    SinkBase() {}
+    SinkBase(ImageBaseType::ConstPointer inputImage) :
+      m_InputImage(inputImage)
+    {}
+    virtual ~SinkBase() {}
+    virtual ImageBaseType::ConstPointer GetInput() const { return m_InputImage; }
+    virtual ImageBaseType::Pointer GetInput() { return const_cast<ImageBaseType*>(m_InputImage.GetPointer()); }
+    virtual void WriteImageInformation() = 0;
+    virtual void Write(const RegionType & streamRegion) = 0;
+    virtual bool CanStreamWrite() = 0;
+    typedef boost::shared_ptr<SinkBase> Pointer;
+  protected:
+    /** The image on which streaming is performed */
+    ImageBaseType::ConstPointer m_InputImage;
+  };
+
+  /** \class Sink
+   *  Wrapper class for each ImageFileWriter
+   *
+   * \ingroup OTBImageIO
+   */
+  template <class TImage>
+  class Sink : public SinkBase
+  {
+  public:
+    Sink() {}
+    Sink(typename TImage::ConstPointer inputImage,
+         const std::string & filename);
+    Sink(typename otb::ImageFileWriter<TImage>::ConstPointer writer);
+
+    virtual ~Sink() {}
+
+    virtual void WriteImageInformation();
+    virtual void Write(const RegionType & streamRegion);
+    virtual bool CanStreamWrite();
+    typedef boost::shared_ptr<Sink> Pointer;
+  private:
+    /** Actual writer for this image */
+    typename otb::ImageFileWriter<TImage>::Pointer m_Writer;
+
+    /** An ImageIO used to actually write data to a file */
+    otb::ImageIOBase::Pointer m_ImageIO;
+  };
+
+  /** The list of inputs and their associated parameters, built using AddInput */
+  typedef std::vector<boost::shared_ptr<SinkBase> > SinkListType;
+  SinkListType m_SinkList;
+
+  std::vector<RegionType> m_StreamRegionList;
+};
+
+} // end of namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbMultiImageFileWriter.txx"
+#endif
+
+#endif // otbMultiImageFileWriter_h
diff --git a/Modules/IO/ImageIO/include/otbMultiImageFileWriter.txx b/Modules/IO/ImageIO/include/otbMultiImageFileWriter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..5e6d6304dc52aa3dc95a497b45264ff2cc15d277
--- /dev/null
+++ b/Modules/IO/ImageIO/include/otbMultiImageFileWriter.txx
@@ -0,0 +1,89 @@
+/*
+ * Copyright (C) CS SI
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbMultiImageFileWriter_txx
+#define otbMultiImageFileWriter_txx
+
+#include "otbMultiImageFileWriter.h"
+#include "otbImageIOFactory.h"
+#include "otbMacro.h"
+
+namespace otb
+{
+
+template <class TImage>
+MultiImageFileWriter::Sink<TImage>
+::Sink(typename TImage::ConstPointer inputImage,
+       const std::string & fileName):
+  SinkBase(dynamic_cast<const ImageBaseType*>(inputImage.GetPointer())),
+  m_Writer(otb::ImageFileWriter<TImage>::New()),
+  m_ImageIO(NULL)
+{
+  m_Writer->SetFileName(fileName);
+  m_Writer->SetInput(inputImage);
+}
+
+template <class TImage>
+MultiImageFileWriter::Sink<TImage>
+::Sink(typename otb::ImageFileWriter<TImage>::ConstPointer writer):
+  SinkBase(dynamic_cast<const ImageBaseType*>(writer->GetInput()->GetPointer())),
+  m_Writer(writer),
+  m_ImageIO(NULL)
+{
+}
+
+template <class TImage>
+bool
+MultiImageFileWriter::Sink<TImage>
+::CanStreamWrite()
+{
+  if (m_ImageIO.IsNull())
+    return false;
+  return m_ImageIO->CanStreamWrite();
+}
+
+template <class TImage>
+void
+MultiImageFileWriter::Sink<TImage>
+::WriteImageInformation()
+{
+  m_Writer->UpdateOutputInformation();
+  m_ImageIO = m_Writer->GetImageIO();
+}
+
+template <class TImage>
+void
+MultiImageFileWriter::Sink<TImage>
+::Write(const RegionType & streamRegion)
+{
+  // Write the image stream
+  itk::ImageIORegion ioRegion(TImage::ImageDimension);
+  for (unsigned int i = 0; i < TImage::ImageDimension; ++i)
+    {
+    ioRegion.SetSize(i, streamRegion.GetSize(i));
+    ioRegion.SetIndex(i, streamRegion.GetIndex(i));
+    }
+  m_ImageIO->SetIORegion(ioRegion);
+  m_Writer->UpdateOutputData(nullptr);
+}
+
+} // end of namespace otb
+
+#endif // otbMultiImageFileWriter_txx
diff --git a/Modules/IO/ImageIO/otb-module.cmake b/Modules/IO/ImageIO/otb-module.cmake
index 552d431d5a41aba8c23ba62004eac3e60a905c80..d592ba8b086549b5a132f64f0867282e9b7a4de8 100644
--- a/Modules/IO/ImageIO/otb-module.cmake
+++ b/Modules/IO/ImageIO/otb-module.cmake
@@ -22,6 +22,7 @@ set(DOCUMENTATION "This module contains classes related to the reading and the
 writing of remote sensing images.")
 
 otb_module(OTBImageIO
+  ENABLE_SHARED
   DEPENDS
     OTBBoostAdapters
     OTBCommon
diff --git a/Modules/IO/ImageIO/src/CMakeLists.txt b/Modules/IO/ImageIO/src/CMakeLists.txt
index 2d580c6541b14e92fb19e500276f5146103257e0..e604942a7cd4b8473e595bd457dfbdac05523938 100644
--- a/Modules/IO/ImageIO/src/CMakeLists.txt
+++ b/Modules/IO/ImageIO/src/CMakeLists.txt
@@ -20,6 +20,7 @@
 
 set(OTBImageIO_SRC
   otbImageIOFactory.cxx
+  otbMultiImageFileWriter.cxx
   )
 
 add_library(OTBImageIO ${OTBImageIO_SRC})
diff --git a/Modules/IO/ImageIO/src/otbMultiImageFileWriter.cxx b/Modules/IO/ImageIO/src/otbMultiImageFileWriter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a3f32a8d36bd404a271c2fd58d61b94488faf981
--- /dev/null
+++ b/Modules/IO/ImageIO/src/otbMultiImageFileWriter.cxx
@@ -0,0 +1,443 @@
+/*
+ * Copyright (C) CS SI
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "otbMultiImageFileWriter.h"
+#include "otbImageIOFactory.h"
+
+namespace otb
+{
+
+MultiImageFileWriter
+::MultiImageFileWriter() :
+ m_NumberOfDivisions(0),
+ m_CurrentDivision(0),
+ m_DivisionProgress(0.0),
+ m_IsObserving(true),
+ m_ObserverID(0)
+{
+  // By default, we use tiled streaming, with automatic tile size
+  // We don't set any parameter, so the memory size is retrieved from the OTB configuration options
+  this->SetAutomaticAdaptativeStreaming();
+  // add a fake output to drive memory estimation
+  this->SetNthOutput(0, FakeOutputType::New());
+}
+
+void
+MultiImageFileWriter
+::SetNumberOfDivisionsStrippedStreaming(unsigned int nbDivisions)
+{
+  typedef NumberOfDivisionsStrippedStreamingManager<FakeOutputType> NumberOfDivisionsStrippedStreamingManagerType;
+  NumberOfDivisionsStrippedStreamingManagerType::Pointer streamingManager =
+    NumberOfDivisionsStrippedStreamingManagerType::New();
+  streamingManager->SetNumberOfDivisions(nbDivisions);
+
+  m_StreamingManager = streamingManager;
+}
+
+void
+MultiImageFileWriter
+::SetNumberOfDivisionsTiledStreaming(unsigned int nbDivisions)
+{
+  typedef NumberOfDivisionsTiledStreamingManager<FakeOutputType> NumberOfDivisionsTiledStreamingManagerType;
+  NumberOfDivisionsTiledStreamingManagerType::Pointer streamingManager =
+    NumberOfDivisionsTiledStreamingManagerType::New();
+  streamingManager->SetNumberOfDivisions(nbDivisions);
+
+  m_StreamingManager = streamingManager;
+}
+
+void
+MultiImageFileWriter
+::SetNumberOfLinesStrippedStreaming(unsigned int nbLinesPerStrip)
+{
+  typedef NumberOfLinesStrippedStreamingManager<FakeOutputType> NumberOfLinesStrippedStreamingManagerType;
+  NumberOfLinesStrippedStreamingManagerType::Pointer streamingManager =
+    NumberOfLinesStrippedStreamingManagerType::New();
+  streamingManager->SetNumberOfLinesPerStrip(nbLinesPerStrip);
+
+  m_StreamingManager = streamingManager;
+}
+
+void
+MultiImageFileWriter
+::SetAutomaticStrippedStreaming(unsigned int availableRAM, double bias)
+{
+  typedef RAMDrivenStrippedStreamingManager<FakeOutputType> RAMDrivenStrippedStreamingManagerType;
+  RAMDrivenStrippedStreamingManagerType::Pointer streamingManager =
+    RAMDrivenStrippedStreamingManagerType::New();
+  streamingManager->SetAvailableRAMInMB(availableRAM);
+  streamingManager->SetBias(bias);
+
+  m_StreamingManager = streamingManager;
+}
+
+void
+MultiImageFileWriter
+::SetTileDimensionTiledStreaming(unsigned int tileDimension)
+{
+  typedef TileDimensionTiledStreamingManager<FakeOutputType> TileDimensionTiledStreamingManagerType;
+  TileDimensionTiledStreamingManagerType::Pointer streamingManager =
+    TileDimensionTiledStreamingManagerType::New();
+  streamingManager->SetTileDimension(tileDimension);
+
+  m_StreamingManager = streamingManager;
+}
+
+void
+MultiImageFileWriter
+::SetAutomaticTiledStreaming(unsigned int availableRAM, double bias)
+{
+  typedef RAMDrivenTiledStreamingManager<FakeOutputType> RAMDrivenTiledStreamingManagerType;
+  RAMDrivenTiledStreamingManagerType::Pointer streamingManager =
+    RAMDrivenTiledStreamingManagerType::New();
+  streamingManager->SetAvailableRAMInMB(availableRAM);
+  streamingManager->SetBias(bias);
+  m_StreamingManager = streamingManager;
+}
+
+void
+MultiImageFileWriter
+::SetAutomaticAdaptativeStreaming(unsigned int availableRAM, double bias)
+{
+  typedef RAMDrivenAdaptativeStreamingManager<FakeOutputType> RAMDrivenAdaptativeStreamingManagerType;
+  RAMDrivenAdaptativeStreamingManagerType::Pointer streamingManager =
+    RAMDrivenAdaptativeStreamingManagerType::New();
+  streamingManager->SetAvailableRAMInMB(availableRAM);
+  streamingManager->SetBias(bias);
+  m_StreamingManager = streamingManager;
+}
+
+void
+MultiImageFileWriter
+::InitializeStreaming()
+{
+//  const ImageBaseType* inputPtr = this->GetInput(0);
+  if(m_SinkList.size() == 0)
+    itkExceptionMacro("At least one input must be connected to the writer\n");
+  const ImageBaseType* inputPtr = m_SinkList[0]->GetInput();
+  if(!inputPtr)
+    itkExceptionMacro("At least one input must be connected to the writer\n");
+
+  /** Control if the ImageIO is CanStreamWrite */
+  m_NumberOfDivisions = 1;
+  bool canStream = true;
+  bool isBuffered = true;
+  for (unsigned int inputIndex = 0; inputIndex < m_SinkList.size(); ++inputIndex)
+    {
+    if (!m_SinkList[inputIndex]->CanStreamWrite())
+      {
+      canStream = false;
+      }
+    if (m_SinkList[inputIndex]->GetInput()->GetBufferedRegion() !=
+        m_SinkList[inputIndex]->GetInput()->GetLargestPossibleRegion())
+      {
+      isBuffered = false;
+      }
+    }
+  if (canStream == false)
+    {
+    otbWarningMacro(
+      << "One of the selected ImageIO does not support streaming.");
+    this->SetNumberOfDivisionsStrippedStreaming(m_NumberOfDivisions);
+    }
+
+  /** Compare the buffered region  with the inputRegion which is the largest
+  * possible region or a user defined region through extended filename
+  * Not sure that if this modification is needed  */
+  else if (isBuffered)
+    {
+    otbMsgDevMacro(<< "Buffered region is the largest possible region, there is"
+      " no need for streaming.");
+    this->SetNumberOfDivisionsStrippedStreaming(m_NumberOfDivisions);
+    }
+  else
+    {
+    /**
+     * Determine of number of pieces to divide the input.  This will be the
+     * first estimated on the fake output, which has the same size as the
+     * first input. Then there is a check that each input can be split into
+     * this number of pieces.
+     */
+    FakeOutputType * fakeOut = static_cast<FakeOutputType *>(
+      this->itk::ProcessObject::GetOutput(0));
+    RegionType region = fakeOut->GetLargestPossibleRegion();
+    m_StreamingManager->PrepareStreaming(fakeOut, region);
+    m_NumberOfDivisions = m_StreamingManager->GetNumberOfSplits();
+    // Check this number of division is compatible with all inputs
+    bool nbDivValid = false;
+    while ( (!nbDivValid) && 1 < m_NumberOfDivisions)
+      {
+      unsigned int smallestNbDiv = m_NumberOfDivisions;
+      for (unsigned int i = 0; i < m_SinkList.size(); ++i)
+        {
+        unsigned int div = m_StreamingManager->GetSplitter()->GetNumberOfSplits(
+          m_SinkList[i]->GetInput()->GetLargestPossibleRegion(),
+          m_NumberOfDivisions);
+        smallestNbDiv = std::min(div, smallestNbDiv);
+        }
+      if (smallestNbDiv == m_NumberOfDivisions)
+        {
+        nbDivValid = true;
+        }
+      else
+        {
+        m_NumberOfDivisions = smallestNbDiv;
+        }
+      }
+    if (m_NumberOfDivisions == 1)
+      {
+      otbWarningMacro("Can't find a common split scheme between all inputs, streaming disabled\n");
+      }
+    otbMsgDebugMacro(<< "Number Of Stream Divisions : " << m_NumberOfDivisions);
+    }
+}
+
+void
+MultiImageFileWriter
+::ResetAllRequestedRegions(ImageBaseType* imagePtr)
+{
+  RegionType nullRegion = imagePtr->GetLargestPossibleRegion();
+  nullRegion.SetSize(0, 0);
+  nullRegion.SetSize(1, 0);
+
+  imagePtr->SetRequestedRegion(nullRegion);
+  if(imagePtr->GetSource())
+    {
+    itk::ProcessObject::DataObjectPointerArray inputs = imagePtr->GetSource()->GetInputs();
+    for( itk::ProcessObject::DataObjectPointerArray::iterator
+        it = inputs.begin();
+        it != inputs.end();
+        ++it )
+      {
+      ImageBaseType * inputImagePtr = dynamic_cast<ImageBaseType*>(it->GetPointer());
+      if(inputImagePtr != NULL)
+        {
+        ResetAllRequestedRegions(inputImagePtr);
+        }
+      }
+    }
+}
+
+void
+MultiImageFileWriter
+::UpdateOutputInformation()
+{
+  for(unsigned int inputIndex = 0; inputIndex < m_SinkList.size(); ++inputIndex)
+    {
+    m_SinkList[inputIndex]->WriteImageInformation();
+    }
+  this->GenerateOutputInformation();
+}
+
+void
+MultiImageFileWriter
+::UpdateOutputData(itk::DataObject * itkNotUsed(output))
+{
+  /**
+   * prevent chasing our tail
+   */
+  if (this->m_Updating)
+    {
+    return;
+    }
+
+  // Initialize streaming
+  this->InitializeStreaming();
+
+  this->SetAbortGenerateData(0);
+  this->SetProgress(0.0);
+  this->m_Updating = true;
+
+  /**
+   * Tell all Observers that the filter is starting
+   */
+  this->InvokeEvent(itk::StartEvent());
+
+  this->UpdateProgress(0);
+  m_CurrentDivision = 0;
+  m_DivisionProgress = 0;
+
+  /** Loop over the number of pieces, set and propagate requested regions then
+   * update pipeline upstream
+   */
+  int numInputs = m_SinkList.size(); //this->GetNumberOfInputs();
+  m_StreamRegionList.resize(numInputs);
+  itkDebugMacro( "Number of streaming divisions: " << m_NumberOfDivisions);
+
+  // Add observer only to first input
+  if(numInputs > 0)
+    {
+    m_IsObserving = false;
+    m_ObserverID = 0;
+
+    typedef itk::MemberCommand<Self> CommandType;
+    typedef CommandType::Pointer CommandPointerType;
+
+    CommandPointerType command = CommandType::New();
+    command->SetCallbackFunction(this, &Self::ObserveSourceFilterProgress);
+
+    itk::ProcessObject* src = this->GetInput(0)->GetSource();
+    m_ObserverID = src->AddObserver(itk::ProgressEvent(), command);
+    m_IsObserving = true;
+    }
+
+  for (m_CurrentDivision = 0; m_CurrentDivision < m_NumberOfDivisions && !this->GetAbortGenerateData();
+      m_CurrentDivision++, m_DivisionProgress = 0, this->UpdateFilterProgress())
+    {
+    // Update all stream regions
+    for(int inputIndex = 0; inputIndex < numInputs; ++inputIndex)
+      {
+      m_StreamRegionList[inputIndex] = GetStreamRegion(inputIndex);
+      }
+
+    // NOTE : this reset was probably designed to work with the next section
+    // Where the final requested region is the "union" between the computed
+    // requested region and the current requested region.
+    
+    // Reset requested regions for all images
+    for(int inputIndex = 0; inputIndex < numInputs; ++inputIndex)
+      {
+      ResetAllRequestedRegions(m_SinkList[inputIndex]->GetInput());
+      }
+
+    for(int inputIndex = 0; inputIndex < numInputs; ++inputIndex)
+      {
+      ImageBaseType::Pointer inputPtr = m_SinkList[inputIndex]->GetInput();
+      RegionType inputRequestedRegion = m_StreamRegionList[inputIndex];
+      const RegionType & currentInputRequestedRegion = inputPtr->GetRequestedRegion();
+      if( currentInputRequestedRegion != inputPtr->GetLargestPossibleRegion()
+        && currentInputRequestedRegion.GetNumberOfPixels() != 0)
+        {
+        IndexType startIndex = currentInputRequestedRegion.GetIndex();
+        IndexType lastIndex = currentInputRequestedRegion.GetUpperIndex();
+        startIndex[0] = std::min(startIndex[0], inputRequestedRegion.GetIndex(0));
+        startIndex[1] = std::min(startIndex[1], inputRequestedRegion.GetIndex(1));
+        lastIndex[0] = std::max(lastIndex[0], inputRequestedRegion.GetUpperIndex()[0]);
+        lastIndex[1] = std::max(lastIndex[1], inputRequestedRegion.GetUpperIndex()[1]);
+        inputRequestedRegion.SetIndex(startIndex);
+        inputRequestedRegion.SetUpperIndex(lastIndex);
+        }
+
+      inputPtr->SetRequestedRegion(inputRequestedRegion);
+      inputPtr->PropagateRequestedRegion();
+      }
+
+    /** Call GenerateData to write streams to files if needed */
+    this->GenerateData();
+    }
+
+  /**
+   * If we ended due to aborting, push the progress up to 1.0 (since
+   * it probably didn't end there)
+   */
+  if (!this->GetAbortGenerateData())
+    {
+    this->UpdateProgress(1.0);
+    }
+
+  // Notify end event observers
+  this->InvokeEvent(itk::EndEvent());
+
+  if (m_IsObserving)
+    {
+    ImageBaseType::Pointer inputPtr = m_SinkList[0]->GetInput(); // const_cast<ImageBaseType *>(this->GetInput(0));
+    itk::ProcessObject* source = inputPtr->GetSource();
+    m_IsObserving = false;
+    source->RemoveObserver(m_ObserverID);
+    }
+
+  /**
+   * Release any inputs if marked for release
+   */
+  this->ReleaseInputs();
+
+  // Mark that we are no longer updating the data in this filter
+  this->m_Updating = false;
+}
+
+
+void
+MultiImageFileWriter
+::GenerateInputRequestedRegion()
+{
+  Superclass::GenerateInputRequestedRegion();
+
+  // Approximate conversion of output requested region into each input,
+  // but this is only to have a consistent pipeline memory estimation.
+  int numInputs = m_SinkList.size(); //this->GetNumberOfInputs();
+
+  FakeOutputType* fakeOut = static_cast<FakeOutputType *>(
+    this->itk::ProcessObject::GetOutput(0));
+
+  if (numInputs)
+    {
+    RegionType refLargest = fakeOut->GetLargestPossibleRegion();
+    RegionType refRequest = fakeOut->GetRequestedRegion();
+    IndexType idxLargest = refLargest.GetIndex();
+    SizeType sizeLargest = refLargest.GetSize();
+    IndexType idxRequest = refRequest.GetIndex();
+    SizeType sizeRequest = refRequest.GetSize();
+    for (int i = 0; i < numInputs; ++i)
+      {
+      ImageBaseType* inputPtr = m_SinkList[i]->GetInput();
+      if(!inputPtr)
+        {
+        return;
+        }
+      RegionType region = inputPtr->GetLargestPossibleRegion();
+      IndexType idx = region.GetIndex();
+      SizeType size = region.GetSize();
+      idx[0] += size[0] * (idxRequest[0] - idxLargest[0]) / sizeLargest[0];
+      idx[1] += size[1] * (idxRequest[1] - idxLargest[1]) / sizeLargest[1];
+      size[0] *= sizeRequest[0] / sizeLargest[0];
+      size[1] *= sizeRequest[1] / sizeLargest[1];
+      region.SetIndex(idx);
+      region.SetSize(size);
+      inputPtr->SetRequestedRegion(region);
+      }
+    }
+}
+
+void
+MultiImageFileWriter
+::GenerateData()
+{
+  int numInputs = m_SinkList.size();
+  for(int inputIndex = 0; inputIndex < numInputs; ++inputIndex)
+    {
+    m_SinkList[inputIndex]->Write(m_StreamRegionList[inputIndex]);
+    }
+}
+
+MultiImageFileWriter::RegionType
+MultiImageFileWriter
+::GetStreamRegion(int inputIndex)
+{
+  const SinkBase::Pointer sink = m_SinkList[inputIndex];
+  RegionType region = sink->GetInput()->GetLargestPossibleRegion();
+
+  m_StreamingManager->GetSplitter()->GetSplit(
+    m_CurrentDivision,
+    m_NumberOfDivisions,
+    region);
+  return region;
+}
+
+} // end of namespace otb
diff --git a/Modules/IO/ImageIO/test/CMakeLists.txt b/Modules/IO/ImageIO/test/CMakeLists.txt
index aac49fcdeab48d245227689b5aeb6461f6d6dd07..110b25410e8bfd8f663a505fa202b4399d213609 100644
--- a/Modules/IO/ImageIO/test/CMakeLists.txt
+++ b/Modules/IO/ImageIO/test/CMakeLists.txt
@@ -76,6 +76,7 @@ otbImageIOFactoryNew.cxx
 otbCompareWritingComplexImage.cxx
 otbImageFileReaderOptBandTest.cxx
 otbImageFileWriterOptBandTest.cxx
+otbMultiImageFileWriterTest.cxx
 )
 
 add_executable(otbImageIOTestDriver ${OTBImageIOTests})
@@ -1339,3 +1340,31 @@ otb_add_test(NAME ioTvImageIOToWriterOptions_OptBandReorgTest COMMAND otbImageIO
   ${TEMP}/QB_Toulouse_Ortho_XS_WriterOptBandReorg.tif?bands=2,:,-3,2:-1
   4
   )
+
+otb_add_test(NAME ioTvMultiImageFileWriter_SameSize
+  COMMAND otbImageIOTestDriver
+  --compare-n-images ${EPSILON_9} 2
+  ${INPUTDATA}/GomaAvant.png
+  ${TEMP}/ioTvMultiImageFileWriter_SameSize1.tif
+  ${INPUTDATA}/GomaApres.png
+  ${TEMP}/ioTvMultiImageFileWriter_SameSize2.tif
+  otbMultiImageFileWriterTest
+  ${INPUTDATA}/GomaAvant.png
+  ${INPUTDATA}/GomaApres.png
+  ${TEMP}/ioTvMultiImageFileWriter_SameSize1.tif
+  ${TEMP}/ioTvMultiImageFileWriter_SameSize2.tif
+  50)
+
+otb_add_test(NAME ioTvMultiImageFileWriter_DiffSize
+  COMMAND otbImageIOTestDriver
+  --compare-n-images ${EPSILON_9} 2
+  ${INPUTDATA}/GomaAvant.png
+  ${TEMP}/ioTvMultiImageFileWriter_DiffSize1.tif
+  ${INPUTDATA}/QB_Toulouse_Ortho_PAN.tif
+  ${TEMP}/ioTvMultiImageFileWriter_DiffSize2.tif
+  otbMultiImageFileWriterTest
+  ${INPUTDATA}/GomaAvant.png
+  ${INPUTDATA}/QB_Toulouse_Ortho_PAN.tif
+  ${TEMP}/ioTvMultiImageFileWriter_DiffSize1.tif
+  ${TEMP}/ioTvMultiImageFileWriter_DiffSize2.tif
+  25)
diff --git a/Modules/IO/ImageIO/test/otbImageIOTestDriver.cxx b/Modules/IO/ImageIO/test/otbImageIOTestDriver.cxx
index fb2aa583aa73653ece1ffb9f523ae6caaec71d84..4a7ff9409a53d34ec7d404c6f6d56cc57c8917bf 100644
--- a/Modules/IO/ImageIO/test/otbImageIOTestDriver.cxx
+++ b/Modules/IO/ImageIO/test/otbImageIOTestDriver.cxx
@@ -154,4 +154,5 @@ void RegisterTests()
   REGISTER_TEST(otbCompareWritingComplexImageTest);
   REGISTER_TEST(otbImageFileReaderOptBandTest);
   REGISTER_TEST(otbImageFileWriterOptBandTest);
+  REGISTER_TEST(otbMultiImageFileWriterTest);
 }
diff --git a/Modules/IO/ImageIO/test/otbMultiImageFileWriterTest.cxx b/Modules/IO/ImageIO/test/otbMultiImageFileWriterTest.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..77d6ede243b19c34694ff18a74959457c060a0d0
--- /dev/null
+++ b/Modules/IO/ImageIO/test/otbMultiImageFileWriterTest.cxx
@@ -0,0 +1,65 @@
+/*
+ * Copyright (C) CS SI
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "otbMultiImageFileWriter.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+
+int otbMultiImageFileWriterTest(int argc, char* argv[])
+{
+  typedef unsigned short PixelType1;
+  typedef otb::Image<PixelType1, 2> ImageType1;
+  typedef otb::ImageFileReader<ImageType1> ReaderType1;
+
+  typedef double PixelType2;
+  typedef otb::Image<PixelType2, 2> ImageType2;
+  typedef otb::ImageFileReader<ImageType2> ReaderType2;
+
+  typedef otb::MultiImageFileWriter WriterType;
+
+  if (argc < 6)
+    {
+    std::cout << "Usage: " << argv[0] << " inputImageFileName1 inputImageFileName2 outputImageFileName1 outputImageFileName2 numberOfLinesPerStrip\n";
+    return EXIT_FAILURE;
+    }
+
+  const char * inputImageFileName1 = argv[1];
+  const char * inputImageFileName2 = argv[2];
+  const std::string outputImageFileName1 = argv[3];
+  const std::string outputImageFileName2 = argv[4];
+  const int numberOfLinesPerStrip = atoi(argv[5]);
+
+  ReaderType1::Pointer reader1 = ReaderType1::New();
+  reader1->SetFileName( inputImageFileName1 );
+
+  ReaderType2::Pointer reader2 = ReaderType2::New();
+  reader2->SetFileName( inputImageFileName2 );
+
+  WriterType::Pointer writer = WriterType::New();
+  writer->AddInputImage( reader1->GetOutput(), outputImageFileName1);
+  writer->AddInputImage( reader2->GetOutput(), outputImageFileName2);
+  writer->SetNumberOfLinesStrippedStreaming( numberOfLinesPerStrip );
+
+  writer->Update();
+
+  std::cout << writer << std::endl;
+
+  return EXIT_SUCCESS;
+}
diff --git a/Modules/Learning/DimensionalityReductionLearning/include/otbAutoencoderModel.h b/Modules/Learning/DimensionalityReductionLearning/include/otbAutoencoderModel.h
index 7bf3231a9b8d372473fd5ed0f598d3f2fd7c7c06..5fe7ec2f27c9f91f21a2ce2417502f5d4beac54d 100644
--- a/Modules/Learning/DimensionalityReductionLearning/include/otbAutoencoderModel.h
+++ b/Modules/Learning/DimensionalityReductionLearning/include/otbAutoencoderModel.h
@@ -33,8 +33,9 @@
 #endif
 #include "otb_shark.h"
 #include <shark/Algorithms/StoppingCriteria/AbstractStoppingCriterion.h>
-#include <shark/Models/FFNet.h>
-#include <shark/Models/Autoencoder.h>
+#include <shark/Models/LinearModel.h>
+#include <shark/Models/ConcatenatedModel.h>
+#include <shark/Models/NeuronLayers.h>
 #if defined(__GNUC__) || defined(__clang__)
 #pragma GCC diagnostic pop
 #endif
@@ -76,9 +77,9 @@ public:
   typedef typename Superclass::ConfidenceListSampleType         ConfidenceListSampleType;
 
   /// Neural network related typedefs
-  typedef shark::Autoencoder<NeuronType,shark::LinearNeuron> OutAutoencoderType;
-  typedef shark::Autoencoder<NeuronType,NeuronType> AutoencoderType;
-  typedef shark::FFNet<NeuronType,shark::LinearNeuron> NetworkType;
+  typedef shark::ConcatenatedModel<shark::RealVector> ModelType;
+  typedef shark::LinearModel<shark::RealVector,NeuronType> LayerType;
+  typedef shark::LinearModel<shark::RealVector, shark::LinearNeuron> OutLayerType;
 
   itkNewMacro(Self);
   itkTypeMacro(AutoencoderModel, DimensionalityReductionModel);
@@ -127,18 +128,16 @@ public:
 
   void Train() override;
 
-  template <class T, class Autoencoder>
+  template <class T>
   void TrainOneLayer(
     shark::AbstractStoppingCriterion<T> & criterion,
-    Autoencoder &,
     unsigned int,
     shark::Data<shark::RealVector> &,
     std::ostream&);
 
-  template <class T, class Autoencoder>
+  template <class T>
   void TrainOneSparseLayer(
     shark::AbstractStoppingCriterion<T> & criterion,
-    Autoencoder &,
     unsigned int,
     shark::Data<shark::RealVector> &,
     std::ostream&);
@@ -166,7 +165,9 @@ protected:
 
 private:
   /** Internal Network */
-  NetworkType m_Net;
+  ModelType m_Encoder;
+  std::vector<LayerType> m_InLayers;
+  OutLayerType m_OutLayer;
   itk::Array<unsigned int> m_NumberOfHiddenNeurons;
   /** Training parameters */
   unsigned int m_NumberOfIterations; // stop the training after a fixed number of iterations
diff --git a/Modules/Learning/DimensionalityReductionLearning/include/otbAutoencoderModel.txx b/Modules/Learning/DimensionalityReductionLearning/include/otbAutoencoderModel.txx
index 33f1c28e247c43f80ac28a1d608b1c15967c6a5e..e5a26e9ee3dc8cbf4918222f4b7b45bc93e925cb 100644
--- a/Modules/Learning/DimensionalityReductionLearning/include/otbAutoencoderModel.txx
+++ b/Modules/Learning/DimensionalityReductionLearning/include/otbAutoencoderModel.txx
@@ -34,18 +34,17 @@
 #include "otbSharkUtils.h"
 //include train function
 #include <shark/ObjectiveFunctions/ErrorFunction.h>
-#include <shark/ObjectiveFunctions/SparseAutoencoderError.h>//the error function performing the regularisation of the hidden neurons
+//~ #include <shark/ObjectiveFunctions/SparseAutoencoderError.h>//the error function performing the regularisation of the hidden neurons
 
 #include <shark/Algorithms/GradientDescent/Rprop.h>// the RProp optimization algorithm
 #include <shark/ObjectiveFunctions/Loss/SquaredLoss.h> // squared loss used for regression
 #include <shark/ObjectiveFunctions/Regularizer.h> //L2 regulariziation
-#include <shark/Models/ImpulseNoiseModel.h> //noise source to corrupt the inputs
-#include <shark/Models/ConcatenatedModel.h>//to concatenate the noise with the model
+//~ #include <shark/Models/ImpulseNoiseModel.h> //noise source to corrupt the inputs
 
 #include <shark/Algorithms/StoppingCriteria/MaxIterations.h> //A simple stopping criterion that stops after a fixed number of iterations
 #include <shark/Algorithms/StoppingCriteria/TrainingProgress.h> //Stops when the algorithm seems to converge, Tracks the progress of the training error over a period of time
 
-#include <shark/Algorithms/GradientDescent/SteepestDescent.h>
+#include <shark/Algorithms/GradientDescent/Adam.h>
 #if defined(__GNUC__) || defined(__clang__)
 #pragma GCC diagnostic pop
 #endif
@@ -83,96 +82,56 @@ AutoencoderModel<TInputValue,NeuronType>
     }
 
   // Initialization of the feed forward neural network
-  std::vector<size_t> layers;
-  layers.push_back(shark::dataDimension(inputSamples));
+  m_Encoder = ModelType();
+  m_InLayers.clear();
+  size_t previousShape = shark::dataDimension(inputSamples);
   for (unsigned int i = 0 ; i < m_NumberOfHiddenNeurons.Size(); ++i)
     {
-    layers.push_back(m_NumberOfHiddenNeurons[i]);
+    m_InLayers.push_back( LayerType(previousShape, m_NumberOfHiddenNeurons[i]) );
+    previousShape = m_NumberOfHiddenNeurons[i];
+    m_Encoder.add(&(m_InLayers.back()), true);
     }
-
   for (unsigned int i = std::max(0,static_cast<int>(m_NumberOfHiddenNeurons.Size()-1)) ; i > 0; --i)
     {
-    layers.push_back(m_NumberOfHiddenNeurons[i-1]);
-    }
-
-  layers.push_back(shark::dataDimension(inputSamples));
-  m_Net.setStructure(layers);
-  shark::initRandomNormal(m_Net,0.1);
-
-  // Training of the first Autoencoder (first and last layer of the FF network)
-  if (m_Epsilon > 0)
-    {
-    shark::TrainingProgress<> criterion(5,m_Epsilon);
-
-    OutAutoencoderType net;
-    // Shark doesn't allow to train a layer using a sparsity term AND a noisy input. 
-    if (m_Noise[0] != 0)
-      {
-      TrainOneLayer(criterion, net, 0, inputSamples, ofs);
-      }
-    else
-      {
-      TrainOneSparseLayer(criterion, net, 0, inputSamples, ofs);
-      }
-    criterion.reset();
+    m_InLayers.push_back( LayerType(previousShape, m_NumberOfHiddenNeurons[i-1]) );
+    previousShape = m_NumberOfHiddenNeurons[i-1];
     }
-  else
-    {
-    shark::MaxIterations<> criterion(m_NumberOfIterations);
+  m_OutLayer = OutLayerType(previousShape, shark::dataDimension(inputSamples));
 
-    OutAutoencoderType net;
-    // Shark doesn't allow to train a layer using a sparsity term AND a noisy input.
-    if (m_Noise[0] != 0)
-      {
-      TrainOneLayer(criterion, net, 0, inputSamples, ofs);
-      otbMsgDevMacro(<< "m_Noise " << m_Noise[0]);
-      }
-    else
-      {
-      TrainOneSparseLayer(criterion, net, 0, inputSamples, ofs);
-      }
-    criterion.reset();
-    }
-
-  // Training of the other autoencoders
-  if (m_Epsilon > 0)
+  // Training of the autoencoders pairwise, starting from the first and last layers
+  for (unsigned int i = 0 ; i < m_NumberOfHiddenNeurons.Size(); ++i)
     {
-    shark::TrainingProgress<> criterion(5,m_Epsilon);
-
-    for (unsigned int i = 1 ; i < m_NumberOfHiddenNeurons.Size(); ++i)
+    if (m_Epsilon > 0)
       {
-      AutoencoderType net;
+      shark::TrainingProgress<> criterion(5,m_Epsilon);
       // Shark doesn't allow to train a layer using a sparsity term AND a noisy input.
       if (m_Noise[i] != 0)
         {
-        TrainOneLayer(criterion, net, i, inputSamples, ofs);
+        TrainOneLayer(criterion, i, inputSamples, ofs);
         }
       else
         {
-        TrainOneSparseLayer(criterion, net, i, inputSamples, ofs);
+        TrainOneSparseLayer(criterion, i, inputSamples, ofs);
         }
       criterion.reset();
       }
-    }
-  else
-    {
-    shark::MaxIterations<> criterion(m_NumberOfIterations);
-
-    for (unsigned int i = 1 ; i < m_NumberOfHiddenNeurons.Size(); ++i)
+    else
       {
-      AutoencoderType net;
+      shark::MaxIterations<> criterion(m_NumberOfIterations);
       // Shark doesn't allow to train a layer using a sparsity term AND a noisy input.
       if (m_Noise[i] != 0)
         {
-        TrainOneLayer(criterion, net, i, inputSamples, ofs);
+        TrainOneLayer(criterion, i, inputSamples, ofs);
         otbMsgDevMacro(<< "m_Noise " << m_Noise[0]);
         }
       else
         {
-        TrainOneSparseLayer( criterion, net, i, inputSamples, ofs);
+        TrainOneSparseLayer( criterion, i, inputSamples, ofs);
         }
       criterion.reset();
       }
+    // encode the samples with the last encoder trained
+    inputSamples = m_InLayers[i](inputSamples);
     }
   if (m_NumberOfIterationsFineTuning > 0)
     {
@@ -183,31 +142,37 @@ AutoencoderModel<TInputValue,NeuronType>
 }
 
 template <class TInputValue, class NeuronType>
-template <class T, class Autoencoder>
+template <class T>
 void
 AutoencoderModel<TInputValue,NeuronType>
 ::TrainOneLayer(
   shark::AbstractStoppingCriterion<T> & criterion,
-  Autoencoder & net,
   unsigned int layer_index,
   shark::Data<shark::RealVector> &samples,
   std::ostream& File)
 {
+  typedef shark::AbstractModel<shark::RealVector,shark::RealVector> BaseModelType;
+  ModelType net;
+  net.add(&(m_InLayers[layer_index]), true);
+  net.add( (layer_index ?
+    (BaseModelType*) &(m_InLayers[m_NumberOfHiddenNeurons.Size()*2 - 1 - layer_index]) :
+    (BaseModelType*) &m_OutLayer) , true);
+
   otbMsgDevMacro(<< "Noise " <<  m_Noise[layer_index]);
   std::size_t inputs = dataDimension(samples);
-  net.setStructure(inputs, m_NumberOfHiddenNeurons[layer_index]);
   initRandomUniform(net,-m_InitFactor*std::sqrt(1.0/inputs),m_InitFactor*std::sqrt(1.0/inputs));
 
-  shark::ImpulseNoiseModel noise(inputs,m_Noise[layer_index],1.0); //set an input pixel with probability m_Noise to 0
-  shark::ConcatenatedModel<shark::RealVector,shark::RealVector> model = noise>> net;
+  //~ shark::ImpulseNoiseModel noise(inputs,m_Noise[layer_index],1.0); //set an input pixel with probability m_Noise to 0
+  //~ shark::ConcatenatedModel<shark::RealVector,shark::RealVector> model = noise>> net;
   shark::LabeledData<shark::RealVector,shark::RealVector> trainSet(samples,samples);//labels identical to inputs
   shark::SquaredLoss<shark::RealVector> loss;
-  shark::ErrorFunction error(trainSet, &model, &loss);
+  //~ shark::ErrorFunction error(trainSet, &model, &loss);
+  shark::ErrorFunction<> error(trainSet, &net, &loss);
 
-  shark::TwoNormRegularizer regularizer(error.numberOfVariables());
+  shark::TwoNormRegularizer<> regularizer(error.numberOfVariables());
   error.setRegularizer(m_Regularization[layer_index],&regularizer);
 
-  shark::IRpropPlusFull optimizer;
+  shark::Adam<> optimizer;
   error.init();
   optimizer.init(error);
 
@@ -230,35 +195,37 @@ AutoencoderModel<TInputValue,NeuronType>
     } while( !criterion.stop( optimizer.solution() ) );
 
   net.setParameterVector(optimizer.solution().point);
-  m_Net.setLayer(layer_index,net.encoderMatrix(),net.hiddenBias());  // Copy the encoder in the FF neural network
-  m_Net.setLayer( m_NumberOfHiddenNeurons.Size()*2 - 1 - layer_index,net.decoderMatrix(),net.outputBias()); // Copy the decoder in the FF neural network
-  samples = net.encode(samples);
 }
 
 template <class TInputValue, class NeuronType>
-template <class T, class Autoencoder>
+template <class T>
 void AutoencoderModel<TInputValue,NeuronType>::TrainOneSparseLayer(
   shark::AbstractStoppingCriterion<T> & criterion,
-  Autoencoder & net,
   unsigned int layer_index,
   shark::Data<shark::RealVector> &samples,
   std::ostream& File)
 {
-  //AutoencoderType net;
-  std::size_t inputs = dataDimension(samples);
-  net.setStructure(inputs, m_NumberOfHiddenNeurons[layer_index]);
+  typedef shark::AbstractModel<shark::RealVector,shark::RealVector> BaseModelType;
+  ModelType net;
+  net.add(&(m_InLayers[layer_index]), true);
+  net.add( (layer_index ?
+    (BaseModelType*) &(m_InLayers[m_NumberOfHiddenNeurons.Size()*2 - 1 - layer_index]) :
+    (BaseModelType*) &m_OutLayer) , true);
 
+  std::size_t inputs = dataDimension(samples);
   shark::initRandomUniform(net,-m_InitFactor*std::sqrt(1.0/inputs),m_InitFactor*std::sqrt(1.0/inputs));
 
   // Idea : set the initials value for the output weights higher than the input weights
 
   shark::LabeledData<shark::RealVector,shark::RealVector> trainSet(samples,samples);//labels identical to inputs
   shark::SquaredLoss<shark::RealVector> loss;
-  shark::SparseAutoencoderError error(trainSet,&net, &loss, m_Rho[layer_index], m_Beta[layer_index]);
-
-  shark::TwoNormRegularizer regularizer(error.numberOfVariables());
+  //~ shark::SparseAutoencoderError error(trainSet,&net, &loss, m_Rho[layer_index], m_Beta[layer_index]);
+  // SparseAutoencoderError doesn't exist anymore, for now use a plain ErrorFunction
+  shark::ErrorFunction<> error(trainSet, &net, &loss);
+  
+  shark::TwoNormRegularizer<> regularizer(error.numberOfVariables());
   error.setRegularizer(m_Regularization[layer_index],&regularizer);
-  shark::IRpropPlusFull optimizer;
+  shark::Adam<> optimizer;
   error.init();
   optimizer.init(error);
 
@@ -279,9 +246,6 @@ void AutoencoderModel<TInputValue,NeuronType>::TrainOneSparseLayer(
     File << "end layer" << std::endl;
     }
   net.setParameterVector(optimizer.solution().point);
-  m_Net.setLayer(layer_index,net.encoderMatrix(),net.hiddenBias());  // Copy the encoder in the FF neural network
-  m_Net.setLayer( m_NumberOfHiddenNeurons.Size()*2 - 1 - layer_index,net.decoderMatrix(),net.outputBias()); // Copy the decoder in the FF neural network
-  samples = net.encode(samples);
 }
 
 template <class TInputValue, class NeuronType>
@@ -293,15 +257,23 @@ AutoencoderModel<TInputValue,NeuronType>
   shark::Data<shark::RealVector> &samples,
   std::ostream& File)
 {
+  // create full network
+  ModelType net;
+  for (auto &layer : m_InLayers)
+    {
+    net.add(&layer, true);
+    }
+  net.add(&m_OutLayer, true);
+  
   //labels identical to inputs
   shark::LabeledData<shark::RealVector,shark::RealVector> trainSet(samples,samples);
   shark::SquaredLoss<shark::RealVector> loss;
 
-  shark::ErrorFunction error(trainSet, &m_Net, &loss);
-  shark::TwoNormRegularizer regularizer(error.numberOfVariables());
+  shark::ErrorFunction<> error(trainSet, &net, &loss);
+  shark::TwoNormRegularizer<> regularizer(error.numberOfVariables());
   error.setRegularizer(m_Regularization[0],&regularizer);
 
-  shark::IRpropPlusFull optimizer;
+  shark::Adam<> optimizer;
   error.init();
   optimizer.init(error);
   otbMsgDevMacro(<<"Error before training : " << optimizer.solution().value);
@@ -326,7 +298,6 @@ AutoencoderModel<TInputValue,NeuronType>
   try
     {
     this->Load(filename);
-    m_Net.name();
     }
   catch(...)
     {
@@ -350,22 +321,15 @@ AutoencoderModel<TInputValue,NeuronType>
 {
   otbMsgDevMacro(<< "saving model ...");
   std::ofstream ofs(filename);
-  ofs << m_Net.name() << std::endl; // the first line of the model file contains a key
+  ofs << "Autoencoder" << std::endl; // the first line of the model file contains a key
+  ofs << (m_InLayers.size() + 1) << std::endl; // second line is the number of encoders/decoders 
   shark::TextOutArchive oa(ofs);
-  oa << m_Net;
-  ofs.close();
-
-  if (this->m_WriteWeights == true)     // output the map vectors in a txt file
+  for (const auto &layer : m_InLayers)
     {
-    std::ofstream otxt(filename+".txt");
-    for (unsigned int i = 0 ; i < m_Net.layerMatrices().size(); ++i)
-      {
-      otxt << "layer " << i << std::endl;
-      otxt << m_Net.layerMatrix(i) << std::endl;
-      otxt << m_Net.bias(i) << std::endl;
-      otxt << std::endl;
-      }
+    oa << layer;
     }
+  oa << m_OutLayer;
+  ofs.close();
 }
 
 template <class TInputValue, class NeuronType>
@@ -373,23 +337,39 @@ void
 AutoencoderModel<TInputValue,NeuronType>
 ::Load(const std::string & filename, const std::string & /*name*/)
 {
-  NetworkType net;
   std::ifstream ifs(filename);
-  char autoencoder[256];
-  ifs.getline(autoencoder,256);
-  std::string autoencoderstr(autoencoder);
-
-  if (autoencoderstr != net.name()){
+  char buffer[256];
+  // check first line
+  ifs.getline(buffer,256);
+  std::string bufferStr(buffer);
+  if (bufferStr != "Autoencoder"){
     itkExceptionMacro(<< "Error opening " << filename.c_str() );
     }
+  // check second line
+  ifs.getline(buffer,256);
+  int nbLevels = boost::lexical_cast<int>(buffer);
+  if (nbLevels < 2 || nbLevels%2 == 1)
+    {
+    itkExceptionMacro(<< "Unexpected number of levels : "<<buffer );
+    }
+  m_InLayers.clear();
+  m_Encoder = ModelType();
   shark::TextInArchive ia(ifs);
-  ia >> m_Net;
+  for (int i=0 ; (i+1) < nbLevels ; i++)
+    {
+    LayerType layer;
+    ia >> layer;
+    m_InLayers.push_back(layer);
+    }
+  ia >> m_OutLayer;
   ifs.close();
 
-  // This gives us the dimension if we keep the encoder and decoder
-  size_t feature_layer_index = m_Net.layerMatrices().size()/2;
-  // number of neurons in the feature layer (second dimension of the first decoder weight matrix)
-  this->SetDimension(m_Net.layerMatrix(feature_layer_index).size2());
+  for (int i=0 ; i < nbLevels/2 ; i++)
+    {
+    m_Encoder.add(&(m_InLayers[i]) ,true);
+    }
+
+  this->SetDimension( m_Encoder.outputShape()[0] );
 }
 
 template <class TInputValue, class NeuronType>
@@ -409,7 +389,7 @@ AutoencoderModel<TInputValue,NeuronType>
   shark::Data<shark::RealVector> data = shark::createDataFromRange(features);
 
   // features layer for a network containing the encoder and decoder part
-  data = m_Net.evalLayer( m_Net.layerMatrices().size()/2-1 ,data);
+  data = m_Encoder(data);
   TargetSampleType target;
   target.SetSize(this->m_Dimension);
 
@@ -435,7 +415,7 @@ AutoencoderModel<TInputValue,NeuronType>
   shark::Data<shark::RealVector> data = shark::createDataFromRange(features);
   TargetSampleType target;
   // features layer for a network containing the encoder and decoder part
-  data = m_Net.evalLayer( m_Net.layerMatrices().size()/2-1 ,data);
+  data = m_Encoder(data);
 
   unsigned int id = startIndex;
   target.SetSize(this->m_Dimension);
diff --git a/Modules/Learning/DimensionalityReductionLearning/include/otbPCAModel.txx b/Modules/Learning/DimensionalityReductionLearning/include/otbPCAModel.txx
index 9f39326a21bc5f1980a49d80ecdaea55b42a450a..a387852fecc386d9c5f2a6c27c7bf39cd7a3649d 100644
--- a/Modules/Learning/DimensionalityReductionLearning/include/otbPCAModel.txx
+++ b/Modules/Learning/DimensionalityReductionLearning/include/otbPCAModel.txx
@@ -137,11 +137,11 @@ PCAModel<TInputValue>::Load(const std::string & filename, const std::string & /*
   ifs.close();
   if (this->m_Dimension ==0)
   {
-    this->m_Dimension = m_Encoder.outputSize();
+    this->m_Dimension = m_Encoder.outputShape()[0];
   }
 
   auto eigenvectors = m_Encoder.matrix();
-  eigenvectors.resize(this->m_Dimension,m_Encoder.inputSize());
+  eigenvectors.resize(this->m_Dimension,m_Encoder.inputShape()[0]);
 
   m_Encoder.setStructure(eigenvectors, m_Encoder.offset() );
 }
diff --git a/Modules/Learning/LearningBase/otb-module.cmake b/Modules/Learning/LearningBase/otb-module.cmake
index afa2a339a1813cf16e5f6ea3700f079a36180dcd..c0af985032de6d4a2acd11988be1b9a177cf8219 100644
--- a/Modules/Learning/LearningBase/otb-module.cmake
+++ b/Modules/Learning/LearningBase/otb-module.cmake
@@ -28,7 +28,11 @@ otb_module(OTBLearningBase
     OTBImageBase
     OTBITK
 
-  TEST_DEPENDS
+    OPTIONAL_DEPENDS
+    OTBShark
+
+    TEST_DEPENDS
+    OTBBoost
     OTBTestKernel
     OTBImageIO
 
diff --git a/Modules/Learning/LearningBase/test/CMakeLists.txt b/Modules/Learning/LearningBase/test/CMakeLists.txt
index d1d16c3e65801e606c6e6903538b65264a4483a6..48e28cc5cad320ffa41eee0659ff6979d0bf4457 100644
--- a/Modules/Learning/LearningBase/test/CMakeLists.txt
+++ b/Modules/Learning/LearningBase/test/CMakeLists.txt
@@ -32,6 +32,10 @@ otbKMeansImageClassificationFilterNew.cxx
 otbMachineLearningModelTemplates.cxx
 )
 
+if(OTB_USE_SHARK)
+  set(OTBLearningBaseTests ${OTBLearningBaseTests} otbSharkUtilsTests.cxx)
+endif()
+
 add_executable(otbLearningBaseTestDriver ${OTBLearningBaseTests})
 target_link_libraries(otbLearningBaseTestDriver ${OTBLearningBase-Test_LIBRARIES})
 otb_module_target_label(otbLearningBaseTestDriver)
@@ -68,3 +72,7 @@ otb_add_test(NAME leTuDecisionTreeNew COMMAND otbLearningBaseTestDriver
 otb_add_test(NAME leTuKMeansImageClassificationFilterNew COMMAND otbLearningBaseTestDriver
   otbKMeansImageClassificationFilterNew)
 
+if(OTB_USE_SHARK)
+  otb_add_test(NAME leTuSharkNormalizeLabels COMMAND otbLearningBaseTestDriver
+    otbSharkNormalizeLabels)
+endif()
diff --git a/Modules/Learning/LearningBase/test/otbLearningBaseTestDriver.cxx b/Modules/Learning/LearningBase/test/otbLearningBaseTestDriver.cxx
index 5b38bf300dd4520c18e198b6e6643848cbdc937c..dc2d36b7943129ec6519ebbc4f194d1dd6078800 100644
--- a/Modules/Learning/LearningBase/test/otbLearningBaseTestDriver.cxx
+++ b/Modules/Learning/LearningBase/test/otbLearningBaseTestDriver.cxx
@@ -29,4 +29,7 @@ void RegisterTests()
   REGISTER_TEST(otbSEMClassifierNew);
   REGISTER_TEST(otbDecisionTreeNew);
   REGISTER_TEST(otbKMeansImageClassificationFilterNew);
+#ifdef OTB_USE_SHARK
+  REGISTER_TEST(otbSharkNormalizeLabels);
+#endif
 }
diff --git a/Modules/Learning/LearningBase/test/otbSharkUtilsTests.cxx b/Modules/Learning/LearningBase/test/otbSharkUtilsTests.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..bc3783cb728b0f5ad0f6b2d43620b18ba7939e30
--- /dev/null
+++ b/Modules/Learning/LearningBase/test/otbSharkUtilsTests.cxx
@@ -0,0 +1,54 @@
+/*
+ * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "itkMacro.h"
+#include "otbSharkUtils.h"
+
+
+int otbSharkNormalizeLabels(int itkNotUsed(argc), char* itkNotUsed(argv) [])
+{
+  std::vector<unsigned int> inLabels = {2, 2, 3, 20, 1};
+  std::vector<unsigned int> expectedDictionary = {2, 3, 20, 1};
+  std::vector<unsigned int> expectedLabels = {0, 0, 1, 2, 3};
+
+  auto newLabels = inLabels;
+  std::vector<unsigned int> labelDict;
+  otb::Shark::NormalizeLabelsAndGetDictionary(newLabels, labelDict);
+
+  if(newLabels != expectedLabels)
+    {
+    std::cout << "Wrong new labels\n";
+    for(size_t i = 0; i<newLabels.size(); ++i)
+      std::cout << "Got " << newLabels[i] << " expected " << expectedLabels[i] << '\n';
+
+    return EXIT_FAILURE;
+    }
+
+  if(labelDict != expectedDictionary)
+    {
+    std::cout << "Wrong dictionary\n";
+    for(size_t i = 0; i<labelDict.size(); ++i)
+      std::cout << "Got " << labelDict[i] << " expected " << expectedDictionary[i] << '\n';
+
+    return EXIT_FAILURE;
+    }
+
+  return EXIT_SUCCESS;
+}
diff --git a/Modules/Learning/Sampling/include/otbSampleAugmentation.h b/Modules/Learning/Sampling/include/otbSampleAugmentation.h
index 84e67ffdd413940e2646ff2f35d6423241ab5f3f..8196b3e9da759e330f6a652a13377b9f609e0d57 100644
--- a/Modules/Learning/Sampling/include/otbSampleAugmentation.h
+++ b/Modules/Learning/Sampling/include/otbSampleAugmentation.h
@@ -21,6 +21,10 @@
 #ifndef otbSampleAugmentation_h
 #define otbSampleAugmentation_h
 
+#ifdef _OPENMP
+ # include <omp.h>
+#endif
+
 #include <vector>
 #include <algorithm>
 #include <random>
@@ -49,7 +53,9 @@ SampleType EstimateStds(const SampleVectorType& samples)
   for(size_t i=0; i<nbSamples; ++i)
     {
     auto norm_factor = 1.0/(i+1);
-#pragma omp parallel for 
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif 
     for(size_t j=0; j<nbComponents; ++j)
       {
       const auto mu = means[j];
@@ -59,7 +65,9 @@ SampleType EstimateStds(const SampleVectorType& samples)
       means[j] = muNew;
       }
     }
+#ifdef _OPENMP
 #pragma omp parallel for
+#endif
   for(size_t j=0; j<nbComponents; ++j)
     {
     stds[j] = std::sqrt(stds[j]/nbSamples);
@@ -77,7 +85,9 @@ void ReplicateSamples(const SampleVectorType& inSamples,
 {
   newSamples.resize(nbSamples);
   size_t imod{0};
+#ifdef _OPENMP
 #pragma omp parallel for
+#endif
   for(size_t i=0; i<nbSamples; ++i)
     {
     if (imod == inSamples.size()) imod = 0;
@@ -108,14 +118,18 @@ void JitterSamples(const SampleVectorType& inSamples,
   // have different stds
   auto stds = EstimateStds(inSamples);
   std::vector<std::normal_distribution<double>> gaussDis(nbComponents);
+#ifdef _OPENMP
 #pragma omp parallel for
+#endif
   for(size_t i=0; i<nbComponents; ++i)
     gaussDis[i] = std::normal_distribution<double>{0.0, stds[i]/stdFactor};
 
   for(size_t i=0; i<nbSamples; ++i)
     {
     newSamples[i] = inSamples[std::rand()%inSamples.size()];
+#ifdef _OPENMP
 #pragma omp parallel for
+#endif
     for(size_t j=0; j<nbComponents; ++j)
       newSamples[i][j] += gaussDis[j](gen);
     }
@@ -157,7 +171,9 @@ void FindKNNIndices(const SampleVectorType& inSamples,
 {
   const auto nbSamples = inSamples.size();
   nnVector.resize(nbSamples);
+  #ifdef _OPENMP
   #pragma omp parallel for
+  #endif
   for(size_t sampleIdx=0; sampleIdx<nbSamples; ++sampleIdx)
     {
     NNIndicesType nns;
@@ -200,7 +216,9 @@ void Smote(const SampleVectorType& inSamples,
   FindKNNIndices(inSamples, nbNeighbors, nnVector);
   // The input samples are selected randomly with replacement
   std::srand(seed);
+  #ifdef _OPENMP
   #pragma omp parallel for
+  #endif
   for(size_t i=0; i<nbSamples; ++i)
     {
     const auto sampleIdx = std::rand()%(inSamples.size());
diff --git a/Modules/Learning/Supervised/include/otbSharkRandomForestsMachineLearningModel.h b/Modules/Learning/Supervised/include/otbSharkRandomForestsMachineLearningModel.h
index 7dd41d6eed22208d512ecf0e5d4eb5116f2bf5c0..9e4a985901c2f7007529c279ce156e371ebbcbd2 100644
--- a/Modules/Learning/Supervised/include/otbSharkRandomForestsMachineLearningModel.h
+++ b/Modules/Learning/Supervised/include/otbSharkRandomForestsMachineLearningModel.h
@@ -34,6 +34,7 @@
 #pragma GCC diagnostic ignored "-Wcast-align"
 #pragma GCC diagnostic ignored "-Wunknown-pragmas"
 #endif
+#include <shark/Models/Classifier.h>
 #include "otb_shark.h"
 #include "shark/Algorithms/Trainers/RFTrainer.h"
 #if defined(__GNUC__) || defined(__clang__)
@@ -134,6 +135,10 @@ public:
   /** If true, margin confidence value will be computed */
   itkSetMacro(ComputeMargin, bool);
 
+  /** If true, class labels will be normalised in [0 ... nbClasses] */
+  itkGetMacro(NormalizeClassLabels, bool);
+  itkSetMacro(NormalizeClassLabels, bool);
+
 protected:
   /** Constructor */
   SharkRandomForestsMachineLearningModel();
@@ -154,8 +159,10 @@ private:
   SharkRandomForestsMachineLearningModel(const Self &); //purposely not implemented
   void operator =(const Self&); //purposely not implemented
 
-  shark::RFClassifier m_RFModel;
-  shark::RFTrainer m_RFTrainer;
+  shark::RFClassifier<unsigned int> m_RFModel;
+  shark::RFTrainer<unsigned int> m_RFTrainer;
+  std::vector<unsigned int> m_ClassDictionary;
+  bool m_NormalizeClassLabels;
 
   unsigned int m_NumberOfTrees;
   unsigned int m_MTry;
diff --git a/Modules/Learning/Supervised/include/otbSharkRandomForestsMachineLearningModel.txx b/Modules/Learning/Supervised/include/otbSharkRandomForestsMachineLearningModel.txx
index 207f1abdd77e4b5cfffd9bc5d104c4b40232f853..72c816069bebddc048a0f8af48f24579a55fa38b 100644
--- a/Modules/Learning/Supervised/include/otbSharkRandomForestsMachineLearningModel.txx
+++ b/Modules/Learning/Supervised/include/otbSharkRandomForestsMachineLearningModel.txx
@@ -32,7 +32,6 @@
 #pragma GCC diagnostic ignored "-Woverloaded-virtual"
 #pragma GCC diagnostic ignored "-Wignored-qualifiers"
 #endif
-#include <shark/Models/Converter.h>
 #if defined(__GNUC__) || defined(__clang__)
 #pragma GCC diagnostic pop
 #endif
@@ -52,6 +51,7 @@ SharkRandomForestsMachineLearningModel<TInputValue,TOutputValue>
   this->m_ConfidenceIndex = true;
   this->m_IsRegressionSupported = false;
   this->m_IsDoPredictBatchMultiThreaded = true;
+  this->m_NormalizeClassLabels = true;
 }
 
 
@@ -76,13 +76,17 @@ SharkRandomForestsMachineLearningModel<TInputValue,TOutputValue>
 
   Shark::ListSampleToSharkVector(this->GetInputListSample(), features);
   Shark::ListSampleToSharkVector(this->GetTargetListSample(), class_labels);
+  if(m_NormalizeClassLabels)
+    {
+    Shark::NormalizeLabelsAndGetDictionary(class_labels, m_ClassDictionary);
+    }
   shark::ClassificationDataset TrainSamples = shark::createLabeledDataFromRange(features,class_labels);
 
   //Set parameters
   m_RFTrainer.setMTry(m_MTry);
   m_RFTrainer.setNTrees(m_NumberOfTrees);
   m_RFTrainer.setNodeSize(m_NodeSize);
-  m_RFTrainer.setOOBratio(m_OobRatio);
+  //  m_RFTrainer.setOOBratio(m_OobRatio);
   m_RFTrainer.train(m_RFModel, TrainSamples);
 
 }
@@ -125,15 +129,20 @@ SharkRandomForestsMachineLearningModel<TInputValue,TOutputValue>
     }
   if (quality != ITK_NULLPTR)
     {
-    shark::RealVector probas = m_RFModel(samples);
+    shark::RealVector probas = m_RFModel.decisionFunction()(samples);
     (*quality) = ComputeConfidence(probas, m_ComputeMargin);
     }
-  shark::ArgMaxConverter<shark::RFClassifier> amc;
-  amc.decisionFunction() = m_RFModel;
-  unsigned int res;
-  amc.eval(samples, res);
+  unsigned int res{0};
+  m_RFModel.eval(samples, res);
   TargetSampleType target;
-  target[0] = static_cast<TOutputValue>(res);
+  if(m_NormalizeClassLabels)
+    {
+    target[0] = m_ClassDictionary[static_cast<TOutputValue>(res)];
+    }
+  else
+    {
+    target[0] = static_cast<TOutputValue>(res);
+    }
   return target;
 }
 
@@ -157,13 +166,13 @@ SharkRandomForestsMachineLearningModel<TInputValue,TOutputValue>
   Shark::ListSampleRangeToSharkVector(input, features,startIndex,size);
   shark::Data<shark::RealVector> inputSamples = shark::createDataFromRange(features);
 
-  #ifdef _OPENMP
+#ifdef _OPENMP
   omp_set_num_threads(itk::MultiThreader::GetGlobalDefaultNumberOfThreads());
-  #endif
+#endif
   
   if(quality != ITK_NULLPTR)
     {
-    shark::Data<shark::RealVector> probas = m_RFModel(inputSamples);
+    shark::Data<shark::RealVector> probas = m_RFModel.decisionFunction()(inputSamples);
     unsigned int id = startIndex;
     for(shark::RealVector && p : probas.elements())
       {
@@ -175,14 +184,19 @@ SharkRandomForestsMachineLearningModel<TInputValue,TOutputValue>
       }
     }
     
-  shark::ArgMaxConverter<shark::RFClassifier> amc;
-  amc.decisionFunction() = m_RFModel;
-  auto prediction = amc(inputSamples);
+  auto prediction = m_RFModel(inputSamples);
   unsigned int id = startIndex;
   for(const auto& p : prediction.elements())
     {
     TargetSampleType target;
-    target[0] = static_cast<TOutputValue>(p);
+    if(m_NormalizeClassLabels)
+      {
+      target[0] = m_ClassDictionary[static_cast<TOutputValue>(p)];
+      }
+    else
+      {
+      target[0] = static_cast<TOutputValue>(p);
+      }
     targets->SetMeasurementVector(id,target);
     ++id;
     }
@@ -199,7 +213,18 @@ SharkRandomForestsMachineLearningModel<TInputValue,TOutputValue>
     itkExceptionMacro(<< "Error opening " << filename.c_str() );
     }
   // Add comment with model file name
-  ofs << "#" << m_RFModel.name() << std::endl;
+  ofs << "#" << m_RFModel.name();
+  if(m_NormalizeClassLabels) ofs  << " with_dictionary";
+  ofs << std::endl;
+  if(m_NormalizeClassLabels)
+    {
+    ofs << m_ClassDictionary.size() << " ";
+    for(const auto& l : m_ClassDictionary)
+      {
+      ofs << l << " ";
+      }
+    ofs << std::endl;
+    }
   shark::TextOutArchive oa(ofs);
   m_RFModel.save(oa,0);
 }
@@ -219,6 +244,10 @@ SharkRandomForestsMachineLearningModel<TInputValue,TOutputValue>
       {
       if( line.find( m_RFModel.name() ) == std::string::npos )
         itkExceptionMacro( "The model file : " + filename + " cannot be read." );
+      if( line.find( "with_dictionary" ) == std::string::npos )
+        {
+        m_NormalizeClassLabels=false;
+        }
       }
     else
       {
@@ -226,6 +255,18 @@ SharkRandomForestsMachineLearningModel<TInputValue,TOutputValue>
       ifs.clear();
       ifs.seekg( 0, std::ios::beg );
       }
+    if(m_NormalizeClassLabels)
+      {
+      size_t nbLabels{0};
+      ifs >> nbLabels;
+      m_ClassDictionary.resize(nbLabels);
+      for(size_t i=0; i<nbLabels; ++i)
+        {
+        unsigned int label;
+        ifs >> label;
+        m_ClassDictionary[i]=label;
+        }
+      }
     shark::TextInArchive ia( ifs );
     m_RFModel.load( ia, 0 );
     }
diff --git a/Modules/Learning/Unsupervised/include/otbSharkKMeansMachineLearningModel.txx b/Modules/Learning/Unsupervised/include/otbSharkKMeansMachineLearningModel.txx
index 9dd43948a719c9305dace0a6366ebfd40e4b3e24..1b08d538c943001279d9401f314d51e21e8dbf88 100644
--- a/Modules/Learning/Unsupervised/include/otbSharkKMeansMachineLearningModel.txx
+++ b/Modules/Learning/Unsupervised/include/otbSharkKMeansMachineLearningModel.txx
@@ -55,6 +55,7 @@ SharkKMeansMachineLearningModel<TInputValue, TOutputValue>
         m_Normalized( false ), m_K(2), m_MaximumNumberOfIterations( 10 )
 {
   // Default set HardClusteringModel
+  this->m_ConfidenceIndex = true;
   m_ClusteringModel = boost::make_shared<ClusteringModelType>( &m_Centroids );
 }
 
@@ -174,7 +175,7 @@ SharkKMeansMachineLearningModel<TInputValue, TOutputValue>
   // Change quality measurement only if SoftClustering or other clustering method is used.
   if( quality != ITK_NULLPTR )
     {
-    for( unsigned int qid = startIndex; qid < size; ++qid )
+    for( unsigned int qid = startIndex; qid < startIndex+size; ++qid )
       {
       quality->SetMeasurementVector( qid, static_cast<ConfidenceValueType>(1.) );
       }
diff --git a/Modules/MPI/MPIVrtWriter/include/otbMPIVrtWriter.h b/Modules/MPI/MPIVrtWriter/include/otbMPIVrtWriter.h
index ecec379e6d0cdaab2a6b39817b61ccd1a533c524..1800dc62ae4a170c70f0039906c79f2b2093e40d 100644
--- a/Modules/MPI/MPIVrtWriter/include/otbMPIVrtWriter.h
+++ b/Modules/MPI/MPIVrtWriter/include/otbMPIVrtWriter.h
@@ -47,232 +47,102 @@
 namespace otb {
 
 /**
- *\brief Write image data to multiple files with MPI processus and add a VRT file.
+ * \class MPIVrtWriter
  *
- * The image is divided into several pieces.
- * Each pieces is distributed to a MPI processus.
- * Each MPI processus write their pieces into a separate
- * file.
- * The master processus writes a VRT file (optional).
+ * \brief Write each tile of an image into a separate Tiff file and join them in a VRT
  *
- *\param img Image
- *\param output Output Filename
- *\param availableRAM Available memory for streaming
- *\param writeVRTFile Activate the VRT file writing
+ * \ingroup OTBMPIVrtWriter
  */
-template <typename TImage> void WriteMPI(TImage *img, const std::string &output, unsigned int availableRAM = 0, bool writeVRTFile=true)
+template <typename TImage>
+class MPIVrtWriter: public itk::ProcessObject
 {
-  typename otb::MPIConfig::Pointer mpiConfig = otb::MPIConfig::Instance();
-
-  unsigned int myRank = mpiConfig->GetMyRank();
-  unsigned int nbProcs = mpiConfig->GetNbProcs();
-
-  typedef otb::ImageFileWriter<TImage>                                           WriterType;
-  typedef otb::NumberOfDivisionsTiledStreamingManager<TImage>                    StreamingManagerType;
-  typedef itk::RegionOfInterestImageFilter<TImage, TImage>                       ExtractFilterType;
-
-  // First, update infomration from image to write
-  img->UpdateOutputInformation();
-
-  // Configure streaming manager
-  typename StreamingManagerType::Pointer streamingManager = StreamingManagerType::New();
-  streamingManager->SetNumberOfDivisions(nbProcs);
-  streamingManager->PrepareStreaming(img,img->GetLargestPossibleRegion());
-  unsigned int numberOfSplits = streamingManager->GetNumberOfSplits();
-
-  // This vector will hold all regions to write for current rank
-  std::vector<typename TImage::RegionType> regionsToWrite;
-
-  // Handle both cases when there are much more (resp. less) region to
-  // write than NbProcs
-  if(myRank < numberOfSplits)
-  {
-    unsigned int splitIdx = myRank;
-    while(splitIdx < numberOfSplits)
-    {
-      typename TImage::RegionType currentRegion = streamingManager->GetSplit(splitIdx);
-      regionsToWrite.push_back(currentRegion);
-      splitIdx+=nbProcs;
-    }
-  }
-
-  // Output prefix
-  std::string extension = itksys::SystemTools::GetFilenameExtension(output);
-  if (extension != ".vrt")
-  {
-
-  // TODO: Maybe remove this
-    if (extension == "")
-	{
-	  // Missing extension
-	  mpiConfig->logInfo("Filename has no extension. Adding <.vrt> extension.");
-	}
-	else
-	{
-	  // Bad extension
-	  mpiConfig->logError("Filename must have .vrt extension!");
-	  mpiConfig->abort(EXIT_FAILURE);
-	}
-  }
-  std::vector<std::string> joins;
-  joins.push_back(itksys::SystemTools::GetFilenamePath(output).append("/"));
-  joins.push_back(itksys::SystemTools::GetFilenameWithoutExtension(output));
-  std::string prefix = itksys::SystemTools::JoinPath(joins);
-
+public:
+  /** Standard class typedefs. */
+  typedef MPIVrtWriter                                      Self;
+  typedef itk::ProcessObject                                Superclass;
+  typedef itk::SmartPointer<Self>                           Pointer;
+  typedef itk::SmartPointer<const Self>                     ConstPointer;
 
-  // Data type
-  std::string dataTypeStr = "Float32";
-  GDALImageIO::Pointer gdalImageIO;
+  typedef TImage InputImageType;
 
-  // Now write all the regions
-  for(typename std::vector<typename TImage::RegionType>::const_iterator it = regionsToWrite.begin();
-      it!=regionsToWrite.end();++it)
-  {
-    typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
-    extractFilter->SetInput(img);
-    extractFilter->SetRegionOfInterest(*it);
-    // Writer
-	  // Output Filename
-    std::stringstream ss;
-    ss<<prefix<<"_"<<it->GetIndex()[0]<<"_"<<it->GetIndex()[1]<<"_"<<it->GetSize()[0]<<"_"<<it->GetSize()[1]<<".tif";
-    typename WriterType::Pointer writer = WriterType::New();
-    writer->SetFileName(ss.str());
-    writer->SetInput(extractFilter->GetOutput());
-    // Datatype
-    gdalImageIO = dynamic_cast<GDALImageIO *>(writer->GetImageIO());
-    if(gdalImageIO.IsNotNull())
-    {
-      dataTypeStr = gdalImageIO->GetGdalPixelTypeAsString();
-    }
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
 
-    if (!availableRAM)
-      {
-      writer->SetNumberOfDivisionsTiledStreaming(0);
-      }
-    else
-      {
-      writer->SetAutomaticAdaptativeStreaming(availableRAM);
-      }
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(MPIVrtWriter, itk::ProcessObject);
 
-    // Pipeline execution
-    try
-    {
-      writer->Update();
-    }
-    catch (itk::ExceptionObject& err)
-    {
-      std::stringstream message;
-      message << "ExceptionObject caught: " << err << std::endl;
-      mpiConfig->logError(message.str());
-      mpiConfig->abort(EXIT_FAILURE);
-    }
-  }
+  using Superclass::SetInput;
+  virtual void SetInput(const InputImageType *input);
 
-  // MPI process synchronization
-  mpiConfig->barrier();
+  /** Get writer only input */
+  const InputImageType* GetInput();
 
-  // Write VRT file
-  try
-  {
-    if (writeVRTFile && (myRank == 0))
-    {
-      // VRT Filename
-      std::stringstream vrtfOut;
-      vrtfOut<< prefix << ".vrt";
+  /** Does the real work. */
+  virtual void Update() override;
 
-      // Data type
-      GDALDataType dataType;
-      dataType = GDALGetDataTypeByName(dataTypeStr.c_str());
+  /** SimpleParallelTiffWriter Methods */
+  virtual void SetFileName(const char* extendedFileName);
+  virtual void SetFileName(std::string extendedFileName);
+  virtual const char* GetFileName () const;
 
-      int imageSizeY = img->GetLargestPossibleRegion().GetSize()[1];
-      int imageSizeX = img->GetLargestPossibleRegion().GetSize()[0];
-      const unsigned int nbBands = img->GetNumberOfComponentsPerPixel();
+  /** Specify the region to write. If left NULL, then the whole image
+   * is written. */
+  void SetIORegion(const itk::ImageIORegion& region);
+  itkGetConstReferenceMacro(IORegion, itk::ImageIORegion);
 
-      // Get VRT driver
-      GDALAllRegister();
-      GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("VRT");
-      if (driver == NULL) {
-         mpiConfig->logError("Error opening VRT driver.");
-         mpiConfig->abort(EXIT_FAILURE);
-      }
+  itkSetMacro(WriteVRT, bool);
+  itkGetMacro(WriteVRT, bool);
 
-      // Create output raster
-      GDALDataset *VRTOutput = driver->Create(vrtfOut.str().c_str(),
-            imageSizeX,
-            imageSizeY,
-            0,
-            dataType,
-            NULL); // No options
-      if (VRTOutput == NULL) {
-        mpiConfig->logError("driver->Create call failed.\n");
-        mpiConfig->abort(EXIT_FAILURE);
-      }
+  itkSetMacro(AvailableRAM, unsigned int);
+  itkGetMacro(AvailableRAM, unsigned int);
 
-      // Set GeoTransform
-      double gt[6];
-      gt[0] = img->GetOrigin()[0] - 0.5*img->GetSignedSpacing()[0];
-      gt[1] = img->GetSignedSpacing()[0];
-      gt[2] = 0.0;
-      gt[3] = img->GetOrigin()[1] - 0.5*img->GetSignedSpacing()[1];
-      gt[4] = 0.0;
-      gt[5] = img->GetSignedSpacing()[1];
-      VRTOutput->SetGeoTransform(gt);
+protected:
+  MPIVrtWriter();
+  virtual ~MPIVrtWriter();
+  void PrintSelf(std::ostream& os, itk::Indent indent) const override;
 
-      // Set projection
-      OGRSpatialReference out_sr;
-      char *wkt = NULL;
-      out_sr.SetFromUserInput(img->GetProjectionRef().c_str());
-      out_sr.exportToWkt(&wkt);
-      VRTOutput->SetProjection(wkt);
+private:
+  MPIVrtWriter(const MPIVrtWriter &) = delete;
+  void operator =(const MPIVrtWriter&) = delete;
 
-      for(unsigned int band = 1; band<=nbBands;++band)
-      {
-        VRTOutput->AddBand(dataType, NULL);
+  unsigned int m_AvailableRAM;
 
-        typename TImage::RegionType currentRegion;
-        for(unsigned int id=0; id < numberOfSplits; ++id)
-        {
-          currentRegion = streamingManager->GetSplit(id);
-          int tileSizeX = currentRegion.GetSize()[0];
-          int tileSizeY = currentRegion.GetSize()[1];
-          int tileIndexX = currentRegion.GetIndex()[0];
-          int tileIndexY = currentRegion.GetIndex()[1];
-          std::stringstream tileFileName;
-          tileFileName <<prefix<<"_"<<tileIndexX<<"_"<<tileIndexY<<"_"<<tileSizeX<<"_"<<tileSizeY<<".tif";
-          std::cout<<tileFileName.str()<<std::endl;
+  itk::ImageIORegion m_IORegion;
 
-          GDALDataset *dataset = (GDALDataset*) GDALOpen(tileFileName.str().c_str(), GA_ReadOnly);
+  std::string m_Filename;
 
-          VRTSourcedRasterBand *VRTBand = dynamic_cast<VRTSourcedRasterBand*> (VRTOutput->GetRasterBand(band));
-          VRTBand->AddSimpleSource(dataset->GetRasterBand(band),
-                                      0, //xOffSrc
-                                      0, //yOffSrc
-                                      tileSizeX, //xSizeSrc
-                                      tileSizeY, //ySizeSrc
-                                      tileIndexX, //xOffDest
-                                      tileIndexY, //yOffDest
-                                      tileSizeX, //xSizeDest
-                                      tileSizeY, //ySizeDest
-                                      "near",
-                                      VRT_NODATA_UNSET);
-        }
+  bool m_WriteVRT;
 
-      }
-
-      // Close
-      CPLFree(wkt);
-      GDALClose(VRTOutput);
-    }
-  }
-  catch (itk::ExceptionObject& err)
-  {
-    std::stringstream message;
-    message << "ExceptionObject caught: " << err << std::endl;
-    mpiConfig->logError(message.str());
-    mpiConfig->abort(EXIT_FAILURE);
-  }
+};
 
+/**
+ *\brief Write image data to multiple files with MPI processus and add a VRT file.
+ *
+ * The image is divided into several pieces.
+ * Each pieces is distributed to a MPI processus.
+ * Each MPI processus write their pieces into a separate
+ * file.
+ * The master processus writes a VRT file (optional).
+ *
+ *\param img Image
+ *\param output Output Filename
+ *\param availableRAM Available memory for streaming
+ *\param writeVRTFile Activate the VRT file writing
+ */
+template <typename TImage> void WriteMPI(TImage *img, const std::string &output, unsigned int availableRAM = 0, bool writeVRTFile=true)
+{
+  typename MPIVrtWriter<TImage>::Pointer writer = MPIVrtWriter<TImage>::New();
+  writer->SetInput(img);
+  writer->SetFileName(output);
+  writer->SetAvailableRAM(availableRAM);
+  writer->SetWriteVRT(writeVRTFile);
+  writer->Update();
 }
 
 } // End namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbMPIVrtWriter.txx"
+#endif
+
 #endif //__otbMPIVrtWriter_h
diff --git a/Modules/MPI/MPIVrtWriter/include/otbMPIVrtWriter.txx b/Modules/MPI/MPIVrtWriter/include/otbMPIVrtWriter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..b57518ed73d613b37e33b86c59ba732ed0e44cc2
--- /dev/null
+++ b/Modules/MPI/MPIVrtWriter/include/otbMPIVrtWriter.txx
@@ -0,0 +1,327 @@
+/*
+ * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbMPIVrtWriter_txx
+#define otbMPIVrtWriter_txx
+
+#include "otbMPIVrtWriter.h"
+#include "otbMacro.h"
+
+namespace otb
+{
+
+template <typename TImage>
+MPIVrtWriter<TImage>::MPIVrtWriter()
+  : m_AvailableRAM(0)
+  , m_IORegion()
+  , m_Filename("")
+  , m_WriteVRT(true)
+{
+}
+
+template <typename TImage>
+MPIVrtWriter<TImage>::~MPIVrtWriter()
+{
+}
+
+template <typename TImage>
+void
+MPIVrtWriter<TImage>::SetInput(const InputImageType *input)
+{
+  this->ProcessObject::SetNthInput(0,const_cast<InputImageType*>(input));
+}
+
+template <typename TImage>
+const TImage*
+MPIVrtWriter<TImage>::GetInput()
+{
+  if (this->GetNumberOfInputs() < 1)
+    {
+    return 0;
+    }
+  return static_cast<const InputImageType*>(this->ProcessObject::GetInput(0));
+}
+
+template <typename TImage>
+void
+MPIVrtWriter<TImage>::SetFileName(const char* extendedFileName)
+{
+  if (m_Filename.compare(extendedFileName) != 0 )
+    {
+    m_Filename = std::string(extendedFileName);
+    this->Modified();
+    }
+}
+
+template <typename TImage>
+void
+MPIVrtWriter<TImage>::SetFileName(std::string extendedFileName)
+{
+  this->SetFileName(extendedFileName.c_str());
+}
+
+template <typename TImage>
+const char*
+MPIVrtWriter<TImage>::GetFileName () const
+{
+  return m_Filename.c_str();
+}
+
+template <typename TImage>
+void
+MPIVrtWriter<TImage>::SetIORegion(const itk::ImageIORegion& region)
+{
+  if (m_IORegion != region)
+    {
+    m_IORegion = region;
+    this->Modified();
+    }
+}
+
+template <typename TImage>
+void
+MPIVrtWriter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+  os << indent << "File Name: "<< m_Filename << std::endl;
+  os << indent << "Available RAM: "<< m_AvailableRAM << std::endl;
+  os << indent << "Write VRT: "<< m_WriteVRT << std::endl;
+}
+
+template <typename TImage>
+void
+MPIVrtWriter<TImage>::Update()
+{
+  typename otb::MPIConfig::Pointer mpiConfig = otb::MPIConfig::Instance();
+
+  unsigned int myRank = mpiConfig->GetMyRank();
+  unsigned int nbProcs = mpiConfig->GetNbProcs();
+
+  typedef otb::ImageFileWriter<TImage>                                           WriterType;
+  typedef otb::NumberOfDivisionsTiledStreamingManager<TImage>                    StreamingManagerType;
+  typedef itk::RegionOfInterestImageFilter<TImage, TImage>                       ExtractFilterType;
+
+  // First, update infomration from image to write
+  UpdateOutputInformation();
+  InputImageType* img = const_cast<InputImageType*>(GetInput());
+  std::string output = GetFileName();
+
+  // Configure streaming manager
+  typename StreamingManagerType::Pointer streamingManager = StreamingManagerType::New();
+  streamingManager->SetNumberOfDivisions(nbProcs);
+  streamingManager->PrepareStreaming(img,img->GetLargestPossibleRegion());
+  unsigned int numberOfSplits = streamingManager->GetNumberOfSplits();
+
+  // This vector will hold all regions to write for current rank
+  std::vector<typename TImage::RegionType> regionsToWrite;
+
+  // Handle both cases when there are much more (resp. less) region to
+  // write than NbProcs
+  if(myRank < numberOfSplits)
+  {
+    unsigned int splitIdx = myRank;
+    while(splitIdx < numberOfSplits)
+    {
+      typename TImage::RegionType currentRegion = streamingManager->GetSplit(splitIdx);
+      regionsToWrite.push_back(currentRegion);
+      splitIdx+=nbProcs;
+    }
+  }
+
+  // Output prefix
+  std::string extension = itksys::SystemTools::GetFilenameExtension(output);
+  if (extension != ".vrt")
+  {
+
+  // TODO: Maybe remove this
+    if (extension == "")
+	{
+	  // Missing extension
+	  mpiConfig->logInfo("Filename has no extension. Adding <.vrt> extension.");
+	}
+	else
+	{
+	  // Bad extension
+	  mpiConfig->logError("Filename must have .vrt extension!");
+	  mpiConfig->abort(EXIT_FAILURE);
+	}
+  }
+  std::vector<std::string> joins;
+  itksys::SystemTools::SplitPath(output, joins);
+  joins.back() = itksys::SystemTools::GetFilenameWithoutExtension(output);
+  std::string prefix = itksys::SystemTools::JoinPath(joins);
+
+  // Data type
+  std::string dataTypeStr = "Float32";
+  GDALImageIO::Pointer gdalImageIO;
+
+  // Now write all the regions
+  for(typename std::vector<typename TImage::RegionType>::const_iterator it = regionsToWrite.begin();
+      it!=regionsToWrite.end();++it)
+  {
+    typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
+    extractFilter->SetInput(img);
+    extractFilter->SetRegionOfInterest(*it);
+    // Writer
+	  // Output Filename
+    std::stringstream ss;
+    ss<<prefix<<"_"<<it->GetIndex()[0]<<"_"<<it->GetIndex()[1]<<"_"<<it->GetSize()[0]<<"_"<<it->GetSize()[1]<<".tif";
+    typename WriterType::Pointer writer = WriterType::New();
+    writer->SetFileName(ss.str());
+    writer->SetInput(extractFilter->GetOutput());
+    // Datatype
+    gdalImageIO = dynamic_cast<GDALImageIO *>(writer->GetImageIO());
+    if(gdalImageIO.IsNotNull())
+    {
+      dataTypeStr = gdalImageIO->GetGdalPixelTypeAsString();
+    }
+
+    if (!m_AvailableRAM)
+      {
+      writer->SetNumberOfDivisionsTiledStreaming(0);
+      }
+    else
+      {
+      writer->SetAutomaticAdaptativeStreaming(m_AvailableRAM);
+      }
+
+    // Pipeline execution
+    try
+    {
+      writer->Update();
+    }
+    catch (itk::ExceptionObject& err)
+    {
+      std::stringstream message;
+      message << "ExceptionObject caught: " << err << std::endl;
+      mpiConfig->logError(message.str());
+      mpiConfig->abort(EXIT_FAILURE);
+    }
+  }
+
+  // MPI process synchronization
+  mpiConfig->barrier();
+
+  // Write VRT file
+  try
+  {
+    if (m_WriteVRT && (myRank == 0))
+    {
+      // VRT Filename
+      std::stringstream vrtfOut;
+      vrtfOut<< prefix << ".vrt";
+
+      // Data type
+      GDALDataType dataType;
+      dataType = GDALGetDataTypeByName(dataTypeStr.c_str());
+
+      int imageSizeY = img->GetLargestPossibleRegion().GetSize()[1];
+      int imageSizeX = img->GetLargestPossibleRegion().GetSize()[0];
+      const unsigned int nbBands = img->GetNumberOfComponentsPerPixel();
+
+      // Get VRT driver
+      GDALAllRegister();
+      GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("VRT");
+      if (driver == NULL) {
+         mpiConfig->logError("Error opening VRT driver.");
+         mpiConfig->abort(EXIT_FAILURE);
+      }
+
+      // Create output raster
+      GDALDataset *VRTOutput = driver->Create(vrtfOut.str().c_str(),
+            imageSizeX,
+            imageSizeY,
+            0,
+            dataType,
+            NULL); // No options
+      if (VRTOutput == NULL) {
+        mpiConfig->logError("driver->Create call failed.\n");
+        mpiConfig->abort(EXIT_FAILURE);
+      }
+
+      // Set GeoTransform
+      double gt[6];
+      gt[0] = img->GetOrigin()[0] - 0.5*img->GetSignedSpacing()[0];
+      gt[1] = img->GetSignedSpacing()[0];
+      gt[2] = 0.0;
+      gt[3] = img->GetOrigin()[1] - 0.5*img->GetSignedSpacing()[1];
+      gt[4] = 0.0;
+      gt[5] = img->GetSignedSpacing()[1];
+      VRTOutput->SetGeoTransform(gt);
+
+      // Set projection
+      OGRSpatialReference out_sr;
+      char *wkt = NULL;
+      out_sr.SetFromUserInput(img->GetProjectionRef().c_str());
+      out_sr.exportToWkt(&wkt);
+      VRTOutput->SetProjection(wkt);
+
+      for(unsigned int band = 1; band<=nbBands;++band)
+      {
+        VRTOutput->AddBand(dataType, NULL);
+
+        typename TImage::RegionType currentRegion;
+        for(unsigned int id=0; id < numberOfSplits; ++id)
+        {
+          currentRegion = streamingManager->GetSplit(id);
+          int tileSizeX = currentRegion.GetSize()[0];
+          int tileSizeY = currentRegion.GetSize()[1];
+          int tileIndexX = currentRegion.GetIndex()[0];
+          int tileIndexY = currentRegion.GetIndex()[1];
+          std::stringstream tileFileName;
+          tileFileName <<prefix<<"_"<<tileIndexX<<"_"<<tileIndexY<<"_"<<tileSizeX<<"_"<<tileSizeY<<".tif";
+          otbDebugMacro(<<tileFileName.str());
+
+          GDALDataset *dataset = (GDALDataset*) GDALOpen(tileFileName.str().c_str(), GA_ReadOnly);
+
+          VRTSourcedRasterBand *VRTBand = dynamic_cast<VRTSourcedRasterBand*> (VRTOutput->GetRasterBand(band));
+          VRTBand->AddSimpleSource(dataset->GetRasterBand(band),
+                                      0, //xOffSrc
+                                      0, //yOffSrc
+                                      tileSizeX, //xSizeSrc
+                                      tileSizeY, //ySizeSrc
+                                      tileIndexX, //xOffDest
+                                      tileIndexY, //yOffDest
+                                      tileSizeX, //xSizeDest
+                                      tileSizeY, //ySizeDest
+                                      "near",
+                                      VRT_NODATA_UNSET);
+        }
+
+      }
+
+      // Close
+      CPLFree(wkt);
+      GDALClose(VRTOutput);
+    }
+  }
+  catch (itk::ExceptionObject& err)
+  {
+    std::stringstream message;
+    message << "ExceptionObject caught: " << err << std::endl;
+    mpiConfig->logError(message.str());
+    mpiConfig->abort(EXIT_FAILURE);
+  }
+}
+
+
+} // end of namespace otb
+
+#endif
diff --git a/Modules/MPI/MPIVrtWriter/otb-module.cmake b/Modules/MPI/MPIVrtWriter/otb-module.cmake
index 0e95d97f7cbe02b27ba37be548a0106b3d70ea9e..f98ab957d207a9bf02f6f13c27002ad054582897 100644
--- a/Modules/MPI/MPIVrtWriter/otb-module.cmake
+++ b/Modules/MPI/MPIVrtWriter/otb-module.cmake
@@ -22,8 +22,9 @@ set(DOCUMENTATION "Provides a template function for MPI writing to a VRT file")
 
 otb_module(OTBMPIVrtWriter
   DEPENDS
-  OTBMPIConfig
-    OTBPanSharpening
+    OTBCommon
+    OTBMPIConfig
+#    OTBPanSharpening
     OTBProjection
     OTBInterpolation
     OTBTestKernel
diff --git a/Modules/ThirdParty/Boost/otb-module-init.cmake b/Modules/ThirdParty/Boost/otb-module-init.cmake
index a5f58041fb0bf5001d123ff2e0f772a26cff6d34..0a07bdfdbc7a9d33a03464925c95667a036bea59 100644
--- a/Modules/ThirdParty/Boost/otb-module-init.cmake
+++ b/Modules/ThirdParty/Boost/otb-module-init.cmake
@@ -31,3 +31,8 @@ if (BUILD_TESTING)
     message(STATUS "Found Boost components: unit_test_framework")
   endif()
 endif() #BUILD_TESTING
+
+if(WIN32)
+  # disable autolinking in boost
+	add_definitions( -DBOOST_ALL_NO_LIB )
+endif()
diff --git a/Modules/ThirdParty/Shark/include/otbSharkUtils.h b/Modules/ThirdParty/Shark/include/otbSharkUtils.h
index de3adf77401d0f131d2bd7d447627829b3df64ff..04c57b6d4e7f5a022b0c4fafa86ac41b134f690c 100644
--- a/Modules/ThirdParty/Shark/include/otbSharkUtils.h
+++ b/Modules/ThirdParty/Shark/include/otbSharkUtils.h
@@ -23,6 +23,7 @@
 
 #include <stdexcept>
 #include <string>
+#include <unordered_map>
 
 #if defined(__GNUC__) || defined(__clang__)
 #pragma GCC diagnostic push
@@ -127,6 +128,27 @@ template <class T> void ListSampleToSharkVector(const T * listSample, std::vecto
   assert(listSample != nullptr);
   ListSampleRangeToSharkVector(listSample,output,0, static_cast<unsigned int>(listSample->Size()));
 }
+
+/** Shark assumes that labels are 0 ... (nbClasses-1). This function modifies the labels contained in the input vector and returns a vector with size = nbClasses which allows the translation from the normalised labels to the new ones oldLabel = dictionary[newLabel].
+*/
+template <typename T> void NormalizeLabelsAndGetDictionary(std::vector<T>& labels, 
+                                                           std::vector<T>& dictionary)
+{
+  std::unordered_map<T, T> dictMap;
+  T labelCount{0};
+  for(const auto& l : labels)
+    {
+    if(dictMap.find(l)==dictMap.end())
+      dictMap.insert({l, labelCount++});
+    }
+  dictionary.resize(labelCount);
+  for(auto& l : labels)
+    {
+    auto newLabel = dictMap[l];
+    dictionary[newLabel] = l;
+    l = newLabel;
+    }
+}
   
 }
 }
diff --git a/Modules/Wrappers/ApplicationEngine/include/otbWrapperOutputImageParameter.h b/Modules/Wrappers/ApplicationEngine/include/otbWrapperOutputImageParameter.h
index b9bc451c3b3b9a59437ac306cf48929e9459d734..eac7df4cadad4e1c8376918f53401439748a6094 100644
--- a/Modules/Wrappers/ApplicationEngine/include/otbWrapperOutputImageParameter.h
+++ b/Modules/Wrappers/ApplicationEngine/include/otbWrapperOutputImageParameter.h
@@ -117,14 +117,8 @@ protected:
   /** Destructor */
   ~OutputImageParameter() override;
 
-  template <class TInputVectorImageType>
-    void SwitchVectorImageWrite();
-
-  template <class TInputVectorImageType>
-    void SwitchRGBImageWrite();
-
-  template <class TInputVectorImageType>
-    void SwitchRGBAImageWrite();
+  template <class TInput>
+    int SwitchInput(TInput *img);
 
   //FloatVectorImageType::Pointer m_Image;
   ImageBaseType::Pointer m_Image;
@@ -132,46 +126,28 @@ protected:
   ImagePixelType         m_PixelType;
   ImagePixelType         m_DefaultPixelType;
 
-  typedef otb::ImageFileWriter<UInt8VectorImageType>  VectorUInt8WriterType;
-  typedef otb::ImageFileWriter<Int16VectorImageType>  VectorInt16WriterType;
-  typedef otb::ImageFileWriter<UInt16VectorImageType> VectorUInt16WriterType;
-  typedef otb::ImageFileWriter<Int32VectorImageType>  VectorInt32WriterType;
-  typedef otb::ImageFileWriter<UInt32VectorImageType> VectorUInt32WriterType;
-  typedef otb::ImageFileWriter<FloatVectorImageType>  VectorFloatWriterType;
-  typedef otb::ImageFileWriter<DoubleVectorImageType> VectorDoubleWriterType;
-
-  typedef otb::ImageFileWriter<UInt8RGBAImageType>  RGBAUInt8WriterType;
-  typedef otb::ImageFileWriter<UInt8RGBImageType>   RGBUInt8WriterType;
-
-  typedef otb::ImageFileWriter<ComplexInt16VectorImageType>  ComplexVectorInt16WriterType;
-  typedef otb::ImageFileWriter<ComplexInt32VectorImageType>  ComplexVectorInt32WriterType;
-  typedef otb::ImageFileWriter<ComplexFloatVectorImageType>  ComplexVectorFloatWriterType;
-  typedef otb::ImageFileWriter<ComplexDoubleVectorImageType> ComplexVectorDoubleWriterType;
-
-  VectorUInt8WriterType::Pointer  m_VectorUInt8Writer;
-  VectorInt16WriterType::Pointer  m_VectorInt16Writer;
-  VectorUInt16WriterType::Pointer m_VectorUInt16Writer;
-  VectorInt32WriterType::Pointer  m_VectorInt32Writer;
-  VectorUInt32WriterType::Pointer m_VectorUInt32Writer;
-  VectorFloatWriterType::Pointer  m_VectorFloatWriter;
-  VectorDoubleWriterType::Pointer m_VectorDoubleWriter;
-
-  RGBUInt8WriterType::Pointer   m_RGBUInt8Writer;
-  RGBAUInt8WriterType::Pointer  m_RGBAUInt8Writer;
-
-  ComplexVectorInt16WriterType::Pointer  m_ComplexVectorInt16Writer;
-  ComplexVectorInt32WriterType::Pointer  m_ComplexVectorInt32Writer;
-  ComplexVectorFloatWriterType::Pointer  m_ComplexVectorFloatWriter;
-  ComplexVectorDoubleWriterType::Pointer m_ComplexVectorDoubleWriter;
-
 private:
   OutputImageParameter(const Parameter &); //purposely not implemented
   void operator =(const Parameter&); //purposely not implemented
 
   unsigned int                  m_RAMValue;
 
+  itk::ProcessObject::Pointer m_Caster;
+
+  itk::ProcessObject::Pointer m_Writer;
+
 }; // End class OutputImage Parameter
 
+// Declare specialisation for UInt8RGBAImageType
+template <>
+int
+OutputImageParameter::SwitchInput(UInt8RGBAImageType *img);
+
+// Declare specialisation for UInt8RGBImageType
+template <>
+int
+OutputImageParameter::SwitchInput(UInt8RGBImageType *img);
+
 } // End namespace Wrapper
 } // End namespace otb
 
diff --git a/Modules/Wrappers/ApplicationEngine/src/otbWrapperApplication.cxx b/Modules/Wrappers/ApplicationEngine/src/otbWrapperApplication.cxx
index 7aec4f87620aa6917040ee282d4c5e0110848d53..31f46f3a019e6fe335ebcae5535b4b6284bc2f4a 100644
--- a/Modules/Wrappers/ApplicationEngine/src/otbWrapperApplication.cxx
+++ b/Modules/Wrappers/ApplicationEngine/src/otbWrapperApplication.cxx
@@ -686,7 +686,6 @@ int Application::ExecuteAndWriteOutput()
 
           if(outputParam!=ITK_NULLPTR)
             {
-            outputParam->InitializeWriters();
             std::string checkReturn = outputParam->CheckFileName(true);
             if (!checkReturn.empty())
               {
@@ -696,6 +695,7 @@ int Application::ExecuteAndWriteOutput()
               {
               outputParam->SetRAMValue(ram);
               }
+            outputParam->InitializeWriters();
             std::ostringstream progressId;
             progressId << "Writing " << outputParam->GetFileName() << "...";
             AddProcess(outputParam->GetWriter(), progressId.str());
diff --git a/Modules/Wrappers/ApplicationEngine/src/otbWrapperComplexOutputImageParameter.cxx b/Modules/Wrappers/ApplicationEngine/src/otbWrapperComplexOutputImageParameter.cxx
index b97171cb71d9b270f07437cf22fa5116bc0b0d7c..4218836ed7eb3091f1cf1a451c165d9bb77a040a 100644
--- a/Modules/Wrappers/ApplicationEngine/src/otbWrapperComplexOutputImageParameter.cxx
+++ b/Modules/Wrappers/ApplicationEngine/src/otbWrapperComplexOutputImageParameter.cxx
@@ -259,20 +259,22 @@ ComplexOutputImageParameter::GetWriter()
     case ComplexImagePixelType_int16:
     {
       writer = m_ComplexVectorInt16Writer;
+      break;
     }
     case ComplexImagePixelType_int32:
     {
       writer = m_ComplexVectorInt32Writer;
+      break;
     }
     case ComplexImagePixelType_float:
     {
       writer = m_ComplexVectorFloatWriter;
-    break;
+      break;
     }
     case ComplexImagePixelType_double:
     {
       writer = m_ComplexVectorDoubleWriter;
-    break;
+      break;
     }
     }
   return writer;
diff --git a/Modules/Wrappers/ApplicationEngine/src/otbWrapperOutputImageParameter.cxx b/Modules/Wrappers/ApplicationEngine/src/otbWrapperOutputImageParameter.cxx
index 09396447765610f1b1f2dd903943b9a2cc116922..9e6e33efdac2c1223830dad9db754e7110777817 100644
--- a/Modules/Wrappers/ApplicationEngine/src/otbWrapperOutputImageParameter.cxx
+++ b/Modules/Wrappers/ApplicationEngine/src/otbWrapperOutputImageParameter.cxx
@@ -42,7 +42,9 @@ namespace Wrapper
 OutputImageParameter::OutputImageParameter()
   : m_PixelType(ImagePixelType_float),
     m_DefaultPixelType(ImagePixelType_float),
-    m_RAMValue(0)
+    m_RAMValue(0),
+    m_Caster(nullptr),
+    m_Writer(nullptr)
 {
   this->SetName("Output Image");
   this->SetKey("out");
@@ -149,35 +151,102 @@ OutputImageParameter::ConvertStringToPixelType(const std::string &value, ImagePi
 
 void OutputImageParameter::InitializeWriters()
 {
-  m_VectorUInt8Writer = VectorUInt8WriterType::New();
-  m_VectorInt16Writer = VectorInt16WriterType::New();
-  m_VectorUInt16Writer = VectorUInt16WriterType::New();
-  m_VectorInt32Writer = VectorInt32WriterType::New();
-  m_VectorUInt32Writer = VectorUInt32WriterType::New();
-  m_VectorFloatWriter = VectorFloatWriterType::New();
-  m_VectorDoubleWriter = VectorDoubleWriterType::New();
-
-  m_RGBUInt8Writer = RGBUInt8WriterType::New();
-  m_RGBAUInt8Writer = RGBAUInt8WriterType::New();
-
-  m_ComplexVectorInt16Writer = ComplexVectorInt16WriterType::New();
-  m_ComplexVectorInt32Writer = ComplexVectorInt32WriterType::New();
-  m_ComplexVectorFloatWriter = ComplexVectorFloatWriterType::New();
-  m_ComplexVectorDoubleWriter = ComplexVectorDoubleWriterType::New();
+  ImageBaseType* imgBase = m_Image.GetPointer();
+  // Guess the image type
+  std::string className(m_Image->GetNameOfClass());
+  if (className == "VectorImage")
+    {
+    UInt8VectorImageType* imgUInt8 = dynamic_cast<UInt8VectorImageType*>(imgBase);
+    if (imgUInt8 && SwitchInput(imgUInt8)) return;
+
+    Int16VectorImageType* imgInt16 = dynamic_cast<Int16VectorImageType*>(imgBase);
+    if (imgInt16 && SwitchInput(imgInt16)) return;
+
+    UInt16VectorImageType* imgUInt16 = dynamic_cast<UInt16VectorImageType*>(imgBase);
+    if (imgUInt16 && SwitchInput(imgUInt16)) return;
+
+    Int32VectorImageType* imgInt32 = dynamic_cast<Int32VectorImageType*>(imgBase);
+    if (imgInt32 && SwitchInput(imgInt32)) return;
+
+    UInt32VectorImageType* imgUInt32 = dynamic_cast<UInt32VectorImageType*>(imgBase);
+    if (imgUInt32 && SwitchInput(imgUInt32)) return;
+
+    FloatVectorImageType* imgFloat = dynamic_cast<FloatVectorImageType*>(imgBase);
+    if (imgFloat && SwitchInput(imgFloat)) return;
+
+    DoubleVectorImageType* imgDouble = dynamic_cast<DoubleVectorImageType*>(imgBase);
+    if (imgDouble && SwitchInput(imgDouble)) return;
+
+    ComplexInt16VectorImageType* imgCInt16 = dynamic_cast<ComplexInt16VectorImageType*>(imgBase);
+    if (imgCInt16 && SwitchInput(imgCInt16)) return;
+
+    ComplexInt32VectorImageType* imgCInt32 = dynamic_cast<ComplexInt32VectorImageType*>(imgBase);
+    if (imgCInt32 && SwitchInput(imgCInt32)) return;
+
+    ComplexFloatVectorImageType* imgCFloat = dynamic_cast<ComplexFloatVectorImageType*>(imgBase);
+    if (imgCFloat && SwitchInput(imgCFloat)) return;
+
+    ComplexDoubleVectorImageType* imgCDouble = dynamic_cast<ComplexDoubleVectorImageType*>(imgBase);
+    if (imgCDouble && SwitchInput(imgCDouble)) return;
+    }
+  else
+    {
+    UInt8ImageType* imgUInt8 = dynamic_cast<UInt8ImageType*>(imgBase);
+    if (imgUInt8 && SwitchInput(imgUInt8)) return;
+
+    Int16ImageType* imgInt16 = dynamic_cast<Int16ImageType*>(imgBase);
+    if (imgInt16 && SwitchInput(imgInt16)) return;
+
+    UInt16ImageType* imgUInt16 = dynamic_cast<UInt16ImageType*>(imgBase);
+    if (imgUInt16 && SwitchInput(imgUInt16)) return;
+
+    Int32ImageType* imgInt32 = dynamic_cast<Int32ImageType*>(imgBase);
+    if (imgInt32 && SwitchInput(imgInt32)) return;
+
+    UInt32ImageType* imgUInt32 = dynamic_cast<UInt32ImageType*>(imgBase);
+    if (imgUInt32 && SwitchInput(imgUInt32)) return;
+
+    FloatImageType* imgFloat = dynamic_cast<FloatImageType*>(imgBase);
+    if (imgFloat && SwitchInput(imgFloat)) return;
+
+    DoubleImageType* imgDouble = dynamic_cast<DoubleImageType*>(imgBase);
+    if (imgDouble && SwitchInput(imgDouble)) return;
+
+    ComplexInt16ImageType* imgCInt16 = dynamic_cast<ComplexInt16ImageType*>(imgBase);
+    if (imgCInt16 && SwitchInput(imgCInt16)) return;
+
+    ComplexInt32ImageType* imgCInt32 = dynamic_cast<ComplexInt32ImageType*>(imgBase);
+    if (imgCInt32 && SwitchInput(imgCInt32)) return;
+
+    ComplexFloatImageType* imgCFloat = dynamic_cast<ComplexFloatImageType*>(imgBase);
+    if (imgCFloat && SwitchInput(imgCFloat)) return;
+
+    ComplexDoubleImageType* imgCDouble = dynamic_cast<ComplexDoubleImageType*>(imgBase);
+    if (imgCDouble && SwitchInput(imgCDouble)) return;
+
+    UInt8RGBImageType* imgRGB = dynamic_cast<UInt8RGBImageType*>(imgBase);
+    if (imgRGB && SwitchInput(imgRGB)) return;
+
+    UInt8RGBAImageType* imgRGBA = dynamic_cast<UInt8RGBAImageType*>(imgBase);
+    if (imgRGBA && SwitchInput(imgRGBA)) return;
+    }
+
+  itkExceptionMacro("Unknown image type");
 }
 
 
 template <typename TInput, typename TOutput> 
-void 
-ClampAndWriteVectorImage( itk::ImageBase<2> * in ,
-                    otb::ImageFileWriter<TOutput> * writer , 
+std::pair<itk::ProcessObject::Pointer,itk::ProcessObject::Pointer> 
+ClampAndWriteVectorImage( TInput * in ,
                     const std::string & filename , 
                     const unsigned int & ramValue )
 {
+  std::pair<itk::ProcessObject::Pointer,itk::ProcessObject::Pointer> ret;
   typedef ClampImageFilter < TInput , TOutput > ClampFilterType;
   typename ClampFilterType::Pointer clampFilter ( ClampFilterType::New() );
 
-  clampFilter->SetInput( dynamic_cast<TInput*>(in));
+  clampFilter->SetInput( in);
+  ret.first = clampFilter.GetPointer();
   
   bool useStandardWriter = true;
 
@@ -194,8 +263,14 @@ ClampAndWriteVectorImage( itk::ImageBase<2> * in ,
 
     if(extension == ".vrt")
       {
-      // Use the WriteMPI function
-      WriteMPI(clampFilter->GetOutput(),filename,ramValue);      
+      // Use the MPIVrtWriter
+      typedef otb::MPIVrtWriter<TOutput> VRTWriterType;
+
+      typename VRTWriterType::Pointer vrtWriter = VRTWriterType::New();
+      vrtWriter->SetInput(clampFilter->GetOutput());
+      vrtWriter->SetFileName(filename);
+      vrtWriter->SetAvailableRAM(ramValue);
+      ret.second = vrtWriter.GetPointer();
       }
     #ifdef OTB_USE_SPTW
     else if (extension == ".tif")
@@ -206,9 +281,8 @@ ClampAndWriteVectorImage( itk::ImageBase<2> * in ,
       typename SPTWriterType::Pointer sptWriter = SPTWriterType::New();
       sptWriter->SetFileName(filename);
       sptWriter->SetInput(clampFilter->GetOutput());
-      sptWriter->SetAutomaticAdaptativeStreaming(ramValue);
       sptWriter->GetStreamingManager()->SetDefaultRAM(ramValue);
-      sptWriter->Update();
+      ret.second = sptWriter.GetPointer();
       }
     
     #endif
@@ -223,114 +297,110 @@ ClampAndWriteVectorImage( itk::ImageBase<2> * in ,
   
   if(useStandardWriter)
     {
+    typename otb::ImageFileWriter<TOutput>::Pointer writer =
+      otb::ImageFileWriter<TOutput>::New();
     writer->SetFileName( filename );
     writer->SetInput(clampFilter->GetOutput());
     writer->GetStreamingManager()->SetDefaultRAM(ramValue);
-    writer->Update();
+    ret.second = writer.GetPointer();
     }
+
+  return ret;
 }
 
 template <class TInput>
-void
-OutputImageParameter::SwitchVectorImageWrite()
-  {
+int
+OutputImageParameter::SwitchInput(TInput *img)
+{
+  if (! img) return 0;
+
+  std::pair<itk::ProcessObject::Pointer,itk::ProcessObject::Pointer> ret;
   switch(m_PixelType )
     {
     case ImagePixelType_uint8:
     {
-    ClampAndWriteVectorImage< TInput , UInt8VectorImageType > (
-      m_Image ,
-      m_VectorUInt8Writer ,
+    ret = ClampAndWriteVectorImage< TInput , UInt8VectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue );
     break;
     }
     case ImagePixelType_int16:
     {
-    ClampAndWriteVectorImage< TInput , Int16VectorImageType > (
-      m_Image ,
-      m_VectorInt16Writer ,
+    ret = ClampAndWriteVectorImage< TInput , Int16VectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue );
     break;
     }
     case ImagePixelType_uint16:
     {
-    ClampAndWriteVectorImage< TInput , UInt16VectorImageType > (
-      m_Image ,
-      m_VectorUInt16Writer ,
+    ret = ClampAndWriteVectorImage< TInput , UInt16VectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue );
     break;
     }
     case ImagePixelType_int32:
     {
-    ClampAndWriteVectorImage< TInput , Int32VectorImageType > (
-      m_Image ,
-      m_VectorInt32Writer ,
+    ret = ClampAndWriteVectorImage< TInput , Int32VectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue );
     break;
     }
     case ImagePixelType_uint32:
     {
-    ClampAndWriteVectorImage< TInput , UInt32VectorImageType > (
-      m_Image ,
-      m_VectorUInt32Writer ,
+    ret = ClampAndWriteVectorImage< TInput , UInt32VectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue );
     break;
     }
     case ImagePixelType_float:
     {
-    ClampAndWriteVectorImage< TInput , FloatVectorImageType > (
-      m_Image ,
-      m_VectorFloatWriter ,
+    ret = ClampAndWriteVectorImage< TInput , FloatVectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue );
     break;
     }
     case ImagePixelType_double:
     {
-    ClampAndWriteVectorImage< TInput , DoubleVectorImageType > (
-      m_Image ,
-      m_VectorDoubleWriter ,
+    ret = ClampAndWriteVectorImage< TInput , DoubleVectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue );
     break;
     }
     case ImagePixelType_cint16:
     {
-    ClampAndWriteVectorImage < TInput , ComplexInt16VectorImageType > (
-      m_Image ,
-      m_ComplexVectorInt16Writer ,
+    ret = ClampAndWriteVectorImage < TInput , ComplexInt16VectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue ); 
     break;
     }
     case ImagePixelType_cint32:
     {
-    ClampAndWriteVectorImage < TInput , ComplexInt32VectorImageType > (
-      m_Image ,
-      m_ComplexVectorInt32Writer ,
+    ret = ClampAndWriteVectorImage < TInput , ComplexInt32VectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue ); 
     break;
     }
     case ImagePixelType_cfloat:
     {
-    ClampAndWriteVectorImage < TInput , ComplexFloatVectorImageType > (
-      m_Image ,
-      m_ComplexVectorFloatWriter ,
+    ret = ClampAndWriteVectorImage < TInput , ComplexFloatVectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue ); 
     break;
     }
     case ImagePixelType_cdouble:
     {
-    ClampAndWriteVectorImage < TInput , ComplexDoubleVectorImageType > (
-      m_Image ,
-      m_ComplexVectorDoubleWriter ,
+    ret = ClampAndWriteVectorImage < TInput , ComplexDoubleVectorImageType > (
+      img ,
       m_FileName ,
       m_RAMValue );
     break;
@@ -338,256 +408,27 @@ OutputImageParameter::SwitchVectorImageWrite()
     default:
       break;
     }
-  }
-
-template <class TInputRGBAImageType>
-void
-OutputImageParameter::SwitchRGBAImageWrite()
-  {
-  if( m_PixelType == ImagePixelType_uint8 )
-    {
-    m_RGBAUInt8Writer->SetFileName( this->GetFileName() );
-    m_RGBAUInt8Writer->SetInput(dynamic_cast<UInt8RGBAImageType*>(m_Image.GetPointer()) );
-    m_RGBAUInt8Writer->GetStreamingManager()->SetDefaultRAM(m_RAMValue);
-    m_RGBAUInt8Writer->Update();
-    }
-   else
-     itkExceptionMacro("Unknown PixelType for RGBA Image. Only uint8 is supported.");
-  }
-
-template <class TInputRGBImageType>
-void
-OutputImageParameter::SwitchRGBImageWrite()
-  {
-   if( m_PixelType == ImagePixelType_uint8 )
-    {
-    m_RGBUInt8Writer->SetFileName( this->GetFileName() );
-    m_RGBUInt8Writer->SetInput(dynamic_cast<UInt8RGBImageType*>(m_Image.GetPointer()) );
-    m_RGBUInt8Writer->GetStreamingManager()->SetDefaultRAM(m_RAMValue);
-    m_RGBUInt8Writer->Update();
-    }
-   else
-     itkExceptionMacro("Unknown PixelType for RGB Image. Only uint8 is supported.");
-  }
+  // Save the caster and writer
+  m_Caster = ret.first;
+  m_Writer = ret.second;
+  return 1;
+}
 
 void
 OutputImageParameter::Write()
 {
-  m_Image->UpdateOutputInformation();
+  m_Writer->Update();
 
-  if (dynamic_cast<UInt8ImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<UInt8ImageType>();
-    }
-  else if (dynamic_cast<Int16ImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<Int16ImageType>();
-    }
-  else if (dynamic_cast<UInt16ImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<UInt16ImageType>();
-    }
-  else if (dynamic_cast<Int32ImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<Int32ImageType>();
-    }
-  else if (dynamic_cast<UInt32ImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<UInt32ImageType>();
-    }
-  else if (dynamic_cast<FloatImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<FloatImageType>();
-    }
-  else if (dynamic_cast<DoubleImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<DoubleImageType>();
-    }
-  else if (dynamic_cast<ComplexInt16ImageType*>(m_Image.GetPointer()) )
-    {
-    SwitchVectorImageWrite<ComplexInt16ImageType>();
-    }
-  else if (dynamic_cast<ComplexInt32ImageType*>(m_Image.GetPointer()) )
-    {
-    SwitchVectorImageWrite<ComplexInt32ImageType>();
-    }
-  else if (dynamic_cast<ComplexFloatImageType*>(m_Image.GetPointer()) )
-    {
-    SwitchVectorImageWrite<ComplexFloatImageType>();
-    }
-  else if (dynamic_cast<ComplexDoubleImageType*>(m_Image.GetPointer()) )
-    {
-    SwitchVectorImageWrite<ComplexDoubleImageType>();
-    }
-  else if (dynamic_cast<UInt8VectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<UInt8VectorImageType>();
-    }
-  else if (dynamic_cast<Int16VectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<Int16VectorImageType>();
-    }
-  else if (dynamic_cast<UInt16VectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<UInt16VectorImageType>();
-    }
-  else if (dynamic_cast<Int32VectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<Int32VectorImageType>();
-    }
-  else if (dynamic_cast<UInt32VectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<UInt32VectorImageType>();
-    }
-  else if (dynamic_cast<FloatVectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<FloatVectorImageType>();
-    }
-  else if (dynamic_cast<DoubleVectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<DoubleVectorImageType>();
-    }
-  else if (dynamic_cast<ComplexInt16VectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<ComplexInt16VectorImageType>();
-    }
-  else if (dynamic_cast<ComplexInt32VectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<ComplexInt32VectorImageType>();
-    }
-  else if (dynamic_cast<ComplexFloatVectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<ComplexFloatVectorImageType>();
-    }
-  else if (dynamic_cast<ComplexDoubleVectorImageType*>(m_Image.GetPointer()))
-    {
-    SwitchVectorImageWrite<ComplexDoubleVectorImageType>();
-    }
-  else if (dynamic_cast<UInt8RGBImageType*>(m_Image.GetPointer()))
-    {
-    SwitchRGBImageWrite<UInt8RGBImageType>();
-    }
-  else if (dynamic_cast<UInt8RGBAImageType*>(m_Image.GetPointer()))
-    {
-    SwitchRGBAImageWrite<UInt8RGBAImageType>();
-    }
-  else
-    {
-    itkExceptionMacro("Unknown image type");
-    }
-  }
+  // Clean internal filters
+  m_Caster = nullptr;
+  m_Writer = nullptr;
+}
 
 
 itk::ProcessObject*
 OutputImageParameter::GetWriter()
 {
-  int type = 1;
-  // 0 : image
-  // 1 : VectorImage
-  // 2 : RGBAImage
-  // 3 : RGBImage
-  itk::ProcessObject* writer = ITK_NULLPTR;
-
-  if (dynamic_cast<UInt8RGBAImageType*> (m_Image.GetPointer()))
-    {
-    type = 2;
-    writer = m_RGBAUInt8Writer;
-    itkWarningMacro("UInt8RGBAImageType will be saved in UInt8 format.");
-    return writer;
-    }
-  else if (dynamic_cast<UInt8RGBImageType*> (m_Image.GetPointer()))
-    {
-    type = 3;
-    writer = m_RGBUInt8Writer;
-    itkWarningMacro("UInt8RGBImageType will be saved in UInt8 format.");
-    return writer;
-    }
-  
-  switch (GetPixelType())
-    {
-    case ImagePixelType_uint8:
-      {
-      switch(type)
-        {
-        case 1:
-          writer = m_VectorUInt8Writer;
-          break;
-        case 2:
-          writer = m_RGBAUInt8Writer;
-          break;
-        default:
-          writer = m_RGBUInt8Writer;
-          break;
-        }
-      break;
-      }
-    case ImagePixelType_int16:
-      {
-      if (type == 1)
-        writer = m_VectorInt16Writer;
-      break;
-      }
-    case ImagePixelType_uint16:
-      {
-      if (type == 1)
-        writer = m_VectorUInt16Writer;
-      break;
-      }
-    case ImagePixelType_int32:
-      {
-      if (type == 1)
-        writer = m_VectorInt32Writer;
-      break;
-      }
-    case ImagePixelType_uint32:
-      {
-      if (type == 1)
-        writer = m_VectorUInt32Writer;
-      break;
-      }
-    case ImagePixelType_float:
-      {
-      if (type == 1)
-        writer = m_VectorFloatWriter;
-      break;
-      }
-    case ImagePixelType_double:
-      {
-      if (type == 1)
-        writer = m_VectorDoubleWriter;
-      break;
-      }
-    case ImagePixelType_cint16:
-      {
-      if( type == 1 )
-        writer = m_ComplexVectorInt16Writer;
-      break;
-      }
-    case ImagePixelType_cint32:
-      {
-      if( type == 1 )
-        writer = m_ComplexVectorInt32Writer;
-      break;
-      }
-    case ImagePixelType_cfloat:
-      {
-      if( type == 1 )
-        writer = m_ComplexVectorFloatWriter;
-      break;
-      }
-    case ImagePixelType_cdouble:
-      {
-      if( type == 1 )
-        writer = m_ComplexVectorDoubleWriter;
-      break;
-      }
-    }
-  if (ITK_NULLPTR == writer)
-    {
-    itkExceptionMacro("Unknown Writer type.");
-    }
-
-  return writer;
+  return m_Writer;
 }
 
 ImageBaseType*
@@ -651,5 +492,45 @@ OutputImageParameter::CheckFileName(bool fixMissingExtension)
   return ret;
 }
 
+// Specialization for UInt8RGBAImageType
+template <>
+int
+OutputImageParameter::SwitchInput(UInt8RGBAImageType *img)
+{
+  if (! img) return 0;
+  if( m_PixelType == ImagePixelType_uint8 )
+    {
+    typename otb::ImageFileWriter<UInt8RGBAImageType>::Pointer writer =
+      otb::ImageFileWriter<UInt8RGBAImageType>::New();
+    writer->SetFileName( this->GetFileName() );
+    writer->SetInput(img);
+    writer->GetStreamingManager()->SetDefaultRAM(m_RAMValue);
+    m_Writer = writer.GetPointer();
+    }
+   else
+     itkExceptionMacro("Unknown PixelType for RGBA Image. Only uint8 is supported.");
+  return 1;
+}
+
+// Specialization for UInt8RGBImageType
+template <>
+int
+OutputImageParameter::SwitchInput(UInt8RGBImageType *img)
+{
+  if (! img) return 0;
+  if( m_PixelType == ImagePixelType_uint8 )
+    {
+    typename otb::ImageFileWriter<UInt8RGBImageType>::Pointer writer =
+      otb::ImageFileWriter<UInt8RGBImageType>::New();
+    writer->SetFileName( this->GetFileName() );
+    writer->SetInput(img);
+    writer->GetStreamingManager()->SetDefaultRAM(m_RAMValue);
+    m_Writer = writer.GetPointer();
+    }
+   else
+     itkExceptionMacro("Unknown PixelType for RGB Image. Only uint8 is supported.");
+  return 1;
+}
+
 }
 }
diff --git a/SuperBuild/CMake/External_shark.cmake b/SuperBuild/CMake/External_shark.cmake
index ce8486db084935352b4266fc384f40be3604a29c..33934d4bb943dce00faceb9d6910534197a2342e 100644
--- a/SuperBuild/CMake/External_shark.cmake
+++ b/SuperBuild/CMake/External_shark.cmake
@@ -30,8 +30,8 @@ ADD_SUPERBUILD_CMAKE_VAR(SHARK BOOST_LIBRARYDIR)
 
 ExternalProject_Add(SHARK
   PREFIX SHARK
-  URL "https://github.com/Shark-ML/Shark/archive/349f29bd71c370e0f88f7fc9aa66fa5c4768fcb0.zip"
-  URL_MD5 d6e4310f943e8dda4a0151612b5c62ce
+  URL "https://github.com/Shark-ML/Shark/archive/2fd55e2b83f0666d05b403b291712668f4b76a13.zip"
+  URL_MD5 863bb5f0d94b01be5292867beb05a0bb
   SOURCE_DIR ${SHARK_SB_SRC}
   BINARY_DIR ${SHARK_SB_BUILD_DIR}
   INSTALL_DIR ${SB_INSTALL_PREFIX}
@@ -45,6 +45,7 @@ ExternalProject_Add(SHARK
   -DENABLE_HDF5:BOOL=OFF
   -DENABLE_CBLAS:BOOL=OFF
   -DENABLE_OPENMP:BOOL=${OTB_USE_OPENMP}
+  -DSHARK_INSTALL_LIB_DIR:STRING=lib/
   ${SHARK_SB_CONFIG}
   CMAKE_COMMAND ${SB_CMAKE_COMMAND}
   LOG_DOWNLOAD 1
diff --git a/SuperBuild/patches/SHARK/shark-2-ext-num-literals-all.diff b/SuperBuild/patches/SHARK/shark-2-ext-num-literals-all.diff
new file mode 100644
index 0000000000000000000000000000000000000000..0b964c1b9ada7aa4409f0f032285a70723caacfe
--- /dev/null
+++ b/SuperBuild/patches/SHARK/shark-2-ext-num-literals-all.diff
@@ -0,0 +1,13 @@
+diff -burN Shark.orig/CMakeLists.txt Shark/CMakeLists.txt
+--- Shark.orig/CMakeLists.txt	2018-02-05 18:04:58.012612932 +0100
++++ Shark/CMakeLists.txt	2018-02-05 18:20:50.032233165 +0100
+@@ -415,6 +415,9 @@
+ #####################################################################
+ #                       General Path settings
+ #####################################################################
++if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
++  add_definitions(-fext-numeric-literals)
++endif()
+ include_directories( ${shark_SOURCE_DIR}/include )
+ include_directories( ${shark_BINARY_DIR}/include )
+ add_subdirectory( include )
diff --git a/SuperBuild/patches/SHARK/shark-2-find-boost-all.diff b/SuperBuild/patches/SHARK/shark-2-find-boost-all.diff
deleted file mode 100644
index a97c1ac4afd1f56118fdba14cf7b993755bb5c00..0000000000000000000000000000000000000000
--- a/SuperBuild/patches/SHARK/shark-2-find-boost-all.diff
+++ /dev/null
@@ -1,16 +0,0 @@
-diff -burN Shark-349f29bd71c370e0f88f7fc9aa66fa5c4768fcb0.orig/CMakeLists.txt Shark-349f29bd71c370e0f88f7fc9aa66fa5c4768fcb0/CMakeLists.txt
---- Shark-349f29bd71c370e0f88f7fc9aa66fa5c4768fcb0.orig/CMakeLists.txt	2017-08-22 11:31:50.472052695 +0200
-+++ Shark-349f29bd71c370e0f88f7fc9aa66fa5c4768fcb0/CMakeLists.txt	2017-08-22 11:32:36.448358789 +0200
-@@ -141,10 +141,8 @@
- 
- find_package( 
- 	Boost 1.48.0 REQUIRED COMPONENTS
--	system date_time filesystem
--	program_options serialization thread
--	unit_test_framework
--)
-+	serialization
-+	)
- 
- if(NOT Boost_FOUND)
- 	message(FATAL_ERROR "Please make sure Boost 1.48.0 is installed on your system")