ENH: add option to enable zero-y intercept

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);