Plugin for element-wise arithmetic operators with seamless compute dispatching to GPU or CPU. More...

Collaboration diagram for ImageMath Plugin:

Detailed Description

Plugin for element-wise arithmetic operators with seamless compute dispatching to GPU or CPU.

Overview of main components of the ImFusionImageMath module. Does not show all classes and some intermediate parts of the inheritence hierarchy are omitted for simplicity.

Overview and quick start

This plugin provides element-wise arithmetic operators for SharedImage, SharedImageSet, and TypedImage.

When to use SharedImageSetArithmetics?

You need element-wise arithmetic operations for our image types. Trivial example: You want to add, subtract, multiply, or divide two images.
All images/extracted blocks have the same dimensions.
You want GPU and CPU support without writing a single shader.
You do not want to manually worry about shift, scale, intensity shifts/scales, or selection.

How to use SharedImageSetArithmetics?

Just include "ImFusion/ImageMath/SharedImageSetArithmetic.h" and you are good to go. You can also use ImageMath expressions as abstract includes in your own shaders.

Basic usage examples

Build-in operators:

In ExprBaseTpl we provide a lot of build-in operators.

// Any right-hand leaf can also be ChannelWiseScalar, Scalar, or Variable
SharedImage res, img1, img2, img3;
 
// Assuming all images have the same type. Otherwise see below.
makeArray(res) = makeArray(img1) + makeArray(img2);
makeArray(res) = makeArray(img1) - makeArray(img2);
makeArray(res) = makeArray(img1) * makeArray(img2);
makeArray(res) = makeArray(img1) / makeArray(img2);
 
 
// We also support comparison operators. The result will be either 0 or 1 such that it can be used as mask by multiplication
makeArray(res) = makeArray(img1) < makeArray(img2);
makeArray(res) = makeArray(img1) <= makeArray(img2);
makeArray(res) = makeArray(img1) > makeArray(img2);
makeArray(res) = makeArray(img1) >= makeArray(img2);
makeArray(res) = makeArray(img1) == makeArray(img2);
 
// If img1 is a 1-channel operator, analogously to the ternary ?: operator we can select either img2 or img3 depending on the values in img1
 
makeArray(res) = makeArray(img1).select(makeArray(img2), makeArray(img3));
 
// Channels can be combined from multiple expressions. Warning: expressions must have compatible image descriptors,
// and the total number of channels may not exceed 4 when working on the GPU.
// The below also demonstrates the evaluateIntoImage operator, which allocates a new image and evaluates the expression into it.
 
std::unique_ptr<SharedImage> multiChannelImage = combineChannels(makeArray(img1),makeArray(img2)).evaluateIntoImage();
 
// Short-cuts. Feel free to extend these in ExprBaseTpl.h
makeArray(res) = makeArray(img1).abs();
makeArray(res) = makeArray(img1).exp();
makeArray(res) = makeArray(img1).log();
makeArray(res) = makeArray(img1).pow2();
makeArray(res) = makeArray(img1).sqrt();
makeArray(res) = makeArray(img1).cast<float>();
makeArray(res) = makeArray(img1).divideIfNotZero(makeArray(img2));
makeArray(res) = makeArray(img1).min(makeArray(img2));
makeArray(res) = makeArray(img1).max(makeArray(img2));
makeArray(res) = makeArray(img1).clamp(makeArray(img2), makeArray(img3));

Special array-like leaves

Coordinates.

// The following writes the first component of the world coordinates on `img1` into `res`
// Image, Pixel, and Texture coordinates are also supported.
makeArray(res) = makeCoordinate(img1, ImageMath::CoordinateType::World).x();
auto worldCoords = makeCoordinates(sis, CoordinateType::World);
auto circularMask = (2.0 * (makeCoordinates(sis, CoordinateType::Texture).xy() - vec2(0.5, 0.5))).length() < 1.0;
auto coords = makeCoordinates(labelSis, CoordinateType::Image);
Eigen::VectorXd upperBound = (makeArray(labelSis) == 1).select(coords, -labelSis.extent()).max();
auto mask = (makeArray(labelSis) == 1);                         // 1 channel binary mask expression
auto coords = makeCoordinates(labelSis, CoordinateType::Image); // 3 channels coordinate expression
double count = mask.sum()[0];                                   // scalar Σ1, evaluation of the mask expression
Eigen::VectorXd summed = mask.select(coords, 0.0).sum();        // 3-vector [Σx, Σy, Σz], evaluates the coords expression where mask is true
                                                                // Zeros are streamed implicitly over 3 channels
Eigen::VectorXd centerOfMass = summed / count;                  // assumes count > 0
auto mask = (makeArray(labelSis) == 1);                         // 1 channel binary mask expression
// append homogeneous 1 to coordinates to form 4 channels; stream the mask to 4 channels to match (optional)
auto zeros4 = makeScalar(0.0).streamChannels(4);                // 4-channel zero vector
auto homogeneousCoords = combineChannels(makeCoordinates(labelSis, CoordinateType::Pixel), 1.0);
auto coordForLabel = mask.select(homogeneousCoords, zeros4);    // equivalent to mask.select(homogeneousCoords, 0.0);
Eigen::VectorXd accum = coordForLabel.sum();                    // [Σx, Σy, Σz, Σ1], evaluates the expression
Eigen::VectorXd centerOfMass = accum.topRows<3>() / accum[3];   // assumes accum[3] > 0 — equivalent to accum.hnormalized().head<3>()
auto intensityForLabel = mask.select(makeArray(sis), 0.0);                          // 1 channel expression
Eigen::VectorXd weightedSum = (coords * intensityForLabel.streamChannels(3)).sum(); // 3-vector, evaluation 1
double massSum = intensityForLabel.sum()[0];                                        // total mass scalar, evaluation 2
Eigen::VectorXd centerOfMass = weightedSum / massSum;                               // assumes massSum > 0
auto homogeneousCoords = combineChannels(makeCoordinates(labelSis, CoordinateType::Pixel), 1.0);  // 4-channel coordinate expression
auto coordTimesMass4 = intensityForLabel.streamChannels(4) * homogeneousCoords;                   // weighted coordinates expression
auto weightedHom = mask.select(coordTimesMass4, 0.0);                                             // implicit channel streaming of 0.0
Eigen::VectorXd accumWeighted = weightedHom.sum();                                                // [Σx·m, Σy·m, Σz·m, Σm]
Eigen::VectorXd centerOfMass = accumWeighted.topRows<3>() / accumWeighted[3];                     // assumes accumWeighted[3] > 0

Combining channels:

auto mergeChannels = combineChannels(makeArray(img1), makeArray(img2), makeArray(img3), 1.0);
auto addCoords = combineChannels(makeArray(img1), makeCoordinates(img1, CoordinateType::World));

Mask/GlMask values:

std::shared_ptr<CroppingMask> m = std::make_shared<CroppingMask>(si.descriptor().dimensions());
m->setRoundness(100);
auto maskLeaf = ImageMath::makeMaskLeaf( ImageMath::makeArray(si), m);
auto derivedExpr = (maskLeaf > 0).select(1.0,2.0);

You can also use ImageMath expressions to create a GlMask, see OperandMask.h for further information.

Deformations:

std::shared_ptr<Deformation> def = si.deformation();
auto defExpr = DeformationLeaf(si,def);
auto defExprX = defExpr.x();

You can also use ImageMath expressions to create a Deformation, see OperandDeformation.h for further information.

Reduction operators:

Regular:

// Compute channel wise reduction
Eigen::VectorXd reducedSum = makeArray(img1).sum();
// Computed squared sum
Eigen::VectorXd l2norms = makeArray(img1).pow2().sum();
// Number of pixel bigger than
Eigen::VectorXd biggerThan5Sum= (makeArray(img1) > 5).sum();

Visitors in max/min reductions:

std::vector<vec4i> pixelCoordsMax; // Fourth component is SharedImage index for SharedImageSet expressions

makeArray(img1).max(pixelCoordsMax);

std::vector

Reductions over channels

We support common reductions (sum, product, min, max, length) over channels.

In the following example, si is e.g. a 3-channel image:

auto grayScale = makeArray(si).channelSum();
auto distanceToPoint = (makeCoordinates(sis, CoordinateType::World) - pt).length();
auto circularMask = makeCoordinates(sis, CoordinateType::World).dot(vec3(1.0,1.0,1.0));

View operators:

We currently offer Views for block access as well as channel swizzling.

// block read/write access. E.g. to compute forward differences
// Note: Currently dimension squeezing is not implemented, i.e. if the block is 2D you still need to match the dimensions.
makeArray(img1).block(vec3i(1, 0, 0), vec3i(width - 1, height, slices)) -=
    makeArray(vol).block(vec3i(0, 0, 0), vec3i(width - 1, height, slices));
 
// channel-wise read access. E.g. to compute total variation
makeArray(res) =
   (makeArray(img1).channelSwizzle({0}).pow2() + makeArray(img1).channelSwizzle({1}).pow2() +
    makeArray(img1).channelSwizzle({2}).pow2())
      .sqrt();
 
// channel-wise write access. E.g. to merge-channels
makeArray(res).channelSwizzle({0}) = makeArray(img1);
makeArray(res).channelSwizzle({1}) = makeArray(img2);
makeArray(res).channelSwizzle({2}) = makeArray(img3);
 
// Convenience functions for common channel swizzles
// This is implemented also for .r(), .g(), .b(), .xy(), .rgba(), .xyzw(), etc ...
 
makeArray(res).x() = makeArray(img1);
makeArray(res).y() = makeArray(img2);
makeArray(res).z() = makeArray(img3);

Broadcasting

We support broadcasting/streaming in several settings:

You can broadcast a single image to all images in a SharedImageSet.

SharedImage si1...; //Single image
SharedImageSet sis...; // SharedImageSet with multiple images
makeArray(sis) *= makeArray(si1);

You can broadcast a single channel to multiple channels.

SharedImage si1...; // Single channel image
SharedImage si2...; // Image with 4 channels.
makeArray(si2) *= makeArray(si1).streamChannels(4);

You can broadcast any dimension using samplers, i.e. boundary conditions. We support clampToEdge (default) or clampToBorder with a specific border color. E.g. you can stream an image of a single slice to all slices of a volume like this:

SharedImage si1...; // Image with dimensions 4,5,1
SharedImage si2...; // Image with dimensions 4,5,6
makeArray(si2) *= makeArray(si1).block(vec3i::Zero(), vec3i(4,5,6));

You can combine channelSwizzle and broadcasting to stream to other slices as well:

SharedImage si1...; // Image with dimensions 4,6,1
SharedImage si2...; // Image with dimensions 4,5,6
makeArray(si2) *= makeArray(si1).block(vec3i::Zero(), vec3i(4,6,5)).channelSwizzle({0,2,1});

You can also define a border color instead of using ClampToEdge. For instance to realize finite differences with zero padding:

SharedImage si2...; // Image with dimensions 4,5,6
// Forward differences in x
makeArray(si2) -= makeArray(si2, false, MagFilter::Nearest, Wrap::ClampToBorder, vec4f::Zero()).block(vec3i(-1,0,0), vec3i(4,6,5));

Compute 3D finite differences:

SharedImage si2...; // Image with dimensions 4,5,6
SharedImage grad...; // Image with dimensions 4,5,6 and 3 channels
 
makeArray(grad) = makeArray(si2).streamChannels(3) - combineChannels(
    makeArray(si2, false, MagFilter::Nearest, Wrap::ClampToBorder, vec4f::Zero()).block(vec3i(-1,0,0), vec3i(4,6,5)),
    makeArray(si2, false, MagFilter::Nearest, Wrap::ClampToBorder, vec4f::Zero()).block(vec3i(0,-1,0), vec3i(4,6,5)),
    makeArray(si2, false, MagFilter::Nearest, Wrap::ClampToBorder, vec4f::Zero()).block(vec3i(0,0,-1), vec3i(4,6,5)));

Evaluating an expression into an empty image

We provide the evaluateIntoImage template to conveniently evaluate an expression into a new image:

std::shared_ptr<SharedImageSet> result = (makeArray(sis0) + makeArray(sis1)).evaluateIntoImage<float>();

Storing an expression

There are several use-cases where you might want to store the intermediate object, representing an expression. In this example we store the expression representing a gradient descent step to reuse it in each iteration.

SharedImageSet x, grad;
 
// These store the expressions (!) not the results
auto xExpr = makeArray(x);
auto stepLengthExpr = makeScalar(1.0);
auto gradStep = xExpr - stepLengthExpr * makeArray(grad);
 
for(int i = 0; i < N; ++i)
{
    // update grad
    ...
    // update stepLength
    stepLengthExpr.val() = ...;
 
    // Perform gradient descent. At this point the expression will be evaluated.
    xExpr = gradStep;
}

Implementation details

This plugin provides element-wise arithmetic operators for SharedImage, SharedImageSet, and TypedImage. The evaluation of an expression is Selection, shift/scale and type aware.

Expression tree

For each expression an binary expression tree gets build. Leaves are values, inner nodes are operators. This either happens at compile time or at runtime (see below for more details).

Expression Leaves

These leaves are currently supported:

R-value leaves: ChannelWiseScalar, Scalar, Variable, Array (including block access), or internal::ChannelSwizzleViewTpl
L-value leaves: Array (including block access), or internal::ChannelSwizzleViewTpl

Compile-Time Expressions

We support compile-time expressions. These expressions need to be fully posed and can be evaluated without the need for any temporary images, i.e. memory, to store intermediate results. In order to enable this we use the curiously recurring template pattern (CRTP) (very similar to how Eigen works). The corresponding base class is ExprBaseTpl. Child classes are for instance the aforementioned expression leaves but also templated implementations of arithmetic operations.

We support the following operators: internal::BinaryOpTpl, internal::UnaryOpTpl, internal::VariableSubstitutionOp, internal::PolyWrapper, internal::DeviceStrategyWrapperTpl, and internal::ChannelSwizzleViewTpl. In almost all cases as a user you do not need these classes but you instead want to use the corresponding methods provided by ExprBaseTpl:

Run-Time Expressions

In order to use expression as arguments to methods or functions without the need for templates we also support run-time expressions. The major difference is that these are not compile time optimized and require temporary image allocation in order to store intermediate results while evaluating an expression.

The base class for polymorphic run-time expressions is ExprBase.

Compile-time expressions can be used in run-time expressions as well. For this purpose we offer a wrapper class internal::PolyWrapper for which we also offer the builder method ExprBaseTpl::polyWrapper().

GlExpr

Finally any expression can be used within OpenGL shaders. For this purpose we have the GlExpr which is a specialized GlAbstractInclude.

#version 400
in vec3 in_texCoord;
 
layout(location = 0) out vec4 out_color;
 
#pragma abstract_include "EXPR"
#ifndef EXPR
// default implementation if no definition was set
vec4 exprVal(vec3 tc, vec4 val, int index)
{
    return vec4(0.0);
}
#endif
 
void main()
{
    vec4 varVal = some_fancy_computation();
 
    out_color = exprVal(in_texCoord, varVal, 0);
}

 
auto compute = [...](auto expr)
{
    GlExpr exprInInclude(exprIn);
    exprInInclude.setIndices({0});
 
    prog.addAbstractInclude(&exprInInclude);
 
    prog.enable();
    exprInInclude.setIncludeArguments(prog);
    ...
    prog.compute(*siOut.gl());
    ...
    prog.disable();
};
 
// use compute with different expressions
 
// siOut will contain si1 + si2 + some_fancy_computation()
compute(makeArray(si1) + makeArray(si2) + Variable());
 
// siOut will contain si1 + si2 + some_fancy_computation()^2
compute(makeArray(si1) + makeArray(si2) + Variable().pow2());
 
// siOut will contain some_fancy_computation() * si1 + si2
compute(Variable() * makeArray(si1) + makeArray(si2));

NOTE: This also supports multiple input and output images. See the documentation of GlExpr.

GPU and CPU support and how is determined where to execute the expression

The DeviceStrategy of the target as well as the source expression will define where the expression gets evaluated. If the DeviceStrategy is set to DeviceStrategy::Auto the expression will be run on the GPU if the target image has a GlImage representation. Otherwise the CPU execution will be used. This can be overridden with ExprBaseTpl::forceCPU() or ExprBaseTpl::forceGPU(). DeviceStrategy::ForceCPU will override DeviceStrategy::ForceGPU and both will override DeviceStrategy::Auto.

Array returns DeviceStrategy::ForceCPU as DeviceStrategy if the image has more than 4 channels or has type double, uint, or int.

Handling mixed types

Similar to the case above, the type is determined using the target image. The main reason is that you cannot build up the expression for any permutation of types because this will become to complex. Thus, as long as all the images involved have the same type, you do not have to provide any type information. I.e. You can write the expression like this.

makeArray(res) = makeArray(img1) + makeArray(img2);

However, if you mix types like in the following example you need to provide the type for all images which are not the same type as the target image. Note that the labelMap is a unsigned short image while img1 is float. In such cases you need to provide the type.

// Extract an intensity range as label
makeArray(labelMap) = makeScalar(0)
    + label
      makeArray<float>(img1).binaryOp<ImgOps::thresholdMaskGreater>(rangeLower)
          makeArray<float>(img1).binaryOp<ImgOps::thresholdMaskLessEqual>(rangeUpper);

Topics
	Core/Low-level API
	Low-level data structures for element-wise arithmetic operators.

	Extending ImageMath
	How to extend the ImageMath plugin.

	Pitfalls
	Common pitfalls.

Files
file	SharedImageSetArithmetic.h
	This header provides arithmetic operators for TypeImage<>, SharedImage, and SharedImageSet.

Namespaces
namespace	ImFusion::ImageMath
	Arithmetic operations on images and arrays.

Classes
class	SharedImageArithmeticAlgorithm
	ImageMath algorithm class. More...

class	GlExpr
	Wrapper for using a ExprBase as GlAbstractInclude. More...

class	ExprBase
	Base class for expressions (non-templated) (needs to be empty in order to enable empty base class optimization!) More...

class	PolyWrapper< ExprT >
	Specialization of ExprBase wrapping a ExprBaseTpl. More...

class	Array< T, ImgType, BlockAccessEnabled >
	Specialization of ExprBase for a image leaf, i.e. More...

class	ChannelWiseScalar< T, Dim >
	Specialization of ExprBaseTpl for a channel wise scalar value. More...

class	ExprBaseTpl< ImplType >
	Base class for compile-time optimized expressions (curiously recurring template pattern (CRTP)) More...

class	Scalar< T >
	Specialization of ExprBaseTpl for a single scalar value. More...

class	StringLeaf
	Specialization of ExprBaseTpl for a simple leaf representing a string. More...

class	Variable
	Specialization of ExprBaseTpl for a variable, e.g. More...

Functions
template<typename T>
auto	makeArray (TypedImage< T > &img, bool ignoreShiftAndScale, MagFilter magFilter, Wrap wrap, const vec4f &borderColor)
	Convenience function to create an array leaf (needed as until c++17 there is no template type deduction from constructor arguments)

template<typename T>
auto	makeArray (SharedImage &img, bool ignoreShiftAndScale, MagFilter magFilter, Wrap wrap, const vec4f &borderColor)
	Convenience function to create an array leaf (needed as until c++17 there is no template type deduction from constructor arguments)

template<typename T>
auto	makeArray (SharedImageSet &img, bool ignoreShiftAndScale, MagFilter magFilter, Wrap wrap, const vec4f &borderColor)
	Convenience function to create an array leaf (needed as until c++17 there is no template type deduction from constructor arguments)

template<typename T>
auto	makeArray (const TypedImage< T > &img, bool ignoreShiftAndScale, MagFilter magFilter, Wrap wrap, const vec4f &borderColor)
	Convenience function to create an array leaf (needed as until c++17 there is no template type deduction from constructor arguments)

template<typename T>
auto	makeArray (const SharedImage &img, bool ignoreShiftAndScale, MagFilter magFilter, Wrap wrap, const vec4f &borderColor)
	Convenience function to create an array leaf (needed as until c++17 there is no template type deduction from constructor arguments)

template<typename T>
auto	makeArray (const SharedImageSet &img, bool ignoreShiftAndScale, MagFilter magFilter, Wrap wrap, const vec4f &borderColor)
	Convenience function to create an array leaf (needed as until c++17 there is no template type deduction from constructor arguments)

template<typename T, int Dim>
ChannelWiseScalar< T, Dim >	makeChannelWiseScalar (Eigen::Matrix< T, Dim, 1 > val)
	Convenience function to create a channel wise scalar leaf (needed as until c++17 there is no template type deduction from constructor arguments)

template<typename T>
Coordinates< const TypedImage< T > >	makeCoordinates (const TypedImage< T > &img, CoordinateType coordType=CoordinateType::World)
	This function creates a suitable Coordinates object.

template<typename T>
Noise< const TypedImage< T > >	makeNoise (const TypedImage< T > &img, NoiseType noiseType=NoiseType::Gaussian, float alpha=0.0, float beta=1.0, int seed=0)
	This function creates a suitable Noise object.

template<typename T>
Scalar< typename std::decay_t< T > >	makeScalar (T &&val)
	Convenience function to create a scalar leaf (needed as until c++17 there is no template type deduction from constructor arguments)

template<typename T>
Scalar< typename std::decay_t< T > & >	makeScalarWrapper (T &val)
	Convenience function to create a scalar leaf wrapping a scalar variable (needed as until c++17 there is no template type deduction from constructor arguments)

Function Documentation

◆ makeCoordinates()

template<typename T>

Coordinates< const TypedImage< T > > makeCoordinates	(	const TypedImage< T > &	img,
		CoordinateType	coordType = CoordinateType::World )

#include <ImFusion/ImageMath/Tpl/Coordinates.h>

This function creates a suitable Coordinates object.

It is needed, as C++ only allows template deduction from constructor arguments in C++17

See also: Coordinates::Coordinates

◆ makeNoise()