octave_packages/image-1.0.15/fftconv2.m

   1 ## Copyright (C) 2004 Stefan van der Walt <stefan@sun.ac.za>
   2 ##
   3 ## This program is free software; redistribution and use in source and
   4 ## binary forms, with or without modification, are permitted provided that
   5 ## the following conditions are met:
   6 ##
   7 ## 1. Redistributions of source code must retain the above copyright
   8 ##    notice, this list of conditions and the following disclaimer.
   9 ## 2. Redistributions in binary form must reproduce the above copyright
  10 ##    notice, this list of conditions and the following disclaimer in the
  11 ##    documentation and/or other materials provided with the distribution.
  12 ##
  13 ## THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  14 ## ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  15 ## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  16 ## ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  17 ## FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  18 ## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  19 ## OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  20 ## HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  21 ## LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  22 ## OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  23 ## SUCH DAMAGE.
  24
  25 ## -*- texinfo -*-
  26 ## @deftypefn {Function File} fftconv2 (@var{a}, @var{b}, @var{shape})
  27 ## @deftypefnx{Function File} fftconv2 (@var{v1}, @var{v2}, @var{a}, @var{shape})
  28 ## Convolve 2 dimensional signals using the FFT.
  29 ##
  30 ## This method is faster but less accurate than @var{conv2} for large @var{a} and @var{b}.
  31 ## It also uses more memory. A small complex component will be
  32 ## introduced even if both @var{a} and @var{b} are real.
  33 ## @seealso{conv2}
  34 ## @end deftypefn
  35
  36 ## Author: Stefan van der Walt <stefan@sun.ac.za>
  37 ## Date: 2004
  38
  39 function X = fftconv2(varargin)
  40     if (nargin < 2)
  41         usage("fftconv2(a,b[,shape]) or fftconv2(v1, v2, a, shape)")
  42     endif
  43
  44     shape = "full";
  45     rowcolumn = 0;
  46
  47     if ((nargin > 2) && ismatrix(varargin{3}) && !ischar(varargin{3}))
  48         ## usage: fftconv2(v1, v2, a[, shape])
  49
  50         rowcolumn = 1;
  51         v1 = varargin{1}(:)';
  52         v2 = varargin{2}(:);
  53         orig_a = varargin{3};
  54
  55         if (nargin == 4) shape = varargin{4}; endif
  56     else
  57         ## usage: fftconv2(a, b[, shape])
  58
  59         a = varargin{1};
  60         b = varargin{2};
  61         if (nargin == 3) shape = varargin{3}; endif
  62
  63     endif
  64
  65     if (rowcolumn)
  66         a = fftconv2(orig_a, v2);
  67         b = v1;
  68     endif
  69
  70     ra = rows(a);
  71     ca = columns(a);
  72     rb = rows(b);
  73     cb = columns(b);
  74
  75     A = fft2(impad(a, [0 cb-1], [0 rb-1]));
  76     B = fft2(impad(b, [0 ca-1], [0 ra-1]));
  77
  78     X = ifft2(A.*B);
  79
  80     if (rowcolumn)
  81         rb = rows(v2);
  82         ra = rows(orig_a);
  83         cb = columns(v1);
  84         ca = columns(orig_a);
  85     endif
  86
  87     if strcmp(shape,"same")
  88         r_top = ceil((rb + 1) / 2);
  89         c_top = ceil((cb + 1) / 2);
  90         X = X(r_top:r_top + ra - 1, c_top:c_top + ca - 1);
  91     elseif strcmp(shape, "valid")
  92         X = X(rb:ra, cb:ca);
  93     endif
  94 endfunction
  95
  96 %!# usage: fftconv2(a,b,[, shape])
  97 %!shared a,b
  98 %! a = repmat(1:10, 5);
  99 %! b = repmat(10:-1:3, 7);
 100 %!assert(norm(fftconv2(a,b)-conv2(a,b)), 0, 1e6*eps)
 101 %!assert(norm(fftconv2(b,a)-conv2(b,a)), 0, 1e6*eps)
 102 %!assert(norm(fftconv2(a,b,'full')-conv2(a,b,'full')), 0, 1e6*eps)
 103 %!assert(norm(fftconv2(b,a,'full')-conv2(b,a,'full')), 0, 1e6*eps)
 104 %!assert(norm(fftconv2(a,b,'same')-conv2(a,b,'same')), 0, 1e6*eps)
 105 %!assert(norm(fftconv2(b,a,'same')-conv2(b,a,'same')), 0, 1e6*eps)
 106 %!assert(isempty(fftconv2(a,b,'valid')));
 107 %!assert(norm(fftconv2(b,a,'valid')-conv2(b,a,'valid')), 0, 1e6*eps)
 108
 109 %!# usage: fftconv2(v1, v2, a[, shape])
 110 %!shared x,y,a
 111 %! x = 1:4; y = 4:-1:1; a = repmat(1:10, 5);
 112 %!assert(norm(fftconv2(x,y,a)-conv2(x,y,a)), 0, 1e6*eps)
 113 %!assert(norm(fftconv2(x,y,a,'full')-conv2(x,y,a,'full')), 0, 1e6*eps)
 114 %!assert(norm(fftconv2(x,y,a,'same')-conv2(x,y,a,'same')), 0, 1e6*eps)
 115 %!assert(norm(fftconv2(x,y,a,'valid')-conv2(x,y,a,'valid')), 0, 1e6*eps)
 116
 117 %!demo
 118 %! ## Draw a cross
 119 %! N = 100;
 120 %! [x,y] = meshgrid(-N:N, -N:N);
 121 %! z = 0*x;
 122 %! z(N,1:2*N+1) = 1; z(1:2*N+1, N) = 1;
 123 %! imshow(z);
 124 %!
 125 %! ## Draw a sinc blob
 126 %! n = floor(N/10);
 127 %! [x,y] = meshgrid(-n:n, -n:n);
 128 %! b = x.^2 + y.^2; b = max(b(:)) - b; b = b / max(b(:));
 129 %! imshow(b);
 130 %!
 131 %! ## Convolve the cross with the blob
 132 %! imshow(real(fftconv2(z, b, 'same')*N))
 133
 134