function [ColoredImage] = getColorExact(BinaryMarkupImage, CombinedImageYUV)

imageHeight = size(CombinedImageYUV, 1);
imageWidth  = size(CombinedImageYUV, 2);
imageSize   = imageHeight * imageWidth;

ColoredImage(:, :, 1) = CombinedImageYUV(:, :, 1); % Y channel is the same

linearIndices = reshape([1 : imageSize], imageHeight, imageWidth); % pixels' linear indices ('r' in the paper)
linearMarkedIndices = find(BinaryMarkupImage);                     % marked pixels' linear indices

neighborhoodRadius = 1; % pixels neighborhoodRadius away are considered inside the neighborhood

weightCount  = 0; % number of non-zero elements in matrix A which contains all the weights ('w_rs' in the paper)
elementCount = 0; % total number of elements in A (= imageWidth * imageHeight eventually)

% column and row indices of all the weights (i.e. columnIndices[0] = column of weight 0)
% since the neighborhood of each pixel can contain at most (2 * neighborhoodRadius + 1) ^ 2 pixels,
% the maximum number of indices is imageSize * (2 * neighborhoodRadius + 1) ^ 2
columnIndices = zeros(imageSize * (2 * neighborhoodRadius + 1) ^ 2, 1);
rowIndices    = zeros(imageSize * (2 * neighborhoodRadius + 1) ^ 2, 1);

weights = zeros(imageSize * (2 * neighborhoodRadius + 1) ^ 2, 1); % the weights themselves
localWeights = zeros(1, (2 * neighborhoodRadius + 1) ^ 2); % weights inside a neighborhood

% for all pixels (in column-first order)
for j = 1 : imageWidth
  for i = 1 : imageHeight
    elementCount = elementCount + 1;
    if (~BinaryMarkupImage(i, j)) % if this pixel is not marked
      localWeightCount = 0;
      % for all pixels in a neighborhood (neighborhoodRadius pixels away)
      for ii = max(1, i - neighborhoodRadius) : min(i + neighborhoodRadius, imageHeight)
        for jj = max(1, j - neighborhoodRadius) : min(j + neighborhoodRadius, imageWidth)
          if (ii ~= i) || (jj ~= j) % if it is not the current pixel
            weightCount = weightCount + 1;            
            rowIndices(weightCount)    = elementCount;          % row index is the linear index of the current pixel
            columnIndices(weightCount) = linearIndices(ii, jj); % column index is the linear index of the current neighbor
            localWeightCount = localWeightCount + 1;
            localWeights(localWeightCount) = CombinedImageYUV(ii, jj, 1); % weight is the Y value of the current neighbor (temporarily)
          end
        end
      end
      currentY = CombinedImageYUV(i, j, 1);
      localWeights(localWeightCount + 1) = currentY; % put the Y value of the current pixel here for convenience, will be discarded eventually      
      mean_r = mean(localWeights(1 : localWeightCount + 1)); % mean Y values in the neighborhood (M_r in the paper)
      variance = mean((localWeights(1 : localWeightCount + 1) - mean_r) .^ 2); % variance, or 'sigma^2' in the paper
      csig = variance * 0.6; % scale the variance to avoid color leaking and haloing artifacts
                             % very important parameter which gets no discussion in the paper
                             % 0.6 is to be used with equation (2), for equation (3), see below
      % uncomment the below if using equation (3), which gives slightly worse results
      %csig = variance * 0.3;
      if (csig < 0.000002) % avoid dividing by 0
        csig = 0.000002;
      end
      % compute the weights (equation (2) and (3) in the paper)
      localWeights(1 : localWeightCount) = exp(-(localWeights(1 : localWeightCount) - currentY) .^ 2 / csig); % equation (2)
      %localWeights(1 : localWeightCount) = 1.0 + (currentY - mean_r) * (localWeights(1 : localWeightCount) - mean_r) / (csig); % equation (3)
      % avoid very small weights which can lead to unstable solutions
      for iii = 1 : localWeightCount
        if localWeights(iii) < 0.01
          localWeights(iii) = 0.01;
        end
      end
      % normalize the weights so that they sum to 1
      localWeights(1 : localWeightCount) = localWeights(1 : localWeightCount) / sum(localWeights(1 : localWeightCount));
      % copy local weights to global weights
      % the weights are negated according to equation (1) in the paper
      weights(weightCount - localWeightCount + 1 : weightCount) = -localWeights(1 : localWeightCount);
    end
    % we're done with the neighbors
    % for marked pixels, all of its neighbors' weights are 0
    % now onto the current pixel itself    
    weightCount = weightCount + 1;
    rowIndices(weightCount) = elementCount; % row index is the linear index of the current pixel
    columnIndices(weightCount) = linearIndices(i, j); % column index is the linear index of the current pixel as well
    weights(weightCount) = 1; % the weight is 1 (from equation (1) in the paper)
  end
end
       
% shrink the arrays to the actual number of weights       
weights       = weights      (1 : weightCount);
columnIndices = columnIndices(1 : weightCount);
rowIndices    = rowIndices   (1 : weightCount);

% create the corresponding matrix and column vector
A = sparse(rowIndices, columnIndices, weights, elementCount, imageSize);
b = zeros(size(A, 1), 1);

for t = 2 : 3 % for U and V channels
  EachChannel = CombinedImageYUV(:, :, t);
  b(linearMarkedIndices) = EachChannel(linearMarkedIndices); % enforce the constraints on marked pixels
  x = A \ b; % find column vector x that minimizes |Ax - b| in the least square sense
  ColoredImage(:, :, t) = reshape(x, imageHeight, imageWidth, 1);
end

ColoredImage = ntsc2rgb(ColoredImage); % convert back to rgb