5.0

5.0 | 1 rating Rate this file 37 Downloads (last 30 days) File Size: 2.87 KB File ID: #39013

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

by

 

12 Nov 2012 (Updated )

Native Fourier implementation, support GPU computation and anisotropic voxel.

| Watch this File

File Information
Description

 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/

Required Products MATLAB
MATLAB release MATLAB 7 (R14)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (1)
12 Mar 2013 UWO  
Updates
14 Nov 2012

Minor update on the file description

16 Nov 2012

Minor help info. update

Contact us