From 08aa77b0ecd87640aa71509d88da1ef0ee152141 Mon Sep 17 00:00:00 2001
From: remi <remi.cresson@irstea.fr>
Date: Mon, 30 Sep 2019 16:04:33 +0200
Subject: [PATCH] ENH: add option to enable zero-y intercept

---
 app/otbHarmonizer.cxx | 26 ++++++++++++++++++++++++--
 1 file changed, 24 insertions(+), 2 deletions(-)

diff --git a/app/otbHarmonizer.cxx b/app/otbHarmonizer.cxx
index 8bfc7a3..1bb1d3a 100644
--- a/app/otbHarmonizer.cxx
+++ b/app/otbHarmonizer.cxx
@@ -296,6 +296,10 @@ private:
     AddParameter(ParameterType_Int, "smoothingradius", "Smoothing radius (in pixels of the reference image)");
     SetDefaultParameterInt         ("smoothingradius", 10);
 
+    AddParameter(ParameterType_Choice, "useoffset", "Zero-y intercept mode in linear transform");
+    AddChoice("useoffset.on", "Yes (y=ax+b)");
+    AddChoice("useoffset.off", "No (y=ax)");
+
     AddRAMParameter();
 
     // Doc example parameter settings
@@ -406,6 +410,15 @@ private:
     offset.SetSize(refImage->GetNumberOfComponentsPerPixel());
     std::stringstream exp;
     exp << "{";
+
+    if (GetParameterInt("useoffset") == 0)
+      {
+      otbAppLogINFO("Using zero-y intercept (y=ax)");
+      }
+    else
+      {
+      otbAppLogINFO("Using linear correction model (y=ax+b)");
+      }
     for (unsigned int band = 0 ; band < refImage->GetNumberOfComponentsPerPixel() ; band++)
     {
       float meanY = m_StatsFilter->GetMeans().at(band)[0][0];
@@ -414,8 +427,17 @@ private:
       float meanX2 = m_StatsFilter->GetMeansOfProducts().at(band)[1][1];
       float mean2X = meanX * meanX;
 
-      float b1 = (meanX * meanY - meanXY) / (mean2X - meanX2);
-      float b0 = meanY - b1 * meanX;
+      float b0, b1;
+      if (GetParameterInt("useoffset") == 1)
+        {
+        b1 = (meanX * meanY - meanXY) / (mean2X - meanX2);
+        b0 = meanY - b1 * meanX;
+        }
+      else
+        {
+        b1 = meanXY / meanX2;
+        b0 = 0;
+        }
 
       otbAppLogINFO("Band " << band << " gain: " << b1 << " bias: " << b0);
 
-- 
GitLab