1 // Copyright (C) 2008 Søren Hauberg <soren@hauberg.org>
2 //
3 // This program is free software; you can redistribute it and/or modify it under
4 // the terms of the GNU General Public License as published by the Free Software
5 // Foundation; either version 3 of the License, or (at your option) any later
6 // version.
7 //
8 // This program is distributed in the hope that it will be useful, but WITHOUT
9 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
11 // details.
12 //
13 // You should have received a copy of the GNU General Public License along with
14 // this program; if not, see <http://www.gnu.org/licenses/>.
15
16 #include <vector>
17
18 #include <octave/oct.h>
19
20 inline
gauss(const std::vector<double> x,const std::vector<double> mu,const double sigma)21 double gauss (const std::vector<double> x, const std::vector<double> mu,
22 const double sigma)
23 {
24 double s = 0;
25 for (size_t i = 0; i < x.size (); i++)
26 {
27 const double d = x[i] - mu[i];
28 s += d*d;
29 }
30 return exp (-0.5*s/(sigma*sigma));
31 }
32
33 template <class MatrixType>
34 octave_value
bilateral(const MatrixType & im,const double sigma_d,const double sigma_r)35 bilateral (const MatrixType &im, const double sigma_d, const double sigma_r)
36 {
37 // Get sizes
38 const octave_idx_type ndims = im.ndims ();
39 const dim_vector size = im.dims ();
40 const octave_idx_type num_planes = (ndims == 2) ? 1 : size (2);
41
42 // Build spatial kernel
43 const int s = std::max ((int)std::round (3*sigma_d), 1);
44 Matrix kernel (2*s+1, 2*s+1);
45 for (octave_idx_type r = 0; r < 2*s+1; r++)
46 {
47 for (octave_idx_type c = 0; c < 2*s+1; c++)
48 {
49 const int dr = r-s;
50 const int dc = c-s;
51 kernel (r,c) = exp (-0.5 * (dr*dr + dc*dc)/(sigma_d*sigma_d));
52 }
53 }
54
55 // Allocate output
56 dim_vector out_size (size);
57 out_size (0) = std::max (size (0) - 2*s, (octave_idx_type)0);
58 out_size (1) = std::max (size (1) - 2*s, (octave_idx_type)0);
59 MatrixType out = MatrixType (out_size);
60
61 // Iterate over every element of 'out'.
62 for (octave_idx_type r = 0; r < out_size (0); r++)
63 {
64 for (octave_idx_type c = 0; c < out_size (1); c++)
65 {
66 OCTAVE_QUIT;
67
68 // For each neighbour
69 std::vector<double> val (num_planes);
70 std::vector<double> sum (num_planes);
71 double k = 0;
72 for (octave_idx_type i = 0; i < num_planes; i++)
73 {
74 val[i] = im (r+s,c+s,i);
75 sum[i] = 0;
76 }
77 for (octave_idx_type kr = 0; kr < 2*s+1; kr++)
78 {
79 for (octave_idx_type kc = 0; kc < 2*s+1; kc++)
80 {
81 std::vector<double> lval (num_planes);
82 for (octave_idx_type i = 0; i < num_planes; i++)
83 lval[i] = im (r+kr, c+kc, i);
84 const double w = kernel (kr, kc) * gauss (val, lval, sigma_r);
85 for (octave_idx_type i = 0; i < num_planes; i++)
86 sum[i] += w * lval[i];
87 k += w;
88 }
89 }
90 for (octave_idx_type i = 0; i < num_planes; i++)
91 out (r, c, i) = sum[i]/k;
92 }
93 }
94
95 return octave_value (out);
96 }
97
98 DEFUN_DLD (__bilateral__, args, , "\
99 -*- texinfo -*-\n\
100 @deftypefn {Loadable Function} __bilateral__(@var{im}, @var{sigma_d}, @var{sigma_r})\n\
101 Performs Gaussian bilateral filtering in the image @var{im}. @var{sigma_d} is the\n\
102 spread of the Gaussian used as closenes function, and @var{sigma_r} is the spread\n\
103 of Gaussian used as similarity function. This function is internal and should NOT\n\
104 be called directly. Instead use @code{imsmooth}.\n\
105 @end deftypefn\n\
106 ")
107 {
108 octave_value_list retval;
109 if (args.length () != 3)
110 print_usage ();
111
112 const octave_idx_type ndims = args (0).ndims ();
113 if (ndims != 2 && ndims != 3)
114 error ("__bilateral__: only 2 and 3 dimensional is supported");
115
116 const double sigma_d = args (1).scalar_value ();
117 const double sigma_r = args (2).scalar_value ();
118
119 // Take action depending on input type
120 if (args (0).is_real_matrix ())
121 {
122 const NDArray im = args(0).array_value ();
123 retval = bilateral<NDArray> (im, sigma_d, sigma_r);
124 }
125 else if (args (0).is_int8_type ())
126 {
127 const int8NDArray im = args (0).int8_array_value ();
128 retval = bilateral<int8NDArray> (im, sigma_d, sigma_r);
129 }
130 else if (args (0).is_int16_type ())
131 {
132 const int16NDArray im = args (0).int16_array_value ();
133 retval = bilateral<int16NDArray> (im, sigma_d, sigma_r);
134 }
135 else if (args (0).is_int32_type ())
136 {
137 const int32NDArray im = args (0).int32_array_value ();
138 retval = bilateral<int32NDArray> (im, sigma_d, sigma_r);
139 }
140 else if (args (0).is_int64_type ())
141 {
142 const int64NDArray im = args (0).int64_array_value ();
143 retval = bilateral<int64NDArray> (im, sigma_d, sigma_r);
144 }
145 else if (args (0).is_uint8_type ())
146 {
147 const uint8NDArray im = args (0).uint8_array_value ();
148 retval = bilateral<uint8NDArray> (im, sigma_d, sigma_r);
149 }
150 else if (args(0).is_uint16_type())
151 {
152 const uint16NDArray im = args (0).uint16_array_value ();
153 retval = bilateral<uint16NDArray> (im, sigma_d, sigma_r);
154 }
155 else if (args (0).is_uint32_type ())
156 {
157 const uint32NDArray im = args (0).uint32_array_value ();
158 retval = bilateral<uint32NDArray> (im, sigma_d, sigma_r);
159 }
160 else if (args (0).is_uint64_type ())
161 {
162 const uint64NDArray im = args (0).uint64_array_value ();
163 retval = bilateral<uint64NDArray> (im, sigma_d, sigma_r);
164 }
165 else
166 error ("__bilateral__: first input should be a real or integer array");
167
168 return retval;
169 }
170