diff --git a/Data/Baseline/OTB/Images/apTvUtQuicklookWithGCP.tif b/Data/Baseline/OTB/Images/apTvUtQuicklookWithGCP.tif
new file mode 100644
index 0000000000000000000000000000000000000000..85a15cbe7ebc4f29122eacedefcf5138a236e4bc
--- /dev/null
+++ b/Data/Baseline/OTB/Images/apTvUtQuicklookWithGCP.tif
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:590e2e3df54c12833852bf26c00961cf24fb7fca69a25077423c8e56316bfb0b
+size 10580
diff --git a/Data/Input/Scene.png b/Data/Input/Scene.png
deleted file mode 100644
index 85e1f78f536c0ac65fc652f28005d96ef63128fe..0000000000000000000000000000000000000000
--- a/Data/Input/Scene.png
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:1e249ea80956b59a9009f49dcf18292804d89a0bc67d71e957fc51f8e81fcf60
-size 55503
diff --git a/Modules/Applications/AppClassification/include/otbVectorPrediction.h b/Modules/Applications/AppClassification/include/otbVectorPrediction.h
index 3e54cb3d0879eb3b40aed1acf3491c0a31794323..7441863f7b632f85f57f8b3da7c015b99477ab93 100644
--- a/Modules/Applications/AppClassification/include/otbVectorPrediction.h
+++ b/Modules/Applications/AppClassification/include/otbVectorPrediction.h
@@ -104,9 +104,12 @@ private:
   /** Normalize a list sample using the statistic file given  */
   typename ListSampleType::Pointer NormalizeListSample(ListSampleType::Pointer input);
 
-  /** Create the output DataSource, in update mode the input layer is buffered and the input
-   * data source is re opened in update mode. */
-  otb::ogr::DataSource::Pointer CreateOutputDataSource(otb::ogr::DataSource::Pointer source, otb::ogr::Layer& layer, bool updateMode);
+  /** Update the output DataSource : the input layer is buffered and the input data source is re opened in update mode. */
+  otb::ogr::DataSource::Pointer ReopenDataSourceInUpdateMode(ogr::DataSource::Pointer source, ogr::Layer& layer,
+                                                             ogr::DataSource::Pointer buffer);
+
+  /** Create the output DataSource. */
+  otb::ogr::DataSource::Pointer CreateOutputDataSource(ogr::DataSource::Pointer source, ogr::Layer& layer);
 
   /** Add a prediction field in the output layer if it does not exist.
    * If computeConfidenceMap evaluates to true a confidence field will be
diff --git a/Modules/Applications/AppClassification/include/otbVectorPrediction.hxx b/Modules/Applications/AppClassification/include/otbVectorPrediction.hxx
index ec99e681467b473f820404ea88dffd55c7ded187..b90b50f476885e0df9081456cf435173656443b6 100644
--- a/Modules/Applications/AppClassification/include/otbVectorPrediction.hxx
+++ b/Modules/Applications/AppClassification/include/otbVectorPrediction.hxx
@@ -147,41 +147,39 @@ typename VectorPrediction<RegressionMode>::ListSampleType::Pointer VectorPredict
 
 
 template <bool                RegressionMode>
-otb::ogr::DataSource::Pointer VectorPrediction<RegressionMode>::CreateOutputDataSource(otb::ogr::DataSource::Pointer source, otb::ogr::Layer& layer,
-                                                                                       bool updateMode)
+otb::ogr::DataSource::Pointer VectorPrediction<RegressionMode>::ReopenDataSourceInUpdateMode(ogr::DataSource::Pointer source, ogr::Layer& layer,
+                                                                                             ogr::DataSource::Pointer buffer)
 {
   ogr::DataSource::Pointer output;
-  ogr::DataSource::Pointer buffer = ogr::DataSource::New();
-  if (updateMode)
-  {
-    // Update mode
-    otbAppLogINFO("Update input vector data.");
-    // fill temporary buffer for the transfer
-    otb::ogr::Layer inputLayer = layer;
-    layer                      = buffer->CopyLayer(inputLayer, std::string("Buffer"));
-    // close input data source
-    source->Clear();
-    // Re-open input data source in update mode
-    output = otb::ogr::DataSource::New(GetParameterString("in"), otb::ogr::DataSource::Modes::Update_LayerUpdate);
-  }
-  else
+  // Update mode
+  otbAppLogINFO("Update input vector data.");
+  // fill temporary buffer for the transfer
+  otb::ogr::Layer inputLayer = layer;
+  layer                      = buffer->CopyLayer(inputLayer, std::string("Buffer"));
+  // close input data source
+  source->Clear();
+  // Re-open input data source in update mode
+  output = otb::ogr::DataSource::New(GetParameterString("in"), otb::ogr::DataSource::Modes::Update_LayerUpdate);
+  return output;
+}
+
+template <bool                RegressionMode>
+otb::ogr::DataSource::Pointer VectorPrediction<RegressionMode>::CreateOutputDataSource(ogr::DataSource::Pointer source, ogr::Layer& layer)
+{
+  ogr::DataSource::Pointer output;
+  // Create new OGRDataSource
+  output                   = ogr::DataSource::New(GetParameterString("out"), ogr::DataSource::Modes::Overwrite);
+  otb::ogr::Layer newLayer = output->CreateLayer(GetParameterString("out"), const_cast<OGRSpatialReference*>(layer.GetSpatialRef()), layer.GetGeomType());
+  // Copy existing fields
+  OGRFeatureDefn& inLayerDefn = layer.GetLayerDefn();
+  for (int k = 0; k < inLayerDefn.GetFieldCount(); k++)
   {
-    // Create new OGRDataSource
-    output                   = ogr::DataSource::New(GetParameterString("out"), ogr::DataSource::Modes::Overwrite);
-    otb::ogr::Layer newLayer = output->CreateLayer(GetParameterString("out"), const_cast<OGRSpatialReference*>(layer.GetSpatialRef()), layer.GetGeomType());
-    // Copy existing fields
-    OGRFeatureDefn& inLayerDefn = layer.GetLayerDefn();
-    for (int k = 0; k < inLayerDefn.GetFieldCount(); k++)
-    {
-      OGRFieldDefn fieldDefn(inLayerDefn.GetFieldDefn(k));
-      newLayer.CreateField(fieldDefn);
-    }
+    OGRFieldDefn fieldDefn(inLayerDefn.GetFieldDefn(k));
+    newLayer.CreateField(fieldDefn);
   }
-
   return output;
 }
 
-
 template <bool RegressionMode>
 void VectorPrediction<RegressionMode>::AddPredictionField(otb::ogr::Layer& outLayer, otb::ogr::Layer const& layer, bool computeConfidenceMap)
 {
@@ -306,7 +304,21 @@ void           VectorPrediction<RegressionMode>::DoExecute()
 
   const bool updateMode = !(IsParameterEnabled("out") && HasValue("out"));
 
-  auto            output   = CreateOutputDataSource(source, layer, updateMode);
+  ogr::DataSource::Pointer buffer;
+  ogr::DataSource::Pointer output;
+
+  if (updateMode)
+  {
+    // in update mode, output is added to input data source.
+    // buffer needs to be allocated here, as its life-cycle is bound to "layer"
+    buffer = ogr::DataSource::New();
+    output = ReopenDataSourceInUpdateMode(source, layer, buffer);
+  }
+  else
+  {
+    output = CreateOutputDataSource(source, layer);
+  }
+
   otb::ogr::Layer outLayer = output->GetLayer(0);
 
   OGRErr errStart = outLayer.ogr().StartTransaction();
diff --git a/Modules/Applications/AppImageUtils/test/CMakeLists.txt b/Modules/Applications/AppImageUtils/test/CMakeLists.txt
index 44eb49e7f9cbc2842df10dac085472d759b266bd..73a1ed86a3dbd831826175285f8f8ef08036b17b 100644
--- a/Modules/Applications/AppImageUtils/test/CMakeLists.txt
+++ b/Modules/Applications/AppImageUtils/test/CMakeLists.txt
@@ -143,7 +143,7 @@ otb_test_application(NAME apTvUtTileFusion
                              -rows 2
                              -out ${TEMP}/apTvUtTileFusion.png uint8
                      VALID   --compare-image ${NOTOL}
-                             ${INPUTDATA}/Scene.png
+                             ${INPUTDATA}/scene.png
                            ${TEMP}/apTvUtTileFusion.png)
 
 
@@ -203,6 +203,20 @@ otb_test_application(NAME apTvUtQuicklookROI1Channel
                              ${TEMP}/apTvUtQuicklookROI1Channel.tif
                      )
 
+otb_test_application(NAME apTvUtQuicklookWithGCP
+                     APP Quicklook
+                     OPTIONS -in ${INPUTDATA}/spot5SubWithGcps.tif
+                             -out ${TEMP}/apTvUtQuicklookWithGCP.tif uint8
+                             -sr 4
+                             -rox 100
+                             -roy 100
+                             -rsx 200
+                             -rsy 200
+                     VALID   --compare-metadata ${NOTOL}
+                             ${BASELINE}/apTvUtQuicklookWithGCP.tif
+                             ${TEMP}/apTvUtQuicklookWithGCP.tif
+                     )
+
 otb_test_application(NAME apTvUtQuicklookSpot5
                      APP Quicklook
                      OPTIONS -in LARGEINPUT{SPOT5/TEHERAN/IMAGERY.TIF}
diff --git a/Modules/Core/ImageBase/test/otbExtractROITestMetaData.cxx b/Modules/Core/ImageBase/test/otbExtractROITestMetaData.cxx
index b19c7e65840c4c9f95be8183b0d60ba7b61f5b0d..d3c3061185ea1bb243669e349bb5d75f2cde7d9c 100644
--- a/Modules/Core/ImageBase/test/otbExtractROITestMetaData.cxx
+++ b/Modules/Core/ImageBase/test/otbExtractROITestMetaData.cxx
@@ -30,8 +30,6 @@
 #include "otbImageFileReader.h"
 #include "otbImageFileWriter.h"
 
-#include "otb_boost_string_header.h"
-
 int otbExtractROITestMetaData(int itkNotUsed(argc), char* argv[])
 {
   typedef float PixelType;
@@ -100,17 +98,20 @@ int otbExtractROITestMetaData(int itkNotUsed(argc), char* argv[])
   reader00->SetFileName(argv[2]);
   reader00->GenerateOutputInformation();
 
-  if (reader00->GetOutput()->GetProjectionRef() != "" || boost::algorithm::istarts_with(reader00->GetOutput()->GetProjectionRef(), "LOCAL_CS"))
+  // The input image should have a sensor model and GCP, but the output images
+  // should only have the sensor model (priority over GCP). This behaviour
+  // must be consistent regardless of the ROI.
+
+  if (reader00->GetOutput()->GetProjectionRef().size())
   {
     std::cout << "The read generated extract from index (0, 0) must NOT contain a ProjectionReference." << std::endl;
     std::cout << "Found ProjectionReference: " << reader00->GetOutput()->GetProjectionRef() << std::endl;
-
     return EXIT_FAILURE;
   }
 
-  if (reader00->GetOutput()->GetGCPCount() == 0)
+  if (reader00->GetOutput()->GetGCPCount())
   {
-    std::cout << "The read generated extract from index (0, 0) must contain a list a GCPs.." << std::endl;
+    std::cout << "The read generated extract from index (0, 0) must NOT contain a list a GCPs.." << std::endl;
     return EXIT_FAILURE;
   }
 
@@ -118,14 +119,14 @@ int otbExtractROITestMetaData(int itkNotUsed(argc), char* argv[])
   reader57->SetFileName(argv[3]);
   reader57->GenerateOutputInformation();
 
-  if (reader57->GetOutput()->GetProjectionRef() != "" || boost::algorithm::istarts_with(reader57->GetOutput()->GetProjectionRef(), "LOCAL_CS"))
+  if (reader57->GetOutput()->GetProjectionRef().size())
   {
     std::cout << "The read generated extract from index (x, y) must NOT contain a ProjectionReference." << std::endl;
     std::cout << "Found ProjectionReference: " << reader57->GetOutput()->GetProjectionRef() << std::endl;
     return EXIT_FAILURE;
   }
 
-  if (reader57->GetOutput()->GetGCPCount() != 0)
+  if (reader57->GetOutput()->GetGCPCount())
   {
     std::cout << "The read generated extract from index (x, y) must NOT contain a list a GCPs.." << std::endl;
     return EXIT_FAILURE;
diff --git a/Modules/Filtering/MathParserX/src/otbParserXPlugins.cxx b/Modules/Filtering/MathParserX/src/otbParserXPlugins.cxx
index 9784759177908412d6b611b7f5ae8ee84252ad23..203f34b18dc2634b36eee50bd3e1ad97d08b8542 100644
--- a/Modules/Filtering/MathParserX/src/otbParserXPlugins.cxx
+++ b/Modules/Filtering/MathParserX/src/otbParserXPlugins.cxx
@@ -34,9 +34,7 @@ namespace otb
 void bands::Eval(mup::ptr_val_type& ret, const mup::ptr_val_type* a_pArg, int a_iArgc)
 {
 
-  if (a_iArgc != 2)
-    return;
-
+  assert (a_iArgc == 2);
   assert(a_pArg[0]->GetType() == 'm');
   assert(a_pArg[1]->GetType() == 'm');
 
diff --git a/Modules/Filtering/MathParserX/test/otbBandMathXImageFilter.cxx b/Modules/Filtering/MathParserX/test/otbBandMathXImageFilter.cxx
index 2c076a902dd5689d22429ae95cf00fda5e399426..df56d155674e609b4107b521383dafb253dc5cc2 100644
--- a/Modules/Filtering/MathParserX/test/otbBandMathXImageFilter.cxx
+++ b/Modules/Filtering/MathParserX/test/otbBandMathXImageFilter.cxx
@@ -31,6 +31,23 @@
 #include "itkImageRegionIteratorWithIndex.h"
 
 
+double relativeError(double a, double b)
+{
+  return std::pow(a-b, 2) / std::fabs(a+b);
+}
+
+template<typename T> typename T::Pointer createTestImage(typename T::RegionType& region,
+                                                         unsigned int components)
+{
+  typename T::Pointer image = T::New();
+  image->SetLargestPossibleRegion(region);
+  image->SetBufferedRegion(region);
+  image->SetRequestedRegion(region);
+  image->SetNumberOfComponentsPerPixel(components);
+  image->Allocate();
+  return image;
+}
+
 int otbBandMathXImageFilter(int itkNotUsed(argc), char* itkNotUsed(argv)[])
 {
 
@@ -39,7 +56,7 @@ int otbBandMathXImageFilter(int itkNotUsed(argc), char* itkNotUsed(argv)[])
   typedef otb::BandMathXImageFilter<ImageType> FilterType;
 
   unsigned int       i;
-  const unsigned int N = 100, D1 = 3, D2 = 1, D3 = 1;
+  const unsigned int N = 100, D1 = 3, D2 = 1, D3 = 1, D4=4;
 
   ImageType::SizeType size;
   size.Fill(N);
@@ -49,39 +66,29 @@ int otbBandMathXImageFilter(int itkNotUsed(argc), char* itkNotUsed(argv)[])
   region.SetSize(size);
   region.SetIndex(index);
 
-  ImageType::Pointer image1 = ImageType::New();
-  ImageType::Pointer image2 = ImageType::New();
-  ImageType::Pointer image3 = ImageType::New();
-
-  image1->SetLargestPossibleRegion(region);
-  image1->SetBufferedRegion(region);
-  image1->SetRequestedRegion(region);
-  image1->SetNumberOfComponentsPerPixel(D1);
-  image1->Allocate();
-
-  image2->SetLargestPossibleRegion(region);
-  image2->SetBufferedRegion(region);
-  image2->SetRequestedRegion(region);
-  image2->SetNumberOfComponentsPerPixel(D2);
-  image2->Allocate();
-
-  image3->SetLargestPossibleRegion(region);
-  image3->SetBufferedRegion(region);
-  image3->SetRequestedRegion(region);
-  image3->SetNumberOfComponentsPerPixel(D3);
-  image3->Allocate();
+  ImageType::Pointer image1 = createTestImage<ImageType>(region, D1);
+  ImageType::Pointer image2 = createTestImage<ImageType>(region, D2);
+  ImageType::Pointer image3 = createTestImage<ImageType>(region, D3);
+  ImageType::Pointer image4 = createTestImage<ImageType>(region, D2);
+  ImageType::Pointer image5 = createTestImage<ImageType>(region, D4);
 
   typedef itk::ImageRegionIteratorWithIndex<ImageType> IteratorType;
   IteratorType                                         it1(image1, region);
   IteratorType                                         it2(image2, region);
   IteratorType                                         it3(image3, region);
+  IteratorType                                         it4(image4, region);
+  IteratorType                                         it5(image5, region);
 
-  ImageType::PixelType val1, val2, val3;
+  ImageType::PixelType val1, val2, val3, val4, val5;
   val1.SetSize(D1);
   val2.SetSize(D2);
   val3.SetSize(D3);
+  val4.SetSize(D2);
+  val5.SetSize(D4);
 
-  for (it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++it3)
+  for (i=0, it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(), it4.GoToBegin(), it5.GoToBegin();
+       !it1.IsAtEnd();
+       ++i, ++it1, ++it2, ++it3, ++it4, ++it5)
   {
     ImageType::IndexType i1 = it1.GetIndex();
     ImageType::IndexType i2 = it2.GetIndex();
@@ -92,10 +99,18 @@ int otbBandMathXImageFilter(int itkNotUsed(argc), char* itkNotUsed(argv)[])
     val1[2] = i1[0] / (i1[1] + 1) + 5;
     val2[0] = i2[0] * i2[1];
     val3[0] = i3[0] + i3[1] * i3[1];
+    val4[0] = 2e-4 * i;
+
+    val5[0] = i;
+    val5[1] = i*2;
+    val5[2] = i*3;
+    val5[3] = i*4;
 
     it1.Set(val1);
     it2.Set(val2);
     it3.Set(val3);
+    it4.Set(val4);
+    it5.Set(val5);
   }
 
 
@@ -105,15 +120,41 @@ int otbBandMathXImageFilter(int itkNotUsed(argc), char* itkNotUsed(argv)[])
   filter->SetNthInput(0, image1);
   filter->SetNthInput(1, image2);
   filter->SetNthInput(2, image3, "canal3");
+  filter->SetNthInput(3, image4);
+  filter->SetNthInput(4, image5);
+
+  filter->SetExpression( // Sub-test 1
+      "vcos(2 * pi * im1) div (2 * pi * (bands(im2,{1,1,1}) pw 2)"
+      "+ {3.38,3.38,3.38} pow {2.1,3.2,4}) mult (vsin(pi * bands(canal3,{1,1,1})) dv 2.0)");
 
-  filter->SetExpression("vcos(2 * pi * im1) div (2 * pi * bands(im2,{1,1,1}) + {3.38,3.38,3.38}) mult vsin(pi * bands(canal3,{1,1,1}))"); // Sub-test 1
   filter->SetExpression("im1b1 / im2b1"); // Sub-test 2 (Edge Effect Handling)
+
+  // Sub-test 3 : MultiplicationByScalar, ElementWisePower, DivisionByScalar
+  filter->SetExpression("((im2 mlt 3) pow (im4 mlt {2.0})) dv {2.0}");
+
+  // Sub-test 4 : PowerByScalar, vtan, vtanh, vsinh, vcosh
+  filter->SetExpression("vcosh(vsinh(vtanh(vtan(im2 pw 2 + im4 pw {2.2}))))");
+
+  // Sub-test 5 : vlog, vlog10, vabs, vexp, vsqrt
+  filter->SetExpression("vlog(vabs(vsqrt(vexp(vlog10(im2 + bands(im1,{3}))))))");
+
+  // Sub-test 6 : mean
+  filter->SetExpression("mean(10, 3.14)");
+
+  // Sub-test 7 : corr
+  filter->SetExpression("corr(bands(im5, {1,2,3,4}), bands(im5, {1,2,4,3}))");
+
   filter->Update();
 
   // if (filter->GetNumberOfOutputs() != 2)
 
   ImageType::Pointer output1 = filter->GetOutput(0);
   ImageType::Pointer output2 = filter->GetOutput(1);
+  ImageType::Pointer output3 = filter->GetOutput(2);
+  ImageType::Pointer output4 = filter->GetOutput(3);
+  ImageType::Pointer output5 = filter->GetOutput(4);
+  ImageType::Pointer output6 = filter->GetOutput(5);
+  ImageType::Pointer output7 = filter->GetOutput(6);
 
   std::cout << "\n---  Standard Use\n";
   std::cout << "Parsed Expression :   " << filter->GetExpression(0) << std::endl;
@@ -138,28 +179,28 @@ int otbBandMathXImageFilter(int itkNotUsed(argc), char* itkNotUsed(argv)[])
     double error1, error2, error3;
 
 
-    double expected1 = std::cos(2 * otb::CONST_PI * px1[0]) / (2 * otb::CONST_PI * px2[0] + 3.38) * std::sin(otb::CONST_PI * px3[0]);
-    double expected2 = std::cos(2 * otb::CONST_PI * px1[1]) / (2 * otb::CONST_PI * px2[0] + 3.38) * std::sin(otb::CONST_PI * px3[0]);
-    double expected3 = std::cos(2 * otb::CONST_PI * px1[2]) / (2 * otb::CONST_PI * px2[0] + 3.38) * std::sin(otb::CONST_PI * px3[0]);
+    double expected1 = std::cos(2 * otb::CONST_PI * px1[0]) / (2 * otb::CONST_PI * px2[0] * px2[0] + pow(3.38,2.1)) * std::sin(otb::CONST_PI * px3[0]) / 2;
+    double expected2 = std::cos(2 * otb::CONST_PI * px1[1]) / (2 * otb::CONST_PI * px2[0] * px2[0] + pow(3.38,3.2)) * std::sin(otb::CONST_PI * px3[0]) / 2;
+    double expected3 = std::cos(2 * otb::CONST_PI * px1[2]) / (2 * otb::CONST_PI * px2[0] * px2[0] + pow(3.38,4)) * std::sin(otb::CONST_PI * px3[0]) / 2;
 
     /*std::cout << "Pixel_1 =  " << it1.Get()[0] << "     Pixel_2 =  " << it2.Get()[0] << "     Pixel_3 =  " << it3.Get()[0]
         << "     Result =  " << itoutput1.Get()[0] << "     Expected =  " << expected1 << std::endl; */
 
 
-    error1 = (result1 - expected1) * (result1 - expected1) / abs(result1 + expected1);
-    error2 = (result2 - expected2) * (result2 - expected2) / abs(result2 + expected2);
-    error3 = (result3 - expected3) * (result3 - expected3) / abs(result3 + expected3);
+    error1 = relativeError(result1, expected1);
+    error2 = relativeError(result2, expected2);
+    error3 = relativeError(result3, expected3);
 
     if ((error1 > 1E-9) || (error2 > 1E-9) || (error3 > 1E-9))
     {
       itkGenericExceptionMacro(<< std::endl
-                               << "Error = " << error1 << "  > 1E-9     -> TEST FAILLED" << std::endl
+                               << "Error = " << error1 << "  > 1E-9     -> TEST FAILED" << std::endl
                                << "Pixel_1 =  " << it1.Get()[0] << "     Pixel_2 =  " << it2.Get()[0] << "     Pixel_3 =  " << it3.Get()[0]
                                << "     Result =  " << result1 << "     Expected =  " << expected1 << std::endl
-                               << "Error = " << error2 << "  > 1E-9     -> TEST FAILLED" << std::endl
+                               << "Error = " << error2 << "  > 1E-9     -> TEST FAILED" << std::endl
                                << "Pixel_1 =  " << it1.Get()[1] << "     Pixel_2 =  " << it2.Get()[0] << "     Pixel_3 =  " << it3.Get()[0]
                                << "     Result =  " << result2 << "     Expected =  " << expected2 << std::endl
-                               << "Error = " << error3 << "  > 1E-9     -> TEST FAILLED" << std::endl
+                               << "Error = " << error3 << "  > 1E-9     -> TEST FAILED" << std::endl
                                << "Pixel_1 =  " << it1.Get()[2] << "     Pixel_2 =  " << it2.Get()[0] << "     Pixel_3 =  " << it3.Get()[0]
                                << "     Result =  " << result3 << "     Expected =  " << expected3 << std::endl);
     }
@@ -187,6 +228,96 @@ int otbBandMathXImageFilter(int itkNotUsed(argc), char* itkNotUsed(argv)[])
                              << "     Expected =  nan\n");
   std::cout << std::endl;
 
+  // Sub-test 3
+
+  IteratorType itoutput3(output3, region);
+  for (i=0, it2.GoToBegin(); !it2.IsAtEnd(); ++i, ++it2, ++itoutput3)
+  {
+    ImageType::IndexType i2 = it2.GetIndex();
+
+    double expected = pow(i2[0] * i2[1] * 3, 2 * 2e-4 * i) / 2.0;
+    double result = itoutput3.Get()[0];
+    double error;
+
+    error = std::fabs(result - expected);
+
+    if (error > 1e-9) {
+      itkGenericExceptionMacro(<< std::endl
+                               << "Result = " << itoutput3.Get()[0]
+                               << " Expected = " << expected
+                               << " Error = " << error << "  > 1E-9     -> TEST FAILED"
+                               << std::endl);
+    }
+  }
+
+  // Sub-test 4
+
+  IteratorType itoutput4(output4, region);
+  for (i=0, it2.GoToBegin(); !it2.IsAtEnd(); ++i, ++it2, ++itoutput4)
+  {
+    ImageType::IndexType i2 = it2.GetIndex();
+
+    double expected = std::cosh(std::sinh(std::tanh(std::tan(i2[0] * i2[1] * i2[0] * i2[1]
+                                                             + pow(2e-4 * i, 2.2)))));
+    double result = itoutput4.Get()[0];
+    double error;
+
+    error = std::fabs(result - expected);
+
+    if (error > 1e-12) {
+      itkGenericExceptionMacro(<< std::endl
+                               << "Result = " << itoutput4.Get()[0]
+                               << " Expected = " << expected
+                               << " Error = " << error << "  > 1E-12     -> TEST FAILED"
+                               << std::endl);
+    }
+  }
+
+  // Sub-test 5
+
+  IteratorType itoutput5(output5, region);
+  for (i=0, it1.GoToBegin(), it2.GoToBegin(); !it2.IsAtEnd(); ++i, ++it1, ++it2, ++itoutput5)
+  {
+    ImageType::IndexType i2 = it2.GetIndex();
+
+    double expected = std::log(std::fabs(std::sqrt(std::exp(std::log10(i2[0] * i2[1] + it1.Get()[2])))));
+    double result = itoutput5.Get()[0];
+    double error;
+
+    error = std::fabs(result - expected);
+
+    if (error > 1e-12) {
+      itkGenericExceptionMacro(<< std::endl
+                               << "Result = " << itoutput5.Get()[0]
+                               << " Expected = " << expected
+                               << " Error = " << error << "  > 1E-12     -> TEST FAILED"
+                               << std::endl);
+    }
+  }
+
+  // Sub-test 6
+  PixelType result = output6->GetPixel( {0, 0} );
+  if (result[0] != 10 || std::abs(result[1] - 3.14) > 1e-12) {
+    itkGenericExceptionMacro(<< std::endl
+                             << " Expected = [10, 3.14]"
+                             << " Got = " << result << "     -> TEST FAILED"
+                             << std::endl);
+  }
+
+  // Sub-test 7
+  IteratorType itoutput7(output7, region);
+  for (i=0, itoutput7.GoToBegin(); !itoutput7.IsAtEnd(); ++i, ++itoutput7)
+  {
+    double result = itoutput7.Get()[0];
+    double error = std::abs(result - 0.8);
+
+    if (error > 1e-12) {
+      itkGenericExceptionMacro(<< std::endl
+                               << " Expected = 0.8"
+                               << " Got = " << result << "     -> TEST FAILED"
+                               << std::endl);
+    }
+  }
 
   return EXIT_SUCCESS;
 }
@@ -210,27 +341,9 @@ int otbBandMathXImageFilterConv(int itkNotUsed(argc), char* argv[])
   region.SetSize(size);
   region.SetIndex(index);
 
-  ImageType::Pointer image1 = ImageType::New();
-  ImageType::Pointer image2 = ImageType::New();
-  ImageType::Pointer image3 = ImageType::New();
-
-  image1->SetLargestPossibleRegion(region);
-  image1->SetBufferedRegion(region);
-  image1->SetRequestedRegion(region);
-  image1->SetNumberOfComponentsPerPixel(D1);
-  image1->Allocate();
-
-  image2->SetLargestPossibleRegion(region);
-  image2->SetBufferedRegion(region);
-  image2->SetRequestedRegion(region);
-  image2->SetNumberOfComponentsPerPixel(D2);
-  image2->Allocate();
-
-  image3->SetLargestPossibleRegion(region);
-  image3->SetBufferedRegion(region);
-  image3->SetRequestedRegion(region);
-  image3->SetNumberOfComponentsPerPixel(D3);
-  image3->Allocate();
+  ImageType::Pointer image1 = createTestImage<ImageType>(region, D1);
+  ImageType::Pointer image2 = createTestImage<ImageType>(region, D2);
+  ImageType::Pointer image3 = createTestImage<ImageType>(region, D3);
 
   typedef itk::ConstNeighborhoodIterator<ImageType> IteratorType;
   IteratorType::RadiusType                          radius;
@@ -386,13 +499,13 @@ int otbBandMathXImageFilterConv(int itkNotUsed(argc), char* argv[])
 
     double expected1 = px1[0], expected2 = px1[1], expected3 = px1[2], expected4 = px1[3], expected5 = px1[4], expected6 = px1[5], expected7 = px1[6];
 
-    error1 = (result1 - expected1) * (result1 - expected1) / (result1 + expected1);
-    error2 = (result2 - expected2) * (result2 - expected2) / (result2 + expected2);
-    error3 = (result3 - expected3) * (result3 - expected3) / (result3 + expected3);
-    error4 = (result4 - expected4) * (result4 - expected4) / (result4 + expected4);
-    error5 = (result5 - expected5) * (result5 - expected5) / (result5 + expected5);
-    error6 = (result6 - expected6) * (result6 - expected6) / (result6 + expected6);
-    error7 = (result7 - expected7) * (result7 - expected7) / (result7 + expected7);
+    error1 = relativeError(result1, expected1);
+    error2 = relativeError(result2, expected2);
+    error3 = relativeError(result3, expected3);
+    error4 = relativeError(result4, expected4);
+    error5 = relativeError(result5, expected5);
+    error6 = relativeError(result6, expected6);
+    error7 = relativeError(result7, expected7);
 
 
     // expression 2
@@ -400,13 +513,13 @@ int otbBandMathXImageFilterConv(int itkNotUsed(argc), char* argv[])
     double result9   = itoutput2.GetCenterPixel()[1];
     double expected8 = px2[0];
     double expected9 = px2[1];
-    double error8    = (result8 - expected8) * (result8 - expected8) / (result8 + expected8);
-    double error9    = (result9 - expected9) * (result9 - expected9) / (result9 + expected9);
+    double error8    = relativeError(result8, expected8);
+    double error9    = relativeError(result9, expected9);
 
     if ((error1 > 1E-9) || (error2 > 1E-9) || (error3 > 1E-9) || (error4 > 1E-9) || (error5 > 1E-9) || (error6 > 1E-9) || (error7 > 1E-9) || (error8 > 1E-9) ||
         (error9 > 1E-9))
     {
-      itkGenericExceptionMacro(<< "TEST FAILLED" << std::endl
+      itkGenericExceptionMacro(<< "TEST FAILED" << std::endl
                                << "Error1 = " << error1 << std::endl
                                << "     Result1 =  " << result1 << "     Expected1 =  " << expected1 << std::endl
                                << "Error2 = " << error2 << std::endl
@@ -475,27 +588,9 @@ int otbBandMathXImageFilterWithIdx(int itkNotUsed(argc), char* argv[])
   spacing[0] = 0.5;
   spacing[1] = 0.5;
 
-  ImageType::Pointer image1 = ImageType::New();
-  ImageType::Pointer image2 = ImageType::New();
-  ImageType::Pointer image3 = ImageType::New();
-
-  image1->SetLargestPossibleRegion(region);
-  image1->SetBufferedRegion(region);
-  image1->SetRequestedRegion(region);
-  image1->SetNumberOfComponentsPerPixel(D1);
-  image1->Allocate();
-
-  image2->SetLargestPossibleRegion(region);
-  image2->SetBufferedRegion(region);
-  image2->SetRequestedRegion(region);
-  image2->SetNumberOfComponentsPerPixel(D2);
-  image2->Allocate();
-
-  image3->SetLargestPossibleRegion(region);
-  image3->SetBufferedRegion(region);
-  image3->SetRequestedRegion(region);
-  image3->SetNumberOfComponentsPerPixel(D3);
-  image3->Allocate();
+  ImageType::Pointer image1 = createTestImage<ImageType>(region, D1);
+  ImageType::Pointer image2 = createTestImage<ImageType>(region, D2);
+  ImageType::Pointer image3 = createTestImage<ImageType>(region, D3);
 
   typedef itk::ImageRegionIteratorWithIndex<ImageType> IteratorType;
   IteratorType                                         it1(image1, region);
@@ -520,11 +615,9 @@ int otbBandMathXImageFilterWithIdx(int itkNotUsed(argc), char* argv[])
     it3.Get()[0] = i3[0] + i3[1] * i3[1];
   }
 
-
   FilterType::Pointer filter = FilterType::New();
   std::cout << "Number Of Threads  :  " << filter->GetNumberOfThreads() << std::endl;
 
-
   filter->SetNthInput(0, image1);
   filter->SetNthInput(1, image2);
   filter->SetNthInput(2, image3);
@@ -561,13 +654,7 @@ int otbBandMathXImageFilterBandsFailures(int itkNotUsed(argc), char* itkNotUsed(
   region.SetSize(size);
   region.SetIndex(index);
 
-  ImageType::Pointer image = ImageType::New();
-
-  image->SetLargestPossibleRegion(region);
-  image->SetBufferedRegion(region);
-  image->SetRequestedRegion(region);
-  image->SetNumberOfComponentsPerPixel(D1);
-  image->Allocate();
+  ImageType::Pointer image = createTestImage<ImageType>(region, D1);
 
   FilterType::Pointer filter = FilterType::New();
 
@@ -586,10 +673,14 @@ int otbBandMathXImageFilterBandsFailures(int itkNotUsed(argc), char* itkNotUsed(
       std::cout << "INFO: Exception thrown as expected : " << e.what() << std::endl;
       filter->ClearExpression();
     }
-    catch (...) {
-      itkGenericExceptionMacro(<< "TEST FAILLED: " << exp
-                               << "should have raise a RangeError exception" << std::endl);
+    catch (const std::exception& e) {
+      itkGenericExceptionMacro(<< "TEST FAILED: " << exp
+                               << "should have raise a RangeError exception."
+                               << e.what() << std::endl);
     }
+
+    filter->ClearExpression(); // clean-up filter for next test string
+
   }
   return EXIT_SUCCESS;
 }
diff --git a/Modules/IO/IOGDAL/src/otbGDALImageIO.cxx b/Modules/IO/IOGDAL/src/otbGDALImageIO.cxx
index aca9771e01c53c6aedf3da4201ef91058ee13584..165b4549ecf9aabd60e1e7484d7a43414123a6f8 100644
--- a/Modules/IO/IOGDAL/src/otbGDALImageIO.cxx
+++ b/Modules/IO/IOGDAL/src/otbGDALImageIO.cxx
@@ -743,14 +743,18 @@ void GDALImageIO::InternalReadImageInformation()
   /* -------------------------------------------------------------------- */
   /* Get the projection coordinate system of the image : ProjectionRef  */
   /* -------------------------------------------------------------------- */
-  if (dataset->GetProjectionRef() != nullptr && !std::string(dataset->GetProjectionRef()).empty())
+  const char* pszProjection = dataset->GetProjectionRef();
+  if (pszProjection != nullptr && !std::string(pszProjection).empty())
   {
     OGRSpatialReferenceH pSR = OSRNewSpatialReference(nullptr);
-
-    const char* pszProjection = nullptr;
-    pszProjection             = dataset->GetProjectionRef();
-
-    if (OSRImportFromWkt(pSR, (char**)(&pszProjection)) == OGRERR_NONE)
+    if (strncmp(pszProjection, "LOCAL_CS",8) == 0)
+      {
+      // skip local coordinate system as they will cause crashed later
+      // In GDAL 3, they begin to do special processing for Transmercator local
+      // coordinate system
+      otbLogMacro(Debug, << "Skipping LOCAL_CS projection")
+      }
+    else if (OSRImportFromWkt(pSR, (char**)(&pszProjection)) == OGRERR_NONE)
     {
       char* pszPrettyWkt = nullptr;
       OSRExportToPrettyWkt(pSR, &pszPrettyWkt, FALSE);
@@ -761,7 +765,7 @@ void GDALImageIO::InternalReadImageInformation()
     }
     else
     {
-      itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey, static_cast<std::string>(dataset->GetProjectionRef()));
+      itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey, static_cast<std::string>(pszProjection));
     }
 
     if (pSR != nullptr)
@@ -1399,111 +1403,120 @@ void GDALImageIO::InternalWriteImageInformation(const void* buffer)
 
   std::string projectionRef;
   itk::ExposeMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey, projectionRef);
+  ImageKeywordlist otb_kwl;
+  itk::ExposeMetaData<ImageKeywordlist>(dict, MetaDataKey::OSSIMKeywordlistKey, otb_kwl);
+  unsigned int gcpCount = 0;
+  itk::ExposeMetaData<unsigned int>(dict, MetaDataKey::GCPCountKey, gcpCount);
+
+  // In OTB we can have simultaneously projection ref, sensor keywordlist, GCPs
+  // but GDAL can't handle both geotransform and GCP (see issue #303). Here is the priority
+  // order we will be using in OTB:
+  // [ProjRef+geotransform] > [SensorKeywordlist+geotransform] > [GCPs]
 
   /* -------------------------------------------------------------------- */
-  /* Set the GCPs                                                          */
+  /* Pre-compute geotransform                                             */
   /* -------------------------------------------------------------------- */
   const double Epsilon = 1E-10;
-  if (projectionRef.empty() && (std::abs(m_Origin[0] - 0.5) > Epsilon || std::abs(m_Origin[1] - 0.5) > Epsilon ||
-                                std::abs(m_Spacing[0] * m_Direction[0][0] - 1.0) > Epsilon || std::abs(m_Spacing[1] * m_Direction[1][1] - 1.0) > Epsilon))
-  {
-    // See issue #303 :
-    // If there is no ProjectionRef, and the GeoTransform is not the identity,
-    // then saving also GCPs is undefined behavior for GDAL, and a WGS84 projection crs
-    // is assigned arbitrarily
-    otbLogMacro(Warning, << "Skipping GCPs saving to prevent GDAL from assigning a WGS84 projref to file (" << m_FileName << ")")
-  }
-  else
-  {
-    unsigned int gcpCount = 0;
-    itk::ExposeMetaData<unsigned int>(dict, MetaDataKey::GCPCountKey, gcpCount);
-
-    if (gcpCount > 0)
-    {
-
-      GDAL_GCP* gdalGcps = new GDAL_GCP[gcpCount];
-
-      for (unsigned int gcpIndex = 0; gcpIndex < gcpCount; ++gcpIndex)
-      {
-        // Build the GCP string in the form of GCP_n
-        std::ostringstream lStream;
-        lStream << MetaDataKey::GCPParametersKey << gcpIndex;
-        std::string key = lStream.str();
-
-        OTB_GCP gcp;
-        itk::ExposeMetaData<OTB_GCP>(dict, key, gcp);
-
-        gdalGcps[gcpIndex].pszId      = const_cast<char*>(gcp.m_Id.c_str());
-        gdalGcps[gcpIndex].pszInfo    = const_cast<char*>(gcp.m_Info.c_str());
-        gdalGcps[gcpIndex].dfGCPPixel = gcp.m_GCPCol;
-        gdalGcps[gcpIndex].dfGCPLine  = gcp.m_GCPRow;
-        gdalGcps[gcpIndex].dfGCPX     = gcp.m_GCPX;
-        gdalGcps[gcpIndex].dfGCPY     = gcp.m_GCPY;
-        gdalGcps[gcpIndex].dfGCPZ     = gcp.m_GCPZ;
-      }
-
-      std::string gcpProjectionRef;
-      itk::ExposeMetaData<std::string>(dict, MetaDataKey::GCPProjectionKey, gcpProjectionRef);
-
-      dataset->SetGCPs(gcpCount, gdalGcps, gcpProjectionRef.c_str());
-
-      delete[] gdalGcps;
-    }
-  }
+  double geoTransform[6];
+  geoTransform[0] = m_Origin[0] - 0.5 * m_Spacing[0] * m_Direction[0][0];
+  geoTransform[3] = m_Origin[1] - 0.5 * m_Spacing[1] * m_Direction[1][1];
+  geoTransform[1] = m_Spacing[0] * m_Direction[0][0];
+  geoTransform[5] = m_Spacing[1] * m_Direction[1][1];
+  // FIXME: Here component 1 and 4 should be replaced by the orientation parameters
+  geoTransform[2] = 0.;
+  geoTransform[4] = 0.;
+  // only write geotransform if it has non-default values
+  bool writeGeotransform =
+    std::abs(geoTransform[0]) > Epsilon ||
+    std::abs(geoTransform[1] - 1.0) > Epsilon ||
+    std::abs(geoTransform[2]) > Epsilon ||
+    std::abs(geoTransform[3]) > Epsilon ||
+    std::abs(geoTransform[4]) > Epsilon ||
+    std::abs(geoTransform[5] - 1.0) > Epsilon;
 
   /* -------------------------------------------------------------------- */
-  /* Set the projection coordinate system of the image : ProjectionRef    */
+  /* Case 1: Set the projection coordinate system of the image            */
   /* -------------------------------------------------------------------- */
   if (!projectionRef.empty())
-  {
+    {
     dataset->SetProjection(projectionRef.c_str());
-  }
-  else
-  {
+    }
+  /* -------------------------------------------------------------------- */
+  /* Case 2: Sensor keywordlist                                           */
+  /* -------------------------------------------------------------------- */
+  else if (otb_kwl.GetSize())
+    {
     /* -------------------------------------------------------------------- */
-    /* Set the RPC coeffs if no projection available (since GDAL 1.10.0)    */
+    /* Set the RPC coeffs (since GDAL 1.10.0)                               */
     /* -------------------------------------------------------------------- */
-    ImageKeywordlist otb_kwl;
-    itk::ExposeMetaData<ImageKeywordlist>(dict, MetaDataKey::OSSIMKeywordlistKey, otb_kwl);
-    if (m_WriteRPCTags && otb_kwl.GetSize() != 0)
-    {
+    if (m_WriteRPCTags)
+      {
       GDALRPCInfo gdalRpcStruct;
       if (otb_kwl.convertToGDALRPC(gdalRpcStruct))
-      {
+        {
+        otbLogMacro(Debug, << "Saving RPC to file (" << m_FileName << ")")
         char** rpcMetadata = RPCInfoToMD(&gdalRpcStruct);
         dataset->SetMetadata(rpcMetadata, "RPC");
         CSLDestroy(rpcMetadata);
+        }
       }
     }
-  }
+  /* -------------------------------------------------------------------- */
+  /* Case 3: Set the GCPs                                                 */
+  /* -------------------------------------------------------------------- */
+  else if (gcpCount)
+    {
+    otbLogMacro(Debug, << "Saving GCPs to file (" << m_FileName << ")")
+    std::vector<GDAL_GCP> gdalGcps(gcpCount);
+    std::vector<OTB_GCP> otbGcps(gcpCount);
+    for (unsigned int gcpIndex = 0; gcpIndex < gcpCount; ++gcpIndex)
+      {
+      // Build the GCP string in the form of GCP_n
+      std::ostringstream lStream;
+      lStream << MetaDataKey::GCPParametersKey << gcpIndex;
+      std::string key = lStream.str();
+
+      OTB_GCP &gcp = otbGcps[gcpIndex];
+      itk::ExposeMetaData<OTB_GCP>(dict, key, gcp);
+
+      gdalGcps[gcpIndex].pszId      = const_cast<char*>(gcp.m_Id.c_str());
+      gdalGcps[gcpIndex].pszInfo    = const_cast<char*>(gcp.m_Info.c_str());
+      gdalGcps[gcpIndex].dfGCPPixel = gcp.m_GCPCol;
+      gdalGcps[gcpIndex].dfGCPLine  = gcp.m_GCPRow;
+      gdalGcps[gcpIndex].dfGCPX     = gcp.m_GCPX;
+      gdalGcps[gcpIndex].dfGCPY     = gcp.m_GCPY;
+      gdalGcps[gcpIndex].dfGCPZ     = gcp.m_GCPZ;
+
+      if (writeGeotransform)
+        {
+        // we need to transform GCP col and row accordingly
+        // WARNING: assume no rotation is there
+        gdalGcps[gcpIndex].dfGCPPixel = (gcp.m_GCPCol - geoTransform[0]) / geoTransform[1];
+        gdalGcps[gcpIndex].dfGCPLine  = (gcp.m_GCPRow - geoTransform[3]) / geoTransform[5];
+        }
+      }
+
+    std::string gcpProjectionRef;
+    itk::ExposeMetaData<std::string>(dict, MetaDataKey::GCPProjectionKey, gcpProjectionRef);
+
+    dataset->SetGCPs(gcpCount, gdalGcps.data(), gcpProjectionRef.c_str());
+
+    // disable geotransform with GCP
+    writeGeotransform = false;
+    }
 
   /* -------------------------------------------------------------------- */
-  /*  Set the six coefficients of affine geoTransform                     */
+  /*  Save geotransform if needed.                                        */
   /* -------------------------------------------------------------------- */
-  if (std::abs(m_Origin[0] - 0.5) > Epsilon || std::abs(m_Origin[1] - 0.5) > Epsilon || std::abs(m_Spacing[0] * m_Direction[0][0] - 1.0) > Epsilon ||
-      std::abs(m_Spacing[1] * m_Direction[1][1] - 1.0) > Epsilon)
-  {
-    // Only set the geotransform if it is not identity (it may erase GCP)
-    itk::VariableLengthVector<double> geoTransform(6);
-    /// Reporting origin and spacing
-    // Beware : GDAL origin is at the corner of the top-left pixel
-    // whereas OTB/ITK origin is at the centre of the top-left pixel
-    geoTransform[0] = m_Origin[0] - 0.5 * m_Spacing[0] * m_Direction[0][0];
-    geoTransform[3] = m_Origin[1] - 0.5 * m_Spacing[1] * m_Direction[1][1];
-    geoTransform[1] = m_Spacing[0] * m_Direction[0][0];
-    geoTransform[5] = m_Spacing[1] * m_Direction[1][1];
-
-    // FIXME: Here component 1 and 4 should be replaced by the orientation parameters
-    geoTransform[2] = 0.;
-    geoTransform[4] = 0.;
-    dataset->SetGeoTransform(const_cast<double*>(geoTransform.GetDataPointer()));
-
-    // Error if writing to a ENVI file with a positive Y spacing
-    if (driverShortName == "ENVI" && geoTransform[5] > 0.)
+  if (writeGeotransform)
     {
+    if ( driverShortName == "ENVI" && geoTransform[5] > 0.)
+      {
+      // Error if writing to a ENVI file with a positive Y spacing
       itkExceptionMacro(<< "Can not write to ENVI file format with a positive Y spacing (" << m_FileName << ")");
+      }
+    dataset->SetGeoTransform(geoTransform);
     }
-  }
 
   /* -------------------------------------------------------------------- */
   /*      Report metadata.                                                */
diff --git a/Modules/Wrappers/ApplicationEngine/include/otbWrapperApplication.h b/Modules/Wrappers/ApplicationEngine/include/otbWrapperApplication.h
index 726055b9cfcc206182e783e945165534efa2751d..045a8352e42a85d5fe228f394244e44392645e82 100644
--- a/Modules/Wrappers/ApplicationEngine/include/otbWrapperApplication.h
+++ b/Modules/Wrappers/ApplicationEngine/include/otbWrapperApplication.h
@@ -219,6 +219,7 @@ public:
    * \li ParameterType_Int
    * \li ParameterType_Bool
    * \li ParameterType_Float
+   * \li ParameterType_Double
    * \li ParameterType_Radius
    * \li ParameterType_Choice
    */
@@ -228,9 +229,18 @@ public:
    *
    * Can be called for types :
    * \li ParameterType_Float
+   * \li ParameterType_Double
    */
   void SetParameterFloat(std::string const& parameter, float value, bool hasUserValueFlag = true);
 
+  /* Set a double precision floating value
+   *
+   * Can be called for types :
+   * \li ParameterType_Float
+   * \li ParameterType_Double
+   */
+  void SetParameterDouble(std::string const& parameter, double value, bool hasUserValueFlag = true);
+
   /* Set a string value
    *
    * Can be called for types :
@@ -253,6 +263,7 @@ public:
    * \li ParameterType_Directory
    * \li ParameterType_Choice
    * \li ParameterType_Float
+   * \li ParameterType_Double
    * \li ParameterType_Int
    * \li ParameterType_Radius
    * \li ParameterType_InputImageParameter
@@ -285,6 +296,7 @@ public:
    * Can be called for types :
    * \li ParameterType_Int
    * \li ParameterType_Float
+   * \li ParameterType_Double
    * \li ParameterType_Radius
    * \li ParameterType_Choice
    */
@@ -295,6 +307,7 @@ public:
    * Can be called for types :
    * \li ParameterType_Int
    * \li ParameterType_Float
+   * \li ParameterType_Double
    * \li ParameterType_Radius
    * \li ParameterType_Choice
    */
@@ -306,6 +319,7 @@ public:
    *
    * Can be called for types :
    * \li ParameterType_Float
+   * \li ParameterType_Double
    */
   void SetDefaultParameterFloat(std::string const& parameter, float value);
 
@@ -313,9 +327,28 @@ public:
    *
    * Can be called for types :
    * \li ParameterType_Float
+   * \li ParameterType_Double
    */
   float GetDefaultParameterFloat(std::string const& parameter);
 
+  /* Set a default double precision floating value, must be used in the
+   * DoInit when setting a value by default
+   * for the parameter
+   *
+   * Can be called for types :
+   * \li ParameterType_Float
+   * \li ParameterType_Double
+   */
+  void SetDefaultParameterDouble(std::string const& parameter, double value);
+
+  /* Get the default double precision floating value of a parameter
+   *
+   * Can be called for types :
+   * \li ParameterType_Float
+   * \li ParameterType_Double
+   */
+  double GetDefaultParameterDouble(std::string const& parameter);
+
   /** Set a default pixel type for an output image parameter
    *
    * \param[in] parameter Name of the output image parameter
@@ -348,7 +381,7 @@ public:
    */
   void SetMaximumParameterIntValue(std::string const& parameter, int value);
 
-  /* Set a minimum int value, must used in the
+  /* Set a minimum float value, must used in the
    * DoInit when setting a value by default
    * for the parameter
    *
@@ -357,7 +390,7 @@ public:
    */
   void SetMinimumParameterFloatValue(std::string const& parameter, float value);
 
-  /* Set a maximum int value, must used in the
+  /* Set a maximum float value, must used in the
    * DoInit when setting a value by default
    * for the parameter
    *
@@ -366,6 +399,26 @@ public:
    */
   void SetMaximumParameterFloatValue(std::string const& parameter, float value);
 
+  /* Set a minimum double precision float value, must used in the
+   * DoInit when setting a value by default
+   * for the parameter
+   *
+   * Can be called for types :
+   * \li ParameterType_Float
+   * \li ParameterType_Double
+   */
+  void SetMinimumParameterDoubleValue(std::string const& parameter, double value);
+
+  /* Set a maximum double precision value, must used in the
+   * DoInit when setting a value by default
+   * for the parameter
+   *
+   * Can be called for types :
+   * \li ParameterType_Float
+   * \li ParameterType_Double
+   */
+  void SetMaximumParameterDoubleValue(std::string const& parameter, double value);
+
 
   /**
    * Enable single selection mode for list view if status is true
@@ -411,6 +464,7 @@ public:
    * \li ParameterType_Int
    * \li ParameterType_Bool
    * \li ParameterType_Float
+   * \li ParameterType_Double
    * \li ParameterType_Radius
    * \li ParameterType_Choice
    */
@@ -420,9 +474,18 @@ public:
    *
    * Can be called for types :
    * \li ParameterType_Float
+   * \li ParameterType_Double
    */
   float GetParameterFloat(std::string const& parameter) const;
 
+  /* Get a double precision floating parameter value
+   *
+   * Can be called for types :
+   * \li ParameterType_Float
+   * \li ParameterType_Double
+   */
+  double GetParameterDouble(std::string const& parameter) const;
+
   /* Get a string parameter value
    *
    * Can be called for types :
@@ -609,6 +672,7 @@ public:
     *
     * Can be called for types :
     * \li ParameterType_Float
+    * \li ParameterType_Double
     * \li ParameterType_Int
     * \li ParameterType_Choice
     * \li ParameterType_Radius
diff --git a/Modules/Wrappers/ApplicationEngine/include/otbWrapperNumericalParameter.h b/Modules/Wrappers/ApplicationEngine/include/otbWrapperNumericalParameter.h
index 5be2d24308ec554479d0d19dc49506a9a95272a7..17abbe7007d00d4037bcb97c9df5234fe581ec45 100644
--- a/Modules/Wrappers/ApplicationEngine/include/otbWrapperNumericalParameter.h
+++ b/Modules/Wrappers/ApplicationEngine/include/otbWrapperNumericalParameter.h
@@ -127,6 +127,15 @@ public:
     return static_cast<float>(*m_Value);
   }
 
+  double ToDouble() const override
+  {
+    if (!HasValue())
+    {
+      itkExceptionMacro("Cannot convert parameter " << GetKey() << " to double (no value).");
+    }
+    return static_cast<double>(*m_Value);
+  }
+
   void FromInt(int value) override
   {
     SetValue(value);
@@ -185,7 +194,7 @@ public:
   itkNewMacro(Self);
   itkTypeMacro(NumericalParameter, Parameter);
 
-  virtual ParameterType GetType() const override
+  ParameterType GetType() const override
   {
     return ParameterType_Float;
   }
@@ -196,6 +205,28 @@ public:
   }
 };
 
+class OTBApplicationEngine_EXPORT DoubleParameter : public NumericalParameter<double>
+{
+public:
+  /** Standard class typedef */
+  typedef DoubleParameter                Self;
+  typedef itk::SmartPointer<Self>       Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  itkNewMacro(Self);
+  itkTypeMacro(NumericalParameter, Parameter);
+
+  ParameterType GetType() const override
+  {
+    return ParameterType_Double;
+  }
+
+  void FromDouble(double value) override
+  {
+    SetValue(value);
+  }
+};
+
 class OTBApplicationEngine_EXPORT IntParameter : public NumericalParameter<int>
 {
 public:
@@ -207,7 +238,7 @@ public:
   itkNewMacro(Self);
   itkTypeMacro(NumericalParameter, Parameter);
 
-  virtual ParameterType GetType() const override
+  ParameterType GetType() const override
   {
     return ParameterType_Int;
   }
@@ -228,7 +259,7 @@ public:
   /** RTTI support */
   itkTypeMacro(RAMParameter, Parameter);
 
-  virtual ParameterType GetType() const override
+  ParameterType GetType() const override
   {
     return ParameterType_RAM;
   }
@@ -260,7 +291,7 @@ public:
     return true;
   }
 
-  virtual ParameterType GetType() const override
+  ParameterType GetType() const override
   {
     return ParameterType_Radius;
   }
diff --git a/Modules/Wrappers/ApplicationEngine/include/otbWrapperParameter.h b/Modules/Wrappers/ApplicationEngine/include/otbWrapperParameter.h
index a82bdfc8e01b734d03c9a43ad4083b86c19f297b..b8ea41480e0914dc1089365aa1292d1739f2dcb5 100644
--- a/Modules/Wrappers/ApplicationEngine/include/otbWrapperParameter.h
+++ b/Modules/Wrappers/ApplicationEngine/include/otbWrapperParameter.h
@@ -147,11 +147,13 @@ public:
    */
   virtual int                      ToInt() const;
   virtual float                    ToFloat() const;
+  virtual double                   ToDouble() const;
   virtual std::string              ToString() const;
   virtual std::vector<std::string> ToStringList() const;
 
   virtual void FromInt(int);
   virtual void FromFloat(float);
+  virtual void FromDouble(double);
   virtual void FromString(const std::string&);
   virtual void FromStringList(const std::vector<std::string>&);
 
diff --git a/Modules/Wrappers/ApplicationEngine/include/otbWrapperTypes.h b/Modules/Wrappers/ApplicationEngine/include/otbWrapperTypes.h
index 2b2f0d0171052e121a0615c85dd66e051f39a357..3a8c602bb70a5f326fdc4343f53193e2882a3a7f 100644
--- a/Modules/Wrappers/ApplicationEngine/include/otbWrapperTypes.h
+++ b/Modules/Wrappers/ApplicationEngine/include/otbWrapperTypes.h
@@ -37,6 +37,7 @@ namespace Wrapper
 typedef enum {
   ParameterType_Int,
   ParameterType_Float,
+  ParameterType_Double,
   ParameterType_String,
   ParameterType_StringList,
   ParameterType_InputFilename,
@@ -64,6 +65,7 @@ namespace
 {
 constexpr char const* parameterTypesStrings[] = {"Int",
                                                  "Float",
+                                                 "Double",
                                                  "String",
                                                  "StringList",
                                                  "InputFilename",
diff --git a/Modules/Wrappers/ApplicationEngine/src/otbWrapperApplication.cxx b/Modules/Wrappers/ApplicationEngine/src/otbWrapperApplication.cxx
index c4b8205deef1071098ccf296611ece0f4a2f26df..9312e6350366651865a71ace60479b7b7d69ad05 100644
--- a/Modules/Wrappers/ApplicationEngine/src/otbWrapperApplication.cxx
+++ b/Modules/Wrappers/ApplicationEngine/src/otbWrapperApplication.cxx
@@ -411,6 +411,12 @@ void Application::SetParameterFloat(std::string const& key, float value, bool ha
   this->SetParameterUserValue(key, hasUserValueFlag);
 }
 
+void Application::SetParameterDouble(std::string const& key, double value, bool hasUserValueFlag)
+{
+  GetParameterByKey(key)->FromDouble(value);
+  this->SetParameterUserValue(key, hasUserValueFlag);
+}
+
 void Application::SetParameterString(std::string const& parameter, std::string value, bool hasUserValueFlag)
 {
   GetParameterByKey(parameter)->FromString(value);
@@ -1097,6 +1103,13 @@ void Application::SetDefaultParameterInt(std::string const& parameter, int value
     if (!hasUserValue)
       paramFloat->SetValue(static_cast<float>(value));
   }
+  else if (dynamic_cast<DoubleParameter*>(param))
+  {
+    DoubleParameter* paramDouble = dynamic_cast<DoubleParameter*>(param);
+    paramDouble->SetDefaultValue(static_cast<double>(value));
+    if (!hasUserValue)
+      paramDouble->SetValue(static_cast<double>(value));
+  }
   else if (dynamic_cast<RAMParameter*>(param))
   {
     RAMParameter* paramRAM = dynamic_cast<RAMParameter*>(param);
@@ -1125,6 +1138,11 @@ int Application::GetDefaultParameterInt(std::string const& parameter)
     FloatParameter* paramFloat = dynamic_cast<FloatParameter*>(param);
     ret                        = paramFloat->GetDefaultValue();
   }
+  else if (dynamic_cast<DoubleParameter*>(param))
+  {
+    DoubleParameter* paramDouble = dynamic_cast<DoubleParameter*>(param);
+    ret                        = paramDouble->GetDefaultValue();
+  }
   else if (dynamic_cast<RAMParameter*>(param))
   {
     RAMParameter* paramRAM = dynamic_cast<RAMParameter*>(param);
@@ -1154,6 +1172,23 @@ float Application::GetDefaultParameterFloat(std::string const& key)
   return param->GetDefaultValue();
 }
 
+void Application::SetDefaultParameterDouble(std::string const& key, double value)
+{
+  auto param = downcast_check<DoubleParameter>(GetParameterByKey(key));
+  param->SetDefaultValue(value);
+
+  if (!param->HasUserValue())
+  {
+    param->SetValue(value);
+  }
+}
+
+double Application::GetDefaultParameterDouble(std::string const& key)
+{
+  auto param = downcast_check<DoubleParameter>(GetParameterByKey(key));
+  return param->GetDefaultValue();
+}
+
 void Application::SetDefaultOutputPixelType(std::string const& key, ImagePixelType type)
 {
   auto param = downcast_check<OutputImageParameter>(GetParameterByKey(key));
@@ -1185,6 +1220,18 @@ void Application::SetMaximumParameterFloatValue(std::string const& key, float va
   param->SetMaximumValue(value);
 }
 
+void Application::SetMinimumParameterDoubleValue(std::string const& key, double value)
+{
+  auto param = downcast_check<DoubleParameter>(GetParameterByKey(key));
+  param->SetMinimumValue(value);
+}
+
+void Application::SetMaximumParameterDoubleValue(std::string const& key, double value)
+{
+  auto param = downcast_check<DoubleParameter>(GetParameterByKey(key));
+  param->SetMaximumValue(value);
+}
+
 void Application::SetListViewSingleSelectionMode(std::string const& key, bool status)
 {
   auto param = downcast_check<ListViewParameter>(GetParameterByKey(key));
@@ -1246,6 +1293,11 @@ float Application::GetParameterFloat(std::string const& key) const
   return GetParameterByKey(key)->ToFloat();
 }
 
+double Application::GetParameterDouble(std::string const& key) const
+{
+  return GetParameterByKey(key)->ToDouble();
+}
+
 std::string Application::GetParameterString(std::string const& key) const
 {
   return GetParameterByKey(key)->ToString();
@@ -1413,6 +1465,13 @@ std::vector<std::pair<std::string, std::string>> Application::GetOutputParameter
           oss << GetParameterFloat(*it);
           keyVal.second = oss.str();
         }
+        else if (type == ParameterType_Double)
+        {
+          std::ostringstream oss;
+          oss << std::setprecision(19);
+          oss << GetParameterDouble(*it);
+          keyVal.second = oss.str();
+        }
         else
         {
           keyVal.second = GetParameterAsString(*it);
diff --git a/Modules/Wrappers/ApplicationEngine/src/otbWrapperInputImageListParameter.cxx b/Modules/Wrappers/ApplicationEngine/src/otbWrapperInputImageListParameter.cxx
index 9dbddc7704c52a9c95ecb7496bda6d71e72117ad..912860fdf134a67d7ba7fe38882820e4632a1629 100644
--- a/Modules/Wrappers/ApplicationEngine/src/otbWrapperInputImageListParameter.cxx
+++ b/Modules/Wrappers/ApplicationEngine/src/otbWrapperInputImageListParameter.cxx
@@ -132,6 +132,8 @@ void InputImageListParameter::SetNthImage(std::size_t i, ImageBaseType* image)
   assert(image != nullptr);
 
   // Check input availability
+  if (image == nullptr)
+    itkExceptionMacro("Image number " << i << " is null. Please check the upstream source.");
   image->UpdateOutputInformation();
 
   // Build parameter.
@@ -142,6 +144,12 @@ void InputImageListParameter::SetNthImage(std::size_t i, ImageBaseType* image)
 /*****************************************************************************/
 void InputImageListParameter::AddImage(ImageBaseType* image)
 {
+
+  // Check input availability
+  if (image == nullptr)
+    itkExceptionMacro("Image is null. Please check the upstream source.");
+  image->UpdateOutputInformation();
+
   AddData(image, [this](auto i) -> auto { return this->FromImage(i); });
 }
 
diff --git a/Modules/Wrappers/ApplicationEngine/src/otbWrapperInputXML.cxx b/Modules/Wrappers/ApplicationEngine/src/otbWrapperInputXML.cxx
index 6178e4dd0876da5034d42821791bef73eb421a66..8f2f3fbf42fb5c46ec95c080c11343df785064e8 100644
--- a/Modules/Wrappers/ApplicationEngine/src/otbWrapperInputXML.cxx
+++ b/Modules/Wrappers/ApplicationEngine/src/otbWrapperInputXML.cxx
@@ -269,6 +269,12 @@ int Read(const std::string& filename, Application::Pointer this_)
       std::stringstream(value) >> floatValue;
       this_->SetParameterFloat(key, floatValue);
     }
+    else if (type == ParameterType_Double)
+    {
+      double doubleValue;
+      std::stringstream(value) >> doubleValue;
+      this_->SetParameterDouble(key, doubleValue);
+    }
     else if (type == ParameterType_StringList || type == ParameterType_ListView)
     {
       this_->SetParameterStringList(key, values);
diff --git a/Modules/Wrappers/ApplicationEngine/src/otbWrapperOutputXML.cxx b/Modules/Wrappers/ApplicationEngine/src/otbWrapperOutputXML.cxx
index 4f6e5a912bec9c9957ff2594aee82af48e38eddd..4b8a60619efd253dd33c94047fad0f12dfe13cc8 100644
--- a/Modules/Wrappers/ApplicationEngine/src/otbWrapperOutputXML.cxx
+++ b/Modules/Wrappers/ApplicationEngine/src/otbWrapperOutputXML.cxx
@@ -250,6 +250,13 @@ void ParseGroup(Application::Pointer app, TiXmlElement* n_App, const std::string
           oss << app->GetParameterFloat(key);
           value = oss.str();
         }
+        else if (type == ParameterType_Double)
+        {
+          std::ostringstream oss;
+          oss << std::setprecision(std::numeric_limits<double>::digits10 + 1);
+          oss << app->GetParameterDouble(key);
+          value = oss.str();
+        }
         else if (type == ParameterType_String || type == ParameterType_InputFilename || type == ParameterType_Directory || type == ParameterType_InputImage ||
                  type == ParameterType_InputVectorData || type == ParameterType_Choice || type == ParameterType_OutputVectorData ||
                  type == ParameterType_OutputFilename || type == ParameterType_Bool)
diff --git a/Modules/Wrappers/ApplicationEngine/src/otbWrapperParameter.cxx b/Modules/Wrappers/ApplicationEngine/src/otbWrapperParameter.cxx
index bd057c792c149b514f2972e01d9dc1c667273c55..bae8e60f378cc590e98ea17020b4c188638e9726 100644
--- a/Modules/Wrappers/ApplicationEngine/src/otbWrapperParameter.cxx
+++ b/Modules/Wrappers/ApplicationEngine/src/otbWrapperParameter.cxx
@@ -229,6 +229,11 @@ float Parameter::ToFloat() const
   TypeError("float");
 }
 
+double Parameter::ToDouble() const
+{
+  TypeError("double");
+}
+
 std::string Parameter::ToString() const
 {
   TypeError("std::string");
@@ -249,6 +254,11 @@ void Parameter::FromFloat(float)
   TypeError("float");
 }
 
+void Parameter::FromDouble(double)
+{
+  TypeError("double");
+}
+
 void Parameter::FromString(const std::string&)
 {
   TypeError("std::string");
diff --git a/Modules/Wrappers/ApplicationEngine/src/otbWrapperParameterGroup.cxx b/Modules/Wrappers/ApplicationEngine/src/otbWrapperParameterGroup.cxx
index 803212e8724ff1ac17517c5203a9618c67304e00..6212dfb2119d12b17f1bb8df773b327d0cc29bc9 100644
--- a/Modules/Wrappers/ApplicationEngine/src/otbWrapperParameterGroup.cxx
+++ b/Modules/Wrappers/ApplicationEngine/src/otbWrapperParameterGroup.cxx
@@ -231,6 +231,11 @@ void ParameterGroup::AddParameter(ParameterType type, std::string paramKey, std:
       newParam = FloatParameter::New();
     }
     break;
+    case ParameterType_Double:
+    {
+      newParam = DoubleParameter::New();
+    }
+    break;
     case ParameterType_String:
     {
       newParam = StringParameter::New();
diff --git a/Modules/Wrappers/ApplicationEngine/test/CMakeLists.txt b/Modules/Wrappers/ApplicationEngine/test/CMakeLists.txt
index 070eed0d442c9345576c3ca30e526339c3ad4256..ffb6ce0cd16834729f93c4325632ed7a8b7f75a9 100644
--- a/Modules/Wrappers/ApplicationEngine/test/CMakeLists.txt
+++ b/Modules/Wrappers/ApplicationEngine/test/CMakeLists.txt
@@ -56,6 +56,10 @@ otb_add_test(NAME owTvFloatParameter COMMAND otbApplicationEngineTestDriver
   otbWrapperFloatParameterTest
   )
 
+otb_add_test(NAME owTvDoubleParameter COMMAND otbApplicationEngineTestDriver
+  otbWrapperDoubleParameterTest
+  )
+
 otb_add_test(NAME owTvIntParameter COMMAND otbApplicationEngineTestDriver
   otbWrapperIntParameterTest
   )
diff --git a/Modules/Wrappers/ApplicationEngine/test/otbApplicationEngineTestDriver.cxx b/Modules/Wrappers/ApplicationEngine/test/otbApplicationEngineTestDriver.cxx
index 9fa82c2c1d4f0f70a369082e52b17d5e9a8b3269..541f3c6d2ba6458361f5d374c5f329ae311203df 100644
--- a/Modules/Wrappers/ApplicationEngine/test/otbApplicationEngineTestDriver.cxx
+++ b/Modules/Wrappers/ApplicationEngine/test/otbApplicationEngineTestDriver.cxx
@@ -24,6 +24,7 @@ void RegisterTests()
 {
   REGISTER_TEST(otbWrapperInputImageParameterTest);
   REGISTER_TEST(otbWrapperFloatParameterTest);
+  REGISTER_TEST(otbWrapperDoubleParameterTest);
   REGISTER_TEST(otbWrapperIntParameterTest);
   REGISTER_TEST(otbWrapperRAMParameterTest);
   REGISTER_TEST(otbWrapperStringParameterTest1);
diff --git a/Modules/Wrappers/ApplicationEngine/test/otbWrapperNumericalParameterTest.cxx b/Modules/Wrappers/ApplicationEngine/test/otbWrapperNumericalParameterTest.cxx
index fc3cb8553755f43678045796d2b82345242af949..a101fdb85de910eb9f7de301c5bd49cc628b4ede 100644
--- a/Modules/Wrappers/ApplicationEngine/test/otbWrapperNumericalParameterTest.cxx
+++ b/Modules/Wrappers/ApplicationEngine/test/otbWrapperNumericalParameterTest.cxx
@@ -23,6 +23,7 @@
 #endif
 
 #include "otbWrapperNumericalParameter.h"
+#include <typeinfo>
 
 using namespace otb::Wrapper;
 
@@ -35,6 +36,15 @@ void assert_equal(const T& a, const T& b)
   }
 }
 
+template <typename T>
+void assert_close(const T& a, const T& b)
+{
+  if (std::abs(a-b) > 1e-7)
+  {
+    itkGenericExceptionMacro("assert_close failed");
+  }
+}
+
 int otbWrapperFloatParameterTest(int, char* [])
 {
   auto param = FloatParameter::New();
@@ -50,6 +60,7 @@ int otbWrapperFloatParameterTest(int, char* [])
     assert_equal(param->GetValue(), val);
     assert_equal(param->ToInt(), int(val));
     assert_equal(param->ToFloat(), val);
+    assert_equal(param->ToDouble(), (double)val);
     assert_equal(param->ToString(), std::string("42.005"));
   }
 
@@ -59,6 +70,7 @@ int otbWrapperFloatParameterTest(int, char* [])
     assert_equal(param->GetValue(), val);
     assert_equal(param->ToInt(), int(val));
     assert_equal(param->ToFloat(), val);
+    assert_equal(param->ToDouble(), (double)val);
     assert_equal(param->ToString(), std::string("-6.5"));
   }
 
@@ -69,6 +81,50 @@ int otbWrapperFloatParameterTest(int, char* [])
     assert_equal(param->GetValue(), val);
     assert_equal(param->ToInt(), int(val));
     assert_equal(param->ToFloat(), val);
+    assert_equal(param->ToDouble(), (double)val);
+    assert_equal(param->ToString(), std::string("-100.01"));
+  }
+
+  return EXIT_SUCCESS;
+}
+
+int otbWrapperDoubleParameterTest(int, char* [])
+{
+  auto param = DoubleParameter::New();
+
+  param->SetKey("mykey");
+  param->SetDescription("My description.");
+
+  assert_equal(param->GetType(), ParameterType_Double);
+
+  { // SetValue
+    const double val = 42.005;
+    param->SetValue(val);
+    assert_equal(param->GetValue(), val);
+    assert_equal(param->ToInt(), int(val));
+    assert_close(param->ToFloat(), (float)val);
+    assert_close(param->ToDouble(), val);
+    assert_equal(param->ToString(), std::string("42.005"));
+  }
+
+  { // FromDouble
+    const double val = -6.5;
+    param->FromDouble(val);
+    assert_equal(param->GetValue(), val);
+    assert_equal(param->ToInt(), int(val));
+    assert_close(param->ToFloat(), (float)val);
+    assert_close(param->ToDouble(), val);
+    assert_equal(param->ToString(), std::string("-6.5"));
+  }
+
+  { // FromString
+    const std::string str = "-100.01";
+    const double       val = -100.01;
+    param->FromString(str);
+    assert_equal(param->GetValue(), val);
+    assert_equal(param->ToInt(), int(val));
+    assert_close(param->ToFloat(), (float)val);
+    assert_close(param->ToDouble(), val);
     assert_equal(param->ToString(), std::string("-100.01"));
   }
 
@@ -90,6 +146,7 @@ int otbWrapperIntParameterTest(int, char* [])
     assert_equal(param->GetValue(), val);
     assert_equal(param->ToInt(), val);
     assert_equal(param->ToFloat(), float(val));
+    assert_equal(param->ToDouble(), double(val));
     assert_equal(param->ToString(), std::string("42"));
   }
 
@@ -100,6 +157,7 @@ int otbWrapperIntParameterTest(int, char* [])
     assert_equal(param->GetValue(), val);
     assert_equal(param->ToInt(), val);
     assert_equal(param->ToFloat(), float(val));
+    assert_equal(param->ToDouble(), double(val));
     assert_equal(param->ToString(), std::string("-100"));
   }
 
diff --git a/Modules/Wrappers/CommandLine/src/otbWrapperCommandLineLauncher.cxx b/Modules/Wrappers/CommandLine/src/otbWrapperCommandLineLauncher.cxx
index c23738b01eaf1cf62da2a2ee9035646c8c390f25..59e1d907509db11616cd6b28a2f4cb16d49cf7dd 100644
--- a/Modules/Wrappers/CommandLine/src/otbWrapperCommandLineLauncher.cxx
+++ b/Modules/Wrappers/CommandLine/src/otbWrapperCommandLineLauncher.cxx
@@ -431,7 +431,7 @@ CommandLineLauncher::ParamResultType CommandLineLauncher::LoadParameters()
           // Multiple values parameters
           m_Application->SetParameterStringList(paramKey, values);
         }
-        else if (type == ParameterType_Choice || type == ParameterType_Float || type == ParameterType_Int || type == ParameterType_Radius ||
+        else if (type == ParameterType_Choice || type == ParameterType_Float || type == ParameterType_Double || type == ParameterType_Int || type == ParameterType_Radius ||
                  type == ParameterType_Directory || type == ParameterType_String || type == ParameterType_InputFilename ||
                  type == ParameterType_OutputFilename || type == ParameterType_InputImage || type == ParameterType_OutputImage ||
                  type == ParameterType_InputVectorData || type == ParameterType_OutputVectorData || type == ParameterType_RAM || type == ParameterType_Bool)
@@ -708,6 +708,10 @@ std::string CommandLineLauncher::DisplayParameterHelp(const Parameter::Pointer&
   {
     oss << "<float>         ";
   }
+  else if (type == ParameterType_Double)
+  {
+    oss << "<double>         ";
+  }
   else if (type == ParameterType_InputFilename || type == ParameterType_OutputFilename || type == ParameterType_Directory || type == ParameterType_InputImage ||
            type == ParameterType_InputVectorData || type == ParameterType_OutputVectorData || type == ParameterType_String || type == ParameterType_Choice ||
            (type == ParameterType_ListView && singleSelectionForListView))
diff --git a/Modules/Wrappers/QGIS/src/otbQgisDescriptor.cxx b/Modules/Wrappers/QGIS/src/otbQgisDescriptor.cxx
index 226d5c9a48244d4a8da89776ff58e293456de341..2b27702fdad6eba4231b52a2afb4026090d1843c 100644
--- a/Modules/Wrappers/QGIS/src/otbQgisDescriptor.cxx
+++ b/Modules/Wrappers/QGIS/src/otbQgisDescriptor.cxx
@@ -57,6 +57,7 @@ int main(int argc, char* argv[])
   parameterTypeToString[ParameterType_Bool]                = "QgsProcessingParameterBoolean";
   parameterTypeToString[ParameterType_Int]                 = "QgsProcessingParameterNumber";
   parameterTypeToString[ParameterType_Float]               = "QgsProcessingParameterNumber";
+  parameterTypeToString[ParameterType_Double]              = "QgsProcessingParameterNumber";
   parameterTypeToString[ParameterType_RAM]                 = "QgsProcessingParameterNumber";
   parameterTypeToString[ParameterType_Radius]              = "QgsProcessingParameterNumber";
   parameterTypeToString[ParameterType_Choice]              = "OTBParameterChoice";
@@ -201,7 +202,7 @@ int main(int argc, char* argv[])
         default_value = param->HasValue() ? appli->GetParameterAsString(name) : "0";
       }
     }
-    else if (type == ParameterType_Float)
+    else if (type == ParameterType_Float || type == ParameterType_Double)
     {
       dFile << "|QgsProcessingParameterNumber.Double";
       default_value = param->HasValue() ? appli->GetParameterAsString(name) : "0";
diff --git a/Modules/Wrappers/SWIG/src/otbApplication.i b/Modules/Wrappers/SWIG/src/otbApplication.i
index 2f4139254868b91935628529664031ff5fff91e1..226e8136cd0c45fce69ab3f6a398ca83fc2baeb3 100644
--- a/Modules/Wrappers/SWIG/src/otbApplication.i
+++ b/Modules/Wrappers/SWIG/src/otbApplication.i
@@ -78,6 +78,7 @@ namespace Wrapper
   {
     ParameterType_Int,
     ParameterType_Float,
+    ParameterType_Double,
     ParameterType_String,
     ParameterType_StringList,
     ParameterType_InputFilename,
@@ -217,6 +218,10 @@ public:
     self.UpdateParameters()
 %}
 
+%pythonappend Application::SetParameterDouble %{
+    self.UpdateParameters()
+%}
+
 %pythonappend Application::SetParameterString %{
     self.UpdateParameters()
 %}
@@ -307,6 +312,7 @@ public:
 
   void SetParameterInt(std::string parameter, int value, bool hasUserValueFlag = true);
   void SetParameterFloat(std::string parameter, float value, bool hasUserValueFlag = true);
+  void SetParameterDouble(std::string parameter, double value, bool hasUserValueFlag = true);
   void SetParameterString(std::string parameter, std::string value, bool hasUserValueFlag = true);
   void SetParameterStringList(std::string parameter, std::vector<std::string> values, bool hasUserValueFlag = true);
 
@@ -316,6 +322,7 @@ public:
 
   int GetParameterInt(std::string parameter);
   float GetParameterFloat(std::string parameter);
+  double GetParameterDouble(std::string parameter);
   std::string GetParameterString(std::string parameter);
   std::vector<std::string> GetParameterStringList(std::string parameter);
   std::string GetParameterAsString(std::string paramKey);
@@ -625,6 +632,7 @@ class ApplicationProxy(object):
         ParameterType_Int : 'ParameterType_Int',
         ParameterType_Radius : 'ParameterType_Radius',
         ParameterType_RAM : 'ParameterType_RAM',
+        ParameterType_Double : 'ParameterType_Double',
         ParameterType_Float : 'ParameterType_Float',
         ParameterType_Choice : 'ParameterType_Choice',
         ParameterType_Group : 'ParameterType_Group',
@@ -658,6 +666,8 @@ class ApplicationProxy(object):
         return self.SetParameterInt(paramKey, value)
       elif paramType in [ParameterType_Float]:
         return self.SetParameterFloat(paramKey, value)
+      elif paramType in [ParameterType_Double]:
+        return self.SetParameterDouble(paramKey, value)
       elif paramType in [ParameterType_Bool]:
         return self.SetParameterString(paramKey, str(value) )
       elif paramType in [ParameterType_Group]:
@@ -692,6 +702,8 @@ class ApplicationProxy(object):
         return self.GetParameterInt(paramKey)
       elif paramType in [ParameterType_Float]:
         return self.GetParameterFloat(paramKey)
+      elif paramType in [ParameterType_Double]:
+        return self.GetParameterDouble(paramKey)
       elif paramType in [ParameterType_Bool]:
         return bool(self.GetParameterInt(paramKey))
       elif paramType in [ParameterType_Group, ParameterType_Choice]: