Simple polygon triangulation
If all you need is a triangulation, there are two very simple algorithms. There are others that are much harder to implement, and have minimal benefit (other than nice asymptotic properties as nVertices becomes very large.
a) subdivision : at each step, find a “chord” (a connection between two vertices that lies entirely INSIDE the polygon) and use it to divide the polygon into two pieces. Recursively triangulate the two pieces.
b) ear-trimming: a special case of a) where only chords cutting off a single triangle on one side are considered. You look for v0,v1,v2 - consecutive vertices, where v0->v1->v2 is convex and v0->v2 is entirely INSIDE the polygon. The naive algorithm is O(n^4), but with care you can reduce this to O(n^2) - this is an exercise in Joe O’Rourke’s book Computational Geometry in C
For “slightly non-planar” triangles (and, in general, for triangles in/near a plane other than the x-y plane, fit a plane to the vertices (I like eigenvectors for this task) and map all the vertices to that plane (there are various ways to do this). Triangulate in that plane, and apply that triangulation to the original points. Define “slightly non-planar” to be any problem for which this works!
I always advise people to implement b) first - and then (after putting it into production) investigate “better” algorithms. Don’t bother to push out a “better” method until customers complain that b) is “too slow”.
Now…Delaunay Triangles. These are nice if you want to avoid skinny triangles - but for “slightly non-planar” polygons, there’s not that much benefit. I recommend leaving that for later. One simple algorithm (but notoriously plagued by numerical issues) starts with an arbitrary triangulation and then swaps edges to improve the triangulation). Most simple examples you will find laying about do not deal with arbitrary polygons, but instead triangulate the convex hull of the vertices. I think OpenSCAD would be satisfied with a “Delaunay-like” improvement of the triangulation, without insisting on the absolute “best” triangulation.
As for holes, that will depend on how you represent the polygons with holes. I suspect I might de-compose a polygon with holes into several polygons without holes (noting the cut edges), triangulate the several polygons, and then stitch them back together at the cut edges). Looks easy enough, but I haven’t ever actually implemented it. It may be that simple ear-clipping can be made to understand holes - I just don’t know, and (again) that might depend on the representation.
--
Kenneth Sloan
KennethRSloan@gmail.com
Vision is the art of seeing what is invisible to others.
On Mar 29, 2016, at 17:33 , Ronaldo rcmpersiano@gmail.com wrote:
Parkinbot wrote
good question ;-)
so you'd need a convexity test first.
-
If convexity is given, you can triangulate by either
a) using one corner as common point (this is trivial)
b) pull up a new triangle at one at the CW or CCW endpoint, until you're
done. This gives you two choices for each new triangle. If both choices
are valid, you get a possible recursion retrace point for sceanrios where
both choices fail. -
For a non-convex polygon, you can decompose the shape into a set convex
shapes and use 1a. So your question could be, how to decompose a
non-convex polygon into a set of convex polygons. -
I think alternatively you could also try 1b and find some additional
test for doing the right choice.
I delayed my answer to your suggestions because I needed to take a look on
the convex decomposition methods. From my search, I found that:
a. the convex decomposition of simple polygons is a task so hard as to
triangulate it; some methods triangulate the polygon to find its convex
decomposition;
b. it is not difficult to find simple polygons whose convex decomposition
is itself a triangulation; the case I have in hands and that stimulated me
to implement a simple polygon triangulation falls exactly in that category;
its boundary is composed by two Bezier curves and two line segments;
c. convex decomposition methods require, in general, some kind of data
structure to reduce their complexity and those data structure are likewise
hard to implement in OpenSCAD language.
The ear-cut methods seem to be easier to implement in OpenSCAD. The simplest
ones are sub-optimal in complexity but require just lists as data structure.
I implemented a very simple version of one of those methods and it seems
that, for less than a thousand vertices, the time to triangulate a polygon
is of order n^2. This seems fair to me by now.
A question which is not well addressed yet is how to deal with almost planar
polygonals. The ear-cut method requires that the planar polygon vertices are
ordered in a CCW circulation. In 3D space, there is no well defined
orientation for circulating a polygon. In a planar polygon in 3D, the normal
to the polygon plane may be used to induce an orientation to the vertex
circulation. If the polygon normal is well chosen, it induces the correct
CCW orientation. If the orientation is not CCW, it is sufficient to invert
the normal. A simple check of the orientation correctness is a first issue.
Almost planar polygonals pose more issues since no normal really exists. My
approach is to find a vector that reduces to the polygonal plane normal when
the polygonal vertices are in a plane. I call it a false normal. A least
square adjustment is an alternative but too expensive. The crude method I
tried it is not fool proof. It takes the barycenters of two disjoint subsets
of the vertices and do a cross product of the vectors from the starting
vertex and the barycenters. Unhappily, this cross product may be near zero
inducing errors later on.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16803.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
OpenSCAD mailing list
Discuss@lists.openscad.org
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Personally I find it impossible to implement some complex algorithm with
OpenSCAD first hand. I would always use a language with a rich library
background and a comfortable IDE with debugging facilities like Matlab and
port the solution to OpenSCAD after having it refactored into a more
functional formulation that still passes all tests.
When implementing my NACA airfoil library in order to freely sweep around
morphing airfoils in 3D I encountered the same triangulation problem and
started to think about how to triangulate non-convex shapes. In the end I
decided to exploit some known symmetry of my airfoil generator and use a
straight forward scheme doing the job without any test. This solution is of
course far from being general.
I think triangulation of a general well-formed (=non self-intersecting, non
singular) N-polygon can be one by a recursive scheme, which I tried to
explain as 1b. Here some more words about it:
We head for a triangulation without adding extra points. I don't have a
proof for it, but I'd guess it will always exist for well-formed N-polygon,
N>3.
Recursion starts with three points in sequence that may be connected into a
valid triangle inside the shape. You can ensure this with a simple test by
use of normals. I guess it will save some work to prescribe an orientation
for the point sequence, but orientation can also be detected by an extra
cycle.
To see, whether a triple is also valid with respect to all other points in
the N-sequence, we have to implement an in_triangle() predicate telling us
whether a given point is inside the triangle and test all remaining points.
So it might need cycling around a bit, until a first valid triple is found.
Cycling also indices ensures the first triangle can be named [1,2,3].
Now we have two possibilities to construct the next triangle [1,3,4] and
[N,1,2]. We test the first one against the remaining N-4 points and if it is
valid, we reenter the recursion scheme. If it is not valid or the recursion
fails, we test the second one. If it is valid we reenter the recursion
scheme. If not, or the recursion fails, we continue cycling until we find
the next valid triangle is found to start the recursion again. (Once the
cycle is done, a counter example is found bringing my guess down, or the
polygon was not well-formed.)
In 2D the predicate can be implemented exploiting the normals of the
triangle sides. If a point is on the same side with respect to each triangle
side, it is inside the triangle.
Concerning ‘almost’ planar polygons living in 3D: If they behave 'well
enough' to be projected to a common plane, one could use any (or with some
more work a good or even best selection of) three points forming a triangle
as reference plane and project all other points to it. A simple coordinate
system transformation reduces the problem to 2D ...
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16808.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
On 2016-03-30 00:33, Ronaldo wrote:
In 3D space, there is no well defined
orientation for circulating a polygon. In a planar polygon in 3D, the
normal
to the polygon plane may be used to induce an orientation to the vertex
circulation. If the polygon normal is well chosen, it induces the
correct
CCW orientation. If the orientation is not CCW, it is sufficient to
invert
the normal. A simple check of the orientation correctness is a first
issue.
This appears like circular logic, unless you have an explicitly stated
normal independent of the boundary vertices. Usually, the normal is
derived from the boundary vertices. If you have only a general, flat
(possibly concave) polygon in 3d space, you can compute its normal with
two possible outomes, but you cannot know what is "correct" unless you
have something else to compare with.
Almost planar polygonals pose more issues since no normal really
exists.
Actually, there are several normals since the face is curved. In
principle, one could compute some add them vectorially to arrive at some
average.
Carsten Arnholm
Thank you, Sloan, for your comments and suggestions. They gave me more
confidence on may track.
Kenneth Sloan wrote
If all you need is a triangulation, there are two very simple algorithms.
There are others that are much harder to implement, and have minimal
benefit (other than nice asymptotic properties as nVertices becomes very
large.
a) subdivision : at each step, find a “chord” (a connection between two
vertices that lies entirely INSIDE the polygon) and use it to divide the
polygon into two pieces. Recursively triangulate the two pieces.
b) ear-trimming: a special case of a) where only chords cutting off a
single triangle on one side are considered. You look for v0,v1,v2 -
consecutive vertices, where v0->v1->v2 is convex and v0->v2 is entirely
INSIDE the polygon. The naive algorithm is O(n^4), but with care you can
reduce this to O(n^2) - this is an exercise in Joe O’Rourke’s book
Computational Geometry in C
I always advise people to implement b) first - and then (after putting it
into production) investigate “better” algorithms. Don’t bother to push
out a “better” method until customers complain that b) is “too slow”.
That was the way I followed. I implemented a very simple version of the
ear-trimming method. It is very easy to express it in OpenSCAD language
(there is no complex data structure) and I am able to sophisticate it if its
lack of efficiency turns to be a concern. I don't expect it will deal with
more than one thousand points in the applications I have in mind.
Kenneth Sloan wrote
For “slightly non-planar” triangles (and, in general, for triangles
in/near a plane other than the x-y plane, fit a plane to the vertices (I
like eigenvectors for this task) and map all the vertices to that plane
(there are various ways to do this). Triangulate in that plane, and apply
that triangulation to the original points. Define “slightly non-planar”
to be any problem for which this works!
I am working on this line now. For symmetrical 3X3 matrices, there is very
simple and efficient methods using eigenvectors.
Regarding Delaunay triangulation: very very narrow triangles may be a
concern as CGAL may reject them. So some care to find "fat" triangle may be
needed. The only concern I have in implementing triangulation edge flipping
is the required triangulation data structure. The triangulation I generate
with the ear-trimming method is a simple list of the three triangle
vertices. If I devise something more manageable for the edge flipping
procedure, I could get better polygon triangulations.
For a while, I don't intend to consider polygons with holes.
Well before that I need to identify eventual auto-intersections. For now,
auto-intersection will crash my code. I need at least identify it and echo
an error message.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16826.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Parkinbot, thank you for your contribution to the thread.
Parkinbot wrote
Personally I find it impossible to implement some complex algorithm with
OpenSCAD first hand. ...
I understand your first point and agree that with the lack of debugging
tools in OpenSCAD it is hard to develop complex codes. However, I never
worked as a professional programmer and my debugging skills are poor in any
language. Besides, OpenSCAD is my first language based on the principles of
functional programming. So, to learn this style of programming I challenged
myself to code directly in OpenSCAD language. It is a hard track, I know.
The ideas you expressed for the triangulation of simple polygon are the core
of the ear-trimming methods, the exact kind of method I am working with. I
use recursion and an in_triangle test. My method distinguish from your
proposal just that I do the test after the recursion where I am looking for
a new ear (three points in sequence forming a triangle contained in the
polygon). A well known theorem assures that any simple polygon has at least
two ears. So, you are right, if any ear is found, the polygon is not simple.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16827.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
cacb wrote
On 2016-03-30 00:33, Ronaldo wrote:
In 3D space, there is no well defined
orientation for circulating a polygon. In a planar polygon in 3D, the
normal
to the polygon plane may be used to induce an orientation to the vertex
circulation. If the polygon normal is well chosen, it induces the
correct
CCW orientation. If the orientation is not CCW, it is sufficient to
invert
the normal. A simple check of the orientation correctness is a first
issue.
This appears like circular logic, unless you have an explicitly stated
normal independent of the boundary vertices. Usually, the normal is
derived from the boundary vertices. If you have only a general, flat
(possibly concave) polygon in 3d space, you can compute its normal with
two possible outomes, but you cannot know what is "correct" unless you
have something else to compare with.
Almost planar polygonals pose more issues since no normal really
exists.
Actually, there are several normals since the face is curved. In
principle, one could compute some add them vectorially to arrive at some
average.
Let me clarify what I said. Whether or not the polygon vertices are in a
plane, we may project them in a appropriate plane chosen by fitting
techniques, for instance. We may use a plane normal to define a 2D
coordinate system in the fitting plane. With this 2D coordinate system the
sequence of polygon vertices may or may not be the CCW orientation the
triangulation methods require. If we take the symmetrical of the normal, the
induced orientation of the vertices will be reversed. So, a correct choice
of the normal orientation is fundamental. It is easier to invert the normal
than invert the sequence.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16828.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
you are right. It is enough (and much easier) to cycle in search for an ear,
exclude the outer point of it from the sequence and recurse with the
shortend sequence until you are done.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16832.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
On 30. mars 2016 22:02, Ronaldo wrote:
Let me clarify what I said. Whether or not the polygon vertices are in a
plane, we may project them in a appropriate plane chosen by fitting
techniques, for instance. We may use a plane normal to define a 2D
coordinate system in the fitting plane.
Maybe I am not understanding you right, but then you have to choose one
of 2 possible normals of that plane, same problem.
With this 2D coordinate system the
sequence of polygon vertices may or may not be the CCW orientation the
triangulation methods require. If we take the symmetrical of the normal, the
induced orientation of the vertices will be reversed. So, a correct choice
of the normal orientation is fundamental. It is easier to invert the normal
than invert the sequence.
At least in the case where the polygon is planar, you can always compute
the correct normal (regarding CCW vs CW) directly from the vertices,
without having to assume anything about the polygon except that the
vertices are listed sequentially around the contour. The plane normal
defines 3 of 4 coefficients in the plane equation, the 4th can be
trivially computed using one vertex.
I don't know for sure, but I am guessing it would work for "reasonably
planar" polygons as well.
Carsten Arnholm
this works, just the polygon generator does not always produce a non
selfintersecting polygon.
N = 20;
p = randpolyg(N); // too easy to be not perfect
q = triag(p);
for (i=[0:len(q)]) color([i,i,i]/N) polygon(q[i]);
translate([0, 0,-1]) polygon(p);
// generate a randomized polygon
function randpolyg(N=10) = [for (i=[0:360/N:359.9])
let (rnd = rands(3,10,2)) [rnd[0]*sin(i), rnd[1]*cos(i)]];
function triag(S) =
len(S) == 3 ?
[S]: is_valid(S)?
concat([[S[0], S[1], S[2]]],triag(exclude_ear(S))):
triag(cycle(S));
function exclude_ear(S) = [for(i=[0:len(S)-1]) if (i!=1) S[i]];
function cycle(S) = let (n = len(S)) [for (i=[0:n-1]) S[(i+1)%n]];
function is_valid(S) = // false, if not ear or any other point of S inside
ear
is_left(S[0], S[1], S[2])? // ear?
let (res = [for(i=[3:len(S)-1]) // test all other points
if(is_left(S[0], S[1], S[i]) &&
is_left(S[1], S[2], S[i]) &&
is_left(S[2], S[0], S[i])) 1])
res==[]:false;
function is_left(a, b, c) = // true if c is left of a--->b
(b[0] - a[0])(c[1] - a[1]) - (b[1] - a[1])(c[0] - a[0]) < 0;
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16837.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Well, for polyhedron calculation (e.g. for a sweep) one might be more
interested in the indices. Usually it is difficult to keep track of them
within a recursion. An easy solution is to pack the index as extra dimension
with a point:
q1 = triag_I(p);
function triag_I(S) = triag_(arm(S));
function arm(S) = [for (i=[0:len(S)-1]) concat(S[i],[i])];
function triag_(S) = len(S) == 3 ? [I(S)]:
is_valid(S)? concat([I(S)],triag_(exclude_ear(S))): triag_(cycle(S));
function I(S) = [S[0][2], S[1][2], S[2][2]];
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16842.html
Sent from the OpenSCAD mailing list archive at Nabble.com.