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

In reply to this post by Ronaldo
> On Mar 27, 2016, at 18:02 PM, Ronaldo <[hidden email]> wrote: > 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. > We use a modified version of libtess2 (https://github.com/memononen/libtess2) I haven’t managed to find a robust general purpose triangulator which supports slightly nonplanar polygons as well as polygons with holes. Even CGAL crashes on some cornercase polygons. libtess2 is the best I’ve found so far. The “Triangle” library is rumored to be very good, but it’s not open source: https://www.cs.cmu.edu/~quake/triangle.html Marius _______________________________________________ OpenSCAD mailing list [hidden email] http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org 
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

In reply to this post by kintel
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

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. 20160328 15:29 GMT03:00 Torsten Paul <[hidden email]>: On 03/28/2016 07:34 PM, Ronaldo wrote: _______________________________________________ OpenSCAD mailing list [hidden email] http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org 
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

In reply to this post by Parkinbot
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 
In reply to this post by Ronaldo
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 ... 
In reply to this post by Ronaldo
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 
In reply to this post by Kenneth Sloan
Thank you, Sloan, for your comments and suggestions. They gave me more confidence on may track.
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. 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. 
In reply to this post by Parkinbot
Parkinbot, thank you for your contribution to the thread.
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. 
In reply to this post by cacb
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. 
In reply to this post by Ronaldo
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.

In reply to this post by Ronaldo
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 
In reply to this post by Parkinbot
this works, just the polygon generator does not always produce a non selfintersecting polygon.

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:

Free forum by Nabble  Edit this page 