A Survey of Gaussian Convolution Algorithms
gaussian_conv_box.c
Go to the documentation of this file.
1 
20 #include "gaussian_conv_box.h"
21 #include <assert.h>
22 #include <math.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include "filter_util.h"
26 
43 static void box_filter(num *dest, long dest_stride,
44  const num *src, long src_stride, long N, long r)
45 {
46  long n;
47  num accum = 0;
48 
49  assert(dest && src && dest != src && N > 0 && r >= 0);
50 
51  /* Initialize the filter on the left boundary by directly computing
52  dest(0) = accum = \sum_{n=-r}^r src(n). */
53  for (n = -r; n <= r; ++n)
54  accum += src[src_stride * extension(N, n)];
55 
56  dest[0] = accum;
57 
58  /* Filter the interior samples. */
59  for (n = 1; n < N; ++n)
60  {
61  /* Update accum: add sample src(n + r) and remove src(n - r - 1). */
62  accum += src[src_stride * extension(N, n + r)]
63  - src[src_stride * extension(N, n - r - 1)];
64  dest[dest_stride * n] = accum;
65  }
66 
67  return;
68 }
69 
105 void box_gaussian_conv(num *dest_data, num *buffer_data, const num *src,
106  long N, long stride, num sigma, int K)
107 {
108  struct
109  {
110  num *data;
111  long stride;
112  } dest, buffer, cur, next;
113  num scale;
114  long r;
115  int step;
116 
117  assert(dest_data && buffer_data && src && dest_data != buffer_data
118  && N > 0 && sigma > 0 && K > 0);
119 
120  /* Compute the box radius according to Wells' formula. */
121  r = (long)(0.5 * sqrt((12.0 * sigma * sigma) / K + 1.0));
122  scale = (num)(1.0 / pow(2*r + 1, K));
123 
124  dest.data = dest_data;
125  dest.stride = stride;
126  buffer.data = buffer_data;
127  buffer.stride = (buffer_data == src) ? stride : 1;
128 
129  /* Here we decide whether dest or buffer should be the first output array.
130  If K is even, then buffer is the better choice so that the result is in
131  dest after K iterations, e.g. for K = 4 (as in the function comment),
132 
133  src -> buffer -> dest -> buffer -> dest.
134 
135  If K is odd, we would like to choose dest, e.g. for K = 3,
136 
137  src -> dest -> buffer -> dest.
138 
139  However, if src and dest point to the same memory (i.e., in-place
140  computation), then we must select buffer as the first output array. */
141  if (buffer_data == src || (dest_data != src && K % 2 == 1))
142  next = dest;
143  else
144  next = buffer;
145 
146  /* Perform the first step of box filtering. */
147  box_filter(next.data, next.stride, src, stride, N, r);
148 
149  /* Perform another (K - 1) steps of box filtering, alternating the roles
150  of the dest and buffer arrays. */
151  for (step = 1; step < K; ++step)
152  {
153  cur = next;
154  next = (cur.data == buffer_data) ? dest : buffer;
155  box_filter(next.data, next.stride, cur.data, cur.stride, N, r);
156  }
157 
158  /* Multiply by the constant scale factor. */
159  if (next.data != dest_data)
160  {
161  long n, i;
162 
163  for (n = i = 0; n < N; ++n, i += stride)
164  dest_data[i] = buffer_data[n] * scale;
165  }
166  else
167  {
168  long i, i_end = stride * N;
169 
170  for (i = 0; i < i_end; i += stride)
171  dest_data[i] *= scale;
172  }
173 
174  return;
175 }
176 
195 void box_gaussian_conv_image(num *dest, num *buffer, const num *src,
196  int width, int height, int num_channels, num sigma, int K)
197 {
198  const long num_pixels = ((long)width) * ((long)height);
199  int x, y, channel;
200 
201  assert(dest && buffer && src && dest != buffer
202  && num_pixels > 0 && sigma > 0 && K > 0);
203 
204  /* Loop over the image channels. */
205  for (channel = 0; channel < num_channels; ++channel)
206  {
207  num *dest_y = dest;
208  const num *src_y = src;
209 
210  /* Filter each column of the channel. */
211  for (y = 0; y < height; ++y)
212  {
213  box_gaussian_conv(dest_y, buffer, src_y,
214  width, 1, sigma, K);
215  dest_y += width;
216  src_y += width;
217  }
218 
219  /* Filter each row of the channel. */
220  for (x = 0; x < width; ++x)
221  box_gaussian_conv(dest + x, buffer, dest + x,
222  height, width, sigma, K);
223 
224  dest += num_pixels;
225  src += num_pixels;
226  }
227 
228  return;
229 }