A Survey of Gaussian Convolution Algorithms
gaussian_conv_ebox.c
Go to the documentation of this file.
1 
20 #include "gaussian_conv_ebox.h"
21 #include <assert.h>
22 #include <math.h>
23 #include <string.h>
24 #include "filter_util.h"
25 
35 void ebox_precomp(ebox_coeffs *c, double sigma, int K)
36 {
37  double alpha;
38 
39  assert(c && sigma > 0 && K > 0);
40 
41  /* Set the parameters according to Gwosdek et al. */
42  c->r = (long)(0.5 * sqrt((12.0 * sigma * sigma) / K + 1.0) - 0.5);
43  alpha = (2 * c->r + 1) * (c->r * (c->r + 1) - 3.0 * sigma * sigma / K)
44  / (6.0 * (sigma * sigma / K - (c->r + 1) * (c->r + 1)));
45  c->c_1 = alpha / (2.0 * (alpha + c->r) + 1);
46  c->c_2 = (1.0 - alpha) / (2.0 * (alpha + c->r) + 1);
47  c->K = K;
48  return;
49 }
50 
69 static void ebox_filter(num *dest, long dest_stride,
70  const num *src, long src_stride, long N, long r, num c_1, num c_2)
71 {
72  long n;
73  num accum = 0;
74 
75  assert(dest && src && dest != src && N > 0 && r >= 0);
76 
77  for (n = -r; n <= r; ++n)
78  accum += src[src_stride * extension(N, n)];
79 
80  dest[0] = accum = c_1 * (src[src_stride * extension(N, r + 1)]
81  + src[src_stride * extension(N, -r - 1)])
82  + (c_1 + c_2) * accum;
83 
84  for (n = 1; n < N; ++n)
85  {
86  accum += c_1 * (src[src_stride * extension(N, n + r + 1)]
87  - src[src_stride * extension(N, n - r - 2)])
88  + c_2 * (src[src_stride * extension(N, n + r)]
89  - src[src_stride * extension(N, n - r - 1)]);
90  dest[dest_stride * n] = accum;
91  }
92 
93  return;
94 }
95 
117 void ebox_gaussian_conv(ebox_coeffs c, num *dest_data, num *buffer_data,
118  const num *src, long N, long stride)
119 {
120  struct
121  {
122  num *data;
123  long stride;
124  } dest, buffer, cur, next;
125  int step;
126 
127  assert(dest_data && buffer_data && src
128  && dest_data != buffer_data && N > 0);
129 
130  dest.data = dest_data;
131  dest.stride = stride;
132  buffer.data = buffer_data;
133  buffer.stride = (buffer_data == src) ? stride : 1;
134 
135  next = (buffer_data == src || (dest_data != src && c.K % 2 == 1))
136  ? dest : buffer;
137  /* Perform the first filtering pass. */
138  ebox_filter(next.data, next.stride, src, stride,
139  N, c.r, c.c_1, c.c_2);
140 
141  /* Perform (K - 1) filtering passes, alternating the roles of the dest
142  and buffer arrays. */
143  for (step = 1; step < c.K; ++step)
144  {
145  cur = next;
146  next = (cur.data == buffer_data) ? dest : buffer;
147  ebox_filter(next.data, next.stride, cur.data, cur.stride,
148  N, c.r, c.c_1, c.c_2);
149  }
150 
151  /* If necessary, copy the result to the destination array. */
152  if (next.data != dest_data)
153  {
154  if (stride == 1)
155  memcpy(dest_data, buffer_data, sizeof(num) * N);
156  else
157  {
158  long n, i;
159 
160  for (n = i = 0; n < N; ++n, i += stride)
161  dest_data[i] = buffer_data[n];
162  }
163  }
164 
165  return;
166 }
167 
186  const num *src, int width, int height, int num_channels)
187 {
188  const long num_pixels = ((long)width) * ((long)height);
189  int x, y, channel;
190 
191  assert(dest && buffer && src && dest != buffer && num_pixels > 0);
192 
193  /* Loop over the image channels. */
194  for (channel = 0; channel < num_channels; ++channel)
195  {
196  num *dest_y = dest;
197  const num *src_y = src;
198 
199  /* Filter each row of the channel. */
200  for (y = 0; y < height; ++y)
201  {
202  ebox_gaussian_conv(c, dest_y, buffer, src_y, width, 1);
203  dest_y += width;
204  src_y += width;
205  }
206 
207  /* Filter each column of the channel. */
208  for (x = 0; x < width; ++x)
209  ebox_gaussian_conv(c, dest + x, buffer, dest + x, height, width);
210 
211  dest += num_pixels;
212  src += num_pixels;
213  }
214 
215  return;
216 }