A Survey of Gaussian Convolution Algorithms
gaussian_conv_sii.c
Go to the documentation of this file.
1 
20 #include "gaussian_conv_sii.h"
21 #include <assert.h>
22 #include <stdio.h>
23 #include "filter_util.h"
24 
25 #ifndef M_PI
26 
27 #define M_PI 3.14159265358979323846264338327950288
28 #endif
29 
41 void sii_precomp(sii_coeffs *c, double sigma, int K)
42 {
43  /* Elboher and Werman's optimal radii and weights. */
44  const double sigma0 = 100.0 / M_PI;
45  static const short radii0[SII_MAX_K - SII_MIN_K + 1][SII_MAX_K] =
46  {{76, 46, 23, 0, 0},
47  {82, 56, 37, 19, 0},
48  {85, 61, 44, 30, 16}};
49  static const float weights0[SII_MAX_K - SII_MIN_K + 1][SII_MAX_K] =
50  {{0.1618f, 0.5502f, 0.9495f, 0, 0},
51  {0.0976f, 0.3376f, 0.6700f, 0.9649f, 0},
52  {0.0739f, 0.2534f, 0.5031f, 0.7596f, 0.9738f}};
53 
54  const int i = K - SII_MIN_K;
55  double sum;
56  int k;
57 
58  assert(c && sigma > 0 && SII_VALID_K(K));
59  c->K = K;
60 
61  for (k = 0, sum = 0; k < K; ++k)
62  {
63  c->radii[k] = (long)(radii0[i][k] * (sigma / sigma0) + 0.5);
64  sum += weights0[i][k] * (2 * c->radii[k] + 1);
65  }
66 
67  for (k = 0; k < K; ++k)
68  c->weights[k] = (num)(weights0[i][k] / sum);
69 
70  return;
71 }
72 
84 long sii_buffer_size(sii_coeffs c, long N)
85 {
86  long pad = c.radii[0] + 1;
87  return N + 2 * pad;
88 }
89 
108 void sii_gaussian_conv(sii_coeffs c, num *dest, num *buffer,
109  const num *src, long N, long stride)
110 {
111  num accum;
112  long pad, n;
113  int k;
114 
115  assert(dest && buffer && src && dest != buffer && src != buffer
116  && N > 0 && stride != 0);
117 
118  pad = c.radii[0] + 1;
119  buffer += pad;
120 
121  /* Compute cumulative sum of src over n = -pad,..., N + pad - 1. */
122  for (n = -pad, accum = 0; n < N + pad; ++n)
123  {
124  accum += src[stride * extension(N, n)];
125  buffer[n] = accum;
126  }
127 
128  /* Compute stacked box filters. */
129  for (n = 0; n < N; ++n, dest += stride)
130  {
131  accum = c.weights[0] * (buffer[n + c.radii[0]]
132  - buffer[n - c.radii[0] - 1]);
133 
134  for (k = 1; k < c.K; ++k)
135  accum += c.weights[k] * (buffer[n + c.radii[k]]
136  - buffer[n - c.radii[k] - 1]);
137 
138  *dest = accum;
139  }
140 
141  return;
142 }
143 
162 void sii_gaussian_conv_image(sii_coeffs c, num *dest, num *buffer,
163  const num *src, int width, int height, int num_channels)
164 {
165  long num_pixels = ((long)width) * ((long)height);
166  int x, y, channel;
167 
168  assert(dest && buffer && src && num_pixels > 0);
169 
170  /* Loop over the image channels. */
171  for (channel = 0; channel < num_channels; ++channel)
172  {
173  num *dest_y = dest;
174  const num *src_y = src;
175 
176  /* Filter each row of the channel. */
177  for (y = 0; y < height; ++y)
178  {
180  dest_y, buffer, src_y, width, 1);
181  dest_y += width;
182  src_y += width;
183  }
184 
185  /* Filter each column of the channel. */
186  for (x = 0; x < width; ++x)
188  dest + x, buffer, dest + x, height, width);
189 
190  dest += num_pixels;
191  src += num_pixels;
192  }
193 
194  return;
195 }