% MATLAB script for generating in/out-gamut mask
% in = white, out = black

function generateGamutMask(img_name, mask_name, display_gamut)

if nargin <= 2
    display_gamut = false;
end

% read and convert the image to YyCxCz color space
Name = ['input/',img_name];
ycxcz_img = rgb2ycxcz4Visual(Name);

% declare the mask image
mask_img = zeros(size(ycxcz_img,1),size(ycxcz_img,2));

% load the printer gamut information
load data/IndigoGamut.mat;
printer_gamut = IndigoGamut;

% find the 3-D convex hull of the gamut
cvxh = convhulln(printer_gamut);

% display the convex hull of the gamut
if display_gamut
    figure();
    trisurf(cvxh,printer_gamut(:,2),printer_gamut(:,3),printer_gamut(:,1));
end

% center of the gamut
% it's just a single point pretty much guaranteed to be 
% inside the gamut of the printer
% all the points will be thought of as shifted so that
% the center becomes the origin (0,0,0)
center = [116/2;0;0];

% for each 3-point set of convex hull we find the corresponding plane
% using homogeneous coordinates representation
numpts = size(cvxh,1);
cvhull_planes = zeros(4,numpts);

for i = 1:numpts
    A = [printer_gamut(cvxh(i,:),:)-repmat(center',[3,1]),[1;1;1]];
    cvhull_planes(:,i) = null(A);
    cvhull_planes(:,i) = -cvhull_planes(:,i)/cvhull_planes(4,i);
    cvhull_planes(:,i) = cvhull_planes(:,i)/norm(cvhull_planes(1:3,i));
end

% we also find the inverse of the matrix constructed by putting together
% the three points each as a column vector
% this is used later
ivs = zeros(numpts,3,3);

for i = 1:numpts
    ivs(i,:,:) = reshape(...
        inv((printer_gamut(cvxh(i,:),:)-repmat(center',[3,1]))'),[1,3,3]);
end

% go through each pixel and determine if the color is inside or outside
% the printer gamut's convex hull
for i = 1:size(ycxcz_img,1)
    if rem(i,20) == 0
        fprintf('Processing %dth row\n',i);
    end
    for j = 1:size(ycxcz_img,2)
        curr_pt = reshape(ycxcz_img(i,j,:),[3,1])-center;
        
        % for each 3-point set of the convex hull
        % first represent the current point as a linear combination
        % of the three points
        % if the (elongated) current point intersects with the
        % plane inside the triangle constructed by the three points
        % all the coefficients for the linear combination will be positive
        % otherwise we conclude that the current point vector and the plane
        % intersects outside the triangle
        for l = 1:numpts
            iP = reshape(ivs(l,:,:),[3,3]);
            a = iP*curr_pt;
            if (a(1) >= 0) && (a(2) >= 0) && (a(3) >= 0)
                
                % if we conclude that the current vector meets the plane
                % inside the triangle, we compute the ratio between the 
                % the point of intersection and the current point
                % if the ratio is larger than 1 (i.e., the current point
                % needs to be stretched longer to meet the plane)
                % the current point is inside the convex hull
                % thus we mark the pixel as white
                alpha = -cvhull_planes(4,l)/...
                    (cvhull_planes(1:3,l)'*curr_pt);
                if alpha >= 1
                    mask_img(i,j) = 1;
                end
                break;
            end
        end
        
    end
end

% save the mask image
imwrite(uint8(255*mask_img),['output/',mask_name]);