obtainShapeFromFunctional

 extract shape information from collective structure

 This does also include roughness information.

 THIS IS NO USER FUNCTION

• computeDistanceFromHrep2D compute the signed distance of points from polytope
• computeNormalFromAngle get normal vetor into direction
• convertCartToPolar convert to polar data and sort for increasing angles

• fitShapeBnd fit a rounded shape to boundary data

0001 function shape = obtainShapeFromFunctional(collectiveStruct, bnd, lambda, ...
0002     vertexIndex)
0003 % extract shape information from collective structure
0004 %
0005 % This does also include roughness information.
0006 %
0007 % THIS IS NO USER FUNCTION
0008 
0009 % The elk-library: convex geometry applied to crystallization modeling.
0010 %   Copyright (C) 2013 Alexander Reinhold
0011 %
0012 % This program is free software: you can redistribute it and/or modify it
0013 %   under the terms of the GNU General Public License as published by the
0014 %   Free Software Foundation, either version 3 of the License, or (at your
0015 %   option) any later version.
0016 %
0017 % This program is distributed in the hope that it will be useful, but
0018 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0019 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0020 %   General Public License for more details.
0021 %
0022 % You should have received a copy of the GNU General Public License along
0023 %   with this program.  If not, see <http://www.gnu.org/licenses/>.
0024 
0025 %% identify precalculated inner base shape
0026 lambdaVector = [collectiveStruct(:).lambda];
0027 baseIndex = find(lambdaVector <= lambda, 1, 'last');
0028 
0029 %% shape information
0030 baseShape = collectiveStruct(baseIndex);
0031 shape.hrep = struct('A', baseShape.hrep.A, 'h', baseShape.hrep.h - ...
0032     (lambda - baseShape.lambda));
0033 shape.lambda = lambda;
0034 shape.edgeLengthVector = baseShape.edgeLengthFunctional(shape.hrep.h);
0035 shape.vrep = baseShape.vrepFunctional(shape.hrep.h);
0036 
0037 
0038 if ~exist('vertexIndex', 'var') || isempty(vertexIndex) || vertexIndex == 0
0039     %% overall distance
0040     shape.distanceVector = computeDistanceFromHrep2D(shape.hrep, shape.vrep, ...
0041         [bnd.cartX; bnd.cartY]);
0042 else
0043     %% evaluating a specific edge
0044     % In this case, we only consider points that belong to a specific
0045     %   vertex of the initial shape. In doing so, we are able to fit this
0046     %   particular vertex.
0047     % The core is the proper selection of points for roughness calculation.
0048     %   Using too few point, especially when the amount increases later on,
0049     %   is a bad choice because roughness will be initially extremely
0050     %   small. A good choice is to identify the sample point right above
0051     %   the vertex, meaning: close to the angular bisecting line of the
0052     %   original vertex. Because the initial vertex vanishes with
0053     %   increasing lambda, this angular bisector is always shited to the
0054     %   currently supporting vertex. We then consider simply the 10% sample
0055     %   points around.
0056     % Each vertex itself features only a limited portion of the whole
0057     %   shape, given the opening angle of the adjacent edges. This opening
0058     %   angle is tracked, also for later evaluation (weighting of the
0059     %   obtained lambda values).
0060     
0061     % The original edges creating the vertex are, with index and normal
0062     %   vectors:
0063     originalEdgeFilter = collectiveStruct(1).incidenceMatrix(:, vertexIndex);
0064     originalEdgeNormalMatrix = collectiveStruct(1).hrep.A(originalEdgeFilter, :)';
0065     
0066     % we generate some angle information: angle of normals, opening angle
0067     %   and angle of the angle bisector (meanAngle)
0068     normalAngleVector = convertCartToPolar(originalEdgeNormalMatrix(1,:), ...
0069         originalEdgeNormalMatrix(2, :));
0070     if diff(normalAngleVector) < pi
0071         meanAngle = mean(normalAngleVector);
0072         shape.angleRatio = (max(normalAngleVector) ...
0073                           - min(normalAngleVector)) / 2 / pi;
0074     else
0075         meanAngle = mod(0.5*max(normalAngleVector) + ...
0076                         0.5*(min(normalAngleVector)+2*pi), 2*pi);
0077         shape.angleRatio = (min(normalAngleVector) ...
0078                           - max(normalAngleVector) + 2*pi) / 2 / pi;
0079     end
0080     
0081     % The supporting vertex for these edges is computed, by:
0082     scalarProductMatrix = shape.vrep.V * originalEdgeNormalMatrix;
0083     vertexIndex = find(...
0084         (abs(scalarProductMatrix(:, 1) - max(scalarProductMatrix(:, 1))) < 1e-10) & ...
0085         (abs(scalarProductMatrix(:, 2) - max(scalarProductMatrix(:, 2))) < 1e-10) );
0086     if numel(vertexIndex) ~= 1
0087         error('elk:boundary:numericalError', ['Supporting vertex for ' ...
0088             'initially adjacent edges is not found. This should not happen.']);
0089     end
0090     supportingVertex = shape.vrep.V(vertexIndex, :)';
0091     
0092     % the polar angle of the bisecting point on the shape is:
0093     bisectPoint = supportingVertex + lambda*computeNormalFromAngle(meanAngle);
0094     bisectAngle = convertCartToPolar(bisectPoint(1), bisectPoint(2));
0095     
0096     % the shape point sitting on the angular bisecting line is:
0097     closenessVector = abs(bnd.polarAngle - bisectAngle);
0098     bisectIndex = find(closenessVector == min(closenessVector), 1);
0099     
0100     % these are the indices of the neighbouring points:
0101     halfRange = round(length(bnd.polarAngle)/10/2);
0102     sampleIndexVector = mod(...
0103         ((bisectIndex-halfRange):(bisectIndex+halfRange))-1, ...
0104         length(bnd.polarAngle)) + 1;
0105     
0106     % get distance from base hrep
0107     shape.distanceVector = computeDistanceFromHrep2D(shape.hrep, shape.vrep, ...
0108         [bnd.cartX(sampleIndexVector); bnd.cartY(sampleIndexVector)]);
0109 end
0110 
0111 %% roughness information from distances
0112 shape.roughness = max([0 (shape.distanceVector - lambda)]) ...
0113     + max([0 (lambda - shape.distanceVector)]);
0114 shape.squareRoughness = sqrt(max([0 (shape.distanceVector - lambda)])^2 ...
0115     + max([0 (lambda - shape.distanceVector)])^2);
0116 shape.meanSquare = mean((shape.distanceVector - lambda).^2);