function newM = restoreBoundary(pointMatrix)
% restore point set that was reduced by edge detection
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>

% The input points are assumed to be in the right order. It is assumed that
%   from straight horizontal, straight vertical or straight diagonal point
%   series only the first and last point was retained. In the final data,
%   two adjacent points must have a distance:
%     max(abs(x - xN), abs(y - yN)) = 1
%   and points must be added when this distance is larger.

distanceMatrix = max(...
    abs(pointMatrix(1, 1:end) - pointMatrix(1, [2:end 1])), ...
    abs(pointMatrix(2, 1:end) - pointMatrix(2, [2:end 1])));
insertIndexVector = find(distanceMatrix > 1);

% prepare filling points
fillIn = cell(1, length(insertIndexVector));
nAdd = 0;
for iInsert = 1:length(insertIndexVector);
    thisIndex = insertIndexVector(iInsert);
    nextIndex = mod(thisIndex, size(pointMatrix, 2)) + 1;
    xDist = pointMatrix(1, thisIndex) - pointMatrix(1, nextIndex);
    yDist = pointMatrix(2, thisIndex) - pointMatrix(2, nextIndex);
    
    % prepare filling - vertical line
    if yDist < 0
        yFill = (pointMatrix(2, thisIndex)+1):...
                (pointMatrix(2, nextIndex)-1);
    elseif yDist > 0
        yFill = (pointMatrix(2, thisIndex)-1):-1:...
                (pointMatrix(2, nextIndex)+1);
    end
    % horizontal line
    if xDist < 0
        xFill = (pointMatrix(1, thisIndex)+1):...
                (pointMatrix(1, nextIndex)-1);
    elseif xDist > 0
        xFill = (pointMatrix(1, thisIndex)-1):-1:...
                (pointMatrix(1, nextIndex)+1);
    end
    
    if xDist == 0
        % vertical line
        dbgType(iInsert) = 'v';
        dbgCount(iInsert) =  abs(yDist)-1;
        nAdd = nAdd + abs(yDist)-1;
        fillIn{iInsert} = [yFill*0 + pointMatrix(1, thisIndex); ...
                           yFill];
        
    elseif yDist == 0
        % horizontal line
        dbgType(iInsert) = 'h';
        dbgCount(iInsert) = abs(xDist)-1;
        nAdd = nAdd + abs(xDist)-1;
        fillIn{iInsert} = [xFill; ...
                           xFill*0 + pointMatrix(2, thisIndex)];
    elseif abs(xDist) == abs(yDist)
        % diagonal line
        dbgType(iInsert) = 'd';
        dbgCount(iInsert) = abs(xDist)-1;
        nAdd = nAdd + abs(xDist)-1;
        fillIn{iInsert} = [xFill; yFill];
%         if dbgCount(iInsert) > 2
%             disp([pointMatrix(:, thisIndex) fillIn{iInsert} pointMatrix(:,nextIndex)]);
%         end
    else
        % none from above.. ..this should be an error
        error('elk:boundary:restoreBoundary', ['The input points are not ' ...
            'continuously neighboured OR have only gaps in diagonal, ' ...
            'horizontal and vertical direction']);
    end
end
% fillIn

newM = nan(size(pointMatrix) + [0 nAdd]);
newM(:, 1:insertIndexVector(1)) = pointMatrix(:, 1:insertIndexVector(1));
thisNewIndex = insertIndexVector(1) + 1;
% thisOldIndex = 1;
for iInsert = 1:length(insertIndexVector)
    % insert stuff
    nAdd = size(fillIn{iInsert}, 2);
    newM(:, (thisNewIndex-1) + (1:nAdd)) = fillIn{iInsert};
    thisNewIndex = thisNewIndex + nAdd;
    % append normal stuff
    if iInsert < length(insertIndexVector)
        nAdd = insertIndexVector(iInsert + 1) - insertIndexVector(iInsert);
    else
        nAdd = size(pointMatrix, 2) - insertIndexVector(iInsert);
    end
    newM(:, (thisNewIndex-1) + (1:nAdd)) = ...
        pointMatrix(:, insertIndexVector(iInsert) + (1:nAdd));
    thisNewIndex = thisNewIndex + nAdd;
end

if any(isnan(newM(:)))
    error('elk:boundary:internalError', 'This should not happen.');
end

newM = unique(newM', 'rows')';