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: 70.096% (+0.003%) from 70.093%
when pulling 54fd38f on dbaston:vrt-expression
into ee04fce 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

[rouault]

That's interesting... from a security point of view.

• can that trigger arbitrary memory allocation ?
• does exprtk do bound checking when accessing the arrays ?

I suspect you could also cause infinite /very long looping ?

[ArashPartow]

@rouault - Yes ExprTk has runtime checks. One of which is explicitly for vector index access.

An example can be found here: https://www.partow.net/programming/exprtk/index.html#simple_example_20

More detailed documentation regarding the checks can be found in Section 24 of the readme: https://www.partow.net/programming/exprtk/readme.html#line_4696

[dbaston]

can that trigger arbitrary memory allocation

Array sizes can be up to 2 billion, so you can cause a large allocation with:

var out[2000000000];
return [out];

does exprtk do bound checking when accessing the arrays ?

Not by default, but enabled in fd94c9e

[rouault]

Array sizes can be up to 2 billion, so you can cause a large allocation with:

could we improve exprtk to offer a way to limit the allocations?

[ArashPartow]

could we improve exprtk to offer a way to limit the allocations?

@rouault: Yes the maximum vector size that can be created in the expression can be set explicitly:

https://www.partow.net/programming/exprtk/readme.html#line_1473

parser_t parser; 
parser.settings().set_max_local_vector_size(12345); 

[ArashPartow]

@rouault - Yes ExprTk has runtime checks. One of which is explicitly for vector index access.

An example can be found here: https://www.partow.net/programming/exprtk/index.html#simple_example_20

More detailed documentation regarding the checks can be found in Section 24 of the readme: https://www.partow.net/programming/exprtk/readme.html#line_4696

[rouault]

The ASAN failures are unrelated to this PR and are avoided in master per fbaa99c. You can rebase on top of that.

Perhaps we could limit the amount of memory to some default value, let's say 10 MB, and if above that, fail with an error like "you need to set the GDAL_VRT_MAX_MEMORY_EXPR configuration option".

[rouault]

you'll need to rebase your PR on top of latest master and run scripts/collect_config_options.py

[rouault]

For the msys2-mingw64 "'file too big" error, see https://github.com/qgis/QGIS/blob/ac7a7fc78eadae3682c57387d2863128f6e65901/src/app/CMakeLists.txt#L493

[dbaston]

@dbaston I've pushed to your branch a commit (cherry-picked from master) that should hopefully fix the mysterious unrelated crash under the ASAN CI config (no idea why it just popped up on your pull request!)

Sorry, I didn't notice this message and apparently clobbered your commit with a force-push. I'll rebase from master again to pick it up.

Is the commit history in this PR important to you? I think the simplest way to remove the vendored exprtk from the history is to squash to a single commit. (Or maybe two commits, with the patch to the vendored muparser kept separate)

[rouault]

Is the commit history in this PR important to you?

Squashing in a single commit is probably going to be the simplest solution . I tend to prefer atomic history in general like "Add function foo()", "use function foo in driver A", "use function foo in driver B", but only if each individual commit builds properly and doesn't break a git bisect section, and when doing week-long developments it can be a headache to reconstruct a clean history (I'd wish github had a function where we could ask it to run CI on each commit to ensure the history is clean). We're not yet at a position to be as demanding as the Linux kernel project is to its contributors :-)

[rouault]

We'll need to document the variables related to EXPRTK and MUPARSER dependencies in https://gdal.org/en/latest/development/building_from_source.html#cmake-package-dependent-options

[rouault]

dbaston#2 : CI: add muparser to conda windows build

[rouault]

dbaston#3: ExprPixelFunc(): fix memory leak in error code path

[rouault]

The PR LGTM. Does @hobu want to give it some additional review ?

[hobu]

Is one of EXPRTK or MUPARSER a hard dependency of GDAL going forward?

[dbaston]

Is one of EXPRTK or MUPARSER a hard dependency of GDAL going forward?

They are both optional. CMake is setup to build with muparser if it's available, and to build with ExprTk only if requested (because of the effect on binary size). I would recommend that standard packages (conda, debian, etc.) build with muparser.

[hobu]

Is the "commonly expected usage" going to require an expression engine? IMO if that's the case, we should just vendor ExprTk and be done with it.

[dbaston]

I think the default would be to have an expression engine, yes. I don't recommend building w/ExprTk by default because it adds ~7MB to a ~16MB library. We could vendor muparser but it seems to be widely available (conda, msys2, debian, fedora, alpine).

[rouault]

Merged!

[rouault]

We now need a gdal raster calc --input=A:in.tif --input=B:in2.tif --output out.tif --calc "A+B" or shorter gdal raster calc A:in.tif B:in2.tif out.tif --calc "A+B"

[dbaston]

That's the plan! (#10067)