accumconncomps(a,b,fun)

function [anew,bnew] = accumconncomps(a,b,fun)

% Construct array with accumulation of connected components 
% 
%     [c,valnew] = accumconncomps(cc,val)
%     [c,valnew] = accumconncomps(cc,val,fun)
%     [c,valnew] = accumconncomps(cc)
%
% accumconncomps creates vectors by accumulating elements in val using
% connected components in cc. Connected components are subsequent,
% identical elements in the vector cc.
%
% The input vectors cc and val must have same size. The output array c 
% contains the values of each connected component 
% (e.g. cc=[1 1 2 2 2 1 2 3 3] returns c=[1 2 1 2 3]). valnew contains the
% aggregated values in val in each connected component.
%
% fun is a function handle that determines the aggregation mode (default:
% @sum). fun must be a function that takes a vector and returns a scalar 
% (e.g. @mean, @var, @(x) max(x)-min(x)).
%
% accumconncomps(subs) is an equal expression to 
% accumconncomps(subs,ones(size(subs)),@sum) and counts the number of
% elements in each connected component.
%
% Example 1:
%
% Sum the values in val according to the connected components in cc.
%
%     cc  = [1 1 2 1 1 1 3 3 3 4 2 2];
%     val = [2.3 1.2 5 3 2 5 3.2 4.5 2 2.2 1.2 2.2];
% 
%     [c,valnew] = accumconncomps(cc,val,@sum)
% 
%     c    =
% 
%          1     2     1     3     4     2
% 
%     valnew =
% 
%          3.5   5.0   10.0  9.7   2.2   3.4
%
%
% Example 2:
%
% What does a distribution of the lengths of connected components in a
% random sequence of zeros and ones look like?
%
%    cc  = round(rand(1000000,1));
%    [c,valnew] = accumconncomps(cc);
%    hist(valnew,20)
%
%
% See also: ACCUMARRAY
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]unibas.ch)
% Date: 10. Sept. 2008


% check input arguments
if nargin == 1
    b   = ones(size(a));
    fun = @sum;
elseif nargin == 2
    fun = @sum;
elseif nargin == 3
    % check if functionhandle is provided
    if ~isa(fun,'function_handle')
        error('the third input argument must be a function handle')
    end
else
    error('wrong number of input arguments')
end

% are vectors provided?
if ~isvector(a)
    error('val and subs must be vectors')
end

% do they have the same size?
siza = size(a);
sizb = size(b);
if siza~= sizb;
    error('val and subs must have same size')
end

% transpose if both vectors are row vectors
if siza(1)<siza(2)
    a = a(:);
    b = b(:);
    flagtranspose = true;
else
    flagtranspose = false;
end

% find beginnings of connected components
ad = [true; diff(a)~=0];

% assign new subs to the vector 
adc = cumsum(ad);

% nr of independent layers
uniqueLayers = adc(end);

% use accumarray to construct new vector
bnew = accumarray(adc,b,[uniqueLayers 1],fun);
anew = a(ad);

% transpose back if row vectors were provided
if flagtranspose
    anew = anew';
    bnew = bnew';
end