fit a rounded shape to boundary data
Syntax: result = fitShapeBnd(bnd)
This function takes a boundary curve (generated by obtainBnd) and fits a
rounded shape. The resulting structure contains the final result but
also the original boundary, characteristic values and intermediate
steps of the fitting process.
The returned structure contains (despite others) the fields:
lambda - the estimated additive sphericity in units of the input data
circumfence - circumfence in units of the input data
meanWidth - mean width (=mean Feret diameter) in units of the input
data
roughness - the maximal deviation of input points from the estimated
shape
shapeStruct - the estimated shape as a structure array for each edge
or vertex, respectively. It has the following fields:
>.edgeAngle - the angle of each edge in radian
>.edgeDistance - the distance of the edge from the origin
>.edgeCart - a matrix containing the start point (first column) and
end point (second column) of the edge
>.vertexOrigin - the center of the vertex circle
>.vertexRadius - the radius of the vertex radius
>.vertexStartAngle - start angle for vertex circle
>.vertexEndAngle - end angle for vertex circle
Optional parameters in the form fitShapeBnd(bnd, 'option', value) are the
following:
'tolAngle' - accuracy for edge angle in radian. Default: 2/180*pi
'minDiffAngle' - minimum difference for the angle of two edges.
Default:, 15/180*pi
'tolDistance' - accuracy for edge distance in pixel. Default: 1
'tolLambda' - accuracy for the fitted radius. Default: 1
'lambdaLocal' - switch for local (1) or global (0) approaches.
Default: 0
'lambdaError' - selected objective method with 'meanSq' for mean
squares and 'sqRough' for // meanSq
The following activate additional filters based on values of the
boundary data (see obtianBnd). Particles are not processed if the
filter criteria is violated.
'minConvexity' - minimum convexity ratio - Default: 0.9
'minDiameter' - minimum diameter - Default: 100
'maxSphericity' - maximum sphericity - Default: 0 (filter off)
'maxSphericity2' - maximum sphericity2 - Default: 0 (filter off)
'minBrightFraction' - minimum bright fraction (bubble detection).
Default: 0 (filter off)
'maxAngleGap' - maximum gap in polar angles of input data
(convexity). Default: 0 (filter off)
See also: obtainBnd, viewBoundary, viewBoundaryShape0001 function result = fitShapeBnd(bnd, varargin) 0002 % fit a rounded shape to boundary data 0003 % 0004 % Syntax: result = fitShapeBnd(bnd) 0005 % 0006 % This function takes a boundary curve (generated by obtainBnd) and fits a 0007 % rounded shape. The resulting structure contains the final result but 0008 % also the original boundary, characteristic values and intermediate 0009 % steps of the fitting process. 0010 % 0011 % The returned structure contains (despite others) the fields: 0012 % lambda - the estimated additive sphericity in units of the input data 0013 % circumfence - circumfence in units of the input data 0014 % meanWidth - mean width (=mean Feret diameter) in units of the input 0015 % data 0016 % roughness - the maximal deviation of input points from the estimated 0017 % shape 0018 % shapeStruct - the estimated shape as a structure array for each edge 0019 % or vertex, respectively. It has the following fields: 0020 % >.edgeAngle - the angle of each edge in radian 0021 % >.edgeDistance - the distance of the edge from the origin 0022 % >.edgeCart - a matrix containing the start point (first column) and 0023 % end point (second column) of the edge 0024 % >.vertexOrigin - the center of the vertex circle 0025 % >.vertexRadius - the radius of the vertex radius 0026 % >.vertexStartAngle - start angle for vertex circle 0027 % >.vertexEndAngle - end angle for vertex circle 0028 % 0029 % Optional parameters in the form fitShapeBnd(bnd, 'option', value) are the 0030 % following: 0031 % 'tolAngle' - accuracy for edge angle in radian. Default: 2/180*pi 0032 % 'minDiffAngle' - minimum difference for the angle of two edges. 0033 % Default:, 15/180*pi 0034 % 'tolDistance' - accuracy for edge distance in pixel. Default: 1 0035 % 'tolLambda' - accuracy for the fitted radius. Default: 1 0036 % 'lambdaLocal' - switch for local (1) or global (0) approaches. 0037 % Default: 0 0038 % 'lambdaError' - selected objective method with 'meanSq' for mean 0039 % squares and 'sqRough' for // meanSq 0040 % The following activate additional filters based on values of the 0041 % boundary data (see obtianBnd). Particles are not processed if the 0042 % filter criteria is violated. 0043 % 'minConvexity' - minimum convexity ratio - Default: 0.9 0044 % 'minDiameter' - minimum diameter - Default: 100 0045 % 'maxSphericity' - maximum sphericity - Default: 0 (filter off) 0046 % 'maxSphericity2' - maximum sphericity2 - Default: 0 (filter off) 0047 % 'minBrightFraction' - minimum bright fraction (bubble detection). 0048 % Default: 0 (filter off) 0049 % 'maxAngleGap' - maximum gap in polar angles of input data 0050 % (convexity). Default: 0 (filter off) 0051 % 0052 % See also: obtainBnd, viewBoundary, viewBoundaryShape 0053 0054 % The elk-library: convex geometry applied to crystallization modeling. 0055 % Copyright (C) 2013 Alexander Reinhold 0056 % 0057 % This program is free software: you can redistribute it and/or modify it 0058 % under the terms of the GNU General Public License as published by the 0059 % Free Software Foundation, either version 3 of the License, or (at your 0060 % option) any later version. 0061 % 0062 % This program is distributed in the hope that it will be useful, but 0063 % WITHOUT ANY WARRANTY; without even the implied warranty of 0064 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0065 % General Public License for more details. 0066 % 0067 % You should have received a copy of the GNU General Public License along 0068 % with this program. If not, see <http://www.gnu.org/licenses/>. 0069 0070 %! DEV: The tolerances were tweaked and denoted by new values: 0071 % new tolerances - angle: 2/180*pi, tolDistance 1, tolLambda 1 0072 % but the server images went missing so that it could not be published 0073 % with these values. The old varlues are: 0074 % old tolerances - angle 0.05; distance 2, tolLambda 2 0075 %/ 0076 optionStruct = mapOptionStruct(varargin, ... 0077 'tolAngle', 2/180*pi, ... 0078 'minDiffAngle', 15/180*pi, ... 0079 'tolDistance', 1, ... 0080 'tolLambda', 1, ... 0081 'lambdaLocal', 0, ... 0082 'lambdaError', 'meanSq', ... % sqHeight // meanSq 0083 'minConvexity', 0.9, ... 0084 'minDiameter', 100, ... 0085 'maxSphericity', 0, ... 0086 'maxSphericity2', 0, ... 0087 'minBrightFraction', 0, ... 0088 'maxAngleGap', 0, ... 0089 'debugHough', 0, ... 0090 'debugEdgeResult', 0, ... 0091 'debugLambda', 0); 0092 0093 %% ignore particles 0094 result = handleIgnoreFilter(bnd, optionStruct); 0095 if ~isempty(result) 0096 return 0097 end 0098 0099 %% identify edges 0100 if optionStruct.debugHough 0101 [hrep debug.edge.hough] = identifyEdges(bnd, optionStruct); 0102 else 0103 hrep = identifyEdges(bnd, optionStruct); 0104 end 0105 debug.edge.hrep = hrep; 0106 0107 %% debug code 0108 if optionStruct.debugEdgeResult 0109 figure,hold on; 0110 viewBoundary(bnd); 0111 viewPolytope(hrep, 'shade', 0, 'edgeColor', [1 0 0]); 0112 title('complete set of edges, lambda = 0') 0113 end 0114 0115 %% Precompute edge reduction 0116 % the following struct contains all information of subsequently reduced 0117 % edges for P(h-lambda). 0118 collectiveStruct = computeCollectiveShapeStruct(hrep, ... 0119 1e-5*optionStruct.tolDistance); 0120 % the next variable denotes where no inner polytope exists anymore: 0121 seperatingLambda = collectiveStruct(end).lambda; 0122 debug.lambda.collectiveStruct = collectiveStruct; 0123 0124 %% Identify lambda 0125 % The shape model is P(h) + lam*B and for P(h) the edge angles are already 0126 % determined by the previous step. The previous step also identified edge 0127 % distances so that we might rewrite the model as (P(hI-lam) + lam*B). 0128 % The roughness is defined by two limiting shapes. A lower bound shape: 0129 % (P(hI-lam) + lam*B) - gL*B 0130 % uses the maximum gL so that no sample point is contained in that shape. 0131 % And an upper bound shape: 0132 % (P(hI-lam) + lam*B) + gU*B 0133 % uses a minimum gU so that all points are contained in that shape. Note 0134 % that these decisions can be made based on the (signed) distance of 0135 % all points from the polytope P(hI-lam). Signed means using negative 0136 % values for inner points and positive values for outer points. 0137 % Rouhness is, however not an ideal objective function. Usually, the 0138 % roughness does not change if lam is chosen in the interval between 0139 % (lam-gL) and (lam+gU) so that both extreme cases are possible. To avoid 0140 % this, we use a squarred roughness: 0141 % squareRoughness = sqrt(gU^2 + gL^2) 0142 % so that lambda should lie in the middle of (lam-gL) and (lam+gU). 0143 0144 %% all vertices at once 0145 % This old code fragment fits the whole shape as one. It was found that 0146 % this method often fits to the smallest lambda that is featured by an 0147 % vertex. As soon as lambda exceeds this value, gU rises. Even if a lot 0148 % more vertices feature a larger lambda value, this will not be taken 0149 % into account. 0150 0151 % getObjectFunction = @(lam)(obtainShapeFromFunctional(... 0152 % collectiveStruct, bnd, lam)); 0153 % objectiveFunction = @(shape)(shape.squareRoughness); 0154 % 0155 % shape = lineSearchGolden(getObjectFunction, objectiveFunction, ... 0156 % 0, seperatingLambda, optionStruct.tolLambda); 0157 0158 if optionStruct.lambdaLocal 0159 %% seperate vertices 0160 % in this approach, we consider each vertex of the initial 0161 % H-representation. 0162 nVertex = size(collectiveStruct(1).vrep.V, 1); 0163 lambdaVector = nan(1, nVertex); 0164 angleRatioVector = nan(1,nVertex); 0165 for iVertex = 1:nVertex 0166 getObjectFunction = @(lam)(obtainShapeFromFunctional(... 0167 collectiveStruct, bnd, lam, iVertex)); 0168 if strcmpi(optionStruct.lambdaError, 'sqHeight') 0169 % Well, square heigt was named square roughness in the early 0170 % stages. 0171 objectiveFunction = @(shape)(shape.squareRoughness); 0172 elseif strcmpi(optionStruct.lambdaError, 'meanSq') 0173 objectiveFunction = @(shape)(shape.meanSquare); 0174 end 0175 0176 shape = lineSearchGolden(getObjectFunction, objectiveFunction, ... 0177 0, seperatingLambda, optionStruct.tolLambda); 0178 0179 lambdaVector(iVertex) = shape.lambda; 0180 % angleRatioVector(iVertex) = shape.angleRatio; 0181 angleRatioVector(iVertex) = 1/nVertex; 0182 debug.lambda.shape(iVertex) = shape; 0183 0184 %! Debug: 0185 % figure,hold on 0186 % h = viewBoundaryShape(collectiveStruct(1));set(h, 'color', [0 1 0]); 0187 % viewBoundary(bnd),viewBoundaryShape(shape) 0188 % plot(collectiveStruct(1).vrep.V(iVertex, 1), ... 0189 % collectiveStruct(1).vrep.V(iVertex, 2), 'm.', 'markerSize', 20); 0190 %/ 0191 end 0192 % we use a mean from the lambdas, weighted by the corresponding opening 0193 % angle: 0194 thisLambda = angleRatioVector * lambdaVector'; 0195 shape = obtainShapeFromFunctional(collectiveStruct, bnd, thisLambda); 0196 % safe local fitting information 0197 debug.lambdaVector = lambdaVector; 0198 debug.angleRatioVector = angleRatioVector; 0199 else 0200 %% global lambda approach 0201 getObjectFunction = @(lam)(obtainShapeFromFunctional(... 0202 collectiveStruct, bnd, lam, 0)); 0203 if strcmpi(optionStruct.lambdaError, 'sqRough') 0204 objectiveFunction = @(shape)(shape.squareRoughness); 0205 elseif strcmpi(optionStruct.lambdaError, 'meanSq') 0206 objectiveFunction = @(shape)(shape.meanSquare); 0207 end 0208 0209 shape = lineSearchGolden(getObjectFunction, objectiveFunction, ... 0210 0, seperatingLambda, optionStruct.tolLambda); 0211 end 0212 0213 result = mergeStruct(shape, bnd); 0214 result.debug = debug; 0215 result.flagIgnore = 0; 0216 result.circumfence = sum(result.edgeLengthVector) + 2*pi*result.lambda; 0217 result.meanWidth = result.circumfence / pi; 0218 result.shapeArea = computeVolumeVrep(result.vrep) + ... 0219 result.lambda*result.meanWidth + pi*result.lambda^2;