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## -*- texinfo -*-
17## @deftypefn {Function File} {@var{E} =} entropyfilt (@var{im})
18## @deftypefnx{Function File} {@var{E} =} entropyfilt (@var{im}, @var{domain})
19## @deftypefnx{Function File} {@var{E} =} entropyfilt (@var{im}, @var{domain}, @var{padding}, @dots{})
20## Computes the local entropy in a neighbourhood around each pixel in an image.
21##
22## The entropy of the elements of the neighbourhood is computed as
23##
24## @example
25## @var{E} = -sum (@var{P} .* log2 (@var{P})
26## @end example
27##
28## where @var{P} is the distribution of the elements of @var{im}. The distribution
29## is approximated using a histogram with @var{nbins} cells. If @var{im} is
30## @code{logical} then two cells are used. For other classes 256 cells
31## are used.
32##
33## When the entropy is computed, zero-valued cells of the histogram are ignored.
34##
35## The neighbourhood is defined by the @var{domain} binary mask. Elements of the
36## mask with a non-zero value are considered part of the neighbourhood. By default
37## a 9 by 9 matrix containing only non-zero values is used.
38##
39## At the border of the image, extrapolation is used. By default symmetric
40## extrapolation is used, but any method supported by the @code{padarray} function
41## can be used. Since extrapolation is used, one can expect a lower entropy near
42## the image border.
43##
44## @seealso{entropy, paddarray, stdfilt}
45## @end deftypefn
46
47function retval = entropyfilt (I, domain = true (9), padding = "symmetric", varargin)
48  ## Check input
49  if (nargin == 0)
50    error ("entropyfilt: not enough input arguments");
51  endif
52
53  if (! isnumeric (I))
54    error ("entropyfilt: I must be numeric");
55  endif
56
57  if (! isnumeric (domain) && ! islogical (domain))
58    error ("entropyfilt: DOMAIN must be a logical matrix");
59  endif
60  domain = logical (domain);
61
62  ## Get number of histogram bins
63  if (islogical (I))
64    nbins = 2;
65  else
66    nbins = 256;
67  endif
68
69  ## Convert to 8 or 16 bit integers if needed
70  ## (accepting single, int8, int16, int32, int64, uint64 is Octave-only)
71  switch (class (I))
72    case {"double", "single", "int16", "int32", "int64", "uint16", "uint32", "uint64"}
73      I = im2uint8 (I); # because this is what Matlab seems to do
74    case {"logical", "int8", "uint8"}
75      ## Do nothing
76    otherwise
77      error ("entropyfilt: cannot handle images of class '%s'", class (I));
78  endswitch
79
80  ## Pad image
81  pad = floor (size (domain)/2);
82  I = padarray (I, pad, padding, varargin {:});
83  even = (round (size (domain)/2) == size (domain)/2);
84  idx = cell (1, ndims (I));
85  for k = 1:ndims (I)
86    idx {k} = (even (k)+1):size (I, k);
87  endfor
88  I = I (idx {:});
89
90  retval = __spatial_filtering__ (I, domain, "entropy", zeros (size (domain)),
91                                  nbins);
92endfunction
93
94%!test
95%! a = log2 (9) * ones (5, 5);
96%! b = -(2*log2 (2/9) + log2 (1/9))/3;
97%! a(1,2:4) = b;
98%! a(5,2:4) = b;
99%! a(2:4,1) = b;
100%! a(2:4,5) = b;
101%! c = -(4*log2 (4/9) + 4*log2 (2/9) + log2 (1/9))/9;
102%! a(1,1) = c;
103%! a(5,1) = c;
104%! a(1,5) = c;
105%! a(5,5) = c;
106%! assert (entropyfilt (uint8 (magic (5)), ones (3, 3)), a, 2*eps);
107
108%!test
109%! assert (entropyfilt (uint8 (ones (10, 10))), zeros (10, 10));
110
111## some (Matlab compatible) tests on simple 2D-images (classes double, uint8, uint16):
112%!test
113%! A = zeros (3,3);
114%! B = ones (3,3);
115%! C = [1 1 1; 2 2 2; 3 3 3];
116%! D = C';
117%! E = ones (3,3);
118%! E(2,2) = 2;
119%! F = 3 .* ones (3,3);
120%! F(2,2) = 1;
121%! G = [-1 2 7; -5 2 8; -7 pi 9];
122%! H = [5 2 8; 1 -3 1; 5 1 0];
123%! Hf = mat2gray(H);
124%! X = uint8(abs(H));
125%! P = [0.2 0.201 0.204; 0.202 0.203 0.205; 0.205 0.206 0.202];
126%! Q = uint16([100 101 103; 100 105 102; 100 102 103]);
127%! R = uint8([1 2 3 4 5; 11 12 13 14 15; 21 22 4 5 6; 5 5 3 2 1; 15 14 14 14 14]);
128%! Aout = zeros (3);
129%! Bout = zeros (3);
130%! Cout = zeros (3);
131%! Dout = zeros (3);
132%! Eout = zeros (3);
133%! Fout = zeros (3);
134%! Gout_1 = -sum([2 7]./9.*log2([2 7]./9));
135%! Gout_2 = -sum([3 6]./9.*log2([3 6]./9));
136%! Gout_3 = -sum([4 5]./9.*log2([4 5]./9));
137%! Gout = [Gout_1 Gout_2 Gout_3; Gout_1 Gout_2 Gout_3; Gout_1 Gout_2 Gout_3];
138%! Hout_5 = -sum([2 7]./9.*log2([2 7]./9)) ;
139%! Hout = [0.8916 0.8256 0.7412; 0.8256 Hout_5 0.6913; 0.7412 0.6913 0.6355];
140%! Hfout_5 =  -sum([3 2 1 1 1 1]./9.*log2([3 2 1 1 1 1]./9));
141%! Hfout = [2.3613 2.3296 2.2252; 2.4571 Hfout_5 2.3090; 2.4805 2.4488 2.3445];
142%! Xout_5 = -sum([1 1 1 1 2 3]./9.*log2([1 1 1 1 2 3]./9));
143%! Xout  = [2.3613 2.3296 2.2252; 2.4571 Xout_5 2.3090; 2.4805 2.4488 2.3445];
144%! Pout_5 = -sum([1 2 6]./9.*log2([1 2 6]./9));
145%! Pout = [1.1137 1.1730 1.2251; 1.1595 Pout_5 1.2774; 1.1556 1.2183 1.2635];
146%! Qout = zeros(3);
147%! Rout = [3.5143 3.5700 3.4871 3.4957 3.4825;
148%!            3.4705 3.5330 3.4341 3.4246 3.3890;
149%!            3.3694 3.4063 3.3279 3.3386 3.3030;
150%!            3.3717 3.4209 3.3396 3.3482 3.3044;
151%!            3.4361 3.5047 3.3999 3.4236 3.3879];
152%! assert (entropyfilt (A), Aout);
153%! assert (entropyfilt (B), Bout);
154%! assert (entropyfilt (C), Cout);
155%! assert (entropyfilt (D), Dout);
156%! assert (entropyfilt (E), Eout);
157%! assert (entropyfilt (F), Fout);
158%! assert (entropyfilt (G), Gout, 1e-4);
159%! assert (entropyfilt (H), Hout, 1e-4);
160%! assert (entropyfilt (Hf), Hfout, 1e-4);
161%! assert (entropyfilt (X), Xout, 1e-4);
162%! assert (entropyfilt (P), Pout, 1e-4);
163%! assert (entropyfilt (Q), Qout);
164%! assert (entropyfilt (R), Rout, 1e-4);
165