diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index f59baac2626b15469be898b1003c153933de6885..8b22dbf873348f64754cfab489f297415c198c90 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -9,7 +9,7 @@
 # GIT_STRATEGY=fetch (https://gitlab.com/gitlab-org/gitlab-runner/issues/3318)
 
 variables:
-  BUILD_IMAGE_REGISTRY: $CI_REGISTRY/gbonnefille/otb-build-env
+  BUILD_IMAGE_REGISTRY: $CI_REGISTRY/orfeotoolbox/otb-build-env
   DOCKER_DRIVER: overlay2
   GIT_DEPTH: "3"
   # Disable automatic checkout to let us fetch LFS before
@@ -44,7 +44,7 @@ fast-build:
   extends: .general
   only: [merge_requests, branches]
   stage: precheck
-  image: $CI_REGISTRY/gpasero/otb/otb-install-ubuntu-native
+  image: $BUILD_IMAGE_REGISTRY/otb-ubuntu-native-develop:latest
   before_script:
     - export GIT_LFS_SKIP_SMUDGE=1
     - git checkout $CI_COMMIT_REF_NAME
@@ -57,6 +57,7 @@ fast-build:
   extends: .general
   only: [merge_requests]
   stage: build
+  dependencies: []
 
 debian-build:
   extends: .common-build
diff --git a/.gitlab/merge_request_templates/request_for_changes.md b/.gitlab/merge_request_templates/request_for_changes.md
index 31bf9cb52c20e5b67e82be1158a46b7bbd7ca268..79acdadae838587114d6aa593ddcf0a78a8ef3bd 100644
--- a/.gitlab/merge_request_templates/request_for_changes.md
+++ b/.gitlab/merge_request_templates/request_for_changes.md
@@ -45,3 +45,4 @@ The copyright owner is *COPYRIGHT OWNER (OR OWNER'S AGENT)* and has signed the O
 - The feature branch is (reasonably) up-to-date with the base branch
 - Dashboard is green
 - Copyright owner has signed the ORFEO ToolBox Contributor License Agreement
+- Optionally, run `git diff develop... -U0 --no-color | clang-format-diff.py -p1 -i` on latest changes and commit
diff --git a/CI/ubuntu-18.04-llvm-nodoc.cmake b/CI/ubuntu-18.04-llvm-nodoc.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..27ab15d0cd7e0e9c9a7802f2831baf5ef5e5f180
--- /dev/null
+++ b/CI/ubuntu-18.04-llvm-nodoc.cmake
@@ -0,0 +1,32 @@
+#
+# Copyright (C) 2005-2019 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.
+#
+
+# Configuration options for ubuntu-18.04-llvm-nodoc
+
+set(site_option
+"opencv_INCLUDE_DIR:PATH=/usr/include
+CMAKE_C_COMPILER:STRING=clang
+CMAKE_CXX_COMPILER:STRING=clang++
+CMAKE_EXE_LINKER_FLAGS:STRING=-fuse-ld=lld
+CMAKE_MODULE_LINKER_FLAGS:STRING=-fuse-ld=lld
+CMAKE_SHARED_LINKER_FLAGS:STRING=-fuse-ld=lld
+CMAKE_C_COMPILER_LAUNCHER:STRING=ccache
+CMAKE_CXX_COMPILER_LAUNCHER:STRING=ccache
+OTB_USE_SHARK:BOOL=OFF")
diff --git a/CI/ubuntu-18.04-llvm.cmake b/CI/ubuntu-18.04-llvm.cmake
index 6480ab2f31b9ddb22b56a15d61ac305a85e85423..9b838da433a06cff9f26281578ad56c0b20c52ec 100644
--- a/CI/ubuntu-18.04-llvm.cmake
+++ b/CI/ubuntu-18.04-llvm.cmake
@@ -49,6 +49,7 @@ OTB_DOXYGEN_ITK_DOXYGEN_URL:STRING=\"https://itk.org/Doxygen413/html\"
   # See otb-devutils/Scripts/tagfile_fix.py
   message(STATUS "Get resources for Doxygen build ...")
   execute_process(COMMAND wget https://www.orfeo-toolbox.org/packages/archives/Doxygen/InsightDoxygenDocTag-4.13.0.gz
-                  COMMAND gzip -d InsightDoxygenDocTag-4.13.0.gz
+                  WORKING_DIRECTORY ${CTEST_BINARY_DIRECTORY})
+  execute_process(COMMAND gzip -d InsightDoxygenDocTag-4.13.0.gz
                   WORKING_DIRECTORY ${CTEST_BINARY_DIRECTORY})
 endif()
diff --git a/Data/Baseline/Examples/Radiometry/ARVIMultiChannelRAndBAndNIRVegetationIndex.tif b/Data/Baseline/Examples/Radiometry/ARVIMultiChannelRAndBAndNIRVegetationIndex.tif
index 0225c9d196f05cec5b484ecdcc45a096f2f3ce03..4361cfbb8aafbcc98f9c8ac99e40739f62f9aaf1 100644
--- a/Data/Baseline/Examples/Radiometry/ARVIMultiChannelRAndBAndNIRVegetationIndex.tif
+++ b/Data/Baseline/Examples/Radiometry/ARVIMultiChannelRAndBAndNIRVegetationIndex.tif
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:42df562a687bdf13d65cacc6da48666b28da72efa79d7d668e607cf82eba4c11
-size 68855
+oid sha256:91667b86ce412524444bde780eee2c7f0baff0c05ecdc5b66a49aca39bc60dda
+size 65842
diff --git a/Data/Baseline/Examples/Radiometry/NDVIRAndNIRVegetationIndex.tif b/Data/Baseline/Examples/Radiometry/NDVIRAndNIRVegetationIndex.tif
deleted file mode 100644
index 0882267d3bfa84e1cf7089f0f6f7aeb50c632151..0000000000000000000000000000000000000000
--- a/Data/Baseline/Examples/Radiometry/NDVIRAndNIRVegetationIndex.tif
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:5727e5525e7597823337b25ca4a86b18f96606d68a8fd8b624c6d3a8796372b4
-size 63578
diff --git a/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_DEMGTIFF.tif b/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_DEMGTIFF.tif
index 9662d6048e1b5bc4eef3ad7c00b24ecfb24db433..3b1b11efef094cc93a80d22849d312da21c7c4f9 100644
--- a/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_DEMGTIFF.tif
+++ b/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_DEMGTIFF.tif
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:604ab2a554cb7f8b8bd2f793da761b944ea1585b267b19063d64fbdd8fe9bd75
-size 2471604
+oid sha256:39a2dd6ff12fd1ad13ce5c3014dac39b546aceb5e5d1fbb8a4bb7e58d06fcbeb
+size 2002472
diff --git a/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_DEMSRTM.tif b/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_DEMSRTM.tif
index 0c871089a6f291752cdfe701b38f59c8b8eb9db3..419e04824bcdb961f1fcc0e54fc1c5ecca21d05f 100644
--- a/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_DEMSRTM.tif
+++ b/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_DEMSRTM.tif
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:d0bf90e6e04a89d3876af4272815b3a47dd3cd2805ed3bdb6900a68ec3fb83a3
-size 2470925
+oid sha256:d684f0d16382f9a7402373b4e17dbc74c8d38b2e2a5ef0c95560374f4ae23a3f
+size 2002472
diff --git a/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_noDEM.tif b/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_noDEM.tif
index 22b97cb747f2ea635594d2d3fe1d3c99a82ad5b6..6a12b74e395f7e78dc669c19ea3463408fb01ebd 100644
--- a/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_noDEM.tif
+++ b/Data/Baseline/OTB/Images/prTvOrthoRectification_sentinel1_noDEM.tif
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:bfcce43b020c8838f1b710dd62cff8781dff095e88e056edf2df8863dc0ae3af
-size 2471252
+oid sha256:1b8c6e80b58829d09ec7ad2b7cca3104c0fd1676a203361bda5493f0060d8fee
+size 2002472
diff --git a/Documentation/Cookbook/rst/C++/UserGuide.rst b/Documentation/Cookbook/rst/C++/UserGuide.rst
index 10304efb75adc6ae0782a52ed781e162e07c37e3..c557ef2e761c22797a3b6ebe000ca350507e2fbd 100644
--- a/Documentation/Cookbook/rst/C++/UserGuide.rst
+++ b/Documentation/Cookbook/rst/C++/UserGuide.rst
@@ -618,8 +618,6 @@ vegetation had been demonstrated, users tended to also use the NDVI to
 quantify the photosynthetic capacity of plant canopies. This, however,
 can be a rather more complex undertaking if not done properly.
 
-* NDVI. See example :ref:`NDVIRAndNIRVegetationIndexImageFilter.cxx`.
-
 * ARVI. See example :ref:`ARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx`.
 
 * AVI. See example :ref:`AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter.cxx`.
diff --git a/Documentation/Cookbook/rst/CompilingOTBFromSource.rst b/Documentation/Cookbook/rst/CompilingOTBFromSource.rst
index 7961cc02b4c01a66b884a02cc678759f70bc5a34..84dccaaf529465b0f31124b78401f8e5f0e3c593 100644
--- a/Documentation/Cookbook/rst/CompilingOTBFromSource.rst
+++ b/Documentation/Cookbook/rst/CompilingOTBFromSource.rst
@@ -23,7 +23,6 @@ majority) are optional and can be activated or not during the build
 process:
 
 .. table:: External libraries used in OTB
-    :widths: 50 20 30
     :align: center
 
     +------------------------------------------------------------------+-----------------------+--------------------------+
diff --git a/Documentation/Cookbook/rst/recipes/featextract.rst b/Documentation/Cookbook/rst/recipes/featextract.rst
index 64afc6814ead8eed8f69d63f7b41af1aaab1f47e..23267e65276dc6d011b2a4967ff2dc73bf335792 100644
--- a/Documentation/Cookbook/rst/recipes/featextract.rst
+++ b/Documentation/Cookbook/rst/recipes/featextract.rst
@@ -111,37 +111,62 @@ The *RadiometricIndices* application has the following input parameters:
    (default value is 1)
 
 -``-list`` the list of available radiometric indices (default value is
-   Vegetation:NDVI)
+   ``Vegetation:NDVI``)
 
-The available radiometric indices to be listed into -list with their
-relevant channels in brackets are:
+Note that band numbering starts at 1.
 
-::
+The available radiometric indices to be listed into ``-list`` with their
+relevant channels in brackets are:
 
-    Vegetation:NDVI - Normalized difference vegetation index (Red, NIR)
-    Vegetation:TNDVI - Transformed normalized difference vegetation index (Red, NIR)
-    Vegetation:RVI - Ratio vegetation index (Red, NIR)
-    Vegetation:SAVI - Soil adjusted vegetation index (Red, NIR)
-    Vegetation:TSAVI - Transformed soil adjusted vegetation index (Red, NIR)
-    Vegetation:MSAVI - Modified soil adjusted vegetation index (Red, NIR)
-    Vegetation:MSAVI2 - Modified soil adjusted vegetation index 2 (Red, NIR)
-    Vegetation:GEMI - Global environment monitoring index (Red, NIR)
-    Vegetation:IPVI - Infrared percentage vegetation index (Red, NIR)
-
-    Water:NDWI - Normalized difference water index (Gao 1996) (NIR, MIR)
-    Water:NDWI2 - Normalized difference water index (Mc Feeters 1996) (Green, NIR)
-    Water:MNDWI - Modified normalized difference water index (Xu 2006) (Green, MIR)
-    Water:NDPI - Normalized difference pond index (Lacaux et al.) (MIR, Green)
-    Water:NDTI - Normalized difference turbidity index (Lacaux et al.) (Red, Green)
-
-    Soil:RI - Redness index (Red, Green)
-    Soil:CI - Color index (Red, Green)
-    Soil:BI - Brightness index (Red, Green)
-    Soil:BI2 - Brightness index 2 (NIR, Red, Green)
++------------------------------+----------------------------------------------------------+---------------+
+|Index                         |Description                                               |Required bands |
++==============================+==========================================================+===============+
+|Vegetation:NDVI               |Normalized difference vegetation index                    |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation:TNDVI              |Transformed normalized difference vegetation index        |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation:RVI                |Ratio vegetation index                                    |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation:SAVI               |Soil adjusted vegetation index                            |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation:TSAVI              |Transformed soil adjusted vegetation index                |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation:MSAVI              |Modified soil adjusted vegetation index                   |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation:MSAVI2             |Modified soil adjusted vegetation index 2                 |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation:GEMI               |Global environment monitoring index                       |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation:IPVI               |Infrared percentage vegetation index                      |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation:LAIFromNDVILog     |Leaf Area Index from log NDVI                             |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation::LAIFromReflLinear |Leaf Area Index from reflectances with linear combination |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Vegetation::LAIFromNDVIFormo  |Leaf Area Index from Formosat 2  TOC                      |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Water:NDWI                    |Normalized difference water index (Gao 1996)              |NIR, MIR       |
++------------------------------+----------------------------------------------------------+---------------+
+|Water:NDWI2                   |Normalized difference water index (Mc Feeters 1996)       |Green, NIR     |
++------------------------------+----------------------------------------------------------+---------------+
+|Water:MNDWI                   |Modified normalized difference water index (Xu 2006)      |Green, MIR     |
++------------------------------+----------------------------------------------------------+---------------+
+|Water:NDTI                    |Normalized difference turbidity index (Lacaux et al.)     |Green, Red     |
++------------------------------+----------------------------------------------------------+---------------+
+|Soil:RI                       |Redness index                                             |Green, Red     |
++------------------------------+----------------------------------------------------------+---------------+
+|Soil:CI                       |Color index                                               |Green, Red     |
++------------------------------+----------------------------------------------------------+---------------+
+|Soil:BI                       |Brightness index                                          |Green, Red     |
++------------------------------+----------------------------------------------------------+---------------+
+|Soil:BI2                      |Brightness index 2                                        |Green, Red, NIR|
++------------------------------+----------------------------------------------------------+---------------+
+|BuiltUp:ISU                   |Built Surfaces Index                                      |Red, NIR       |
++------------------------------+----------------------------------------------------------+---------------+
 
 The application can be used as follows, which would produce an output image
-containing 3 bands, respectively with the Vegetation:NDVI, Vegetation:RVI and
-Vegetation:IPVI radiometric indices in this exact order:
+containing 3 bands, respectively with the ``Vegetation:NDVI``, ``Vegetation:RVI`` and
+``Vegetation:IPVI`` radiometric indices in this exact order:
 
 ::
 
@@ -154,7 +179,7 @@ Vegetation:IPVI radiometric indices in this exact order:
                                               Vegetation:IPVI 
 
 or as follows, which would produce a single band output image with the
-Water:NDWI2 radiometric index:
+``Water:NDWI2`` radiometric index:
 
 ::
 
diff --git a/Examples/OBIA/RadiometricAttributesLabelMapFilterExample.cxx b/Examples/OBIA/RadiometricAttributesLabelMapFilterExample.cxx
index 9991cf7db9e6627f3341c92a2908829044d33509..1ab6443c048376eafa4f9cbd31dabafba345711c 100644
--- a/Examples/OBIA/RadiometricAttributesLabelMapFilterExample.cxx
+++ b/Examples/OBIA/RadiometricAttributesLabelMapFilterExample.cxx
@@ -107,7 +107,7 @@ int main(int argc, char* argv[])
   typedef otb::BandsStatisticsAttributesLabelMapFilter<LabelMapType, VectorImageType> RadiometricLabelMapFilterType;
   typedef otb::AttributesMapOpeningLabelMapFilter<LabelMapType>                       OpeningLabelMapFilterType;
   typedef itk::LabelMapToBinaryImageFilter<LabelMapType, MaskImageType>               LabelMapToBinaryImageFilterType;
-  typedef itk::UnaryFunctorImageFilter<VectorImageType, ImageType,otb::Functor::NDVI<PixelType,PixelType,PixelType> > NDVIImageFilterType;
+  typedef itk::UnaryFunctorImageFilter<VectorImageType, ImageType, otb::Functor::NDVI<PixelType, PixelType>> NDVIImageFilterType;
   typedef otb::ImageToVectorImageCastFilter<ImageType, VectorImageType>   ImageToVectorImageCastFilterType;
 
   ReaderType::Pointer reader = ReaderType::New();
@@ -166,8 +166,8 @@ int main(int argc, char* argv[])
   //  In our case, statistics are computed on the NDVI coefficient on each label object.
   NDVIImageFilterType::Pointer ndviImageFilter = NDVIImageFilterType::New();
 
-  ndviImageFilter->GetFunctor().SetRedIndex(3);
-  ndviImageFilter->GetFunctor().SetNIRIndex(4);
+  ndviImageFilter->GetFunctor().SetBandIndex(CommonBandNames::RED, 3);
+  ndviImageFilter->GetFunctor().SetBandIndex(CommonBandNames::NIR, 4);
   ndviImageFilter->SetInput(vreader->GetOutput());
 
   ImageToVectorImageCastFilterType::Pointer ndviVectorImageFilter = ImageToVectorImageCastFilterType::New();
diff --git a/Examples/Radiometry/ARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx b/Examples/Radiometry/ARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx
index 7f051c3e438c87c6dae83164058571601eccb65c..6845750ac2c420db14567e6aa6dc0c85b720072e 100644
--- a/Examples/Radiometry/ARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx
+++ b/Examples/Radiometry/ARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx
@@ -26,8 +26,7 @@
                                                         Output/pretty_ARVIMultiChannelRAndBAndNIRVegetationIndex.png \
                                                         1 \
                                                         3 \
-                                                        2 \
-                                                        0.6
+                                                        2
 */
 
 
@@ -87,11 +86,11 @@
 
 int main(int argc, char* argv[])
 {
-  if (argc < 8)
+  if (argc < 7)
   {
     std::cerr << "Missing Parameters " << std::endl;
     std::cerr << "Usage: " << argv[0];
-    std::cerr << " inputImage , outputImage , prettyInput , prettyOutput , redChannel , blueChannel , nirChannel , gama" << std::endl;
+    std::cerr << " inputImage , outputImage , prettyInput , prettyOutput , redChannel , blueChannel , nirChannel" << std::endl;
     return 1;
   }
 
@@ -114,7 +113,7 @@ int main(int argc, char* argv[])
   // Note that we also can use other functors which operate with the
   // Red, Blue and Nir channels such as EVI, ARVI and TSARVI.
 
-  typedef otb::Functor::ARVI<InputPixelType, InputPixelType, InputPixelType, OutputPixelType> FunctorType;
+  typedef otb::Functor::ARVI<InputPixelType, OutputPixelType> FunctorType;
 
   // The
   // \doxygen{itk}{UnaryFunctorImageFilter}
@@ -134,17 +133,9 @@ int main(int argc, char* argv[])
   writer->SetFileName(argv[2]);
 
   // The three used index bands (red, blue and NIR) are declared.
-
-  filter->GetFunctor().SetRedIndex(::atoi(argv[5]));
-  filter->GetFunctor().SetBlueIndex(::atoi(argv[6]));
-  filter->GetFunctor().SetNIRIndex(::atoi(argv[7]));
-
-  // The $\gamma$ parameter is set. The
-  // \doxygen{otb::Functor}{ARVI}
-  // class sets the default value of $\gamma$ to $0.5$.  This parameter
-  // is used to reduce the atmospheric effect on a global scale.
-
-  filter->GetFunctor().SetGamma(::atof(argv[8]));
+  filter->GetFunctor().SetBandIndex(CommonBandNames::RED, ::atoi(argv[5]));
+  filter->GetFunctor().SetBandIndex(CommonBandNames::BLUE, ::atoi(argv[6]));
+  filter->GetFunctor().SetBandIndex(CommonBandNames::NIR, ::atoi(argv[7]));
 
   // The filter input is linked to the reader output and
   // the filter output is linked to the writer input.
diff --git a/Examples/Radiometry/AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter.cxx b/Examples/Radiometry/AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter.cxx
index 68653b4e0d5bb5a4a17e240bb6068bb8f94264e1..d42f8cd56e1227d88f1fe1b8bd956538a4bb80aa 100644
--- a/Examples/Radiometry/AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter.cxx
+++ b/Examples/Radiometry/AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter.cxx
@@ -26,10 +26,7 @@
                                                        Output/pretty_AVIMultiChannelRAndGAndNIRVegetationIndex.png \
                                                        3 \
                                                        2 \
-                                                       4 \
-                                                       660 \
-                                                       560 \
-                                                       830
+                                                       4
 */
 
 
@@ -76,12 +73,11 @@
 
 int main(int argc, char* argv[])
 {
-  if (argc < 11)
+  if (argc < 8)
   {
     std::cerr << "Missing Parameters " << std::endl;
     std::cerr << "Usage: " << argv[0];
-    std::cerr << " inputImage , outputImage , prettyInput , prettyOutput , redChannel , greenChannel , nirChannel ,";
-    std::cerr << " lambdaR, lambdaG, lambdaNIR " << std::endl;
+    std::cerr << " inputImage , outputImage , prettyInput , prettyOutput , redChannel , greenChannel , nirChannel ," << std::endl;
     return 1;
   }
 
@@ -102,7 +98,7 @@ int main(int argc, char* argv[])
   // The AVI (Angular Vegetation Index) is
   // instantiated using the image pixel types as template parameters.
 
-  typedef otb::Functor::AVI<InputPixelType, InputPixelType, InputPixelType, OutputPixelType> FunctorType;
+  typedef otb::Functor::AVI<InputPixelType, OutputPixelType> FunctorType;
 
   // The
   // \doxygen{itk}{UnaryFunctorImageFilter}
@@ -123,18 +119,9 @@ int main(int argc, char* argv[])
 
   // The three used index bands (red, green and NIR) are declared.
 
-  filter->GetFunctor().SetRedIndex(::atoi(argv[5]));
-  filter->GetFunctor().SetGreenIndex(::atoi(argv[6]));
-  filter->GetFunctor().SetNIRIndex(::atoi(argv[7]));
-
-  // The $\lambda$ R, G and NIR parameters are set. The
-  // \doxygen{otb::Functor}{AVI}
-  // class sets the default values of $\lambda$ to $660$, $560$ and
-  // $830$.
-
-  filter->GetFunctor().SetLambdaR(::atof(argv[8]));
-  filter->GetFunctor().SetLambdaG(::atof(argv[9]));
-  filter->GetFunctor().SetLambdaNir(::atof(argv[10]));
+  filter->GetFunctor().SetBandIndex(CommonBandNames::RED, ::atoi(argv[5]));
+  filter->GetFunctor().SetBandIndex(CommonBandNames::GREEN, ::atoi(argv[6]));
+  filter->GetFunctor().SetBandIndex(CommonBandNames::NIR, ::atoi(argv[7]));
 
   // The filter input is linked to the reader output and
   // the filter output is linked to the writer input.
diff --git a/Examples/Radiometry/CMakeLists.txt b/Examples/Radiometry/CMakeLists.txt
index c33eed313333af0b2a675d1ef787e596ebe4a3d0..48af7c1d8fb654efb198553c3b648b932073c558 100644
--- a/Examples/Radiometry/CMakeLists.txt
+++ b/Examples/Radiometry/CMakeLists.txt
@@ -32,10 +32,6 @@ endif()
 add_executable(AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter.cxx)
 target_link_libraries(AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter ${OTB_LIBRARIES})
 
-add_executable(NDVIRAndNIRVegetationIndexImageFilter NDVIRAndNIRVegetationIndexImageFilter.cxx)
-target_link_libraries(NDVIRAndNIRVegetationIndexImageFilter ${OTB_LIBRARIES})
-
-
 if(BUILD_TESTING)
   add_subdirectory(test)
 endif()
diff --git a/Examples/Radiometry/NDVIRAndNIRVegetationIndexImageFilter.cxx b/Examples/Radiometry/NDVIRAndNIRVegetationIndexImageFilter.cxx
deleted file mode 100644
index d95aae28065d70a953109885e2f428fddc14f5a9..0000000000000000000000000000000000000000
--- a/Examples/Radiometry/NDVIRAndNIRVegetationIndexImageFilter.cxx
+++ /dev/null
@@ -1,207 +0,0 @@
-/*
- * Copyright (C) 2005-2019 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.
- */
-
-
-/* Example usage:
-./NDVIRAndNIRVegetationIndexImageFilter Input/NDVI_2.hdr \
-                                        Input/NDVI_3.hdr \
-                                        Output/NDVIRAndNIRVegetationIndex.tif \
-                                        Output/pretty_Red.png \
-                                        Output/pretty_NIR.png \
-                                        Output/pretty_NDVIRAndNIRVegetationIndex.png
-*/
-
-
-// \index{otb::VegetationIndicesFunctor}
-// \index{otb::VegetationIndicesFunctor!header}
-//
-// The following example illustrates the use of the
-// \doxygen{itk}{BinaryFunctorImageFilter} with the use of the Normalized
-// Difference Vegatation Index (NDVI).
-// NDVI computes the difference between the NIR channel, noted $L_{NIR}$, and the red channel,
-// noted $L_{r}$ radiances reflected from the surface and transmitted through the atmosphere:
-//
-// \begin{equation}
-// \mathbf{NDVI} = \frac{L_{NIR}-L_{r}}{L_{NIR}+L_{r}}
-// \end{equation}
-//
-//  \relatedClasses
-//  \begin{itemize}
-//  \item \subdoxygen{otb}{Functor}{RVI}
-//  \item \subdoxygen{otb}{Functor}{PVI}
-//  \item \subdoxygen{otb}{Functor}{SAVI}
-//  \item \subdoxygen{otb}{Functor}{TSAVI}
-//  \item \subdoxygen{otb}{Functor}{MSAVI}
-//  \item \subdoxygen{otb}{Functor}{GEMI}
-//  \item \subdoxygen{otb}{Functor}{WDVI}
-//  \item \subdoxygen{otb}{Functor}{IPVI}
-//  \item \subdoxygen{otb}{Functor}{TNDVI}
-//  \end{itemize}
-//
-// Let's look at the minimal code required to use this algorithm.
-
-#include "itkMacro.h"
-#include "otbImage.h"
-#include "otbImageFileReader.h"
-#include "otbImageFileWriter.h"
-#include "itkBinaryFunctorImageFilter.h"
-#include "otbVegetationIndicesFunctor.h"
-#include "itkRescaleIntensityImageFilter.h"
-
-int main(int argc, char* argv[])
-{
-  if (argc < 6)
-  {
-    std::cerr << "Missing Parameters " << std::endl;
-    std::cerr << "Usage: " << argv[0];
-    std::cerr << " inputImage1 , inputImage2 , outputImage , prettyinputImage1 , prettyinputImage2 , prettyOutput" << std::endl;
-    return 1;
-  }
-
-  // The image types are now defined using pixel types the
-  // dimension. Input and output images are defined as \doxygen{otb}{Image}.
-
-  const unsigned int                             Dimension = 2;
-  typedef double                                 InputPixelType;
-  typedef float                                  OutputPixelType;
-  typedef otb::Image<InputPixelType, Dimension>  InputRImageType;
-  typedef otb::Image<InputPixelType, Dimension>  InputNIRImageType;
-  typedef otb::Image<OutputPixelType, Dimension> OutputImageType;
-
-  // We instantiate reader and writer types
-  typedef otb::ImageFileReader<InputRImageType>   RReaderType;
-  typedef otb::ImageFileReader<InputNIRImageType> NIRReaderType;
-  typedef otb::ImageFileWriter<OutputImageType>   WriterType;
-
-  // The NDVI (Normalized Difference Vegetation Index) is instantiated using
-  // the images pixel type as template parameters. It is
-  // implemented as a functor class which will be passed as a
-  // parameter to an \doxygen{itk}{BinaryFunctorImageFilter}.
-
-  typedef otb::Functor::NDVI<InputPixelType, InputPixelType, OutputPixelType> FunctorType;
-
-  // The \doxygen{itk}{BinaryFunctorImageFilter} type is instantiated using the images
-  // types and the NDVI functor as template parameters.
-
-  typedef itk::BinaryFunctorImageFilter<InputRImageType, InputNIRImageType, OutputImageType, FunctorType> NDVIImageFilterType;
-
-  // Instantiating object
-  NDVIImageFilterType::Pointer filter            = NDVIImageFilterType::New();
-  RReaderType::Pointer                 readerR   = RReaderType::New();
-  NIRReaderType::Pointer               readerNIR = NIRReaderType::New();
-  WriterType::Pointer                  writer    = WriterType::New();
-
-  //  Now the input images are set and a name is given to the output image.
-
-  readerR->SetFileName(argv[1]);
-  readerNIR->SetFileName(argv[2]);
-  writer->SetFileName(argv[3]);
-
-  // We set the processing pipeline: filter inputs are linked to
-  // the reader output and the filter output is linked to the writer
-  // input.
-
-  filter->SetInput1(readerR->GetOutput());
-  filter->SetInput2(readerNIR->GetOutput());
-
-  writer->SetInput(filter->GetOutput());
-
-  // Invocation of the \code{Update()} method on the writer triggers the
-  // execution of the pipeline.  It is recommended to place \code{update()} calls in a
-  // \code{try/catch} block in case errors occur and exceptions are thrown.
-
-  try
-  {
-    writer->Update();
-  }
-  catch (itk::ExceptionObject& excep)
-  {
-    std::cerr << "Exception caught !" << std::endl;
-    std::cerr << excep << std::endl;
-  }
-  catch (...)
-  {
-    std::cout << "Unknown exception !" << std::endl;
-    return EXIT_FAILURE;
-  }
-
-  // Pretty image creation for the printing
-  typedef otb::Image<unsigned char, Dimension>                                       OutputPrettyImageType;
-  typedef otb::ImageFileWriter<OutputPrettyImageType>                                WriterPrettyType;
-  typedef itk::RescaleIntensityImageFilter<OutputImageType, OutputPrettyImageType>   RescalerType;
-  typedef itk::RescaleIntensityImageFilter<InputRImageType, OutputPrettyImageType>   RescalerRType;
-  typedef itk::RescaleIntensityImageFilter<InputNIRImageType, OutputPrettyImageType> RescalerNIRType;
-
-  RescalerType::Pointer     rescaler     = RescalerType::New();
-  WriterPrettyType::Pointer prettyWriter = WriterPrettyType::New();
-  rescaler->SetInput(filter->GetOutput());
-  rescaler->SetOutputMinimum(0);
-  rescaler->SetOutputMaximum(255);
-  prettyWriter->SetFileName(argv[6]);
-  prettyWriter->SetInput(rescaler->GetOutput());
-
-  RescalerRType::Pointer    rescalerR       = RescalerRType::New();
-  RescalerNIRType::Pointer  rescalerNIR     = RescalerNIRType::New();
-  WriterPrettyType::Pointer prettyWriterR   = WriterPrettyType::New();
-  WriterPrettyType::Pointer prettyWriterNIR = WriterPrettyType::New();
-  rescalerR->SetInput(readerR->GetOutput());
-  rescalerR->SetOutputMinimum(0);
-  rescalerR->SetOutputMaximum(255);
-  prettyWriterR->SetFileName(argv[4]);
-  prettyWriterR->SetInput(rescalerR->GetOutput());
-
-  rescalerNIR->SetInput(readerNIR->GetOutput());
-  rescalerNIR->SetOutputMinimum(0);
-  rescalerNIR->SetOutputMaximum(255);
-  prettyWriterNIR->SetFileName(argv[5]);
-  prettyWriterNIR->SetInput(rescalerNIR->GetOutput());
-
-  try
-  {
-    prettyWriter->Update();
-    prettyWriterNIR->Update();
-    prettyWriterR->Update();
-  }
-  catch (itk::ExceptionObject& excep)
-  {
-    std::cerr << "Exception caught !" << std::endl;
-    std::cerr << excep << std::endl;
-  }
-  catch (...)
-  {
-    std::cout << "Unknown exception !" << std::endl;
-    return EXIT_FAILURE;
-  }
-
-  // Let's now run this example using as input the images
-  // \code{NDVI\_3.hdr} and  \code{NDVI\_4.hdr} (images kindly and free of charge given by SISA and CNES)
-  // provided in the directory \code{Examples/Data}.
-  //
-  //
-  // \begin{figure} \center
-  // \includegraphics[width=0.24\textwidth]{pretty_Red.eps}
-  // \includegraphics[width=0.24\textwidth]{pretty_NIR.eps}
-  // \includegraphics[width=0.24\textwidth]{pretty_NDVIRAndNIRVegetationIndex.eps}
-  // \itkcaption[ARVI Example]{NDVI input images on the left (Red channel and NIR channel), on the right the result of the algorithm.}
-  // \label{fig:NDVIRAndNIRIndex}
-  // \end{figure}
-
-  return EXIT_SUCCESS;
-}
diff --git a/Examples/Radiometry/test/CMakeLists.txt b/Examples/Radiometry/test/CMakeLists.txt
index b30d3aa7ce82dc01e623f5bbd7efb4d843133498..4bd6befaff4c0b1de63209fcbf8f48bd55990365 100644
--- a/Examples/Radiometry/test/CMakeLists.txt
+++ b/Examples/Radiometry/test/CMakeLists.txt
@@ -35,7 +35,6 @@ otb_add_test(NAME raTeARVIMultiChannelRAndBAndNIRVegetationIndexImageFilterTest
     1
     2
     3
-    0.6 # Gamma parameter
 )
 
 # -------            AVIMultiChannelRAndGAndNIRVegetationIndexImageFilterTest   ------------------------------
@@ -52,24 +51,6 @@ otb_add_test(NAME raTeAVIMultiChannelRAndGAndNIRVegetationIndexImageFilterTest C
     3
     2
     4 # indices of the channels
-    660.
-    560.
-    830. # lambdaR, lambdaG, lambdaNir
-)
-
-# -------            NDVIRAndNIRVegetationIndexImageFilter   ------------------------------
-
-otb_add_test(NAME raTeNDVIRAndNIRVegetationIndexImageFilterTest COMMAND ${OTB_TEST_DRIVER}
-  --compare-image ${NOTOL}
-    ${BASELINE}/NDVIRAndNIRVegetationIndex.tif
-    ${TEMP}/NDVIRAndNIRVegetationIndex.tif
-  Execute $<TARGET_FILE:NDVIRAndNIRVegetationIndexImageFilter>
-    ${INPUTDATA}/poupees_sub_c1.png
-    ${INPUTDATA}/poupees_sub_c2.png
-    ${TEMP}/NDVIRAndNIRVegetationIndex.tif
-    ${TEMP}/NDVIRAndNIRVegetationIndex2.tif
-    ${TEMP}/NDVIRAndNIRVegetationIndex3.tif
-    ${TEMP}/NDVIRAndNIRVegetationIndex4.tif
 )
 
 if(OTBOpticalCalibration_LOADED)
@@ -103,4 +84,4 @@ otb_add_test(NAME raTeAtmosphericCorrectionSequencementTest COMMAND ${OTB_TEST_D
     2       # Radius;
     0.020   # pixel spacing in kilometers
 )
-endif()
\ No newline at end of file
+endif()
diff --git a/Examples/Simulation/LAIFromNDVIImageTransform.cxx b/Examples/Simulation/LAIFromNDVIImageTransform.cxx
index ec6dce3b1ce40748987273f04334b077785945a1..133e5dd6444ddfd50625cc421e9642dc5558bd98 100644
--- a/Examples/Simulation/LAIFromNDVIImageTransform.cxx
+++ b/Examples/Simulation/LAIFromNDVIImageTransform.cxx
@@ -64,8 +64,7 @@ int main(int argc, char* argv[])
   // Filter type is a generic \doxygen{itk}{UnaryFunctorImageFilter} using Formosat2 specific LAI
   //  \doxygen{otb}{LAIFromNDVIFormosat2Functor}.
 
-  typedef otb::Functor::LAIFromNDVIFormosat2Functor<InputImageType::InternalPixelType, InputImageType::InternalPixelType, OutputImageType::PixelType>
-                                                                                                 FunctorType;
+  typedef otb::Functor::LAIFromNDVIFormosat2Functor<InputImageType::InternalPixelType, OutputImageType::PixelType> FunctorType;
   typedef itk::UnaryFunctorImageFilter<InputImageType, OutputImageType, FunctorType> LAIFRomNDVIImageFilterType;
 
   // Instantiating object
@@ -96,8 +95,8 @@ int main(int argc, char* argv[])
   //
   unsigned int redChannel = static_cast<unsigned int>(atoi(argv[5]));
   unsigned int nirChannel = static_cast<unsigned int>(atoi(argv[6]));
-  filter->GetFunctor().SetRedIndex(redChannel);
-  filter->GetFunctor().SetNIRIndex(nirChannel);
+  filter->GetFunctor().SetBandIndex(CommonBandNames::RED, redChannel);
+  filter->GetFunctor().SetBandIndex(CommonBandNames::NIR, nirChannel);
 
   //  The invocation of the \code{Update()} method triggers the
   //  execution of the pipeline.
diff --git a/Modules/Applications/AppIndices/app/otbRadiometricIndices.cxx b/Modules/Applications/AppIndices/app/otbRadiometricIndices.cxx
index 4d13ff7279378f64dcc20f7e5175b5e005900751..082106fd43fca3df6e443e7cea0009cfc873cd8e 100644
--- a/Modules/Applications/AppIndices/app/otbRadiometricIndices.cxx
+++ b/Modules/Applications/AppIndices/app/otbRadiometricIndices.cxx
@@ -20,16 +20,14 @@
 
 #include "otbWrapperApplication.h"
 #include "otbWrapperApplicationFactory.h"
+#include "otbWrapperNumericalParameter.h"
 
-#include "itkUnaryFunctorImageFilter.h"
+#include "otbVegetationIndicesFunctor.h"
 #include "otbWaterIndicesFunctor.h"
 #include "otbBuiltUpIndicesFunctor.h"
 #include "otbSoilIndicesFunctor.h"
-
-#include "otbImageList.h"
-#include "otbImageListToVectorImageFilter.h"
-
-#include "otbWrapperNumericalParameter.h"
+#include "otbIndicesStackFunctor.h"
+#include "otbFunctorImageFilter.h"
 
 namespace otb
 {
@@ -50,78 +48,21 @@ public:
 
   itkTypeMacro(RadiometricIndices, otb::Wrapper::Application);
 
-  /** Output  containers typedef */
-  typedef ObjectList<itk::ProcessObject>                                    FilterListType;
-  typedef ImageList<FloatImageType>                                         ImageListType;
-  typedef ImageListToVectorImageFilter<ImageListType, FloatVectorImageType> ImageListToVectorImageFilterType;
-
-  /** Radiometric water indices functors typedef */
-  typedef Functor::SRWI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType>  SRWIFunctorType;
-  typedef Functor::NDWI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType>  NDWIFunctorType;
-  typedef Functor::NDWI2<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> NDWI2FunctorType;
-  typedef Functor::MNDWI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> MNDWIFunctorType;
-  typedef Functor::NDPI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType>  NDPIFunctorType;
-  typedef Functor::NDTI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType>  NDTIFunctorType;
-  //typedef Functor::WaterSqrtSpectralAngleFunctor<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> WaterSqrtSpectralAngleFunctor;
-
-  /** Radiometric vegetation indices functors typedef */
-  typedef Functor::NDVI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> NDVIFunctor;
-  typedef Functor::TNDVI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> TNDVIFunctor;
-  typedef Functor::RVI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> RVIFunctor;
-  typedef Functor::SAVI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> SAVIFunctor;
-  typedef Functor::TSAVI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> TSAVIFunctor;
-  typedef Functor::MSAVI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> MSAVIFunctor;
-  typedef Functor::MSAVI2<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> MSAVI2Functor;
-  typedef Functor::GEMI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> GEMIFunctor;
-  typedef Functor::IPVI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> IPVIFunctor;
-  typedef Functor::LAIFromNDVILogarithmic<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> LAIFromNDVILogFunctor;
-  typedef Functor::LAIFromReflectancesLinear<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> LAIFromReflLinearFunctor;
-  typedef Functor::LAIFromNDVIFormosat2Functor<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> LAIFromNDVIFormoFunctor;
-
-  /** Radiometric soil indices functors typedef */
-  typedef Functor::IR<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> IRFunctor;
-  typedef Functor::IC<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> ICFunctor;
-  typedef Functor::IB<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> IBFunctor;
-  typedef Functor::IB2<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType, FloatImageType::PixelType> IB2Functor;
-
-  /** Radiometric built up indices functors typedef */
-  typedef Functor::NDBI<FloatVectorImageType::InternalPixelType, FloatVectorImageType::InternalPixelType, FloatImageType::PixelType> NDBIFunctor;
-
-  /** Radiometric indices filters typedef */
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, NDWIFunctorType>  NDWIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, NDWI2FunctorType> NDWI2FilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, MNDWIFunctorType> MNDWIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, NDPIFunctorType>  NDPIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, NDTIFunctorType>  NDTIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, SRWIFunctorType>  SRWIFilterType;
-
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, NDVIFunctor>              NDVIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, TNDVIFunctor>             TNDVIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, RVIFunctor>               RVIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, SAVIFunctor>              SAVIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, TSAVIFunctor>             TSAVIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, MSAVIFunctor>             MSAVIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, MSAVI2Functor>            MSAVI2FilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, GEMIFunctor>              GEMIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, IPVIFunctor>              IPVIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, LAIFromNDVILogFunctor>    LAIFromNDVILogFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, LAIFromReflLinearFunctor> LAIFromReflLinearFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, LAIFromNDVIFormoFunctor>  LAIFromNDVIFormoFilterType;
-
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, IRFunctor>                  RIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, ICFunctor>                  CIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, IBFunctor>                  BIFilterType;
-  typedef itk::UnaryFunctorImageFilter<FloatVectorImageType, FloatImageType, IB2Functor>                 BI2FilterType;
-
-  struct indiceSpec
+  using InputType  = FloatVectorImageType::InternalPixelType;
+  using OutputType = FloatImageType::PixelType;
+
+  using RadiometricIndexType    = otb::Functor::RadiometricIndex<InputType, OutputType>;
+  using IndicesStackFunctorType = otb::Functor::IndicesStackFunctor<RadiometricIndexType>;
+
+  class indiceSpec
   {
+  public:
+    indiceSpec(std::string k, std::string i, RadiometricIndexType* ind) : key(k), item(i), indice(ind)
+    {
+    }
     std::string key;
     std::string item;
-    std::string description;
-    std::string type;
-    std::string chan1;
-    std::string chan2;
-    std::string chan3;
+    std::unique_ptr<RadiometricIndexType> indice;
   };
 
 
@@ -153,46 +94,52 @@ private:
     AddParameter(ParameterType_Int,  "channels.blue",  "Blue Channel");
     SetParameterDescription("channels.blue", "Blue channel index");
     SetDefaultParameterInt("channels.blue", 1);
+    SetMinimumParameterIntValue("channels.blue", 1);
+
     AddParameter(ParameterType_Int,  "channels.green",  "Green Channel");
     SetParameterDescription("channels.green", "Green channel index");
     SetDefaultParameterInt("channels.green", 1);
+    SetMinimumParameterIntValue("channels.green", 1);
+
     AddParameter(ParameterType_Int,  "channels.red",  "Red Channel");
     SetParameterDescription("channels.red", "Red channel index");
     SetDefaultParameterInt("channels.red", 1);
+    SetMinimumParameterIntValue("channels.red", 1);
+
     AddParameter(ParameterType_Int,  "channels.nir",  "NIR Channel");
     SetParameterDescription("channels.nir", "NIR channel index");
     SetDefaultParameterInt("channels.nir", 1);
+    SetMinimumParameterIntValue("channels.nir", 1);
+
     AddParameter(ParameterType_Int,  "channels.mir",  "Mir Channel");
     SetParameterDescription("channels.mir", "Mir channel index");
     SetDefaultParameterInt("channels.mir", 1);
-    //AddParameter(ParameterType_Int,  "channels.rho860",  "Rho860 Channel");
-    //SetParameterDescription("channels.rho860", "860nm band channel index");
-    //SetDefaultParameterInt("channels.rho860", 1);
-    //AddParameter(ParameterType_Int,  "channels.rho1240",  "Rho1240 Channel");
-    //SetParameterDescription("channels.rho1240", "1240nm band channel index");
-    //SetDefaultParameterInt("channels.rho1240", 1);
+    SetMinimumParameterIntValue("channels.mir", 1);
 
     AddParameter(ParameterType_ListView,  "list", "Available Radiometric Indices");
     SetParameterDescription("list",
-        "List of available radiometric indices with their relevant channels in brackets:\n\n"
-        "* Vegetation:NDVI - Normalized difference vegetation index (Red, NIR)\n"
-        "* Vegetation:TNDVI - Transformed normalized difference vegetation index (Red, NIR)\n"
-        "* Vegetation:RVI - Ratio vegetation index (Red, NIR)\n"
-        "* Vegetation:SAVI - Soil adjusted vegetation index (Red, NIR)\n"
-        "* Vegetation:TSAVI - Transformed soil adjusted vegetation index (Red, NIR)\n"
-        "* Vegetation:MSAVI - Modified soil adjusted vegetation index (Red, NIR)\n"
-        "* Vegetation:MSAVI2 - Modified soil adjusted vegetation index 2 (Red, NIR)\n"
-        "* Vegetation:GEMI - Global environment monitoring index (Red, NIR)\n"
-        "* Vegetation:IPVI - Infrared percentage vegetation index (Red, NIR)\n"
-        "* Water:NDWI - Normalized difference water index (Gao 1996) (NIR, MIR)\n"
-        "* Water:NDWI2 - Normalized difference water index (Mc Feeters 1996) (Green, NIR)\n"
-        "* Water:MNDWI - Modified normalized difference water index (Xu 2006) (Green, MIR)\n"
-        "* Water:NDPI - Normalized difference pond index (Lacaux et al.) (MIR, Green)\n"
-        "* Water:NDTI - Normalized difference turbidity index (Lacaux et al.) (Red, Green)\n"
-        "* Soil:RI - Redness index (Red, Green)\n"
-        "* Soil:CI - Color index (Red, Green)\n"
-        "* Soil:BI - Brightness index (Red, Green)\n"
-        "* Soil:BI2 - Brightness index 2 (NIR, Red, Green)");
+                            "List of available radiometric indices with their relevant channels in brackets:\n\n"
+                            "* Vegetation:NDVI - Normalized difference vegetation index (Red, NIR)\n"
+                            "* Vegetation:TNDVI - Transformed normalized difference vegetation index (Red, NIR)\n"
+                            "* Vegetation:RVI - Ratio vegetation index (Red, NIR)\n"
+                            "* Vegetation:SAVI - Soil adjusted vegetation index (Red, NIR)\n"
+                            "* Vegetation:TSAVI - Transformed soil adjusted vegetation index (Red, NIR)\n"
+                            "* Vegetation:MSAVI - Modified soil adjusted vegetation index (Red, NIR)\n"
+                            "* Vegetation:MSAVI2 - Modified soil adjusted vegetation index 2 (Red, NIR)\n"
+                            "* Vegetation:GEMI - Global environment monitoring index (Red, NIR)\n"
+                            "* Vegetation:IPVI - Infrared percentage vegetation index (Red, NIR)\n"
+                            "* Vegetation:LAIFromNDVILog - Leaf Area Index from log NDVI (Red, NIR)\n"
+                            "* Vegetation::LAIFromReflLinear - Leaf Area Index from reflectances with linear combination (Red, NIR)\n"
+                            "* Vegetation::LAIFromNDVIFormo - Leaf Area Index from Formosat 2  TOC (Red, NIR)\n"
+                            "* Water:NDWI - Normalized difference water index (Gao 1996) (NIR, MIR)\n"
+                            "* Water:NDWI2 - Normalized difference water index (Mc Feeters 1996) (Green, NIR)\n"
+                            "* Water:MNDWI - Modified normalized difference water index (Xu 2006) (Green, MIR)\n"
+                            "* Water:NDTI - Normalized difference turbidity index (Lacaux et al.) (Red, Green)\n"
+                            "* Soil:RI - Redness index (Red, Green)\n"
+                            "* Soil:CI - Color index (Red, Green)\n"
+                            "* Soil:BI - Brightness index (Red, Green)\n"
+                            "* Soil:BI2 - Brightness index 2 (NIR, Red, Green)\n"
+                            "* BuiltUp:ISU - Built Surfaces Index (NIR,Red) ");
 
     AddRAMParameter();
 
@@ -205,392 +152,118 @@ private:
 
     m_Map.clear();
 
-    // Add Available choices
-    indiceSpec s_NDVI;
-    s_NDVI.key   = "list.ndvi";
-    s_NDVI.item  = "Vegetation:NDVI";
-    s_NDVI.description = "";
-    s_NDVI.type  = "NDVI";
-    s_NDVI.chan1 = "red";
-    s_NDVI.chan2 = "nir";
-    s_NDVI.chan3 = "";
-    m_Map.push_back(s_NDVI);
-
-    indiceSpec s_TNDVI;
-    s_TNDVI.key   = "list.tndvi";
-    s_TNDVI.item  = "Vegetation:TNDVI";
-    s_TNDVI.description = "";
-    s_TNDVI.type  = "TNDVI";
-    s_TNDVI.chan1 = "red";
-    s_TNDVI.chan2 = "nir";
-    s_TNDVI.chan3 = "";
-    m_Map.push_back(s_TNDVI);
-
-    indiceSpec s_RVI;
-    s_RVI.key   = "list.rvi";
-    s_RVI.item  = "Vegetation:RVI";
-    s_RVI.description = "";
-    s_RVI.type  = "RVI";
-    s_RVI.chan1 = "red";
-    s_RVI.chan2 = "nir";
-    s_RVI.chan3 = "";
-    m_Map.push_back(s_RVI);
-
-    indiceSpec s_SAVI;
-    s_SAVI.key   = "list.savi";
-    s_SAVI.item  = "Vegetation:SAVI";
-    s_SAVI.description = "";
-    s_SAVI.type  = "SAVI";
-    s_SAVI.chan1 = "red";
-    s_SAVI.chan2 = "nir";
-    s_SAVI.chan3 = "";
-    m_Map.push_back(s_SAVI);
-
-    indiceSpec s_TSAVI;
-    s_TSAVI.key   = "list.tsavi";
-    s_TSAVI.item  = "Vegetation:TSAVI";
-    s_TSAVI.description = "";
-    s_TSAVI.type  = "TSAVI";
-    s_TSAVI.chan1 = "red";
-    s_TSAVI.chan2 = "nir";
-    s_TSAVI.chan3 = "";
-    m_Map.push_back(s_TSAVI);
-
-    indiceSpec s_MSAVI;
-    s_MSAVI.key   = "list.msavi";
-    s_MSAVI.item  = "Vegetation:MSAVI";
-    s_MSAVI.description = "";
-    s_MSAVI.type  = "MSAVI";
-    s_MSAVI.chan1 = "red";
-    s_MSAVI.chan2 = "nir";
-    s_MSAVI.chan3 = "";
-    m_Map.push_back(s_MSAVI);
-
-    indiceSpec s_MSAVI2;
-    s_MSAVI2.key   = "list.msavi2";
-    s_MSAVI2.item  = "Vegetation:MSAVI2";
-    s_MSAVI2.description = "";
-    s_MSAVI2.type  = "MSAVI2";
-    s_MSAVI2.chan1 = "red";
-    s_MSAVI2.chan2 = "nir";
-    s_MSAVI2.chan3 = "";
-    m_Map.push_back(s_MSAVI2);
-
-    indiceSpec s_GEMI;
-    s_GEMI.key   = "list.gemi";
-    s_GEMI.item  = "Vegetation:GEMI";
-    s_GEMI.description = "";
-    s_GEMI.type  = "GEMI";
-    s_GEMI.chan1 = "red";
-    s_GEMI.chan2 = "nir";
-    s_GEMI.chan3 = "";
-    m_Map.push_back(s_GEMI);
-
-    indiceSpec s_IPVI;
-    s_IPVI.key   = "list.ipvi";
-    s_IPVI.item  = "Vegetation:IPVI";
-    s_IPVI.description = "";
-    s_IPVI.type  = "IPVI";
-    s_IPVI.chan1 = "red";
-    s_IPVI.chan2 = "nir";
-    s_IPVI.chan3 = "";
-    m_Map.push_back(s_IPVI);
-
-    indiceSpec s_LAIFromNDVILog;
-    s_LAIFromNDVILog.key   = "list.laindvilog";
-    s_LAIFromNDVILog.item  = "Vegetation:LAIFromNDVILog";
-    s_LAIFromNDVILog.description = "";
-    s_LAIFromNDVILog.type  = "LAIFromNDVILog";
-    s_LAIFromNDVILog.chan1 = "red";
-    s_LAIFromNDVILog.chan2 = "nir";
-    s_LAIFromNDVILog.chan3 = "";
-    m_Map.push_back(s_LAIFromNDVILog);
-
-    indiceSpec s_LAIFromReflLinear;
-    s_LAIFromReflLinear.key   = "list.lairefl";
-    s_LAIFromReflLinear.item  = "Vegetation:LAIFromReflLinear";
-    s_LAIFromReflLinear.description = "";
-    s_LAIFromReflLinear.type  = "LAIFromReflLinear";
-    s_LAIFromReflLinear.chan1 = "red";
-    s_LAIFromReflLinear.chan2 = "nir";
-    s_LAIFromReflLinear.chan3 = "";
-    m_Map.push_back(s_LAIFromReflLinear);
-
-    indiceSpec s_LAIFromNDVIFormo;
-    s_LAIFromNDVIFormo.key   = "list.laindviformo";
-    s_LAIFromNDVIFormo.item  = "Vegetation:LAIFromNDVIFormo";
-    s_LAIFromNDVIFormo.description = "";
-    s_LAIFromNDVIFormo.type  = "LAIFromNDVIFormo";
-    s_LAIFromNDVIFormo.chan1 = "red";
-    s_LAIFromNDVIFormo.chan2 = "nir";
-    s_LAIFromNDVIFormo.chan3 = "";
-    m_Map.push_back(s_LAIFromNDVIFormo);
-
-    indiceSpec s_NDWI;
-    s_NDWI.key   = "list.ndwi";
-    s_NDWI.item  = "Water:NDWI";
-    s_NDWI.description = "";
-    s_NDWI.type  = "NDWI";
-    s_NDWI.chan1 = "nir";
-    s_NDWI.chan2 = "mir";
-    s_NDWI.chan3 = "";
-    m_Map.push_back(s_NDWI);
-
-    indiceSpec s_NDWI2;
-    s_NDWI2.key   = "list.ndwi2";
-    s_NDWI2.item  = "Water:NDWI2";
-    s_NDWI2.description = "";
-    s_NDWI2.type  = "NDWI2";
-    s_NDWI2.chan1 = "green";
-    s_NDWI2.chan2 = "nir";
-    s_NDWI2.chan3 = "";
-    m_Map.push_back(s_NDWI2);
-
-    indiceSpec s_MNDWI;
-    s_MNDWI.key   = "list.mndwi";
-    s_MNDWI.item  = "Water:MNDWI";
-    s_MNDWI.description = "";
-    s_MNDWI.type  = "MNDWI";
-    s_MNDWI.chan1 = "green";
-    s_MNDWI.chan2 = "mir";
-    s_MNDWI.chan3 = "";
-    m_Map.push_back(s_MNDWI);
-
-    indiceSpec s_NDPI;
-    s_NDPI.key   = "list.ndpi";
-    s_NDPI.item  = "Water:NDPI";
-    s_NDPI.description = "";
-    s_NDPI.type  = "NDPI";
-    s_NDPI.chan1 = "mir";
-    s_NDPI.chan2 = "green";
-    s_NDPI.chan3 = "";
-    m_Map.push_back(s_NDPI);
-
-    indiceSpec s_NDTI;
-    s_NDTI.key   = "list.ndti";
-    s_NDTI.item  = "Water:NDTI";
-    s_NDTI.description = "";
-    s_NDTI.type  = "NDTI";
-    s_NDTI.chan1 = "red";
-    s_NDTI.chan2 = "green";
-    s_NDTI.chan3 = "";
-    m_Map.push_back(s_NDTI);
-
-    indiceSpec s_SRWI;
-    s_SRWI.key   = "list.srwi";
-    s_SRWI.item  = "Water:SRWI";
-    s_SRWI.description = "";
-    s_SRWI.type  = "SRWI";
-    s_SRWI.chan1 = "rho860";
-    s_SRWI.chan2 = "rho1240";
-    s_SRWI.chan3 = "";
-    //m_Map.push_back(s_SRWI);
-
-    indiceSpec s_RI;
-    s_RI.key   = "list.ri";
-    s_RI.item  = "Soil:RI";
-    s_RI.description = "";
-    s_RI.type  = "RI";
-    s_RI.chan1 = "red";
-    s_RI.chan2 = "green";
-    s_RI.chan3 = "";
-    m_Map.push_back(s_RI);
-
-    indiceSpec s_CI;
-    s_CI.key   = "list.ci";
-    s_CI.item  = "Soil:CI";
-    s_CI.description = "";
-    s_CI.type  = "CI";
-    s_CI.chan1 = "red";
-    s_CI.chan2 = "green";
-    s_CI.chan3 = "";
-    m_Map.push_back(s_CI);
-
-    indiceSpec s_BI;
-    s_BI.key   = "list.bi";
-    s_BI.item  = "Soil:BI";
-    s_BI.description = "";
-    s_BI.type  = "BI";
-    s_BI.chan1 = "red";
-    s_BI.chan2 = "green";
-    s_BI.chan3 = "";
-    m_Map.push_back(s_BI);
-
-    indiceSpec s_BI2;
-    s_BI2.key   = "list.bi2";
-    s_BI2.item  = "Soil:BI2";
-    s_BI2.description = "";
-    s_BI2.type  = "BI2";
-    s_BI2.chan1 = "nir";
-    s_BI2.chan2 = "red";
-    s_BI2.chan3 = "green";
-    m_Map.push_back(s_BI2);
+    m_Map.push_back({"list.ndvi", "Vegetation:NDVI", new otb::Functor::NDVI<InputType, OutputType>()});
+    m_Map.push_back({"list.tndvi", "Vegetation:TNDVI", new otb::Functor::TNDVI<InputType, OutputType>()});
+    m_Map.push_back({"list.rdvi", "Vegetation:RVI", new otb::Functor::RVI<InputType, OutputType>()});
+    m_Map.push_back({"list.savi", "Vegetation:SAVI", new otb::Functor::SAVI<InputType, OutputType>()});
+    m_Map.push_back({"list.tsavi", "Vegetation:TSAVI", new otb::Functor::TSAVI<InputType, OutputType>()});
+    m_Map.push_back({"list.msavi", "Vegetation:MSAVI", new otb::Functor::MSAVI<InputType, OutputType>()});
+    m_Map.push_back({"list.msavi2", "Vegetation:MSAVI2", new otb::Functor::MSAVI2<InputType, OutputType>()});
+    m_Map.push_back({"list.gemi", "Vegetation:GEMI", new otb::Functor::GEMI<InputType, OutputType>()});
+    m_Map.push_back({"list.ipvi", "Vegetation:IPVI", new otb::Functor::IPVI<InputType, OutputType>()});
+    m_Map.push_back({"list.laindvilog", "Vegetation:LAIFromNDVILog", new otb::Functor::LAIFromNDVILogarithmic<InputType, OutputType>()});
+    m_Map.push_back({"list.lairefl", "Vegetation:LAIFromReflLinear", new otb::Functor::LAIFromReflectancesLinear<InputType, OutputType>()});
+    m_Map.push_back({"list.laindviformo", "Vegetation:LAIFromNDVIFormo", new otb::Functor::LAIFromNDVIFormosat2Functor<InputType, OutputType>()});
+    m_Map.push_back({"list.ndwi", "Water:NDWI", new otb::Functor::NDWI<InputType, OutputType>()});
+    m_Map.push_back({"list.ndwi2", "Water:NDWI2", new otb::Functor::NDWI2<InputType, OutputType>()});
+    m_Map.push_back({"list.mndwi", "Water:MNDWI", new otb::Functor::MNDWI<InputType, OutputType>()});
+    m_Map.push_back({"list.ndti", "Water:NDTI", new otb::Functor::NDTI<InputType, OutputType>()});
+    m_Map.push_back({"list.ri", "Soil:RI", new otb::Functor::RI<InputType, OutputType>()});
+    m_Map.push_back({"list.ci", "Soil:CI", new otb::Functor::CI<InputType, OutputType>()});
+    m_Map.push_back({"list.bi", "Soil:BI", new otb::Functor::BI<InputType, OutputType>()});
+    m_Map.push_back({"list.bi2", "Soil:BI2", new otb::Functor::BI2<InputType, OutputType>()});
+    m_Map.push_back({"list.isu", "BuiltUp:ISU", new otb::Functor::ISU<InputType, OutputType>()});
 
     ClearChoices("list");
+
     for ( unsigned int i=0; i<m_Map.size(); i++ )
       {
       AddChoice(m_Map[i].key, m_Map[i].item);
-      //SetParameterDescription(m_Map[i].item, m_Map[i].description);
       }
   }
 
-  void DoUpdateParameters() override
+  // Compute required bands for selected indices
+  std::set<CommonBandNames> GetRequiredBands()
   {
-    //Nothing to do here
-  }
+    std::set<CommonBandNames> required;
 
-#define otbRadiometricWaterIndicesMacro( type )                           \
-    {                                                                     \
-    type##FilterType::Pointer l_##type##Filter = type##FilterType::New(); \
-    std::ostringstream oss;                                               \
-    oss<<"channels."<<m_Map[GetSelectedItems("list")[idx]].chan1;         \
-    l_##type##Filter->GetFunctor().SetIndex1(this->GetParameterInt(oss.str()));\
-    oss.str("");                                                          \
-    oss<<"channels."<<m_Map[GetSelectedItems("list")[idx]].chan2;         \
-    l_##type##Filter->GetFunctor().SetIndex2(this->GetParameterInt(oss.str()));\
-    l_##type##Filter->SetInput(inImage);                                  \
-    m_FilterList->PushBack( l_##type##Filter );                           \
-    m_ImageList->PushBack( l_##type##Filter->GetOutput() );               \
-    otbAppLogINFO(<< m_Map[GetSelectedItems("list")[idx]].item << " added.");\
-    }
-
-#define otbRadiometricVegetationIndicesMacro( type )                      \
-    {                                                                     \
-    type##FilterType::Pointer l_##type##Filter = type##FilterType::New(); \
-    std::ostringstream oss;                                               \
-    oss<<"channels."<<m_Map[GetSelectedItems("list")[idx]].chan1;         \
-    l_##type##Filter->GetFunctor().SetRedIndex(this->GetParameterInt(oss.str())); \
-    oss.str("");                                                          \
-    oss<<"channels."<<m_Map[GetSelectedItems("list")[idx]].chan2;         \
-    l_##type##Filter->GetFunctor().SetNIRIndex(this->GetParameterInt(oss.str())); \
-    l_##type##Filter->SetInput(inImage);                                  \
-    m_FilterList->PushBack( l_##type##Filter );                           \
-    m_ImageList->PushBack( l_##type##Filter->GetOutput() );               \
-    otbAppLogINFO(<<m_Map[GetSelectedItems("list")[idx]].item<<" added.");\
+    for (unsigned int idx = 0; idx < GetSelectedItems("list").size(); ++idx)
+    {
+      auto requiredForCurrentIndice = m_Map[GetSelectedItems("list")[idx]].indice->GetRequiredBands();
+      required.insert(requiredForCurrentIndice.begin(), requiredForCurrentIndice.end());
     }
+    return required;
+  }
 
-#define otbRadiometricSoilIndicesMacro( type )                            \
-    {                                                                     \
-    type##FilterType::Pointer l_##type##Filter = type##FilterType::New(); \
-    std::ostringstream oss;                                               \
-    oss<<"channels."<<m_Map[GetSelectedItems("list")[idx]].chan1;         \
-    l_##type##Filter->GetFunctor().SetRedIndex(this->GetParameterInt(oss.str()));\
-    oss.str("");                                                          \
-    oss<<"channels."<<m_Map[GetSelectedItems("list")[idx]].chan2;         \
-    l_##type##Filter->GetFunctor().SetGreenIndex(this->GetParameterInt(oss.str()));\
-    l_##type##Filter->SetInput(inImage);                                  \
-    m_FilterList->PushBack( l_##type##Filter );                           \
-    m_ImageList->PushBack( l_##type##Filter->GetOutput() );               \
-    otbAppLogINFO(<< m_Map[GetSelectedItems("list")[idx]].item << " added.");\
-    }
+  void DoUpdateParameters() override
+  {
+    //Nothing to do here
+  }
 
   void DoExecute() override
   {
+    // Retrieve number of bands of input image
+    unsigned int nbChan = GetParameterImage("in")->GetNumberOfComponentsPerPixel();
 
-    int nbChan = GetParameterImage("in")->GetNumberOfComponentsPerPixel();
+    // Derive required bands from selected indices
+    auto requiredBands = GetRequiredBands();
 
-    if (   (this->GetParameterInt("channels.blue")  <= nbChan)
-        && (this->GetParameterInt("channels.green") <= nbChan)
-        && (this->GetParameterInt("channels.red")   <= nbChan)
-        && (this->GetParameterInt("channels.nir")   <= nbChan)
-        && (this->GetParameterInt("channels.mir")   <= nbChan))
-      {
+    // Map to store association between bands and indices
+    std::map<CommonBandNames, size_t> bandIndicesMap;
 
-      m_FilterList = FilterListType::New();
-      m_ImageList  = ImageListType::New();
-      m_Concatener = ImageListToVectorImageFilterType::New();
-
-      FloatVectorImageType* inImage = GetParameterImage("in");
+    // Lambda that will:
+    // - Check if band is required,
+    // - Check band index range,
+    // - Populate the bandIndicesMap
+    auto bandChecker = [this, requiredBands, nbChan](std::map<CommonBandNames, size_t>& indicesMap, const CommonBandNames& band, const std::string& key) {
+      if (requiredBands.find(band) != requiredBands.end())
+      {
+        unsigned int idx = this->GetParameterInt(key);
 
-      for (unsigned int idx = 0; idx < GetSelectedItems("list").size(); ++idx)
+        if (idx > nbChan)
         {
-
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:NDVI")
-          otbRadiometricVegetationIndicesMacro(NDVI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:TNDVI")
-          otbRadiometricVegetationIndicesMacro(TNDVI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:RVI")
-          otbRadiometricVegetationIndicesMacro(RVI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:SAVI")
-          otbRadiometricVegetationIndicesMacro(SAVI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:TSAVI")
-          otbRadiometricVegetationIndicesMacro(TSAVI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:MSAVI")
-          otbRadiometricVegetationIndicesMacro(MSAVI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:MSAVI2")
-          otbRadiometricVegetationIndicesMacro(MSAVI2);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:GEMI")
-          otbRadiometricVegetationIndicesMacro(GEMI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:IPVI")
-          otbRadiometricVegetationIndicesMacro(IPVI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:LAIFromNDVILog")
-          otbRadiometricVegetationIndicesMacro(LAIFromNDVILog);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:LAIFromReflLinear")
-          otbRadiometricVegetationIndicesMacro(LAIFromReflLinear);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Vegetation:LAIFromNDVIFormo")
-          otbRadiometricVegetationIndicesMacro(LAIFromNDVIFormo);
-
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Water:NDWI")
-          otbRadiometricWaterIndicesMacro(NDWI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Water:NDWI2")
-          otbRadiometricWaterIndicesMacro(NDWI2);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Water:MNDWI")
-          otbRadiometricWaterIndicesMacro(MNDWI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Water:NDPI")
-          otbRadiometricWaterIndicesMacro(NDPI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Water:NDTI")
-          otbRadiometricWaterIndicesMacro(NDTI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Water:SRWI")
-          otbRadiometricWaterIndicesMacro(SRWI);
-
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Soil:RI")
-          otbRadiometricSoilIndicesMacro(RI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Soil:CI")
-          otbRadiometricSoilIndicesMacro(CI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Soil:BI")
-          otbRadiometricSoilIndicesMacro(BI);
-        if (m_Map[GetSelectedItems("list")[idx]].item == "Soil:BI2")
-          {
-          BI2FilterType::Pointer l_BI2Filter = BI2FilterType::New();
-          std::ostringstream oss;
-          oss<<"channels."<<m_Map[GetSelectedItems("list")[idx]].chan1;
-          l_BI2Filter->GetFunctor().SetNIRIndex(this->GetParameterInt(oss.str()));
-          oss.str("");
-          oss<<"channels."<<m_Map[GetSelectedItems("list")[idx]].chan2;
-          l_BI2Filter->GetFunctor().SetRedIndex(this->GetParameterInt(oss.str()));
-          oss.str("");
-          oss<<"channels."<<m_Map[GetSelectedItems("list")[idx]].chan3;
-          l_BI2Filter->GetFunctor().SetGreenIndex(this->GetParameterInt(oss.str()));
-          l_BI2Filter->SetInput(inImage);
-          m_FilterList->PushBack( l_BI2Filter );
-          m_ImageList->PushBack( l_BI2Filter->GetOutput() );
-          otbAppLogINFO(<< m_Map[GetSelectedItems("list")[idx]].item << " added.");
-          }
-
+          otbAppLogFATAL(<< "Index for band " << key << " exceeds the number of channels in image (" << nbChan << ")");
         }
-
-      if( m_ImageList->Size() == 0 )
+        else
         {
-        itkExceptionMacro(<< "No indices selected...");
+          indicesMap[band] = idx;
         }
+      }
+    };
 
-      m_Concatener->SetInput( m_ImageList );
-      m_Concatener->UpdateOutputInformation();
+    // Call lambda for each possible band
+    bandChecker(bandIndicesMap, CommonBandNames::BLUE, "channels.blue");
+    bandChecker(bandIndicesMap, CommonBandNames::GREEN, "channels.green");
+    bandChecker(bandIndicesMap, CommonBandNames::RED, "channels.red");
+    bandChecker(bandIndicesMap, CommonBandNames::NIR, "channels.nir");
+    bandChecker(bandIndicesMap, CommonBandNames::MIR, "channels.mir");
 
-      SetParameterOutputImage("out", m_Concatener->GetOutput());
-      }
-    else
-      {
-      itkExceptionMacro(<< "At least one needed channel has an invalid index");
+    std::vector<RadiometricIndexType*> indices;
+
+    // Find selected indices
+    for (unsigned int idx = 0; idx < GetSelectedItems("list").size(); ++idx)
+    {
+      // Retrieve the indice instance
+      indices.push_back(m_Map[GetSelectedItems("list")[idx]].indice.get());
+
+      // And set bands using the band map
+      indices.back()->SetBandsIndices(bandIndicesMap);
       }
 
-  }
+      // Build a composite indices functor to compute all indices at
+      // once
+      auto compositeFunctor = IndicesStackFunctorType(indices);
 
-  FilterListType::Pointer                   m_FilterList;
-  ImageListType::Pointer                    m_ImageList;
-  ImageListToVectorImageFilterType::Pointer m_Concatener;
-  std::vector<indiceSpec>                   m_Map;
+      // Build and plug functor filter
+      auto filter = NewFunctorFilter(compositeFunctor);
+      filter->SetInputs(GetParameterImage("in"));
+      SetParameterOutputImage("out", filter->GetOutput());
+
+      // Call register pipeline to allow streaming and garbage collection
+      RegisterPipeline();
+  }
 
+  std::vector<indiceSpec> m_Map;
 };
 
 }
diff --git a/Modules/Applications/AppIndices/otb-module.cmake b/Modules/Applications/AppIndices/otb-module.cmake
index 6133056e60beb13dd6e2f8ed203274ab1ae3a751..3c100b67faa40cd069df1c1ae1064023ddcc74f4 100644
--- a/Modules/Applications/AppIndices/otb-module.cmake
+++ b/Modules/Applications/AppIndices/otb-module.cmake
@@ -24,7 +24,7 @@ otb_module(OTBAppIndices
   DEPENDS
     OTBIndices
     OTBApplicationEngine
-    OTBObjectList
+    OTBFunctor
   TEST_DEPENDS
     OTBTestKernel
     OTBCommandLine
diff --git a/Modules/Detection/UrbanArea/include/otbUrbanAreaDetectionImageFilter.h b/Modules/Detection/UrbanArea/include/otbUrbanAreaDetectionImageFilter.h
index cfee272c5486c76220c4e513ea43b8647e518dbe..c178c9acc56ab65b4118729618177488b95368cd 100644
--- a/Modules/Detection/UrbanArea/include/otbUrbanAreaDetectionImageFilter.h
+++ b/Modules/Detection/UrbanArea/include/otbUrbanAreaDetectionImageFilter.h
@@ -51,19 +51,25 @@ template<class TInput, class TOutput = double>
 class RadiometricNonWaterNonVegetationIndexFunctor
 {
 public:
-  typedef Functor::NDVI<double, double, double>  VegetationFunctorType;
-  typedef Functor::NDWI2<double, double, double> WaterFunctorType;
+  typedef Functor::NDVI<double, double>          VegetationFunctorType;
+  typedef Functor::NDWI2<double, double>         WaterFunctorType;
   typedef TOutput                                ValueType;
 
   VegetationFunctorType GetVegetationFunctor(){ return m_VegetationFunctor; }
   WaterFunctorType GetWaterFunctor(){ return m_WaterFunctor; }
 
-  void SetRedIndex(int id){ m_VegetationFunctor.SetRedIndex(id); }
-  void SetGreenIndex(int id){ m_WaterFunctor.SetGIndex(id); }
+  void SetRedIndex(int id)
+  {
+    m_VegetationFunctor.SetBandIndex(CommonBandNames::RED, id);
+  }
+  void SetGreenIndex(int id)
+  {
+    m_WaterFunctor.SetBandIndex(CommonBandNames::GREEN, id);
+  }
   void SetNIRIndex(int id)
   {
-    m_VegetationFunctor.SetNIRIndex(id);
-    m_WaterFunctor.SetNIRIndex(id);
+    m_VegetationFunctor.SetBandIndex(CommonBandNames::NIR, id);
+    m_WaterFunctor.SetBandIndex(CommonBandNames::NIR, id);
   }
 
   RadiometricNonWaterNonVegetationIndexFunctor(){}
diff --git a/Modules/Filtering/Projection/test/CMakeLists.txt b/Modules/Filtering/Projection/test/CMakeLists.txt
index e2e346b3c1805dd0124080886dde1d0356cdd5ee..adbca8d880be8c086e3e33a8d86c4424ebd1b0e9 100644
--- a/Modules/Filtering/Projection/test/CMakeLists.txt
+++ b/Modules/Filtering/Projection/test/CMakeLists.txt
@@ -188,6 +188,16 @@ RAPIDEYE/level1B/2008-12-25T005918_RE3_1B-NAC_397971_12345_band3.ntf
 SENTINEL1/S1A_S6_SLC__1SSV_20150619T195043/measurement/s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff
 )
 
+set(TOLERANCE_RATIO
+0
+0
+0
+0
+0
+0
+0.0002
+)
+
 set(IMG_TYPE
 "pleiades-1"
 "wv2-1"
@@ -303,6 +313,8 @@ foreach(current_img ${IMG_TEST_ORTHO})
 	list(GET RESOL ${IMGNB} current_resol         )
 	list(GET GRIDSPACING ${IMGNB} current_grid_spacing)
 	list(GET ISCOMPLEX ${IMGNB} current_is_compex)
+  list(GET TOLERANCE_RATIO ${IMGNB} current_tolerance_ratio)
+
 	math(EXPR IMGNB "${IMGNB} + 1")
 
 	set( MODENB 0)
@@ -313,6 +325,7 @@ foreach(current_img ${IMG_TEST_ORTHO})
 	  otb_add_test(NAME prTvOrthoRectification_${current_imgtype}_${current_mode} COMMAND otbProjectionTestDriver
 		  --compare-image ${EPSILON_4}  ${BASELINE}/prTvOrthoRectification_${current_imgtype}_${current_mode}.tif
 		  ${TEMP}/prTvOrthoRectification_${current_imgtype}_${current_mode}.tif
+      --tolerance-ratio ${current_tolerance_ratio}
 		  otbOrthoRectificationFilter
 		  LARGEINPUT{${current_img}?&geom=${INPUTDATA}/${current_geomgcp}.geom}
 		  ${TEMP}/prTvOrthoRectification_${current_imgtype}_${current_mode}.tif
diff --git a/Modules/Core/Metadata/include/otbBandName.h b/Modules/Radiometry/Indices/include/otbBandName.h
similarity index 83%
rename from Modules/Core/Metadata/include/otbBandName.h
rename to Modules/Radiometry/Indices/include/otbBandName.h
index 6bd25a144e5ac7a0506694035e591195a478e92c..118086139976cff411152658b054e26c59e93a44 100644
--- a/Modules/Core/Metadata/include/otbBandName.h
+++ b/Modules/Radiometry/Indices/include/otbBandName.h
@@ -27,11 +27,30 @@ namespace BandName
 {
 
 /**
-* Provides a way to identify bands when passing the parameters
-* to the radiometric functors.*
-*/
-enum BandName {BLUE, GREEN, RED, NIR, MIR};
-enum LandsatTMBandNames {TM1, TM2, TM3, TM4, TM5, TM60, TM61, TM62, TM7};
+ * Provides a way to identify bands when passing the parameters
+ * to the radiometric functors.*
+ */
+enum class CommonBandNames
+{
+  BLUE,
+  GREEN,
+  RED,
+  NIR,
+  MIR,
+  MAX
+};
+enum LandsatTMBandNames
+{
+  TM1,
+  TM2,
+  TM3,
+  TM4,
+  TM5,
+  TM60,
+  TM61,
+  TM62,
+  TM7
+};
 
 // Note for landsat equivalence
 // http://landsat.gsfc.nasa.gov/news/news-archive/sci_0017.html
@@ -57,8 +76,8 @@ enum LandsatTMBandNames {TM1, TM2, TM3, TM4, TM5, TM60, TM61, TM62, TM7};
 // 6              10.40-12.50 microm                Thermal IR
 // 7                2.08-2.35 microm                Mid-IR
 
-}
+} // namespace BandName
 
-}
+} // namespace otb
 
 #endif
diff --git a/Modules/Radiometry/Indices/include/otbBuiltUpIndicesFunctor.h b/Modules/Radiometry/Indices/include/otbBuiltUpIndicesFunctor.h
index ac8fff8d7e364ae00790595caaaeadada28ce211..426548fd5fc40aedf09bed790d64ccb4e2431d08 100644
--- a/Modules/Radiometry/Indices/include/otbBuiltUpIndicesFunctor.h
+++ b/Modules/Radiometry/Indices/include/otbBuiltUpIndicesFunctor.h
@@ -21,161 +21,12 @@
 #ifndef otbBuiltUpIndicesFunctor_h
 #define otbBuiltUpIndicesFunctor_h
 
-#include "otbVegetationIndicesFunctor.h"
-#include <string>
+#include "otbRadiometricIndex.h"
 
 namespace otb
 {
 namespace Functor
 {
-/**
-   * \class TM4AndTM5IndexBase
-   * \brief Base class for TM4 And TM5 channels of Land Sat
-   * (equivalent to Red and NIR of SPOT5)
-   *
-   *  Implement operators for UnaryFunctorImageFilter templated with a
-   *  VectorImage and BinaryFunctorImageFilter templated with single
-   *  images.
-   *  Subclasses should NOT overload operators, they must  re-implement
-   *  the Evaluate() method.
-   *
-   * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template<class TInput1, class TInput2, class TOutput>
-class TM4AndTM5IndexBase
-{
-public:
-  /// Vector pixel type used to support both vector images and multiple
-  /// input images
-  typedef itk::VariableLengthVector<TInput1> InputVectorType;
-
-  //operators !=
-  bool operator !=(const TM4AndTM5IndexBase&) const
-  {
-    return true;
-  }
-  //operator ==
-  bool operator ==(const TM4AndTM5IndexBase& other) const
-  {
-    return !(*this != other);
-  }
-
-  // Operator on vector pixel type
-  inline TOutput operator ()(const InputVectorType& inputVector) const
-  {
-    return this->Evaluate(inputVector[m_TM4Index - 1], static_cast<TInput2>(inputVector[m_TM5Index - 1]));
-  }
-
-  // Binary operator
-  inline TOutput operator ()(const TInput1& tm4, const TInput2& tm5) const
-  {
-    return this->Evaluate(tm4, tm5);
-  }
-  /// Constructor
-  TM4AndTM5IndexBase() : m_TM4Index(4), m_TM5Index(5) {}
-  /// Desctructor
-  virtual ~TM4AndTM5IndexBase() {}
-
-  /// Set TM4 Index
-  void SetIndex1(unsigned int channel)
-  {
-    m_TM4Index = channel;
-  }
-  /// Get TM4 Index
-  unsigned int GetIndex1() const
-  {
-    return m_TM4Index;
-  }
-  /// Set TM5 Index
-  void SetIndex2(unsigned int channel)
-  {
-    m_TM5Index = channel;
-  }
-  /// Get TM5 Index
-  unsigned int GetIndex2() const
-  {
-    return m_TM5Index;
-  }
-
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::RED)
-      {
-      m_TM4Index = channel;
-      }
-    if (band == BandName::NIR)
-      {
-      m_TM5Index = channel;
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::RED)
-      {
-      return m_TM4Index;
-      }
-    if (band == BandName::NIR)
-      {
-      return m_TM5Index;
-      }
-  }
-
-  /** Return the index name */
-  virtual std::string GetName() const = 0;
-
-protected:
-  // This method must be reimplemented in subclasses to actually
-  // compute the index value
-  virtual TOutput Evaluate(const TInput1& tm4, const TInput2& tm5) const = 0;
-
-private:
-  unsigned int m_TM4Index;
-  unsigned int m_TM5Index;
-};
-
-/** \class NDBI
- *  \brief This functor computes the Normalized Difference Built Up Index (NDBI)
- *
- *  [Zha 2003]
- *
- *  \ingroup Functor
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template <class TInput1, class TInput2, class TOutput>
-class NDBI : public TM4AndTM5IndexBase<TInput1, TInput2, TOutput>
-{
-public:
-  /** Return the index name */
-  std::string GetName() const override
-  {
-    return "NDBI";
-  }
-
-  /// Constructor
-  NDBI() {}
-  /// Desctructor
-  ~NDBI() override {}
-  // Operator on r and nir single pixel values
-protected:
-  inline TOutput Evaluate(const TInput1& pTM4, const TInput2& pTM5) const override
-  {
-    double dTM4 = static_cast<double>(pTM4);
-    double dTM5 = static_cast<double>(pTM5);
-    if (dTM5 + dTM4 == 0)
-      {
-      return static_cast<TOutput>(0.);
-      }
-
-    return (static_cast<TOutput>((dTM5 - dTM4) / (dTM5 + dTM4)));
-  }
-};
-
 /** \class ISU
  *  \brief This functor computes the Index surfaces built (ISU)
  *
@@ -186,56 +37,29 @@ protected:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class ISU : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class ISU : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  ISU() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "ISU";
   }
 
-  /// Constructor
-  ISU() : m_A(100.), m_B(25.) {}
-  /// Desctructor
-  ~ISU() override {}
-
-  /** Set/Get A correction */
-  void SetA(const double pA)
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    m_A = pA;
-  }
-  double GetA(void) const
-  {
-    return (m_A);
-  }
-  /** Set/Get B correction */
-  void SetB(const double pB)
-  {
-    m_B = pB;
-  }
-  double GetB(void) const
-  {
-    return (m_B);
-  }
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
 
-protected:
-  inline TOutput Evaluate(const TInput1& pRed, const TInput2& pNIR) const override
-  {
-    double dRed = static_cast<double>(pRed);
-    double dNIR = static_cast<double>(pNIR);
-    if (dNIR == 0)
-      {
+    if (nir == 0)
+    {
       return static_cast<TOutput>(0.);
       }
 
-    return (static_cast<TOutput>(m_A - (m_B * dRed) / (dNIR)));
+      return (static_cast<TOutput>(A - (B * red) / nir));
   }
 
-private:
-  double m_A;
-  double m_B;
+  static constexpr double A = 100.;
+  static constexpr double B = 25.;
 };
 
 } // namespace Functor
diff --git a/Modules/Radiometry/Indices/include/otbIndicesStackFunctor.h b/Modules/Radiometry/Indices/include/otbIndicesStackFunctor.h
new file mode 100644
index 0000000000000000000000000000000000000000..257f66a84dc7a581a14ad2e1133c2aedf142a7f2
--- /dev/null
+++ b/Modules/Radiometry/Indices/include/otbIndicesStackFunctor.h
@@ -0,0 +1,101 @@
+/*
+ * Copyright (C) 2005-2019 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 otbIndicesStackFunctor_h
+#define otbIndicesStackFunctor_h
+
+#include <vector>
+#include <stdexcept>
+
+namespace otb
+{
+
+namespace Functor
+{
+/**
+ * \class IndicesStackFunctor
+ * \brief A class to compute a stack of radiometric indices
+ *
+ * This functor can be built from a vector of TIndice*. its operator()
+ * will apply each functor of this vector to the input pixel, and
+ * return a VariableLengthVector containing the list resulting
+ * values. It can be used with otb::FunctorImageFilter
+ *
+ * \sa FunctorImageFilter
+ *
+ * \ingroup OTBIndices
+ */
+template <typename TIndice>
+class IndicesStackFunctor
+{
+public:
+  /// Read input / output types from TIndice
+  using IndiceType = TIndice;
+  using PixelType  = typename IndiceType::PixelType;
+  // Output will be a VariableLengthVector of values return by
+  // radiometric indices
+  using OutputType = itk::VariableLengthVector<typename IndiceType::OutputType>;
+
+  /**
+   * \param indices A std::vector<IndiceType*> for indices to compute
+   * the indice stack
+   * \throw std::runtime_error if indices is empty
+   */
+  IndicesStackFunctor(const std::vector<IndiceType*>& indices) : m_Indices(std::move(indices))
+  {
+    if (indices.empty())
+    {
+      throw std::runtime_error("Can not build IndicesStackFunctor from an empty list of indices.");
+    }
+  }
+
+  /**
+   * \param input A itk::VariableLengthVector<TInput> holding the
+   * pixel values for each band
+   * \return A VariableLengthVector<TInput::OutputType> holding all
+   * the indices values
+   */
+  void operator()(OutputType& out, const PixelType& in) const
+  {
+    size_t idx = 0;
+    for (auto indice : m_Indices)
+    {
+      out[idx] = (*indice)(in);
+      ++idx;
+    }
+  }
+  /**
+   * \return the size of the indices list (to be used by FunctorImgeFilter)
+   */
+  size_t OutputSize(...) const
+  {
+    return m_Indices.size();
+  }
+
+private:
+  /// The list of indices to use
+  std::vector<IndiceType*> m_Indices;
+};
+
+} // End namespace Functor
+
+} // End namespace otb
+
+#endif
diff --git a/Modules/Radiometry/Indices/include/otbLandsatTMIndices.h b/Modules/Radiometry/Indices/include/otbLandsatTMIndices.h
index 8dfa9cb4649682049b04f36f20fe789c739724de..aeb68d40248b1d9ae364638a060f0e60fd327b48 100644
--- a/Modules/Radiometry/Indices/include/otbLandsatTMIndices.h
+++ b/Modules/Radiometry/Indices/include/otbLandsatTMIndices.h
@@ -22,8 +22,8 @@
 #define otbLandsatTMIndices_h
 
 #include "otbMath.h"
-#include "itkVariableLengthVector.h"
 #include "otbBandName.h"
+#include "itkVariableLengthVector.h"
 #include "otbFuzzyVariable.h"
 #include <vector>
 #include <algorithm>
diff --git a/Modules/Radiometry/Indices/include/otbNDVIDataNodeFeatureFunction.h b/Modules/Radiometry/Indices/include/otbNDVIDataNodeFeatureFunction.h
index 58120ad259a7666af506d883193f75299d85735a..9d17e0198a0af4e365c7af5747f8edd91e1bb1ad 100644
--- a/Modules/Radiometry/Indices/include/otbNDVIDataNodeFeatureFunction.h
+++ b/Modules/Radiometry/Indices/include/otbNDVIDataNodeFeatureFunction.h
@@ -83,8 +83,7 @@ public:
   typedef PolyLineImageConstIterator<InputImageType, LineType>
                                                       ImageLineIteratorType;
 
-  typedef Functor::NDVI<ScalarRealType, ScalarRealType, PrecisionType>
-                                                      NDVIFunctorType;
+  typedef Functor::NDVI<ScalarRealType, ScalarRealType> NDVIFunctorType;
 
   typedef std::vector<PrecisionType>                  OutputType;
 
@@ -97,22 +96,22 @@ public:
   //TODO replace by metadata parsing
   unsigned int GetREDChannelIndex() const
   {
-    return m_NDVIFunctor.GetRedIndex()+1;
+    return m_NDVIFunctor.GetBandIndex(CommonBandNames::RED);
   }
 
   void SetREDChannelIndex(unsigned int id)
   {
-    m_NDVIFunctor.SetRedIndex(id-1);
+    m_NDVIFunctor.SetBandIndex(CommonBandNames::RED, id);
   }
 
   unsigned int GetNIRChannelIndex() const
   {
-    return m_NDVIFunctor.GetNIRIndex()+1;
+    return m_NDVIFunctor.GetBandIndex(CommonBandNames::NIR);
   }
 
   void SetNIRChannelIndex(unsigned int id)
   {
-    m_NDVIFunctor.SetNIRIndex(id-1);
+    m_NDVIFunctor.SetBandIndex(CommonBandNames::NIR, id);
   }
 
 protected:
diff --git a/Modules/Radiometry/Indices/include/otbNDVIDataNodeFeatureFunction.hxx b/Modules/Radiometry/Indices/include/otbNDVIDataNodeFeatureFunction.hxx
index 910bc13c1f3baf276218f7af94ceb33a9fc0edb7..f6f33564d5c29ed8cec88bf73b4a6445df28c076 100644
--- a/Modules/Radiometry/Indices/include/otbNDVIDataNodeFeatureFunction.hxx
+++ b/Modules/Radiometry/Indices/include/otbNDVIDataNodeFeatureFunction.hxx
@@ -106,8 +106,7 @@ typename NDVIDataNodeFeatureFunction<TImage, TCoordRep, TPrecision>
     if(this->IsInsideBuffer(lineIt.GetIndex()))
       {
       PixelType pixel = this->GetInputImage()->GetPixel(lineIt.GetIndex());
-      if(m_NDVIFunctor(pixel [this->GetREDChannelIndex() - 1],
-                       pixel [this->GetNIRChannelIndex() - 1]) >= this->GetNDVIThreshold())
+      if(m_NDVIFunctor(pixel) >= this->GetNDVIThreshold())
         {
         nbValidPixel += 1;
         }
diff --git a/Modules/Radiometry/Indices/include/otbRadiometricIndex.h b/Modules/Radiometry/Indices/include/otbRadiometricIndex.h
new file mode 100644
index 0000000000000000000000000000000000000000..d65213a2ec63b3a6d40e575b01b5f29b5847e53c
--- /dev/null
+++ b/Modules/Radiometry/Indices/include/otbRadiometricIndex.h
@@ -0,0 +1,217 @@
+/*
+ * Copyright (C) 2005-2019 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 otbRadiometricIndex_h
+#define otbRadiometricIndex_h
+
+#include "itkVariableLengthVector.h"
+#include "otbBandName.h"
+#include <array>
+#include <set>
+#include <string>
+#include <map>
+#include <stdexcept>
+
+using namespace otb::BandName;
+
+namespace otb
+{
+namespace Functor
+{
+/**
+ * \class RadiometricIndex
+ * \brief Base class for all radiometric indices
+ *
+ * This class is the base class for all radiometric indices.
+ *
+ * It offers services to:
+ * - Indicate which band are required among CommonBandNames enum
+ * - Set indices of each required band
+ * - Compute the indice response to a pixel by subclassing the pure
+ * virtual operator()
+ *
+ * This class is designed for performance on the critical path. For
+ * best performances use the Value() method when implementing
+ * operator() to avoid branches.
+ *
+ * \ingroup OTBIndices
+ */
+template <typename TInput, typename TOutput>
+class RadiometricIndex
+{
+public:
+  /// Types for input/output
+  using InputType  = TInput;
+  using PixelType  = itk::VariableLengthVector<InputType>;
+  using OutputType = TOutput;
+
+  /// Enum Among which bands are used
+  using BandNameType = CommonBandNames;
+
+  /// The number of bands, derived from the Enum MAX value
+  static constexpr size_t NumberOfBands = static_cast<size_t>(BandNameType::MAX);
+
+  static constexpr double Epsilon = 0.0000001;
+
+  /**
+   * \param requiredBands the set<CommonBandNames> of required bands
+   * \throw runtime_error if requiredBands contains CommonBandNames::MAX
+   */
+  RadiometricIndex(const std::set<BandNameType>& requiredBands) : m_RequiredBands(), m_BandIndices()
+  {
+    if (requiredBands.find(BandNameType::MAX) != requiredBands.end())
+    {
+      throw std::runtime_error("TBandNameEnum::MAX can not be used as a required band");
+    }
+
+    // Fill the required bands array
+    m_RequiredBands.fill(false);
+    m_BandIndices.fill(0);
+
+    for (auto b : requiredBands)
+    {
+      m_RequiredBands[static_cast<size_t>(b)] = true;
+    }
+  }
+
+  /**
+   * \return a set<CommandBandName> containing the required bands for
+   * this indice.
+   */
+  std::set<BandNameType> GetRequiredBands() const
+  {
+    std::set<BandNameType> resp;
+    for (size_t i = 0; i < NumberOfBands; ++i)
+    {
+      if (m_RequiredBands[i])
+      {
+        resp.insert(static_cast<BandNameType>(i));
+      }
+    }
+
+    return resp;
+  }
+
+  /**
+   * \param band The band to set (value in CommandBandName)
+   * \param index The index of the band to set (starts at 1 for first band)
+   * \throw runtime_error if band is CommandBandName::MAX
+   */
+  void SetBandIndex(BandNameType band, size_t index)
+  {
+    if (band == BandNameType::MAX)
+    {
+      throw std::runtime_error("Can not set index for CommandBandName::MAX");
+    }
+    m_BandIndices[static_cast<size_t>(band)] = index;
+  }
+
+  /**
+   * \param indicesMap a std::map<CommandBandName,size_t> containing all
+   * bands indices to set  (starts at 1 for first band)
+   * \throw runtime_error if indicesMap contains CommandBandName::MAX
+   */
+  void SetBandsIndices(const std::map<BandNameType, size_t>& indicesMap)
+  {
+    for (auto it : indicesMap)
+    {
+      SetBandIndex(it.first, it.second);
+    }
+  }
+
+  /**
+   * \param band The band for which to retrieve indice
+   * \return The indices of the band
+   * \throw runtime_error if band is CommandBandName::MAX
+   */
+  size_t GetBandIndex(BandNameType band) const
+  {
+    if (band == BandNameType::MAX)
+    {
+      throw std::runtime_error("Can not get index for CommandBandName::MAX");
+    }
+    return m_BandIndices[static_cast<size_t>(band)];
+  }
+
+  /**
+   * Astract method which will compute the radiometric indice
+   * \param input A itk::VariableLengthVector<TInput> holding the
+   * pixel values for each band
+   * \return The indice value as TOutput  (starts at 1 for first band)
+   */
+  virtual TOutput operator()(const itk::VariableLengthVector<TInput>& input) const = 0;
+
+protected:
+  /**
+   * Helper method to retrieve index for band name. With respect to
+   * the public method, this method will not throw an exception if
+   * CommandBandName::MAX is used as a parameter. Since it is meant for
+   * internal use in the critical path and not for client code, it
+   * will only assert that band is not CommandBandName::MAX in debug
+   * mode.
+   *
+   * \param band The band for which to retrieve indice
+   * \return The indices of the band
+   */
+  size_t UncheckedBandIndex(BandNameType band) const
+  {
+    assert(band != BandNameType::MAX && "Can not retrieve index for band CommandBandName::MAX");
+    return m_BandIndices[static_cast<size_t>(band)];
+  }
+
+  /**
+   * Helper method to parse input  itk::VariableLengthVector<TInput>
+   * and get the corresponding band value.
+   * For instance:
+   * \snippet auto red   = this->Value(CommonBandNamess::RED,input);
+   *
+   * As this function is on the critical performance path, no checks
+   * are made to see wether this band is really required for this
+   * indice. However an assertion will be raised in debug mode.
+   *
+   * \param band The band for which to retrieve the value
+   * \param input A itk::VariableLengthVector<TInput> holding the
+   * pixel values for each band
+   * \return The value of the band as double
+   *
+   */
+  double Value(BandNameType band, const itk::VariableLengthVector<TInput>& input) const
+  {
+    assert(m_RequiredBands[band] && "Retrieving value for a band that is not in the required bands list");
+    return static_cast<double>(input[UncheckedBandIndex(band) - 1]);
+  }
+
+private:
+  // Explicitely disable default constructor
+  RadiometricIndex() = delete;
+
+  /// An array storing the required status for each band
+  using RequiredBandsContainer = std::array<bool, NumberOfBands>;
+  RequiredBandsContainer m_RequiredBands;
+
+  /// An array storing the indice for each band
+  using BandIndicesContainer = std::array<size_t, NumberOfBands>;
+  BandIndicesContainer m_BandIndices;
+};
+
+} // namespace Functor
+} // End namespace otb
+
+#endif
diff --git a/Modules/Radiometry/Indices/include/otbSoilIndicesFunctor.h b/Modules/Radiometry/Indices/include/otbSoilIndicesFunctor.h
index 1f8ee072431c31fbcd476db730c772ca2e80fec1..aaedbc7c0efe62610a54a97817c87e52ff26dd1d 100644
--- a/Modules/Radiometry/Indices/include/otbSoilIndicesFunctor.h
+++ b/Modules/Radiometry/Indices/include/otbSoilIndicesFunctor.h
@@ -22,255 +22,14 @@
 #define otbSoilIndicesFunctor_h
 
 #include "otbMath.h"
-#include "itkVariableLengthVector.h"
-#include "otbBandName.h"
-#include <string>
+#include "otbRadiometricIndex.h"
 
 namespace otb
 {
 namespace Functor
 {
-/**
- * \class GAndRIndexBase
- *
- * \brief Base class for Green And Red channels of Spot Images
- *  XS1 corresponds to the green channel
- *  XS2 corresponds to the red channel
- *  XS3 corresponds to the Nir channel
- *  XS4 corresponds to the Mir channel (for Spot 4 & 5)
- *  Implement operators for UnaryFunctorImageFilter templated with a
- *  VectorImage and BinaryFunctorImageFilter templated with single
- *  images.
- *  Subclasses should NOT overload operators, they must  re-implement
- *  the Evaluate() method.
- *
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
-*/
-template<class TInput1, class TInput2, class TOutput>
-class GAndRIndexBase
-{
-public:
-  /// Vector pixel type used to support both vector images and multiple
-  /// input images
-  typedef itk::VariableLengthVector<TInput1> InputVectorType;
-
-  //operators !=
-  bool operator !=(const GAndRIndexBase&) const
-  {
-    return true;
-  }
-  //operator ==
-  bool operator ==(const GAndRIndexBase& other) const
-  {
-    return !(*this != other);
-  }
-
-  // Operator on vector pixel type
-  inline TOutput operator ()(const InputVectorType& inputVector) const
-  {
-    return this->Evaluate(inputVector[m_GreenIndex - 1], static_cast<TInput2>(inputVector[m_RedIndex - 1]));
-  }
-
-  // Binary operator
-  inline TOutput operator ()(const TInput1& g, const TInput2& r) const
-  {
-    return this->Evaluate(g, r);
-  }
-  /// Constructor
-  GAndRIndexBase() : m_EpsilonToBeConsideredAsZero(0.0000001), m_GreenIndex(1), m_RedIndex(2) {}
-  /// Desctructor
-  virtual ~GAndRIndexBase() {}
-
-  /// Set Green Index
-  void SetGreenIndex(unsigned int channel)
-  {
-    m_GreenIndex = channel;
-  }
-  /// Get Green Index
-  unsigned int GetGreenIndex() const
-  {
-    return m_GreenIndex;
-  }
-  /// Set Red Index
-  void SetRedIndex(unsigned int channel)
-  {
-    m_RedIndex = channel;
-  }
-  /// Get Red Index
-  unsigned int GetRedIndex() const
-  {
-    return m_RedIndex;
-  }
-
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::RED)
-      {
-      m_RedIndex = channel;
-      }
-    if (band == BandName::GREEN)
-      {
-      m_GreenIndex = channel;
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::RED)
-      {
-      return m_RedIndex;
-      }
-    if (band == BandName::GREEN)
-      {
-      return m_GreenIndex;
-      }
-  }
-
-  /** Return the index name */
-  virtual std::string GetName() const = 0;
-
-protected:
-  // This method must be reimplemented in subclasses to actually
-  // compute the index value
-  virtual TOutput Evaluate(const TInput1& g, const TInput2& r) const = 0;
-  double m_EpsilonToBeConsideredAsZero;
-
-private:
-  unsigned int m_GreenIndex;
-  unsigned int m_RedIndex;
-};
-
-/**
- * \class GAndRAndNirIndexBase
- * \brief Base class for Green And Red And NIR channels of Spot Images
- *
- *
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template<class TInput1, class TInput2, class TInput3, class TOutput>
-class GAndRAndNirIndexBase
-{
-public:
-  /// Vector pixel type used to support both vector images and multiple
-  /// input images
-  typedef itk::VariableLengthVector<TInput1> InputVectorType;
-
-  //operators !=
-  bool operator !=(const GAndRAndNirIndexBase&) const
-  {
-    return true;
-  }
-  //operator ==
-  bool operator ==(const GAndRAndNirIndexBase& other) const
-  {
-    return !(*this != other);
-  }
-
-  // Operator on vector pixel type
-  inline TOutput operator ()(const InputVectorType& inputVector) const
-  {
-    return this->Evaluate(static_cast<TInput1>(inputVector[m_GreenIndex - 1]),
-                          static_cast<TInput2>(inputVector[m_RedIndex - 1]),
-                          static_cast<TInput3>(inputVector[m_NIRIndex - 1]));
-  }
-
-  // Binary operator
-  inline TOutput operator ()(const TInput1& g, const TInput2& r, const TInput2& nir) const
-  {
-    return this->Evaluate(g, r, nir);
-  }
-  /// Constructor
-  GAndRAndNirIndexBase() : m_EpsilonToBeConsideredAsZero(0.0000001), m_GreenIndex(1), m_RedIndex(2),  m_NIRIndex(3) {}
-  /// Desctructor
-  virtual ~GAndRAndNirIndexBase() {}
-
-  /// Set Green Index
-  void SetGreenIndex(unsigned int channel)
-  {
-    m_GreenIndex = channel;
-  }
-  /// Get Green Index
-  unsigned int GetGreenIndex() const
-  {
-    return m_GreenIndex;
-  }
-  /// Set Red Index
-  void SetRedIndex(unsigned int channel)
-  {
-    m_RedIndex = channel;
-  }
-  /// Get Red Index
-  unsigned int GetRedIndex() const
-  {
-    return m_RedIndex;
-  }
-  /// Set Nir Index
-  void SetNIRIndex(unsigned int channel)
-  {
-    m_NIRIndex = channel;
-  }
-  /// Get Nir Index
-  unsigned int GetNIRIndex() const
-  {
-    return m_NIRIndex;
-  }
-
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::RED)
-      {
-      m_RedIndex = channel;
-      }
-    if (band == BandName::GREEN)
-      {
-      m_GreenIndex = channel;
-      }
-    if (band == BandName::NIR)
-      {
-      m_NIRIndex = channel;
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::RED)
-      {
-      return m_RedIndex;
-      }
-    if (band == BandName::GREEN)
-      {
-      return m_GreenIndex;
-      }
-    if (band == BandName::NIR)
-      {
-      return m_NIRIndex;
-      }
-  }
-
-  /** Return the index name */
-  virtual std::string GetName() const = 0;
-
-protected:
-  // This method must be reimplemented in subclasses to actually
-  // compute the index value
-  virtual TOutput Evaluate(const TInput1& g, const TInput2& r, const TInput2& nir) const = 0;
-
-  double m_EpsilonToBeConsideredAsZero;
-
-private:
-  unsigned int m_GreenIndex;
-  unsigned int m_RedIndex;
-  unsigned int m_NIRIndex;
-};
-
-/** \class IR
- *  \brief This functor computes the Redness Index (IR)
+/** \class RI
+ *  \brief This functor computes the Redness Index (RI)
  *
  *  [Pouget et al., "Caracteristiques spectrales des surfaces sableuses
  *   de la region cotiere nord-ouest de l'Egypte: application aux donnees
@@ -283,36 +42,29 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class IR : public GAndRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class RI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  RI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::GREEN})
   {
-    return "IR";
   }
 
-  /// Constructor
-  IR() {}
-  /// Desctructor
-  ~IR() override {}
-  // Operator on r and nir single pixel values
-protected:
-  inline TOutput Evaluate(const TInput1& pGreen, const TInput2& pRed) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    double dGreen = static_cast<double>(pGreen);
-    double dRed = static_cast<double>(pRed);
-    if (std::abs(dGreen) < this->m_EpsilonToBeConsideredAsZero)
-      {
+    auto green = this->Value(CommonBandNames::GREEN, input);
+    auto red   = this->Value(CommonBandNames::RED, input);
+
+    if (std::abs(green) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       return static_cast<TOutput>(0.);
       }
 
-    return static_cast<TOutput>(dRed * dRed / (dGreen * dGreen * dGreen));
+      return static_cast<TOutput>(red * red / (green * green * green));
   }
 };
 
-/** \class IC
+/** \class CI
  *  \brief This functor computes the Color Index (IC)
  *
  *  [Pouget et al., "Caracteristiques spectrales des surfaces sableuses
@@ -326,37 +78,30 @@ protected:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class IC : public GAndRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class CI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  CI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::GREEN})
   {
-    return "IC";
   }
 
-  /// Constructor
-  IC() {}
-  /// Desctructor
-  ~IC() override {}
-  // Operator on r and nir single pixel values
-protected:
-  inline TOutput Evaluate(const TInput1& pGreen, const TInput2& pRed) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    double dGreen = static_cast<double>(pGreen);
-    double dRed = static_cast<double>(pRed);
-    if (std::abs(dGreen + dRed) < this->m_EpsilonToBeConsideredAsZero)
-      {
+    auto green = this->Value(CommonBandNames::GREEN, input);
+    auto red   = this->Value(CommonBandNames::RED, input);
+
+    if (std::abs(green + red) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       return static_cast<TOutput>(0.);
       }
 
-    return (static_cast<TOutput>((dRed - dGreen) / (dRed + dGreen)));
+      return (static_cast<TOutput>((red - green) / (red + green)));
   }
 };
 
-/** \class IB
- *  \brief This functor computes the Brilliance Index (IB)
+/** \class BI
+ *  \brief This functor computes the Brilliance Index (BI)
  *
  *  [ ]
  *
@@ -365,33 +110,25 @@ protected:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class IB : public GAndRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class BI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  BI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::GREEN})
   {
-    return "IB";
   }
 
-  /// Constructor
-  IB() {}
-  /// Desctructor
-  ~IB() override {}
-  // Operator on r and nir single pixel values
-protected:
-  inline TOutput Evaluate(const TInput1& pGreen, const TInput2& pRed) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    double dGreen = static_cast<double>(pGreen);
-    double dRed = static_cast<double>(pRed);
+    auto green = this->Value(CommonBandNames::GREEN, input);
+    auto red   = this->Value(CommonBandNames::RED, input);
 
-    return (static_cast<TOutput>(std::sqrt((dRed * dRed + dGreen * dGreen) / 2.)));
+    return (static_cast<TOutput>(std::sqrt((red * red + green * green) / 2.)));
   }
 };
 
-/** \class IB2
- *  \brief This functor computes the Brilliance Index (IB2)
+/** \class BI2
+ *  \brief This functor computes the Brilliance Index (BI2)
  *
  *  [ ]
  *
@@ -400,29 +137,21 @@ protected:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TInput3, class TOutput>
-class IB2 : public GAndRAndNirIndexBase<TInput1, TInput2, TInput3, TOutput>
+template <class TInput, class TOutput>
+class BI2 : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  BI2() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::GREEN, CommonBandNames::NIR})
   {
-    return "IB2";
   }
 
-  /// Constructor
-  IB2() {}
-  /// Desctructor
-  ~IB2() override {}
-  // Operator on r and nir single pixel values
-protected:
-  inline TOutput Evaluate(const TInput1& pGreen, const TInput2& pRed, const TInput2& pNir) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    double dGreen = static_cast<double>(pGreen);
-    double dRed = static_cast<double>(pRed);
-    double dNir = static_cast<double>(pNir);
+    auto green = this->Value(CommonBandNames::GREEN, input);
+    auto red   = this->Value(CommonBandNames::RED, input);
+    auto nir   = this->Value(CommonBandNames::NIR, input);
 
-    return (static_cast<TOutput>(std::sqrt((dRed * dRed + dGreen * dGreen + dNir * dNir) / 3.)));
+    return (static_cast<TOutput>(std::sqrt((red * red + green * green + nir * nir) / 3.)));
   }
 };
 
diff --git a/Modules/Radiometry/Indices/include/otbVegetationIndicesFunctor.h b/Modules/Radiometry/Indices/include/otbVegetationIndicesFunctor.h
index 7131e2cdc08cf575a0070e846b02ae172b0c516d..63bc7299527ad45a5ab829f6a69d53f6cab3ec94 100644
--- a/Modules/Radiometry/Indices/include/otbVegetationIndicesFunctor.h
+++ b/Modules/Radiometry/Indices/include/otbVegetationIndicesFunctor.h
@@ -22,383 +22,13 @@
 #define otbVegetationIndicesFunctor_h
 
 #include "otbMath.h"
-#include "itkVariableLengthVector.h"
-#include "otbBandName.h"
-#include <string>
+#include "otbRadiometricIndex.h"
 
 namespace otb
 {
 
 namespace Functor
 {
-
-/**
-   * \class RAndNIRIndexBase
-   * \brief Base class for R And NIR based Index
-   *
-   *  Implement operators for UnaryFunctorImageFilter templated with a
-   *  VectorImage and BinaryFunctorImageFilter templated with single
-   *  images.
-   *  Subclasses should NOT overload operators, they must  re-implement
-   *  the Evaluate() method.
-   *
-   * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template<class TInput1, class TInput2, class TOutput>
-class RAndNIRIndexBase
-{
-public:
-  /// Vector pixel type used to support both vector images and multiple
-  /// input images
-  typedef itk::VariableLengthVector<TInput1> InputVectorType;
-
-  //operators !=
-  bool operator !=(const RAndNIRIndexBase&) const
-  {
-    return true;
-  }
-  //operator ==
-  bool operator ==(const RAndNIRIndexBase& other) const
-  {
-    return !(*this != other);
-  }
-
-  // Operator on vector pixel type
-  inline TOutput operator ()(const InputVectorType& inputVector) const
-  {
-    return this->Evaluate(inputVector[m_RedIndex - 1], static_cast<TInput2>(inputVector[m_NIRIndex - 1]));
-  }
-
-  // Binary operator
-  inline TOutput operator ()(const TInput1& r, const TInput2& nir) const
-  {
-    return this->Evaluate(r, nir);
-  }
-  /// Constructor
-  RAndNIRIndexBase() :  m_EpsilonToBeConsideredAsZero(0.0000001), m_RedIndex(3), m_NIRIndex(4) {}
-  /// Desctructor
-  virtual ~RAndNIRIndexBase() {}
-
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::RED)
-      {
-      m_RedIndex = channel;
-      }
-    if (band == BandName::NIR)
-      {
-      m_NIRIndex = channel;
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::RED)
-      {
-      return m_RedIndex;
-      }
-    if (band == BandName::NIR)
-      {
-      return m_NIRIndex;
-      }
-  }
-
-  /// Set Red Index
-  void SetRedIndex(unsigned int channel)
-  {
-    m_RedIndex = channel;
-  }
-  /// Get Red Index
-  unsigned int GetRedIndex() const
-  {
-    return m_RedIndex;
-  }
-  /// Set NIR Index
-  void SetNIRIndex(unsigned int channel)
-  {
-    m_NIRIndex = channel;
-  }
-  /// Get NIR Index
-  unsigned int GetNIRIndex() const
-  {
-    return m_NIRIndex;
-  }
-
-  /** Return the index name */
-  virtual std::string GetName() const = 0;
-
-protected:
-  // This method must be reimplemented in subclasses to actually
-  // compute the index value
-  virtual TOutput Evaluate(const TInput1& r, const TInput2& nir) const = 0;
-  double m_EpsilonToBeConsideredAsZero;
-
-private:
-  unsigned int m_RedIndex;
-  unsigned int m_NIRIndex;
-};
-
-/**
- * \class RAndBAndNIRIndexBase
- * \brief base class for R, B And NIR based Index
- *  Implement operators for UnaryFunctorImageFilter templated with a
- *  VectorImage and BinaryFunctorImageFilter templated with single
- *  images.
- *  Subclasses should NOT overload operators, they must  re-implement
- *  the Evaluate() method.
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template<class TInput1, class TInput2, class TInput3, class TOutput>
-class RAndBAndNIRIndexBase
-{
-public:
-  /// Vector pixel type used to support both vector images and multiple
-  /// input images
-  typedef itk::VariableLengthVector<TInput1> InputVectorType;
-
-  //operators !=
-  bool operator !=(const RAndBAndNIRIndexBase&) const
-  {
-    return true;
-  }
-
-  //operator ==
-  bool operator ==(const RAndBAndNIRIndexBase& other) const
-  {
-    return !(*this != other);
-  }
-
-  // Operator on vector pixel type
-  inline TOutput operator ()(const InputVectorType& inputVector)
-  {
-    return this->Evaluate(inputVector[m_RedIndex - 1],
-                          static_cast<TInput2>(inputVector[m_BlueIndex - 1]),
-                          static_cast<TInput3>(inputVector[m_NIRIndex - 1]));
-  }
-  // Binary operator
-  inline TOutput operator ()(const TInput1& r, const TInput2& b, const TInput2& nir)
-  {
-    return this->Evaluate(r, b, nir);
-  }
-  /// Constructor
-  RAndBAndNIRIndexBase() : m_EpsilonToBeConsideredAsZero(0.0000001), m_RedIndex(3), m_BlueIndex(1), m_NIRIndex(4) {}
-  /// Desctructor
-  virtual ~RAndBAndNIRIndexBase() {}
-
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::RED)
-      {
-      m_RedIndex = channel;
-      }
-    if (band == BandName::BLUE)
-      {
-      m_BlueIndex = channel;
-      }
-    if (band == BandName::NIR)
-      {
-      m_NIRIndex = channel;
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::RED)
-      {
-      return m_RedIndex;
-      }
-    if (band == BandName::BLUE)
-      {
-      return m_BlueIndex;
-      }
-    if (band == BandName::NIR)
-      {
-      return m_NIRIndex;
-      }
-  }
-
-  /// Set Red Index
-  void SetRedIndex(unsigned int channel)
-  {
-    m_RedIndex = channel;
-  }
-  /// Get Red Index
-  unsigned int GetRedIndex() const
-  {
-    return m_RedIndex;
-  }
-  /// Set Blue Index
-  void SetBlueIndex(unsigned int channel)
-  {
-    m_BlueIndex = channel;
-  }
-  /// Get Blue Index
-  unsigned int GetBlueIndex() const
-  {
-    return m_BlueIndex;
-  }
-
-  /// Set NIR Index
-  void SetNIRIndex(unsigned int channel)
-  {
-    m_NIRIndex = channel;
-  }
-  /// Get NIR Index
-  unsigned int GetNIRIndex() const
-  {
-    return m_NIRIndex;
-  }
-
-  /** Return the index name */
-  virtual std::string GetName() const = 0;
-
-protected:
-  // This method must be reimplemented in subclasses to actually
-  // compute the index value
-  virtual TOutput Evaluate(const TInput1& r, const TInput2& b, const TInput3& nir) const = 0;
-  double m_EpsilonToBeConsideredAsZero;
-
-private:
-  unsigned int m_RedIndex;
-  unsigned int m_BlueIndex;
-  unsigned int m_NIRIndex;
-};
-
-/**
- * \class RAndGAndNIRIndexBase
- * \brief base class for R, G And NIR based Index
- *  Implement operators for UnaryFunctorImageFilter templated with a
- *  VectorImage and BinaryFunctorImageFilter templated with single
- *  images.
- *  Subclasses should NOT overload operators, they must  re-implement
- *  the Evaluate() method.
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template<class TInput1, class TInput2, class TInput3, class TOutput>
-class RAndGAndNIRIndexBase
-{
-public:
-  /// Vector pixel type used to support both vector images and multiple
-  /// input images
-  typedef itk::VariableLengthVector<TInput1> InputVectorType;
-
-  //operators !=
-  bool operator !=(const RAndGAndNIRIndexBase&) const
-  {
-    return true;
-  }
-  //operator ==
-  bool operator ==(const RAndGAndNIRIndexBase& other) const
-  {
-    return !(*this != other);
-  }
-
-  // Operator on vector pixel type
-  inline TOutput operator ()(const InputVectorType& inputVector)
-  {
-    return this->Evaluate(inputVector[m_RedIndex - 1],
-                          static_cast<TInput2>(inputVector[m_GreenIndex - 1]),
-                          static_cast<TInput3>(inputVector[m_NIRIndex - 1]));
-  }
-
-  // Binary operator
-  inline TOutput operator ()(const TInput1& r, const TInput2& g, const TInput2& nir)
-  {
-    return this->Evaluate(r, g, nir);
-  }
-  /// Constructor
-  RAndGAndNIRIndexBase() : m_EpsilonToBeConsideredAsZero(0.0000001), m_RedIndex(3), m_GreenIndex(2), m_NIRIndex(4)  {}
-  /// Desctructor
-  virtual ~RAndGAndNIRIndexBase() {}
-
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::RED)
-      {
-      m_RedIndex = channel;
-      }
-    if (band == BandName::GREEN)
-      {
-      m_GreenIndex = channel;
-      }
-    if (band == BandName::NIR)
-      {
-      m_NIRIndex = channel;
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::RED)
-      {
-      return m_RedIndex;
-      }
-    if (band == BandName::GREEN)
-      {
-      return m_GreenIndex;
-      }
-    if (band == BandName::NIR)
-      {
-      return m_NIRIndex;
-      }
-  }
-
-  /// Set Red Index
-  void SetRedIndex(unsigned int channel)
-  {
-    m_RedIndex = channel;
-  }
-  /// Get Red Index
-  unsigned int GetRedIndex() const
-  {
-    return m_RedIndex;
-  }
-  /// Set Green Index
-  void SetGreenIndex(unsigned int channel)
-  {
-    m_GreenIndex = channel;
-  }
-  /// Get Green Index
-  unsigned int GetGreenIndex() const
-  {
-    return m_GreenIndex;
-  }
-
-  /// Set NIR Index
-  void SetNIRIndex(unsigned int channel)
-  {
-    m_NIRIndex = channel;
-  }
-  /// Get NIR Index
-  unsigned int GetNIRIndex() const
-  {
-    return m_NIRIndex;
-  }
-
-  /** Return the index name */
-  virtual std::string GetName() const = 0;
-
-protected:
-  // This method must be reimplemented in subclasses to actually
-  // compute the index value
-  virtual TOutput Evaluate(const TInput1& r, const TInput2& g, const TInput3& nir) const = 0;
-  double m_EpsilonToBeConsideredAsZero;
-
-private:
-  unsigned int m_RedIndex;
-  unsigned int m_GreenIndex;
-  unsigned int m_NIRIndex;
-};
-
 /** \class NDVI
  *  \brief This functor computes the Normalized Difference Vegetation Index (NDVI)
  *
@@ -409,33 +39,31 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class NDVI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class NDVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
+  NDVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
+  {
+  }
 
-  /** Return the index name */
-  std::string GetName() const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    return "NDVI";
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
+
+    return static_cast<TOutput>(Compute(red, nir));
   }
 
-  /// Constructor
-  NDVI() {}
-  /// Desctructor
-  ~NDVI() override {}
-  // Operator on r and nir single pixel values
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
+  // This static compute will be used in indices derived from NDVI
+  static double Compute(const double& red, const double& nir)
   {
-    double dr = static_cast<double>(r);
-    double dnir = static_cast<double>(nir);
-    if (std::abs(dnir + dr) < this->m_EpsilonToBeConsideredAsZero)
-      {
-      return static_cast<TOutput>(0.);
+    if (std::abs(nir + red) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
+      return 0.;
       }
 
-    return (static_cast<TOutput>((dnir - dr) / (dnir + dr)));
+      return (nir - red) / (nir + red);
   }
 };
 
@@ -449,29 +77,24 @@ protected:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class RVI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class RVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-
-  /** Return the index name */
-  std::string GetName() const override
+  RVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "RVI";
   }
 
-  RVI() {}
-  ~RVI() override {}
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    double dr = static_cast<double>(r);
-    double dnir = static_cast<double>(nir);
-    if (std::abs(dr)  < this->m_EpsilonToBeConsideredAsZero)
-      {
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
+
+    if (std::abs(red) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       return static_cast<TOutput>(0.);
       }
-    return (static_cast<TOutput>(dnir / dr));
+      return (static_cast<TOutput>(nir / red));
   }
 };
 
@@ -484,57 +107,31 @@ protected:
  * C. L. Wiegand, A. J. Richardson, D. E. Escobar, and A. H. Gerbermann,
  * "Vegetation Indices in Crop Assessments", REMOTE SENS. ENVIRON. 35:105-119 (1991)
  *
- *  \ingroup Functor2
+ *  \ingroup Functor
  * \ingroup Radiometry
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class PVI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class PVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  PVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "PVI";
   }
 
-  PVI() : m_A(0.90893), m_B(7.46216), m_Coeff(0.74) {}
-  ~PVI() override {}
-  /** Set/Get A and B parameters */
-  void SetA(const double A)
-  {
-    m_A = A;
-    m_Coeff = 1. / (std::sqrt(m_A * m_A + 1.));
-  }
-  double GetA(void) const
-  {
-    return (m_A);
-  }
-  void SetB(const double B)
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    m_B = B;
-  }
-  double GetB(void) const
-  {
-    return (m_B);
-  }
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
-  {
-    double dnir = static_cast<double>(nir);
-    double dr = static_cast<double>(r);
-    return (static_cast<TOutput>((dnir - m_A * dr - m_B) * m_Coeff));
-  }
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
 
-private:
+    return (static_cast<TOutput>((nir - A * red - B) * C));
+  }
 
   /** A and B parameters */
-  double m_A;
-  double m_B;
-  /** Denominator, pre-calculed when the A variable is set */
-  double m_Coeff;
-
+  static constexpr double A = 0.90893;
+  static constexpr double B = 7.46216;
+  static constexpr double C = 9.74;
 };
 
 /** \class SAVI
@@ -547,48 +144,28 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class SAVI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class SAVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-
-  /** Return the index name */
-  std::string GetName() const override
+  SAVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "SAVI";
   }
 
-  SAVI() : m_L(0.5) {}
-  ~SAVI() override {}
-
-  /** Set/Get L correction */
-  void SetL(const double L)
-  {
-    m_L = L;
-  }
-  double GetL(void) const
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    return (m_L);
-  }
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
 
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
-  {
-    double dnir = static_cast<double>(nir);
-    double dr = static_cast<double>(r);
-    double denominator = dnir + dr + m_L;
-    if (std::abs(denominator)  < this->m_EpsilonToBeConsideredAsZero)
-      {
+    if (std::abs(nir + red + L) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       return static_cast<TOutput>(0.);
       }
-    return (static_cast<TOutput>(((dnir - dr) * (1 + m_L)) / denominator));
+      return (static_cast<TOutput>(((nir - red) * (1 + L)) / (nir + red + L)));
   }
 
-private:
-
   /** L correction */
-  double m_L;
-
+  static constexpr double L = 0.5;
 };
 
 /** \class TSAVI
@@ -601,68 +178,33 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class TSAVI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class TSAVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-
-  /** Return the index name */
-  std::string GetName() const override
+  TSAVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "TSAVI";
   }
 
-  TSAVI() : m_A(0.7), m_S(0.9), m_X(0.08) {}
-  ~TSAVI() override {}
-
-  /** Set/Get S and A parameters */
-  void SetS(const double S)
-  {
-    m_S = S;
-  }
-  double GetS(void) const
-  {
-    return (m_S);
-  }
-  void SetA(const double A)
-  {
-    m_A = A;
-  }
-  double GetA(void) const
-  {
-    return (m_A);
-  }
-  /** Set/Get X parameter */
-  void SetX(const double X)
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    m_X = X;
-  }
-  double GetX(void) const
-  {
-    return (m_X);
-  }
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
 
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
-  {
-    double dnir = static_cast<double>(nir);
-    double dr = static_cast<double>(r);
-    double denominator = m_A * dnir + dr + m_X * (1. + m_A * m_A);
-    if (std::abs(denominator) < this->m_EpsilonToBeConsideredAsZero)
-      {
+    double denominator = A * nir + red + X * (1. + A * A);
+
+    if (std::abs(denominator) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       return static_cast<TOutput>(0.);
       }
-    return (static_cast<TOutput>((m_A * (dnir - m_A * dr - m_S)) / denominator));
+      return (static_cast<TOutput>((A * (nir - A * red - S)) / denominator));
   }
 
-private:
-
   /** A and S parameters */
-  double m_A;
-  double m_S;
+  static constexpr double A = 0.7;
+  static constexpr double S = 0.9;
   /** X parameter */
-  double m_X;
-
+  static constexpr double X = 0.08;
 };
 
 /** \class WDVI
@@ -675,41 +217,30 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class WDVI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class WDVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  /// Constructor
+  WDVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "WDVI";
   }
 
-  /// Constructor
-  WDVI() : m_S(0.4) {}
-  /// Desctructor
-  ~WDVI() override {}
-  // Operator on r and nir single pixel values
-/** Set/Get Slop of soil line */
-  void SetS(const double s)
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    m_S = s;
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
+
+    return static_cast<TOutput>(Compute(red, nir));
   }
-  double GetS(void) const
+
+  static double Compute(const double& red, const double& nir)
   {
-    return (m_S);
+    return (nir - S * red);
   }
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
-  {
-    double dr = static_cast<double>(r);
-    double dnir = static_cast<double>(nir);
 
-    return (dnir - m_S * dr);
-  }
-private:
   /** Slope of soil line */
-  double m_S;
+  static constexpr double S = 0.4;
 };
 
 /** \class MSAVI
@@ -723,69 +254,37 @@ private:
  * \ingroup OTBIndices
  */
 
-template <class TInput1, class TInput2, class TOutput>
-class MSAVI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class MSAVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  MSAVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "MSAVI";
   }
 
-  typedef NDVI<TInput1, TInput2, TOutput> NDVIFunctorType;
-  typedef SAVI<TInput1, TInput2, TOutput> SAVIFunctorType;
-  typedef WDVI<TInput1, TInput2, TOutput> WDVIFunctorType;
-  MSAVI() : m_S(0.4)
-  {
-    m_WDVIfunctor.SetS(m_S);
-  }
-  ~MSAVI() override {}
-/** Set/Get Slop of soil line */
-  void SetS(const double s)
-  {
-    m_S = s;
-    m_WDVIfunctor.SetS(m_S);
-  }
-  double GetS(void) const
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    return (m_S);
-  }
-  NDVIFunctorType GetNDVI(void) const
-  {
-    return (m_NDVIfunctor);
-  }
-  WDVIFunctorType GetWDVI(void) const
-  {
-    return (m_WDVIfunctor);
-  }
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
 
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
-  {
-    double dnir = static_cast<double>(nir);
-    double dr = static_cast<double>(r);
+    double ndvi = NDVI<TInput, TOutput>::Compute(red, nir);
+    double wdvi = WDVI<TInput, TOutput>::Compute(red, nir);
 
-    double dNDVI = this->GetNDVI() (r, nir);
-    double dWDVI = this->GetWDVI() (r, nir);
-    double dL = 1 - 2 * m_S * dNDVI * dWDVI;
+    double L = 1 - 2 * S * ndvi * wdvi;
 
-    double denominator = dnir + dr + dL;
+    double denominator = nir + red + L;
 
-    if (std::abs(denominator)  < this->m_EpsilonToBeConsideredAsZero)
-      {
+    if (std::abs(denominator) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       return static_cast<TOutput>(0.);
       }
 
-    return (static_cast<TOutput>(((dnir - dr) * (1 + dL)) / denominator));
+      return (static_cast<TOutput>(((nir - red) * (1 + L)) / denominator));
   }
 
 private:
   /** Slope of soil line */
-  double                m_S;
-  NDVIFunctorType m_NDVIfunctor;
-  WDVIFunctorType       m_WDVIfunctor;
-
+  static constexpr double S = 0.4;
 };
 
 /** \class MSAVI2
@@ -798,30 +297,25 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class MSAVI2 : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class MSAVI2 : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  MSAVI2() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "MSAVI2";
   }
 
-  MSAVI2() {}
-  ~MSAVI2() override {}
-
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    double dnir = static_cast<double>(nir);
-    double dr = static_cast<double>(r);
-    double sqrt_value = (2 * dnir + 1) * (2 * dnir + 1) - 8 * (dnir - dr);
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
+
+    double sqrt_value = (2 * nir + 1) * (2 * nir + 1) - 8 * (nir - red);
     if (sqrt_value < 0.)
       {
       return static_cast<TOutput>(0.);
       }
-    return (static_cast<TOutput>((2 * dnir + 1 - std::sqrt(sqrt_value)) / 2.));
+      return (static_cast<TOutput>((2 * nir + 1 - std::sqrt(sqrt_value)) / 2.));
   }
 
 };
@@ -836,44 +330,39 @@ protected:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class GEMI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class GEMI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  GEMI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "GEMI";
   }
 
-  GEMI() {}
-  ~GEMI() override {}
-
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    double dnir = static_cast<double>(nir);
-    double dr = static_cast<double>(r);
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
 
-    double dnu;
-    double dnumerateur_nu;
-    double ddenominateur_nu = dnir + dr + 0.5;
-    if (std::abs(ddenominateur_nu)  < this->m_EpsilonToBeConsideredAsZero)
-      {
-      dnu = 0;
+    double nu;
+    double num_nu;
+    double denom_nu = nir + red + 0.5;
+
+    if (std::abs(denom_nu) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
+      nu = 0;
       }
     else
       {
-      dnumerateur_nu = 2 * (dnir * dnir - dr * dr) + 1.5 * dnir + 0.5 * dr;
-      dnu = dnumerateur_nu / ddenominateur_nu;
+        num_nu = 2 * (nir * nir - red * red) + 1.5 * nir + 0.5 * red;
+        nu     = num_nu / denom_nu;
       }
 
-    double ddenominateur_GEMI = 1 - dr;
-    if (std::abs(ddenominateur_GEMI)  < this->m_EpsilonToBeConsideredAsZero)
+      double denom_GEMI = 1 - red;
+      if (std::abs(denom_GEMI) < RadiometricIndex<TInput, TOutput>::Epsilon)
       {
       return static_cast<TOutput>(0.);
       }
-    return (static_cast<TOutput>((dnu * (1 - 0.25 * dnu) - (dr - 0.125)) / ddenominateur_GEMI));
+      return (static_cast<TOutput>((nu * (1 - 0.25 * nu) - (red - 0.125)) / denom_GEMI));
   }
 
 };
@@ -890,87 +379,54 @@ protected:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TInput3, class TOutput>
-class AVI : public RAndGAndNIRIndexBase<TInput1, TInput2, TInput3, TOutput>
+template <class TInput, class TOutput>
+class AVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  AVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::GREEN, CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "AVI";
   }
 
-  AVI() : m_LambdaG(560.), m_LambdaR(660.), m_LambdaNir(830.) {}
-  ~AVI() override {}
-/** Set/Get Lambda red parameter*/
-  void SetLambdaR(const double lr)
-  {
-    m_LambdaR = lr;
-  }
-  double GetLambdaR(void) const
-  {
-    return (m_LambdaR);
-  }
-/** Set/Get Lambda green parameter */
-  void SetLambdaG(const double lg)
-  {
-    m_LambdaG = lg;
-  }
-  double GetLambdaG(void) const
-  {
-    return (m_LambdaG);
-  }
-/** Set/Get Lambda red parameter */
-  void SetLambdaNir(const double lnir)
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    m_LambdaNir = lnir;
-  }
-  double GetLambdaNir(void) const
-  {
-    return (m_LambdaNir);
-  }
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& g, const TInput3& nir) const override
-  {
-    double dr = static_cast<double>(r);
-    double dg = static_cast<double>(g);
-    double dnir = static_cast<double>(nir);
+    auto green = this->Value(CommonBandNames::GREEN, input);
+    auto red   = this->Value(CommonBandNames::RED, input);
+    auto nir   = this->Value(CommonBandNames::NIR, input);
 
-    double dfact1 = (m_LambdaNir - m_LambdaR) / m_LambdaR;
-    double dfact2 = (m_LambdaR - m_LambdaG) / m_LambdaR;
+    constexpr double dfact1 = (LambdaNir - LambdaR) / LambdaR;
+    constexpr double dfact2 = (LambdaR - LambdaG) / LambdaR;
     double dterm1;
     double dterm2;
-    if (std::abs(dnir - dr)  < this->m_EpsilonToBeConsideredAsZero)
-      {
+    if (std::abs(nir - red) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       dterm1 = 0;
       }
     else
       {
-      dterm1 = std::atan(dfact1 / (dnir - dr));
+        dterm1 = std::atan(dfact1 / (nir - red));
       }
 
-    if (std::abs(dg - dr)  < this->m_EpsilonToBeConsideredAsZero)
+      if (std::abs(green - red) < RadiometricIndex<TInput, TOutput>::Epsilon)
       {
       dterm2 = 0;
       }
     else
       {
-      dterm2 = std::atan(dfact2 / (dg - dr));
+        dterm2 = std::atan(dfact2 / (green - red));
       }
 
     return static_cast<TOutput>(dterm1 + dterm2);
 
   }
-private:
 
   /**  Central wavelength of the green channel (=Lambda1) */
-  double m_LambdaG;
+  static constexpr double LambdaG = 560;
 
   /**  Central wavelength of the red channel (=Lambda2) */
-  double m_LambdaR;
+  static constexpr double LambdaR = 660;
 
   /**  Central wavelength of the nir channel (=Lambda3) */
-  double m_LambdaNir;
+  static constexpr double LambdaNir = 830;
 };
 
 /** \class ARVI
@@ -985,134 +441,31 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TInput3, class TOutput>
-class ARVI : public RAndBAndNIRIndexBase<TInput1, TInput2, TInput3, TOutput>
-{
-public:
-  /** Return the index name */
-  std::string GetName() const override
-  {
-    return "ARVI";
-  }
-
-  ARVI() : m_Gamma(0.5) {}
-  ~ARVI() override {}
-
-  /** Set/Get Gamma parameter */
-  void SetGamma(const double gamma)
-  {
-    m_Gamma = gamma;
-  }
-  double GetGamma(void) const
-  {
-    return (m_Gamma);
-  }
-
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& b, const TInput3& nir) const override
-  {
-    double dr = static_cast<double>(r);
-    double db = static_cast<double>(b);
-    double dnir = static_cast<double>(nir);
-    double RHOrb = dr - m_Gamma * (db - dr);
-    double denominator = dnir + RHOrb;
-    if (std::abs(denominator)  < this->m_EpsilonToBeConsideredAsZero)
-      {
-      return static_cast<TOutput>(0.);
-      }
-    return (static_cast<TOutput>((dnir - RHOrb) / denominator));
-  }
-
-private:
-
-  /** Gamma parameter */
-  double m_Gamma;
-};
-
-/** \class TSARVI
- *  \brief This functor computes the Transformed Soil Atmospherical Resistant Vegetation Index (TSARVI)
- *
- *  [Yoram J. Kaufman and Didier Tanre, 1992]
- *
- *  \ingroup Functor
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template <class TInput1, class TInput2, class TInput3, class TOutput>
-class TSARVI : public RAndBAndNIRIndexBase<TInput1, TInput2, TInput3, TOutput>
+template <class TInput, class TOutput>
+class ARVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  ARVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::BLUE, CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "TSARVI";
   }
 
-  TSARVI() : m_A(0.0), m_B(0.0), m_X(0.08), m_Gamma(0.5) {}
-  ~TSARVI() override {}
-
-  /** Set/Get A and B parameters */
-  void SetA(const double A)
-  {
-    m_A = A;
-  }
-  double GetA(void) const
-  {
-    return (m_A);
-  }
-  void SetB(const double B)
-  {
-    m_B = B;
-  }
-  double GetB(void) const
-  {
-    return (m_B);
-  }
-  /** Set/Get X parameter */
-  void SetX(const double X)
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    m_X = X;
-  }
-  double GetX(void) const
-  {
-    return (m_X);
-  }
-  /** Set/Get the gamma parameter */
-  void SetGamma(const double gamma)
-  {
-    m_Gamma = gamma;
-  }
-  double GetGamma(void) const
-  {
-    return (m_Gamma);
-  }
+    auto blue = this->Value(CommonBandNames::BLUE, input);
+    auto red  = this->Value(CommonBandNames::RED, input);
+    auto nir  = this->Value(CommonBandNames::NIR, input);
 
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& b, const TInput3& nir) const override
-  {
-    double dr = static_cast<double>(r);
-    double db = static_cast<double>(b);
-    double dnir = static_cast<double>(nir);
-    double dRB = dr - m_Gamma * (db - dr);
-    double denominator = dRB + m_A * dnir - m_A * m_B + m_X * (1. + m_A * m_A);
-    if (std::abs(denominator)  < this->m_EpsilonToBeConsideredAsZero)
-      {
+    double RHOrb       = red - Gamma * (blue - red);
+    double denominator = nir + RHOrb;
+    if (std::abs(denominator) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       return static_cast<TOutput>(0.);
       }
-    return (static_cast<TOutput>((m_A * (dnir - m_A * dRB - m_B)) / denominator));
+      return (static_cast<TOutput>((nir - RHOrb) / denominator));
   }
 
-private:
-
-  /** A and B parameters */
-  double m_A;
-  double m_B;
-  /** X parameter */
-  double m_X;
   /** Gamma parameter */
-  double m_Gamma;
-
+  static constexpr double Gamma = 0.5;
 };
 
 /** \class EVI
@@ -1127,81 +480,39 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TInput3, class TOutput>
-class EVI : public RAndBAndNIRIndexBase<TInput1, TInput2, TInput3, TOutput>
+template <class TInput, class TOutput>
+class EVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  EVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::BLUE, CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "EVI";
   }
 
-  EVI() : m_G(2.5), m_C1(6.0), m_C2(7.5), m_L(1.0) {}
-  ~EVI() override {}
-/** Set/Get G parameter */
-  void SetG(const double g)
-  {
-    m_G = g;
-  }
-  double GetG(void) const
-  {
-    return (m_G);
-  }
-  /** Set/Get C1 parameter */
-  void SetC1(const double c1)
-  {
-    m_C1 = c1;
-  }
-  double GetC1(void) const
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    return (m_C1);
-  }
-  /** Set/Get C2 parameter */
-  void SetC2(const double c2)
-  {
-    m_C2 = c2;
-  }
-  double GetC2(void) const
-  {
-    return (m_C2);
-  }
-  /** Set/Get L parameter */
-  void SetL(const double l)
-  {
-    m_L = l;
-  }
-  double GetL(void) const
-  {
-    return (m_L);
-  }
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& b, const TInput3& nir) const override
-  {
-    double dr = static_cast<double>(r);
-    double db = static_cast<double>(b);
-    double dnir = static_cast<double>(nir);
-    double denominator = dnir + m_C1 * dr - m_C2 * db + m_L;
-    if (std::abs(denominator) < this->m_EpsilonToBeConsideredAsZero)
-      {
+    auto blue = this->Value(CommonBandNames::BLUE, input);
+    auto red  = this->Value(CommonBandNames::RED, input);
+    auto nir  = this->Value(CommonBandNames::NIR, input);
+
+    double denominator = nir + C1 * red - C2 * blue + L;
+    if (std::abs(denominator) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       return (static_cast<TOutput>(0.));
       }
-    return (static_cast<TOutput>(m_G * (dnir - dr) / denominator));
+      return (static_cast<TOutput>(G * (nir - red) / denominator));
   }
 
-private:
-
   /** Gain factor */
-  double m_G;
+  static constexpr double G = 2.5;
 
   /** Coefficient of the aerosol resistance term */
-  double m_C1;
+  static constexpr double C1 = 6.0;
 
   /** Coefficient of the aerosol resistance term */
-  double m_C2;
+  static constexpr double C2 = 7.5;
 
   /** Canopy background adjustment */
-  double m_L;
+  static constexpr double L = 1.0;
 };
 
 /** \class IPVI
@@ -1214,31 +525,26 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class IPVI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class IPVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  IPVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "IPVI";
   }
 
-  IPVI() {}
-  ~IPVI() override {}
-
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    double dr = static_cast<double>(r);
-    double dnir = static_cast<double>(nir);
-    if (std::abs(dnir + dr)  < this->m_EpsilonToBeConsideredAsZero)
-      {
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
+
+    if (std::abs(nir + red) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
       return static_cast<TOutput>(0.);
       }
     else
       {
-      return (static_cast<TOutput>(dnir / (dnir + dr)));
+        return (static_cast<TOutput>(nir / (nir + red)));
       }
   }
 };
@@ -1253,40 +559,30 @@ protected:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class TNDVI : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class TNDVI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  TNDVI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-    return "TNDVI";
   }
 
-  typedef NDVI<TInput1, TInput2, TOutput> NDVIFunctorType;
-  TNDVI() {}
-  ~TNDVI() override {}
-
-  NDVIFunctorType GetNDVI(void) const
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    return (m_NDVIfunctor);
-  }
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
 
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
-  {
-    double dval = this->GetNDVI() (r, nir) + 0.5;
-    if (dval < 0)
-      {
+    double val = NDVI<TInput, TOutput>::Compute(red, nir) + 0.5;
+
+    if (val < 0)
+    {
       return  (static_cast<TOutput>(0));
       }
     else
       {
-      return (static_cast<TOutput>(std::sqrt(dval)));
+        return (static_cast<TOutput>(std::sqrt(val)));
       }
   }
-private:
-  NDVIFunctorType m_NDVIfunctor;
 };
 
 /** \class LAIFromNDVILogarithmic
@@ -1305,69 +601,62 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class LAIFromNDVILogarithmic : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class LAIFromNDVILogarithmic : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  LAIFromNDVILogarithmic()
+    : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR}), m_NdviSoil(0.1), m_NdviInf(0.89), m_ExtinctionCoefficient(0.71)
   {
-    return "LAIFromNDVILogarithmic";
   }
 
-  typedef NDVI<TInput1, TInput2, TOutput> NDVIFunctorType;
-  LAIFromNDVILogarithmic() : m_NdviSoil(0.10), m_NdviInf(0.89), m_ExtinctionCoefficient(0.71) {}
-  ~LAIFromNDVILogarithmic() override {}
-
-  NDVIFunctorType GetNDVI(void) const
-  {
-    return (m_NDVIfunctor);
-  }
-
-  void SetNdviSoil(const double val)
+  void SetNdviSoil(const double& val)
   {
     m_NdviSoil = val;
   }
-  double GetNdviSoil(void) const
+
+  const double& GetNdviSoil() const
   {
-    return (m_NdviSoil);
+    return m_NdviSoil;
   }
 
-  void SetNdviInf(const double val)
+  void SetNdviInf(const double& val)
   {
     m_NdviInf = val;
   }
-  double GetNdviInf(void) const
+
+  const double& GetNdviInf() const
   {
-    return (m_NdviInf);
+    return m_NdviInf;
   }
 
-  void SetExtinctionCoefficient(const double val)
+  void SetExtinctionCoefficient(const double& val)
   {
     m_ExtinctionCoefficient = val;
   }
-  double GetExtinctionCoefficient(void) const
+
+  const double& GetExtionctionCoefficient() const
   {
-    return (m_ExtinctionCoefficient);
+    return m_ExtinctionCoefficient;
   }
 
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    double dval = this->GetNDVI() (r, nir);
-    if (dval < 0)
-      {
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
+
+    double val = NDVI<TInput, TOutput>::Compute(red, nir);
+
+    if (val < 0)
+    {
       return  (static_cast<TOutput>(0));
       }
     else
       {
-      return (static_cast<TOutput>(
-    -(1.0/m_ExtinctionCoefficient)*std::log((dval- m_NdviInf)/(m_NdviSoil-m_NdviInf))
-    ));
+        return static_cast<TOutput>(-(1.0 / m_ExtinctionCoefficient) * std::log((val - m_NdviInf) / (m_NdviSoil - m_NdviInf)));
       }
   }
-private:
-  NDVIFunctorType m_NDVIfunctor;
+
   double m_NdviSoil;
   double m_NdviInf;
   double m_ExtinctionCoefficient;
@@ -1391,50 +680,42 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class LAIFromReflectancesLinear : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class LAIFromReflectancesLinear : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
-  {
-    return "LAIFromReflectancesLinear";
-  }
-
-  typedef NDVI<TInput1, TInput2, TOutput> NDVIFunctorType;
-  LAIFromReflectancesLinear() : m_RedCoef(-17.91), m_NirCoef(12.26) {}
-  ~LAIFromReflectancesLinear() override {}
-
-  NDVIFunctorType GetReflectances(void) const
+  LAIFromReflectancesLinear() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR}), m_RedCoef(-17.91), m_NirCoef(12.26)
   {
-    return (m_NDVIfunctor);
   }
 
-  void SetRedCoef(const double val)
+  void SetRedCoef(const double& val)
   {
     m_RedCoef = val;
   }
-  double GetRedCoef(void) const
+
+  const double& GetRedCoef() const
   {
-    return (m_RedCoef);
+    return m_RedCoef;
   }
 
-  void SetNirCoef(const double val)
+  void SetNirCoef(const double& val)
   {
     m_NirCoef = val;
   }
-  double GetNirCoef(void) const
+
+  const double& GetNirCoef() const
   {
-    return (m_NirCoef);
+    return m_NirCoef;
   }
 
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-      return (static_cast<TOutput>(m_RedCoef*r+m_NirCoef*nir));
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
+
+    return (static_cast<TOutput>(m_RedCoef * red + m_NirCoef * nir));
   }
-private:
-  NDVIFunctorType m_NDVIfunctor;
+
   double m_RedCoef;
   double m_NirCoef;
 };
@@ -1459,40 +740,29 @@ private:
   */
 
 
-  template <class TInput1, class TInput2, class TOutput>
-  class LAIFromNDVIFormosat2Functor : public RAndNIRIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class LAIFromNDVIFormosat2Functor : public RadiometricIndex<TInput, TOutput>
+{
+public:
+  LAIFromNDVIFormosat2Functor() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::NIR})
   {
-  public:
+  }
 
-    /** Return the index name */
-    std::string GetName() const override
-    {
-      return "LAIFromNDVIFormosat2Functor";
-    }
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
+  {
+    auto red = this->Value(CommonBandNames::RED, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
 
-    /// Constructor
-    LAIFromNDVIFormosat2Functor() {}
-    /// Desctructor
-    ~LAIFromNDVIFormosat2Functor() override {}
-    // Operator on r and nir single pixel values
-  protected:
-    inline TOutput Evaluate(const TInput1& r, const TInput2& nir) const override
+    if (std::abs(nir + red) < RadiometricIndex<TInput, TOutput>::Epsilon)
     {
-    double a = 0.1519;
-    double b = 3.9443;
-    double c = 0.13;
-
-      double dr = static_cast<double>(r);
-      double dnir = static_cast<double>(nir);
-      if (std::abs(dnir + dr) < this->m_EpsilonToBeConsideredAsZero)
-        {
-        return static_cast<TOutput>(0.);
-        }
-
-      return  static_cast<TOutput>(a*(std::exp(static_cast<double>(dnir-dr)/static_cast<double>(dr+dnir)*b)-std::exp(c*b)));
-    };
-
+      return static_cast<TOutput>(0.);
+    }
+    return static_cast<TOutput>(A * (std::exp((nir - red) / (red + nir) * B) - std::exp(C * B)));
+  }
 
+  static constexpr double A = 0.1519;
+  static constexpr double B = 3.9443;
+  static constexpr double C = 0.13;
 };
 
 
diff --git a/Modules/Radiometry/Indices/include/otbWaterIndicesFunctor.h b/Modules/Radiometry/Indices/include/otbWaterIndicesFunctor.h
index e68d884d721fbdbe9fe4cbf67b2f4c034c8e66c3..19acf0fb19bd6681451935869625f2b1267d0665 100644
--- a/Modules/Radiometry/Indices/include/otbWaterIndicesFunctor.h
+++ b/Modules/Radiometry/Indices/include/otbWaterIndicesFunctor.h
@@ -22,168 +22,12 @@
 #define otbWaterIndicesFunctor_h
 
 #include "otbMath.h"
-#include "itkVariableLengthVector.h"
-#include "otbSqrtSpectralAngleFunctor.h"
-#include "otbBandName.h"
-#include <string>
+#include "otbRadiometricIndex.h"
 
 namespace otb
 {
 namespace Functor
 {
-/**
-   * \class WaterIndexBase
-   * \brief Base class
-   *
-   *  Implement operators for UnaryFunctorImageFilter templated with a
-   *  VectorImage and BinaryFunctorImageFilter templated with single
-   *  images.
-   *  Subclasses should NOT overload operators, they must  re-implement
-   *  the Evaluate() method.
-   *
-   * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template<class TInput1, class TInput2, class TOutput>
-class WaterIndexBase
-{
-public:
-  /// Vector pixel type used to support both vector images and multiple
-  /// input images
-  typedef itk::VariableLengthVector<TInput1> InputVectorType;
-
-  //operators !=
-  bool operator !=(const WaterIndexBase&) const
-  {
-    return true;
-  }
-  //operator ==
-  bool operator ==(const WaterIndexBase& other) const
-  {
-    return !(*this != other);
-  }
-
-  // Operator on vector pixel type
-  inline TOutput operator ()(const InputVectorType& inputVector) const
-  {
-    return this->Evaluate(inputVector[m_Index1 - 1], static_cast<TInput2>(inputVector[m_Index2 - 1]));
-  }
-
-  // Binary operator
-  inline TOutput operator ()(const TInput1& id1, const TInput2& id2) const
-  {
-    return this->Evaluate(id1, id2);
-  }
-  /// Constructor
-  WaterIndexBase() {}
-  /// Desctructor
-  virtual ~WaterIndexBase() {}
-
-  /// Set Index 1
-  void SetIndex1(unsigned int channel)
-  {
-    m_Index1 = channel;
-  }
-  /// Get Index 1
-  unsigned int GetIndex1() const
-  {
-    return m_Index1;
-  }
-  /// Set Index 2
-  void SetIndex2(unsigned int channel)
-  {
-    m_Index2 = channel;
-  }
-  /// Get Index 2
-  unsigned int GetIndex2() const
-  {
-    return m_Index2;
-  }
-
-  /** Return the index name */
-  virtual std::string GetName() const = 0;
-
-protected:
-  // This method must be reimplemented in subclasses to actually
-  // compute the index value
-  virtual TOutput Evaluate(const TInput1& id1, const TInput2& id2) const = 0;
-
-private:
-  unsigned int m_Index1;
-  unsigned int m_Index2;
-};
-
-/** \class WaterIndexFunctor
- *  \brief This functor will be used for most of water index functors.
- *
- *  \ingroup Functor
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template <class TInput1, class TInput2, class TOutput>
-class WaterIndexFunctor : public WaterIndexBase<TInput1, TInput2, TOutput>
-{
-public:
-  /** Return the index name */
-  std::string GetName() const override
-  {
-    return "WaterIndexFunctor";
-  }
-
-  WaterIndexFunctor() {}
-  ~WaterIndexFunctor() override {}
-protected:
-  inline TOutput Evaluate(const TInput1& id1, const TInput2& id2) const override
-  {
-    double dindex1 = static_cast<double>(id1);
-    double dindex2 = static_cast<double>(id2);
-    double ddenom = dindex1 + dindex2;
-    if (ddenom == 0)
-      {
-      return static_cast<TOutput>(0.);
-      }
-    return (static_cast<TOutput>((dindex1 - dindex2) / ddenom));
-  }
-};
-
-/** \class SRWI
- *  \brief This functor computes the Simple Ratio Water Index (SRWI)
- *  \brief For MODIS bands 860 & 1240
- *
- *   [Zarco-Tejada 2001]
- *
- *  \ingroup Functor
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template <class TInput1, class TInput2, class TOutput>
-class SRWI : public WaterIndexBase<TInput1, TInput2, TOutput>
-{
-public:
-  /** Return the index name */
-  virtual std::string GetName() const
-  {
-    return "SRWI";
-  }
-
-  SRWI() {}
-  virtual ~SRWI() {}
-protected:
-  inline TOutput Evaluate(const TInput1& rho860, const TInput2& rho1240) const
-  {
-    double drho860 = static_cast<double>(rho860);
-    double drho1240 = static_cast<double>(rho1240);
-    if (drho1240 == 0)
-      {
-      return static_cast<TOutput>(0.);
-      }
-    return (static_cast<TOutput>(drho860 / drho1240));
-  }
-};
-
 /** \class NDWI
  *  \brief This functor computes the Normalized Difference Water Index (NDWI)
  *  \brief Also called :
@@ -198,79 +42,26 @@ protected:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class NDWI : public WaterIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class NDWI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  NDWI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::NIR, CommonBandNames::MIR})
   {
-    return "NDWI";
   }
 
-  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
-  /// Constructor
-  NDWI() {}
-  /// Desctructor
-  ~NDWI() override {}
-  WIFunctorType GetWIFunctor(void) const
-  {
-    return (m_WIFunctor);
-  }
-  /// Set Index NIR
-  void SetNIRIndex(unsigned int channel)
-  {
-    this->SetIndex1(channel);
-  }
-  /// Get Index NIR
-  unsigned int GetNIRIndex() const
-  {
-    return this->GetIndex1();
-  }
-  /// Set Index MIR
-  void SetMIRIndex(unsigned int channel)
-  {
-    this->SetIndex2(channel);
-  }
-  /// Get Index MIR
-  unsigned int GetMIRIndex() const
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    return this->GetIndex2();
-  }
+    auto mir = this->Value(CommonBandNames::MIR, input);
+    auto nir = this->Value(CommonBandNames::NIR, input);
 
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::NIR)
-      {
-      this->SetIndex1(channel);
-      }
-    if (band == BandName::MIR)
-      {
-      this->SetIndex2(channel);
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::NIR)
-      {
-      return this->GetIndex1();
+    if (std::abs(nir + mir) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
+      return 0.;
       }
-    if (band == BandName::MIR)
-      {
-      return this->GetIndex2();
-      }
-  }
 
-protected:
-  inline TOutput Evaluate(const TInput1& nir, const TInput2& mir) const override
-  {
-    return (static_cast<TOutput>(GetWIFunctor() (nir, mir)));
+      return (nir - mir) / (nir + mir);
   }
-private:
-  // Water Index Classic Functor
-  WIFunctorType m_WIFunctor;
 };
 
 /** \class NDWI2
@@ -283,79 +74,26 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class NDWI2 : public WaterIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class NDWI2 : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  std::string GetName() const override
+  NDWI2() : RadiometricIndex<TInput, TOutput>({CommonBandNames::NIR, CommonBandNames::GREEN})
   {
-    return "NDWI2";
   }
 
-  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
-  /// Constructor
-  NDWI2() {}
-  /// Desctructor
-  ~NDWI2() override {}
-  WIFunctorType GetWIFunctor(void) const
-  {
-    return (m_WIFunctor);
-  }
-  /// Set Index G
-  void SetGIndex(unsigned int channel)
-  {
-    this->SetIndex1(channel);
-  }
-  /// Get Index G
-  unsigned int GetGIndex() const
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    return this->GetIndex1();
-  }
-  /// Set Index NIR
-  void SetNIRIndex(unsigned int channel)
-  {
-    this->SetIndex2(channel);
-  }
-  /// Get Index NIR
-  unsigned int GetNIRIndex() const
-  {
-    return this->GetIndex2();
-  }
+    auto green = this->Value(CommonBandNames::GREEN, input);
+    auto nir   = this->Value(CommonBandNames::NIR, input);
 
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::GREEN)
-      {
-      this->SetIndex1(channel);
-      }
-    if (band == BandName::NIR)
-      {
-      this->SetIndex2(channel);
+    if (std::abs(nir + green) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
+      return 0.;
       }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::GREEN)
-      {
-      return this->GetIndex1();
-      }
-    if (band == BandName::NIR)
-      {
-      return this->GetIndex2();
-      }
-  }
 
-protected:
-  inline TOutput Evaluate(const TInput1& g, const TInput2& nir) const override
-  {
-    return (static_cast<TOutput>(GetWIFunctor() (g, nir)));
+      return (green - nir) / (green + nir);
   }
-private:
-  // Water Index Classic Functor
-  WIFunctorType m_WIFunctor;
 };
 
 /** \class MNDWI
@@ -363,88 +101,7 @@ private:
  *
  *  [Xu & al., 2006 ]
  *
- *  \ingroup Functor
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template <class TInput1, class TInput2, class TOutput>
-class MNDWI : public WaterIndexBase<TInput1, TInput2, TOutput>
-{
-public:
-  /** Return the index name */
-  virtual std::string GetName() const
-  {
-    return "MNDWI";
-  }
-
-  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
-  /// Constructor
-  MNDWI() {}
-  /// Desctructor
-  virtual ~MNDWI() {}
-  WIFunctorType GetWIFunctor(void) const
-  {
-    return (m_WIFunctor);
-  }
-  /// Set Index G
-  void SetGIndex(unsigned int channel)
-  {
-    this->SetIndex1(channel);
-  }
-  /// Get Index G
-  unsigned int GetGIndex() const
-  {
-    return this->GetIndex1();
-  }
-  /// Set Index MIR
-  void SetMIRIndex(unsigned int channel)
-  {
-    this->SetIndex2(channel);
-  }
-  /// Get Index MIR
-  unsigned int GetMIRIndex() const
-  {
-    return this->GetIndex2();
-  }
-
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::GREEN)
-      {
-      this->SetIndex1(channel);
-      }
-    if (band == BandName::MIR)
-      {
-      this->SetIndex2(channel);
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::GREEN)
-      {
-      return this->GetIndex1();
-      }
-    if (band == BandName::MIR)
-      {
-      return this->GetIndex2();
-      }
-  }
-
-protected:
-  inline TOutput Evaluate(const TInput1& g, const TInput2& mir) const
-  {
-    return (static_cast<TOutput>(GetWIFunctor() (g, mir)));
-  }
-private:
-  // Water Index Classic Functor
-  WIFunctorType m_WIFunctor;
-};
-
-/** \class NDPI
- *  \brief This functor computes the Normalized Difference Pond Index (NDPI)
+ * Similar to Normalized Difference Pond Index (NDPI)
  *
  *  [J.P Lacaux & al., 2006 ]
  *
@@ -453,79 +110,26 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class NDPI : public WaterIndexBase<TInput1, TInput2, TOutput>
+template <class TInput, class TOutput>
+class MNDWI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  virtual std::string GetName() const
+  MNDWI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::MIR, CommonBandNames::GREEN})
   {
-    return "NDPI";
   }
 
-  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
-  /// Constructor
-  NDPI() {}
-  /// Desctructor
-  virtual ~NDPI() {}
-  WIFunctorType GetWIFunctor(void) const
-  {
-    return (m_WIFunctor);
-  }
-  /// Set Index MIR
-  void SetMIRIndex(unsigned int channel)
-  {
-    this->SetIndex1(channel);
-  }
-  /// Get Index MIR
-  unsigned int GetMIRIndex() const
-  {
-    return this->GetIndex1();
-  }
-  /// Set Index G
-  void SetGIndex(unsigned int channel)
-  {
-    this->SetIndex2(channel);
-  }
-  /// Get Index G
-  unsigned int GetGIndex() const
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    return this->GetIndex2();
-  }
+    auto green = this->Value(CommonBandNames::GREEN, input);
+    auto mir   = this->Value(CommonBandNames::MIR, input);
 
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::MIR)
-      {
-      this->SetIndex1(channel);
-      }
-    if (band == BandName::GREEN)
-      {
-      this->SetIndex2(channel);
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::MIR)
-      {
-      return this->GetIndex1();
-      }
-    if (band == BandName::GREEN)
-      {
-      return this->GetIndex2();
+    if (std::abs(mir + green) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
+      return 0.;
       }
-  }
 
-protected:
-  inline TOutput Evaluate(const TInput1& mir, const TInput2& g) const
-  {
-    return (static_cast<TOutput>(GetWIFunctor() (mir, g)));
+      return (green - mir) / (green + mir);
   }
-private:
-  // Water Index Classic Functor
-  WIFunctorType m_WIFunctor;
 };
 
 /** \class NDTI
@@ -538,221 +142,26 @@ private:
  *
  * \ingroup OTBIndices
  */
-template <class TInput1, class TInput2, class TOutput>
-class NDTI : public WaterIndexBase<TInput1, TInput2, TOutput>
-{
-public:
-  /** Return the index name */
-  virtual std::string GetName() const
-  {
-    return "NDTI";
-  }
-
-  typedef WaterIndexFunctor<TInput1, TInput2, TOutput> WIFunctorType;
-  /// Constructor
-  NDTI() {}
-  /// Desctructor
-  virtual ~NDTI() {}
-  WIFunctorType GetWIFunctor(void) const
-  {
-    return (m_WIFunctor);
-  }
-  // FIXME why now using Red and Green fully spelled as everywhere
-  //else ???
-  /// Set Index R
-  void SetRIndex(unsigned int channel)
-  {
-    this->SetIndex1(channel);
-  }
-  /// Get Index R
-  unsigned int GetRIndex() const
-  {
-    return this->GetIndex1();
-  }
-  /// Set Index G
-  void SetGIndex(unsigned int channel)
-  {
-    this->SetIndex2(channel);
-  }
-  /// Get Index G
-  unsigned int GetGIndex() const
-  {
-    return this->GetIndex2();
-  }
-
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::RED)
-      {
-      this->SetIndex1(channel);
-      }
-    if (band == BandName::GREEN)
-      {
-      this->SetIndex2(channel);
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::RED)
-      {
-      return this->GetIndex1();
-      }
-    if (band == BandName::GREEN)
-      {
-      return this->GetIndex2();
-      }
-  }
-
-protected:
-  inline TOutput Evaluate(const TInput1& r, const TInput2& g) const
-  {
-    return (static_cast<TOutput>(GetWIFunctor() (r, g)));
-  }
-private:
-  // Water Index Classic Functor
-  WIFunctorType m_WIFunctor;
-};
-
-/** \class WaterSqrtSpectralAngleFunctor
- *  \brief This functor uses a spectral angle with a particular reference pixel.
- *
- *
- *  \ingroup Functor
- * \ingroup Radiometry
- *
- * \ingroup OTBIndices
- */
-template <class TInputVectorPixel, class TOutputPixel>
-class WaterSqrtSpectralAngleFunctor : public SqrtSpectralAngleFunctor<TInputVectorPixel, TOutputPixel>
+template <class TInput, class TOutput>
+class NDTI : public RadiometricIndex<TInput, TOutput>
 {
 public:
-  /** Return the index name */
-  virtual std::string GetName() const
+  NDTI() : RadiometricIndex<TInput, TOutput>({CommonBandNames::RED, CommonBandNames::GREEN})
   {
-    return "WaterSqrtSpectralAngleFunctor";
-  }
-  typedef WaterSqrtSpectralAngleFunctor                             Self;
-  typedef SqrtSpectralAngleFunctor<TInputVectorPixel, TOutputPixel> Superclass;
-  typedef TInputVectorPixel                                         InputVectorPixelType;
-  WaterSqrtSpectralAngleFunctor()
-  {
-
-    //Set the channels indices
-    m_BlueIndex = 0;
-    m_GreenIndex = 1;
-    m_RedIndex = 2;
-    m_NIRIndex = 3;
-
-    //Set reference water value
-    InputVectorPixelType reference;
-    reference.SetSize(4);
-    reference[0] = 136.0; reference[1] = 132.0; reference[2] = 47.0; reference[3] = 24.0;
-    this->SetReferenceWaterPixel(reference);
   }
-  ~WaterSqrtSpectralAngleFunctor() override {}
 
-  /** Set Reference Pixel */
-  void SetReferenceWaterPixel(InputVectorPixelType ref)
+  TOutput operator()(const itk::VariableLengthVector<TInput>& input) const override
   {
-    if (ref.GetSize() != 4)
-      {
-      }
-    InputVectorPixelType reference;
-    reference.SetSize(4);
-    reference[m_BlueIndex] = ref[0]; reference[m_GreenIndex] = ref[1]; reference[m_RedIndex] = ref[2];
-    reference[m_NIRIndex] = ref[3];
-    this->SetReferencePixel(reference);
-
-  }
+    auto green = this->Value(CommonBandNames::GREEN, input);
+    auto red   = this->Value(CommonBandNames::RED, input);
 
-  /** Getters and setters */
-  void SetBlueChannel(unsigned int channel)
-  {
-    m_BlueIndex = channel;
-  }
-  unsigned int GetBlueChannel() const
-  {
-    return m_BlueIndex;
-  }
-  void SetGreenChannel(unsigned int channel)
-  {
-    m_GreenIndex = channel;
-  }
-  unsigned int GetGreenChannel() const
-  {
-    return m_GreenIndex;
-  }
-  void SetRedChannel(unsigned int channel)
-  {
-    m_RedIndex = channel;
-  }
-  unsigned int GetRedChannel() const
-  {
-    return m_RedIndex;
-  }
-  void SetNIRChannel(unsigned int channel)
-  {
-    m_NIRIndex = channel;
-  }
-  unsigned int GetNIRChannel() const
-  {
-    return m_NIRIndex;
-  }
-
-  /** Set index, generic method */
-  void SetIndex(BandName::BandName band, unsigned int channel)
-  {
-    if (band == BandName::RED)
-      {
-      m_RedIndex = channel;
-      }
-    if (band == BandName::GREEN)
-      {
-      m_GreenIndex = channel;
-      }
-    if (band == BandName::BLUE)
-      {
-      m_BlueIndex = channel;
-      }
-    if (band == BandName::NIR)
-      {
-      m_NIRIndex = channel;
-      }
-  }
-  /** Get index, generic method */
-  unsigned int GetIndex(BandName::BandName band) const
-  {
-    if (band == BandName::RED)
-      {
-      return m_RedIndex;
+    if (std::abs(red + green) < RadiometricIndex<TInput, TOutput>::Epsilon)
+    {
+      return 0.;
       }
-    if (band == BandName::GREEN)
-      {
-      return m_GreenIndex;
-      }
-    if (band == BandName::BLUE)
-      {
-      return m_BlueIndex;
-      }
-    if (band == BandName::NIR)
-      {
-      return m_NIRIndex;
-      }
-  }
 
-protected:
-  inline TOutputPixel Evaluate(const TInputVectorPixel& inPix) const override
-  {
-    return static_cast<TOutputPixel>(Superclass::Evaluate(inPix));
+      return (red - green) / (green + red);
   }
-
-  /** Channels */
-  int m_BlueIndex;
-  int m_GreenIndex;
-  int m_RedIndex;
-  int m_NIRIndex;
 };
 
 } // namespace Functor
diff --git a/Modules/Radiometry/Indices/include/otbWaterSqrtSpectralAngleImageFilter.h b/Modules/Radiometry/Indices/include/otbWaterSqrtSpectralAngleImageFilter.h
index 0ac02c7d6fec990ecbf17b7f981dad2e00603300..adbe0dee543242d0156ac2aafb22444f32b52c60 100644
--- a/Modules/Radiometry/Indices/include/otbWaterSqrtSpectralAngleImageFilter.h
+++ b/Modules/Radiometry/Indices/include/otbWaterSqrtSpectralAngleImageFilter.h
@@ -21,12 +21,116 @@
 #ifndef otbWaterSqrtSpectralAngleImageFilter_h
 #define otbWaterSqrtSpectralAngleImageFilter_h
 
-#include "otbWaterIndicesFunctor.h"
+#include "otbSqrtSpectralAngleFunctor.h"
 #include "itkUnaryFunctorImageFilter.h"
 
 namespace otb
 {
 
+namespace Functor
+{
+/** \class WaterSqrtSpectralAngleFunctor
+ *  \brief This functor uses a spectral angle with a particular reference pixel.
+ *
+ *
+ *  \ingroup Functor
+ * \ingroup Radiometry
+ *
+ * \ingroup OTBIndices
+ */
+template <class TInputVectorPixel, class TOutputPixel>
+class WaterSqrtSpectralAngleFunctor : public SqrtSpectralAngleFunctor<TInputVectorPixel, TOutputPixel>
+{
+public:
+  typedef WaterSqrtSpectralAngleFunctor                             Self;
+  typedef SqrtSpectralAngleFunctor<TInputVectorPixel, TOutputPixel> Superclass;
+  typedef TInputVectorPixel                                         InputVectorPixelType;
+  WaterSqrtSpectralAngleFunctor()
+  {
+
+    // Set the channels indices
+    m_BlueIndex  = 0;
+    m_GreenIndex = 1;
+    m_RedIndex   = 2;
+    m_NIRIndex   = 3;
+
+    // Set reference water value
+    InputVectorPixelType reference;
+    reference.SetSize(4);
+    reference[0] = 136.0;
+    reference[1] = 132.0;
+    reference[2] = 47.0;
+    reference[3] = 24.0;
+    this->SetReferenceWaterPixel(reference);
+  }
+  ~WaterSqrtSpectralAngleFunctor() override
+  {
+  }
+
+  /** Set Reference Pixel */
+  void SetReferenceWaterPixel(InputVectorPixelType ref)
+  {
+    if (ref.GetSize() != 4)
+    {
+    }
+    InputVectorPixelType reference;
+    reference.SetSize(4);
+    reference[m_BlueIndex]  = ref[0];
+    reference[m_GreenIndex] = ref[1];
+    reference[m_RedIndex]   = ref[2];
+    reference[m_NIRIndex]   = ref[3];
+    this->SetReferencePixel(reference);
+  }
+
+  /** Getters and setters */
+  void SetBlueChannel(unsigned int channel)
+  {
+    m_BlueIndex = channel;
+  }
+  unsigned int GetBlueChannel() const
+  {
+    return m_BlueIndex;
+  }
+  void SetGreenChannel(unsigned int channel)
+  {
+    m_GreenIndex = channel;
+  }
+  unsigned int GetGreenChannel() const
+  {
+    return m_GreenIndex;
+  }
+  void SetRedChannel(unsigned int channel)
+  {
+    m_RedIndex = channel;
+  }
+  unsigned int GetRedChannel() const
+  {
+    return m_RedIndex;
+  }
+  void SetNIRChannel(unsigned int channel)
+  {
+    m_NIRIndex = channel;
+  }
+  unsigned int GetNIRChannel() const
+  {
+    return m_NIRIndex;
+  }
+
+protected:
+  inline TOutputPixel Evaluate(const TInputVectorPixel& inPix) const override
+  {
+    return static_cast<TOutputPixel>(Superclass::Evaluate(inPix));
+  }
+
+  /** Channels */
+  int m_BlueIndex;
+  int m_GreenIndex;
+  int m_RedIndex;
+  int m_NIRIndex;
+};
+} // End namespace Functor
+
+
 /** \class WaterSqrtSpectralAngleImageFilter
  *  \brief Compute a radiometric water indice
  *
@@ -46,7 +150,6 @@ namespace otb
  *
  * \ingroup OTBIndices
  */
-
 template <class TInputVectorImage, class TOutputImage,
     class TFunction = Functor::WaterSqrtSpectralAngleFunctor <
         typename TInputVectorImage::PixelType,
diff --git a/Modules/Radiometry/Indices/test/CMakeLists.txt b/Modules/Radiometry/Indices/test/CMakeLists.txt
index 3e5c8354670345f66bd113cf0c312d76dbcc88c9..79e4a828212956c5a5a7aa28977d4767115aa396 100644
--- a/Modules/Radiometry/Indices/test/CMakeLists.txt
+++ b/Modules/Radiometry/Indices/test/CMakeLists.txt
@@ -22,13 +22,12 @@ otb_module_test()
 
 set(OTBIndicesTests
 otbIndicesTestDriver.cxx
+otbRadiometricIndicesTest.cxx
 otbNDVIDataNodeFeatureFunction.cxx
 otbLandsatTMIndexNDSITest.cxx
 otbLandsatTMIndexBIOTest.cxx
-otbLAIFromReflectancesLinearFunctorTest.cxx
 otbLandsatTMIndexMIR2Test.cxx
 otbLandsatTMIndexNDVITest.cxx
-otbLAIFromNDVILogarithmicFunctorTest.cxx
 otbLandsatTMIndexVisTest.cxx
 otbWaterSqrtSpectralAngleImageFilter.cxx
 otbLandsatTMIndexBrightTest.cxx
@@ -102,14 +101,6 @@ otb_add_test(NAME raTvLandsatTMIndexBIOTest COMMAND otbIndicesTestDriver
   21  #TM7
   )
 
-otb_add_test(NAME raTvLAIFromReflectancesLinearFunctorTest COMMAND otbIndicesTestDriver
-  otbLAIFromReflectancesLinear
-  3   # red
-  4   # nir
-  -18   # red coef
-  13   # nir coef
-  )
-
 otb_add_test(NAME raTvLandsatTMIndexMIR2Test COMMAND otbIndicesTestDriver
   otbLandsatTMIndexMIR2
   3   #TM1
@@ -134,16 +125,6 @@ otb_add_test(NAME raTvLandsatTMIndexNDVITest COMMAND otbIndicesTestDriver
   21  #TM7
   )
 
-
-otb_add_test(NAME raTvLAIFromNDVILogarithmicFunctorTest COMMAND otbIndicesTestDriver
-  otbLAIFromNDVILogarithmic
-  3   # red
-  4   # nir
-  0.12   # ndvi soil
-  0.91   # ndvi infinity
-  0.70   # extinction coefficient
-  )
-
 otb_add_test(NAME raTvLandsatTMIndexVisTest COMMAND otbIndicesTestDriver
   otbLandsatTMIndexVis
   3   #TM1
@@ -340,3 +321,21 @@ otb_add_test(NAME raTvLandsatTMThickCloudTest COMMAND otbIndicesTestDriver
   ${TEMP}/raTvLandsatTMThickCloudTest_cloudImage.tif
   )
 
+
+otb_add_test(NAME raTvRadiometricIndexBaseClassTest COMMAND otbIndicesTestDriver
+                  otbRadiometricIndexTest)
+
+otb_add_test(NAME raTvVegetationIndicesTest COMMAND otbIndicesTestDriver
+                  otbVegetationIndicesTest)
+
+otb_add_test(NAME raTvWaterIndicesTest COMMAND otbIndicesTestDriver
+                  otbWaterIndicesTest)
+
+otb_add_test(NAME raTvBuiltUpIndicesTest COMMAND otbIndicesTestDriver
+                  otbBuiltUpIndicesTest)
+
+otb_add_test(NAME raTvSoilIndicesTest COMMAND otbIndicesTestDriver
+                  otbSoilIndicesTest)
+
+otb_add_test(NAME raTvIndicesStackFunctorTest COMMAND otbIndicesTestDriver
+                  otbIndicesStackFunctorTest)
diff --git a/Modules/Radiometry/Indices/test/otbIndicesTestDriver.cxx b/Modules/Radiometry/Indices/test/otbIndicesTestDriver.cxx
index d0b41d61e7d623c2690731cbc478ddb24253cdb8..884ebeabb10bac520f40d10336bd96aec25b5b93 100644
--- a/Modules/Radiometry/Indices/test/otbIndicesTestDriver.cxx
+++ b/Modules/Radiometry/Indices/test/otbIndicesTestDriver.cxx
@@ -25,10 +25,8 @@ void RegisterTests()
   REGISTER_TEST(otbNDVIDataNodeFeatureFunction);
   REGISTER_TEST(otbLandsatTMIndexNDSI);
   REGISTER_TEST(otbLandsatTMIndexBIO);
-  REGISTER_TEST(otbLAIFromReflectancesLinear);
   REGISTER_TEST(otbLandsatTMIndexMIR2);
   REGISTER_TEST(otbLandsatTMIndexNDVI);
-  REGISTER_TEST(otbLAIFromNDVILogarithmic);
   REGISTER_TEST(otbLandsatTMIndexVis);
   REGISTER_TEST(otbWaterSqrtSpectralAngleImageFilter);
   REGISTER_TEST(otbLandsatTMIndexBright);
@@ -46,4 +44,10 @@ void RegisterTests()
   REGISTER_TEST(otbLandsatTMKernelSpectralRulesWithImage);
   REGISTER_TEST(otbLandsatTMIndexNDBSI);
   REGISTER_TEST(otbLandsatTMThickCloudTest);
+  REGISTER_TEST(otbVegetationIndicesTest);
+  REGISTER_TEST(otbWaterIndicesTest);
+  REGISTER_TEST(otbBuiltUpIndicesTest);
+  REGISTER_TEST(otbSoilIndicesTest);
+  REGISTER_TEST(otbRadiometricIndexTest);
+  REGISTER_TEST(otbIndicesStackFunctorTest);
 }
diff --git a/Modules/Radiometry/Indices/test/otbLAIFromNDVILogarithmicFunctorTest.cxx b/Modules/Radiometry/Indices/test/otbLAIFromNDVILogarithmicFunctorTest.cxx
deleted file mode 100644
index 09a2f5ff71f4393797077b67c485e3b3777e5a45..0000000000000000000000000000000000000000
--- a/Modules/Radiometry/Indices/test/otbLAIFromNDVILogarithmicFunctorTest.cxx
+++ /dev/null
@@ -1,58 +0,0 @@
-/*
- * Copyright (C) 2005-2019 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 "otbVegetationIndicesFunctor.h"
-
-int otbLAIFromNDVILogarithmic(int itkNotUsed(argc), char * argv[])
-{
-  typedef double                           PixelType;
-
-  typedef otb::Functor::LAIFromNDVILogarithmic<PixelType, PixelType, PixelType> FunctorType;
-
-  FunctorType laiFunct = FunctorType();
-
-  double redValue = (::atof(argv[1]));
-  double nirValue = (::atof(argv[2]));
-  double ndviSoil(::atof(argv[3]));
-  double ndviInf(::atof(argv[4]));
-  double extCoef(::atof(argv[5]));
-
-  double ndvi = (nirValue-redValue)/(nirValue+redValue);
-  double goodResult = -1/extCoef*std::log((ndvi-ndviInf)/(ndviSoil-ndviInf));
-
-  laiFunct.SetNdviInf(ndviInf);
-  laiFunct.SetNdviSoil(ndviSoil);
-  laiFunct.SetExtinctionCoefficient(extCoef);
-
-  laiFunct.SetRedIndex(1);
-  laiFunct.SetNIRIndex(2);
-
-  itk::VariableLengthVector<PixelType> pixel;
-  pixel.Reserve(2);
-  pixel[0] = redValue;
-  pixel[1] = nirValue;
-
-  double result = laiFunct(pixel);
-
-  if( result!=goodResult ) return EXIT_FAILURE;
-
-  return EXIT_SUCCESS;
-}
diff --git a/Modules/Radiometry/Indices/test/otbLAIFromReflectancesLinearFunctorTest.cxx b/Modules/Radiometry/Indices/test/otbLAIFromReflectancesLinearFunctorTest.cxx
deleted file mode 100644
index 388b36a054b9054d3b362c29e148ee2bbd01e81f..0000000000000000000000000000000000000000
--- a/Modules/Radiometry/Indices/test/otbLAIFromReflectancesLinearFunctorTest.cxx
+++ /dev/null
@@ -1,55 +0,0 @@
-/*
- * Copyright (C) 2005-2019 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 "otbVegetationIndicesFunctor.h"
-
-int otbLAIFromReflectancesLinear(int itkNotUsed(argc), char * argv[])
-{
-  typedef double                           PixelType;
-
-  typedef otb::Functor::LAIFromReflectancesLinear<PixelType, PixelType, PixelType> FunctorType;
-
-  FunctorType laiFunct = FunctorType();
-
-  double redValue = (::atof(argv[1]));
-  double nirValue = (::atof(argv[2]));
-  double redCoef(::atof(argv[3]));
-  double nirCoef(::atof(argv[4]));
-
-  double goodResult = redCoef*redValue+nirCoef*nirValue;
-
-  laiFunct.SetRedCoef(redCoef);
-  laiFunct.SetNirCoef(nirCoef);
-
-  laiFunct.SetRedIndex(1);
-  laiFunct.SetNIRIndex(2);
-
-  itk::VariableLengthVector<PixelType> pixel;
-  pixel.Reserve(2);
-  pixel[0] = redValue;
-  pixel[1] = nirValue;
-
-  double result = laiFunct(pixel);
-
-  if( result!=goodResult ) return EXIT_FAILURE;
-
-  return EXIT_SUCCESS;
-}
diff --git a/Modules/Radiometry/Indices/test/otbRadiometricIndicesTest.cxx b/Modules/Radiometry/Indices/test/otbRadiometricIndicesTest.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..58f76ed98817ba7a2f33e291eb85f31a6e5edccd
--- /dev/null
+++ b/Modules/Radiometry/Indices/test/otbRadiometricIndicesTest.cxx
@@ -0,0 +1,279 @@
+/*
+ * Copyright (C) 2005-2019 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 "otbVegetationIndicesFunctor.h"
+#include "otbWaterIndicesFunctor.h"
+#include "otbBuiltUpIndicesFunctor.h"
+#include "otbSoilIndicesFunctor.h"
+#include "otbIndicesStackFunctor.h"
+
+#include <iomanip>
+
+template <typename T>
+itk::VariableLengthVector<T> build_pixel(const std::initializer_list<T>& il)
+{
+  itk::VariableLengthVector<T> res(il.size());
+  size_t                       idx = 0;
+
+  for (auto v : il)
+  {
+    res[idx] = v;
+    ++idx;
+  }
+  return res;
+}
+
+template <class TIndice>
+bool CheckResult(const std::string& testName, std::map<typename TIndice::BandNameType, size_t> bandMap,
+                 const std::initializer_list<typename TIndice::InputType>& input, const typename TIndice::OutputType& expected)
+{
+  TIndice indice;
+
+  indice.SetBandsIndices(bandMap);
+
+  auto pixel = build_pixel(input);
+
+  typename TIndice::OutputType v = indice(pixel);
+
+  if (std::abs(expected - v) > TIndice::Epsilon)
+  {
+    std::cerr << std::setprecision(10);
+    std::cerr << testName << "\t- failed: expected " << expected << ", got " << v << std::endl;
+    return false;
+  }
+  else
+  {
+    return true;
+  }
+}
+
+
+using namespace otb::Functor;
+
+int otbVegetationIndicesTest(int, char ** const)
+{
+
+  const std::map<CommonBandNames, size_t> bandMap = {
+      {CommonBandNames::BLUE, 1}, {CommonBandNames::GREEN, 2}, {CommonBandNames::RED, 3}, {CommonBandNames::NIR, 4}};
+
+  // Syntax: CheckResult<Indice Class>("test_name",bandMap,{red_value,nir_value},expected_result)
+  bool res = CheckResult<NDVI<int, double>>("ndvi_null ", bandMap, {0, 0, 0, 0}, 0.);
+  res      = res & CheckResult<NDVI<int, double>>("ndvi_pixel", bandMap, {0, 0, 1, 2}, 0.3333333);
+  res      = res & CheckResult<RVI<int, double>>("rvi_null", bandMap, {0, 0, 0, 0}, 0.);
+  res      = res & CheckResult<RVI<int, double>>("rvi_pixel", bandMap, {0, 0, 2, 1}, 0.5);
+  res      = res & CheckResult<PVI<int, double>>("pvi_pixel", bandMap, {0, 0, 1, 2}, -62.0544166);
+  res      = res & CheckResult<SAVI<double, double>>("savi_null", bandMap, {0, 0, 0, -0.5}, 0.);
+  res      = res & CheckResult<SAVI<int, double>>("savi_pixel", bandMap, {0, 0, 1, 2}, 0.42857142857);
+  res      = res & CheckResult<TSAVI<double, double>>("tsavi_null", bandMap, {0, 0, -0.1192, 0.}, 0.);
+  res      = res & CheckResult<TSAVI<int, double>>("tsavi_pixel", bandMap, {0, 0, 1, 2}, 0.1111463957);
+  res      = res & CheckResult<WDVI<int, double>>("wdvi_pixel", bandMap, {0, 0, 1, 2}, 1.6);
+  res      = res & CheckResult<MSAVI<int, double>>("msavi_pixel", bandMap, {0, 0, 1, 2}, 0.4402985075);
+  res      = res & CheckResult<MSAVI2<int, double>>("msavi2_pixel", bandMap, {0, 0, 1, 2}, 0.4384471872);
+  res      = res & CheckResult<GEMI<int, double>>("gemi_pixel", bandMap, {1, 4, 3, 2}, 2.0625);
+  res      = res & CheckResult<AVI<int, double>>("avi_pixel", bandMap, {0, 0, 1, 2}, 0.1017245527);
+  res      = res & CheckResult<ARVI<int, double>>("arvi_pixel", bandMap, {0, 0, 1, 2}, 0.1428571429);
+  res      = res & CheckResult<EVI<int, double>>("evi_pixel", bandMap, {0, 0, 1, 2}, 0.2777777778);
+  res      = res & CheckResult<IPVI<int, double>>("ipvi_pixel", bandMap, {0, 0, 1, 2}, 0.6666666667);
+  res      = res & CheckResult<LAIFromNDVILogarithmic<int, double>>("lailog_pixel", bandMap, {0, 0, 1, 2}, 0.4930511672);
+  res      = res & CheckResult<LAIFromReflectancesLinear<int, double>>("lailog_pixel", bandMap, {0, 0, 1, 2}, 6.61);
+  res      = res & CheckResult<LAIFromNDVIFormosat2Functor<int, double>>("laifrom_pixel", bandMap, {0, 0, 1, 2}, 0.3120010659);
+
+  if (res)
+  {
+    return EXIT_SUCCESS;
+  }
+  else
+  {
+    return EXIT_FAILURE;
+  }
+}
+
+int otbWaterIndicesTest(int, char ** const)
+{
+  const std::map<CommonBandNames, size_t> bandMap = {
+      {CommonBandNames::BLUE, 1}, {CommonBandNames::GREEN, 2}, {CommonBandNames::RED, 3}, {CommonBandNames::NIR, 4}, {CommonBandNames::MIR, 5}};
+
+  // Syntax: CheckResult<Indice Class>("test_name",bandMap,{red_value,nir_value},expected_result)
+  bool res = CheckResult<NDWI<int, double>>("ndwi_null ", bandMap, {0, 0, 0, 0, 0}, 0.);
+  res      = res & CheckResult<NDWI<int, double>>("ndwi_pixel", bandMap, {1, 2, 3, 4, 5}, -0.1111111111);
+  res      = res & CheckResult<NDWI2<int, double>>("ndwi2_null", bandMap, {0, 0, 0, 0, 0}, 0.);
+  res      = res & CheckResult<NDWI2<int, double>>("ndwi2_pixel", bandMap, {1, 2, 3, 4, 5}, -0.3333333333);
+  res      = res & CheckResult<MNDWI<int, double>>("mndwi_null", bandMap, {0, 0, 0, 0, 0}, 0.);
+  res      = res & CheckResult<MNDWI<int, double>>("mndwi_pixel", bandMap, {1, 2, 3, 4, 5}, -0.4285714286);
+  res      = res & CheckResult<NDTI<int, double>>("ndti_null", bandMap, {0, 0, 0, 0, 0}, 0.);
+  res      = res & CheckResult<NDTI<int, double>>("ndti_pixel", bandMap, {1, 2, 3, 4, 5}, 0.2);
+
+  if (res)
+  {
+    return EXIT_SUCCESS;
+  }
+  else
+  {
+    return EXIT_FAILURE;
+  }
+}
+
+int otbSoilIndicesTest(int, char ** const)
+{
+  const std::map<CommonBandNames, size_t> bandMap = {
+      {CommonBandNames::BLUE, 1}, {CommonBandNames::GREEN, 2}, {CommonBandNames::RED, 3}, {CommonBandNames::NIR, 4}, {CommonBandNames::MIR, 5}};
+
+  // Syntax: CheckResult<Indice Class>("test_name",bandMap,{red_value,nir_value},expected_result)
+  bool res = CheckResult<CI<int, double>>("ci_null ", bandMap, {0, 0, 0, 0, 0}, 0.);
+  res      = res & CheckResult<CI<int, double>>("ci_pixel", bandMap, {1, 2, 3, 4, 5}, 0.2);
+  res      = res & CheckResult<BI<int, double>>("bi_pixel", bandMap, {1, 2, 3, 4, 5}, 2.549509757);
+  res      = res & CheckResult<BI2<int, double>>("bi2_pixel", bandMap, {1, 2, 3, 4, 5}, 3.109126351);
+
+  if (res)
+  {
+    return EXIT_SUCCESS;
+  }
+  else
+  {
+    return EXIT_FAILURE;
+  }
+}
+
+int otbBuiltUpIndicesTest(int, char ** const)
+{
+  const std::map<CommonBandNames, size_t> bandMap = {
+      {CommonBandNames::BLUE, 1}, {CommonBandNames::GREEN, 2}, {CommonBandNames::RED, 3}, {CommonBandNames::NIR, 4}, {CommonBandNames::MIR, 5}};
+
+  // Syntax: CheckResult<Indice Class>("test_name",bandMap,{red_value,nir_value},expected_result)
+  bool res = CheckResult<ISU<int, double>>("isu_null", bandMap, {0, 0, 0, 0, 0}, 0.);
+  res      = res & CheckResult<ISU<int, double>>("isu_pixel", bandMap, {1, 2, 3, 4, 5}, 81.25);
+
+  if (res)
+  {
+    return EXIT_SUCCESS;
+  }
+  else
+  {
+    return EXIT_FAILURE;
+  }
+}
+
+
+int otbRadiometricIndexTest(int, char ** const)
+{
+  auto ndvi = NDVI<double, double>();
+
+  auto requiredBands = ndvi.GetRequiredBands();
+
+  bool success = true;
+
+  if (requiredBands.size() != 2 || requiredBands.find(CommonBandNames::RED) == requiredBands.end() ||
+      requiredBands.find(CommonBandNames::NIR) == requiredBands.end())
+  {
+    std::cerr << "Required bands is not {RED,NIR} for NDVI" << std::endl;
+    success = false;
+  }
+
+  ndvi.SetBandIndex(CommonBandNames::RED, 10);
+
+  if (ndvi.GetBandIndex(CommonBandNames::RED) != 10)
+  {
+    std::cerr << "Could not Set/Get band index properly" << std::endl;
+    success = false;
+  }
+
+  const std::map<CommonBandNames, size_t> bandMap = {{CommonBandNames::RED, 100}, {CommonBandNames::NIR, 200}};
+
+  ndvi.SetBandsIndices(bandMap);
+
+  if (ndvi.GetBandIndex(CommonBandNames::RED) != 100 || ndvi.GetBandIndex(CommonBandNames::NIR) != 200)
+  {
+    std::cerr << "Could not set all band indices at once with SetBandIndices" << std::endl;
+    success = false;
+  }
+
+  try
+  {
+    ndvi.SetBandIndex(CommonBandNames::MAX, 1);
+    std::cerr << "Calling SetBandIndices with ::MAX should raise a runtime_error exception." << std::endl;
+    success = false;
+  }
+  catch (const std::runtime_error& e)
+  {
+  }
+
+  if (success)
+  {
+    return EXIT_SUCCESS;
+  }
+  else
+  {
+    return EXIT_FAILURE;
+  }
+}
+
+int otbIndicesStackFunctorTest(int, char ** const)
+{
+  using IndicesType      = RadiometricIndex<double, int>;
+  using StackFunctorType = IndicesStackFunctor<IndicesType>;
+
+  auto ndvi = NDVI<double, int>();
+  auto ndwi = NDWI<double, int>();
+
+  std::vector<IndicesType*> indices = {&ndvi, &ndwi};
+
+  auto stack = StackFunctorType(indices);
+
+  bool success = true;
+
+  if (stack.OutputSize() != 2)
+  {
+    std::cerr << "Size of output pixel for stack functor should be 2" << std::endl;
+    success = false;
+  }
+
+  const std::map<CommonBandNames, size_t> bandMap = {
+      {CommonBandNames::BLUE, 1}, {CommonBandNames::GREEN, 2}, {CommonBandNames::RED, 3}, {CommonBandNames::NIR, 4}, {CommonBandNames::MIR, 5}};
+
+  ndvi.SetBandsIndices(bandMap);
+  ndwi.SetBandsIndices(bandMap);
+
+  StackFunctorType::OutputType out(2);
+
+  auto in = build_pixel<double>({1, 2, 3, 4, 5});
+
+  stack(out, in);
+
+  if (out[0] != ndvi(in))
+  {
+    std::cerr << "First output band should correspond to ndvi" << std::endl;
+    success = false;
+  }
+
+  if (out[1] != ndwi(in))
+  {
+    std::cerr << "Second output band should correspond to ndwi" << std::endl;
+    success = false;
+  }
+
+  if (success)
+  {
+    return EXIT_SUCCESS;
+  }
+  else
+  {
+    return EXIT_FAILURE;
+  }
+}
diff --git a/Modules/Radiometry/Simulation/include/otbReduceSpectralResponseClassifierRAndNIR.h b/Modules/Radiometry/Simulation/include/otbReduceSpectralResponseClassifierRAndNIR.h
index ec0ef3d8264028aaa656959abc693b19f4eec0ae..ea82717f97b18260a1b5e281937b879903607081 100644
--- a/Modules/Radiometry/Simulation/include/otbReduceSpectralResponseClassifierRAndNIR.h
+++ b/Modules/Radiometry/Simulation/include/otbReduceSpectralResponseClassifierRAndNIR.h
@@ -52,47 +52,46 @@ namespace otb
  *
  * \ingroup OTBSimulation
  */
-  template <class TReduceSpectralResponse , class TFunction = Functor::NDVI< typename TReduceSpectralResponse::ValuePrecisionType, typename TReduceSpectralResponse::ValuePrecisionType,
-  typename TReduceSpectralResponse::ValuePrecisionType > >
-      class ReduceSpectralResponseClassifierRAndNIR
-  : public itk::DataObject
-      {
-        //friend class
-        public:
-          /** Standard class typedefs */
-          typedef ReduceSpectralResponseClassifierRAndNIR Self;
-          typedef itk::DataObject Superclass;
-          typedef itk::SmartPointer<Self> Pointer;
-          typedef itk::SmartPointer<const Self> ConstPointer;
-
-          /** Template parameters typedef */
-          typedef TReduceSpectralResponse InputReduceSpectralResponseType;
-          typedef TFunction   FunctorType;
-          typedef typename TReduceSpectralResponse::Pointer InputReduceSpectralResponsePointerType;
-          typedef typename InputReduceSpectralResponseType::ValuePrecisionType ValuePrecisionType;
-
-
-          /** Standard macros */
-          itkNewMacro(Self);
-          itkTypeMacro(ReduceSpectralResponseClassifierRAndNIR, DataObject);
-
-          itkGetConstObjectMacro(InputReduceSpectralResponse, InputReduceSpectralResponseType);
-          itkSetObjectMacro(InputReduceSpectralResponse, InputReduceSpectralResponseType);
-
-          itkGetConstMacro(RBandNumber, unsigned int);
-          itkSetMacro(RBandNumber, unsigned int);
-
-          itkGetConstMacro(NIRBandNumber, unsigned int);
-          itkSetMacro(NIRBandNumber, unsigned int);
-
-          /** Get the functor object.  The functor is returned by reference.
-           * (Functors do not have to derive from itk::LightObject, so they do
-           * not necessarily have a reference count. So we cannot return a
-           * SmartPointer.) */
-          FunctorType& GetFunctor()
-          {
-            return m_Functor;
-          };
+template <class TReduceSpectralResponse,
+          class TFunction = Functor::NDVI<typename TReduceSpectralResponse::ValuePrecisionType, typename TReduceSpectralResponse::ValuePrecisionType>>
+class ReduceSpectralResponseClassifierRAndNIR : public itk::DataObject
+{
+  // friend class
+public:
+  /** Standard class typedefs */
+  typedef ReduceSpectralResponseClassifierRAndNIR Self;
+  typedef itk::DataObject                         Superclass;
+  typedef itk::SmartPointer<Self>                 Pointer;
+  typedef itk::SmartPointer<const Self>           ConstPointer;
+
+  /** Template parameters typedef */
+  typedef TReduceSpectralResponse                                      InputReduceSpectralResponseType;
+  typedef TFunction                                                    FunctorType;
+  typedef typename TReduceSpectralResponse::Pointer                    InputReduceSpectralResponsePointerType;
+  typedef typename InputReduceSpectralResponseType::ValuePrecisionType ValuePrecisionType;
+
+
+  /** Standard macros */
+  itkNewMacro(Self);
+  itkTypeMacro(ReduceSpectralResponseClassifierRAndNIR, DataObject);
+
+  itkGetConstObjectMacro(InputReduceSpectralResponse, InputReduceSpectralResponseType);
+  itkSetObjectMacro(InputReduceSpectralResponse, InputReduceSpectralResponseType);
+
+  itkGetConstMacro(RBandNumber, unsigned int);
+  itkSetMacro(RBandNumber, unsigned int);
+
+  itkGetConstMacro(NIRBandNumber, unsigned int);
+  itkSetMacro(NIRBandNumber, unsigned int);
+
+  /** Get the functor object.  The functor is returned by reference.
+   * (Functors do not have to derive from itk::LightObject, so they do
+   * not necessarily have a reference count. So we cannot return a
+   * SmartPointer.) */
+  FunctorType& GetFunctor()
+  {
+    return m_Functor;
+  };
 
   /** Set the functor object.  This replaces the current Functor with a
            * copy of the specified Functor. This allows the user to specify a
diff --git a/Modules/Radiometry/Simulation/include/otbReduceSpectralResponseClassifierRAndNIR.hxx b/Modules/Radiometry/Simulation/include/otbReduceSpectralResponseClassifierRAndNIR.hxx
index 4d97fb2b292b4d8f62359b9497f19f4908de680c..bd43ed41d45b7c073210ead4443c3b1643dc35a2 100644
--- a/Modules/Radiometry/Simulation/include/otbReduceSpectralResponseClassifierRAndNIR.hxx
+++ b/Modules/Radiometry/Simulation/include/otbReduceSpectralResponseClassifierRAndNIR.hxx
@@ -50,7 +50,13 @@ namespace otb
       ReduceSpectralResponseClassifierRAndNIR<TReduceSpectralResponse , TFunction>
   ::operator()()
   {
-      return m_Functor((*m_InputReduceSpectralResponse)(m_RBandNumber), (*m_InputReduceSpectralResponse)(m_NIRBandNumber));
+    itk::VariableLengthVector<ValuePrecisionType> sr(2);
+    sr[0]=(*m_InputReduceSpectralResponse)(m_RBandNumber);
+    sr[1]=(*m_InputReduceSpectralResponse)(m_NIRBandNumber);
+    m_Functor.SetBandIndex(CommonBandNames::RED,1);
+    m_Functor.SetBandIndex(CommonBandNames::NIR,2);
+
+    return m_Functor(sr);
   }
 
 
diff --git a/Modules/Radiometry/Simulation/test/otbReduceSpectralResponseClassifierRAndNIR.cxx b/Modules/Radiometry/Simulation/test/otbReduceSpectralResponseClassifierRAndNIR.cxx
index 87fda203e2af13b01998c836307f0376f3199b7d..a0e0d61e169419bdc67eabee18589f0353b4f44f 100644
--- a/Modules/Radiometry/Simulation/test/otbReduceSpectralResponseClassifierRAndNIR.cxx
+++ b/Modules/Radiometry/Simulation/test/otbReduceSpectralResponseClassifierRAndNIR.cxx
@@ -45,7 +45,7 @@ int otbReduceSpectralResponseClassifierRAndNIR(int argc, char * argv[])
   typedef otb::ReduceSpectralResponse < ResponseType, SatRSRType>  ReduceResponseType;
   typedef ReduceResponseType::Pointer  ReduceResponseTypePointerType;
 
-  typedef otb::Functor::NDVI<double, double, double >               TFunctionType;
+  typedef otb::Functor::NDVI<double, double>                                               TFunctionType;
   typedef otb::ReduceSpectralResponseClassifierRAndNIR <ReduceResponseType, TFunctionType> ReduceSpectralResponseClassifierRAndNIRType;
   typedef ReduceSpectralResponseClassifierRAndNIRType::Pointer  ReduceSpectralResponseClassifierRAndNIRPointerType;
 
diff --git a/Modules/ThirdParty/OssimPlugins/src/ossim/ossimSarSensorModel.cpp b/Modules/ThirdParty/OssimPlugins/src/ossim/ossimSarSensorModel.cpp
index e57304fb7751a2ff0e330741cc179a2cb75c81d9..5819e87018a93914463a6dab2304002ce6563e1f 100644
--- a/Modules/ThirdParty/OssimPlugins/src/ossim/ossimSarSensorModel.cpp
+++ b/Modules/ThirdParty/OssimPlugins/src/ossim/ossimSarSensorModel.cpp
@@ -1058,29 +1058,39 @@ bool ossimSarSensorModel::worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeT
 
       theAzimuthTimeOffset = count > 0 ? cumulAzimuthTime / count : DurationType(0);
 
-      // Then, fix the range time
-      count=0;
+      // Patch_locS1 : do not fix the range time
+      // This fix takes GCPs as inputs to calculate an offset in range dimension. 
+      // However for recent S1 products (after March 2017 with IPF >= 2.82), GCPs contain wrong values 
+      // for lat/lon coordinates => the range time offset is wrong and shift the localization.
+      // To avoid wrong offset, fixRangeTimeWithGCPs is always set to false
+      bool fixRangeTimeWithGCPs = false;
+
+      if (fixRangeTimeWithGCPs)
+	{
+	  // Then, fix the range time
+	  count=0;
 
-      for(std::vector<GCPRecordType>::const_iterator gcpIt = theGCPRecords.begin(); gcpIt!=theGCPRecords.end();++gcpIt)
-      {
-         ossimDpt estimatedImPt;
-         TimeType estimatedAzimuthTime;
-         double   estimatedRangeTime;
+	  for(std::vector<GCPRecordType>::const_iterator gcpIt = theGCPRecords.begin(); gcpIt!=theGCPRecords.end();++gcpIt)
+	  {
+	     ossimDpt estimatedImPt;
+	     TimeType estimatedAzimuthTime;
+	     double   estimatedRangeTime;
 
-         ossimEcefPoint sensorPos;
-         ossimEcefVector sensorVel;
+	     ossimEcefPoint sensorPos;
+	     ossimEcefVector sensorVel;
 
-         // Estimate times
-         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime, sensorPos, sensorVel);
+	     // Estimate times
+	     const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime, sensorPos, sensorVel);
 
-         if(s1)
-         {
-            cumulRangeTime+=-estimatedRangeTime+gcpIt->slantRangeTime;
-            ++count;
-         }
-      }
+	     if(s1)
+	     {
+	        cumulRangeTime+=-estimatedRangeTime+gcpIt->slantRangeTime;
+	        ++count;
+	     }
+	  }
 
-      theRangeTimeOffset = count > 0 ? cumulRangeTime/count : 0;
+	  theRangeTimeOffset = count > 0 ? cumulRangeTime/count : 0;
+	}
    }
 
    void get(
diff --git a/Modules/Wrappers/QtWidget/src/otbWrapperQtWidgetSpinBoxes.cxx b/Modules/Wrappers/QtWidget/src/otbWrapperQtWidgetSpinBoxes.cxx
index 4af61804517d82518ce13d78ec021a86b044fee5..9984222ae6a1fd7240dc544e2d9a811d15e245ce 100644
--- a/Modules/Wrappers/QtWidget/src/otbWrapperQtWidgetSpinBoxes.cxx
+++ b/Modules/Wrappers/QtWidget/src/otbWrapperQtWidgetSpinBoxes.cxx
@@ -103,7 +103,8 @@ void QtWidgetSpinBox::SetValueNoSignal(int value)
 int QtWidgetSpinBox::valueFromText(const QString &text) const
 {
   bool ok;
-  int result = QLocale::system().toInt(text, &ok);
+  // Force C locale because OTB gui is not i18n
+  int result = QLocale::c().toInt(text, &ok);
   if (ok)
   {
     return result;
@@ -162,7 +163,8 @@ void QtWidgetDoubleSpinBox::SetValueNoSignal(double value)
 double QtWidgetDoubleSpinBox::valueFromText(const QString &text) const
 {
   bool ok;
-  double result = QLocale::system().toDouble(text, &ok);
+  // Force C locale because OTB gui is not i18n
+  double result = QLocale::c().toDouble(text, &ok);
   if (ok)
   {
     return result;
@@ -180,7 +182,9 @@ QString QtWidgetDoubleSpinBox::textFromValue(double value) const
   // which leads to ugly trailing zeros for small values (e.g 1.50000)
   // We use std::ostringstream because QString::arg formatting support is too limited
   std::ostringstream oss;
-  oss.imbue(std::locale("")); // use system's locale for formatting
+
+  // Force C locale because OTB gui is not i18n
+  oss.imbue(std::locale::classic());
 
   // Set precision to the number of decimal digits that can be represented without change.
   // Use float precision because OTB parameter is float
@@ -190,9 +194,8 @@ QString QtWidgetDoubleSpinBox::textFromValue(double value) const
 
   // Add a trailing dot if the number is integer,
   // so that int and float parameters are more visually different.
-  // For now this is done for all locales, even though not all locales use this
-  // convention for formatting decimals...
-  const char dot = std::use_facet<std::numpunct<char>>(std::locale("")).decimal_point();
+  // This is an ok convention as long as we stay in C or english locale
+  const char dot = std::use_facet<std::numpunct<char>>(std::locale::classic()).decimal_point();
   if (oss.str().find(dot) == std::string::npos)
   {
       oss << dot;