VRT Pixel Functions: Add function to evaluate arbitrary expression

[dbaston]

What does this PR do?

Adds a C++ pixel function called "expression" that can evaluate an arbitrary expression (or indeed, a mini-program) using the exprtk library. This is a single MIT-licensed header that is easily integrated into GDAL. It appears to be quite widely used. However, given that a major purpose of this is to extend the utility of VRT to users that don't want to allow functions defined in Python (possibly because of security concerns), I'm a little uncomfortable with the scope of expressions allowed by this library. Some of these, such as file I/O, can be disabled (and I've done so.)

I've also looked at:

• tinyexpr, which seems easily integrated but lacks conditionals
• muparser, a little harder to integrate (includes multiple header/object files) and supports a reasonable variety of functions including a ternary conditional operator

I need to do more work to investigate details such as nodata handling, but wanted to share the concept sooner rather than later.

An alternative implementation would expose exprtk as an additional pixel function language alongside Python, instead of hooking into the existing C++ pixel function machinery.

An example VRT that is allowed is:

<VRTDataset rasterXSize="1" rasterYSize="1">
  <VRTRasterBand dataType="Float64" band="1" subClass="VRTDerivedRasterBand">
    <PixelFunctionType>expression</PixelFunctionType>
    <PixelFunctionArguments expression="(NIR-R)/(NIR+R)" />
    <SimpleSource name="NIR">
      <SourceFilename relativeToVRT="0">source_0.tif</SourceFilename>
      <SourceBand>1</SourceBand>
    </SimpleSource>
    <SimpleSource name="R">
      <SourceFilename relativeToVRT="0">source_1.tif</SourceFilename>
      <SourceBand>1</SourceBand>
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

[rouault]

That's super cool, and so compact as an integration!

Any hints on the performance ? In particular it would be cool if exprtk could use SIMD to vectorize computations, but I'm probably asking too much. Regarding performance, I'd suspect the biggest issue to be the GDALCopyWords() done for each pixel. Creating a temporary array with all the double values and copying in one-go to the output data type would certainly be beneficial (this remark applies to existing pixel functions).

Can exptrtk can be used with non-double variables ? (in case that would make things faster for integer only computations)

Thinking out loud: for ultimate performance, we should probabl use some LLVM-based solution where the expression would be just-in-time compiled to the target architecture, a bit like numba does.

[rouault]

Maybe the following to avoid namespace collisions in case GDAL is integrated with another software also using exprtk?

Suggested change

	#include <exprtk.hpp>
	#define exprtk gdal_exprtk
	#include <exprtk.hpp>

[dbaston]

Any hints on the performance ?

Haven't looked into this yet. In particular, comparing to Python would be interesting.

Creating a temporary array with all the double values and copying in one-go to the output data type would certainly be beneficial (this remark applies to existing pixel functions).

Maybe they could be rewritten as templates...

Can exptrtk can be used with non-double variables ? (in case that would make things faster for integer only computations)

The docs note "Support for various numeric types (float, double, long double, MPFR/GMP)"

Thinking out loud: for ultimate performance, we should probabl use some LLVM-based solution where the expression would be just-in-time compiled to the target architecture, a bit like numba does.

I also came across mathpresso that does something along these lines.

[coveralls]

coverage: 72.812% (-0.8%) from 73.626%
when pulling 43a7c75 on dbaston:vrt-expression
into c578e9d on OSGeo:master.

[dbaston]

@rouault A limitation of the pixel function integration is that we can only output a single band. exprtk would happily accept an expression like "var q[2] := {B2, B3}; B1 * q;" but we have no way to return the second band to the user. Do you think looking into a VRTProcessedDataset integration might make sense instead?

On performance, I've done basic tests like summing up a 24-band netCDF. Performance is roughly equivalent to Python with perf showing the following for exprtk:

  23.81%  gdal_translate  libgdal.so.36.3.11.0        [.] ExprPixelFunc
  15.50%  gdal_translate  libgdal.so.36.3.11.0        [.] exprtk::details::vec_add_op<double>::process
  12.32%  gdal_translate  libgdal.so.36.3.11.0        [.] (anonymous namespace)::GDALCopyWordsFromT<short>
  11.21%  gdal_translate  libgdal.so.36.3.11.0        [.] GDALCopyWords64
   6.31%  gdal_translate  libnetcdf.so.19             [.] ncx_getn_short_short
   4.06%  gdal_translate  [unknown]                   [k] 0xffffffff98a0109c
   3.31%  gdal_translate  libc.so.6                   [.] __memset_avx2_unaligned_erms
   2.47%  gdal_translate  [unknown]                   [k] 0xffffffff98a00158

and numpy:

  16.49%  gdal_translate  libc.so.6                                          [.] __memmove_avx_unaligned_erms
  13.34%  gdal_translate  _multiarray_umath.cpython-310-x86_64-linux-gnu.so  [.] DOUBLE_add_FMA3__AVX2
  12.76%  gdal_translate  libgdal.so.36.3.11.0                               [.] (anonymous namespace)::GDALCopyWordsFromT<short>
   6.90%  gdal_translate  [unknown]                                          [k] 0xffffffff98a00158
   6.82%  gdal_translate  libnetcdf.so.19                                    [.] ncx_getn_short_short
   4.28%  gdal_translate  [unknown]                                          [k] 0xffffffff98a0109c
   4.23%  gdal_translate  ld-linux-x86-64.so.2                               [.] do_lookup_x
   3.13%  gdal_translate  libc.so.6                                          [.] __memset_avx2_unaligned_erms
   1.98%  gdal_translate  libpython3.10.so.1.0                               [.] _PyEval_EvalFrameDefault

[rouault]

Do you think looking into a VRTProcessedDataset integration might make sense instead?

It accepts only a single (potentially multiband) input dataset. That said the Input can be a VRTDataset, so you can assemble several single-band datasets as a virtual multiple one. Hard to tell which way is the preferred one. Would having it available in both VRTDerivedRasterBand and VRTProcessedDataset be overkill ? Having it available in VRTDerivedRasterBand makes it probably easier/intutive to write formulas that are different per output-band, even if exprtk allows to return a tuple.

[dbaston]

Would having it available in both VRTDerivedRasterBand and VRTProcessedDataset be overkill ?

Not necessarily. I'll look into VRTProcessedDataset some more.

Having it available in VRTDerivedRasterBand makes it probably easier/intutive to write formulas that are different per output-band

On the other hand, this requires repeating all of the SimpleSource definitions (and re-reading the input pixels) for every output band...unless I'm missing something.

[rouault]

On the other hand, this requires repeating all of the SimpleSource definitions (and re-reading the input pixels) for every output band.

that's true