function [ROIS] = dicomrt_rasterizeVoi(VOI, ct, ct_xmesh, ct_ymesh, ct_zmesh)
% See also: dicomrt_build3dVOI()

% Check input
% % % Commented for ease of downsample
% % % [ct_temp,type_ct,label,PatientPosition]=dicomrt_checkinput(ct);
% % % ct=dicomrt_varfilter(ct_temp);
[VOI_temp]=dicomrt_checkinput(VOI);
VOI=dicomrt_varfilter(VOI_temp);

% TODO check VOI has to be 2D and along z
% TODO(?) check x,y,zmesh has to be monotonic and unique

% assistant variables
xdim = length(ct_xmesh);
ydim = length(ct_ymesh);
zdim = length(ct_zmesh);

voi2use = 1:size(VOI,1);

[XMESH YMESH] = meshgrid(ct_xmesh, ct_ymesh);

% Progress bar data
h = waitbar(0, 'Rasterizing', 'Name', 'dicomrt_rasterizeVoi');
for kk = 1:length(voi2use)
    waitbar(kk./length(voi2use), h, ['Rasterizing: ', VOI{kk,1}]);
        
    % For each voi create a matrix with just 1 and 0
    rasterVOI = zeros(xdim, ydim, zdim, 'int8');
    num_curves = numel(VOI{kk,2});
    for k = 1:num_curves
        
%         if k == 6 %% manual fix for hippocampus
%             continue
%         end
        
        % the N x 3 curve matrix
        curve = VOI{kk,2}{k};
        % skip empty curves
        numPoints = size(curve,1);
        if numPoints < 3
            continue;
        end
        % z coordinate of current VOI
        voi_z = VOI{kk,2}{k}(1,3);
        % find the appropriate slice
        % TODO can simplify this once I'm more confident
        dist_z = abs(ct_zmesh-voi_z);
        slice = find(dist_z == min(dist_z));
%         slice = find(abs(ct_zmesh-voi_z) < dz);
        if isempty(slice)
            warning('Can not find matching slice');
            error('Closest slice found at %f cm', min(abs(ct_zmesh-voi_z)));
        elseif numel(slice) > 1
            warning('Multiple matching slices found');
%         elseif dist_z(slice) >= abs(ct_zmesh(1)-ct_zmesh(2)) % grozomah
%             error('Closest slice found at %f cm', min(abs(ct_zmesh-voi_z)));
        end
        % binary mask of current curve
        if exist('InPolygon', 'file') == 3
            % favor the mex function for performance
            BW = InPolygon(XMESH,YMESH,curve(:,1),curve(:,2));
        else
            BW = inpolygon(XMESH,YMESH,curve(:,1),curve(:,2));
        end
        
        % merge with existing curve
        rasterVOI(:,:,slice) = int8(xor(BW, rasterVOI(:,:,slice)));
    end
    
    % Create contour polygon on X and Y plane
    Xcurves = cell(0);
    for slice = 1:size(rasterVOI, 1)
        B = bwboundaries(squeeze(rasterVOI(slice,:,:)));
        if ~isempty(B)
            for ii = 1:numel(B)
                % Here x, y, z vertices means coordinate along 1st, 2nd and 3rd dimension,
                % respectively as in line Xcurves{end+1} = cat(2, x, y, z);
                x_vertices = ct_xmesh(slice) * ones(size(B{ii}, 1), 1);
                y_vertices = interp1(1:numel(ct_ymesh), ct_ymesh, B{ii}(:,1));
                z_vertices = interp1(1:numel(ct_zmesh), ct_zmesh, B{ii}(:,2));
                Xcurves{end+1} = cat(2, x_vertices, y_vertices, z_vertices);
            end
        end
    end
    
    Ycurves = cell(0);
    for slice = 1:size(rasterVOI, 2)
        B = bwboundaries(squeeze(rasterVOI(:,slice,:)));
        if ~isempty(B)
            for ii = 1:numel(B)
                y_vertices = ct_ymesh(slice) * ones(size(B{ii}, 1), 1);
                x_vertices = interp1(1:numel(ct_xmesh), ct_xmesh, B{ii}(:,1));
                z_vertices = interp1(1:numel(ct_zmesh), ct_zmesh, B{ii}(:,2));
                Ycurves{end+1} = cat(2, x_vertices, y_vertices, z_vertices);
            end
        end
    end
    
    % Save the ROIs struct
    % ROIS{kk}.BW = rasterVOI;
    ROIS{kk}.name = VOI{kk,1};
    ROIS{kk}.Xcurves = Xcurves;
    ROIS{kk}.Ycurves = Ycurves;
    % legacy nomenclature: num_curves and curves both refer to Z dimension
    ROIS{kk}.num_curves = num_curves;
    for k = 1:num_curves
        ROIS{kk}.curves{k} = VOI{kk,2}{k};
    end
    ROIS{kk}.ind = int32(find(rasterVOI));
    ROIS{kk}.dim = [xdim ydim zdim];
    ROIS{kk}.pixdim = [abs(ct_xmesh(2)-ct_xmesh(1)) abs(ct_ymesh(2)-ct_ymesh(1)) abs(ct_zmesh(2)-ct_zmesh(1))];
    ROIS{kk}.start = [ct_xmesh(1) ct_ymesh(1) ct_zmesh(1)];
    ROIS{kk}.start_ind = '(1,1,1)';
end
delete(h);