[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

Image Processing VIGRA

Section Contents

In this chapter we'll use VIGRA's methods for some applications of Image Processing.

Calling Conventions

VIGRA's image processing functions follow a uniform calling convention: The argument list start with the input images or arrays, followed by the output images or arrays, followed by the function's parameters (if any). Some functions additionally accept an option object that allows more fine-grained control of the function's actions and must be passed as the last argument. Most functions assume that the output arrays already have the appropriate shape.

All functions working on arrays expect their arguments to be passed as vigra::MultiArrayView instances. Functions that only support 2-dimensional images usually contain the term "Image" in their name, whereas functions that act on arbitrary many dimensions usually contain the term "Multi" in their name.
Examples:

// determine the connected components in a binary image, using the 8-neighborhood
MultiArray<2, UInt8> image(width, height);
MultiArray<2, UInt32> labels(width, height);
... // fill image
labelImage(image, labels, true);
// smooth a 3D array with a gaussian filter with sigma=2.0
MultiArray<3, float> volume(Shape3(300, 200, 100)),
smoothed(Shape3(300, 200, 100));
... // fill volume
gaussianSmoothMultiArray(volume, smoothed, 2.0);
// compute the determinant of a 5x5 matrix
MultiArray<2, float> matrix(Shape2(5, 5));
... // fill matrix with data
float det = linalg::determinant(matrix);

For historical reasons, VIGRA also supports two alternative APIs in terms of iterators. These APIs used to be considerably faster, but meanwhile compilers and processors have improved to the point where the much simpler MultiArrayView API no longer imposes a significant abstraction penalty. While there are no plans to remove support for the old APIs, they should not be used in new code.

Inverting an Image

Inverting an (gray scale) image is quite easy. We just need to subtract every pixel's value from white (255). This simple task doesn't need an explicit function call at all, but is best solved with a arithmetic expression implemented in namespace vigra::multi_math. To avoid possible overload ambiguities, you must explicitly activate array arithmetic via the command using namespace vigra::multi_math before use. To invert imageArray and overwrite its original contents, you write:

using namespace vigra::multi_math;
imageArray = 255-imageArray;

See here for a complete example: invert_tutorial.cxx

This is the result:

lenna_small.gif
input file
lenna_inverted.gif
inverted output file

Image Blending

In this example, we have two input images and want to blend them into one another. In the combined output image every pixel value is the mean of the two appropriate original pixels. This is also best solved with array arithmetic:

using namespace vigra::multi_math;
exportArray = 0.5*imageArray1 + 0.5*imageArray2;

Since it is not guaranteed that the two input images have the same shape, we first determine the maximum possible shape of the blended image, which equals the minimum size along each axis. With the help of subarray-method we just blend the appropriate parts of the two images. These parts (subimages) are aligned around the centers of the original images.

Here's the code: dissolve.cxx

And here are the results:

lenna_color_small.gif
input file 1
oi_small.jpg
input file 2
dissolved_color.gif
dissolved output file

Creating a Composite Image

Let's come to a little gimmick. Given one input image we want to create a composite image of 4 images reflected with respect to each other. The result resembles the effect of a kaleidoscope. Two of VIGRA's functions are sufficient for this purpose: subarray(p,q) and reflectImage(). Input and output images of reflectImage() are specified by MultiArrayViews. The third parameter specifies the desired reflection axis. The axis can either be horizontal, vertical or both (as in this example):

reflectImage(inputArray, outputArray, horizontal | vertical);

Here's the code: composite.cxx

And here are the results:

lenna_color_small.gif
input file
lenna_composite_color.gif
composite output file

Smoothing

2-dimensional Convolution

There are many different ways to smooth an image. Before we use VIGRA's methods, we want to write a smoothing code of our own. The idea is to choose each pixel in turn and replace it with the mean of itself and the pixels in 5x5 window around it. To calculate the mean in a window, we can just devide the sum of the pixel values within the corresponding subarray by their number. MultiArrayView provides two useful methods for doing this: sum and size. In our code we iterate over every pixel, construct the surrounding 5x5 window via subarray, and write the average of the window into the corresponding output pixel. Near the borders of the image we truncate the window appropriately so that it remains inside the image, and only take the average over the actually existing neighbours of the pixel.

See the code: smooth_explicitly.cxx

The results:

lenna_small.gif
input file
lenna_smoothed.gif
smoothed output file

The technical term for this kind of operation is convolution. VIGRA provides convolveImage as a comfortable way to perform 2-dimensional convolutions with arbitrary filters. You may use it as follows:

convolveImage(inputImage, resultImage, filter);

The filter of convolution kernel is given as argument object by kernel2d(). To implement the above smoothing by taking averages in 3x3 windows, you need an averaging kernel with radius 1. Kernel truncation near the image borders is performed when the filter's border treatment mode is set to BORDER_TREATMENT_CLIP:

Kernel2D<double> filter;
filter.initAveraging(1);
filter.setBorderTreatment(BORDER_TREATMENT_CLIP);

By default, VIGRA's convolution functions use BORDER_TREATMENT_REFLECT (i.e. the image is virtually enlarged by reflecting the pixel values about the border), which usually leads to superior results. The strength of smoothing can be controlled by increasing the filter radius.

Another improvement over simple averaging can be achieved when one takes a weighted average such that pixels near the center have more influence on the result. A popular choice here is the 5x5 binomial filter. VIGRA allows to specify arbitrary filter shapes and coefficients via the Kernel2D::initExplicitly():

Kernel2D<float> filter;
// specify filter shape (lower right corner is inclusive here!)
filter.initExplicitly(Shape2(-2,-2), Shape2(2,2));
// specify filter coefficients
filter = 1.0/256.0, 4.0/256.0, 6.0/256.0, 4.0/256.0, 1.0/256.0,
4.0/256.0, 16.0/256.0, 24.0/256.0, 16.0/256.0, 4.0/256.0,
6.0/256.0, 24.0/256.0, 36.0/256.0, 24.0/256.0, 6.0/256.0,
4.0/256.0, 16.0/256.0, 24.0/256.0, 16.0/256.0, 4.0/256.0,
1.0/256.0, 4.0/256.0, 6.0/256.0, 4.0/256.0, 1.0/256.0;
// apply filter
convolveImage(inputImage, resultImage, filter);

initExplicitly() receives the upper left and lower right corners of the filter window. Note that the lower right corner here is included in the window, in contrast to MultiArray::subarray() where the end point is not included.

The filter weights are provided in a comma separated list. Normally, the sum of the coefficients should to be 1 in order to preserve the average intensity of the image. You must provide either as many coefficients as needed for the given filter size, or exactly one value which will be used for all filter coefficients. Thus, the 3x3 averaging filter can also be created like this:

Kernel2D<double> filter;
filter.initExplicitly(Shape2(-1,-1), Shape2(1,1)) = 1.0/9.0;

For various theoretical and practical reasons, the Gaussian filter is the best choice in most situations. Its coefficients are chosen according to a Gaussian (i.e. bell-shaped) function with given standard deviation. The kernel class has a convenient initGaussian(std_dev) method that creates the appropriate coefficients:

filter.initGaussian(1.5);
convolveImage(inputImage, resultImage, filter);

A complete example using these possibilities can be found in smooth_convolve.cxx.


Separable Convolution in 2D and nD Images

When filtering is implemented with 2-dimensional windows as in the previous section, we need as many multiplications per pixel as there are coefficients in the filter. Fortunately, many important filters (including averaging and Gaussian smoothing) have the property of beeing separable, which allows a much more efficient implementation in terms of 1-dimensional windows. A 2-dimensional filter is separable if its coefficients $f_{ij}$ can be expressend as an outer product of two 1-dimensional filters $h_i$ and $c_j$:

\[ f_{ij} = h_i \cdot c_j \]

For example, the 3x3 averaging filter (with coefficients 1/9) is obtained as the outer product of two 3x1 filters (with coefficients 1/3):

\[ \left( \begin{array}{ccc} \frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\[1ex] \frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\[1ex] \frac{1}{9} & \frac{1}{9} & \frac{1}{9} \end{array} \right) = \left( \begin{array}{c} \frac{1}{3} \\[1ex] \frac{1}{3} \\[1ex] \frac{1}{3} \end{array} \right) \cdot \left( \begin{array}{ccc} \frac{1}{3} & \frac{1}{3} & \frac{1}{3} \end{array} \right) \]

The convolution with separable filters can be implemented by two consecutive 1-dimensional convolutions: first, one filters all rows of the image with the horizontal filter, and then all columns of the result with the vertical filter. Instead of the (n x m) operations required for a 2-dimensional window, we now only need (n + m) operations for the two 1-dimensional ones. Already for a 5x5 window, this reduces the number of operations from 25 to 10, and the difference becomes even bigger with increasing window size.

To construct and apply 1-dimensional filters, VIGRA provides the class vigra::Kernel1D and the functions separableConvolveX() resp. separableConvolveY(). To compute a 2D Gaussian filter we use the following code:

Kernel1D<double> filter;
filter.initGaussian(1.5);
MultiArray<2, float> tmpImage(inputImage.shape());
separateConvolveX(inputImage, tmpImage, filter);
separateConvolveY(tmpImage, resultImage, filter);

Note that we need an intermediate image to hold the result of the horizontal filtering. The same result is more conveniently achieved by the functions convolveImage() and gaussianSmoothing() (see smooth_convolve.cxx for a working example):

// apply 'filter' to both the x- and y-axis
// (calls separateConvolveX() and separateConvolveY() internally)
convolveImage(inputImage, resultImage, filter, filter);
// smooth image with Gaussian filter with sigma=1.5
// (calls convolveImage() with Gaussian filter internally)
gaussianSmoothing(inputImage, resultImage, 1.5);

It is, of course, also possible to apply different filters in the x- and y-directions. This is especially useful for derivative filters which are commonly used to compute image features, for example gaussianGradient() and gaussianGradientMagnitude(). For more information see Convolution Filters.

Separable filters are also the key for efficient convolution of higher-dimensional images and arrays: An n-dimensional filter is simply implemented by n consecutive 1-dimensional filter applications, regardless of the size of n. This is the basis for VIGRA's multi-dimensional filter functions. For example, Gaussian smoothing in arbitrary many dimensions is implemented in gaussianSmoothMultiArray():

MultiArray<3, UInt8> inputArray(Shape3(100, 100, 100));
... // fill inputArray with data
MultiArray<3, float> resultArray(inputArray.shape());
// perform isotropic Gaussian smoothing at scale 1.5
gaussianSmoothMultiArray(inputArray, resultArray, 1.5);

More information about VIGRA's multi-dimensional convolution funcions can be found in the reference manual under Convolution Filters.

Parallel Execution of Gaussian Filters

The computation of Gaussian filters and their derivatives can be accelerated significantly when rectangular blocks of a large image as processed in parallel. This is easily achieved in VIGRA by passing the option object vigra::BlockwiseConvolutionOptions to the convolution functions:

// create a big array
MultiArray<3, UInt8> inputArray(Shape3(1000, 1000, 100));
... // fill inputArray with data
MultiArray<3, float> resultArray(inputArray.shape());
// perform isotropic Gaussian smoothing at scale 1.5 in parallel
gaussianSmoothMultiArray(inputArray, resultArray, 1.5,
BlockwiseConvolutionOptions<3>());

This call will spawn the standard number of threads for the present platform (as returned by std::thread::hardware_concurrency()) and distributes the work across these threads in blocks with a suitable default shape. You can customize the number of threads and the block shape via the option object:

gaussianSmoothMultiArray(inputArray, resultArray, 1.5,
BlockwiseConvolutionOptions<3>().numThreads(6)
.blockShape(Shape3(128, 128, 100)));

The same works for Gaussian derivative filters such as gaussianGradientMultiArray(), gaussianGradientMagnitude(), and hessianOfGaussianMultiArray(). Refer to section Convolution Filters for more details.

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1 (Fri May 19 2017)