diff --git a/autotest/gdrivers/vrtprocesseddataset.py b/autotest/gdrivers/vrtprocesseddataset.py new file mode 100755 index 000000000000..7ebec7613d0d --- /dev/null +++ b/autotest/gdrivers/vrtprocesseddataset.py @@ -0,0 +1,1258 @@ +#!/usr/bin/env pytest +############################################################################### +# $Id$ +# +# Project: GDAL/OGR Test Suite +# Purpose: Test VRTProcessedDataset support. +# Author: Even Rouault +# +############################################################################### +# Copyright (c) 2024, Even Rouault +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the "Software"), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included +# in all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. +############################################################################### + +import struct + +import gdaltest +import pytest + +from osgeo import gdal + +############################################################################### +# Test error cases in general VRTProcessedDataset XML structure + + +def test_vrtprocesseddataset_errors(tmp_vsimem): + + with pytest.raises(Exception, match="Input element missing"): + gdal.Open( + """ + + """ + ) + + with pytest.raises( + Exception, + match="Input element should have a SourceFilename or VRTDataset element", + ): + gdal.Open( + """ + + + """ + ) + + with pytest.raises(Exception): # "No such file or directory'", but O/S dependent + gdal.Open( + """ + + + """ + ) + + with pytest.raises( + Exception, + match="Missing one of rasterXSize, rasterYSize or bands on VRTDataset", + ): + gdal.Open( + """ + + + """ + ) + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 10, 5, 3) + src_ds.GetRasterBand(1).Fill(1) + src_ds.GetRasterBand(2).Fill(2) + src_ds.GetRasterBand(3).Fill(3) + src_ds.Close() + + with pytest.raises(Exception, match="ProcessingSteps element missing"): + gdal.Open( + f""" + + {src_filename} + + + """ + ) + + with pytest.raises( + Exception, match="Inconsistent declared VRT dimensions with input dataset" + ): + gdal.Open( + f""" + + {src_filename} + + + """ + ) + + with pytest.raises( + Exception, match="Inconsistent declared VRT dimensions with input dataset" + ): + gdal.Open( + f""" + + {src_filename} + + + """ + ) + + with pytest.raises(Exception, match="At least one step should be defined"): + gdal.Open( + f""" + + {src_filename} + + + + """ + ) + + +############################################################################### +# Test nominal cases of BandAffineCombination algorithm + + +def test_vrtprocesseddataset_affine_combination_nominal(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 2, 1, 3) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1, b"\x01\x03") + src_ds.GetRasterBand(2).WriteRaster(0, 0, 2, 1, b"\x02\x06") + src_ds.GetRasterBand(3).WriteRaster(0, 0, 2, 1, b"\x03\x03") + src_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + BandAffineCombination + 10,0,1,0 + 20,0,0,1 + 30,1,0,0 + 15 + 32 + + + + """ + ) + assert ds.RasterXSize == 2 + assert ds.RasterYSize == 1 + assert ds.RasterCount == 3 + assert ds.GetSpatialRef() is None + assert ds.GetGeoTransform(can_return_null=True) is None + assert ds.GetRasterBand(1).DataType == gdal.GDT_Byte + assert struct.unpack("B" * 2, ds.GetRasterBand(1).ReadRaster()) == (15, 10 + 6) + assert struct.unpack("B" * 2, ds.GetRasterBand(2).ReadRaster()) == (20 + 3, 20 + 3) + assert struct.unpack("B" * 2, ds.GetRasterBand(3).ReadRaster()) == (30 + 1, 32) + + +############################################################################### +# Test several steps in a VRTProcessedDataset + + +def test_vrtprocesseddataset_several_steps(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 10, 5, 3) + src_ds.GetRasterBand(1).Fill(1) + src_ds.GetRasterBand(2).Fill(2) + src_ds.GetRasterBand(3).Fill(3) + src_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + BandAffineCombination + 0,0,1,0 + 0,0,0,1 + 0,1,0,0 + + + BandAffineCombination + 0,0,1,0 + 0,0,0,1 + 0,1,0,0 + + + BandAffineCombination + 0,0,1,0 + 0,0,0,1 + 0,1,0,0 + + + + """ + ) + assert ds.RasterXSize == 10 + assert ds.RasterYSize == 5 + assert ds.RasterCount == 3 + assert ds.GetSpatialRef() is None + assert ds.GetGeoTransform(can_return_null=True) is None + assert ds.GetRasterBand(1).DataType == gdal.GDT_Byte + assert ds.GetRasterBand(1).ComputeRasterMinMax(False) == (1, 1) + assert ds.GetRasterBand(2).ComputeRasterMinMax(False) == (2, 2) + assert ds.GetRasterBand(3).ComputeRasterMinMax(False) == (3, 3) + + +############################################################################### +# Test nominal cases of BandAffineCombination algorithm with nodata + + +def test_vrtprocesseddataset_affine_combination_nodata(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 2, 1, 2) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1, b"\x01\x02") + src_ds.GetRasterBand(1).SetNoDataValue(1) + src_ds.GetRasterBand(2).WriteRaster(0, 0, 2, 1, b"\x03\x03") + src_ds.GetRasterBand(2).SetNoDataValue(1) + src_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + BandAffineCombination + 0,1,1 + 0,1,-1 + + + + """ + ) + assert ds.GetRasterBand(1).DataType == gdal.GDT_Byte + assert struct.unpack("B" * 2, ds.GetRasterBand(1).ReadRaster()) == (1, 5) + # 0 should actually be 3-2=1, but this is the nodata value hence the replacement value + assert struct.unpack("B" * 2, ds.GetRasterBand(2).ReadRaster()) == (1, 0) + + +def test_vrtprocesseddataset_affine_combination_nodata_as_parameter(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 2, 1, 2) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1, b"\x01\x02") + src_ds.GetRasterBand(2).WriteRaster(0, 0, 2, 1, b"\x03\x03") + src_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + BandAffineCombination + 0,1,1 + 256,1,-1 + 1 + 255 + Byte + + + + """ + ) + assert ds.GetRasterBand(1).DataType == gdal.GDT_Byte + assert struct.unpack("B" * 2, ds.GetRasterBand(1).ReadRaster()) == (255, 5) + # 254 should actually be 256+1*2+(-1)*3=255, but this is the nodata value hence the replacement value + assert struct.unpack("B" * 2, ds.GetRasterBand(2).ReadRaster()) == (255, 254) + + +############################################################################### +# Test replacement_nodata logic of BandAffineCombination + + +def test_vrtprocesseddataset_affine_combination_replacement_nodata(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 2, 1, 2) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1, b"\x01\x02") + src_ds.GetRasterBand(2).WriteRaster(0, 0, 2, 1, b"\x03\x03") + src_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + BandAffineCombination + 0,1,1 + 256,1,-1 + 1 + 255 + 128 + + + + """ + ) + assert ds.GetRasterBand(1).DataType == gdal.GDT_Byte + assert struct.unpack("B" * 2, ds.GetRasterBand(1).ReadRaster()) == (255, 5) + # 254 should actually be 256+1*2+(-1)*3=255, but this is the nodata value hence the replacement value + assert struct.unpack("B" * 2, ds.GetRasterBand(2).ReadRaster()) == (255, 128) + + +############################################################################### +# Test error cases of BandAffineCombination algorithm + + +def test_vrtprocesseddataset_affine_combination_errors(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 10, 5, 3) + src_ds.GetRasterBand(1).Fill(1) + src_ds.GetRasterBand(2).Fill(2) + src_ds.GetRasterBand(3).Fill(3) + src_ds.Close() + + with pytest.raises( + Exception, + match="Step 'Affine combination of band values' lacks required Argument 'coefficients_{band}'", + ): + gdal.Open( + f""" + + {src_filename} + + + + BandAffineCombination + + + + """ + ) + + with pytest.raises( + Exception, match="Argument coefficients_1 has 3 values, whereas 4 are expected" + ): + gdal.Open( + f""" + + {src_filename} + + + + BandAffineCombination + 10,0,1 + + + + """ + ) + + with pytest.raises(Exception, match="Argument coefficients_3 is missing"): + gdal.Open( + f""" + + {src_filename} + + + + BandAffineCombination + 10,0,1,0 + 10,0,1,0 + 10,0,1,0 + + + + """ + ) + + with pytest.raises( + Exception, + match="Final step expect 3 bands, but only 1 coefficient_XX are provided", + ): + gdal.Open( + f""" + + {src_filename} + + + + BandAffineCombination + 10,0,1,0 + + + + """ + ) + + +############################################################################### +# Test nominal cases of LUT algorithm + + +def test_vrtprocesseddataset_lut_nominal(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 3, 1, 2) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 3, 1, b"\x01\x02\x03") + src_ds.GetRasterBand(2).WriteRaster(0, 0, 3, 1, b"\x01\x02\x03") + src_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + LUT + 1.5:10,2.5:20 + 1.5:100,2.5:200 + + + + """ + ) + assert struct.unpack("B" * 3, ds.GetRasterBand(1).ReadRaster()) == (10, 15, 20) + assert struct.unpack("B" * 3, ds.GetRasterBand(2).ReadRaster()) == (100, 150, 200) + + +############################################################################### +# Test nominal cases of LUT algorithm with nodata coming from input dataset + + +def test_vrtprocesseddataset_lut_nodata(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 4, 1, 2) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 4, 1, b"\x00\x01\x02\x03") + src_ds.GetRasterBand(1).SetNoDataValue(0) + src_ds.GetRasterBand(2).WriteRaster(0, 0, 4, 1, b"\x00\x01\x02\x03") + src_ds.GetRasterBand(2).SetNoDataValue(0) + src_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + LUT + 1.5:10,2.5:20 + 1.5:100,2.5:200 + + + + """ + ) + assert struct.unpack("B" * 4, ds.GetRasterBand(1).ReadRaster()) == (0, 10, 15, 20) + assert struct.unpack("B" * 4, ds.GetRasterBand(2).ReadRaster()) == ( + 0, + 100, + 150, + 200, + ) + + +############################################################################### +# Test nominal cases of LUT algorithm with nodata set as a parameter + + +def test_vrtprocesseddataset_lut_nodata_as_parameter(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 4, 1, 2) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 4, 1, b"\x00\x01\x02\x03") + src_ds.GetRasterBand(2).WriteRaster(0, 0, 4, 1, b"\x00\x01\x02\x03") + src_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + LUT + 1.5:10,2.5:20 + 1.5:100,2.5:200 + 0 + 1 + + + + """ + ) + assert struct.unpack("B" * 4, ds.GetRasterBand(1).ReadRaster()) == (1, 10, 15, 20) + assert struct.unpack("B" * 4, ds.GetRasterBand(2).ReadRaster()) == ( + 1, + 100, + 150, + 200, + ) + + +############################################################################### +# Test error cases of LUT algorithm + + +def test_vrtprocesseddataset_lut_errors(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 3, 1, 2) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 3, 1, b"\x01\x02\x03") + src_ds.GetRasterBand(2).WriteRaster(0, 0, 3, 1, b"\x01\x02\x03") + src_ds.Close() + + with pytest.raises(Exception, match="Step 'nr 1' lacks required Argument"): + gdal.Open( + f""" + + {src_filename} + + + + LUT + + + + """ + ) + + with pytest.raises(Exception, match="Invalid value for argument 'lut_1'"): + gdal.Open( + f""" + + {src_filename} + + + + LUT + 1.5:10,2.5 + + + + """ + ) + + with pytest.raises(Exception, match="Invalid band in argument 'lut_3'"): + gdal.Open( + f""" + + {src_filename} + + + + LUT + 1.5:10,2.5:20 + 1.5:10,2.5:20 + + + + """ + ) + + with pytest.raises(Exception, match="Missing lut_XX element"): + gdal.Open( + f""" + + {src_filename} + + + + LUT + 1.5:10,2.5:20 + + + + """ + ) + + +############################################################################### +# Test nominal case of Dehazing algorithm + + +def test_vrtprocesseddataset_dehazing_nominal(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 6, 1, 2) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 6, 1, b"\x01\x02\x03\xff\x01\x01") + src_ds.GetRasterBand(2).WriteRaster(0, 0, 6, 1, b"\x01\x02\x03\xff\x01\x01") + src_ds.GetRasterBand(1).SetNoDataValue(255) + src_ds.GetRasterBand(2).SetNoDataValue(255) + src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + src_ds.Close() + + gain_filename = str(tmp_vsimem / "gain.tif") + gain_ds = gdal.GetDriverByName("GTiff").Create(gain_filename, 6, 1, 2) + gain_ds.GetRasterBand(1).WriteRaster(0, 0, 6, 1, b"\x02\x04\x06\x01\xfe\x01") + gain_ds.GetRasterBand(2).WriteRaster(0, 0, 6, 1, b"\x03\x05\x07\x01\xfe\x01") + gain_ds.GetRasterBand(1).SetNoDataValue(254) + gain_ds.GetRasterBand(2).SetNoDataValue(254) + gain_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + gain_ds.Close() + + offset_filename = str(tmp_vsimem / "offset.tif") + offset_ds = gdal.GetDriverByName("GTiff").Create(offset_filename, 6, 1, 2) + offset_ds.GetRasterBand(1).WriteRaster(0, 0, 6, 1, b"\x01\x02\x03\x01\x01\xfd") + offset_ds.GetRasterBand(2).WriteRaster(0, 0, 6, 1, b"\x02\x03\x04\x01\x01\xfd") + offset_ds.GetRasterBand(1).SetNoDataValue(253) + offset_ds.GetRasterBand(2).SetNoDataValue(253) + offset_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + offset_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + Dehazing + {gain_filename} + 1 + {gain_filename} + 2 + {offset_filename} + 1 + {offset_filename} + 2 + 2 + 16 + + + + """ + ) + assert struct.unpack("B" * 6, ds.GetRasterBand(1).ReadRaster()) == ( + 2, + 6, + 15, + 255, + 255, + 255, + ) + assert struct.unpack("B" * 6, ds.GetRasterBand(2).ReadRaster()) == ( + 2, + 7, + 16, + 255, + 255, + 255, + ) + + +############################################################################### +# Test nominal case of Dehazing algorithm where gain and offset have a lower +# resolution than the input dataset + + +def test_vrtprocesseddataset_dehazing_different_resolution(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 6, 2, 1) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 6, 2, b"\x01\x01\x02\x02\x03\x03" * 2) + src_ds.SetGeoTransform([0, 0.5, 0, 0, 0, 0.5]) + src_ds.Close() + + gain_filename = str(tmp_vsimem / "gain.tif") + gain_ds = gdal.GetDriverByName("GTiff").Create(gain_filename, 3, 1, 1) + gain_ds.GetRasterBand(1).WriteRaster(0, 0, 3, 1, b"\x02\x04\x06") + gain_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + gain_ds.Close() + + offset_filename = str(tmp_vsimem / "offset.tif") + offset_ds = gdal.GetDriverByName("GTiff").Create(offset_filename, 3, 1, 1) + offset_ds.GetRasterBand(1).WriteRaster(0, 0, 3, 1, b"\x01\x02\x03") + offset_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + offset_ds.Close() + + ds = gdal.Open( + f""" + + {src_filename} + + + + Dehazing + {gain_filename} + 1 + {offset_filename} + 1 + + + + """ + ) + assert struct.unpack("B" * 12, ds.GetRasterBand(1).ReadRaster()) == ( + 1, + 2, + 6, + 8, + 15, + 15, + 1, + 2, + 6, + 8, + 15, + 15, + ) + + +############################################################################### +# Test error cases of Dehazing algorithm + + +def test_vrtprocesseddataset_dehazing_error(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 3, 1, 1) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 3, 1, b"\x01\x02\x03") + src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + src_ds.Close() + + with pytest.raises( + Exception, + match="Step 'nr 1' lacks required Argument 'offset_dataset_band_{band}'", + ): + gdal.Open( + f""" + + {src_filename} + + + + Dehazing + {src_filename} + 1 + + + + """ + ) + + with pytest.raises( + Exception, + match="Invalid band in argument 'gain_dataset_filename_2'", + ): + gdal.Open( + f""" + + {src_filename} + + + + Dehazing + {src_filename} + 1 + {src_filename} + 1 + + + + """ + ) + + with pytest.raises( + Exception, + match="Invalid band in argument 'gain_dataset_band_2'", + ): + gdal.Open( + f""" + + {src_filename} + + + + Dehazing + {src_filename} + 1 + {src_filename} + 1 + + + + """ + ) + + with pytest.raises( + Exception, + match="Invalid band in argument 'offset_dataset_filename_2'", + ): + gdal.Open( + f""" + + {src_filename} + + + + Dehazing + {src_filename} + 1 + {src_filename} + 1 + + + + """ + ) + + with pytest.raises( + Exception, + match="Invalid band in argument 'offset_dataset_band_2'", + ): + gdal.Open( + f""" + + {src_filename} + + + + Dehazing + {src_filename} + 1 + {src_filename} + 1 + + + + """ + ) + + with pytest.raises( + Exception, + match=r"Invalid band number \(2\) for a gain dataset", + ): + gdal.Open( + f""" + + {src_filename} + + + + Dehazing + {src_filename} + 2 + {src_filename} + 1 + + + + """ + ) + + with pytest.raises(Exception): # "No such file or directory'", but O/S dependent + gdal.Open( + f""" + + {src_filename} + + + + Dehazing + invalid + 1 + {src_filename} + 1 + + + + """ + ) + + nogt_filename = str(tmp_vsimem / "nogt.tif") + ds = gdal.GetDriverByName("GTiff").Create(nogt_filename, 1, 1, 1) + ds.Close() + + with pytest.raises(Exception, match="lacks a geotransform"): + gdal.Open( + f""" + + {src_filename} + + + + Dehazing + {nogt_filename} + 1 + {nogt_filename} + 1 + + + + """ + ) + + +############################################################################### +# Test nominal cases of Trimming algorithm + + +def test_vrtprocesseddataset_trimming_nominal(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 6, 1, 4) + + R = 100.0 + G = 150.0 + B = 200.0 + NIR = 100.0 + + src_ds.GetRasterBand(1).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, int(R), 150, 200, 0, 0, 0) + ) + src_ds.GetRasterBand(2).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, int(G), 200, 100, 0, 0, 0) + ) + src_ds.GetRasterBand(3).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, int(B), 100, 150, 0, 0, 0) + ) + src_ds.GetRasterBand(4).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, int(NIR), 150, 200, 0, 0, 0) + ) + src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + src_ds.Close() + + trimming_filename = str(tmp_vsimem / "trimming.tif") + trimming_ds = gdal.GetDriverByName("GTiff").Create(trimming_filename, 6, 1, 1) + + localMaxRGB = 205.0 + + trimming_ds.GetRasterBand(1).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, int(localMaxRGB), 210, 220, 0, 0, 0) + ) + trimming_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + trimming_ds.Close() + + top_rgb = 200.0 + tone_ceil = 190.0 + top_margin = 0.1 + + ds = gdal.Open( + f""" + + {src_filename} + + + + Trimming + {trimming_filename} + {top_rgb} + {tone_ceil} + {top_margin} + + + + """ + ) + + # Do algorithm at hand + + # Extract local saturation value from trimming image + reducedRGB = min((1.0 - top_margin) * top_rgb / localMaxRGB, 1) + + # RGB bands specific process + maxRGB = max(R, G, B) + toneMaxRGB = min(tone_ceil / maxRGB, 1) + toneR = min(tone_ceil / R, 1) + toneG = min(tone_ceil / G, 1) + toneB = min(tone_ceil / B, 1) + outputR = min(reducedRGB * R * toneR / toneMaxRGB, top_rgb) + outputG = min(reducedRGB * G * toneG / toneMaxRGB, top_rgb) + outputB = min(reducedRGB * B * toneB / toneMaxRGB, top_rgb) + + # Other bands processing (NIR, ...): only apply RGB reduction factor + outputNIR = reducedRGB * NIR + + # print(outputR, outputG, outputB, outputNIR) + + assert ( + round(outputR) + == struct.unpack("B", ds.GetRasterBand(1).ReadRaster(0, 0, 1, 1))[0] + ) + assert ( + round(outputG) + == struct.unpack("B", ds.GetRasterBand(2).ReadRaster(0, 0, 1, 1))[0] + ) + assert ( + round(outputB) + == struct.unpack("B", ds.GetRasterBand(3).ReadRaster(0, 0, 1, 1))[0] + ) + assert ( + round(outputNIR) + == struct.unpack("B", ds.GetRasterBand(4).ReadRaster(0, 0, 1, 1))[0] + ) + + assert struct.unpack("B" * 6, ds.GetRasterBand(1).ReadRaster()) == ( + 92, # round(outputR) + 135, + 164, + 0, + 0, + 0, + ) + assert struct.unpack("B" * 6, ds.GetRasterBand(2).ReadRaster()) == ( + 139, # round(outputG) + 171, + 86, + 0, + 0, + 0, + ) + assert struct.unpack("B" * 6, ds.GetRasterBand(3).ReadRaster()) == ( + 176, # round(outputB) + 90, + 129, + 0, + 0, + 0, + ) + assert struct.unpack("B" * 6, ds.GetRasterBand(4).ReadRaster()) == ( + 88, # round(outputNIR) + 129, + 164, + 0, + 0, + 0, + ) + + +############################################################################### +# Test error cases of Trimming algorithm + + +def test_vrtprocesseddataset_trimming_errors(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 6, 1, 4) + src_ds.GetRasterBand(1).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, 100, 150, 200, 0, 0, 0) + ) + src_ds.GetRasterBand(2).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, 150, 200, 100, 0, 0, 0) + ) + src_ds.GetRasterBand(3).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, 200, 100, 150, 0, 0, 0) + ) + src_ds.GetRasterBand(4).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, 100, 150, 200, 0, 0, 0) + ) + src_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + src_ds.Close() + + trimming_filename = str(tmp_vsimem / "trimming.tif") + trimming_ds = gdal.GetDriverByName("GTiff").Create(trimming_filename, 6, 1, 1) + trimming_ds.GetRasterBand(1).WriteRaster( + 0, 0, 6, 1, struct.pack("B" * 6, 200, 210, 220, 0, 0, 0) + ) + trimming_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + trimming_ds.Close() + + trimming_two_bands_filename = str(tmp_vsimem / "trimming_two_bands.tif") + trimming_ds = gdal.GetDriverByName("GTiff").Create( + trimming_two_bands_filename, 6, 1, 2 + ) + trimming_ds.SetGeoTransform([0, 1, 0, 0, 0, 1]) + trimming_ds.Close() + + with pytest.raises(Exception): + gdal.Open( + f""" + + {src_filename} + + + + Trimming + invalid + 200 + 190 + 0.1 + + + + """ + ) + + for val in (0, 5): + with pytest.raises(Exception, match="Invalid band in argument 'red_band'"): + gdal.Open( + f""" + + {src_filename} + + + + Trimming + {trimming_filename} + {val} + 200 + 190 + 0.1 + + + + """ + ) + + for val in (0, 5): + with pytest.raises(Exception, match="Invalid band in argument 'green_band'"): + gdal.Open( + f""" + + {src_filename} + + + + Trimming + {trimming_filename} + {val} + 200 + 190 + 0.1 + + + + """ + ) + + for val in (0, 5): + with pytest.raises(Exception, match="Invalid band in argument 'blue_band'"): + gdal.Open( + f""" + + {src_filename} + + + + Trimming + {trimming_filename} + {val} + 200 + 190 + 0.1 + + + + """ + ) + + for (red_band, green_band, blue_band) in [(1, 1, 3), (3, 2, 3), (1, 3, 3)]: + with pytest.raises( + Exception, + match="red_band, green_band and blue_band must have distinct values", + ): + gdal.Open( + f""" + + {src_filename} + + + + Trimming + {trimming_filename} + {red_band} + {green_band} + {blue_band} + 200 + 190 + 0.1 + + + + """ + ) + + with pytest.raises(Exception, match="Trimming dataset should have a single band"): + gdal.Open( + f""" + + {src_filename} + + + + Trimming + {trimming_two_bands_filename} + 200 + 190 + 0.1 + + + + """ + ) + + +############################################################################### +# Test that serialization (for example due to statistics computation) properly +# works + + +def test_vrtprocesseddataset_serialize(tmp_vsimem): + + src_filename = str(tmp_vsimem / "src.tif") + src_ds = gdal.GetDriverByName("GTiff").Create(src_filename, 2, 1, 1) + src_ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1, b"\x01\x02") + src_ds.Close() + + vrt_filename = str(tmp_vsimem / "the.vrt") + content = f""" + + + {src_filename} + + + + BandAffineCombination + 10,1 + + + + """ + with gdaltest.tempfile(vrt_filename, content): + ds = gdal.Open(vrt_filename) + assert struct.unpack("B" * 2, ds.GetRasterBand(1).ReadRaster()) == (11, 12) + assert ds.GetRasterBand(1).GetStatistics(False, False) == [0.0, 0.0, 0.0, -1.0] + ds.GetRasterBand(1).ComputeStatistics(False) + ds.Close() + + ds = gdal.Open(vrt_filename) + assert struct.unpack("B" * 2, ds.GetRasterBand(1).ReadRaster()) == (11, 12) + assert ds.GetRasterBand(1).GetStatistics(False, False) == [ + 11.0, + 12.0, + 11.5, + 0.5, + ] diff --git a/doc/source/drivers/raster/vrt.rst b/doc/source/drivers/raster/vrt.rst index 1c860c439ef2..4aa8da0f06a9 100644 --- a/doc/source/drivers/raster/vrt.rst +++ b/doc/source/drivers/raster/vrt.rst @@ -1907,6 +1907,22 @@ See the dedicated :ref:`vrt_multidimensional` page. vrt_multidimensional +Processed dataset VRT +--------------------- + +.. versionadded:: 3.9 + +A VRT processed dataset is a specific variant of the :ref:`raster.vrt` format, +to apply chained processing steps that apply to several bands at the same time. + +See the dedicated :ref:`vrt_processed_dataset` page. + +.. toctree:: + :maxdepth: 1 + :hidden: + + vrt_processed_dataset + vrt:// connection string ------------------------ diff --git a/doc/source/drivers/raster/vrt_processed_dataset.rst b/doc/source/drivers/raster/vrt_processed_dataset.rst new file mode 100644 index 000000000000..9f00f7f13d74 --- /dev/null +++ b/doc/source/drivers/raster/vrt_processed_dataset.rst @@ -0,0 +1,261 @@ +.. _vrt_processed_dataset: + +================================================================================ +VRT processed dataset +================================================================================ + +.. versionadded:: 3.9 + +A VRT processed dataset is a specific variant of the :ref:`raster.vrt` format, +to apply chained processing steps that apply to several bands at the same time. + +The following built-in algorithms are introduced, and may typically be applied +in the following order: + +- Dehazing: remove haze effects by applying (subsampled) gain and offset + auxiliary datasets. + +- BandAffineCombination: perform an affine transformation combination of bands. + +- Trimming: apply local thresholding of saturation + +- LUT: apply a look-up table (band per band) + +More algorithms can be registered are run-time with the :cpp:func:`GDALVRTRegisterProcessedDatasetFunc` +function` + +Here's an example of such a file to apply various correction to a R,G,B,NIR dataset: + +.. code-block:: xml + + + + source.tif + + + + + Dehazing + + true + + gains.tif + gains.tif + gains.tif + gains.tif + 1 + 2 + 3 + 4 + + offsets.tif + offsets.tif + offsets.tif + offsets.tif + 1 + 2 + 3 + 4 + + 0 + 1 + 10000 + + + + BandAffineCombination + 0,1.2,-0.2,0.0,0.0 + 0,-0.03,1.03,0.0,0.0 + 0,0.0,0.0,1.0,0.0 + 0,0.0,0.0,0.0,1.0 + + 1 + 10000 + + + + Trimming + true + trimming.tif + 10000 + 0 + 10000 + + + + LUT + + 0:0,10000.0:255 + + + 0:0,10000.0:255 + + + 0:0,10000.0:255 + + + 0:0,10000.0:255 + + + + + + Red + + + Green + + + Blue + + + + + +.vrt format +----------- + +The ``VRTDataset`` root element must have a ``subClass="VRTProcessedDataset"`` attribute. + +The following child elements of ``VRTDataset`` may be defined: ``SRS``, ``GeoTransform``, ``Metadata``. If they are not explicitly set, they are inferred from the input dataset. + +``VRTRasterBand`` elements may be explicitly defined, in particular if the data type of the virtual dataset after all processing steps is different from the input one, or if the number of output bands is different from the number of input bands. If there is no explicit ``VRTRasterBand`` element, the number and data types of input bands are used implicitly. When explicitly defined, ``VRTRasterBand`` elements must have a ``subClass="VRTProcessedRasterBand"`` attribute. +` +It must also have the 2 following child elements: + +- ``Input``, which must have one and only one of the following ``SourceFilename`` or ``VRTDataset`` as child elements, to define the input dataset to which to apply the processing steps. + +- ``ProcessingSteps``, with at least one child ``Step`` element. + +Each ``Step`` must have a ``Algorithm`` child element, and an optional ``name`` attribute. +The value of ``Algorithm`` must be a registered VRTProcessedDataset function. At time of writing, the following 4 algorithms are defined: ``Dehazing``, ``BandAffineCombination``, ``Trimming`` and ``LUT``. + +A ``Step`` will generally have one or several ``Argument`` child elements, some of them being required, others optional. Consult the documentation of each algorithm. + +Dehazing algorithm +------------------ + +Remove haze effects by applying (subsampled) gain and offset auxiliary datasets. + +The gain and offset auxiliary datasets must have a georeferencing consistent of +the input dataset, but may have a different resolution. + +The formula applied by that algorithm is: ``output_value = clamp(input_value * gain - offset, min, max)`` + +The following required arguments must be specified: + +- ``gain_dataset_filename_{band}``: Filename to the gain dataset, where {band} must be replaced by 1 to the number of input bands. + +- ``gain_dataset_band_{band}``: Band number corresponding to ``gain_dataset_filename_{band}``, where {band} must be replaced by 1 to the number of input bands. + +- ``offset_dataset_filename_{band}``: Filename to the offset dataset, where {band} must be replaced by 1 to the number of input bands. + +- ``offset_dataset_band_{band}``: Band number corresponding to ``offset_dataset_filename_{band}``, where {band} must be replaced by 1 to the number of input bands. + + +The following optional arguments may be specified: + +- ``relative_filenames``: Whether gain and offset filenames are relative to the VRT. Allowed values are ``true`` and ``false``. Defaults to ``false`` + +- ``min``: Clamp minimum value, applied before writing the output value. + +- ``max``: Clamp maximum value, applied before writing the output value. + +- ``nodata``: Override the input nodata value coming from the previous step (or the input dataset for the first step). + +- ``gain_nodata``: Override the nodata value coming from the gain dataset(s). + +- ``offset_nodata``: Override the nodata value coming from the offset dataset(s). + + +BandAffineCombination algorithm +------------------------------- + +Perform an affine transformation combination of bands. + +The following required argument must be specified: + +- ``coefficients_{band}``: Comma-separated coefficients for combining bands where {band} must be replaced by 1 to the number of output bands. The number of coefficients in each argument must be 1 + number_of_input_bands, where the first coefficient is a constant, the second coefficient is the weight of the first input band, the third coefficient is the weight of the second input band, etc. + + +The following optional arguments may be specified: + +- ``src_nodata``: Override the input nodata value coming from the previous step (or the input dataset for the first step). + +- ``dst_nodata``: Set the output nodata value. + +- ``replacement_nodata``: Value to substitute to a valid computed value that would be equal to dst_nodata. + +- ``dst_intended_datatype``: Intended datatype of output (which might be different than the working data type). Used to infer an appropriate value for replacement_nodata when it is not specified. + +- ``min``: Clamp minimum value, applied before writing the output value. + +- ``max``: Clamp maximum value, applied before writing the output value. + + +Trimming algorithm +------------------ + +Apply local thresholding of saturation, with a special processing of the R,G,B bands compared to other bands. + +The pseudo algorithm used for each pixel is: + +.. code-block:: + + // Extract local saturation value from trimming image + localMaxRGB = value from TrimmingImage + reducedRGB = min ( (1-top_margin)*top_rgb/localMaxRGB ; 1) + + // RGB bands specific process + RGB[] = get red, green, blue components of input buffer + maxRGB = max(RGB[]) + toneMaxRGB = min ( toneCeil/maxRGB ; 1) + toneBand[] = min ( toneCeil/RGB[] ; 1) + + output_value_RGB[] = min ( reducedRGB*RGB[]*toneBand[] / toneMaxRGB ; topRGB) + + // Other bands processing (NIR, ...): only apply RGB reduction factor + Trimmed(OtherBands[]) = reducedRGB * OtherBands[] + + +The following required arguments must be specified: + +- ``trimming_dataset_filename``: Filename of the trimming dataset. It must have one single band. It must have a georeferencing consistent of the input dataset, but may have a different resolution. + +- ``top_rgb``: Maximum saturating RGB output value. + +- ``tone_ceil``: Maximum threshold beyond which we give up saturation. + +- ``top_margin``: Margin to allow for dynamics in brighest areas (between 0 and 1, should be close to 0) + + +The following optional arguments may be specified: + +- ``relative_filenames``: Whether the trimming dataset filename is relative to the VRT. Allowed values are ``true`` and ``false``. Defaults to ``false`` + +- ``red_band``: Index (one-based) of the red band. Defaults to 1. + +- ``green_band``: Index (one-based) of the green band. Defaults to 1. + +- ``blue_band``: Index (one-based) of the blue band. Defaults to 1. + +- ``nodata``: Override the input nodata value coming from the previous step (or the input dataset for the first step). + +- ``trimming_nodata``: Override the nodata value coming from the trimming dataset. + + +LUT +--- + +Apply a look-up table (band per band), typically to get from UInt16 to Byte data types. + +The following required argument must be specified: + +- ``lut_{band}``: List of the form ``[src value 1]:[dest value 1],[src value 2]:[dest value 2],....``. {band} must be replaced by 1 to the number of bands. + + +The following optional arguments may be specified: + +- ``src_nodata``: Override the input nodata value coming from the previous step (or the input dataset for the first step). + +- ``dst_nodata``: Set the output nodata value. diff --git a/frmts/vrt/CMakeLists.txt b/frmts/vrt/CMakeLists.txt index ba770561f569..447d48eaedbe 100644 --- a/frmts/vrt/CMakeLists.txt +++ b/frmts/vrt/CMakeLists.txt @@ -14,6 +14,8 @@ add_gdal_driver( vrtdataset.cpp pixelfunctions.cpp vrtpansharpened.cpp + vrtprocesseddataset.cpp + vrtprocesseddatasetfunctions.cpp vrtmultidim.cpp gdaltileindexdataset.cpp STRONG_CXX_WFLAGS) diff --git a/frmts/vrt/data/gdalvrt.xsd b/frmts/vrt/data/gdalvrt.xsd index 5e204df9ab98..89ce1496a4b7 100644 --- a/frmts/vrt/data/gdalvrt.xsd +++ b/frmts/vrt/data/gdalvrt.xsd @@ -30,30 +30,76 @@ ****************************************************************************/ --> - - - - - - - - - - - - - - - - - - - - - - + + + Root element + + + + + + + + + + + + May be repeated + + + + + May be repeated + + + + + + Allowed only if subClass="VRTWarpedDataset" + + + + + Allowed only if subClass="VRTPansharpenedDataset" + + + + + Allowed only if subClass="VRTProcessedDataset" + + + + + Allowed only if subClass="VRTProcessedDataset" + + + + + only for multidimensional dataset + + + + + + + + + + + + + + + + + Added in GDAL 3.9 + + + + + @@ -158,6 +204,51 @@ + + + + + + + + + + + + + + + + + + Processing step of a VRTPansharpenedDataset + + + + + Builtin allowed names are BandAffineCombination, LUT, Dehazing, Trimming. More algorithms can be registered at run-time. + + + + + + + + + + Argument of a processing function + + + + + + Allowed names are specific of each processing function + + + + + + @@ -234,6 +325,7 @@ + diff --git a/frmts/vrt/vrtdataset.cpp b/frmts/vrt/vrtdataset.cpp index 9331480d9bea..168a2638f876 100644 --- a/frmts/vrt/vrtdataset.cpp +++ b/frmts/vrt/vrtdataset.cpp @@ -97,10 +97,7 @@ VRTDataset::~VRTDataset() { VRTDataset::FlushCache(true); - if (m_poSRS) - m_poSRS->Release(); - if (m_poGCP_SRS) - m_poGCP_SRS->Release(); + if (m_nGCPCount > 0) { GDALDeinitGCPs(m_nGCPCount, m_pasGCPList); @@ -155,6 +152,17 @@ CPLErr VRTPansharpenedDataset::FlushCache(bool bAtClosing) /* FlushCache() */ /************************************************************************/ +CPLErr VRTProcessedDataset::FlushCache(bool bAtClosing) + +{ + return VRTFlushCacheStruct::FlushCache(*this, + bAtClosing); +} + +/************************************************************************/ +/* FlushCache() */ +/************************************************************************/ + template CPLErr VRTFlushCacheStruct::FlushCache(T &obj, bool bAtClosing) { @@ -318,7 +326,7 @@ CPLXMLNode *VRTDataset::SerializeToXML(const char *pszVRTPathIn) if (m_nGCPCount > 0) { GDALSerializeGCPListToXML(psDSTree, m_pasGCPList, m_nGCPCount, - m_poGCP_SRS); + m_poGCP_SRS.get()); } /* -------------------------------------------------------------------- */ @@ -410,10 +418,18 @@ CPLXMLNode *CPL_STDCALL VRTSerializeToXML(VRTDatasetH hDataset, /************************************************************************/ VRTRasterBand *VRTDataset::InitBand(const char *pszSubclass, int nBand, - bool bAllowPansharpened) + bool bAllowPansharpenedOrProcessed) { VRTRasterBand *poBand = nullptr; - if (EQUAL(pszSubclass, "VRTSourcedRasterBand")) + if (auto poProcessedDS = dynamic_cast(this)) + { + if (bAllowPansharpenedOrProcessed && + EQUAL(pszSubclass, "VRTProcessedRasterBand")) + { + poBand = new VRTProcessedRasterBand(poProcessedDS, nBand); + } + } + else if (EQUAL(pszSubclass, "VRTSourcedRasterBand")) poBand = new VRTSourcedRasterBand(this, nBand); else if (EQUAL(pszSubclass, "VRTDerivedRasterBand")) poBand = new VRTDerivedRasterBand(this, nBand); @@ -422,13 +438,17 @@ VRTRasterBand *VRTDataset::InitBand(const char *pszSubclass, int nBand, else if (EQUAL(pszSubclass, "VRTWarpedRasterBand") && dynamic_cast(this) != nullptr) poBand = new VRTWarpedRasterBand(this, nBand); - else if (bAllowPansharpened && + else if (bAllowPansharpenedOrProcessed && EQUAL(pszSubclass, "VRTPansharpenedRasterBand") && dynamic_cast(this) != nullptr) poBand = new VRTPansharpenedRasterBand(this, nBand); - else + + if (!poBand) + { CPLError(CE_Failure, CPLE_AppDefined, "VRTRasterBand of unrecognized subclass '%s'.", pszSubclass); + } + return poBand; } @@ -448,9 +468,7 @@ CPLErr VRTDataset::XMLInit(const CPLXMLNode *psTree, const char *pszVRTPathIn) const CPLXMLNode *psSRSNode = CPLGetXMLNode(psTree, "SRS"); if (psSRSNode) { - if (m_poSRS) - m_poSRS->Release(); - m_poSRS = new OGRSpatialReference(); + m_poSRS.reset(new OGRSpatialReference()); m_poSRS->SetFromUserInput( CPLGetXMLValue(psSRSNode, nullptr, ""), OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get()); @@ -505,8 +523,10 @@ CPLErr VRTDataset::XMLInit(const CPLXMLNode *psTree, const char *pszVRTPathIn) /* -------------------------------------------------------------------- */ if (const CPLXMLNode *psGCPList = CPLGetXMLNode(psTree, "GCPList")) { + OGRSpatialReference *poSRS = nullptr; GDALDeserializeGCPListFromXML(psGCPList, &m_pasGCPList, &m_nGCPCount, - &m_poGCP_SRS); + &poSRS); + m_poGCP_SRS.reset(poSRS); } /* -------------------------------------------------------------------- */ @@ -563,6 +583,13 @@ CPLErr VRTDataset::XMLInit(const CPLXMLNode *psTree, const char *pszVRTPathIn) { const char *pszSubclass = CPLGetXMLValue(psChild, "subclass", "VRTSourcedRasterBand"); + if (dynamic_cast(this) && + !EQUAL(pszSubclass, "VRTProcessedRasterBand")) + { + CPLError(CE_Failure, CPLE_NotSupported, + "Only subClass=VRTProcessedRasterBand supported"); + return CE_Failure; + } VRTRasterBand *poBand = InitBand(pszSubclass, l_nBands + 1, true); if (poBand != nullptr && @@ -642,15 +669,13 @@ CPLErr VRTDataset::SetGCPs(int nGCPCountIn, const GDAL_GCP *pasGCPListIn, const OGRSpatialReference *poGCP_SRS) { - if (m_poGCP_SRS) - m_poGCP_SRS->Release(); if (m_nGCPCount > 0) { GDALDeinitGCPs(m_nGCPCount, m_pasGCPList); CPLFree(m_pasGCPList); } - m_poGCP_SRS = poGCP_SRS ? poGCP_SRS->Clone() : nullptr; + m_poGCP_SRS.reset(poGCP_SRS ? poGCP_SRS->Clone() : nullptr); m_nGCPCount = nGCPCountIn; @@ -668,12 +693,7 @@ CPLErr VRTDataset::SetGCPs(int nGCPCountIn, const GDAL_GCP *pasGCPListIn, CPLErr VRTDataset::SetSpatialRef(const OGRSpatialReference *poSRS) { - if (m_poSRS) - m_poSRS->Release(); - if (poSRS) - m_poSRS = poSRS->Clone(); - else - m_poSRS = nullptr; + m_poSRS.reset(poSRS ? poSRS->Clone() : nullptr); SetNeedsFlush(); @@ -874,8 +894,7 @@ GDALDataset *VRTDataset::Open(GDALOpenInfo *poOpenInfo) /* -------------------------------------------------------------------- */ /* Turn the XML representation into a VRTDataset. */ /* -------------------------------------------------------------------- */ - VRTDataset *poDS = static_cast( - OpenXML(pszXML, pszVRTPath, poOpenInfo->eAccess)); + VRTDataset *poDS = OpenXML(pszXML, pszVRTPath, poOpenInfo->eAccess); if (poDS != nullptr) poDS->m_bNeedsFlush = false; @@ -1588,8 +1607,8 @@ GDALDataset *VRTDataset::OpenVRTProtocol(const char *pszSpec) /* of the dataset. */ /************************************************************************/ -GDALDataset *VRTDataset::OpenXML(const char *pszXML, const char *pszVRTPath, - GDALAccess eAccessIn) +VRTDataset *VRTDataset::OpenXML(const char *pszXML, const char *pszVRTPath, + GDALAccess eAccessIn) { /* -------------------------------------------------------------------- */ @@ -1610,8 +1629,10 @@ GDALDataset *VRTDataset::OpenXML(const char *pszXML, const char *pszVRTPath, const bool bIsPansharpened = strcmp(pszSubClass, "VRTPansharpenedDataset") == 0; + const bool bIsProcessed = strcmp(pszSubClass, "VRTProcessedDataset") == 0; - if (!bIsPansharpened && CPLGetXMLNode(psRoot, "Group") == nullptr && + if (!bIsPansharpened && !bIsProcessed && + CPLGetXMLNode(psRoot, "Group") == nullptr && (CPLGetXMLNode(psRoot, "rasterXSize") == nullptr || CPLGetXMLNode(psRoot, "rasterYSize") == nullptr || CPLGetXMLNode(psRoot, "VRTRasterBand") == nullptr)) @@ -1628,7 +1649,8 @@ GDALDataset *VRTDataset::OpenXML(const char *pszXML, const char *pszVRTPath, const int nXSize = atoi(CPLGetXMLValue(psRoot, "rasterXSize", "0")); const int nYSize = atoi(CPLGetXMLValue(psRoot, "rasterYSize", "0")); - if (!bIsPansharpened && CPLGetXMLNode(psRoot, "VRTRasterBand") != nullptr && + if (!bIsPansharpened && !bIsProcessed && + CPLGetXMLNode(psRoot, "VRTRasterBand") != nullptr && !GDALCheckDatasetDimensions(nXSize, nYSize)) { return nullptr; @@ -1639,6 +1661,8 @@ GDALDataset *VRTDataset::OpenXML(const char *pszXML, const char *pszVRTPath, poDS = new VRTWarpedDataset(nXSize, nYSize); else if (bIsPansharpened) poDS = new VRTPansharpenedDataset(nXSize, nYSize); + else if (bIsProcessed) + poDS = new VRTProcessedDataset(nXSize, nYSize); else { poDS = new VRTDataset(nXSize, nYSize); @@ -2874,4 +2898,94 @@ void VRTDataset::ClearStatistics() GDALDataset::ClearStatistics(); } +/************************************************************************/ +/* BuildSourceFilename() */ +/************************************************************************/ + +/* static */ +std::string VRTDataset::BuildSourceFilename(const char *pszFilename, + const char *pszVRTPath, + bool bRelativeToVRT) +{ + std::string osSrcDSName; + if (pszVRTPath != nullptr && bRelativeToVRT) + { + // Try subdatasetinfo API first + // Note: this will become the only branch when subdatasetinfo will become + // available for NITF_IM, RASTERLITE and TILEDB + const auto oSubDSInfo{GDALGetSubdatasetInfo(pszFilename)}; + if (oSubDSInfo && !oSubDSInfo->GetPathComponent().empty()) + { + auto path{oSubDSInfo->GetPathComponent()}; + osSrcDSName = oSubDSInfo->ModifyPathComponent( + CPLProjectRelativeFilename(pszVRTPath, path.c_str())); + GDALDestroySubdatasetInfo(oSubDSInfo); + } + else + { + bool bDone = false; + for (const char *pszSyntax : VRTDataset::apszSpecialSyntax) + { + CPLString osPrefix(pszSyntax); + osPrefix.resize(strchr(pszSyntax, ':') - pszSyntax + 1); + if (pszSyntax[osPrefix.size()] == '"') + osPrefix += '"'; + if (EQUALN(pszFilename, osPrefix, osPrefix.size())) + { + if (STARTS_WITH_CI(pszSyntax + osPrefix.size(), "{ANY}")) + { + const char *pszLastPart = strrchr(pszFilename, ':') + 1; + // CSV:z:/foo.xyz + if ((pszLastPart[0] == '/' || pszLastPart[0] == '\\') && + pszLastPart - pszFilename >= 3 && + pszLastPart[-3] == ':') + { + pszLastPart -= 2; + } + CPLString osPrefixFilename = pszFilename; + osPrefixFilename.resize(pszLastPart - pszFilename); + osSrcDSName = + osPrefixFilename + + CPLProjectRelativeFilename(pszVRTPath, pszLastPart); + bDone = true; + } + else if (STARTS_WITH_CI(pszSyntax + osPrefix.size(), + "{FILENAME}")) + { + CPLString osFilename(pszFilename + osPrefix.size()); + size_t nPos = 0; + if (osFilename.size() >= 3 && osFilename[1] == ':' && + (osFilename[2] == '\\' || osFilename[2] == '/')) + nPos = 2; + nPos = osFilename.find( + pszSyntax[osPrefix.size() + strlen("{FILENAME}")], + nPos); + if (nPos != std::string::npos) + { + const CPLString osSuffix = osFilename.substr(nPos); + osFilename.resize(nPos); + osSrcDSName = osPrefix + + CPLProjectRelativeFilename( + pszVRTPath, osFilename) + + osSuffix; + bDone = true; + } + } + break; + } + } + if (!bDone) + { + osSrcDSName = + CPLProjectRelativeFilename(pszVRTPath, pszFilename); + } + } + } + else + { + osSrcDSName = pszFilename; + } + return osSrcDSName; +} + /*! @endcond */ diff --git a/frmts/vrt/vrtdataset.h b/frmts/vrt/vrtdataset.h index 6c123027779a..65c68ee9b75e 100644 --- a/frmts/vrt/vrtdataset.h +++ b/frmts/vrt/vrtdataset.h @@ -47,6 +47,7 @@ #include CPLErr GDALRegisterDefaultPixelFunc(); +void GDALVRTRegisterDefaultProcessedDatasetFuncs(); CPLString VRTSerializeNoData(double dfVal, GDALDataType eDataType, int nPrecision); @@ -198,6 +199,7 @@ template struct VRTFlushCacheStruct class VRTWarpedDataset; class VRTPansharpenedDataset; +class VRTProcessedDataset; class VRTGroup; class CPL_DLL VRTDataset CPL_NON_FINAL : public GDALDataset @@ -206,17 +208,15 @@ class CPL_DLL VRTDataset CPL_NON_FINAL : public GDALDataset friend struct VRTFlushCacheStruct; friend struct VRTFlushCacheStruct; friend struct VRTFlushCacheStruct; + friend struct VRTFlushCacheStruct; friend class VRTSourcedRasterBand; + friend class VRTSimpleSource; friend VRTDatasetH CPL_STDCALL VRTCreate(int nXSize, int nYSize); - OGRSpatialReference *m_poSRS = nullptr; - - int m_bGeoTransformSet = false; - double m_adfGeoTransform[6]; - int m_nGCPCount = 0; GDAL_GCP *m_pasGCPList = nullptr; - OGRSpatialReference *m_poGCP_SRS = nullptr; + std::unique_ptr + m_poGCP_SRS{}; bool m_bNeedsFlush = false; bool m_bWritable = true; @@ -244,8 +244,13 @@ class CPL_DLL VRTDataset CPL_NON_FINAL : public GDALDataset VRTSource::WorkingState m_oWorkingState{}; + static constexpr const char *const apszSpecialSyntax[] = { + "NITF_IM:{ANY}:{FILENAME}", "PDF:{ANY}:{FILENAME}", + "RASTERLITE:{FILENAME},{ANY}", "TILEDB:\"{FILENAME}\":{ANY}", + "TILEDB:{FILENAME}:{ANY}"}; + VRTRasterBand *InitBand(const char *pszSubclass, int nBand, - bool bAllowPansharpened); + bool bAllowPansharpenedOrProcessed); static GDALDataset *OpenVRTProtocol(const char *pszSpec); bool AddVirtualOverview(int nOvFactor, const char *pszResampling); @@ -260,6 +265,11 @@ class CPL_DLL VRTDataset CPL_NON_FINAL : public GDALDataset int m_nBlockXSize = 0; int m_nBlockYSize = 0; + std::unique_ptr m_poSRS{}; + + int m_bGeoTransformSet = false; + double m_adfGeoTransform[6]; + virtual int CloseDependentDatasets() override; public: @@ -283,7 +293,7 @@ class CPL_DLL VRTDataset CPL_NON_FINAL : public GDALDataset const OGRSpatialReference *GetSpatialRef() const override { - return m_poSRS; + return m_poSRS.get(); } CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override; @@ -300,7 +310,7 @@ class CPL_DLL VRTDataset CPL_NON_FINAL : public GDALDataset virtual int GetGCPCount() override; const OGRSpatialReference *GetGCPSpatialRef() const override { - return m_poGCP_SRS; + return m_poGCP_SRS.get(); } virtual const GDAL_GCP *GetGCPs() override; using GDALDataset::SetGCPs; @@ -366,8 +376,8 @@ class CPL_DLL VRTDataset CPL_NON_FINAL : public GDALDataset static int Identify(GDALOpenInfo *); static GDALDataset *Open(GDALOpenInfo *); - static GDALDataset *OpenXML(const char *, const char * = nullptr, - GDALAccess eAccess = GA_ReadOnly); + static VRTDataset *OpenXML(const char *, const char * = nullptr, + GDALAccess eAccess = GA_ReadOnly); static GDALDataset *Create(const char *pszName, int nXSize, int nYSize, int nBands, GDALDataType eType, char **papszOptions); @@ -376,6 +386,10 @@ class CPL_DLL VRTDataset CPL_NON_FINAL : public GDALDataset CSLConstList papszRootGroupOptions, CSLConstList papszOptions); static CPLErr Delete(const char *pszFilename); + + static std::string BuildSourceFilename(const char *pszFilename, + const char *pszVRTPath, + bool bRelativeToVRT); }; /************************************************************************/ @@ -509,6 +523,130 @@ class VRTPansharpenedDataset final : public VRTDataset } }; +/************************************************************************/ +/* VRTPansharpenedDataset */ +/************************************************************************/ + +/** Specialized implementation of VRTDataset that chains several processing + * steps applied on all bands at a time. + * + * @since 3.9 + */ +class VRTProcessedDataset final : public VRTDataset +{ + public: + VRTProcessedDataset(int nXSize, int nYSize); + ~VRTProcessedDataset() override; + + virtual CPLErr FlushCache(bool bAtClosing) override; + + virtual CPLErr XMLInit(const CPLXMLNode *, const char *) override; + virtual CPLXMLNode *SerializeToXML(const char *pszVRTPath) override; + + void GetBlockSize(int *, int *) const; + + // GByte whose initialization constructor does nothing +#ifdef __GNUC__ +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Weffc++" +#endif + struct NoInitByte + { + GByte value; + // cppcheck-suppress uninitMemberVar + NoInitByte() + { + // do nothing + /* coverity[uninit_member] */ + } + inline operator GByte() const + { + return value; + } + }; +#ifdef __GNUC__ +#pragma GCC diagnostic pop +#endif + + private: + friend class VRTProcessedRasterBand; + + //! Data for a processing step. + struct Step + { + //! Algorithm name + std::string osAlgorithm{}; + + //! Arguments to pass to the processing function. + CPLStringList aosArguments{}; + + //! Data type of the input buffer. + GDALDataType eInDT = GDT_Unknown; + + //! Data type of the output buffer. + GDALDataType eOutDT = GDT_Unknown; + + //! Number of input bands. + int nInBands = 0; + + //! Number of output bands. + int nOutBands = 0; + + //! Nodata values (nInBands) of the input bands. + std::vector adfInNoData{}; + + //! Nodata values (nOutBands) of the output bands. + std::vector adfOutNoData{}; + + //! Working data structure (private data of the implementation of the function) + VRTPDWorkingDataPtr pWorkingData = nullptr; + + // NOTE: if adding a new member, edit the move constructor and + // assignment operators! + + Step() = default; + ~Step(); + Step(Step &&); + Step &operator=(Step &&); + + private: + Step(const Step &) = delete; + Step &operator=(const Step &) = delete; + void deinit(); + }; + + //! Directory of the VRT + std::string m_osVRTPath{}; + + //! Source dataset + std::unique_ptr m_poSrcDS{}; + + //! Processing steps. + std::vector m_aoSteps{}; + + //! Backup XML tree passed to XMLInit() + CPLXMLTreeCloser m_oXMLTree{nullptr}; + + //! Overview datasets (dynamically generated from the ones of m_poSrcDS) + std::vector> m_apoOverviewDatasets{}; + + //! Input buffer of a processing step + std::vector m_abyInput{}; + + //! Output buffer of a processing step + std::vector m_abyOutput{}; + + CPLErr Init(const CPLXMLNode *, const char *, + const VRTProcessedDataset *poParentDS, + GDALDataset *poParentSrcDS, int iOvrLevel); + + bool ParseStep(const CPLXMLNode *psStep, bool bIsFinalStep, + GDALDataType &eCurrentDT, int &nCurrentBandCount, + std::vector &adfInNoData, + std::vector &adfOutNoData); + bool ProcessRegion(int nXOff, int nYOff, int nBufXSize, int nBufYSize); +}; + /************************************************************************/ /* VRTRasterBand */ /* */ @@ -822,6 +960,7 @@ class CPL_DLL VRTWarpedRasterBand final : public VRTRasterBand virtual int GetOverviewCount() override; virtual GDALRasterBand *GetOverview(int) override; }; + /************************************************************************/ /* VRTPansharpenedRasterBand */ /************************************************************************/ @@ -865,6 +1004,26 @@ class VRTPansharpenedRasterBand final : public VRTRasterBand } }; +/************************************************************************/ +/* VRTProcessedRasterBand */ +/************************************************************************/ + +class VRTProcessedRasterBand final : public VRTRasterBand +{ + public: + VRTProcessedRasterBand(VRTProcessedDataset *poDS, int nBand, + GDALDataType eDataType = GDT_Unknown); + + virtual CPLErr IReadBlock(int, int, void *) override; + + virtual int GetOverviewCount() override; + virtual GDALRasterBand *GetOverview(int) override; + + virtual CPLXMLNode *SerializeToXML(const char *pszVRTPath, + bool &bHasWarnedAboutRAMUsage, + size_t &nAccRAMUsage) override; +}; + /************************************************************************/ /* VRTDerivedRasterBand */ /************************************************************************/ diff --git a/frmts/vrt/vrtdriver.cpp b/frmts/vrt/vrtdriver.cpp index 2883944b4e13..99fa8ca7c4da 100644 --- a/frmts/vrt/vrtdriver.cpp +++ b/frmts/vrt/vrtdriver.cpp @@ -507,6 +507,9 @@ void GDALRegister_VRT() // First register the pixel functions GDALRegisterDefaultPixelFunc(); + // Registr functions for VRTProcessedDataset + GDALVRTRegisterDefaultProcessedDatasetFuncs(); + VRTDriver *poDriver = new VRTDriver(); poDriver->SetDescription("VRT"); diff --git a/frmts/vrt/vrtprocesseddataset.cpp b/frmts/vrt/vrtprocesseddataset.cpp new file mode 100644 index 000000000000..06be6006bbc7 --- /dev/null +++ b/frmts/vrt/vrtprocesseddataset.cpp @@ -0,0 +1,1341 @@ +/****************************************************************************** + * + * Project: Virtual GDAL Datasets + * Purpose: Implementation of VRTProcessedDataset. + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2024, Even Rouault + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +#include "cpl_minixml.h" +#include "cpl_string.h" +#include "vrtdataset.h" + +#include +#include +#include +#include + +/************************************************************************/ +/* VRTProcessedDatasetFunc */ +/************************************************************************/ + +//! Structure holding information for a VRTProcessedDataset function. +struct VRTProcessedDatasetFunc +{ + //! Processing function name + std::string osFuncName{}; + + //! User data to provide to pfnInit, pfnFree, pfnProcess callbacks. + void *pUserData = nullptr; + + //! Whether XML metadata has been specified + bool bMetadataSpecified = false; + + //! Map of (constant argument name, constant value) + std::map oMapConstantArguments{}; + + //! Set of builtin argument names (e.g "offset", "scale", "nodata") + std::set oSetBuiltinArguments{}; + + //! Arguments defined in the VRT + struct OtherArgument + { + std::string osType{}; + bool bRequired = false; + }; + std::map oOtherArguments{}; + + //! Requested input data type. + GDALDataType eRequestedInputDT = GDT_Unknown; + + //! List of supported input datatypes. Empty if no restriction. + std::vector aeSupportedInputDT{}; + + //! List of supported input band counts. Empty if no restriction. + std::vector anSupportedInputBandCount{}; + + //! Optional initialization function + GDALVRTProcessedDatasetFuncInit pfnInit = nullptr; + + //! Optional free function + GDALVRTProcessedDatasetFuncFree pfnFree = nullptr; + + //! Required processing function + GDALVRTProcessedDatasetFuncProcess pfnProcess = nullptr; +}; + +/************************************************************************/ +/* GetGlobalMapProcessedDatasetFunc() */ +/************************************************************************/ + +/** Return the registry of VRTProcessedDatasetFunc functions */ +static std::map & +GetGlobalMapProcessedDatasetFunc() +{ + static std::map goMap; + return goMap; +} + +/************************************************************************/ +/* Step::~Step() */ +/************************************************************************/ + +/*! @cond Doxygen_Suppress */ + +/** Step destructor */ +VRTProcessedDataset::Step::~Step() +{ + deinit(); +} + +/************************************************************************/ +/* Step::deinit() */ +/************************************************************************/ + +/** Free pWorkingData */ +void VRTProcessedDataset::Step::deinit() +{ + if (pWorkingData) + { + const auto &oMapFunctions = GetGlobalMapProcessedDatasetFunc(); + const auto oIterFunc = oMapFunctions.find(osAlgorithm); + if (oIterFunc != oMapFunctions.end()) + { + if (oIterFunc->second.pfnFree) + { + oIterFunc->second.pfnFree(osAlgorithm.c_str(), + oIterFunc->second.pUserData, + pWorkingData); + } + } + else + { + CPLAssert(false); + } + pWorkingData = nullptr; + } +} + +/************************************************************************/ +/* Step::Step(Step&& other) */ +/************************************************************************/ + +/** Move constructor */ +VRTProcessedDataset::Step::Step(Step &&other) + : osAlgorithm(std::move(other.osAlgorithm)), + aosArguments(std::move(other.aosArguments)), eInDT(other.eInDT), + eOutDT(other.eOutDT), nInBands(other.nInBands), + nOutBands(other.nOutBands), adfInNoData(other.adfInNoData), + adfOutNoData(other.adfOutNoData), pWorkingData(other.pWorkingData) +{ + other.pWorkingData = nullptr; +} + +/************************************************************************/ +/* Step operator=(Step&& other) */ +/************************************************************************/ + +/** Move assignment operator */ +VRTProcessedDataset::Step &VRTProcessedDataset::Step::operator=(Step &&other) +{ + if (&other != this) + { + deinit(); + osAlgorithm = std::move(other.osAlgorithm); + aosArguments = std::move(other.aosArguments); + eInDT = other.eInDT; + eOutDT = other.eOutDT; + nInBands = other.nInBands; + nOutBands = other.nOutBands; + adfInNoData = std::move(other.adfInNoData); + adfOutNoData = std::move(other.adfOutNoData); + std::swap(pWorkingData, other.pWorkingData); + } + return *this; +} + +/************************************************************************/ +/* VRTProcessedDataset() */ +/************************************************************************/ + +/** Constructor */ +VRTProcessedDataset::VRTProcessedDataset(int nXSize, int nYSize) + : VRTDataset(nXSize, nYSize) +{ +} + +/************************************************************************/ +/* ~VRTProcessedDataset() */ +/************************************************************************/ + +VRTProcessedDataset::~VRTProcessedDataset() + +{ + VRTProcessedDataset::FlushCache(true); + VRTProcessedDataset::CloseDependentDatasets(); +} + +/************************************************************************/ +/* XMLInit() */ +/************************************************************************/ + +/** Instantiate object from XML tree */ +CPLErr VRTProcessedDataset::XMLInit(const CPLXMLNode *psTree, + const char *pszVRTPathIn) + +{ + if (Init(psTree, pszVRTPathIn, nullptr, nullptr, -1) != CE_None) + return CE_Failure; + + const auto poSrcFirstBand = m_poSrcDS->GetRasterBand(1); + const int nOvrCount = poSrcFirstBand->GetOverviewCount(); + for (int i = 0; i < nOvrCount; ++i) + { + auto poOvrDS = std::make_unique(0, 0); + if (poOvrDS->Init(psTree, pszVRTPathIn, this, m_poSrcDS.get(), i) != + CE_None) + break; + m_apoOverviewDatasets.emplace_back(std::move(poOvrDS)); + } + + return CE_None; +} + +/** Instantiate object from XML tree */ +CPLErr VRTProcessedDataset::Init(const CPLXMLNode *psTree, + const char *pszVRTPathIn, + const VRTProcessedDataset *poParentDS, + GDALDataset *poParentSrcDS, int iOvrLevel) + +{ + const CPLXMLNode *psInput = CPLGetXMLNode(psTree, "Input"); + if (!psInput) + { + CPLError(CE_Failure, CPLE_AppDefined, "Input element missing"); + return CE_Failure; + } + + if (pszVRTPathIn) + m_osVRTPath = pszVRTPathIn; + + if (poParentSrcDS) + { + m_poSrcDS.reset( + GDALCreateOverviewDataset(poParentSrcDS, iOvrLevel, true)); + } + else if (const CPLXMLNode *psSourceFileNameNode = + CPLGetXMLNode(psInput, "SourceFilename")) + { + const bool bRelativeToVRT = CPL_TO_BOOL( + atoi(CPLGetXMLValue(psSourceFileNameNode, "relativetoVRT", "0"))); + const std::string osFilename = VRTDataset::BuildSourceFilename( + CPLGetXMLValue(psInput, "SourceFilename", ""), pszVRTPathIn, + bRelativeToVRT); + m_poSrcDS.reset(GDALDataset::Open( + osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr, + nullptr, nullptr)); + } + else if (const CPLXMLNode *psVRTDataset = + CPLGetXMLNode(psInput, "VRTDataset")) + { + CPLXMLNode sVRTDatasetTmp = *psVRTDataset; + sVRTDatasetTmp.psNext = nullptr; + char *pszXML = CPLSerializeXMLTree(&sVRTDatasetTmp); + m_poSrcDS.reset(VRTDataset::OpenXML(pszXML, pszVRTPathIn, GA_ReadOnly)); + CPLFree(pszXML); + } + else + { + CPLError( + CE_Failure, CPLE_AppDefined, + "Input element should have a SourceFilename or VRTDataset element"); + return CE_Failure; + } + + if (!m_poSrcDS) + return CE_Failure; + + if (nRasterXSize == 0 && nRasterYSize == 0) + { + nRasterXSize = m_poSrcDS->GetRasterXSize(); + nRasterYSize = m_poSrcDS->GetRasterYSize(); + } + else if (nRasterXSize != m_poSrcDS->GetRasterXSize() || + nRasterYSize != m_poSrcDS->GetRasterYSize()) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Inconsistent declared VRT dimensions with input dataset"); + return CE_Failure; + } + + if (m_poSrcDS->GetRasterCount() == 0) + return CE_Failure; + + // Inherit SRS from source if not explicitly defined in VRT + if (!CPLGetXMLNode(psTree, "SRS")) + { + const OGRSpatialReference *poSRS = m_poSrcDS->GetSpatialRef(); + if (poSRS) + { + m_poSRS.reset(poSRS->Clone()); + } + } + + // Inherit GeoTransform from source if not explicitly defined in VRT + if (iOvrLevel < 0 && !CPLGetXMLNode(psTree, "GeoTransform")) + { + if (m_poSrcDS->GetGeoTransform(m_adfGeoTransform) == CE_None) + m_bGeoTransformSet = true; + } + + /* -------------------------------------------------------------------- */ + /* Initialize blocksize before calling sub-init so that the */ + /* band initializers can get it from the dataset object when */ + /* they are created. */ + /* -------------------------------------------------------------------- */ + + const auto poSrcFirstBand = m_poSrcDS->GetRasterBand(1); + poSrcFirstBand->GetBlockSize(&m_nBlockXSize, &m_nBlockYSize); + if (const char *pszBlockXSize = + CPLGetXMLValue(psTree, "BlockXSize", nullptr)) + m_nBlockXSize = atoi(pszBlockXSize); + if (const char *pszBlockYSize = + CPLGetXMLValue(psTree, "BlockYSize", nullptr)) + m_nBlockYSize = atoi(pszBlockYSize); + + // Initialize all the general VRT stuff. + if (VRTDataset::XMLInit(psTree, pszVRTPathIn) != CE_None) + { + return CE_Failure; + } + + // Use geotransform from parent for overviews + if (iOvrLevel >= 0 && poParentDS->m_bGeoTransformSet) + { + m_bGeoTransformSet = true; + m_adfGeoTransform[0] = poParentDS->m_adfGeoTransform[0]; + m_adfGeoTransform[1] = poParentDS->m_adfGeoTransform[1]; + m_adfGeoTransform[2] = poParentDS->m_adfGeoTransform[2]; + m_adfGeoTransform[3] = poParentDS->m_adfGeoTransform[3]; + m_adfGeoTransform[4] = poParentDS->m_adfGeoTransform[4]; + m_adfGeoTransform[5] = poParentDS->m_adfGeoTransform[5]; + + m_adfGeoTransform[1] *= + static_cast(poParentDS->GetRasterXSize()) / nRasterXSize; + m_adfGeoTransform[2] *= + static_cast(poParentDS->GetRasterYSize()) / nRasterYSize; + m_adfGeoTransform[4] *= + static_cast(poParentDS->GetRasterXSize()) / nRasterXSize; + m_adfGeoTransform[5] *= + static_cast(poParentDS->GetRasterYSize()) / nRasterYSize; + } + + // Create bands automatically from source dataset if not explicitly defined + // in VRT. + if (!CPLGetXMLNode(psTree, "VRTRasterBand")) + { + for (int i = 0; i < m_poSrcDS->GetRasterCount(); ++i) + { + const auto poSrcBand = m_poSrcDS->GetRasterBand(i + 1); + auto poBand = new VRTProcessedRasterBand( + this, i + 1, poSrcBand->GetRasterDataType()); + poBand->CopyCommonInfoFrom(poSrcBand); + SetBand(i + 1, poBand); + } + } + + const CPLXMLNode *psProcessingSteps = + CPLGetXMLNode(psTree, "ProcessingSteps"); + if (!psProcessingSteps) + { + CPLError(CE_Failure, CPLE_AppDefined, + "ProcessingSteps element missing"); + return CE_Failure; + } + + const auto eInDT = poSrcFirstBand->GetRasterDataType(); + for (int i = 1; i < m_poSrcDS->GetRasterCount(); ++i) + { + const auto eDT = m_poSrcDS->GetRasterBand(i + 1)->GetRasterDataType(); + if (eDT != eInDT) + { + CPLError(CE_Warning, CPLE_AppDefined, + "Not all bands of the input dataset have the same data " + "type. The data type of the first band will be used as " + "the reference one."); + break; + } + } + + GDALDataType eCurrentDT = eInDT; + int nCurrentBandCount = m_poSrcDS->GetRasterCount(); + + std::vector adfNoData; + for (int i = 1; i <= nCurrentBandCount; ++i) + { + int bHasVal = FALSE; + const double dfVal = + m_poSrcDS->GetRasterBand(i)->GetNoDataValue(&bHasVal); + adfNoData.emplace_back( + bHasVal ? dfVal : std::numeric_limits::quiet_NaN()); + } + + int nStepCount = 0; + for (const CPLXMLNode *psStep = psProcessingSteps->psChild; psStep; + psStep = psStep->psNext) + { + if (psStep->eType == CXT_Element && + strcmp(psStep->pszValue, "Step") == 0) + { + ++nStepCount; + } + } + + int iStep = 0; + for (const CPLXMLNode *psStep = psProcessingSteps->psChild; psStep; + psStep = psStep->psNext) + { + if (psStep->eType == CXT_Element && + strcmp(psStep->pszValue, "Step") == 0) + { + ++iStep; + const bool bIsFinalStep = (iStep == nStepCount); + std::vector adfOutNoData; + if (bIsFinalStep) + { + // Initialize adfOutNoData with nodata value of *output* bands + // for final step + for (int i = 1; i <= nBands; ++i) + { + int bHasVal = FALSE; + const double dfVal = + GetRasterBand(i)->GetNoDataValue(&bHasVal); + adfOutNoData.emplace_back( + bHasVal ? dfVal + : std::numeric_limits::quiet_NaN()); + } + } + if (!ParseStep(psStep, bIsFinalStep, eCurrentDT, nCurrentBandCount, + adfNoData, adfOutNoData)) + return CE_Failure; + adfNoData = std::move(adfOutNoData); + } + } + + if (m_aoSteps.empty()) + { + CPLError(CE_Failure, CPLE_AppDefined, + "At least one step should be defined"); + return CE_Failure; + } + + if (nCurrentBandCount != nBands) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Number of output bands of last step is not consistent with " + "number of VRTProcessedRasterBand's"); + return CE_Failure; + } + + if (nBands > 1) + SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE"); + + m_oXMLTree.reset(CPLCloneXMLTree(psTree)); + + return CE_None; +} + +/************************************************************************/ +/* ParseStep() */ +/************************************************************************/ + +/** Parse the current Step node and create a corresponding entry in m_aoSteps. + * + * @param psStep Step node + * @param bIsFinalStep Whether this is the final step. + * @param[in,out] eCurrentDT Input data type for this step. + * Updated to output data type at end of method. + * @param[in,out] nCurrentBandCount Input band count for this step. + * Updated to output band cout at end of + * method. + * @param adfInNoData Input nodata values + * @param[in,out] adfOutNoData Output nodata values, to be filled by this + * method. When bIsFinalStep, this is also an + * input parameter. + * @return true on success. + */ +bool VRTProcessedDataset::ParseStep(const CPLXMLNode *psStep, bool bIsFinalStep, + GDALDataType &eCurrentDT, + int &nCurrentBandCount, + std::vector &adfInNoData, + std::vector &adfOutNoData) +{ + const char *pszStepName = CPLGetXMLValue( + psStep, "name", CPLSPrintf("nr %d", 1 + int(m_aoSteps.size()))); + const char *pszAlgorithm = CPLGetXMLValue(psStep, "Algorithm", nullptr); + if (!pszAlgorithm) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Step '%s' lacks a Algorithm element", pszStepName); + return false; + } + + const auto &oMapFunctions = GetGlobalMapProcessedDatasetFunc(); + const auto oIterFunc = oMapFunctions.find(pszAlgorithm); + if (oIterFunc == oMapFunctions.end()) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Step '%s' uses unregistered algorithm '%s'", pszStepName, + pszAlgorithm); + return false; + } + + const auto &oFunc = oIterFunc->second; + + if (!oFunc.aeSupportedInputDT.empty()) + { + if (std::find(oFunc.aeSupportedInputDT.begin(), + oFunc.aeSupportedInputDT.end(), + eCurrentDT) == oFunc.aeSupportedInputDT.end()) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Step '%s' (using algorithm '%s') does not " + "support input data type = '%s'", + pszStepName, pszAlgorithm, + GDALGetDataTypeName(eCurrentDT)); + return false; + } + } + + if (!oFunc.anSupportedInputBandCount.empty()) + { + if (std::find(oFunc.anSupportedInputBandCount.begin(), + oFunc.anSupportedInputBandCount.end(), + nCurrentBandCount) == + oFunc.anSupportedInputBandCount.end()) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Step '%s' (using algorithm '%s') does not " + "support input band count = %d", + pszStepName, pszAlgorithm, nCurrentBandCount); + return false; + } + } + + Step oStep; + oStep.osAlgorithm = pszAlgorithm; + oStep.eInDT = oFunc.eRequestedInputDT != GDT_Unknown + ? oFunc.eRequestedInputDT + : eCurrentDT; + oStep.nInBands = nCurrentBandCount; + + // Unless modified by pfnInit... + oStep.eOutDT = oStep.eInDT; + + oStep.adfInNoData = adfInNoData; + oStep.adfOutNoData = bIsFinalStep ? adfOutNoData : adfInNoData; + + // Deal with constant arguments + for (const auto &nameValuePair : oFunc.oMapConstantArguments) + { + oStep.aosArguments.AddNameValue(nameValuePair.first.c_str(), + nameValuePair.second.c_str()); + } + + // Deal with built-in arguments + if (oFunc.oSetBuiltinArguments.find("nodata") != + oFunc.oSetBuiltinArguments.end()) + { + int bHasVal = false; + const auto poSrcFirstBand = m_poSrcDS->GetRasterBand(1); + const double dfVal = poSrcFirstBand->GetNoDataValue(&bHasVal); + if (bHasVal) + { + oStep.aosArguments.AddNameValue("nodata", + CPLSPrintf("%.18g", dfVal)); + } + } + + if (oFunc.oSetBuiltinArguments.find("offset_{band}") != + oFunc.oSetBuiltinArguments.end()) + { + for (int i = 1; i <= m_poSrcDS->GetRasterCount(); ++i) + { + int bHasVal = false; + const double dfVal = GetRasterBand(i)->GetOffset(&bHasVal); + oStep.aosArguments.AddNameValue( + CPLSPrintf("offset_%d", i), + CPLSPrintf("%.18g", bHasVal ? dfVal : 0.0)); + } + } + + if (oFunc.oSetBuiltinArguments.find("scale_{band}") != + oFunc.oSetBuiltinArguments.end()) + { + for (int i = 1; i <= m_poSrcDS->GetRasterCount(); ++i) + { + int bHasVal = false; + const double dfVal = GetRasterBand(i)->GetScale(&bHasVal); + oStep.aosArguments.AddNameValue( + CPLSPrintf("scale_%d", i), + CPLSPrintf("%.18g", bHasVal ? dfVal : 1.0)); + } + } + + // Parse arguments specified in VRT + std::set oFoundArguments; + + for (const CPLXMLNode *psStepChild = psStep->psChild; psStepChild; + psStepChild = psStepChild->psNext) + { + if (psStepChild->eType == CXT_Element && + strcmp(psStepChild->pszValue, "Argument") == 0) + { + const char *pszParamName = + CPLGetXMLValue(psStepChild, "name", nullptr); + if (!pszParamName) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Step '%s' has a Argument without a name attribute", + pszStepName); + return false; + } + const char *pszValue = CPLGetXMLValue(psStepChild, nullptr, ""); + auto oOtherArgIter = + oFunc.oOtherArguments.find(CPLString(pszParamName).tolower()); + if (!oFunc.oOtherArguments.empty() && + oOtherArgIter == oFunc.oOtherArguments.end()) + { + // If we got a parameter name like 'coefficients_1', + // try to fetch the generic 'coefficients_{band}' + std::string osParamName(pszParamName); + const auto nPos = osParamName.rfind('_'); + if (nPos != std::string::npos) + { + osParamName.resize(nPos + 1); + osParamName += "{band}"; + oOtherArgIter = oFunc.oOtherArguments.find( + CPLString(osParamName).tolower()); + } + } + if (oOtherArgIter != oFunc.oOtherArguments.end()) + { + oFoundArguments.insert(oOtherArgIter->first); + + const std::string &osType = oOtherArgIter->second.osType; + if (osType == "boolean") + { + if (!EQUAL(pszValue, "true") && !EQUAL(pszValue, "false")) + { + CPLError(CE_Failure, CPLE_NotSupported, + "Step '%s' has a Argument '%s' whose " + "value '%s' is not a boolean", + pszStepName, pszParamName, pszValue); + return false; + } + } + else if (osType == "integer") + { + if (CPLGetValueType(pszValue) != CPL_VALUE_INTEGER) + { + CPLError(CE_Failure, CPLE_NotSupported, + "Step '%s' has a Argument '%s' whose " + "value '%s' is not a integer", + pszStepName, pszParamName, pszValue); + return false; + } + } + else if (osType == "double") + { + const auto eType = CPLGetValueType(pszValue); + if (eType != CPL_VALUE_INTEGER && eType != CPL_VALUE_REAL) + { + CPLError(CE_Failure, CPLE_NotSupported, + "Step '%s' has a Argument '%s' whose " + "value '%s' is not a double", + pszStepName, pszParamName, pszValue); + return false; + } + } + else if (osType == "double_list") + { + const CPLStringList aosTokens( + CSLTokenizeString2(pszValue, ",", 0)); + for (int i = 0; i < aosTokens.size(); ++i) + { + const auto eType = CPLGetValueType(aosTokens[i]); + if (eType != CPL_VALUE_INTEGER && + eType != CPL_VALUE_REAL) + { + CPLError(CE_Failure, CPLE_NotSupported, + "Step '%s' has a Argument '%s' " + "whose value '%s' is not a " + "comma-separated list of doubles", + pszStepName, pszParamName, pszValue); + return false; + } + } + } + else if (osType != "string") + { + CPLDebug("VRT", "Unhandled argument type '%s'", + osType.c_str()); + CPLAssert(0); + } + } + else if (oFunc.bMetadataSpecified && + oFunc.oSetBuiltinArguments.find( + CPLString(pszParamName).tolower()) == + oFunc.oSetBuiltinArguments.end() && + oFunc.oMapConstantArguments.find( + CPLString(pszParamName).tolower()) == + oFunc.oMapConstantArguments.end()) + { + CPLError(CE_Warning, CPLE_NotSupported, + "Step '%s' has a Argument '%s' which is not " + "supported", + pszStepName, pszParamName); + } + + oStep.aosArguments.AddNameValue(pszParamName, pszValue); + } + } + + // Check that required arguments have been specified + for (const auto &oIter : oFunc.oOtherArguments) + { + if (oIter.second.bRequired && + oFoundArguments.find(oIter.first) == oFoundArguments.end()) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Step '%s' lacks required Argument '%s'", pszStepName, + oIter.first.c_str()); + return false; + } + } + + if (oFunc.pfnInit) + { + double *padfOutNoData = nullptr; + if (bIsFinalStep) + { + oStep.nOutBands = nBands; + padfOutNoData = + static_cast(CPLMalloc(nBands * sizeof(double))); + CPLAssert(adfOutNoData.size() == static_cast(nBands)); + memcpy(padfOutNoData, adfOutNoData.data(), nBands * sizeof(double)); + } + else + { + oStep.nOutBands = 0; + } + + if (oFunc.pfnInit(pszAlgorithm, oFunc.pUserData, + oStep.aosArguments.List(), oStep.nInBands, + oStep.eInDT, adfInNoData.data(), &(oStep.nOutBands), + &(oStep.eOutDT), &padfOutNoData, m_osVRTPath.c_str(), + &(oStep.pWorkingData)) != CE_None) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Step '%s' (using algorithm '%s') init() function " + "failed", + pszStepName, pszAlgorithm); + CPLFree(padfOutNoData); + return false; + } + + // Input nodata values may have been modified by pfnInit() + oStep.adfInNoData = adfInNoData; + + if (padfOutNoData) + { + adfOutNoData = + std::vector(padfOutNoData, padfOutNoData + nBands); + } + else + { + adfOutNoData = std::vector( + oStep.nOutBands, std::numeric_limits::quiet_NaN()); + } + CPLFree(padfOutNoData); + + oStep.adfOutNoData = adfOutNoData; + } + else + { + oStep.nOutBands = oStep.nInBands; + adfOutNoData = oStep.adfOutNoData; + } + + eCurrentDT = oStep.eOutDT; + nCurrentBandCount = oStep.nOutBands; + + m_aoSteps.emplace_back(std::move(oStep)); + + return true; +} + +/************************************************************************/ +/* SerializeToXML() */ +/************************************************************************/ + +CPLXMLNode *VRTProcessedDataset::SerializeToXML(const char *pszVRTPathIn) + +{ + CPLXMLNode *psTree = CPLCloneXMLTree(m_oXMLTree.get()); + if (psTree == nullptr) + return psTree; + + /* -------------------------------------------------------------------- */ + /* Remove VRTRasterBand nodes from the original tree and find the */ + /* last child. */ + /* -------------------------------------------------------------------- */ + CPLXMLNode *psLastChild = psTree->psChild; + CPLXMLNode *psPrevChild = nullptr; + while (psLastChild) + { + CPLXMLNode *psNextChild = psLastChild->psNext; + if (psLastChild->eType == CXT_Element && + strcmp(psLastChild->pszValue, "VRTRasterBand") == 0) + { + if (psPrevChild) + psPrevChild->psNext = psNextChild; + else + psTree->psChild = psNextChild; + psLastChild->psNext = nullptr; + CPLDestroyXMLNode(psLastChild); + psLastChild = psPrevChild ? psPrevChild : psTree->psChild; + } + else if (!psNextChild) + { + break; + } + else + { + psPrevChild = psLastChild; + psLastChild = psNextChild; + } + } + CPLAssert(psLastChild); // we have at least Input + + /* -------------------------------------------------------------------- */ + /* Serialize bands. */ + /* -------------------------------------------------------------------- */ + bool bHasWarnedAboutRAMUsage = false; + size_t nAccRAMUsage = 0; + for (int iBand = 0; iBand < nBands; iBand++) + { + CPLXMLNode *psBandTree = + static_cast(papoBands[iBand]) + ->SerializeToXML(pszVRTPathIn, bHasWarnedAboutRAMUsage, + nAccRAMUsage); + + if (psBandTree != nullptr) + { + psLastChild->psNext = psBandTree; + psLastChild = psBandTree; + } + } + + return psTree; +} + +/************************************************************************/ +/* SerializeToXML() */ +/************************************************************************/ + +CPLXMLNode * +VRTProcessedRasterBand::SerializeToXML(const char *pszVRTPathIn, + bool &bHasWarnedAboutRAMUsage, + size_t &nAccRAMUsage) + +{ + CPLXMLNode *psTree = VRTRasterBand::SerializeToXML( + pszVRTPathIn, bHasWarnedAboutRAMUsage, nAccRAMUsage); + + /* -------------------------------------------------------------------- */ + /* Set subclass. */ + /* -------------------------------------------------------------------- */ + CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"), + CXT_Text, "VRTProcessedRasterBand"); + + return psTree; +} + +/************************************************************************/ +/* GetBlockSize() */ +/************************************************************************/ + +/** Return block size */ +void VRTProcessedDataset::GetBlockSize(int *pnBlockXSize, + int *pnBlockYSize) const + +{ + *pnBlockXSize = m_nBlockXSize; + *pnBlockYSize = m_nBlockYSize; +} + +/************************************************************************/ +/* ProcessRegion() */ +/************************************************************************/ + +/** Compute pixel values for the specified region. + * + * The output is stored in m_abyInput in a pixel-interleaved way. + */ +bool VRTProcessedDataset::ProcessRegion(int nXOff, int nYOff, int nBufXSize, + int nBufYSize) +{ + + CPLAssert(!m_aoSteps.empty()); + + const int nFirstBandCount = m_aoSteps.front().nInBands; + CPLAssert(nFirstBandCount == m_poSrcDS->GetRasterCount()); + const GDALDataType eFirstDT = m_aoSteps.front().eInDT; + const int nFirstDTSize = GDALGetDataTypeSizeBytes(eFirstDT); + auto &abyInput = m_abyInput; + auto &abyOutput = m_abyOutput; + try + { + abyInput.resize(static_cast(nBufXSize) * nBufYSize * + nFirstBandCount * nFirstDTSize); + } + catch (const std::bad_alloc &) + { + CPLError(CE_Failure, CPLE_OutOfMemory, + "Out of memory allocating working buffer"); + return false; + } + + if (m_poSrcDS->RasterIO(GF_Read, nXOff, nYOff, nBufXSize, nBufYSize, + abyInput.data(), nBufXSize, nBufYSize, eFirstDT, + nFirstBandCount, nullptr, + nFirstDTSize * nFirstBandCount, + nFirstDTSize * nFirstBandCount * nBufXSize, + nFirstDTSize, nullptr) != CE_None) + { + return false; + } + + const double dfSrcXOff = nXOff; + const double dfSrcYOff = nYOff; + const double dfSrcXSize = nBufXSize; + const double dfSrcYSize = nBufYSize; + + double adfSrcGT[6]; + if (m_poSrcDS->GetGeoTransform(adfSrcGT) != CE_None) + { + adfSrcGT[0] = 0; + adfSrcGT[1] = 1; + adfSrcGT[2] = 0; + adfSrcGT[3] = 0; + adfSrcGT[4] = 0; + adfSrcGT[5] = 1; + } + + GDALDataType eLastDT = eFirstDT; + const auto &oMapFunctions = GetGlobalMapProcessedDatasetFunc(); + for (const auto &oStep : m_aoSteps) + { + const auto oIterFunc = oMapFunctions.find(oStep.osAlgorithm); + CPLAssert(oIterFunc != oMapFunctions.end()); + + // Data type adaptation + if (eLastDT != oStep.eInDT) + { + try + { + abyOutput.resize(static_cast(nBufXSize) * nBufYSize * + oStep.nInBands * + GDALGetDataTypeSizeBytes(oStep.eInDT)); + } + catch (const std::bad_alloc &) + { + CPLError(CE_Failure, CPLE_OutOfMemory, + "Out of memory allocating working buffer"); + return false; + } + + GDALCopyWords64(abyInput.data(), eLastDT, + GDALGetDataTypeSizeBytes(eLastDT), abyOutput.data(), + oStep.eInDT, GDALGetDataTypeSizeBytes(oStep.eInDT), + static_cast(nBufXSize) * nBufYSize * + oStep.nInBands); + + std::swap(abyInput, abyOutput); + } + + try + { + abyOutput.resize(static_cast(nBufXSize) * nBufYSize * + oStep.nOutBands * + GDALGetDataTypeSizeBytes(oStep.eOutDT)); + } + catch (const std::bad_alloc &) + { + CPLError(CE_Failure, CPLE_OutOfMemory, + "Out of memory allocating working buffer"); + return false; + } + + const auto &oFunc = oIterFunc->second; + if (oFunc.pfnProcess( + oStep.osAlgorithm.c_str(), oFunc.pUserData, oStep.pWorkingData, + oStep.aosArguments.List(), nBufXSize, nBufYSize, + abyInput.data(), abyInput.size(), oStep.eInDT, oStep.nInBands, + oStep.adfInNoData.data(), abyOutput.data(), abyOutput.size(), + oStep.eOutDT, oStep.nOutBands, oStep.adfOutNoData.data(), + dfSrcXOff, dfSrcYOff, dfSrcXSize, dfSrcYSize, adfSrcGT, + m_osVRTPath.c_str(), + /*papszExtra=*/nullptr) != CE_None) + { + return false; + } + + std::swap(abyInput, abyOutput); + eLastDT = oStep.eOutDT; + } + + return true; +} + +/************************************************************************/ +/* VRTProcessedRasterBand() */ +/************************************************************************/ + +/** Constructor */ +VRTProcessedRasterBand::VRTProcessedRasterBand(VRTProcessedDataset *poDSIn, + int nBandIn, + GDALDataType eDataTypeIn) +{ + Initialize(poDSIn->GetRasterXSize(), poDSIn->GetRasterYSize()); + + poDS = poDSIn; + nBand = nBandIn; + eAccess = GA_Update; + eDataType = eDataTypeIn; + + poDSIn->GetBlockSize(&nBlockXSize, &nBlockYSize); +} + +/************************************************************************/ +/* GetOverviewCount() */ +/************************************************************************/ + +int VRTProcessedRasterBand::GetOverviewCount() +{ + auto poVRTDS = cpl::down_cast(poDS); + return static_cast(poVRTDS->m_apoOverviewDatasets.size()); +} + +/************************************************************************/ +/* GetOverview() */ +/************************************************************************/ + +GDALRasterBand *VRTProcessedRasterBand::GetOverview(int iOvr) +{ + auto poVRTDS = cpl::down_cast(poDS); + if (iOvr < 0 || + iOvr >= static_cast(poVRTDS->m_apoOverviewDatasets.size())) + return nullptr; + return poVRTDS->m_apoOverviewDatasets[iOvr]->GetRasterBand(nBand); +} + +/************************************************************************/ +/* IReadBlock() */ +/************************************************************************/ + +CPLErr VRTProcessedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, + void *pImage) + +{ + auto poVRTDS = cpl::down_cast(poDS); + + int nBufXSize = 0; + int nBufYSize = 0; + GetActualBlockSize(nBlockXOff, nBlockYOff, &nBufXSize, &nBufYSize); + + const int nXPixelOff = nBlockXOff * nBlockXSize; + const int nYPixelOff = nBlockYOff * nBlockYSize; + if (!poVRTDS->ProcessRegion(nXPixelOff, nYPixelOff, nBufXSize, nBufYSize)) + { + return CE_Failure; + } + + const int nOutBands = poVRTDS->m_aoSteps.back().nOutBands; + CPLAssert(nOutBands == poVRTDS->GetRasterCount()); + const auto eLastDT = poVRTDS->m_aoSteps.back().eOutDT; + const int nLastDTSize = GDALGetDataTypeSizeBytes(eLastDT); + const int nDTSize = GDALGetDataTypeSizeBytes(eDataType); + + // Dispatch final output buffer to cached blocks of output bands + for (int iDstBand = 0; iDstBand < nOutBands; ++iDstBand) + { + GDALRasterBlock *poBlock = nullptr; + GByte *pDst; + if (iDstBand + 1 == nBand) + { + pDst = static_cast(pImage); + } + else + { + auto poOtherBand = poVRTDS->papoBands[iDstBand]; + poBlock = poOtherBand->TryGetLockedBlockRef(nBlockXOff, nBlockYOff); + if (poBlock) + { + poBlock->DropLock(); + continue; + } + poBlock = poOtherBand->GetLockedBlockRef( + nBlockXOff, nBlockYOff, /* bJustInitialized = */ true); + if (!poBlock) + continue; + pDst = static_cast(poBlock->GetDataRef()); + } + for (int iY = 0; iY < nBufYSize; ++iY) + { + GDALCopyWords(poVRTDS->m_abyInput.data() + + (iDstBand + static_cast(iY) * nBufXSize * + nOutBands) * + nLastDTSize, + eLastDT, nLastDTSize * nOutBands, + pDst + + static_cast(iY) * nBlockXSize * nDTSize, + eDataType, nDTSize, nBufXSize); + } + if (poBlock) + poBlock->DropLock(); + } + + return CE_None; +} + +/*! @endcond */ + +/************************************************************************/ +/* GDALVRTRegisterProcessedDatasetFunc() */ +/************************************************************************/ + +/** Register a function to be used by VRTProcessedDataset. + + An example of content for pszXMLMetadata is: + \verbatim + + + + + + + + \endverbatim + + @param pszFuncName Function name. Must be unique and not null. + @param pUserData User data. May be nullptr. Must remain valid during the + lifetime of GDAL. + @param pszXMLMetadata XML metadata describing the function arguments. May be + nullptr if there are no arguments. + @param eRequestedInputDT If the pfnProcess callback only supports a single + data type, it should be specified in this parameter. + Otherwise set it to GDT_Unknown. + @param paeSupportedInputDT List of supported input data types. May be nullptr + if all are supported or if eRequestedInputDT is + set to a non GDT_Unknown value. + @param nSupportedInputDTSize Size of paeSupportedInputDT + @param panSupportedInputBandCount List of supported band count. May be nullptr + if any source band count is supported. + @param nSupportedInputBandCountSize Size of panSupportedInputBandCount + @param pfnInit Initialization function called when a VRTProcessedDataset + step uses the register function. This initialization function + will return the output data type, output band count and + potentially initialize a working structure, typically parsing + arguments. May be nullptr. + If not specified, it will be assumed that the input and output + data types are the same, and that the input number of bands + and output number of bands are the same. + @param pfnFree Free function that will free the working structure allocated + by pfnInit. May be nullptr. + @param pfnProcess Processing function called to compute pixel values. Must + not be nullptr. + @param papszOptions Unused currently. Must be nullptr. + @return CE_None in case of success, error otherwise. + @since 3.9 + */ +CPLErr GDALVRTRegisterProcessedDatasetFunc( + const char *pszFuncName, void *pUserData, const char *pszXMLMetadata, + GDALDataType eRequestedInputDT, const GDALDataType *paeSupportedInputDT, + size_t nSupportedInputDTSize, const int *panSupportedInputBandCount, + size_t nSupportedInputBandCountSize, + GDALVRTProcessedDatasetFuncInit pfnInit, + GDALVRTProcessedDatasetFuncFree pfnFree, + GDALVRTProcessedDatasetFuncProcess pfnProcess, + CPL_UNUSED CSLConstList papszOptions) +{ + if (pszFuncName == nullptr || pszFuncName[0] == '\0') + { + CPLError(CE_Failure, CPLE_AppDefined, + "pszFuncName should be non-empty"); + return CE_Failure; + } + + auto &oMap = GetGlobalMapProcessedDatasetFunc(); + if (oMap.find(pszFuncName) != oMap.end()) + { + CPLError(CE_Failure, CPLE_AppDefined, "%s already registered", + pszFuncName); + return CE_Failure; + } + + if (!pfnProcess) + { + CPLError(CE_Failure, CPLE_AppDefined, "pfnProcess should not be null"); + return CE_Failure; + } + + VRTProcessedDatasetFunc oFunc; + oFunc.osFuncName = pszFuncName; + oFunc.pUserData = pUserData; + if (pszXMLMetadata) + { + oFunc.bMetadataSpecified = true; + auto psTree = CPLXMLTreeCloser(CPLParseXMLString(pszXMLMetadata)); + if (!psTree) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Cannot parse pszXMLMetadata=%s for %s", pszXMLMetadata, + pszFuncName); + return CE_Failure; + } + const CPLXMLNode *psRoot = CPLGetXMLNode( + psTree.get(), "=ProcessedDatasetFunctionArgumentsList"); + if (!psRoot) + { + CPLError(CE_Failure, CPLE_AppDefined, + "No root ProcessedDatasetFunctionArgumentsList element in " + "pszXMLMetadata=%s for %s", + pszXMLMetadata, pszFuncName); + return CE_Failure; + } + for (const CPLXMLNode *psIter = psRoot->psChild; psIter; + psIter = psIter->psNext) + { + if (psIter->eType == CXT_Element && + strcmp(psIter->pszValue, "Argument") == 0) + { + const char *pszName = CPLGetXMLValue(psIter, "name", nullptr); + if (!pszName) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Missing Argument.name attribute in " + "pszXMLMetadata=%s for %s", + pszXMLMetadata, pszFuncName); + return CE_Failure; + } + const char *pszType = CPLGetXMLValue(psIter, "type", nullptr); + if (!pszType) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Missing Argument.type attribute in " + "pszXMLMetadata=%s for %s", + pszXMLMetadata, pszFuncName); + return CE_Failure; + } + if (strcmp(pszType, "constant") == 0) + { + const char *pszValue = + CPLGetXMLValue(psIter, "value", nullptr); + if (!pszValue) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Missing Argument.value attribute in " + "pszXMLMetadata=%s for %s", + pszXMLMetadata, pszFuncName); + return CE_Failure; + } + oFunc.oMapConstantArguments[CPLString(pszName).tolower()] = + pszValue; + } + else if (strcmp(pszType, "builtin") == 0) + { + if (EQUAL(pszName, "nodata") || + EQUAL(pszName, "offset_{band}") || + EQUAL(pszName, "scale_{band}")) + { + oFunc.oSetBuiltinArguments.insert( + CPLString(pszName).tolower()); + } + else + { + CPLError(CE_Failure, CPLE_NotSupported, + "Unsupported builtin parameter name %s in " + "pszXMLMetadata=%s for %s. Only nodata, " + "offset_{band} and scale_{band} are supported", + pszName, pszXMLMetadata, pszFuncName); + return CE_Failure; + } + } + else if (strcmp(pszType, "boolean") == 0 || + strcmp(pszType, "string") == 0 || + strcmp(pszType, "integer") == 0 || + strcmp(pszType, "double") == 0 || + strcmp(pszType, "double_list") == 0) + { + VRTProcessedDatasetFunc::OtherArgument otherArgument; + otherArgument.bRequired = CPLTestBool( + CPLGetXMLValue(psIter, "required", "false")); + otherArgument.osType = pszType; + oFunc.oOtherArguments[CPLString(pszName).tolower()] = + std::move(otherArgument); + } + else + { + CPLError(CE_Failure, CPLE_NotSupported, + "Unsupported type for parameter %s in " + "pszXMLMetadata=%s for %s. Only boolean, string, " + "integer, double and double_list are supported", + pszName, pszXMLMetadata, pszFuncName); + return CE_Failure; + } + } + } + } + oFunc.eRequestedInputDT = eRequestedInputDT; + if (nSupportedInputDTSize) + { + oFunc.aeSupportedInputDT.insert( + oFunc.aeSupportedInputDT.end(), paeSupportedInputDT, + paeSupportedInputDT + nSupportedInputDTSize); + } + if (nSupportedInputBandCountSize) + { + oFunc.anSupportedInputBandCount.insert( + oFunc.anSupportedInputBandCount.end(), panSupportedInputBandCount, + panSupportedInputBandCount + nSupportedInputBandCountSize); + } + oFunc.pfnInit = pfnInit; + oFunc.pfnFree = pfnFree; + oFunc.pfnProcess = pfnProcess; + + oMap[pszFuncName] = std::move(oFunc); + + return CE_None; +} diff --git a/frmts/vrt/vrtprocesseddatasetfunctions.cpp b/frmts/vrt/vrtprocesseddatasetfunctions.cpp new file mode 100644 index 000000000000..7c41b5b31eec --- /dev/null +++ b/frmts/vrt/vrtprocesseddatasetfunctions.cpp @@ -0,0 +1,1673 @@ +/****************************************************************************** + * + * Project: Virtual GDAL Datasets + * Purpose: Implementation of VRTProcessedDataset processing functions + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2024, Even Rouault + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +#include "cpl_minixml.h" +#include "cpl_string.h" +#include "vrtdataset.h" + +#include +#include +#include +#include +#include + +/************************************************************************/ +/* GetDstValue() */ +/************************************************************************/ + +/** Return a destination value given an initial value, the destination no data + * value and its replacement value + */ +static inline double GetDstValue(double dfVal, double dfDstNoData, + double dfReplacementDstNodata, + GDALDataType eIntendedDstDT, + bool bDstIntendedDTIsInteger) +{ + if (bDstIntendedDTIsInteger && std::round(dfVal) == dfDstNoData) + { + return dfReplacementDstNodata; + } + else if (eIntendedDstDT == GDT_Float32 && + static_cast(dfVal) == static_cast(dfDstNoData)) + { + return dfReplacementDstNodata; + } + else if (eIntendedDstDT == GDT_Float64 && dfVal == dfDstNoData) + { + return dfReplacementDstNodata; + } + else + { + return dfVal; + } +} + +/************************************************************************/ +/* GetReplacementNoDataValue() */ +/************************************************************************/ + +/** Return a value that is close to the dfDstNoData but not identical to it. + * + * To be used to avoid the valid result of a computation to be the nodata value. + */ +static double GetReplacementNoDataValue(double dfDstNoData, + GDALDataType eIntendedDstDT) +{ + const bool bDstIntendedDTIsInteger = GDALDataTypeIsInteger(eIntendedDstDT); + double dfReplacementDstNodata; + if (bDstIntendedDTIsInteger) + { + if ((eIntendedDstDT == GDT_Byte || eIntendedDstDT == GDT_UInt16 || + eIntendedDstDT == GDT_UInt32) && + dfDstNoData == 0) + { + dfReplacementDstNodata = 1; + } + else if (eIntendedDstDT == GDT_Int8 && + dfDstNoData == std::numeric_limits::min()) + { + dfReplacementDstNodata = std::numeric_limits::min() + 1; + } + else if (eIntendedDstDT == GDT_Int16 && + dfDstNoData == std::numeric_limits::min()) + { + dfReplacementDstNodata = std::numeric_limits::min() + 1; + } + else if (eIntendedDstDT == GDT_Int32 && + dfDstNoData == std::numeric_limits::min()) + { + dfReplacementDstNodata = std::numeric_limits::min() + 1; + } + else + { + dfReplacementDstNodata = dfDstNoData - 1; + } + } + else + { + if (dfDstNoData == 0) + { + dfReplacementDstNodata = std::numeric_limits::min(); + } + else + { + // Relative difference between 2 consecutive float32 + dfReplacementDstNodata = dfDstNoData * (1.0 + 1.0 / (1 << 23)); + } + } + return dfReplacementDstNodata; +} + +/************************************************************************/ +/* BandAffineCombinationData */ +/************************************************************************/ + +namespace +{ +/** Working structure for 'BandAffineCombination' builtin function. */ +struct BandAffineCombinationData +{ + static constexpr const char *const EXPECTED_SIGNATURE = + "BandAffineCombination"; + //! Signature (to make sure callback functions are called with the right argument) + const std::string m_osSignature = EXPECTED_SIGNATURE; + + /** Replacement nodata value */ + std::vector m_adfReplacementDstNodata{}; + + /** Intended destination data type. */ + GDALDataType m_eIntendedDstDT = GDT_Float64; + + /** Affine transformation coefficients. + * m_aadfCoefficients[i][0] is the constant term for the i(th) dst band + * m_aadfCoefficients[i][j] is the weight of the j(th) src band for the + * i(th) dst vand. + * Said otherwise dst[i] = m_aadfCoefficients[i][0] + + * sum(m_aadfCoefficients[i][j + 1] * src[j] for j in 0...nSrcBands-1) + */ + std::vector> m_aadfCoefficients{}; + + //! Minimum clamping value. + double m_dfClampMin = std::numeric_limits::quiet_NaN(); + + //! Maximum clamping value. + double m_dfClampMax = std::numeric_limits::quiet_NaN(); +}; +} // namespace + +/************************************************************************/ +/* SetOutputValuesForInNoDataAndOutNoData() */ +/************************************************************************/ + +static std::vector SetOutputValuesForInNoDataAndOutNoData( + int nInBands, double *padfInNoData, int *pnOutBands, + double **ppadfOutNoData, bool bSrcNodataSpecified, double dfSrcNoData, + bool bDstNodataSpecified, double dfDstNoData, bool bIsFinalStep) +{ + if (bSrcNodataSpecified) + { + std::vector adfNoData(nInBands, dfSrcNoData); + memcpy(padfInNoData, adfNoData.data(), + adfNoData.size() * sizeof(double)); + } + + std::vector adfDstNoData; + if (bDstNodataSpecified) + { + adfDstNoData.resize(*pnOutBands, dfDstNoData); + } + else if (bIsFinalStep) + { + adfDstNoData = + std::vector(*ppadfOutNoData, *ppadfOutNoData + *pnOutBands); + } + else + { + adfDstNoData = + std::vector(padfInNoData, padfInNoData + nInBands); + adfDstNoData.resize(*pnOutBands, *padfInNoData); + } + + if (*ppadfOutNoData == nullptr) + { + *ppadfOutNoData = + static_cast(CPLMalloc(*pnOutBands * sizeof(double))); + } + memcpy(*ppadfOutNoData, adfDstNoData.data(), *pnOutBands * sizeof(double)); + + return adfDstNoData; +} + +/************************************************************************/ +/* BandAffineCombinationInit() */ +/************************************************************************/ + +/** Init function for 'BandAffineCombination' builtin function. */ +static CPLErr BandAffineCombinationInit( + const char * /*pszFuncName*/, void * /*pUserData*/, + CSLConstList papszFunctionArgs, int nInBands, GDALDataType eInDT, + double *padfInNoData, int *pnOutBands, GDALDataType *peOutDT, + double **ppadfOutNoData, const char * /* pszVRTPath */, + VRTPDWorkingDataPtr *ppWorkingData) +{ + CPLAssert(eInDT == GDT_Float64); + + *peOutDT = eInDT; + *ppWorkingData = nullptr; + + auto data = std::make_unique(); + + std::map> oMapCoefficients{}; + double dfSrcNoData = std::numeric_limits::quiet_NaN(); + bool bSrcNodataSpecified = false; + double dfDstNoData = std::numeric_limits::quiet_NaN(); + bool bDstNodataSpecified = false; + double dfReplacementDstNodata = std::numeric_limits::quiet_NaN(); + bool bReplacementDstNodataSpecified = false; + for (CSLConstList papszIter = papszFunctionArgs; papszIter && *papszIter; + ++papszIter) + { + char *pszKey = nullptr; + const char *pszValue = CPLParseNameValue(*papszIter, &pszKey); + if (pszKey && pszValue) + { + if (EQUAL(pszKey, "src_nodata")) + { + bSrcNodataSpecified = true; + dfSrcNoData = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "dst_nodata")) + { + bDstNodataSpecified = true; + dfDstNoData = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "replacement_nodata")) + { + bReplacementDstNodataSpecified = true; + dfReplacementDstNodata = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "dst_intended_datatype")) + { + for (GDALDataType eDT = GDT_Byte; eDT < GDT_TypeCount; + eDT = static_cast(eDT + 1)) + { + if (EQUAL(GDALGetDataTypeName(eDT), pszValue)) + { + data->m_eIntendedDstDT = eDT; + break; + } + } + } + else if (STARTS_WITH_CI(pszKey, "coefficients_")) + { + const int nTargetBand = atoi(pszKey + strlen("coefficients_")); + if (nTargetBand <= 0 || nTargetBand > 65536) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band in argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + const CPLStringList aosTokens( + CSLTokenizeString2(pszValue, ",", 0)); + if (aosTokens.size() != 1 + nInBands) + { + CPLError( + CE_Failure, CPLE_AppDefined, + "Argument %s has %d values, whereas %d are expected", + pszKey, aosTokens.size(), 1 + nInBands); + CPLFree(pszKey); + return CE_Failure; + } + std::vector adfValues; + for (int i = 0; i < aosTokens.size(); ++i) + { + adfValues.push_back(CPLAtof(aosTokens[i])); + } + oMapCoefficients[nTargetBand - 1] = std::move(adfValues); + } + else if (EQUAL(pszKey, "min")) + { + data->m_dfClampMin = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "max")) + { + data->m_dfClampMax = CPLAtof(pszValue); + } + else + { + CPLError(CE_Warning, CPLE_AppDefined, + "Unrecognized argument name %s. Ignored", pszKey); + } + } + CPLFree(pszKey); + } + + const bool bIsFinalStep = *pnOutBands != 0; + if (bIsFinalStep) + { + if (*pnOutBands != static_cast(oMapCoefficients.size())) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Final step expect %d bands, but only %d coefficient_XX " + "are provided", + *pnOutBands, static_cast(oMapCoefficients.size())); + return CE_Failure; + } + } + else + { + *pnOutBands = static_cast(oMapCoefficients.size()); + } + + const std::vector adfDstNoData = + SetOutputValuesForInNoDataAndOutNoData( + nInBands, padfInNoData, pnOutBands, ppadfOutNoData, + bSrcNodataSpecified, dfSrcNoData, bDstNodataSpecified, dfDstNoData, + bIsFinalStep); + + if (bReplacementDstNodataSpecified) + { + data->m_adfReplacementDstNodata.resize(*pnOutBands, + dfReplacementDstNodata); + } + else + { + for (double dfVal : adfDstNoData) + { + data->m_adfReplacementDstNodata.emplace_back( + GetReplacementNoDataValue(dfVal, data->m_eIntendedDstDT)); + } + } + + // Check we have a set of coefficient for all output bands and + // convert the map to a vector + for (auto &oIter : oMapCoefficients) + { + const int iExpected = static_cast(data->m_aadfCoefficients.size()); + if (oIter.first != iExpected) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Argument coefficients_%d is missing", iExpected + 1); + return CE_Failure; + } + data->m_aadfCoefficients.emplace_back(std::move(oIter.second)); + } + *ppWorkingData = data.release(); + return CE_None; +} + +/************************************************************************/ +/* BandAffineCombinationFree() */ +/************************************************************************/ + +/** Free function for 'BandAffineCombination' builtin function. */ +static void BandAffineCombinationFree(const char * /*pszFuncName*/, + void * /*pUserData*/, + VRTPDWorkingDataPtr pWorkingData) +{ + BandAffineCombinationData *data = + static_cast(pWorkingData); + CPLAssert(data->m_osSignature == + BandAffineCombinationData::EXPECTED_SIGNATURE); + CPL_IGNORE_RET_VAL(data->m_osSignature); + delete data; +} + +/************************************************************************/ +/* BandAffineCombinationProcess() */ +/************************************************************************/ + +/** Processing function for 'BandAffineCombination' builtin function. */ +static CPLErr BandAffineCombinationProcess( + const char * /*pszFuncName*/, void * /*pUserData*/, + VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs*/, + int nBufXSize, int nBufYSize, const void *pInBuffer, size_t nInBufferSize, + GDALDataType eInDT, int nInBands, const double *CPL_RESTRICT padfInNoData, + void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands, + const double *CPL_RESTRICT padfOutNoData, double /*dfSrcXOff*/, + double /*dfSrcYOff*/, double /*dfSrcXSize*/, double /*dfSrcYSize*/, + const double /*adfSrcGT*/[], const char * /* pszVRTPath */, + CSLConstList /*papszExtra*/) +{ + const size_t nElts = static_cast(nBufXSize) * nBufYSize; + + CPL_IGNORE_RET_VAL(eInDT); + CPLAssert(eInDT == GDT_Float64); + CPL_IGNORE_RET_VAL(eOutDT); + CPLAssert(eOutDT == GDT_Float64); + CPL_IGNORE_RET_VAL(nInBufferSize); + CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double)); + CPL_IGNORE_RET_VAL(nOutBufferSize); + CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double)); + + const BandAffineCombinationData *data = + static_cast(pWorkingData); + CPLAssert(data->m_osSignature == + BandAffineCombinationData::EXPECTED_SIGNATURE); + const double *CPL_RESTRICT padfSrc = static_cast(pInBuffer); + double *CPL_RESTRICT padfDst = static_cast(pOutBuffer); + const bool bDstIntendedDTIsInteger = + GDALDataTypeIsInteger(data->m_eIntendedDstDT); + const double dfClampMin = data->m_dfClampMin; + const double dfClampMax = data->m_dfClampMax; + for (size_t i = 0; i < nElts; ++i) + { + for (int iDst = 0; iDst < nOutBands; ++iDst) + { + const auto &adfCoefficients = data->m_aadfCoefficients[iDst]; + double dfVal = adfCoefficients[0]; + bool bSetNoData = false; + for (int iSrc = 0; iSrc < nInBands; ++iSrc) + { + // written this way to work with a NaN value + if (!(padfSrc[iSrc] != padfInNoData[iSrc])) + { + bSetNoData = true; + break; + } + dfVal += adfCoefficients[iSrc + 1] * padfSrc[iSrc]; + } + if (bSetNoData) + { + *padfDst = padfOutNoData[iDst]; + } + else + { + double dfDstVal = GetDstValue( + dfVal, padfOutNoData[iDst], + data->m_adfReplacementDstNodata[iDst], + data->m_eIntendedDstDT, bDstIntendedDTIsInteger); + if (dfDstVal < dfClampMin) + dfDstVal = dfClampMin; + if (dfDstVal > dfClampMax) + dfDstVal = dfClampMax; + *padfDst = dfDstVal; + } + ++padfDst; + } + padfSrc += nInBands; + } + + return CE_None; +} + +/************************************************************************/ +/* LUTData */ +/************************************************************************/ + +namespace +{ +/** Working structure for 'LUT' builtin function. */ +struct LUTData +{ + static constexpr const char *const EXPECTED_SIGNATURE = "LUT"; + //! Signature (to make sure callback functions are called with the right argument) + const std::string m_osSignature = EXPECTED_SIGNATURE; + + //! m_aadfLUTInputs[i][j] is the j(th) input value for that LUT of band i. + std::vector> m_aadfLUTInputs{}; + + //! m_aadfLUTOutputs[i][j] is the j(th) output value for that LUT of band i. + std::vector> m_aadfLUTOutputs{}; + + /************************************************************************/ + /* LookupValue() */ + /************************************************************************/ + + double LookupValue(int iBand, double dfInput) const + { + const auto &adfInput = m_aadfLUTInputs[iBand]; + const auto &afdOutput = m_aadfLUTOutputs[iBand]; + + // Find the index of the first element in the LUT input array that + // is not smaller than the input value. + int i = static_cast( + std::lower_bound(adfInput.data(), adfInput.data() + adfInput.size(), + dfInput) - + adfInput.data()); + + if (i == 0) + return afdOutput[0]; + + // If the index is beyond the end of the LUT input array, the input + // value is larger than all the values in the array. + if (i == static_cast(adfInput.size())) + return afdOutput.back(); + + if (adfInput[i] == dfInput) + return afdOutput[i]; + + // Otherwise, interpolate. + return afdOutput[i - 1] + (dfInput - adfInput[i - 1]) * + ((afdOutput[i] - afdOutput[i - 1]) / + (adfInput[i] - adfInput[i - 1])); + } +}; +} // namespace + +/************************************************************************/ +/* LUTInit() */ +/************************************************************************/ + +/** Init function for 'LUT' builtin function. */ +static CPLErr LUTInit(const char * /*pszFuncName*/, void * /*pUserData*/, + CSLConstList papszFunctionArgs, int nInBands, + GDALDataType eInDT, double *padfInNoData, int *pnOutBands, + GDALDataType *peOutDT, double **ppadfOutNoData, + const char * /* pszVRTPath */, + VRTPDWorkingDataPtr *ppWorkingData) +{ + CPLAssert(eInDT == GDT_Float64); + + const bool bIsFinalStep = *pnOutBands != 0; + *peOutDT = eInDT; + *ppWorkingData = nullptr; + + if (!bIsFinalStep) + { + *pnOutBands = nInBands; + } + + auto data = std::make_unique(); + + double dfSrcNoData = std::numeric_limits::quiet_NaN(); + bool bSrcNodataSpecified = false; + double dfDstNoData = std::numeric_limits::quiet_NaN(); + bool bDstNodataSpecified = false; + + std::map, std::vector>> oMap{}; + for (CSLConstList papszIter = papszFunctionArgs; papszIter && *papszIter; + ++papszIter) + { + char *pszKey = nullptr; + const char *pszValue = CPLParseNameValue(*papszIter, &pszKey); + if (pszKey && pszValue) + { + if (EQUAL(pszKey, "src_nodata")) + { + bSrcNodataSpecified = true; + dfSrcNoData = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "dst_nodata")) + { + bDstNodataSpecified = true; + dfDstNoData = CPLAtof(pszValue); + } + else if (STARTS_WITH_CI(pszKey, "lut_")) + { + const int nBand = atoi(pszKey + strlen("lut_")); + if (nBand <= 0 || nBand > nInBands) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band in argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + const CPLStringList aosTokens( + CSLTokenizeString2(pszValue, ",", 0)); + std::vector adfInputValues; + std::vector adfOutputValues; + for (int i = 0; i < aosTokens.size(); ++i) + { + const CPLStringList aosTokens2( + CSLTokenizeString2(aosTokens[i], ":", 0)); + if (aosTokens2.size() != 2) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid value for argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + adfInputValues.push_back(CPLAtof(aosTokens2[0])); + adfOutputValues.push_back(CPLAtof(aosTokens2[1])); + } + oMap[nBand - 1] = std::pair(std::move(adfInputValues), + std::move(adfOutputValues)); + } + else + { + CPLError(CE_Warning, CPLE_AppDefined, + "Unrecognized argument name %s. Ignored", pszKey); + } + } + CPLFree(pszKey); + } + + SetOutputValuesForInNoDataAndOutNoData( + nInBands, padfInNoData, pnOutBands, ppadfOutNoData, bSrcNodataSpecified, + dfSrcNoData, bDstNodataSpecified, dfDstNoData, bIsFinalStep); + + int iExpected = 0; + // Check we have values for all bands and convert to vector + for (auto &oIter : oMap) + { + if (oIter.first != iExpected) + { + CPLError(CE_Failure, CPLE_AppDefined, "Argument lut_%d is missing", + iExpected + 1); + return CE_Failure; + } + ++iExpected; + data->m_aadfLUTInputs.emplace_back(std::move(oIter.second.first)); + data->m_aadfLUTOutputs.emplace_back(std::move(oIter.second.second)); + } + + if (static_cast(oMap.size()) < *pnOutBands) + { + CPLError(CE_Failure, CPLE_AppDefined, "Missing lut_XX element(s)"); + return CE_Failure; + } + + *ppWorkingData = data.release(); + return CE_None; +} + +/************************************************************************/ +/* LUTFree() */ +/************************************************************************/ + +/** Free function for 'LUT' builtin function. */ +static void LUTFree(const char * /*pszFuncName*/, void * /*pUserData*/, + VRTPDWorkingDataPtr pWorkingData) +{ + LUTData *data = static_cast(pWorkingData); + CPLAssert(data->m_osSignature == LUTData::EXPECTED_SIGNATURE); + CPL_IGNORE_RET_VAL(data->m_osSignature); + delete data; +} + +/************************************************************************/ +/* LUTProcess() */ +/************************************************************************/ + +/** Processing function for 'LUT' builtin function. */ +static CPLErr +LUTProcess(const char * /*pszFuncName*/, void * /*pUserData*/, + VRTPDWorkingDataPtr pWorkingData, + CSLConstList /* papszFunctionArgs*/, int nBufXSize, int nBufYSize, + const void *pInBuffer, size_t nInBufferSize, GDALDataType eInDT, + int nInBands, const double *CPL_RESTRICT padfInNoData, + void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, + int nOutBands, const double *CPL_RESTRICT padfOutNoData, + double /*dfSrcXOff*/, double /*dfSrcYOff*/, double /*dfSrcXSize*/, + double /*dfSrcYSize*/, const double /*adfSrcGT*/[], + const char * /* pszVRTPath */, CSLConstList /*papszExtra*/) +{ + const size_t nElts = static_cast(nBufXSize) * nBufYSize; + + CPL_IGNORE_RET_VAL(eInDT); + CPLAssert(eInDT == GDT_Float64); + CPL_IGNORE_RET_VAL(eOutDT); + CPLAssert(eOutDT == GDT_Float64); + CPL_IGNORE_RET_VAL(nInBufferSize); + CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double)); + CPL_IGNORE_RET_VAL(nOutBufferSize); + CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double)); + CPLAssert(nInBands == nOutBands); + CPL_IGNORE_RET_VAL(nOutBands); + + const LUTData *data = static_cast(pWorkingData); + CPLAssert(data->m_osSignature == LUTData::EXPECTED_SIGNATURE); + const double *CPL_RESTRICT padfSrc = static_cast(pInBuffer); + double *CPL_RESTRICT padfDst = static_cast(pOutBuffer); + for (size_t i = 0; i < nElts; ++i) + { + for (int iBand = 0; iBand < nInBands; ++iBand) + { + // written this way to work with a NaN value + if (!(*padfSrc != padfInNoData[iBand])) + *padfDst = padfOutNoData[iBand]; + else + *padfDst = data->LookupValue(iBand, *padfSrc); + ++padfSrc; + ++padfDst; + } + } + + return CE_None; +} + +/************************************************************************/ +/* DehazingData */ +/************************************************************************/ + +namespace +{ +/** Working structure for 'Dehazing' builtin function. */ +struct DehazingData +{ + static constexpr const char *const EXPECTED_SIGNATURE = "Dehazing"; + //! Signature (to make sure callback functions are called with the right argument) + const std::string m_osSignature = EXPECTED_SIGNATURE; + + //! Nodata value for gain dataset(s) + double m_dfGainNodata = std::numeric_limits::quiet_NaN(); + + //! Nodata value for offset dataset(s) + double m_dfOffsetNodata = std::numeric_limits::quiet_NaN(); + + //! Minimum clamping value. + double m_dfClampMin = std::numeric_limits::quiet_NaN(); + + //! Maximum clamping value. + double m_dfClampMax = std::numeric_limits::quiet_NaN(); + + //! Map from gain/offset dataset name to datasets + std::map> m_oDatasetMap{}; + + //! Vector of size nInBands that point to the raster band from which to read gains. + std::vector m_oGainBands{}; + + //! Vector of size nInBands that point to the raster band from which to read offsets. + std::vector m_oOffsetBands{}; + + //! Working buffer that contain gain values. + std::vector m_abyGainBuffer{}; + + //! Working buffer that contain offset values. + std::vector m_abyOffsetBuffer{}; +}; +} // namespace + +/************************************************************************/ +/* CheckAllBands() */ +/************************************************************************/ + +/** Return true if the key of oMap is the sequence of all integers between + * 0 and nExpectedBandCount-1. + */ +template +static bool CheckAllBands(const std::map &oMap, int nExpectedBandCount) +{ + int iExpected = 0; + for (const auto &kv : oMap) + { + if (kv.first != iExpected) + return false; + ++iExpected; + } + return iExpected == nExpectedBandCount; +} + +/************************************************************************/ +/* DehazingInit() */ +/************************************************************************/ + +/** Init function for 'Dehazing' builtin function. */ +static CPLErr DehazingInit(const char * /*pszFuncName*/, void * /*pUserData*/, + CSLConstList papszFunctionArgs, int nInBands, + GDALDataType eInDT, double *padfInNoData, + int *pnOutBands, GDALDataType *peOutDT, + double **ppadfOutNoData, const char *pszVRTPath, + VRTPDWorkingDataPtr *ppWorkingData) +{ + CPLAssert(eInDT == GDT_Float64); + + const bool bIsFinalStep = *pnOutBands != 0; + *peOutDT = eInDT; + *ppWorkingData = nullptr; + + if (!bIsFinalStep) + { + *pnOutBands = nInBands; + } + + auto data = std::make_unique(); + + bool bNodataSpecified = false; + double dfNoData = std::numeric_limits::quiet_NaN(); + + bool bGainNodataSpecified = false; + bool bOffsetNodataSpecified = false; + + std::map oGainDatasetNameMap; + std::map oGainDatasetBandMap; + + std::map oOffsetDatasetNameMap; + std::map oOffsetDatasetBandMap; + + bool bRelativeFilenames = false; + + for (CSLConstList papszIter = papszFunctionArgs; papszIter && *papszIter; + ++papszIter) + { + char *pszKey = nullptr; + const char *pszValue = CPLParseNameValue(*papszIter, &pszKey); + if (pszKey && pszValue) + { + if (EQUAL(pszKey, "relative_filenames")) + { + bRelativeFilenames = CPLTestBool(pszValue); + } + else if (EQUAL(pszKey, "nodata")) + { + bNodataSpecified = true; + dfNoData = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "gain_nodata")) + { + bGainNodataSpecified = true; + data->m_dfGainNodata = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "offset_nodata")) + { + bOffsetNodataSpecified = true; + data->m_dfOffsetNodata = CPLAtof(pszValue); + } + else if (STARTS_WITH_CI(pszKey, "gain_dataset_filename_")) + { + const int nBand = + atoi(pszKey + strlen("gain_dataset_filename_")); + if (nBand <= 0 || nBand > nInBands) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band in argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + oGainDatasetNameMap[nBand - 1] = pszValue; + } + else if (STARTS_WITH_CI(pszKey, "gain_dataset_band_")) + { + const int nBand = atoi(pszKey + strlen("gain_dataset_band_")); + if (nBand <= 0 || nBand > nInBands) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band in argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + oGainDatasetBandMap[nBand - 1] = atoi(pszValue); + } + else if (STARTS_WITH_CI(pszKey, "offset_dataset_filename_")) + { + const int nBand = + atoi(pszKey + strlen("offset_dataset_filename_")); + if (nBand <= 0 || nBand > nInBands) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band in argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + oOffsetDatasetNameMap[nBand - 1] = pszValue; + } + else if (STARTS_WITH_CI(pszKey, "offset_dataset_band_")) + { + const int nBand = atoi(pszKey + strlen("offset_dataset_band_")); + if (nBand <= 0 || nBand > nInBands) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band in argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + oOffsetDatasetBandMap[nBand - 1] = atoi(pszValue); + } + else if (EQUAL(pszKey, "min")) + { + data->m_dfClampMin = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "max")) + { + data->m_dfClampMax = CPLAtof(pszValue); + } + else + { + CPLError(CE_Warning, CPLE_AppDefined, + "Unrecognized argument name %s. Ignored", pszKey); + } + } + CPLFree(pszKey); + } + + if (!CheckAllBands(oGainDatasetNameMap, nInBands)) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Missing gain_dataset_filename_XX element(s)"); + return CE_Failure; + } + if (!CheckAllBands(oGainDatasetBandMap, nInBands)) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Missing gain_dataset_band_XX element(s)"); + return CE_Failure; + } + if (!CheckAllBands(oOffsetDatasetNameMap, nInBands)) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Missing offset_dataset_filename_XX element(s)"); + return CE_Failure; + } + if (!CheckAllBands(oOffsetDatasetBandMap, nInBands)) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Missing offset_dataset_band_XX element(s)"); + return CE_Failure; + } + + data->m_oGainBands.resize(nInBands); + data->m_oOffsetBands.resize(nInBands); + + constexpr int IDX_GAIN = 0; + constexpr int IDX_OFFSET = 1; + for (int i : {IDX_GAIN, IDX_OFFSET}) + { + const auto &oMapNames = + (i == IDX_GAIN) ? oGainDatasetNameMap : oOffsetDatasetNameMap; + const auto &oMapBands = + (i == IDX_GAIN) ? oGainDatasetBandMap : oOffsetDatasetBandMap; + for (const auto &kv : oMapNames) + { + const int nInBandIdx = kv.first; + const auto osFilename = VRTDataset::BuildSourceFilename( + kv.second.c_str(), pszVRTPath, bRelativeFilenames); + auto oIter = data->m_oDatasetMap.find(osFilename); + if (oIter == data->m_oDatasetMap.end()) + { + auto poDS = std::unique_ptr(GDALDataset::Open( + osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, + nullptr, nullptr, nullptr)); + if (!poDS) + return CE_Failure; + double adfAuxGT[6]; + if (poDS->GetGeoTransform(adfAuxGT) != CE_None) + { + CPLError(CE_Failure, CPLE_AppDefined, + "%s lacks a geotransform", osFilename.c_str()); + return CE_Failure; + } + oIter = data->m_oDatasetMap + .insert(std::pair(osFilename, std::move(poDS))) + .first; + } + auto poDS = oIter->second.get(); + const auto oIterBand = oMapBands.find(nInBandIdx); + CPLAssert(oIterBand != oMapBands.end()); + const int nAuxBand = oIterBand->second; + if (nAuxBand <= 0 || nAuxBand > poDS->GetRasterCount()) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band number (%d) for a %s dataset", nAuxBand, + (i == IDX_GAIN) ? "gain" : "offset"); + return CE_Failure; + } + auto poAuxBand = poDS->GetRasterBand(nAuxBand); + int bAuxBandHasNoData = false; + const double dfAuxNoData = + poAuxBand->GetNoDataValue(&bAuxBandHasNoData); + if (i == IDX_GAIN) + { + data->m_oGainBands[nInBandIdx] = poAuxBand; + if (!bGainNodataSpecified && bAuxBandHasNoData) + data->m_dfGainNodata = dfAuxNoData; + } + else + { + data->m_oOffsetBands[nInBandIdx] = poAuxBand; + if (!bOffsetNodataSpecified && bAuxBandHasNoData) + data->m_dfOffsetNodata = dfAuxNoData; + } + } + } + + SetOutputValuesForInNoDataAndOutNoData( + nInBands, padfInNoData, pnOutBands, ppadfOutNoData, bNodataSpecified, + dfNoData, bNodataSpecified, dfNoData, bIsFinalStep); + + *ppWorkingData = data.release(); + return CE_None; +} + +/************************************************************************/ +/* DehazingFree() */ +/************************************************************************/ + +/** Free function for 'Dehazing' builtin function. */ +static void DehazingFree(const char * /*pszFuncName*/, void * /*pUserData*/, + VRTPDWorkingDataPtr pWorkingData) +{ + DehazingData *data = static_cast(pWorkingData); + CPLAssert(data->m_osSignature == DehazingData::EXPECTED_SIGNATURE); + CPL_IGNORE_RET_VAL(data->m_osSignature); + delete data; +} + +/************************************************************************/ +/* LoadAuxData() */ +/************************************************************************/ + +// Load auxiliary corresponding offset, gain or trimming data. +static bool LoadAuxData(double dfULX, double dfULY, double dfLRX, double dfLRY, + size_t nElts, int nBufXSize, int nBufYSize, + const char *pszAuxType, GDALRasterBand *poAuxBand, + std::vector &abyBuffer) +{ + double adfAuxGT[6]; + double adfAuxInvGT[6]; + + // Compute pixel/line coordinates from the georeferenced extent + CPL_IGNORE_RET_VAL(poAuxBand->GetDataset()->GetGeoTransform( + adfAuxGT)); // return code already tested + CPL_IGNORE_RET_VAL(GDALInvGeoTransform(adfAuxGT, adfAuxInvGT)); + const double dfULPixel = + adfAuxInvGT[0] + adfAuxInvGT[1] * dfULX + adfAuxInvGT[2] * dfULY; + const double dfULLine = + adfAuxInvGT[3] + adfAuxInvGT[4] * dfULX + adfAuxInvGT[5] * dfULY; + const double dfLRPixel = + adfAuxInvGT[0] + adfAuxInvGT[1] * dfLRX + adfAuxInvGT[2] * dfLRY; + const double dfLRLine = + adfAuxInvGT[3] + adfAuxInvGT[4] * dfLRX + adfAuxInvGT[5] * dfLRY; + if (dfULPixel >= dfLRPixel || dfULLine >= dfLRLine) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Unexpected computed %s pixel/line", pszAuxType); + return false; + } + if (dfULPixel < -1 || dfLRPixel > poAuxBand->GetXSize() || dfULLine < -1 || + dfLRLine > poAuxBand->GetYSize()) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Unexpected computed %s pixel/line", pszAuxType); + return false; + } + + const int nAuxXOff = std::max(0, static_cast(std::round(dfULPixel))); + const int nAuxYOff = std::max(0, static_cast(std::round(dfULLine))); + const int nAuxX2Off = std::min(poAuxBand->GetXSize(), + static_cast(std::round(dfLRPixel))); + const int nAuxY2Off = + std::min(poAuxBand->GetYSize(), static_cast(std::round(dfLRLine))); + + try + { + abyBuffer.resize(nElts * sizeof(float)); + } + catch (const std::bad_alloc &) + { + CPLError(CE_Failure, CPLE_OutOfMemory, + "Out of memory allocating working buffer"); + return false; + } + GDALRasterIOExtraArg sExtraArg; + INIT_RASTERIO_EXTRA_ARG(sExtraArg); + sExtraArg.bFloatingPointWindowValidity = true; + CPL_IGNORE_RET_VAL(sExtraArg.eResampleAlg); + sExtraArg.eResampleAlg = GRIORA_Bilinear; + sExtraArg.dfXOff = std::max(0.0, dfULPixel); + sExtraArg.dfYOff = std::max(0.0, dfULLine); + sExtraArg.dfXSize = std::min(poAuxBand->GetXSize(), dfLRPixel) - + std::max(0.0, dfULPixel); + sExtraArg.dfYSize = std::min(poAuxBand->GetYSize(), dfLRLine) - + std::max(0.0, dfULLine); + return (poAuxBand->RasterIO( + GF_Read, nAuxXOff, nAuxYOff, std::max(1, nAuxX2Off - nAuxXOff), + std::max(1, nAuxY2Off - nAuxYOff), abyBuffer.data(), nBufXSize, + nBufYSize, GDT_Float32, 0, 0, &sExtraArg) == CE_None); +} + +/************************************************************************/ +/* DehazingProcess() */ +/************************************************************************/ + +/** Processing function for 'Dehazing' builtin function. */ +static CPLErr DehazingProcess( + const char * /*pszFuncName*/, void * /*pUserData*/, + VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs*/, + int nBufXSize, int nBufYSize, const void *pInBuffer, size_t nInBufferSize, + GDALDataType eInDT, int nInBands, const double *CPL_RESTRICT padfInNoData, + void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands, + const double *CPL_RESTRICT padfOutNoData, double dfSrcXOff, + double dfSrcYOff, double dfSrcXSize, double dfSrcYSize, + const double adfSrcGT[], const char * /* pszVRTPath */, + CSLConstList /*papszExtra*/) +{ + const size_t nElts = static_cast(nBufXSize) * nBufYSize; + + CPL_IGNORE_RET_VAL(eInDT); + CPLAssert(eInDT == GDT_Float64); + CPL_IGNORE_RET_VAL(eOutDT); + CPLAssert(eOutDT == GDT_Float64); + CPL_IGNORE_RET_VAL(nInBufferSize); + CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double)); + CPL_IGNORE_RET_VAL(nOutBufferSize); + CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double)); + CPLAssert(nInBands == nOutBands); + CPL_IGNORE_RET_VAL(nOutBands); + + DehazingData *data = static_cast(pWorkingData); + CPLAssert(data->m_osSignature == DehazingData::EXPECTED_SIGNATURE); + const double *CPL_RESTRICT padfSrc = static_cast(pInBuffer); + double *CPL_RESTRICT padfDst = static_cast(pOutBuffer); + + // Compute georeferenced extent of input region + const double dfULX = + adfSrcGT[0] + adfSrcGT[1] * dfSrcXOff + adfSrcGT[2] * dfSrcYOff; + const double dfULY = + adfSrcGT[3] + adfSrcGT[4] * dfSrcXOff + adfSrcGT[5] * dfSrcYOff; + const double dfLRX = adfSrcGT[0] + adfSrcGT[1] * (dfSrcXOff + dfSrcXSize) + + adfSrcGT[2] * (dfSrcYOff + dfSrcYSize); + const double dfLRY = adfSrcGT[3] + adfSrcGT[4] * (dfSrcXOff + dfSrcXSize) + + adfSrcGT[5] * (dfSrcYOff + dfSrcYSize); + + auto &abyOffsetBuffer = data->m_abyGainBuffer; + auto &abyGainBuffer = data->m_abyOffsetBuffer; + + for (int iBand = 0; iBand < nInBands; ++iBand) + { + if (!LoadAuxData(dfULX, dfULY, dfLRX, dfLRY, nElts, nBufXSize, + nBufYSize, "gain", data->m_oGainBands[iBand], + abyGainBuffer) || + !LoadAuxData(dfULX, dfULY, dfLRX, dfLRY, nElts, nBufXSize, + nBufYSize, "offset", data->m_oOffsetBands[iBand], + abyOffsetBuffer)) + { + return CE_Failure; + } + + const double *CPL_RESTRICT padfSrcThisBand = padfSrc + iBand; + double *CPL_RESTRICT padfDstThisBand = padfDst + iBand; + const float *pafGain = + reinterpret_cast(abyGainBuffer.data()); + const float *pafOffset = + reinterpret_cast(abyOffsetBuffer.data()); + const double dfSrcNodata = padfInNoData[iBand]; + const double dfDstNodata = padfOutNoData[iBand]; + const double dfGainNodata = data->m_dfGainNodata; + const double dfOffsetNodata = data->m_dfOffsetNodata; + const double dfClampMin = data->m_dfClampMin; + const double dfClampMax = data->m_dfClampMax; + for (size_t i = 0; i < nElts; ++i) + { + const double dfSrcVal = *padfSrcThisBand; + // written this way to work with a NaN value + if (!(dfSrcVal != dfSrcNodata)) + { + *padfDstThisBand = dfDstNodata; + } + else + { + const double dfGain = pafGain[i]; + const double dfOffset = pafOffset[i]; + if (!(dfGain != dfGainNodata) || !(dfOffset != dfOffsetNodata)) + { + *padfDstThisBand = dfDstNodata; + } + else + { + double dfDehazed = dfSrcVal * dfGain - dfOffset; + if (dfDehazed < dfClampMin) + dfDehazed = dfClampMin; + if (dfDehazed > dfClampMax) + dfDehazed = dfClampMax; + + *padfDstThisBand = dfDehazed; + } + } + padfSrcThisBand += nInBands; + padfDstThisBand += nInBands; + } + } + + return CE_None; +} + +/************************************************************************/ +/* TrimmingData */ +/************************************************************************/ + +namespace +{ +/** Working structure for 'Trimming' builtin function. */ +struct TrimmingData +{ + static constexpr const char *const EXPECTED_SIGNATURE = "Trimming"; + //! Signature (to make sure callback functions are called with the right argument) + const std::string m_osSignature = EXPECTED_SIGNATURE; + + //! Nodata value for trimming dataset + double m_dfTrimmingNodata = std::numeric_limits::quiet_NaN(); + + //! Maximum saturating RGB output value. + double m_dfTopRGB = 0; + + //! Maximum threshold beyond which we give up saturation + double m_dfToneCeil = 0; + + //! Margin to allow for dynamics in brighest areas (in [0,1] range) + double m_dfTopMargin = 0; + + //! Index (zero-based) of input/output red band. + int m_nRedBand = 1 - 1; + + //! Index (zero-based) of input/output green band. + int m_nGreenBand = 2 - 1; + + //! Index (zero-based) of input/output blue band. + int m_nBlueBand = 3 - 1; + + //! Trimming dataset + std::unique_ptr m_poTrimmingDS{}; + + //! Trimming raster band. + GDALRasterBand *m_poTrimmingBand = nullptr; + + //! Working buffer that contain trimming values. + std::vector m_abyTrimmingBuffer{}; +}; +} // namespace + +/************************************************************************/ +/* TrimmingInit() */ +/************************************************************************/ + +/** Init function for 'Trimming' builtin function. */ +static CPLErr TrimmingInit(const char * /*pszFuncName*/, void * /*pUserData*/, + CSLConstList papszFunctionArgs, int nInBands, + GDALDataType eInDT, double *padfInNoData, + int *pnOutBands, GDALDataType *peOutDT, + double **ppadfOutNoData, const char *pszVRTPath, + VRTPDWorkingDataPtr *ppWorkingData) +{ + CPLAssert(eInDT == GDT_Float64); + + const bool bIsFinalStep = *pnOutBands != 0; + *peOutDT = eInDT; + *ppWorkingData = nullptr; + + if (!bIsFinalStep) + { + *pnOutBands = nInBands; + } + + auto data = std::make_unique(); + + bool bNodataSpecified = false; + double dfNoData = std::numeric_limits::quiet_NaN(); + std::string osTrimmingFilename; + bool bTrimmingNodataSpecified = false; + bool bRelativeFilenames = false; + + for (CSLConstList papszIter = papszFunctionArgs; papszIter && *papszIter; + ++papszIter) + { + char *pszKey = nullptr; + const char *pszValue = CPLParseNameValue(*papszIter, &pszKey); + if (pszKey && pszValue) + { + if (EQUAL(pszKey, "relative_filenames")) + { + bRelativeFilenames = CPLTestBool(pszValue); + } + else if (EQUAL(pszKey, "nodata")) + { + bNodataSpecified = true; + dfNoData = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "trimming_nodata")) + { + bTrimmingNodataSpecified = true; + data->m_dfTrimmingNodata = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "trimming_dataset_filename")) + { + osTrimmingFilename = pszValue; + } + else if (EQUAL(pszKey, "red_band")) + { + const int nBand = atoi(pszValue) - 1; + if (nBand < 0 || nBand >= nInBands) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band in argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + data->m_nRedBand = nBand; + } + else if (EQUAL(pszKey, "green_band")) + { + const int nBand = atoi(pszValue) - 1; + if (nBand < 0 || nBand >= nInBands) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band in argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + data->m_nGreenBand = nBand; + } + else if (EQUAL(pszKey, "blue_band")) + { + const int nBand = atoi(pszValue) - 1; + if (nBand < 0 || nBand >= nInBands) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Invalid band in argument '%s'", pszKey); + CPLFree(pszKey); + return CE_Failure; + } + data->m_nBlueBand = nBand; + } + else if (EQUAL(pszKey, "top_rgb")) + { + data->m_dfTopRGB = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "tone_ceil")) + { + data->m_dfToneCeil = CPLAtof(pszValue); + } + else if (EQUAL(pszKey, "top_margin")) + { + data->m_dfTopMargin = CPLAtof(pszValue); + } + else + { + CPLError(CE_Warning, CPLE_AppDefined, + "Unrecognized argument name %s. Ignored", pszKey); + } + } + CPLFree(pszKey); + } + + if (data->m_nRedBand == data->m_nGreenBand || + data->m_nRedBand == data->m_nBlueBand || + data->m_nGreenBand == data->m_nBlueBand) + { + CPLError( + CE_Failure, CPLE_NotSupported, + "red_band, green_band and blue_band must have distinct values"); + return CE_Failure; + } + + const auto osFilename = VRTDataset::BuildSourceFilename( + osTrimmingFilename.c_str(), pszVRTPath, bRelativeFilenames); + data->m_poTrimmingDS.reset(GDALDataset::Open( + osFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr, + nullptr, nullptr)); + if (!data->m_poTrimmingDS) + return CE_Failure; + if (data->m_poTrimmingDS->GetRasterCount() != 1) + { + CPLError(CE_Failure, CPLE_NotSupported, + "Trimming dataset should have a single band"); + return CE_Failure; + } + data->m_poTrimmingBand = data->m_poTrimmingDS->GetRasterBand(1); + + double adfAuxGT[6]; + if (data->m_poTrimmingDS->GetGeoTransform(adfAuxGT) != CE_None) + { + CPLError(CE_Failure, CPLE_AppDefined, "%s lacks a geotransform", + osFilename.c_str()); + return CE_Failure; + } + int bAuxBandHasNoData = false; + const double dfAuxNoData = + data->m_poTrimmingBand->GetNoDataValue(&bAuxBandHasNoData); + if (!bTrimmingNodataSpecified && bAuxBandHasNoData) + data->m_dfTrimmingNodata = dfAuxNoData; + + SetOutputValuesForInNoDataAndOutNoData( + nInBands, padfInNoData, pnOutBands, ppadfOutNoData, bNodataSpecified, + dfNoData, bNodataSpecified, dfNoData, bIsFinalStep); + + *ppWorkingData = data.release(); + return CE_None; +} + +/************************************************************************/ +/* TrimmingFree() */ +/************************************************************************/ + +/** Free function for 'Trimming' builtin function. */ +static void TrimmingFree(const char * /*pszFuncName*/, void * /*pUserData*/, + VRTPDWorkingDataPtr pWorkingData) +{ + TrimmingData *data = static_cast(pWorkingData); + CPLAssert(data->m_osSignature == TrimmingData::EXPECTED_SIGNATURE); + CPL_IGNORE_RET_VAL(data->m_osSignature); + delete data; +} + +/************************************************************************/ +/* TrimmingProcess() */ +/************************************************************************/ + +/** Processing function for 'Trimming' builtin function. */ +static CPLErr TrimmingProcess( + const char * /*pszFuncName*/, void * /*pUserData*/, + VRTPDWorkingDataPtr pWorkingData, CSLConstList /* papszFunctionArgs*/, + int nBufXSize, int nBufYSize, const void *pInBuffer, size_t nInBufferSize, + GDALDataType eInDT, int nInBands, const double *CPL_RESTRICT padfInNoData, + void *pOutBuffer, size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands, + const double *CPL_RESTRICT padfOutNoData, double dfSrcXOff, + double dfSrcYOff, double dfSrcXSize, double dfSrcYSize, + const double adfSrcGT[], const char * /* pszVRTPath */, + CSLConstList /*papszExtra*/) +{ + const size_t nElts = static_cast(nBufXSize) * nBufYSize; + + CPL_IGNORE_RET_VAL(eInDT); + CPLAssert(eInDT == GDT_Float64); + CPL_IGNORE_RET_VAL(eOutDT); + CPLAssert(eOutDT == GDT_Float64); + CPL_IGNORE_RET_VAL(nInBufferSize); + CPLAssert(nInBufferSize == nElts * nInBands * sizeof(double)); + CPL_IGNORE_RET_VAL(nOutBufferSize); + CPLAssert(nOutBufferSize == nElts * nOutBands * sizeof(double)); + CPLAssert(nInBands == nOutBands); + CPL_IGNORE_RET_VAL(nOutBands); + + TrimmingData *data = static_cast(pWorkingData); + CPLAssert(data->m_osSignature == TrimmingData::EXPECTED_SIGNATURE); + const double *CPL_RESTRICT padfSrc = static_cast(pInBuffer); + double *CPL_RESTRICT padfDst = static_cast(pOutBuffer); + + // Compute georeferenced extent of input region + const double dfULX = + adfSrcGT[0] + adfSrcGT[1] * dfSrcXOff + adfSrcGT[2] * dfSrcYOff; + const double dfULY = + adfSrcGT[3] + adfSrcGT[4] * dfSrcXOff + adfSrcGT[5] * dfSrcYOff; + const double dfLRX = adfSrcGT[0] + adfSrcGT[1] * (dfSrcXOff + dfSrcXSize) + + adfSrcGT[2] * (dfSrcYOff + dfSrcYSize); + const double dfLRY = adfSrcGT[3] + adfSrcGT[4] * (dfSrcXOff + dfSrcXSize) + + adfSrcGT[5] * (dfSrcYOff + dfSrcYSize); + + if (!LoadAuxData(dfULX, dfULY, dfLRX, dfLRY, nElts, nBufXSize, nBufYSize, + "trimming", data->m_poTrimmingBand, + data->m_abyTrimmingBuffer)) + { + return CE_Failure; + } + + const float *pafTrimming = + reinterpret_cast(data->m_abyTrimmingBuffer.data()); + const int nRedBand = data->m_nRedBand; + const int nGreenBand = data->m_nGreenBand; + const int nBlueBand = data->m_nBlueBand; + const double dfTopMargin = data->m_dfTopMargin; + const double dfTopRGB = data->m_dfTopRGB; + const double dfToneCeil = data->m_dfToneCeil; +#if !defined(trimming_non_optimized_version) + const double dfInvToneCeil = 1.0 / dfToneCeil; +#endif + const bool bRGBBandsAreFirst = + std::max(std::max(nRedBand, nGreenBand), nBlueBand) <= 2; + const double dfNoDataTrimming = data->m_dfTrimmingNodata; + const double dfNoDataRed = padfInNoData[nRedBand]; + const double dfNoDataGreen = padfInNoData[nGreenBand]; + const double dfNoDataBlue = padfInNoData[nBlueBand]; + for (size_t i = 0; i < nElts; ++i) + { + // Extract local saturation value from trimming image + const double dfLocalMaxRGB = pafTrimming[i]; + const double dfReducedRGB = + std::min((1.0 - dfTopMargin) * dfTopRGB / dfLocalMaxRGB, 1.0); + + const double dfRed = padfSrc[nRedBand]; + const double dfGreen = padfSrc[nGreenBand]; + const double dfBlue = padfSrc[nBlueBand]; + bool bNoDataPixel = false; + if ((dfLocalMaxRGB != dfNoDataTrimming) && (dfRed != dfNoDataRed) && + (dfGreen != dfNoDataGreen) && (dfBlue != dfNoDataBlue)) + { + // RGB bands specific process + const double dfMaxRGB = std::max(std::max(dfRed, dfGreen), dfBlue); +#if !defined(trimming_non_optimized_version) + const double dfRedTimesToneRed = std::min(dfRed, dfToneCeil); + const double dfGreenTimesToneGreen = std::min(dfGreen, dfToneCeil); + const double dfBlueTimesToneBlue = std::min(dfBlue, dfToneCeil); + const double dfInvToneMaxRGB = + std::max(dfMaxRGB * dfInvToneCeil, 1.0); + const double dfReducedRGBTimesInvToneMaxRGB = + dfReducedRGB * dfInvToneMaxRGB; + padfDst[nRedBand] = std::min( + dfRedTimesToneRed * dfReducedRGBTimesInvToneMaxRGB, dfTopRGB); + padfDst[nGreenBand] = + std::min(dfGreenTimesToneGreen * dfReducedRGBTimesInvToneMaxRGB, + dfTopRGB); + padfDst[nBlueBand] = std::min( + dfBlueTimesToneBlue * dfReducedRGBTimesInvToneMaxRGB, dfTopRGB); +#else + // Original formulas. Slightly less optimized than the above ones. + const double dfToneMaxRGB = std::min(dfToneCeil / dfMaxRGB, 1.0); + const double dfToneRed = std::min(dfToneCeil / dfRed, 1.0); + const double dfToneGreen = std::min(dfToneCeil / dfGreen, 1.0); + const double dfToneBlue = std::min(dfToneCeil / dfBlue, 1.0); + padfDst[nRedBand] = std::min( + dfReducedRGB * dfRed * dfToneRed / dfToneMaxRGB, dfTopRGB); + padfDst[nGreenBand] = std::min( + dfReducedRGB * dfGreen * dfToneGreen / dfToneMaxRGB, dfTopRGB); + padfDst[nBlueBand] = std::min( + dfReducedRGB * dfBlue * dfToneBlue / dfToneMaxRGB, dfTopRGB); +#endif + + // Other bands processing (NIR, ...): only apply RGB reduction factor + if (bRGBBandsAreFirst) + { + // optimization + for (int iBand = 3; iBand < nInBands; ++iBand) + { + if (padfSrc[iBand] != padfInNoData[iBand]) + { + padfDst[iBand] = dfReducedRGB * padfSrc[iBand]; + } + else + { + bNoDataPixel = true; + break; + } + } + } + else + { + for (int iBand = 0; iBand < nInBands; ++iBand) + { + if (iBand != nRedBand && iBand != nGreenBand && + iBand != nBlueBand) + { + if (padfSrc[iBand] != padfInNoData[iBand]) + { + padfDst[iBand] = dfReducedRGB * padfSrc[iBand]; + } + else + { + bNoDataPixel = true; + break; + } + } + } + } + } + else + { + bNoDataPixel = true; + } + if (bNoDataPixel) + { + for (int iBand = 0; iBand < nInBands; ++iBand) + { + padfDst[iBand] = padfOutNoData[iBand]; + } + } + + padfSrc += nInBands; + padfDst += nInBands; + } + + return CE_None; +} + +/************************************************************************/ +/* GDALVRTRegisterDefaultProcessedDatasetFuncs() */ +/************************************************************************/ + +/** Register builtin functions that can be used in a VRTProcessedDataset. + */ +void GDALVRTRegisterDefaultProcessedDatasetFuncs() +{ + GDALVRTRegisterProcessedDatasetFunc( + "BandAffineCombination", nullptr, + "" + " " + " " + " " + " " + " " + " " + " " + "", + GDT_Float64, nullptr, 0, nullptr, 0, BandAffineCombinationInit, + BandAffineCombinationFree, BandAffineCombinationProcess, nullptr); + + GDALVRTRegisterProcessedDatasetFunc( + "LUT", nullptr, + "" + " " + " " + " " + "", + GDT_Float64, nullptr, 0, nullptr, 0, LUTInit, LUTFree, LUTProcess, + nullptr); + + GDALVRTRegisterProcessedDatasetFunc( + "Dehazing", nullptr, + "" + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + "", + GDT_Float64, nullptr, 0, nullptr, 0, DehazingInit, DehazingFree, + DehazingProcess, nullptr); + + GDALVRTRegisterProcessedDatasetFunc( + "Trimming", nullptr, + "" + " " + " " + " " + " " + " " + " " + " " + " " + " " + " " + "", + GDT_Float64, nullptr, 0, nullptr, 0, TrimmingInit, TrimmingFree, + TrimmingProcess, nullptr); +} diff --git a/frmts/vrt/vrtsources.cpp b/frmts/vrt/vrtsources.cpp index 856927a1300f..97bbd43de8af 100644 --- a/frmts/vrt/vrtsources.cpp +++ b/frmts/vrt/vrtsources.cpp @@ -285,11 +285,6 @@ void VRTSimpleSource::GetDstWindow(double &dfDstXOff, double &dfDstYOff, /* SerializeToXML() */ /************************************************************************/ -static const char *const apszSpecialSyntax[] = { - "NITF_IM:{ANY}:{FILENAME}", "PDF:{ANY}:{FILENAME}", - "RASTERLITE:{FILENAME},{ANY}", "TILEDB:\"{FILENAME}\":{ANY}", - "TILEDB:{FILENAME}:{ANY}"}; - static bool IsSlowSource(const char *pszSrcName) { return strstr(pszSrcName, "/vsicurl/http") != nullptr || @@ -350,11 +345,8 @@ CPLXMLNode *VRTSimpleSource::SerializeToXML(const char *pszVRTPath) } else { - for (size_t i = 0; - i < sizeof(apszSpecialSyntax) / sizeof(apszSpecialSyntax[0]); - ++i) + for (const char *pszSyntax : VRTDataset::apszSpecialSyntax) { - const char *const pszSyntax = apszSpecialSyntax[i]; CPLString osPrefix(pszSyntax); osPrefix.resize(strchr(pszSyntax, ':') - pszSyntax + 1); if (pszSyntax[osPrefix.size()] == '"') @@ -552,86 +544,8 @@ VRTSimpleSource::XMLInit(const CPLXMLNode *psSrc, const char *pszVRTPath, m_nExplicitSharedStatus = CPLTestBool(pszShared); } - if (pszVRTPath != nullptr && m_bRelativeToVRTOri) - { - // Try subdatasetinfo API first - // Note: this will become the only branch when subdatasetinfo will become - // available for NITF_IM, RASTERLITE and TILEDB - const auto oSubDSInfo{GDALGetSubdatasetInfo(pszFilename)}; - if (oSubDSInfo && !oSubDSInfo->GetPathComponent().empty()) - { - auto path{oSubDSInfo->GetPathComponent()}; - m_osSrcDSName = oSubDSInfo->ModifyPathComponent( - CPLProjectRelativeFilename(pszVRTPath, path.c_str())); - GDALDestroySubdatasetInfo(oSubDSInfo); - } - else - { - bool bDone = false; - for (size_t i = 0; - i < sizeof(apszSpecialSyntax) / sizeof(apszSpecialSyntax[0]); - ++i) - { - const char *pszSyntax = apszSpecialSyntax[i]; - CPLString osPrefix(pszSyntax); - osPrefix.resize(strchr(pszSyntax, ':') - pszSyntax + 1); - if (pszSyntax[osPrefix.size()] == '"') - osPrefix += '"'; - if (EQUALN(pszFilename, osPrefix, osPrefix.size())) - { - if (STARTS_WITH_CI(pszSyntax + osPrefix.size(), "{ANY}")) - { - const char *pszLastPart = strrchr(pszFilename, ':') + 1; - // CSV:z:/foo.xyz - if ((pszLastPart[0] == '/' || pszLastPart[0] == '\\') && - pszLastPart - pszFilename >= 3 && - pszLastPart[-3] == ':') - { - pszLastPart -= 2; - } - CPLString osPrefixFilename = pszFilename; - osPrefixFilename.resize(pszLastPart - pszFilename); - m_osSrcDSName = - osPrefixFilename + - CPLProjectRelativeFilename(pszVRTPath, pszLastPart); - bDone = true; - } - else if (STARTS_WITH_CI(pszSyntax + osPrefix.size(), - "{FILENAME}")) - { - CPLString osFilename(pszFilename + osPrefix.size()); - size_t nPos = 0; - if (osFilename.size() >= 3 && osFilename[1] == ':' && - (osFilename[2] == '\\' || osFilename[2] == '/')) - nPos = 2; - nPos = osFilename.find( - pszSyntax[osPrefix.size() + strlen("{FILENAME}")], - nPos); - if (nPos != std::string::npos) - { - const CPLString osSuffix = osFilename.substr(nPos); - osFilename.resize(nPos); - m_osSrcDSName = osPrefix + - CPLProjectRelativeFilename( - pszVRTPath, osFilename) + - osSuffix; - bDone = true; - } - } - break; - } - } - if (!bDone) - { - m_osSrcDSName = - CPLProjectRelativeFilename(pszVRTPath, pszFilename); - } - } - } - else - { - m_osSrcDSName = pszFilename; - } + m_osSrcDSName = VRTDataset::BuildSourceFilename( + pszFilename, pszVRTPath, CPL_TO_BOOL(m_bRelativeToVRTOri)); const char *pszSourceBand = CPLGetXMLValue(psSrc, "SourceBand", "1"); m_bGetMaskBand = false; diff --git a/gcore/gdal.h b/gcore/gdal.h index e98afdc0dd52..96baa8e55292 100644 --- a/gcore/gdal.h +++ b/gcore/gdal.h @@ -1621,6 +1621,115 @@ CPLErr CPL_DLL CPL_STDCALL GDALAddDerivedBandPixelFuncWithArgs( const char *pszName, GDALDerivedPixelFuncWithArgs pfnPixelFunc, const char *pszMetadata); +/** Generic pointer for the working structure of VRTProcessedDataset + * function. */ +typedef void *VRTPDWorkingDataPtr; + +/** Initialization function to pass to GDALVRTRegisterProcessedDatasetFunc. + * + * This initialization function is called for each step of a VRTProcessedDataset + * that uses the related algorithm. + * The initialization function returns the output data type, output band count + * and potentially initializes a working structure, typically parsing arguments. + * + * @param pszFuncName Function name. Must be unique and not null. + * @param pUserData User data. May be nullptr. Must remain valid during the + * lifetime of GDAL. + * @param papszFunctionArgs Function arguments as a list of key=value pairs. + * @param nInBands Number of input bands. + * @param eInDT Input data type. + * @param[in,out] padfInNoData Array of nInBands values for the input nodata + * value. The init function may also override them. + * @param[in,out] pnOutBands Pointer whose value must be set to the number of + * output bands. This will be set to 0 by the caller + * when calling the function, unless this is the + * final step, in which case it will be initialized + * with the number of expected output bands. + * @param[out] peOutDT Pointer whose value must be set to the output + * data type. + * @param[in,out] ppadfOutNoData Pointer to an array of *pnOutBands values + * for the output nodata value that the + * function must set. + * For non-final steps, *ppadfOutNoData + * will be nullptr and it is the responsibility + * of the function to CPLMalloc()'ate it. + * If this is the final step, it will be + * already allocated and initialized with the + * expected nodata values from the output + * dataset (if the init function need to + * reallocate it, it must use CPLRealloc()) + * @param pszVRTPath Directory of the VRT + * @param[out] ppWorkingData Pointer whose value must be set to a working + * structure, or nullptr. + * @return CE_None in case of success, error otherwise. + * @since GDAL 3.9 */ +typedef CPLErr (*GDALVRTProcessedDatasetFuncInit)( + const char *pszFuncName, void *pUserData, CSLConstList papszFunctionArgs, + int nInBands, GDALDataType eInDT, double *padfInNoData, int *pnOutBands, + GDALDataType *peOutDT, double **ppadfOutNoData, const char *pszVRTPath, + VRTPDWorkingDataPtr *ppWorkingData); + +/** Free function to pass to GDALVRTRegisterProcessedDatasetFunc. + * + * @param pszFuncName Function name. Must be unique and not null. + * @param pUserData User data. May be nullptr. Must remain valid during the + * lifetime of GDAL. + * @param pWorkingData Value of the *ppWorkingData output parameter of + * GDALVRTProcessedDatasetFuncInit. + * @since GDAL 3.9 + */ +typedef void (*GDALVRTProcessedDatasetFuncFree)( + const char *pszFuncName, void *pUserData, VRTPDWorkingDataPtr pWorkingData); + +/** Processing function to pass to GDALVRTRegisterProcessedDatasetFunc. + * @param pszFuncName Function name. Must be unique and not null. + * @param pUserData User data. May be nullptr. Must remain valid during the + * lifetime of GDAL. + * @param pWorkingData Value of the *ppWorkingData output parameter of + * GDALVRTProcessedDatasetFuncInit. + * @param papszFunctionArgs Function arguments as a list of key=value pairs. + * @param nBufXSize Width in pixels of pInBuffer and pOutBuffer + * @param nBufYSize Height in pixels of pInBuffer and pOutBuffer + * @param pInBuffer Input buffer. It is pixel-interleaved + * (i.e. R00,G00,B00,R01,G01,B01, etc.) + * @param nInBufferSize Size in bytes of pInBuffer + * @param eInDT Data type of pInBuffer + * @param nInBands Number of bands in pInBuffer. + * @param padfInNoData Input nodata values. + * @param pOutBuffer Output buffer. It is pixel-interleaved + * (i.e. R00,G00,B00,R01,G01,B01, etc.) + * @param nOutBufferSize Size in bytes of pOutBuffer + * @param eOutDT Data type of pOutBuffer + * @param nOutBands Number of bands in pOutBuffer. + * @param padfOutNoData Input nodata values. + * @param dfSrcXOff Source X coordinate in pixel of the top-left of the region + * @param dfSrcYOff Source Y coordinate in pixel of the top-left of the region + * @param dfSrcXSize Width in pixels of the region + * @param dfSrcYSize Height in pixels of the region + * @param adfSrcGT Source geotransform + * @param pszVRTPath Directory of the VRT + * @param papszExtra Extra arguments (unused for now) + * @since GDAL 3.9 + */ +typedef CPLErr (*GDALVRTProcessedDatasetFuncProcess)( + const char *pszFuncName, void *pUserData, VRTPDWorkingDataPtr pWorkingData, + CSLConstList papszFunctionArgs, int nBufXSize, int nBufYSize, + const void *pInBuffer, size_t nInBufferSize, GDALDataType eInDT, + int nInBands, const double *padfInNoData, void *pOutBuffer, + size_t nOutBufferSize, GDALDataType eOutDT, int nOutBands, + const double *padfOutNoData, double dfSrcXOff, double dfSrcYOff, + double dfSrcXSize, double dfSrcYSize, const double adfSrcGT[/*6*/], + const char *pszVRTPath, CSLConstList papszExtra); + +CPLErr CPL_DLL GDALVRTRegisterProcessedDatasetFunc( + const char *pszFuncName, void *pUserData, const char *pszXMLMetadata, + GDALDataType eRequestedInputDT, const GDALDataType *paeSupportedInputDT, + size_t nSupportedInputDTSize, const int *panSupportedInputBandCount, + size_t nSupportedInputBandCountSize, + GDALVRTProcessedDatasetFuncInit pfnInit, + GDALVRTProcessedDatasetFuncFree pfnFree, + GDALVRTProcessedDatasetFuncProcess pfnProcess, CSLConstList papszOptions); + GDALRasterBandH CPL_DLL CPL_STDCALL GDALGetMaskBand(GDALRasterBandH hBand); int CPL_DLL CPL_STDCALL GDALGetMaskFlags(GDALRasterBandH hBand); CPLErr CPL_DLL CPL_STDCALL GDALCreateMaskBand(GDALRasterBandH hBand,