function [heartNormalField] = computeNormalFieldSingleFit(X,mask)
    % mask : 2D mask...
    % X: boundary points
    imageSize = zeros(1,3);
    imageSize(1:2) = size(mask);
    % this assums that there are boundary points for all Z slices.
    imageSize(3) = ceil(max(X(:,3)));
    heartNormalField = zeros([imageSize(1),imageSize(2),2,imageSize(3)]);
    tic
    [C,R] = fitCircle(X(:,1:2));
    for pt = 1:length(X)
        startPoint = X(pt,:);
        dr = [C(1) - startPoint(1), C(2) - startPoint(2)];
        drNorm = norm(dr);
        assert(drNorm > 0)
        dr = dr./drNorm;
        heartNormalField(startPoint(2),startPoint(1),:,startPoint(3)) = dr;
    end
    heartNormalField = squeeze(heartNormalField);
    fprintf('Computed single Fit heartNormalfield in %s seconds\n', duration(0,0,toc));
end