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 <octave/oct.h>
17
18 template <class MT>
19 MT
custom_gaussian_smoothing(const MT & I,const Matrix & lambda1,const Matrix & lambda2,const Matrix & theta)20 custom_gaussian_smoothing (const MT &I, const Matrix &lambda1, const Matrix &lambda2,
21 const Matrix &theta)
22 {
23 const octave_idx_type rows = I.rows ();
24 const octave_idx_type cols = I.columns ();
25
26 // Allocate output
27 MT J (I.dims ());
28
29 // Iterate over every element of 'I'
30 for (octave_idx_type row = 0; row < rows; row++)
31 {
32 for (octave_idx_type col = 0; col < cols; col++)
33 {
34 // Extract parameters
35 const double v1 = lambda1 (row, col);
36 const double v2 = lambda2 (row, col);
37 const double t = theta (row, col);
38
39 // Should we perform any filtering?
40 if (std::min (v1, v2) > 0)
41 {
42 // Compute inverse covariance matrix, C^-1 = [a, b; b, c]
43 const double iv1 = 1.0/v1;
44 const double iv2 = 1.0/v2;
45 const double ct = cos (t);
46 const double st = sin (t);
47 const double ct2 = ct*ct;
48 const double st2 = st*st;
49 const double ctst = ct*st;
50 const double a = ct2*iv2 + st2*iv1;
51 const double b = (iv2-iv1)*ctst;
52 const double c = st2*iv2 + ct2*iv1;
53
54 // Compute bounding box of the filter
55 const double k = 3.0; // The maximally allowed Mahalanobis' distance
56 const double sqrtv1 = sqrt (v1);
57 const double sqrtv2 = sqrt (v2);
58 const octave_idx_type rur = abs (k*(ct*sqrtv2 - st*sqrtv1)); // 'rur' means 'row-upper-right'
59 const octave_idx_type cur = abs (k*(st*sqrtv2 + ct*sqrtv1));
60 const octave_idx_type rlr = abs (k*(ct*sqrtv2 + st*sqrtv1));
61 const octave_idx_type clr = abs (k*(st*sqrtv2 - ct*sqrtv1));
62 const octave_idx_type rul = abs (k*(-ct*sqrtv2 - st*sqrtv1));
63 const octave_idx_type cul = abs (k*(-st*sqrtv2 + ct*sqrtv1));
64 const octave_idx_type rll = abs (k*(-ct*sqrtv2 + st*sqrtv1));
65 const octave_idx_type cll = abs (k*(-st*sqrtv2 - ct*sqrtv1));
66 const octave_idx_type r_delta = std::max (std::max (rur, rlr), std::max (rul, rll));
67 const octave_idx_type c_delta = std::max (std::max (cur, clr), std::max (cul, cll));;
68 // The bounding box is now (row-r_delta):(row+r_delta)x(col-c_delta):(col+c_delta).
69 // We, however, represent the bounding box in a local coordinate system around (row, col).
70 const octave_idx_type r1 = std::max (row-r_delta, octave_idx_type (0)) - row;
71 const octave_idx_type r2 = std::min (row+r_delta, rows-1) - row;
72 const octave_idx_type c1 = std::max (col-c_delta, octave_idx_type (0)) - col;
73 const octave_idx_type c2 = std::min (col+c_delta, cols-1) - col;
74
75 // Perform the actual filtering
76 double sum = 0;
77 double wsum = 0; // for normalisation
78 for (octave_idx_type rl = r1; rl <= r2; rl++)
79 {
80 for (octave_idx_type cl = c1; cl <= c2; cl++)
81 {
82 // Compute Mahalanobis' distance
83 const double dsquare = rl*(a*rl + b*cl) + cl*(b*rl + c*cl);
84
85 // We only do the filtering in an elliptical window
86 if (dsquare > k*k) continue;
87
88 // Update filter values
89 const double w = exp (-0.5*dsquare);
90 wsum += w;
91 sum += w*(double)I.elem (row + rl, col + cl);
92 } // End: cl
93 } // End: rl
94
95 // Compute final result
96 J (row, col) = sum/wsum;
97 }
98 else // No filtering is performed
99 {
100 J.elem (row, col) = I.elem (row, col);
101 }
102 } // End: column iteration
103 } // End: row iteration
104
105 // Return
106 return J;
107 }
108
109 DEFUN_DLD (__custom_gaussian_smoothing__, args, ,"\
110 -*- texinfo -*-\n\
111 @deftypefn {Loadable Function} {@var{J} =} __custom_gaussian_smooting__ (@var{I}, @var{lambda1}, @var{lambda2}, @var{theta})\n\
112 Performs Gaussian smoothing on the image @var{I}. In pixel @math{(r,c)} the \n\
113 Eigenvalues of the Gaussian is @var{lambda1}@math{(r,c)} and @var{lambda2}@math{(r,c)}.\n\
114 The Gaussian is rotated with the angle given in @var{theta}@math{(r,c)}.\n\
115 \n\
116 @strong{Warning:} this function should @i{never} be called directly! The user\n\
117 interface to this function is available in @code{imsmooth}.\n\
118 @seealso{imsmooth}\n\
119 @end deftypefn\n\
120 ")
121 {
122 // Handle Input
123 octave_value_list retval;
124 const int nargin = args.length ();
125
126 if (nargin != 4)
127 print_usage ();
128
129 const Matrix lambda1 = args (1).matrix_value ();
130 const Matrix lambda2 = args (2).matrix_value ();
131 const Matrix theta = args (3).matrix_value ();
132
133 const octave_idx_type rows = args (0).rows();
134 const octave_idx_type cols = args (0).columns();
135 if (lambda1.rows () != rows || lambda1.columns () != cols
136 || lambda2.rows () != rows || lambda2.columns () != cols
137 || theta.rows () != rows || theta.columns () != cols)
138 error ("__custom_gaussian_smoothing__: size mismatch");
139
140 // Take action depending on input type
141 //octave_value J;
142 if (args(0).is_real_matrix())
143 {
144 const Matrix I = args(0).matrix_value();
145 retval.append (custom_gaussian_smoothing<Matrix>(I, lambda1, lambda2, theta));
146 }
147 else if (args(0).is_int8_type())
148 {
149 const int8NDArray I = args(0).int8_array_value();
150 retval.append (custom_gaussian_smoothing<int8NDArray>(I, lambda1, lambda2, theta));
151 }
152 else if (args(0).is_int16_type())
153 {
154 const int16NDArray I = args(0).int16_array_value();
155 retval.append (custom_gaussian_smoothing<int16NDArray>(I, lambda1, lambda2, theta));
156 }
157 else if (args(0).is_int32_type())
158 {
159 const int32NDArray I = args(0).int32_array_value();
160 retval.append (custom_gaussian_smoothing<int32NDArray>(I, lambda1, lambda2, theta));
161 }
162 else if (args(0).is_int64_type())
163 {
164 const int64NDArray I = args(0).int64_array_value();
165 retval.append (custom_gaussian_smoothing<int64NDArray>(I, lambda1, lambda2, theta));
166 }
167 else if (args(0).is_uint8_type())
168 {
169 const uint8NDArray I = args(0).uint8_array_value();
170 retval.append (custom_gaussian_smoothing<uint8NDArray>(I, lambda1, lambda2, theta));
171 }
172 else if (args(0).is_uint16_type())
173 {
174 const uint16NDArray I = args(0).uint16_array_value();
175 retval.append (custom_gaussian_smoothing<uint16NDArray>(I, lambda1, lambda2, theta));
176 }
177 else if (args(0).is_uint32_type())
178 {
179 const uint32NDArray I = args(0).uint32_array_value();
180 retval.append (custom_gaussian_smoothing<uint32NDArray>(I, lambda1, lambda2, theta));
181 }
182 else if (args(0).is_uint64_type())
183 {
184 const uint64NDArray I = args(0).uint64_array_value();
185 retval.append (custom_gaussian_smoothing<uint64NDArray>(I, lambda1, lambda2, theta));
186 }
187 else
188 error("__custom_gaussian_smoothing__: first input should be a real or integer array");
189
190 return retval;
191 }
192