## Copyright (C)
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
##

## -*- texinfo -*-
## @deftypefn  {Function File} {} impyramid (@var{im})
## @deftypefnx {Function File} {} impyramid (@var{im}, @var{direction})
## @var{direction} can be 'reduce' or 'expand'
##
## Create image which is one level up or down in the Gaussian pyramid
##
## The 'reduce' stage is done by low-pass filtering and subsampling of 1:2
## in each axis. If the size of the original image is [M N], the
## size of the reduced image is [ceil((M+1)/2) ceil((N+1)/2)].
##
## The 'expand' stage is done by upsampling the image (2:1 in each axis),
## and then low-pass filtering. If the size of the original image is [M N], the
## size of the exapnded image is [2M 2N].
## 
## Note that image processing pyramids are upside down, so 'reduce' is
## is going one level DOWN in the pyramid, while 'expand' is going 
## one level UP in the pyramid...
##
## The low-pass filter is according to Burt & Adelson paper [1]:
## W(i,j) = w(i)w(j).   
## w = [0.25-alpha/2 0.25 alpha 0.25 0.25-alpha/2]  with alpha=0.375
## 
## [1] Burt and Adelson, "The Laplacian Pyramid as a Compact Image Code",
## IEEE Transactions on Communications, vol. COM-31, no. 4,
## April 1983, pp. 532-540.
## {http://persci.mit.edu/pub_pdfs/pyramid83.pdf}
##
## @example
## @group
## impyramid (im);             # return reduced image (create one level down image)
## impyramid (im, 'reduce');   # return reduced image (same as impyramid (im))
## impyramid (im, 'expand');   # return expanded image (one level up)
## @end group
## @end example
##
## @seealso{imresize, imfilter}
## @end deftypefn

function imp = impyramid (im, direction = "reduce")

  if (nargin < 1 || nargin > 2)
    print_usage ();
  elseif (! isimage (im) || (! isrgb (im) && ! isgray (im)))
    error ("impyramid: IM must be a grayscale or RGB image.")
  elseif (! ischar (direction))
    error ("impyramid: direction must be a string with direction definition")
  endif

  direction = tolower (direction);
  
  ## low pass filter to be used
  alpha = 0.375;
  filt = [0.25-alpha/2 0.25 alpha 0.25 0.25-alpha/2];
  
  r = rows (im);
  c = columns (im);
  
  switch (direction)
    case "reduce"
      ##  perform horizontal low pass filtering  
      im1 = imfilter (im, filt);  
      ##  perform vertical low pass filtering  
      im2 = imfilter (im1, filt');
	
      ##  subsampling	
      if (isgray(im))
        imp = im2(1:2:r,1:2:c);
      else
        imp = im2(1:2:r,1:2:c,:);
      end
    case "expand"
	    if (isgray(im))
      ## creating image with dimensions that are twice 
        im1 = zeros (r*2,c*2, class(im));
    	## put original image in odd pixels
        im1 (1:2:2*r,1:2:2*c) = im;
    	else
      ## creating image with dimensions that are twice 
        im1 = zeros (r*2,c*2, 3, class(im));
    	## put original image in odd pixels
        im1 (1:2:2*r,1:2:2*c,:) = im;
	    endif
	
      ##  perform horizontal low pass filtering  	
      im2 = 2*imfilter (im1, filt);
      ##  perform vertical low pass filtering  
      imp = 2*imfilter (im2, filt');   
    otherwise
      error ("impyramid: direction must be 'reduce' or 'expand'")
  endswitch
endfunction

%!test
%! in = [116  227  153   69  146  194   59  130  139  106
%!         2   47  137  249   90   75   16   24  158   44
%!       155   68   46   84  166  156   69  204   32  152
%!        71  221  137  230  210  153  192  115   30  118
%!       107  143  108   52   51   73  101   21  175   90
%!        54  158  143   77   26  168  113  229  165  225
%!         9   47  133  135  130  207  236   43   19   73]
%!
%! reduced = [  51   93   89   65   69
%!              69  122  140  110   93
%!              77  123  112  124  113
%!              31   74   85  103   71];
%!
%! expand = [ 88  132  160  154  132  108   94  102  120  138  138  100   66   74   96  112  116  104   78   44
%!            62   98  128  142  146  154  154  140  126  126  122   86   54   58   82  114  132  112   74   38
%!            36   54   74  100  130  168  184  156  118  104   92   64   40   44   66  100  122  104   66   32
%!            66   68   64   76   98  130  154  148  132  122  108   80   60   78  104  106   98   98   86   50
%!           104  106   88   78   78   96  122  144  154  154  140  112   98  124  144  110   74   92  106   68
%!           102  130  134  120  108  126  154  174  180  172  156  142  138  146  140   96   60   84  106   68
%!            88  140  170  158  140  156  180  188  180  164  152  154  156  140  112   82   66   84   96   60
%!            90  136  164  154  134  132  138  136  130  122  120  130  134  108   82   86  100  104   92   54
%!            92  126  142  136  116   96   80   74   72   82   94  106  106   88   78  108  138  132  102   56
%!            80  116  140  138  122   96   68   52   52   80  110  114  112  118  128  148  164  164  140   80
%!            58   98  132  140  130  110   82   62   62  102  142  144  138  154  168  164  156  170  162   94
%!            36   68  100  120  130  122  106   92   96  134  174  182  172  156  136  116  104  122  124   74
%!            16   34   58   86  108  114  110  106  112  138  170  184  172  126   74   48   44   60   68   42
%!             6   14   28   46   62   68   68   66   70   84  102  112  104   70   32   16   14   24   30   18];
%!
%! assert (impyramid (uint8(in)), uint8(reduced))
%! assert (impyramid (uint8(in), 'reduce'), uint8(reduced))
%! assert (impyramid (uint8(in), 'expand'), uint8(expand))
