1 /**
2 * \file gaussian_conv_sii.h
3 * \brief Gaussian convolution using stacked integral images
4 * \author Pascal Getreuer <getreuer@cmla.ens-cachan.fr>
5 * http://www.ipol.im/pub/art/2013/87/
6 *
7 * This code implements stacked integral images Gaussian convolution approximation
8 * introduced by Bhatia, Snyder, and Bilbro (http://dx.doi.org/10.1109/ROBOT.2010.5509400)
9 * and refined by Elboher and Werman (http://arxiv.org/abs/1107.4958).
10 *
11 * Modified for integration into PhotoFlow
12 *
13 * Copyright (c) 2012-2013, Pascal Getreuer
14 * All rights reserved.
15 *
16 * This program is free software: you can redistribute it and/or modify it
17 * under, at your option, the terms of the GNU General Public License as
18 * published by the Free Software Foundation, either version 3 of the
19 * License, or (at your option) any later version, or the terms of the
20 * simplified BSD license.
21 *
22 * You should have received a copy of these licenses along with this program.
23 * If not, see <http://www.gnu.org/licenses/> and
24 * <http://www.opensource.org/licenses/bsd-license.html>.
25 */
26
27 /**
28 * \defgroup sii_gaussian Stacked integral images Gaussian convolution
29 * \brief A sum of box filter responses, rather than a cascade.
30 *
31 * This code implements stacked integral images Gaussian convolution
32 * approximation introduced by Bhatia, Snyder, and Bilbro and refined by
33 * Elboher and Werman. The Gaussian is approximated as
34 * \f[ u_n = \sum_{k=1}^K w_k (s_{n+r_k} - s_{n-r_k-1}), \quad
35 s_n = \sum_{m\le n} f_n, \f]
36 * where the kth term of the sum is effectively a box filter of radius
37 * \f$ r_k \f$.
38 *
39 * The process to use these functions is the following:
40 * -# sii_precomp() to precompute coefficients for the convolution
41 * -# sii_gaussian_conv() or sii_gaussian_conv_image() to perform
42 * the convolution itself (may be called multiple times if desired)
43 *
44 * The function sii_buffer_size() should be used to determine the minimum
45 * required buffer size.
46 *
47 * \par Example
48 \code
49 sii_coeffs c;
50 num *buffer;
51
52 sii_precomp(&c, sigma, K);
53 buffer = (num *)malloc(sizeof(num) * sii_buffer_size(c, N));
54 sii_gaussian_conv(c, dest, buffer, src, N, stride);
55 free(buffer);
56 \endcode
57 *
58 * \par References
59 * - A. Bhatia, W.E. Snyder, G. Bilbro, "Stacked Integral Image," IEEE
60 * International Conference on Robotics and Automation (ICRA),
61 * pp. 1530-1535, 2010. http://dx.doi.org/10.1109/ROBOT.2010.5509400
62 * - E. Elboher, M. Werman, "Efficient and Accurate Gaussian Image Filtering
63 * Using Running Sums," Computing Research Repository,
64 * vol. abs/1107.4958, 2011. http://arxiv.org/abs/1107.4958
65 *
66 * \{
67 */
68
69 #ifndef GAUSSIAN_CONV_SII_H
70 #define GAUSSIAN_CONV_SII_H
71
72 //#include "num.h"
73 //#include "gaussblur.hh"
74
75 /** \brief Minimum SII filter order */
76 #define SII_MIN_K 3
77 /** \brief Maximum SII filter order */
78 #define SII_MAX_K 5
79 /** \brief Test whether a given K value is a valid SII filter order */
80 #define SII_VALID_K(K) (SII_MIN_K <= (K) && (K) <= SII_MAX_K)
81
82 /** \brief Parameters for stacked integral images Gaussian approximation */
83 typedef struct sii_coeffs_
84 {
85 double weights[SII_MAX_K]; /**< Box weights */
86 long radii[SII_MAX_K]; /**< Box radii */
87 int K; /**< Number of boxes */
88 } sii_coeffs;
89
90 void sii_precomp(sii_coeffs *c, double sigma, int K);
91 long sii_buffer_size(sii_coeffs c, long N);
92
93
94
95 /**
96 * \brief 2D Gaussian convolution SII approximation
97 * \param c sii_coeffs created by sii_precomp()
98 * \param dest output convolved data
99 * \param buffer array with space for sii_buffer_size() samples
100 * \param src image to be convolved, overwritten if src = dest
101 * \param width image width
102 * \param height image height
103 * \param num_channels number of image channels
104 *
105 * Similar to sii_gaussian_conv(), this routine approximates 2D Gaussian
106 * convolution with stacked integral images. The buffer array must have space
107 * for at least sii_buffer_size(c,max(width,height)) samples.
108 *
109 * The convolution can be performed in-place by setting `src` = `dest` (the
110 * source array is overwritten with the result). However, the buffer array
111 * must be distinct from `src` and `dest`.
112 */
113 template<class num>
sii_gaussian_conv_image(sii_coeffs c,num * dest,num * buffer,const num * src,int width,int height,int num_channels)114 void sii_gaussian_conv_image(sii_coeffs c, num *dest, num *buffer,
115 const num *src, int width, int height, int num_channels)
116 {
117 long num_pixels = ((long)width) * ((long)height);
118 int x, y, channel;
119
120 assert(dest && buffer && src && num_pixels > 0);
121
122 /* Loop over the image channels. */
123 for (channel = 0; channel < num_channels; ++channel)
124 {
125 num *dest_y = dest;
126 const num *src_y = src;
127
128 /* Filter each row of the channel. */
129 for (y = 0; y < height; ++y)
130 {
131 //sii_gaussian_conv(c,
132 // dest_y, buffer, src_y, width, 1);
133 dest_y += width;
134 src_y += width;
135 }
136
137 /* Filter each column of the channel. */
138 for (x = 0; x < width; ++x)
139 //sii_gaussian_conv(c,
140 // dest + x, buffer, dest + x, height, width);
141
142 dest += num_pixels;
143 src += num_pixels;
144 }
145
146 return;
147 }
148
149 /** \} */
150 #endif /* GAUSSIAN_CONV_SII_H */
151