1 /**
2 * \file gaussian_conv_sii.c
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 #include "gaussian_conv_sii.hh"
28 #include <assert.h>
29 #include <stdio.h>
30 //#include "filter_util.h"
31
32 #ifndef M_PI
33 /** \brief The constant pi */
34 #define M_PI 3.14159265358979323846264338327950288
35 #endif
36
37 /**
38 * \brief Precompute filter coefficients for SII Gaussian convolution
39 * \param c sii_coeffs pointer to hold precomputed coefficients
40 * \param sigma Gaussian standard deviation
41 * \param K number of boxes = 3, 4, or 5
42 * \return 1 on success, 0 on failure
43 *
44 * This routine reads Elboher and Werman's optimal SII radii and weights for
45 * reference standard deviation \f$ \sigma_0 = 100/\pi \f$ from a table and
46 * scales them to the specified value of sigma.
47 */
sii_precomp(sii_coeffs * c,double sigma,int K)48 void sii_precomp(sii_coeffs *c, double sigma, int K)
49 {
50 /* Elboher and Werman's optimal radii and weights. */
51 const double sigma0 = 100.0 / M_PI;
52 static const short radii0[SII_MAX_K - SII_MIN_K + 1][SII_MAX_K] =
53 {{76, 46, 23, 0, 0},
54 {82, 56, 37, 19, 0},
55 {85, 61, 44, 30, 16}};
56 static const float weights0[SII_MAX_K - SII_MIN_K + 1][SII_MAX_K] =
57 {{0.1618f, 0.5502f, 0.9495f, 0, 0},
58 {0.0976f, 0.3376f, 0.6700f, 0.9649f, 0},
59 {0.0739f, 0.2534f, 0.5031f, 0.7596f, 0.9738f}};
60
61 const int i = K - SII_MIN_K;
62 double sum;
63 int k;
64
65 assert(c && sigma > 0 && SII_VALID_K(K));
66 c->K = K;
67
68 for (k = 0, sum = 0; k < K; ++k)
69 {
70 c->radii[k] = (long)(radii0[i][k] * (sigma / sigma0) + 0.5);
71 sum += weights0[i][k] * (2 * c->radii[k] + 1);
72 }
73
74 for (k = 0; k < K; ++k)
75 c->weights[k] = (double)(weights0[i][k] / sum);
76
77 return;
78 }
79
80 /**
81 * \brief Determines the buffer size needed for SII Gaussian convolution
82 * \param c sii_coeffs created by sii_precomp()
83 * \param N number of samples
84 * \return required buffer size in units of num samples
85 *
86 * This routine determines the minimum size of the buffer needed for use in
87 * sii_gaussian_conv() or sii_gaussian_conv_image(). This size is the length
88 * of the signal (or in 2D, max(width, height)) plus the twice largest box
89 * radius, for padding.
90 */
sii_buffer_size(sii_coeffs c,long N)91 long sii_buffer_size(sii_coeffs c, long N)
92 {
93 long pad = c.radii[0] + 1;
94 return N + 2 * pad;
95 }
96
97