Convolution support and kernels in MArray expressions
[Multidimensional Dynamic Arrays]

Classes

• struct ltl::cderiv1
• struct ltl::cderiv2
• struct ltl::cderiv3
• struct ltl::cderiv4
• struct ltl::fderiv1h
• struct ltl::fderiv2h
• struct ltl::fderiv3h
• struct ltl::fderiv4h
• struct ltl::fderiv1h2
• struct ltl::fderiv2h2
• struct ltl::fderiv3h2
• struct ltl::fderiv4h2
• struct ltl::Laplacian2D
• struct ltl::Laplacian2D4
• struct ltl::Laplacian3D
• struct ltl::Laplacian3D4
• struct ltl::grad2D
• struct ltl::kernel_return_type< Iter, grad2D >
• struct ltl::grad3D
• struct ltl::kernel_return_type< Iter, grad3D >

Files

• file convolve.h

---

Detailed Description

Convolution support and kernels in MArray expressions

Convolution kernels for MArrays implementing various smoothing kernels and finite differences methods that can be used in any MArray expression. Currently, convoltions are implemented for expressions of up to rank 3.

Convolution operations are used by calling

   convolve( Expr, Kernel )

within any MArray expression. Expr itself may be any valid MArray expression (partial reductions very likely do not produce correct results right now, though).

Declaring and implementing a convolution kernel.

As an example of a declaration of a convolution kernel, let us consider the o(h^2) first derivative along the dimension dim using central differences:

DECLARE_CONVOLUTION_OPERATOR_DIM(cderiv1,A)
    return OFFSET1(A,1,dim) - OFFSET1(A,-1,dim);
END_CONVOLUTION_OPERATOR;

The macros OFFSET0(), OFFSET(i,dim), OFFSET1(i), OFFSET2(i,j), OFFSET3(i,j,k) allow access to elements of an expression relative to the current position. The offsets can be along any one particular dimension, or along dimensions one, two, and/or three (currently, convolutions are only supported in expression up to rank 3). OFFSET0() is short for the current element (element at offset 0).

The extent of the convolution kernel (the maximum offsets used in the kernel) are determined automatically for you before the operation is executed. The extent can be different in each dimension and in each direction within a dimension. The convolution will only be evaluated for elements of the operand expression for which the kernel lies completely within the domain of the expression. The elements at the boundary will not be accessed.

This kernel can then be used in expressions:

    MArray<float,2> B, C;
    ...
    B = 1/(2*h) * convolve( C*C/2, cderiv1(1) );

Which would assign the first derivative along the first dimension of C*C/2 to B. Note the parameter in the constructor of the kernel cderiv1(1). This is the dimension along which we want to take the derivative.

The DECLARE_CONVOLUTION_OPERATOR_DIM macro creates a kernel with one integer member, the one storing the dimension that is given in the constructor call above. If no argument is needed, the DECLARE_CONVOLUTION_OPERATOR version of the macro should be used. In that case, nothing is stored within the kernel struct.

Various finite-differences kernels are provided by ltl out of the box (first to fourth derivatives based on central, forward, or backward differences, Laplacians, ...). For a complete list and some more documentation, see the file convolve.h.

Note that convolutions can provide higher-order tensors than their operands. An example is the gradient operator (there are no convenience macros yet, so definition has to be provided by hand:

struct grad3D
{
      template<typename Iter>
      inline FVector<typename Iter::value_type,3> eval(Iter& A) const
      {
         FVector<typename Iter::value_type,3> g;
         g(1) = OFFSET(A,1,1) - OFFSET(A,-1,1);
         g(2) = OFFSET(A,1,2) - OFFSET(A,-1,2);
         g(3) = OFFSET(A,1,3) - OFFSET(A,-1,3);
         return g;
      }
};
template <typename Iter>
struct kernel_return_type<Iter,grad3D>
{
   typedef FVector<typename Iter::value_type,3>  value_type;
};

Note that the kernel_return_type trait has to be in namespace ltl.

With this definition of the 3-D gradient operator we can write

MArray<float,2> A(10,10);
MArray<FVector<float,2>,2> B(10,10);
A = 0.0f;
A(5,5) = 1.0f;
B = convolve(A,grad2D());

As a more realistic example for the use of these kernels, consider solving a single timestep in the acoustic wave propagation equation in 3D. Let P_tm1, P_t, and P_tp1 be the pressure field at times tm1=t-1, t, and tp1=t+1, and c be the speed of sound.

Without convolution expressions, this would be implemented like this:

   Range X(2,N-1), Y(2,N-1), Z(2,N-1);

   P_tp1(X,Y,Z) = (2-6*c(X,Y,Z)) * P_t(X,Y,Z)
                + c(X,Y,Z)*(P_t(X-1,Y,Z) + P_t(X+1,Y,Z) + P_t(X,Y-1,Z) + P_t(X,Y+1,Z) + P_t(X,Y,Z-1) + P_t(X,Y,Z+1))
                - P_tm1(X,Y,Z);

With convolution expressions, this becomes:

  P_tp1 = 2 * P_t + c * convolve(P_t, Laplacian3D) - P_tm1;

Of course, standard convolution kernels can be implemented in N-Dimensions as well. As an example, here's a 1-D Gaussian:

struct GaussianKernel1D
{
      GaussianKernel1D (const double s, const double k) :
            sigma_ (s), extent_ (k)
      {
         norm_ = 0.0;
         for (int i = -extent_; i<=extent_; ++i)
            norm_ += exp (-0.5*(double (i*i))/(sigma_*sigma_));
      }

      template<typename Iter>
      inline typename Iter::value_type eval (Iter& A) const
      {
         double r = 0.0;
         for (int i = -extent_; i<=extent_; ++i)
            r += exp (-0.5*(double (i*i))/(sigma_*sigma_))*OFFSET1(A, i);

         return typename Iter::value_type (r/norm_);
      }

      double sigma_, norm_;
      int extent_;
};

template <typename Iter>
struct kernel_return_type<Iter,GaussianKernel1D>
{
   typedef typename Iter::value_type  value_type;
};