diff --git a/app/otbHarmonizer.cxx b/app/otbHarmonizer.cxx index 8bfc7a34f3e32898b62dc94a99cd4f363fe121f4..1bb1d3ad308ab736c395ff323ecf939bc62f61be 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);