This is a presentation of the work that we have done in the field of 1D, 2D, 3D and 4D spatial filtering. We have focused on the convolution with Gaussian filter. The results were achieved with our in-house developed LinearFilters library .

Motivation

The motivation for fast and accurate computation of Gaussian filtering was to be able to compute certain optical flow algorithms efficiently. Such group of algorithms employ several Gabor filters (filtering banks) what can be turned into computation of many convolutions with Gaussian filters [1]. The filtering has always been the major slowdown factor of the approach.

The overall goal is to equip such methods with fast and accurate implementation of algorithms which should be able to compute truly arbitrary, i.e. anisotropic with generaly-oriented main axis, Gaussian (and Gabor as well) filtering for 3D and time lapse 3D images. The LinearFilters library should serve this purpose.

Introduction

Naturally, the first improvement in terms of time consumption was to switch from the finite-impulse-response (FIR) Gaussian into the infinite-impulse-response (IIR) version [2,3]. Appart from a huge speed gain, this switch introduced a decent decrease in the accuracy of the filter as demonstrated in the results section.

Nevertheless, the popularity in the literature as well as the number of direct advantages of the IIR approach outweighed. These are:

• simple filter design with fixed number of coefficients irrelevant to Sigma size,
  
• memory efficiency with no temporary images required,
• improved time consumption which does not vary with Sigma size.
  

IIR Gauss filter sweeps the filtered image twice in every direction. I(x) and O(x) are the pixel values at coordinate x in the input and output image, respectively. A(S), B(S), C(S) and D(S) are filter coefficients dependent solely on the Sigma S, which is constant during the filtering operation, of the desired Gaussian.

Illustration of memory blocks accessed during the 3D separable convolution. Note that in the FIR scheme a temporary image T(x) is required since input image is always considered read-only.

Time consumption (in seconds) of 1D, 2D and 3D separable convolutions on a sample image. The image was 5000x5000x1 pixels for 1D and 2D tests, 500x500x100 for 3D test, voxel values were float, image size was then 95.4MB. Computer was Intel(R) Pentium(R) 4 at 2.8GHz with 2GB DDR400 RAM on Linux. Note that the FIR filter's support depends on the Sigma, hence the time demand increases.

Another feature of this library, that is required for the oriented anisotropic Gaussian filtering, is the availability of routines that run 1D convolution along user-specified directions. The demanded filter width, the Sigma, needs to be adjusted in this case because the sampling unit changes. The unit length changes to the length of the direction-specifying vector as illustrated in the next image.

The IIR Gaussian along the (2,1) direction, visualized with arrows, in the 2D image. Three applications of filter coefficients on the image is shown with pixels denoted with corresponding letters and indices.

We support the following configurations for the specification of a 1D convolution direction (currently, only 3D directions are supported):

• (x,1,0) with x being real number,
• (x,y,1) with x and y being real numbers,
• (x,y,z) with x,y and z being integers.

The disadvantage of the first two options is the fact that convolution may "fall off" the voxel grid when x or y is not integer. Under such circumstances, the linear, or bilinear, interpolations are used both for "reading" and "writing" values from and to the images. This introduces additional speed penalty as well as some loss of accuracy.

Last but not least, the following paradigms were kept in mind during the programming:

• utilize modern processors' pipeline processing by not putting if-statements inside loops,
• utilize the availability of n-way caches by accessing small memory amounts at small number of different physical addresses,
• disassemble the process of the whole image filtering and compute certain filtering stages in parallel,
• don't use assembler language and don't do low-level optimalizations in the code.

We believe that, by following these requirements, we can achieve a good performace on most of modern PC-based computers.

Methods

Basically, a Gaussian filter in 3D can be defined using the xy-plane rotation angle, tilting angle around newly positioned y-axis and three Sigmas. These input parameters can be embedded into the covariance matrix C which uniquely represents given Gaussian:

The convolution kernel of a desired Gaussian filter given by covariance matrix C. The spatial coordinate is given with x, filter parameters are encoded in the C.

This matrix is central to the two approaches that appeared quite recently. The first one, due to Lampert and Wirjadi [4], computes a cascade of 3 carefully oriented directional 1D convolutions. The directions and appropriate Sigmas are derived from the triangular decomposition of Cholesky type of the matrix C. The second one, due to Lam and Shi [5], benefits of the cascade of 1D convolutions as well. However, this approach increases the number of convolutions to 6 on behalf an additional constraint on their directions. This constraint pushes convolutions to always stay on the voxel grid what vanishes the need for interpolations. The first approach is optimal in terms of the number of memory accesses and interpolation operations required [4]. The second one seems to aproximate Gaussian better while being only up to 1.5times slower, both will be shown in the Results section.

The Lampert and Wirjadi approach for nD: decompose the input covariance matrix C specifically what results in the effective matrix V and convolve along V's columns with corresponding Sigmas.

Extension of Lam and Shi for nD: select m directional vectors v_i (the matrix A), compute Sigmas corresponding to them and convolve along them. The selection of matrix A can be tricky.

The 6-convolutions approach requires to establish 6 directions prior the computation of the convolution. Unfortunately, we have only the implicit definition of the matrix E, i.e. the Gaussian separated along chosen directions, because the transition matrix A has no inverse. Thus, the process of selecting the appropriate directions for a given Gaussian is more or less a process of guessing. In the library, we have pinpointed, based on our experience, a reasonable set of possible directions among to choose from. On the other hand, the availability of more possible matrices A determining Gaussian filters at roughly the same level of quality, when compared to the original filter, makes a room for an interesting concept "sharing common directions". This may prove advantageous when considering a bank of filters.

Bibliography

[1] A. Bernardino and J. Santos-Victor: Fast IIR Isotropic 2-D Complex Gabor Filters With Boundary Initialization. IEEE Transactions on Image Processing vol. 15(11), pp. 3338-3348, 2006.

[2] I. T. Young and L. J. van Vliet: Recursive implementation of the Gaussian filter. Signal processing, vol. 44(2), pp. 139-151, 1995.

[3] I.T. Young and L.J. van Vliet and M. van Ginkel: Recursive Gabor Filtering. Signal processing, vol. 50(11), pp. 2798-2805, 2005.

[4] C.H. Lampert and O. Wirjadi: An Optimal Nonorthogonal Separation of the Anisotropic Gaussian Convolution Filter. IEEE Transactions on Image Processing vol. 15(11), pp. 3501-3513, 2006.

[5] S.Y.M. Lam and B.E. Shi: Recursive anisotropic 2-D Gaussian filtering based on a triple-axis decomposition. IEEE Trans. on Image Processing, vol. 16(7), pp. 1925-30, 2007.

Results

First of all, the accuracy of FIR and IIR implementations is demonstrated for 1D and 2D sample impulse responses. One may notice that there is relatively small difference in the responses for small sigmas what follows the observation in [3].

The comparison of 1D impulse responses between the estimated, FIR and IIR responses.

The comparison of 2D impulse responses between the estimated, FIR and IIR responses.

Another visual inspection is given in the next three images. In these, two narrow threshold intervals were selected and applied on the impulse responses of the true Gaussian and FIR and IIR convolution results. The first approach is rather squared whereas the second one gets rather broader.

Two fixed threshold intervals for sample impulse response of the approach of Lampert and Wirjadi.

Two fixed threshold intervals for sample impulse response of a true Gaussian.

Two fixed threshold intervals for sample impulse response of the extension of the method of Lam and Shi.

The error measure NSE was adopted from [5] for the purpose of quantitative evaluation. NSE gives the average squared magnitude of errors with respect to the total energy of the test image. We filtered a test image with both approaches as well as with true Gaussian over many 3D Gaussian orientations and over 2 selections of Sigmas: Sigma along main axis was 2 and 5, the other two sigmas were both set to 1 and 2, respectively. The error distributions between each approach and the expected result, the true Gaussian, are depicted.

Lampert and Wirjadi's NSE error distribution.

Lam and Shi's NSE error distribution.

Similarily to the previous NSE error distributions, we took averages over all orientations and summarized them in the table. The last three lines are for the test image whereas the rest is for filtered impulse response.

Finally, the time demands were constant for any desired filter configuration. This is a result of the fixed number of convolutions and the fixed size of filter's support (the number of its coefficients) which disregards the size of desired Sigmas.

We achieved 76ms for Lampert's (3 1D convolutions, two with interpolations) and 110ms for Shi's extension (6 1D convolutions, no interpolations) on average. This was for 3D image of 130x130x90px with float representation (stored using 5.8MB of memory) for voxel intensities on Intel(R) Core2(TM) Quad based PC-type computer.

Download

We provide source codes for the entire project . This is a ZIP archive (0.49MB) designed originaly for Linux operating systems. Beside source codes of the i3dcore, i4dcore and LinearFilters libraries, it contains Makefiles for each library and a global Makefile to easily compile everything including a sample program for 2D oriented filtering. Documentation of filtering functions can be found here.

The included i3dcore was tailored to handle only TIFF, JPEG, TGA, METAIO and I3D file formats. It is expected that libtiff, libjpeg and libz be installed in the operating system already. The original i3dcore can also benefit from ICS (via libics) and Bio-Formats -supported image file formats. These are less common software packages probably not installed on your system. Hence, for the sake of simplicity, we've decided not to support them in this release of i3dcore. However, anyone interested in the topic is welcome to contant us for further details on how to enable these.

A binary versions for Linux (ZIP, 2.0MB) and Windows (TM) (ZIP, 3.3MB) can be downloaded too. A simple example on how to compile your programs with these libraries is included.

7th August 2009

The new version of the library has been released.

We no longer provide separate download boundle of just the LinearFilters library (the links in the text above lead to the version from 25th March 2009). The optical flow download page is now the correct and valid source where this library, among others, can be found. The page also summarizes a small how-to compile and install the library.

 25th March 2009

The published version here already implements the improved extension of the approach of Lam and Shi, the function ApplyANIGaussU2(). Also this release includes first implementations of spatial Gabor filtering.

Please, take a note of the licence below before downloading.

Example of use

The following is a sample snapshot of C++ code that computes convolution with a general Gaussian filter using our library:

#include <i3d/image3d.h>
#include <of/ofLinearFilter.h>

using namespace i3d;

int main(void) {
//read input image
Image3d<float> i_orig("3Dcell.tif");
//create output image
Image3d<float> i_res;

//do Gaussian filtering in the direction of vector (1,1,1),
//the sigma in this direction is 4, the other two sigmas are 2
//
//Lampert and Wirjadi
ApplyANIGaussL<float>(i_orig,i_res,DegToRad(45),DegToRad(45),4,2,2);
//
//extension of Lam and Shi 
ApplyANIGaussU2<float>(i_orig,i_res,DegToRad(45),DegToRad(45),4,2,2);

//save output image
i_res.SaveImage("3Dcell-filtered.tif",IMG_TIFF);

return(0);
} 

Acknowledgement

The research on the topic is supported by the scientific project NPVII (2B06052) by the Ministry of Education of the Czech Republic.

Licence

The library is part of the OpticalFlow boundle in which every part is available under the GNU GPL license. Hence, the LinearFilters library is available under the GNU GPL license.

Conference Paper

This work has also been published in the Proceedings of the 8th WSEAS International Conference on Signal Processing, Computational Geometry and Artificial Vision (ISCGAV'08). Click here to download the paper.