triangulateCurvePair

TRIANGULATECURVEPAIR Compute triangulation between a pair of 3D open curves.

   [V, F] = testTriangulateCurvePair(CURVE1, CURVE2)

   Example
     % triangulate a surface patch between two open curves
     % create two sample curves
     t = linspace(0, 2*pi, 100);
     curve1 = [linspace(0, 10, 100)' 5 * sin(t') zeros(100,1)];
     curve2 = [linspace(0, 10, 100)' 2 * sin(t') 2*ones(100,1)];
     figure; hold on;
     drawPolyline3d(curve1, 'b');
     drawPolyline3d(curve2, 'g');
     view(3); axis equal;
     [v, f] = triangulateCurvePair(curve1, curve2);
     figure; drawMesh(v, f, 'linestyle', 'none');
     view(3); axis equal

   See also
     meshes3D, triangulatePolygonPair, meshSurfaceArea

0001 function [vertices, facets] = triangulateCurvePair(curve1, curve2)
0002 %TRIANGULATECURVEPAIR Compute triangulation between a pair of 3D open curves.
0003 %
0004 %   [V, F] = testTriangulateCurvePair(CURVE1, CURVE2)
0005 %
0006 %   Example
0007 %     % triangulate a surface patch between two open curves
0008 %     % create two sample curves
0009 %     t = linspace(0, 2*pi, 100);
0010 %     curve1 = [linspace(0, 10, 100)' 5 * sin(t') zeros(100,1)];
0011 %     curve2 = [linspace(0, 10, 100)' 2 * sin(t') 2*ones(100,1)];
0012 %     figure; hold on;
0013 %     drawPolyline3d(curve1, 'b');
0014 %     drawPolyline3d(curve2, 'g');
0015 %     view(3); axis equal;
0016 %     [v, f] = triangulateCurvePair(curve1, curve2);
0017 %     figure; drawMesh(v, f, 'linestyle', 'none');
0018 %     view(3); axis equal
0019 %
0020 %   See also
0021 %     meshes3D, triangulatePolygonPair, meshSurfaceArea
0022 
0023 % ------
0024 % Author: David Legland
0025 % E-mail: david.legland@inrae.fr
0026 % Created: 2017-05-18, using Matlab 9.1.0.441655 (R2016b)
0027 % Copyright 2017-2024 INRA - Cepia Software Platform
0028 
0029 %% Memory allocation
0030 
0031 % number of vertices on each curve
0032 n1 = size(curve1, 1);
0033 n2 = size(curve2, 1);
0034 
0035 % allocate the array of facets (each edge of each curve provides a facet)
0036 nFacets = n1 + n2 - 2;
0037 facets = zeros(nFacets, 3);
0038 
0039 % look for the closest ends of two curves and reverse the second curve if
0040 % needed
0041 p1 = curve1(1, :);
0042 if distancePoints(p1, curve2(1,:)) > distancePoints(p1, curve2(end,:))
0043     curve2 = curve2(end:-1:1,:);
0044 end
0045 currentIndex1 = 1;
0046 currentIndex2 = 1;
0047 
0048 % concatenate vertex coordinates for creating mesh
0049 vertices = [curve1 ; curve2];
0050 
0051 
0052 
0053 %% Main iteration
0054 % For each diagonal, consider the two possible facets (one for each 'next'
0055 % vertex on each curve), each create current facet according to the closest
0056 % one.
0057 % Then update current diagonal for next iteration.
0058 
0059 for i = 1:nFacets
0060     nextIndex1 = mod(currentIndex1, n1) + 1;
0061     nextIndex2 = mod(currentIndex2, n2) + 1;
0062     
0063     % compute lengths of diagonals
0064     dist1 = distancePoints(curve1(currentIndex1, :), curve2(nextIndex2,:));
0065     dist2 = distancePoints(curve1(nextIndex1, :), curve2(currentIndex2,:));
0066     
0067     if dist1 < dist2
0068         % keep current vertex of curve1, use next vertex on curve2
0069         facet = [currentIndex1 currentIndex2+n1 nextIndex2+n1];
0070         currentIndex2 = nextIndex2;
0071     else
0072         % keep current vertex of curve2, use next vertex on curve1
0073         facet = [currentIndex1 currentIndex2+n1 nextIndex1];
0074         currentIndex1 = nextIndex1;
0075     end
0076     
0077     % create the facet
0078     facets(i, :) = facet;
0079 end
0080