GDALDataset::ComputeInterBandCovarianceMatrix(): optimize for a all-invalid block

Author: rouault

autotest/gcore/gdal_stats.py

@@ -1501,6 +1501,35 @@ def test_stats_ComputeInterBandCovarianceMatrix_nan_result():
 ###############################################################################
 
 
+def test_stats_ComputeInterBandCovarianceMatrix_whole_block_nodata():
+
+    ds = gdal.GetDriverByName("MEM").Create("", 2, 3, 2, gdal.GDT_Float32)
+    ds.GetRasterBand(1).WriteRaster(
+        0, 0, 2, 1, struct.pack("f" * 2, math.nan, math.nan)
+    )
+    ds.GetRasterBand(1).WriteRaster(0, 1, 2, 1, struct.pack("f" * 2, 0, 1))
+    ds.GetRasterBand(1).WriteRaster(
+        0, 2, 2, 1, struct.pack("f" * 2, math.nan, math.nan)
+    )
+    ds.GetRasterBand(2).WriteRaster(
+        0, 0, 2, 1, struct.pack("f" * 2, math.nan, math.nan)
+    )
+    ds.GetRasterBand(2).WriteRaster(0, 1, 2, 1, struct.pack("f" * 2, 1, 0))
+    ds.GetRasterBand(2).WriteRaster(
+        0, 2, 2, 1, struct.pack("f" * 2, math.nan, math.nan)
+    )
+
+    expected_cov_matrix = [[0.5, -0.5], [-0.5, 0.5]]
+
+    cov_matrix = ds.ComputeInterBandCovarianceMatrix()
+    assert list(chain.from_iterable(cov_matrix)) == pytest.approx(
+        list(chain.from_iterable(expected_cov_matrix))
+    )
+
+
+###############################################################################
+
+
 def test_stats_ComputeInterBandCovarianceMatrix_failed_to_compute_stats():
 
     ds = gdal.GetDriverByName("MEM").Create("", 2, 1, 1, gdal.GDT_Float32)

gcore/gdaldataset.cpp

@@ -12953,6 +12953,7 @@ ComputeInterBandCovarianceMatrixInternal(GDALDataset *poDS,
             return eErr;
 
         // Compute the mean of all bands for this block
+        bool bAllBandsAreAllNodata = false;
         bool bNoBandHasNodata = false;
         for (int i = 0; i < nBandCount; ++i)
         {
@@ -12975,19 +12976,24 @@ ComputeInterBandCovarianceMatrixInternal(GDALDataset *poDS,
             }
             adfCurBlockMean[i] = nCount > 0 ? dfSum / nCount : ZERO;
             anCurBlockCount[i] = nCount;
+            bAllBandsAreAllNodata =
+                (i == 0 || bAllBandsAreAllNodata) && (nCount == 0);
             bNoBandHasNodata = (i == 0 || bNoBandHasNodata) &&
                                (nCount == nThisBlockPixelCount);
         }
 
         // Modify the pixel values to shift them by minus the mean
-        for (int i = 0; i < nBandCount; ++i)
+        if (!bAllBandsAreAllNodata)
         {
-            T *padfI =
-                adfCurBlockPixelsAllBands.data() + i * nThisBlockPixelCount;
-            const T dfMeanI = adfCurBlockMean[i];
-            for (size_t iPixel = 0; iPixel < nThisBlockPixelCount; ++iPixel)
+            for (int i = 0; i < nBandCount; ++i)
             {
-                padfI[iPixel] -= dfMeanI;
+                T *padfI =
+                    adfCurBlockPixelsAllBands.data() + i * nThisBlockPixelCount;
+                const T dfMeanI = adfCurBlockMean[i];
+                for (size_t iPixel = 0; iPixel < nThisBlockPixelCount; ++iPixel)
+                {
+                    padfI[iPixel] -= dfMeanI;
+                }
             }
         }
 
@@ -13057,7 +13063,20 @@ ComputeInterBandCovarianceMatrixInternal(GDALDataset *poDS,
             }
         };
 
-        if (bNoBandHasNodata)
+        if (bAllBandsAreAllNodata)
+        {
+            // Optimized code path where all values in the current block
+            // are invalid
+
+            for (int i = 0; i < nBandCount; ++i)
+            {
+                for (int j = i; j < nBandCount; ++j)
+                {
+                    UpdateGlobalValues(i, j, 0, ZERO);
+                }
+            }
+        }
+        else if (bNoBandHasNodata)
         {
             // Optimized code path where there are no invalid value in the
             // current block