function [ imageO, alphaO ] = imblend ( imageA, alphaA, imageB, alphaB )
%IMBLEND Combine image A and B with alpha channels
% Usage:
%     [ X ] = imblend ( imageA, alphaA, imageB, alphaB )
% Input:
%     imageA, imageB = MxNx3 CDATA
%     alphaA, alphaB = scalar or MxN alpha map
%     imageA in front, imageB in the back
% Output:
%     X = X
%
% See also imagescX
%
% Copyleft 1999-, XMO.

% check input
if any(alphaA(:) < 0) || any(alphaA(:) > 1) || any(alphaB(:) < 0) || any(alphaB(:) > 1)
    error('alpha should be within [0 1]');
end

if size(imageA, 3) ~= 3 && size(imageB, 3) ~= 3
    error('unknown input image format');
end

if size(imageA) ~= size(imageB)
    error('image dimension mismatch');
end

if isscalar(alphaA)
elseif size(alphaA, 1) == size(imageA, 1) && size(alphaA, 2) == size(imageA, 2) && size(alphaA, 3) == 1
else
    error('inappropriate alphaA');
end
if isscalar(alphaB)
elseif size(alphaB, 1) == size(imageB, 1) && size(alphaB, 2) == size(imageB, 2) && size(alphaB, 3) == 1
else
    error('inappropriate alphaB');
end

% return directly if completely transparent
if all(alphaA(:) == 0) && all(alphaB(:) == 0)
    imageO = zeros(size(imageA));
    return;
end

% calculate alpha
alphaO = alphaA + alphaB .* (1 - alphaA);

% calculate image
imageO = zeros(size(imageA), class(imageA));
for i = 1:3
    imageO(:,:,i) = (imageA(:,:,i) .* alphaA + imageB(:,:,i) .* alphaB .* (1 - alphaA)) ./ alphaO;
end