File Exchange

image thumbnail

Efficient three-dimensional (3D) Gaussian smoothing using convolution via frequency domain

version 1.2 (2.87 KB) by

Native Fourier implementation, support GPU computation and anisotropic voxel.

5 Downloads

Updated

View License

 Three-dimensional Gaussian smoothing in the frequency domain, native
 frequency domain implementation. Smoothing is achieved by replacing
 the spatial domain convolution with Fourier coefficient
 multiplication. This is also an excellent example to implement your
 own filters using native Fourier expression.

 R = gauss3filter(I);
 R = gauss3filter(I, sigma);
 R = gauss3filter(I, sigma, pixelspacing);
 In a spatial domain representation, R = convn(I, f(x,y,z));
 The Gaussian kernel f(x,y,z) is different depending on the function
 inputs, see the description below. No image padding is provided, pay
 attention to the Fourier wrap-around artifacts.

 Anisotropic smoothing is partly supported, anisotropic voxel size is
 fully supported. Suband_1.5 frequency oversampling is employed to reduce
 numerical erros when sigma is less than the voxel length. Please refer to
 following paper for the Subband_x frequency oversampling technique:
   Max W. K. Law and Albert C. S. Chung, "Efficient Implementation for Spherical Flux Computation and Its Application to Vascular Segmentation",
   IEEE Transactions on Image Processing, 2009, Volume 18(3), 596–612

   R = gauss3filter(I);
   Smooth the image using isotropic smoothing with sigma = 1 voxel-length,
       f(x,y,z) = (2*pi)^(-3/2) * exp(-(x.^2/2 - y.^2/2 - z.^2/2));

   R = gauss3filter(I, sigma);
   If sigma is a scalar, it smooths the image using isotropic smoothing with
   sigma voxel-length,
       f(x,y,z) = (2*pi)^(-3/2)/(sigma^3) * exp(-(x.^2/sigma^2/2 - y.^2/sigma^2/2 - z.^2/sigma^2/2));
   If sigma is a 3D vector, i.e. sigma = [sigma_x sigma_y sigma_z], it
   smooths the image using anisotropic smoothing (oriented anisotropic
   Gaussian is not supported),
       f(x,y,z) = (2*pi)^(-3/2)/sigma(1)/sigma(2)/sigma(3) * exp(-(x.^2/sigma(1)^2/2 - y.^2/sigma(2)^2/2 - z.^2/sigma(3)^2/2));

   R = gauss3filter(I, sigma, pixelspacing);
   If sigma is a scalar, smooth the image using isotropic smoothing with
   sigma physical-length. pixelspacing is a 3D vector. It defines the size
   of a voxel in physical-length,
       f(x,y,z) = (2*pi)^(-3/2)/(sigma^3) * exp(-((x*pixelspacing(1)).^2/sigma^2/2 - (y*pixelspacing(2)).^2/sigma^2/2 - (z*pixelspacing(3)).^2/sigma^2/2));
   If sigma is a 3D vector, sigma = [sigma_x sigma_y sigma_z],
       f(x,y,z) = (2*pi)^(-3/2)/sigma(1)/sigma(2)/sigma(3) * exp(-((x*pixelspacing(1)).^2/sigma(1)^2/2 - (y*pixelspacing(2)).^2/sigma(2)^2/2 - (z*pixelspacing(3)).^2/sigma(3)^2/2));

 Remarks
   The outputs of gauss3filter(I), gauss3filter(I, 1) and
   gauss3filter(I, 1, [1 1 1]) are identical.
 
   To enable GPU computation (Matlab 2012a or later, CUDA 1.3 GPU are required), use
   R = gauss3filter(gpuArray(I), sigma, pixelspacing).

   The Gaussian kernel in the frequency domain is
   exp(-2*pi*pi* (u.^2 *sigma1 + v.^2 *sigma2 + w.^2 * sigma3));

   Please kindly cite the following paper if you use this program, or any code
   extended from this program.
       Max W. K. Law and Albert C. S. Chung, "Efficient Implementation for Spherical Flux Computation and Its Application to Vascular Segmentation”,
       IEEE Transactions on Image Processing, 2009, Volume 18(3), 596–612

 Author: Max W.K. Law
 Email: max.w.k.law@gmail.com
 Page: http://www.cse.ust.hk/~maxlawwk/

Comments and Ratings (1)

UWO

UWO (view profile)

Updates

1.2

Minor help info. update

1.1

Minor update on the file description

MATLAB Release
MATLAB 7 (R14)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video