Silhouette outline of polygon surface

I'm trying to get the outline of a triangulated 3D surface projected onto the XY plane. Here's an example. First, let's get a triangulated faces/vertices surface into FV:

xyz = [5, 2, 6, 0, 3;  0, 2, 4, 2, 1;  5, 6, 6, 7, 9]';
tri = [1, 2, 4; 1, 5, 4; 1, 2, 5; 2, 3, 5; 2, 4, 3; 3, 4, 5] ;
FV = struct('faces',tri,'vertices',xyz);
figure, patch(FV,'FaceColor','c','FaceAlpha',0.95), hold on

Now, if this surface were convex, I could simply take the convex hull of the XY coordinates of all the outline outXY points:

k = convhull(FV.vertices(:,1:2));
outXY = FV.vertices(k,1:2);
plot(outXY(:,1),outXY(:,2),'.r-','LineWidth',6), axis equal

However, of course, I want to work with non-convex surfaces :) Notice the little patch of white where the convex hull didn't follow the "outline" of my blue patch?

I can imagine a relatively simple hack, where I would just make a solid colour patch similar to this figure, take a screenshot using getframe() then use bwboundaries() on a segmentation of the figure window image.

I'm just trying to work out if there's a nifty trick that would provide an exact result. Votes will be given for good ideas about the approach, even if no code is involved. For instance, the freeBoundary method of delaunay triangulations might be related here, but I can't work out how to reduce my 3d triangulation into a 2d triangulation (which will have overlapping triangles) in a way that freeBoundary() is actually useful.

Note, for the moment, let's just say that I'm not interested in internal holes, such that a doughnut surface would simply return the outline of the outer circle (similar to the 'noholes' option of bwboundaries).

Thanks, Sven.

[Cedric]

I can't run your code (no DELAUNAYTRIANGULATION, and struc-related error when using DELAUNAY).. you want something like what the following illustrates..

 x = [5, 2, 6, 0, 3] ;  y = [0, 2, 4, 2, 1] ;  z = [5, 6, 6, 7, 9] ;
 tri = [1, 2, 4; 1, 5, 4; 1, 2, 5; 2, 3, 5; 2, 4, 3; 3, 4, 5] ;
 trisurf(tri, x, y, z) ;
 z_XY = zeros(size(z)) ;
 hold on ;
 trisurf(tri, x, y, z_XY, 'FaceColor', 'k') ;

but you want the projection to be expressed through a minimal set of triangles, and not just a visual trick like what I did here (?)

[Sven]

Oooh, yeah, lucky me with a current license :) Sorry, the delaunayTriangulation() is a 2013a function. I'll find something to replace it. Actually, your example here is ideal and more "minimal" than my MWE... I'm gonna use it!

And yes, I'm afraid that I'm looking for the numbers rather than the visual trick. My hack would be to use the visual trick, then use getframe, then use bwboundaries.

[Cedric]

Well, I have no solution at this point but at least I understand the question now ;-) While crying on my 2012b license, I'll keep thinking a bit about how to find numbers!

[Sean de Wolski]

I have a few thoughts too, but they'll wait for the 1000mi plane ride tomorrow...

[Sven]

Well, I've got one solution which is a brute force union of all the triangles in the surface, using Polygon Clip, which is mexified so quite fast. I would venture to guess that there's a more optimal way to do this that doesn't involve touching all the triangles one at a time, so I'm still open to suggestions :) A good thing is that using Polygon Clip you get to track holes for free!

% Get a full list of triangular facets as a 3-by-2-by-N matrix
facets = FV.vertices(:,1:2)';
facets = permute(reshape(facets(:,FV.faces'), 2, 3, []),[2 1 3]);
% Turn each facet into a structure acceptable by PolygonClip
facetsCell = num2cell(facets,1);
facetsStruct = squeeze(struct('x',facetsCell(:,1,:),'y',facetsCell(:,2,:)));
% Start with the first facet, then UNIONIZE! (power to the people)
P = PolygonClip(facetsStruct(1),facetsStruct(2),3);
for i = 2:length(facetsStruct)
    P = PolygonClip(P,facetsStruct(i),3);
end

outXY = [P(1).x(:) P(1).y(:)];
outXY = outXY([1:end 1],:);
figure, patch(FV, 'FaceColor','g'), view(3), camlight, hold on
minZ = min(FV.vertices(:,3));
plot3(outXY(:,1), outXY(:,2), ones(size(outXY,1),1)*minZ, '.r-'), axis image

[Cedric]

Another brute-force solution that I had in mind was to compute intersects of all edges of the projection (which isn't too complicated), add these points to the set of vertices, and keep the envelope of the set of all points.

Answer 1

[Sean de Wolski]

Here is what I came up with:

%%Sean's try
% Here's the general approach: 
%
% # Calculate the Delaunay Triangulation in three dimensions to create
%   tetrahedra forming the shape.
% # Use the x and y coordinates of this triangulation to build two
%   dimensional polygons irrespective of Z.
% # Calculate the union of all of these two dimensional polygons to get the minimal
%   bounding polygon.
%
% 03/21/2013

%%Original Data

xyz = [5, 2, 6, 0, 3;  0, 2, 4, 2, 1;  5, 6, 6, 7, 9]';
tri = [1, 2, 4; 1, 5, 4; 1, 2, 5; 2, 3, 5; 2, 4, 3; 3, 4, 5] ;
FV = struct('faces',tri,'vertices',xyz);

%%Delaunay Triangulation

tri3  = triangulation(tri,xyz);

%%Engine:

% Disable warning since we don't need anything depending on CW v. CCW
warning('off','map:polygon:noExternalContours');

% Initialize final polygon vectors
px = [];
py = [];

% Unionize each polygon with the group using |polybool()| from Mapping
% Toolbox
for ii = 1:size(tri3.ConnectivityList,1)
    [px, py] = polybool('union',...
        tri3.Points(tri3.ConnectivityList(ii,:),1),...
        tri3.Points(tri3.ConnectivityList(ii,:),2),...
        px,py);        

end

%%Visualize and Compare:

figure;
patch(FV,'FaceColor','c','FaceAlpha',0.95);
line(px,py,zeros(size(px)),'color','r','LineWidth',3)
view(2); % How did we do?
title('Woohoo!');

[Sean de Wolski]

And I guess there is not really any reason to do the Delaunay triangulation calculations...

[Sven]

Yep, this looks like the same approach but it uses native (or at least toolboxed) polygon boolean-ing.

I suspect that it will be slightly faster than PolygonClip since (I assume) it's also quite optimised/mexified, but at least this way doesn't require packing and unpacking structures.

Nice answer Sean. By the way, doesn't it feel a little... unsatisfying... that we already knew a great deal about the connectedness of these triangles, thus somehow couldn't use this knowledge to, say, start from the the vertices on the convex hull and then traverse (somehow) along these vertices?

All I can think of at the moment is this:

1. Get the vertices on the convex hull (these must be part of the final set)
2. Get the sets of these CH vertices that are connected to each other (kind of like bwconncomp, but for faces/vertices)
3. Make a separate polygon from each of these sets
4. Any triangles that have all of their vertices ON or IN any of these polygons can be skipped in the boolean step (they won't contribute to the outline of the shadow)
5. Begin the boolean starting with these separate CH polygons and iteratively adding the rest of the triangles

[Sean de Wolski]

Thanks!

1) This is not necessarily true, if three or more points are collinear and on the convex hull.

2) This might be one possible optimization. I think you could use the freeBoundary method of the triangulation object to identify these triangles. I didn't do a time test but it's quite posisble that the overhead of doing the triangulation outweighs the boolean operations.

Also, it looks like you can feed nan-delimited vectors into polybool() which might allow for the removal of the loop.

4) This would basically mean using the output from freeBoundary to remove non-surface polygons.

5) Certainly possible if we could get to the above.

How many triangles do you have and how long is this taking? If it's just from a theoretical perspective that we can't do better, I agree it's a bit irking. What lead me to polybool was the thought that there might be some projection algorithm to do this but I couldn't find anything.

I think for research, we may want to look into z-buffering algorithms, since that's essentially what's going on here.

</ random stream of thoughts >

[Sean de Wolski]

Also, as far as having all of this knowledge about the structure, I was originally wondering about using this to do a constrained Delaunay triangulation forcing us to only get the triangles that make up this z-buffered shape.

After a playing with it and thinking about it more, I became more and more convinced that it would never work.

[Sven]

1) Yes, you're right. I considered this and thought it might simply choose non-connected CH points where connected points might exist... but I guess it might actually pick the "wrong" vertex that happens to still be connected to a different part of the CH that might screw things up... checking would be required.

2) True... but I couldn't wrap my head around how to use FB

3) You missed 3)

4) But the thing is, all vertices on a (manifold, not-screwed-up) 3D triangular patch are surface polygons. The fact that they may not be surface polygons in a 2D projection is something different, and not revealed by the freeBoundary() of the 3D triangulation. And to "do" the 2D triangulation, you can't just do a delaunayTriangulation of the projected vertices (since it will triangulate everything up to a convex hull even though the 3D patch may not have been convex).

5) Prossibly...

By now it's theoretical... boolean approach is sufficient for my application. I was reducing a big patch to 10000 vertices, and it was taking ~2secs using PolygonClip. I just realised that I have the Mapping Toolbox. I will try the nan-separated route, and report the timings later.

And you're right... every application of this I searched for online was dealing with rendering and Z-buffers.

My particular application may be interesting to you... imagine you've got patches of the skeleton and patches of some organs. Now imagine "punching" from different angles... I'm checking organ exposure to blunt trauma in pediatrics. Don't worry, purely research... no punching kids involved.