restoreBoundary

 restore point set that was reduced by edge detection

 THIS IS NO USER FUNCTION

• obtainBnd create boundary data structure from raw data

0001 function newM = restoreBoundary(pointMatrix)
0002 % restore point set that was reduced by edge detection
0003 %
0004 % THIS IS NO USER FUNCTION
0005 
0006 % The elk-library: convex geometry applied to crystallization modeling.
0007 %   Copyright (C) 2012 Alexander Reinhold
0008 %
0009 % This program is free software: you can redistribute it and/or modify it
0010 %   under the terms of the GNU General Public License as published by the
0011 %   Free Software Foundation, either version 3 of the License, or (at your
0012 %   option) any later version.
0013 %
0014 % This program is distributed in the hope that it will be useful, but
0015 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0016 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0017 %   General Public License for more details.
0018 %
0019 % You should have received a copy of the GNU General Public License along
0020 %   with this program.  If not, see <http://www.gnu.org/licenses/>
0021 
0022 % The input points are assumed to be in the right order. It is assumed that
0023 %   from straight horizontal, straight vertical or straight diagonal point
0024 %   series only the first and last point was retained. In the final data,
0025 %   two adjacent points must have a distance:
0026 %     max(abs(x - xN), abs(y - yN)) = 1
0027 %   and points must be added when this distance is larger.
0028 
0029 distanceMatrix = max(...
0030     abs(pointMatrix(1, 1:end) - pointMatrix(1, [2:end 1])), ...
0031     abs(pointMatrix(2, 1:end) - pointMatrix(2, [2:end 1])));
0032 insertIndexVector = find(distanceMatrix > 1);
0033 
0034 % prepare filling points
0035 fillIn = cell(1, length(insertIndexVector));
0036 nAdd = 0;
0037 for iInsert = 1:length(insertIndexVector);
0038     thisIndex = insertIndexVector(iInsert);
0039     nextIndex = mod(thisIndex, size(pointMatrix, 2)) + 1;
0040     xDist = pointMatrix(1, thisIndex) - pointMatrix(1, nextIndex);
0041     yDist = pointMatrix(2, thisIndex) - pointMatrix(2, nextIndex);
0042     
0043     % prepare filling - vertical line
0044     if yDist < 0
0045         yFill = (pointMatrix(2, thisIndex)+1):...
0046                 (pointMatrix(2, nextIndex)-1);
0047     elseif yDist > 0
0048         yFill = (pointMatrix(2, thisIndex)-1):-1:...
0049                 (pointMatrix(2, nextIndex)+1);
0050     end
0051     % horizontal line
0052     if xDist < 0
0053         xFill = (pointMatrix(1, thisIndex)+1):...
0054                 (pointMatrix(1, nextIndex)-1);
0055     elseif xDist > 0
0056         xFill = (pointMatrix(1, thisIndex)-1):-1:...
0057                 (pointMatrix(1, nextIndex)+1);
0058     end
0059     
0060     if xDist == 0
0061         % vertical line
0062         dbgType(iInsert) = 'v';
0063         dbgCount(iInsert) =  abs(yDist)-1;
0064         nAdd = nAdd + abs(yDist)-1;
0065         fillIn{iInsert} = [yFill*0 + pointMatrix(1, thisIndex); ...
0066                            yFill];
0067         
0068     elseif yDist == 0
0069         % horizontal line
0070         dbgType(iInsert) = 'h';
0071         dbgCount(iInsert) = abs(xDist)-1;
0072         nAdd = nAdd + abs(xDist)-1;
0073         fillIn{iInsert} = [xFill; ...
0074                            xFill*0 + pointMatrix(2, thisIndex)];
0075     elseif abs(xDist) == abs(yDist)
0076         % diagonal line
0077         dbgType(iInsert) = 'd';
0078         dbgCount(iInsert) = abs(xDist)-1;
0079         nAdd = nAdd + abs(xDist)-1;
0080         fillIn{iInsert} = [xFill; yFill];
0081 %         if dbgCount(iInsert) > 2
0082 %             disp([pointMatrix(:, thisIndex) fillIn{iInsert} pointMatrix(:,nextIndex)]);
0083 %         end
0084     else
0085         % none from above.. ..this should be an error
0086         error('elk:boundary:restoreBoundary', ['The input points are not ' ...
0087             'continuously neighboured OR have only gaps in diagonal, ' ...
0088             'horizontal and vertical direction']);
0089     end
0090 end
0091 % fillIn
0092 
0093 newM = nan(size(pointMatrix) + [0 nAdd]);
0094 newM(:, 1:insertIndexVector(1)) = pointMatrix(:, 1:insertIndexVector(1));
0095 thisNewIndex = insertIndexVector(1) + 1;
0096 % thisOldIndex = 1;
0097 for iInsert = 1:length(insertIndexVector)
0098     % insert stuff
0099     nAdd = size(fillIn{iInsert}, 2);
0100     newM(:, (thisNewIndex-1) + (1:nAdd)) = fillIn{iInsert};
0101     thisNewIndex = thisNewIndex + nAdd;
0102     % append normal stuff
0103     if iInsert < length(insertIndexVector)
0104         nAdd = insertIndexVector(iInsert + 1) - insertIndexVector(iInsert);
0105     else
0106         nAdd = size(pointMatrix, 2) - insertIndexVector(iInsert);
0107     end
0108     newM(:, (thisNewIndex-1) + (1:nAdd)) = ...
0109         pointMatrix(:, insertIndexVector(iInsert) + (1:nAdd));
0110     thisNewIndex = thisNewIndex + nAdd;
0111 end
0112 
0113 if any(isnan(newM(:)))
0114     error('elk:boundary:internalError', 'This should not happen.');
0115 end
0116 
0117 newM = unique(newM', 'rows')';