A Survey of Gaussian Convolution Algorithms

An accurate recursive filter approximation, computed inplace. More...
An accurate recursive filter approximation, computed inplace.
This code implements the recursive filter approximation of Gaussian convolution proposed by Vliet, Young, and Verbeek. The Gaussian is approximated by a cascade of a causal filter and an anticausal filter,
The process to use these functions is the following:
Data Structures  
struct  vyv_coeffs_ 
Coefficients for VlietYoungVerbeek Gaussian approximation. More...  
Macros  
#define  VYV_MIN_K 3 
Minimum valid VYV filter order.  
#define  VYV_MAX_K 5 
Maximum valid VYV filter order.  
#define  VYV_VALID_K(K) (VYV_MIN_K <= (K) && (K) <= VYV_MAX_K) 
Test whether a given K value is a valid VYV filter order.  
Typedefs  
typedef struct vyv_coeffs_  vyv_coeffs 
Coefficients for VlietYoungVerbeek Gaussian approximation. More...  
Functions  
static double  variance (const complex *poles0, int K, double q) 
Compute the variance of the impulse response. More...  
static double  dq_variance (const complex *poles0, int K, double q) 
Derivative of variance with respect to q. More...  
static double  compute_q (const complex *poles0, int K, double sigma, double q0) 
Compute q for a desired sigma using Newton's method. More...  
static void  expand_pole_product (double *c, const complex *poles, int K) 
Expand pole product. More...  
void  vyv_precomp (vyv_coeffs *c, double sigma, int K, num tol) 
Precomputations for VlietYoungVerbeek Gaussian approximation. More...  
void  vyv_gaussian_conv (vyv_coeffs c, num *dest, const num *src, long N, long stride) 
Gaussian convolution VlietYoungVerbeek approximation. More...  
void  vyv_gaussian_conv_image (vyv_coeffs c, num *dest, const num *src, int width, int height, int num_channels) 
2D Gaussian convolution VlietYoungVerbeek approximation More...  
typedef struct vyv_coeffs_ vyv_coeffs 
Coefficients for VlietYoungVerbeek Gaussian approximation.
The vyv_coeffs struct stores the coefficients for the recursive filter of order K. This struct allows to precompute these filter coefficients separately from actually performing the filtering so that filtering may be performed multiple times using the same precomputed coefficients.
This coefficients struct is precomputed by vyv_precomp() and then used by vyv_gaussian_conv() or vyv_gaussian_conv_image().

static 
Compute q for a desired sigma using Newton's method.
poles0  unscaled pole locations 
K  number of poles 
sigma  the desired sigma 
q0  initial estimate of q 
This routine uses Newton's method to solve for the value of q so that the filter achieves the specified variance,
where the are the unscaled pole locations.
Definition at line 99 of file gaussian_conv_vyv.c.

static 
Derivative of variance with respect to q.
poles0  unscaled pole locations 
q  rescaling parameter 
K  number of poles 
This function is used by compute_q() in solving for q.
Definition at line 66 of file gaussian_conv_vyv.c.

static 
Expand pole product.
c  resulting filter coefficients 
poles  pole locations 
K  number of poles 
This routine expands the product to obtain the filter coefficients:
Definition at line 124 of file gaussian_conv_vyv.c.

static 
Compute the variance of the impulse response.
poles0  unscaled pole locations 
q  rescaling parameter 
K  number of poles 
Definition at line 40 of file gaussian_conv_vyv.c.
void vyv_gaussian_conv  (  vyv_coeffs  c, 
num *  dest,  
const num *  src,  
long  N,  
long  stride  
) 
Gaussian convolution VlietYoungVerbeek approximation.
c  vyv_coeffs created by vyv_precomp() 
dest  output convolved data 
src  data to be convolved, modified inplace if src = dest 
N  number of samples 
stride  stride between successive samples 
This routine performs VlietYoungVerbeek recursive filtering Gaussian convolution approximation. The Gaussian is approximated by a causal filter followed by an anticausal filter,
The array c.filter
holds the filter coefficients , which are precomputed by vyv_precomp().
The convolution can be performed inplace by setting src
= dest
(the source array is overwritten with the result).
Definition at line 259 of file gaussian_conv_vyv.c.
void vyv_gaussian_conv_image  (  vyv_coeffs  c, 
num *  dest,  
const num *  src,  
int  width,  
int  height,  
int  num_channels  
) 
2D Gaussian convolution VlietYoungVerbeek approximation
c  vyv_coeffs created by vyv_precomp() 
dest  destination image data 
src  source image data, modified inplace if src = dest 
width  image width 
height  image height 
num_channels  number of image channels 
Similar to vyv_gaussian_conv(), this routine approximates 2D Gaussian convolution with VlietYoungVerbeek recursive filtering.
The convolution can be performed inplace by setting src
= dest
(the source array is overwritten with the result).
Definition at line 397 of file gaussian_conv_vyv.c.
void vyv_precomp  (  vyv_coeffs *  c, 
double  sigma,  
int  K,  
num  tol  
) 
Precomputations for VlietYoungVerbeek Gaussian approximation.
c  coefficients structure to populate 
sigma  the standard deviation of the Gaussian in pixels 
K  filter order = 3, 4, or 5 
tol  accuracy for initialization on the left boundary 
This routine precomputes the filter coefficients for VlietYoungVerbeek approximate Gaussian convolution for use in vyv_gaussian_conv() or vyv_gaussian_conv_image(). The filter coefficients are computed by the following steps:
where the are the unscaled pole locations.
For handling the right boundary, the routine precomputes the inverse to the linear system
The inverse is stored in matrix c>M
, ordered such that
Definition at line 181 of file gaussian_conv_vyv.c.