function boundary_edge = triangulation_node_to_boundary_edge ( triangle_node )

%*****************************************************************************80
%
%% triangulation_node_to_boundary_edge() returns triangulation boundary edges.
%
%  Discussion:
%
%    This routine is given a triangulation, an abstract list of triples
%    of nodes.  It is assumed that the nodes in each triangle are listed
%    in a counterclockwise order, although the routine should work 
%    if the nodes are consistently listed in a clockwise order as well.
%
%    It is assumed that each edge of the triangulation is either 
%    * an INTERIOR edge, which is listed twice, once with positive
%      orientation and once with negative orientation, or;
%    * a BOUNDARY edge, which will occur only once.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 September 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer TRIANGLE_NODE(TRIANGLE_NUM,3), the nodes that make up the
%    triangles.  These should be listed in counterclockwise order.
%
%  Output:
%
%    integer BOUNDARY_EDGE(BOUNDARY_EDGE_NUM,2), the pairs of nodes
%    that make up the boundary edges.
%
  triangle_num = size ( triangle_node, 1 );
%
%  Create a list of all edges.
%
  edge_num = 3 * triangle_num;
  edge = zeros ( edge_num, 3 );

  edge(               1:  triangle_num,1:2) = triangle_node(1:triangle_num,1:2);
  edge(  triangle_num+1:2*triangle_num,1:2) = triangle_node(1:triangle_num,2:3);
  edge(2*triangle_num+1:3*triangle_num,1)   = triangle_node(1:triangle_num,3);
  edge(2*triangle_num+1:3*triangle_num,2)   = triangle_node(1:triangle_num,1);
%
%  For sorting, we need to reorder some edges.  
%  But for recovery later, we need the original ordering.
%  So add a third item which is 1 or 0 depending on whether we reordered.
%
  edge(1:3*triangle_num,3) = ( edge(1:3*triangle_num,1) < edge(1:3*triangle_num,2) );
%
%  In each row, force the smaller entry to appear first.
%  The peculiar syntax is required because min() prefers to 
%  minimize along a column.
%
  e1(1:edge_num) = min ( edge(1:edge_num,1:2), [], 2 );
  e2(1:edge_num) = max ( edge(1:edge_num,1:2), [], 2 );

  edge(1:edge_num,1) = e1(1:edge_num);
  edge(1:edge_num,2) = e2(1:edge_num);
%
%  Now ascending sort the array.
%
  edge = sortrows ( edge );
%
%  The boundary edges are the elements that occur just once in EDGE.
%
  boundary_edge = zeros ( edge_num, 2 );

  e = 0;
  be = 0;
  
  while ( true )
%
%  Only one edge left.
%
    if ( e == edge_num - 1 )
      e = e + 1;
      be = be + 1;
      if ( edge(e,3) )
        boundary_edge(be,1:2) = edge(e,1:2);
      else
        boundary_edge(be,1) = edge(e,2);
        boundary_edge(be,2) = edge(e,1);
      end
      break;
    end
%
%  We can compare edges E+1 and E+2.
%
    if ( edge(e+1,1:2) == edge(e+2,1:2) )
      e = e + 2;
    else
      e = e + 1;
      be = be + 1;
      if ( edge(e,3) )
        boundary_edge(be,1:2) = edge(e,1:2);
      else
        boundary_edge(be,1) = edge(e,2);
        boundary_edge(be,2) = edge(e,1);
      end
    end
    
  end
%
%  Reduce the storage for BOUNDARY_EDGE.
%
  boundary_edge = boundary_edge(1:be,1:2);

  return
end

