function c = subtractiveMix( a, s, bkgnd )
%c = subtractiveMix( a, s )
%   s is a parameter between 0 and 1.
%   a is an N*M*3 array

    cdim = length(size(a));
    mdim = cdim-1;
    if false && (nargin > 2)
        if cdim==3
            a(:,end+1,:) = permute( repmat( bkgnd, size(a,1), 1 ), [1,3,2] );
        else
            a(end+1,:) = bkgnd;
        end
    end
    
    
    min_a = min( a, [], cdim );
    a = 1 - (1-a)*(1-s);
    c = reshape( prod( a, mdim ), [], 3 );
    want_min_c = reshape( prod( min_a, mdim ), [], 1 );
    min_c = reshape( min( c, [], 2 ), [], 1 );
    adjustneeded = (min_c < 1) & (min_c > want_min_c);
    if any(adjustneeded)
        adjustments = (1-want_min_c(adjustneeded))./(1-min_c(adjustneeded));
        for i=1:3
            c(adjustneeded,i) = 1 - (1 - c(adjustneeded,i)) .* adjustments;
        end
    end
    
    if false && (nargin > 2)
        c = backgroundMix( c, bkgnd );
    end
end

function c = backgroundMix( c, bkgnd )
    opacity = max(c,[],2) - min(c,[],2);
    for i=1:3
        c(:,i) = opacity .* c(:,i) + (1-opacity) * bkgnd(i);
    end
end
