
1234
... 6

I am searching for a simple polygon triangulation algorithm to implement in OpenSCAD language. Yes, I know that OpenSCAD is not a general purpose language. But all methods I found for this problem are essentially iterative and require complex data structure hard to implement even in a general purpose functional language. The robust methods are complex to implement. Even the Seidel algorithm, which is allegedly the easier to code, relies on a data structure that required more that 120 lines of code of a C++ implementation I found.
Besides, all methods I found are restricted to polygons in the plane. Which method is used to triangulate the faces of polyhedron in OpenSCAD preview? It seems to be more general.
Any help is welcome.


good question ;)
so you'd need a convexity test first.
1. 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.
2. For a nonconvex polygon, you can decompose the shape into a set convex shapes and use 1a. So your question could be, how to decompose a nonconvex polygon into a set of convex polygons.
3. I think alternatively you could also try 1b and find some additional test for doing the right choice.


Maybe contacting the author, he could change the licence, or give a special permission, considering OpenSCAD is non commercial.
jpmendes


On 03/28/2016 05:19 PM, jpmendes wrote:
> Maybe contacting the author, he could change the licence, or
> give a special permission, considering OpenSCAD is non commercial.
>
The first would be nice, the second one can't work. A special
license for OpenSCAD would make it nondistributable.
ciao,
Torsten.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
 Torsten


Hi, Kintel. Thank you for the references. I will take a look on libtess2.
I visited the Triangle page you mentioned and downloaded the triangle.c code. It is intended to perform Delaunay triangulations. Relying on its ability to compute constrained Delaunay triangulation may be it could be used to triangulate a simple planar polygon. It could be an alternative to me if I had triangulation data structures implemented over OpenSCAD lists. That was in my future plans but this might be a huge work.
I already implemented an eartrim kind of method to triangulate simple almost planar polygons. It is not very efficient however. And very thin triangles is unavoidable.
Regarding the Triangle copyright, the text bellow is found in the triangle.c file. I am not sure if it would restrict its use in OpenSCAD code.
/* (triangle.c) */
/* */
/* Version 1.6 */
/* July 28, 2005 */
/* */
/* Copyright 1993, 1995, 1997, 1998, 2002, 2005 */
/* Jonathan Richard Shewchuk */
/* 2360 Woolsey #H */
/* Berkeley, California 947051927 */
/* jrs@cs.berkeley.edu */
/* */
/* This program may be freely redistributed under the condition that the */
/* copyright notices (including this entire header and the copyright */
/* notice printed when the `h' switch is selected) are not removed, and */
/* no compensation is received. Private, research, and institutional */
/* use is free. You may distribute modified versions of this code UNDER */
/* THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE */
/* SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE */
/* AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR */
/* NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution of this code as */
/* part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT */
/* WITH THE AUTHOR. (If you are not directly supplying this code to a */
/* customer, and you are instead telling them how they can obtain it for */
/* free, then you are not required to make any arrangement with me.) */


On 03/28/2016 07:34 PM, Ronaldo wrote:
> Regarding the Triangle copyright, the text bellow is found in
> the triangle.c file. I am not sure if it would restrict its use
> in OpenSCAD code.
>
I'm not a lawyer, so I can't give official advice about this,
but in my understanding "Private, research, and institutional
use is free." is a big conflict with GPL.
ciao,
Torsten.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
 Torsten


On 03/29/2016 01:24 AM, Ronaldo Persiano wrote:
> I agree but between "Private, research, and institutional use is free." and
> "Distribution of this code as part of a commercial system is permissible
> ONLY BY DIRECT ARRANGEMENT WITH THE AUTHOR" there is a large space for a deal.
>
Yes, there might be room for a deal, but that will not help as long
as there is any additional GPL code used. *Any* special agreement
just for OpenSCAD is not compatible with GPL.
The only way to use that code is if it's released with a license
that is compatible with GPL (e.g. BSD/MIT is compatible, a long
list is available at http://www.gnu.org/licenses/licenselist.html).
Note that this is not so much about the OpenSCAD code itself, this
could be mainly up to Marius to decide, but we use a number of GPL
libraries and we do have to respect their decision.
ciao,
Torsten.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
 Torsten


Parkinbot wrote
good question ;)
so you'd need a convexity test first.
1. 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.
2. For a nonconvex polygon, you can decompose the shape into a set convex shapes and use 1a. So your question could be, how to decompose a nonconvex polygon into a set of convex polygons.
3. 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 earcut methods seem to be easier to implement in OpenSCAD. The simplest ones are suboptimal 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 earcut 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.


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) eartrimming: 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 nonplanar” triangles (and, in general, for triangles in/near a plane other than the xy 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 nonplanar” 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 nonplanar” 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 “Delaunaylike” 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 decompose 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 earclipping can be made to understand holes  I just don’t know, and (again) that might depend on the representation.

Kenneth Sloan
[hidden email]
Vision is the art of seeing what is invisible to others.
> On Mar 29, 2016, at 17:33 , Ronaldo < [hidden email]> wrote:
>
> Parkinbot wrote
>> good question ;)
>>
>> so you'd need a convexity test first.
>> 1. 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.
>>
>> 2. For a nonconvex polygon, you can decompose the shape into a set convex
>> shapes and use 1a. So your question could be, how to decompose a
>> nonconvex polygon into a set of convex polygons.
>>
>> 3. 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 earcut methods seem to be easier to implement in OpenSCAD. The simplest
> ones are suboptimal 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 earcut 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/Simplepolygontriangulationtp16755p16803.html> Sent from the OpenSCAD mailing list archive at Nabble.com.
>
> _______________________________________________
> OpenSCAD mailing list
> [hidden email]
> http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org_______________________________________________
OpenSCAD mailing list
[hidden email]
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 nonconvex 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 wellformed (=non selfintersecting, non singular) Npolygon 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 wellformed Npolygon, 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 Nsequence, 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 N4 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 wellformed.)
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 ...


On 20160330 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
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


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) eartrimming: 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 eartrimming 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 nonplanar” triangles (and, in general, for triangles in/near a plane other than the xy 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 nonplanar” 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 eartrimming 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 autointersections. For now, autointersection will crash my code. I need at least identify it and echo an error message.


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 eartrimming 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.


cacb wrote
On 20160330 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.


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.


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
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


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:n1]) 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;


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]];

1234
... 6
