Efficient algorithm for plotting edges detect in a triangular mesh

I often need to represent geometries by parsing STL as triangular meshes. Often plotting the geometry as a surf or trisurf plot does not do well to describe certain features, especially when in my applications, the geometries can become extremely large. Often the surfaces are very simple but have distinct holes/cuts/edges that I'd like to highlight. To do this I just detect all the edges and plot the eges that aren't shared with other triangles.

In the following example, I detect the edges of a peaks surface and plot as a 3D line:

% Generate Triangular Mesh
N = 1e3;
[x,y] = meshgrid(1:N,1:N);
con = delaunay(x,y);
z = peaks(N);

% Identify Mesh Edges
M = [con(:,[1,2]) ; con(:,[2,3]); con(:,[3,1])];   % Find all edges in mesh, note internal edges are repeated
[u,~,n] = unique(sort(M,2),'rows'); % determine uniqueness of edges
counts = accumarray(n(:), 1);   % determine counts for each unique edge
O = u(counts==1,:); % extract edges that only occurred once
I = O(:,1:2)';
[x0,y0,z0] = deal(x(I),y(I),z(I));

% Plot Results
figure()
trisurf(con,x,y,z)

figure()
plot3(x0,y0,z0)

Aside from potentially being able to improve performance, I'd like to be able to plot the edge as a single object such that it could be labeled in the legend, and automatically detected if a surface contained several boundaries or holes.

Answer 1

[DGM]

If you are reading your data from an STL and can represent it as a triangulation object, then this can be simplified tremendously. Otherwise, you're going to have to assemble your unordered edge segments to form a continuous curve. Unless you do that, you're going to generate thousands of separate plot objects, which is horribly slow.

unzip peaks21.stl.zip % for the forum

% say you're reading an STL
T = stlread('peaks21.stl');

% get the free boundary edges
% depending on the version and scenario, this may be 
% an unordered list of edge segments, or it may be 
% ordered lists of edge segments describing a set of
% closed curves all concatenated together.
E = freeBoundary(T);

% reorder or split them as necessary to form
% a set of ordered closed curves.
boundaries = fixboundaries(E);

% display the surface using patch() or trisurf() or whatever
patch('faces',T.ConnectivityList,'vertices',T.Points, ...
    'facecolor','w','edgecolor','k');
view(3); camlight; view(-15,33)
axis equal; grid on; hold on
xlabel('X'); ylabel('Y'); zlabel('Z')

% display the boundaries
hp = gobjects(numel(boundaries),1);
for k = 1:numel(boundaries)
    thisb = T.Points(boundaries{k}([1:end 1],1),:);
    hp(1) = plot3(thisb(:,1),thisb(:,2),thisb(:,3), ...
        'b','linewidth',2);
end

This example only has one boundary, but the code is written to sort the boundary list and plot all the boundaries. There is only one plot object for each closed boundary curve.