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