VRT: add norm_diff pixel function (#8081)

OSGeo · Jul 12, 2023 · a6d03b2 · a6d03b2
1 parent 5ad5939
commit a6d03b2
Show file tree

Hide file tree

Showing 4 changed files with 93 additions and 3 deletions.
diff --git a/autotest/gcore/data/vrt/pixfun_norm_diff_r.vrt b/autotest/gcore/data/vrt/pixfun_norm_diff_r.vrt
@@ -0,0 +1,18 @@
+<VRTDataset rasterXSize="5" rasterYSize="6">
+  <VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand">
+    <Description>Normalized Difference</Description>
+    <PixelFunctionType>norm_diff</PixelFunctionType>
+    <SimpleSource>
+      <SourceFilename relativeToVRT="1">../int32.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SrcRect xOff="0" yOff="0" xSize="5" ySize="6"/>
+      <DstRect xOff="0" yOff="0" xSize="5" ySize="6"/>
+    </SimpleSource>
+    <SimpleSource>
+      <SourceFilename relativeToVRT="1">../float32.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SrcRect xOff="10" yOff="10" xSize="5" ySize="6"/>
+      <DstRect xOff="0" yOff="0" xSize="5" ySize="6"/>
+    </SimpleSource>
+  </VRTRasterBand>
+</VRTDataset>
diff --git a/autotest/gcore/pixfun.py b/autotest/gcore/pixfun.py
@@ -1045,6 +1045,29 @@ def test_pixfun_pow():
     assert numpy.allclose(data, refdata**3.14)
 
 
+###############################################################################
+# Verify the normalized difference of 2 (real) datasets.
+
+
+def test_pixfun_norm_diff_r():
+
+    filename = "data/vrt/pixfun_norm_diff_r.vrt"
+    ds = gdal.OpenShared(filename, gdal.GA_ReadOnly)
+    assert ds is not None, 'Unable to open "%s" dataset.' % filename
+    data = ds.GetRasterBand(1).ReadAsArray()
+
+    reffilename = "data/int32.tif"
+    refds = gdal.Open(reffilename)
+    assert refds is not None, 'Unable to open "%s" dataset.' % reffilename
+    refdata1 = refds.GetRasterBand(1).ReadAsArray(0, 0, 5, 6)
+
+    reffilename = "data/float32.tif"
+    refds = gdal.Open(reffilename)
+    assert refds is not None, 'Unable to open "%s" dataset.' % reffilename
+    refdata2 = refds.GetRasterBand(1).ReadAsArray(10, 10, 5, 6)
+    assert numpy.allclose(data, (refdata1 - refdata2) / (refdata1 + refdata2))
+
+
 ###############################################################################
 # Verify linear pixel interpolation
 

diff --git a/doc/source/drivers/raster/vrt.rst b/doc/source/drivers/raster/vrt.rst
@@ -890,6 +890,10 @@ GDAL provides a set of default pixel functions that can be used without writing
      - 2
      - -
      - divide one raster band by another (``b1 / b2``)
+   * - **norm_diff**
+     - 2
+     - -
+     - computes the normalized difference between two raster bands: ``(b1 - b2)/(b1 + b2)``
    * - **exp**
      - 1
      - ``base`` (optional), ``fact`` (optional)

diff --git a/frmts/vrt/pixelfunctions.cpp b/frmts/vrt/pixelfunctions.cpp
@@ -1471,6 +1471,49 @@ static CPLErr ScalePixelFunc(void **papoSources, int nSources, void *pData,
     return CE_None;
 }
 
+static CPLErr NormDiffPixelFunc(void **papoSources, int nSources, void *pData,
+                                int nXSize, int nYSize, GDALDataType eSrcType,
+                                GDALDataType eBufType, int nPixelSpace,
+                                int nLineSpace)
+{
+    /* ---- Init ---- */
+    if (nSources != 2)
+        return CE_Failure;
+
+    if (GDALDataTypeIsComplex(eSrcType))
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "norm_diff cannot by applied to complex data types");
+        return CE_Failure;
+    }
+
+    /* ---- Set pixels ---- */
+    size_t ii = 0;
+    for (int iLine = 0; iLine < nYSize; ++iLine)
+    {
+        for (int iCol = 0; iCol < nXSize; ++iCol, ++ii)
+        {
+            const double dfLeftVal = GetSrcVal(papoSources[0], eSrcType, ii);
+            const double dfRightVal = GetSrcVal(papoSources[1], eSrcType, ii);
+
+            const double dfDenom = (dfLeftVal + dfRightVal);
+
+            const double dfPixVal =
+                dfDenom == 0 ? std::numeric_limits<double>::infinity()
+                             : (dfLeftVal - dfRightVal) / dfDenom;
+
+            GDALCopyWords(&dfPixVal, GDT_Float64, 0,
+                          static_cast<GByte *>(pData) +
+                              static_cast<GSpacing>(nLineSpace) * iLine +
+                              iCol * nPixelSpace,
+                          eBufType, nPixelSpace, 1);
+        }
+    }
+
+    /* ---- Return success ---- */
+    return CE_None;
+}  // NormDiffPixelFunc
+
 /************************************************************************/
 /*                     GDALRegisterDefaultPixelFunc()                   */
 /************************************************************************/
@@ -1495,10 +1538,11 @@ static CPLErr ScalePixelFunc(void **papoSources, int nSources, void *pData,
  * - "sum": sum 2 or more raster bands
  * - "diff": computes the difference between 2 raster bands (b1 - b2)
  * - "mul": multiply 2 or more raster bands
- * - "div": divide one rasted band by another (b1 / b2).
- *          Note: no check is performed on zero division
+ * - "div": divide one raster band by another (b1 / b2).
+ * - "norm_diff": computes the normalized difference between two raster bands:
+ *                ``(b1 - b2)/(b1 + b2)``.
  * - "cmul": multiply the first band for the complex conjugate of the second
- * - "inv": inverse (1./x). Note: no check is performed on zero division
+ * - "inv": inverse (1./x).
  * - "intensity": computes the intensity Re(x*conj(x)) of a single raster band
  *                (real or complex)
  * - "sqrt": perform the square root of a single raster band (real only)
@@ -1584,6 +1628,7 @@ CPLErr GDALRegisterDefaultPixelFunc()
                                         pszReplaceNoDataPixelFuncMetadata);
     GDALAddDerivedBandPixelFuncWithArgs("scale", ScalePixelFunc,
                                         pszScalePixelFuncMetadata);
+    GDALAddDerivedBandPixelFunc("norm_diff", NormDiffPixelFunc);
 
     return CE_None;
 }