function [X, Y, Z] = computePatchDataFromPointCloud(dataMatrix)
% compute x, y, z data for patch plots of irregular grids
%
% THIS IS (NOT YET) A USER FUNTION

%% input handling
% even though the algorithm is only required for 2D smapling points, I
%   already implement it as "could be nD".
nDim = size(dataMatrix, 1) - 1;

%% tesselation
% the principle relatively simple: the data points are centered around zero
%   and then elevated into 3D by using z-coordinates that reside on a
%   paraboloid. The mesh, created by convhulln can then be projected back
%   to 2D and represents a tesselation of the input points.

% scaling
centerPoint = mean(dataMatrix(1:nDim, :), 2);
positionMatrix = bsxfun(@minus, dataMatrix(1:nDim, :), centerPoint);
scaleFactor = std(positionMatrix, [], 2);
positionMatrix = bsxfun(@times, positionMatrix, 1./scaleFactor);
% elevate dimension
positionMatrix = [positionMatrix; ...
                  sum(positionMatrix.^2, 1)];
% figure,plot3(positionMatrix(1,:), positionMatrix(2,:), positionMatrix(3,:), 'b.')
              
% triangulate
patchIndexMatrix = convhulln(positionMatrix', {'Qt', 'Qs'});
nPatch = size(patchIndexMatrix, 1);
convhullIndexVector = convhulln(positionMatrix(1:nDim, :)', {'Qt', 'Qs'});
convhullIndexVector = unique(convhullIndexVector(:));

X = nan(nDim+1, nPatch);
Y = X; Z = X;

for iPatch = 1:nPatch
    % we have no patch if ALL points belong to the convex hull points
    incidenceMatrix = bsxfun(@eq, patchIndexMatrix(iPatch, :), ...
        convhullIndexVector(:));
    if all(any(incidenceMatrix, 2))
        continue
    end
    
    X(:, iPatch) = dataMatrix(1, patchIndexMatrix(iPatch, :));
    Y(:, iPatch) = dataMatrix(2, patchIndexMatrix(iPatch, :));
    Z(:, iPatch) = dataMatrix(3, patchIndexMatrix(iPatch, :));
end
