triangulatePolygonPair

 Compute triangulation between a pair of 3D closed curves.

   [V, F] = triangulatePolygonPair(POLY1, POLY2)

   [V, F] = triangulatePolygonPair(..., 'recenter', FLAG)
   Where FLAG is a boolean, specifies whether the second curve should be
   translated to have the same centroid as the first curve. This can
   improve mathcing of vertices. Default is true.


   Example
     % triangulate a surface patch between two ellipses
     % create two sample curves
     poly1 = ellipseToPolygon([50 50 40 20 0], 36);
     poly2 = ellipseToPolygon([50 50 40 20 60], 36);
     poly1 = poly1(1:end-1,:);
     poly2 = poly2(1:end-1,:);
     % transform to 3D polygons / curves
     curve1 = [poly1 10*ones(size(poly1, 1), 1)];
     curve2 = [poly2 20*ones(size(poly2, 1), 1)];
     % draw as 3D curves
     figure(1); clf; hold on;
     drawPolygon3d(curve1, 'b'); drawPoint3d(curve1, 'bo');
     drawPolygon3d(curve2, 'g'); drawPoint3d(curve2, 'go');
     view(3); axis equal;
     [vertices, faces] = triangulatePolygonPair(curve1, curve2);
     % display the resulting mesh
     figure(2); clf; hold on;
     drawMesh(vertices, faces);
     drawPolygon3d(curve1, 'color', 'b', 'linewidth', 2);
     drawPolygon3d(curve2, 'color', 'g', 'linewidth', 2);
     view(3); axis equal;

   See also
     meshes3D, triangulateCurvePair, meshSurfaceArea

0001 function [vertices, faces] = triangulatePolygonPair(poly1, poly2, varargin)
0002 % Compute triangulation between a pair of 3D closed curves.
0003 %
0004 %   [V, F] = triangulatePolygonPair(POLY1, POLY2)
0005 %
0006 %   [V, F] = triangulatePolygonPair(..., 'recenter', FLAG)
0007 %   Where FLAG is a boolean, specifies whether the second curve should be
0008 %   translated to have the same centroid as the first curve. This can
0009 %   improve mathcing of vertices. Default is true.
0010 %
0011 %
0012 %   Example
0013 %     % triangulate a surface patch between two ellipses
0014 %     % create two sample curves
0015 %     poly1 = ellipseToPolygon([50 50 40 20 0], 36);
0016 %     poly2 = ellipseToPolygon([50 50 40 20 60], 36);
0017 %     poly1 = poly1(1:end-1,:);
0018 %     poly2 = poly2(1:end-1,:);
0019 %     % transform to 3D polygons / curves
0020 %     curve1 = [poly1 10*ones(size(poly1, 1), 1)];
0021 %     curve2 = [poly2 20*ones(size(poly2, 1), 1)];
0022 %     % draw as 3D curves
0023 %     figure(1); clf; hold on;
0024 %     drawPolygon3d(curve1, 'b'); drawPoint3d(curve1, 'bo');
0025 %     drawPolygon3d(curve2, 'g'); drawPoint3d(curve2, 'go');
0026 %     view(3); axis equal;
0027 %     [vertices, faces] = triangulatePolygonPair(curve1, curve2);
0028 %     % display the resulting mesh
0029 %     figure(2); clf; hold on;
0030 %     drawMesh(vertices, faces);
0031 %     drawPolygon3d(curve1, 'color', 'b', 'linewidth', 2);
0032 %     drawPolygon3d(curve2, 'color', 'g', 'linewidth', 2);
0033 %     view(3); axis equal;
0034 %
0035 %   See also
0036 %     meshes3D, triangulateCurvePair, meshSurfaceArea
0037 %
0038  
0039 % ------
0040 % Author: David Legland
0041 % e-mail: david.legland@inra.fr
0042 % Created: 2017-05-18,    using Matlab 9.1.0.441655 (R2016b)
0043 % Copyright 2017 INRA - Cepia Software Platform.
0044 
0045 
0046 %% Settings
0047 
0048 recenterFlag = true;
0049 while length(varargin) > 1
0050     pname = varargin{1};
0051     if strcmpi(pname, 'recenter')
0052         recenterFlag = varargin{2};
0053     else
0054         error('Unknown parameter name: %s', pname);
0055     end
0056     varargin(1:2) = [];
0057 end
0058 
0059 
0060 %% Memory allocation
0061 
0062 % concatenate vertex coordinates for creating mesh
0063 vertices = [poly1 ; poly2];
0064 
0065 % number of vertices on each polygon
0066 n1 = size(poly1, 1);
0067 n2 = size(poly2, 1);
0068 
0069 % allocate the array of facets (each edge of each polygon provides a facet)
0070 nFaces = n1 + n2;
0071 faces = zeros(nFaces, 3);
0072 
0073 
0074 % Translate the second polygon such that the centroids of the bounding
0075 % boxes coincide. This is expected to improve the matching of the two
0076 % curves.
0077 if recenterFlag
0078     box1 = boundingBox3d(poly1);
0079     box2 = boundingBox3d(poly2);
0080     center1 = (box1(2:2:end) + box1(1:2:end-1)) / 2;
0081     center2 = (box2(2:2:end) + box2(1:2:end-1)) / 2;
0082     vecTrans = center1 - center2;
0083     trans = createTranslation3d(vecTrans);
0084     poly2 = transformPoint3d(poly2, trans);
0085 end
0086 
0087 
0088 %% Init iteration
0089 
0090 % find the pair of points with smallest distance.
0091 % This will be the current diagonal.
0092 [dists, inds] = minDistancePoints(poly1, poly2);
0093 [dummy, ind1] = min(dists); %#ok<ASGLU>
0094 ind2 = inds(ind1);
0095 
0096 % consider two consecutive vertices on each polygon
0097 currentIndex1 = ind1;
0098 currentIndex2 = ind2;
0099 
0100 
0101 %% Main iteration
0102 % For each diagonal, consider the two possible facets (one for each 'next'
0103 % vertex on each polygon), each create current facet according to the
0104 % closest one.
0105 % Then update current diagonal for next iteration.
0106 
0107 for iFace = 1:nFaces
0108     nextIndex1 = mod(currentIndex1, n1) + 1;
0109     nextIndex2 = mod(currentIndex2, n2) + 1;
0110     
0111     % compute lengths of diagonals
0112     dist1 = distancePoints(poly1(currentIndex1, :), poly2(nextIndex2,:));
0113     dist2 = distancePoints(poly1(nextIndex1, :), poly2(currentIndex2,:));
0114     
0115     if dist1 < dist2
0116         % keep current vertex of curve1, use next vertex on curve2
0117         face = [currentIndex1 currentIndex2+n1 nextIndex2+n1];
0118         currentIndex2 = nextIndex2;
0119     else
0120         % keep current vertex of curve2, use next vertex on curve1
0121         face = [currentIndex1 currentIndex2+n1 nextIndex1];
0122         currentIndex1 = nextIndex1;
0123     end
0124     
0125     % create the facet
0126     faces(iFace, :) = face;
0127 end
0128