function blended = contrast_blend(img, N)
% contrast_blend
% INPUT
%   img - 4D matrix where each color image is stacked in descending order
%   of exposure time in the 4th dimension
%       for example: img(:,:,:,1) =  longest exposure time image
%                   img(:,:,:,2) = next longest exposure time image
%                   ...
%                   img(:,:,:,N) = shortest exposure time image
%   N - number of images to use in the matrix
% OUTPUT
%   blended - a contrast blended image

blended = img(:,:,:,1);
blended = reshape(blended, size(img1,1), size(img1,2), size(img1,3));
for i = 2:N
    img1 = blended;
    img2 = reshape(img(:,:,:,i), size(img1,1), size(img1,2), size(img1,3));

    gray1 = sum(img1, 3);
    gray1 = gray1/max(gray1(:));

    img1 = double(img1);
    img2 = double(img2);

    %blur function
    sigma = 3;
    width = 3*sigma;
    support = -width:width;
    gauss2D = exp(-(support/sigma).^2/2);
    gauss2D = gauss2D/ sum(gauss2D);
    smooth = conv2(conv2(gray1, gauss2D, 'same'), gauss2D', 'same');
    gray2 = 1-smooth;

    %weighting function
    blended(:,:,1) = floor(smooth.*img2(:,:,1) + (1-smooth(:,:)).*img1(:,:,1));
    blended(:,:,2) = floor(smooth.*img2(:,:,2) + (1-smooth).*img1(:,:,2));
    blended(:,:,3) = floor(smooth.*img2(:,:,3) + (1-smooth).*img1(:,:,3));

end
blended = uint8(blended);
imshow(blended);
end
