function [b1,b2] = blkproc(varargin)
%BLKPROC Implement distinct block processing for image.
%   B = BLKPROC(A,[M N],FUN) processes the image A by applying
%   the function FUN to each distinct M-by-N block of A, padding
%   A with zeros if necessary. FUN can be an inline function, a
%   text string containing the name of a function, or a string
%   containing an expression. FUN should operate on an M-by-N
%   block, X, and return a matrix, vector, or scalar Y:
%
%        Y = FUN(X)
%
%   BLKPROC does not require that Y be the same size as X.
%   However, B is the same size as A only if Y is the same size
%   as X. 
%
%   B = BLKPROC(A,[M N],FUN,P1,P2,...) passes the additional
%   parameters P1,P2,..., to FUN.
%
%   B = BLKPROC(A,[M N],[MBORDER NBORDER],FUN,...) defines an
%   overlapping border around the blocks.  BLKPROC extends the
%   original M-by-N blocks by MBORDER on the top and bottom, and
%   NBORDER on the left and right, resulting in blocks of size
%   (M+2*MBORDER)-by-(N+2*NBORDER). BLKPROC pads the border with
%   zeros, if necessary, on the edges of A. FUN should operate
%   on the extended block.
%
%   B = BLKPROC(A,'indexed',...) processes A as an indexed image,
%   padding with zeros if the class of A is uint8, or ones if the
%   class of A is double.
%
%   Class Support
%   -------------
%   The input image A can be of any class supported by FUN. The
%   class of B depends on the class of the output from FUN.
%
%   Example
%   -------
%   This example uses BLKPROC to set the pixels in each 8-by-8
%   block to the standard deviation of the elements in that
%   block.
%   
%       I = imread('alumgrns.tif');
%       I2 = blkproc(I,[8 8],'std2(x)*ones(size(x))');
%       imshow(I)
%       figure, imshow(I2,[])
%
%   See also COLFILT, NLFILTER.

%   Clay M. Thompson 10-6-92
%   Revised Steven L. Eddins, May 1996
%   Copyright 1993-1998 The MathWorks, Inc.  All Rights Reserved.
%   $Revision: 5.13 $  $Date: 1997/11/24 15:33:58 $

[a, block, border, fun, params, padval] = parse_inputs(varargin{:});

% Expand A: Add border, pad if size(a) is not divisible by block.
[ma,na] = size(a);
mpad = rem(ma,block(1)); if mpad>0, mpad = block(1)-mpad; end
npad = rem(na,block(2)); if npad>0, npad = block(2)-npad; end
if (isa(a, 'uint8'))
    if (padval == 1)
        aa = repmat(uint8(1), ma+mpad+2*border(1),na+npad+2*border(2));
    else
        aa = repmat(uint8(0), ma+mpad+2*border(1),na+npad+2*border(2));
    end
else
    if (padval == 1)
        aa = ones(ma+mpad+2*border(1),na+npad+2*border(2));
    else
        aa = zeros(ma+mpad+2*border(1),na+npad+2*border(2));
    end
end
aa(border(1)+(1:ma),border(2)+(1:na)) = a;

%
% Process first block.
%
m = block(1) + 2*border(1);
n = block(2) + 2*border(2);
mblocks = (ma+mpad)/block(1);
nblocks = (na+npad)/block(2);
arows = 1:m; acols = 1:n;
x = aa(arows, acols);
%firstBlock = feval(fun,x,params{:});
[firstBlock1, firstBlock2] = feval(fun,x,params{:});
%if (isempty(firstBlock))
if ( isempty(firstBlock1) | isempty(firstBlock2) )
  style = 'e'; % empty
  b1 = [];
  b2 = [];
%elseif (all(size(firstBlock) == size(x)))
elseif (all(size(firstBlock1) == size(x)))
  style = 's'; % same
  % Preallocate output.
  if (isa(firstBlock1, 'uint8'))
      b1 = repmat(uint8(0), ma+mpad, na+npad);
  else
      b1 = zeros(ma+mpad, na+npad);
  end
  if (isa(firstBlock2, 'uint8'))
      b2 = repmat(uint8(0), ma+mpad, na+npad);
  else
      b2 = zeros(ma+mpad, na+npad);
  end

  brows = 1:block(1);
  bcols = 1:block(2);
  mb = block(1);
  nb = block(2);
  xrows = brows + border(1);
  xcols = bcols + border(2);
  b1(brows, bcols) = firstBlock1(xrows, xcols);
  b2(brows, bcols) = firstBlock2(xrows, xcols);

elseif (all(size(firstBlock1) == (size(x)-2*border)))
  style = 'b'; % border
  % Preallocate output.
  if (isa(firstBlock1, 'uint8'))
      b1 = repmat(uint8(0), ma+mpad, na+npad);
  else
      b1 = zeros(ma+mpad, na+npad);
  end
  if (isa(firstBlock2, 'uint8'))
      b2 = repmat(uint8(0), ma+mpad, na+npad);
  else
      b2 = zeros(ma+mpad, na+npad);
  end
  brows = 1:block(1);
  bcols = 1:block(2);
  b1(brows, bcols) = firstBlock1;
  b2(brows, bcols) = firstBlock2;
  mb = block(1);
  nb = block(2);
else
  style = 'd'; % different
  [P,Q] = size(firstBlock1);
  brows = 1:P;
  bcols = 1:Q;
  mb = P;
  nb = Q;
  if (isa(firstBlock1, 'uint8'))
      b1 = repmat(uint8(0), mblocks*P, nblocks*Q);
  else
      b1 = zeros(mblocks*P, nblocks*Q);
  end
  if (isa(firstBlock2, 'uint8'))
      b2 = repmat(uint8(0), mblocks*P, nblocks*Q);
  else
      b2 = zeros(mblocks*P, nblocks*Q);
  end
  b1(brows, bcols) = firstBlock1;
  b2(brows, bcols) = firstBlock2;
end

[rr,cc] = meshgrid(0:(mblocks-1), 0:(nblocks-1));
rr = rr(:);
cc = cc(:);
mma = block(1);
nna = block(2);
for k = 2:length(rr)
  x = aa(rr(k)*mma+arows,cc(k)*nna+acols);
  [c1, c2] = feval(fun,x,params{:});
  if (style == 's')
    b1(rr(k)*mb+brows,cc(k)*nb+bcols) = c1(xrows,xcols);
    b2(rr(k)*mb+brows,cc(k)*nb+bcols) = c2(xrows,xcols);
  elseif (style == 'b')
    b1(rr(k)*mb+brows,cc(k)*nb+bcols) = c1;
    b2(rr(k)*mb+brows,cc(k)*nb+bcols) = c2;
  elseif (style == 'd')
    b1(rr(k)*mb+brows,cc(k)*nb+bcols) = c1;
    b2(rr(k)*mb+brows,cc(k)*nb+bcols) = c2;
  end
end

if ((style == 's') | (style == 'b'))
  b1 = b1(1:ma,1:na);
  b2 = b2(1:ma,1:na);
end

%%%
%%% Function parse_inputs
%%%
function [a, block, border, fun, params, padval] = parse_inputs(varargin)

switch nargin
case 0
    error('Too few inputs to BLKPROC');
    
case 1
    error('Too few inputs to BLKPROC');

case 2
    error('Too few inputs to BLKPROC');

case 3
    % BLKPROC(A, [m n], 'fun')
    a = varargin{1};
    block = varargin{2};
    border = [0 0];
    fun = varargin{3};
    params = cell(0,0);
    padval = 0;
    
case 4
    if (strcmp(varargin{2}, 'indexed'))
        % BLKPROC(X, 'indexed', [m n], 'fun')
        a = varargin{1};
        block = varargin{3};
        border = [0 0];
        fun = varargin{4};
        params = cell(0,0);
        padval = 1;
        
    elseif (isstr(varargin{3}) | isa(varargin{3},'inline'))
        % BLKPROC(A, [m n], 'fun', P1)
        a = varargin{1};
        block = varargin{2};
        border = [0 0];
        fun = varargin{3};
        params = varargin(4);
        padval = 0;
        
    else
        % BLKPROC(A, [m n], [mb nb], 'fun')
        a = varargin{1};
        block = varargin{2};
        border = varargin{3};
        fun = varargin{4};
        params = cell(0,0);
        padval = 0;
        
    end
    
otherwise
    if (strcmp(varargin{2}, 'indexed'))
        if (isstr(varargin{4}) | isa(varargin{4}))
            % BLKPROC(A, 'indexed', [m n], 'fun', P1, ...)
            a = varargin{1};
            block = varargin{3};
            border = [0 0];
            fun = varargin{4};
            params = varargin(5:end);
            padval = 1;
            
        else
            % BLKPROC(A, 'indexed', [m n], [mb nb], 'fun', P1, ...)
            a = varargin{1};
            block = varargin{3};
            border = varargin{4};
            fun = varargin{5};
            params = varargin(6:end);
            padval = 1;
            
        end
        
    else
        if (isstr(varargin{3}) | isa(varargin{3},'inline'))
            % BLKPROC(A, [m n], 'fun', P1, ...)
            a = varargin{1};
            block = varargin{2};
            border = [0 0];
            fun = varargin{3};
            params = varargin(4:end);
            padval = 0;
            
        else
            % BLKPROC(A, [m n], [mb nb], 'fun', P1, ...)
            a = varargin{1};
            block = varargin{2};
            border = varargin{3};
            fun = varargin{4};
            params = varargin(5:end);
            padval = 0;
            
        end
        
    end
end
    
fun = fcnchk(fun, length(params));

if (isa(a,'uint8'))
    padval = 0;
end

