1## Copyright (C) 2015 Carnë Draug <carandraug@octave.org>
2##
3## This program is free software; you can redistribute it and/or
4## modify it under the terms of the GNU General Public License as
5## published by the Free Software Foundation; either version 3 of the
6## License, or (at your option) any later version.
7##
8## This program is distributed in the hope that it will be useful, but
9## WITHOUT ANY WARRANTY; without even the implied warranty of
10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11## General Public License for more details.
12##
13## You should have received a copy of the GNU General Public License
14## along with this program; if not, see
15## <http://www.gnu.org/licenses/>.
16
17## -*- texinfo -*-
18## @deftypefn {Function File} {} edgetaper (@var{img}, @var{psf})
19## Blur border (edges) of image to prevent ringing artifacts.
20##
21## @emph{Warning}: this function is not @sc{Matlab} compatible and
22## is likely to change in the future.
23##
24## @end deftypefn
25
26function imgt = edgetaper (img, psf)
27  if (nargin != 2)
28    print_usage ();
29  elseif (! isnumeric (img))
30    error ("edgetaper: IMG must be numeric")
31  elseif (! isnumeric (psf))
32    error ("edgetaper: PSF must be numeric")
33  endif
34
35  img_size = size (img);
36  psf_size = size (psf);
37  n = max (numel (img_size), numel (psf_size));
38  img_size = postpad (img_size(:), n, 1);
39  psf_size = postpad (psf_size(:), n, 1);
40
41  ## beware of psf_size = [16 16 1 8] paired with img_size [512 512 1 50]
42  ## which are valid, singleton dimensions do not count for this check
43  if (any ((psf_size > (img_size / 2)) & (psf_size > 1)))
44    error ("edgetaper: PSF must be smaller than half of IMG dimensions")
45  endif
46  psf = psf ./ sum (psf(:));  # we use it for blurring so the sum must be 1
47
48  ## FIXME  this function is not Matlab compatible.  I have no clue what
49  ##        Matlab is doing but is definitely not what they claim on the
50  ##        documentation.  There are no references for the function and
51  ##        the documentation is sparse and incorrect.
52  ##
53  ##        I have found papers that compare their method for reducing
54  ##        boundary artifacts against Matlab's edgetaper but even they
55  ##        do not comment it.
56  ##
57  ##        It an implementation of edgetaper that is close to Matlab's
58  ##        documentation for the function.  If anyone has patience, please
59  ##        fix this.
60  ##
61  ##        Some questions about it:
62  ##
63  ##          1. the autocorrelation of the PSF is twice the size of the PSF
64  ##          but it will still be way smaller than the image.  How can it be
65  ##          used to make a weighted mean between the blurred and the original
66  ##          image? We pretty much split it and only use it on the borders
67  ##          but looks like we end up with an image too blurred.
68  ##
69  ##          2. how do they blur the image? I will guess they pad the
70  ##          image but how?
71  ##
72  ##  Note: always test this function with input as double precision
73  ##        since it returns the same class as input and may hide our
74  ##        differences.
75
76  blurred = fftconvn (img, psf, "same");
77  xpsf = normalized_autocorrelation (psf);
78
79  ## The following will expand a ND matrix into a larger size
80  ## repeating the center elements, e.g.,
81  ##
82  ##              1 2 2 2 3
83  ##  1 2 3       4 5 5 5 6
84  ##  4 5 6   =>  4 5 5 5 6
85  ##  7 8 9       4 5 5 5 6
86  ##              7 8 8 8 9
87
88  ## note the xpsf will always have odd sizes (2*psf_size +1)
89  xdims = ndims (xpsf);
90  xpsf_size = size (xpsf)(:);
91  idim = arrayfun (@(c, n, f) [1:c repmat(c, [1 n]) c:f], ceil (xpsf_size /2),
92                   img_size(1:xdims) - xpsf_size -1, xpsf_size,
93                  "UniformOutput", false);
94  subs = cell (xdims, 1);
95  [subs{:}] = ndgrid (idim{:});
96  inds = sub2ind (xpsf_size, subs{:});
97  weights = xpsf(inds);
98
99  imgt = (img .* weights) + (blurred .* (1 - weights));
100
101  imgt = cast (imgt, class (img));
102endfunction
103
104function acn = normalized_autocorrelation (psf)
105  idx = arrayfun (@colon, size (psf), repmat (-1, [1 ndims(psf)]),
106                  repmat (1, [1 ndims(psf)]), "UniformOutput", false);
107  ac = convn (psf, conj (psf(idx{:})));
108  acn = ac / max (ac(:));
109endfunction
110
111%!assert (class (edgetaper (rand (100), rand (16))), "double")
112%!assert (class (edgetaper (randi (255, 100, "uint8"), rand (16))), "uint8")
113
114