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