computeDistanceFromHrep2D

 compute the signed distance of points from polytope

 THIS IS NO USER FUNCTION

• convertCartToPolar convert to polar data and sort for increasing angles

• obtainShapeFromFunctional extract shape information from collective structure

0001 function distanceVector = computeDistanceFromHrep2D(hrep, vrep, pointMatrix)
0002 % compute the signed distance of points from polytope
0003 %
0004 % THIS IS NO USER FUNCTION
0005 
0006 % The elk-library: convex geometry applied to crystallization modeling.
0007 %   Copyright (C) 2013 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 %% catch empty hrep (origin is hrep)
0023 if isempty(hrep) || isempty(hrep.A)
0024     [tmp distanceVector] = convertCartToPolar(pointMatrix(1,:), pointMatrix(2,:));
0025     return
0026 end
0027 
0028 nEdge = size(hrep.A, 1);
0029 
0030 distanceMatrix = bsxfun(@minus, hrep.A * pointMatrix, hrep.h);
0031 [distanceVector edgeIndexVector] = max(distanceMatrix, [], 1);
0032 % negative or zero values are already appropriate because they imply that
0033 %   all points are contained in the polytope. However, if the points are
0034 %   outside the situation is more complex.
0035 % The current distance vector provides only the maximum distance from
0036 %   halfspaces and the corresponding closest point in this halfspace might
0037 %   not be in the polytope.
0038 % We can observe that the closest point of the polytope to a given point
0039 %   (xi) is
0040 %     1) on the edge with maximum distance, if it is perpendicular above
0041 %        the edge:
0042 %          xi <is in> gamma*edgeNormal + lam*edgePointOne ...
0043 %                     + (1-lam)*edgePointTwo
0044 %             <with> 0<=gamma, 0<=lam<=1
0045 %     2) on one of the vertices of the edge with maximum distance for which
0046 %        the adjacent edge angle must be considered:
0047 %          xi <is in> edgePoint + gammaOne*edgeNormal ...
0048 %                               + gammaTwo*adjacentEdgeNormal
0049 %             <with> 0<=gammaOne, 0<=gammaTwo
0050 %   It is important that no other edge must be considered. The points (xi)
0051 %   can be no closer to the polytope than the maximum distance that was
0052 %   already choosen. Therefore, facets with smaller distances are not
0053 %   actively constraining the closest point -- except the two adjacent
0054 %   edges that do define the two adjacent vertices.
0055 
0056 %% correlate vertices and edges
0057 % the result of this step is knowing which vertices are situated clockwise
0058 %   to each edge and which anticlockwise to each edge.
0059 
0060 [edgeAngleVector ~] = convertCartToPolar(hrep.A(:,1)', hrep.A(:,2)');
0061 [tmp sortedEdgeIndexVector] = sort(edgeAngleVector, 2, 'descend');
0062 % the incidence matrix tells in every column (=vertex) to which edges it
0063 %   belongs to:
0064 incidenceMatrix = abs(bsxfun(@minus, hrep.A * vrep.V', hrep.h)) < elkZerotol;
0065 clockwiseVertexMatrix = nan(2, length(edgeAngleVector));
0066 for iEdge = 1:nEdge
0067     thisEdgeIndex = sortedEdgeIndexVector(iEdge);
0068     nextEdgeIndex = sortedEdgeIndexVector(mod(iEdge, nEdge) + 1);
0069     thisVertexIndex = find(incidenceMatrix(thisEdgeIndex, :));
0070     if incidenceMatrix(nextEdgeIndex, thisVertexIndex(1))
0071         clockwiseVertexMatrix(:, thisEdgeIndex) = vrep.V(thisVertexIndex(1), :)';
0072     elseif incidenceMatrix(nextEdgeIndex, thisVertexIndex(2))
0073         clockwiseVertexMatrix(:, thisEdgeIndex) = vrep.V(thisVertexIndex(2), :)';
0074     else
0075         error('elk:boundary:internalError', ['Subsequent vertex not found.' ...
0076             'This should not happen.']);
0077     end
0078 end
0079 counterclockwiseVertexMatrix = circshift(clockwiseVertexMatrix, [0 -1]);
0080 
0081 %% clockwise and anticlockwise boundary normals and distances
0082 normalMatrix = hrep.A';
0083 clockwiseNormalMatrix = [0 1;-1 0] * normalMatrix;
0084 counterclockwiseNormalMatrix = -1*clockwiseNormalMatrix;
0085 clockwiseDistanceVector = sum(...
0086     clockwiseNormalMatrix .* clockwiseVertexMatrix, 1);
0087 counterclockwiseDistanceVector = sum(...
0088     counterclockwiseNormalMatrix .* counterclockwiseVertexMatrix, 1);
0089 
0090 for iSample = find(distanceVector > 0)
0091     if  clockwiseNormalMatrix(:, edgeIndexVector(iSample))'...
0092        *pointMatrix(:, iSample) > clockwiseDistanceVector(edgeIndexVector(iSample))'
0093         % clockwise vertex is closest
0094         distanceVector(iSample) = norm(pointMatrix(:, iSample) - ...
0095             clockwiseVertexMatrix(:, edgeIndexVector(iSample)));
0096     elseif counterclockwiseNormalMatrix(:, edgeIndexVector(iSample))'...
0097        *pointMatrix(:, iSample) > counterclockwiseDistanceVector(edgeIndexVector(iSample))'
0098         % anticlockwise vertex is closest
0099         distanceVector(iSample) = norm(pointMatrix(:, iSample) - ...
0100             counterclockwiseVertexMatrix(:, edgeIndexVector(iSample)));
0101     else
0102         % entered distance is already correct - nothing to do
0103     end
0104 end