Simple polygon triangulation
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.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
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.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16759.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
On Mar 27, 2016, at 18:02 PM, Ronaldo rcmpersiano@gmail.com 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 non-planar polygons as well as polygons with holes. Even CGAL crashes on some corner-case 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
Maybe contacting the author, he could change the licence, or give a special
permission, considering OpenSCAD is non commercial.
jpmendes
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16764.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
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 non-distributable.
ciao,
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 ear-trim 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 94705-1927
/
/ 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.)
*/
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16771.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
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.
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.
2016-03-28 15:29 GMT-03:00 Torsten Paul Torsten.Paul@gmx.de:
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
Discuss@lists.openscad.org
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/license-list.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.
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.
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.
I see your point. So what you with a almost planar 2D polygon P living
anywhere in 3D is:
- choose any three non-colinear points of P and construct a plane E
- find a map T() that maps the base vectors of E to the base {x, y}
- calculate P'=T(P) and drop the z-dimension (projection)
- cycle around P' and count for each triple of consequent points whether
the third point is left with respect to the first two. If this count is
lower than len(P)/2 use the orientation CW otherwise CCW.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16844.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
On 2016-03-31 13:11, Parkinbot wrote:
I see your point. So what you with a almost planar 2D polygon P living
anywhere in 3D is:
- choose any three non-colinear points of P and construct a plane E
This will not always work. It will give you a plane, but you don't know
which of the 2 possible normals that plane will have, since you cannot
assume the polygon is convex. Look up "Newell's method". It will give
you the polygon normal corresponding to CCW directly for a planar
polygon, including concave polygons. As mentioned, a slightly non-planar
might work too, but I have not checked.
Carsten Arnholm
Your point is adressed in step 4 not in step 1.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16850.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Parkinbot wrote
this works, just the polygon generator does not always produce a non
selfintersecting polygon.
Surprisingly short solution, Parkinbot. However, in my tests, it does not
deal with more than a few dozen vertices. You recurring cycling is the
villain. And to avoid it you should consider that not only the first three
vertices in the sequence is a candidate to be an ear. The following code
changes yours addressing those issues. I provide also a detection of
non-simple polygons to avoid crashing. It is a lot slow yet, taking 16sec
for 200 vertices excluding the display.
seed = floor(rands(0,10000,1)[0]);
//seed = 4628;
echo(seed);
N = 400;
p = randpolyg(N); // too easy to be not perfect
q = triag(p);
if( len(q[len(q)-1]) > 3) echo("Non-simple polygon");
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) = let (rnd = rands(3,10,2N,seed))
[for (i=[0:360/N:359.9]) [rnd[2Ni/360]sin(i),
rnd[2Ni/360+1]*cos(i)]];
function triag(S) =
len(S) <= 3 ?
[S]:
let( ear = valid_ears(S)[0], l = len(S) )
ear== undef? [] : // case S is a non-simple polygon
concat([[S[(ear+l-1)%l], S[ear], S[(ear+1)%l]]],
triag(exclude_vertex(ear,S), k) );
function exclude_vertex(k,S) =
concat( k<=0? []: [for(i=[0 :k-1]) S[i] ] ,
k>=len(S)-1? []: [for(i=[k+1:len(S)-1]) S[i] ] );
function valid_ears(S) =
[ for( i=[0:len(S)-1]) if( is_valid(i,S) ) i ];
function is_valid(k,S) = // false, if not ear or any other point of S
inside ear
let( l = len(S) )
is_left(S[(k+l-1)%l], S[k], S[(k+1)%l],0)? // is S[k] an ear?
let ( res = [ for(i=[0:l-1]) // test all other points
if( (i!=k) && (i!=(k+l-1)%l) && (i!=(k+1)%l) &&
is_left(S[(k+l-1)%l], S[k], S[i],1) &&
is_left(S[k], S[(k+1)%l], S[i],1) &&
is_left(S[(k+1)%l], S[(k+l-1)%l], S[i],1)) 1]
)
res==[]:
false;
function is_left(a, b, c, ind) = // true if c is left of a--->b
ind==0?
(b[0] - a[0])(c[1] - a[1]) - (b[1] - a[1])(c[0] - a[0]) < 0 :
(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-tp16755p16852.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
As I wrote: it is the random polygon generator that freaks out and delivers
malformed polygons, not the triangulation.
As the code cycles around the part of the polygon not being triangulated
yet, it is indeed enough to just regard the first three points.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16853.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Parkinbot wrote
I see your point. So what you do with an almost planar 2D polygon P living
anywhere in 3D is:
- choose any three non-colinear points of P and construct a plane E
- find a map T() that maps the base vectors of E to the base {x, y}
- calculate P'=T(P) and drop the z-dimension (projection)
- cycle around P' and count for each triple of consequent points whether
the third point is left with respect to the first two. If this count is
lower than len(P)/2 use the orientation CW otherwise CCW.
I am trying another approach to step 4. Find the indices of the extreme
vertices of P': those with minimum x, minimum y, maximum x, maximum y, in
that order. If they are cycling ordered (that is they follows the order of
the CCW orientation) you got it. As an example, if the indices are [ 4, 0,
0, 3] it is CCW, if they are [3,0,0,4] it is CW.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16854.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Parkinbot wrote
- choose any three non-colinear points of P and construct a plane E
I would rather prefer a plane fitting. Plane E is a good choice for that
three points but may be very bad for the remaining.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16855.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
I agree. There are lots of ad hoc ways of almost solving this problem - eventually everyone figures out that you need to do a proper fit of a plane to all the points. If you compute eigenvectors for the covariance matrix, you can even use the eigenvalues to determine whether or not the points lie “almost” in a plane. And then, you might as well use the eigenvector associated with the largest eigenvalue as your x-axis.
--
Kenneth Sloan
KennethRSloan@gmail.com
Vision is the art of seeing what is invisible to others.
On Mar 31, 2016, at 13:03 , Ronaldo rcmpersiano@gmail.com wrote:
Parkinbot wrote
- choose any three non-colinear points of P and construct a plane E
I would rather prefer a plane fitting. Plane E is a good choice for that
three points but may be very bad for the remaining.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16855.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
Shrewd!
2016-03-31 15:09 GMT-03:00 Kenneth Sloan kennethrsloan@gmail.com:
I agree. There are lots of ad hoc ways of almost solving this problem -
eventually everyone figures out that you need to do a proper fit of a plane
to all the points. If you compute eigenvectors for the covariance matrix,
you can even use the eigenvalues to determine whether or not the points lie
“almost” in a plane. And then, you might as well use the eigenvector
associated with the largest eigenvalue as your x-axis.
--
Kenneth Sloan
KennethRSloan@gmail.com
Vision is the art of seeing what is invisible to others.
On Mar 31, 2016, at 13:03 , Ronaldo rcmpersiano@gmail.com wrote:
Parkinbot wrote
- choose any three non-colinear points of P and construct a plane E
I would rather prefer a plane fitting. Plane E is a good choice for that
three points but may be very bad for the remaining.
--
View this message in context:
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
Parkinbot wrote
As I wrote: it is the random polygon generator that freaks out and
delivers malformed polygons, not the triangulation.
As the code cycles around the part of the polygon not being triangulated
yet, it is indeed enough to just regard the first three points.
I agree. But my objection is exactly the cycling that increases the
inefficiency and, I suppose, quickly grows the recursion stack.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16859.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Meanwhile I've tested a polygon with 900 points and it took 10s. So, we are
talking about efficiency now.
Letting OpenSCADs poor data manipulation ability and unnessecary copies when
cycling aside, it is the algorithm that takes O(n²) steps, because the ear
cut method relies on testing all the rest of the points in the polygon with
each new candidate.
You can cut that down by implementing the second simple scheme Kenneth
mentioned, because recursively cutting a polygon into two, reduces the
complexity by splitting the set of points to be tested by a division factor
and not only by 1. The method is similar to Quicksort.
You start to split with a guess:
- put floor(n/2) points into a first set, and the rest into a second set.
(Alternatively a random split can be used) - test if all points (except the first one) in the first set are left with
respect to the vector that connects the first points of the sets. - if the test is true, recurse both sets separately until done. If not,
cycle your split (or do alternatively another random guess) and continue
with 2.
There are other schemes mentioned in
https://en.wikipedia.org/wiki/Polygon_triangulation#Ear_clipping_method
which are even faster, but wasn't it you who wanted an implementation with a
small footage?
- Rudolf -
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16863.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
I suspect your algorithm may be O(n^4) and not O(n^2). Testing should resolve this.
Getting it to O(n^3) is easy; getting to O(n^2) is a bit harder.
As long as you are doing recursion, instead of iteration, stepping up to the general subdivision by a chord is probably worthwhile. If you do that, try to bias the selection of the chord to divide the polygon into two more-or-less equal pieces.
Kenneth Sloan
KennethRSloan@gmail.com
Vision is the art of seeing what is invisible to others.
On Mar 31, 2016, at 14:48 , Parkinbot rudolf@parkinbot.com wrote:
Meanwhile I've tested a polygon with 900 points and it took 10s. So, we are
talking about efficiency now.
Letting OpenSCADs poor data manipulation ability and unnessecary copies when
cycling aside, it is the algorithm that takes O(n²) steps, because the ear
cut method relies on testing all the rest of the points in the polygon with
each new candidate.
You can cut that down by implementing the second simple scheme Kenneth
mentioned, because recursively cutting a polygon into two, reduces the
complexity by splitting the set of points to be tested by a division factor
and not only by 1. The method is similar to Quicksort.
You start to split with a guess:
- put floor(n/2) points into a first set, and the rest into a second set.
(Alternatively a random split can be used) - test if all points (except the first one) in the first set are left with
respect to the vector that connects the first points of the sets. - if the test is true, recurse both sets separately until done. If not,
cycle your split (or do alternatively another random guess) and continue
with 2.
There are other schemes mentioned in
https://en.wikipedia.org/wiki/Polygon_triangulation#Ear_clipping_method
which are even faster, but wasn't it you who wanted an implementation with a
small footage?
- Rudolf -
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16863.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
Kenneth Sloan wrote
I suspect your algorithm may be O(n^4) and not O(n^2). Testing should
resolve this.
You are right, it might be even worse ;-). But OpenSCAD is not well enough
equipped for implementing iterative solutions. You have random read access
to vector elements, but a random write forces you to do a full copy of the
whole structure. You can't break a for loop, have no while loop, can't even
sum up anything within a for loop - the only way: do it recusively. The full
curse of functional programming :-)
I might give it a try for the second approach tomorrow, but I wouldn't spend
any minute to tweek an algorithm written in OpenSCAD for speed.
- Rudolf -
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16865.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
I am not expert on this but I estimated that the method we are discussing is
O(n^3). The valid_ears() search is O(n^2), due to two nested for loops. The
recursive triangulation has depth n and call valid_ears() once for each
found triangle. I wouldn't blame OpenSCAD language for the low efficiency;
an implementation of this algorithm in C will be faster but it will have
the same asymptotic complexity.
I have tested it with polygons up to 1600 vertices ( a 313 sec run without
display) and it behaved as O(n^2) which is poor for so small cases. The test
polygon I have used has only four ears even in the intermediate steps which
seems to be the worst case for it. This is the starting CCW polygon:
N = 1600;
M = 5;
function test_polygon(N,M,c,r,a) =
concat(
[ for(i=[ a:-2a/N:-a]) [ c - 1.5rcos(i), rsin(i) ] ],
[ for(i=[ -a: 2a/M: a]) [ -c + 1.5rcos(i), rsin(i) ] ] );
p = test_polygon(N,M, 30, 15, 60);
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16869.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Parkinbot wrote:
OpenSCAD is not well enough equipped for implementing iterative solutions.
You have random read access
to vector elements, but a random write forces you to do a full copy of the
whole structure. You can't break a for loop, have no while loop, can't even
sum up anything within a for loop - the only way: do it recusively. The
full
curse of functional programming :-)
Thanks for pointing this out. General purpose functional languages have a
better set of control structures, and can efficiently update a single
element of an array. It would be reasonable to add array update to
OpenSCAD. We just need a good syntax. I'll think about this.
How do you update an array element without mutable data?
On 1 April 2016 at 11:19, doug moen doug@moens.org wrote:
Parkinbot wrote:
OpenSCAD is not well enough equipped for implementing iterative
solutions. You have random read access
to vector elements, but a random write forces you to do a full copy of the
whole structure. You can't break a for loop, have no while loop, can't
even
sum up anything within a for loop - the only way: do it recusively. The
full
curse of functional programming :-)
Thanks for pointing this out. General purpose functional languages have a
better set of control structures, and can efficiently update a single
element of an array. It would be reasonable to add array update to
OpenSCAD. We just need a good syntax. I'll think about this.
OpenSCAD mailing list
Discuss@lists.openscad.org
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Abstractly, we could provide a function upd(a,i,x) where a is an array, i
is an index, and x is the new value to store in a[i]. The result is a copy
of a in which a[i] is set to x.
There are various ways to implement this efficiently. There is a field of
study called "functional data structures" which provides solutions. In the
OpenSCAD implementation, values are reference counted. A simple approach is
to check if a has a reference count of 1: If so, we destructively update a,
and return it. If not, we return a copy of a with element i changed to x.
That would be good enough.
On 1 April 2016 at 06:31, nop head nop.head@gmail.com wrote:
How do you update an array element without mutable data?
On 1 April 2016 at 11:19, doug moen doug@moens.org wrote:
Parkinbot wrote:
OpenSCAD is not well enough equipped for implementing iterative
solutions. You have random read access
to vector elements, but a random write forces you to do a full copy of
the
whole structure. You can't break a for loop, have no while loop, can't
even
sum up anything within a for loop - the only way: do it recusively. The
full
curse of functional programming :-)
Thanks for pointing this out. General purpose functional languages have a
better set of control structures, and can efficiently update a single
element of an array. It would be reasonable to add array update to
OpenSCAD. We just need a good syntax. I'll think about this.
OpenSCAD mailing list
Discuss@lists.openscad.org
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
I see, so fine to mutate variables as long as no other references exist.
Unless it is a literal or a temporary expression result, won't it always be
bound to at least a local variable?
When calling a function and passing a variable as an argument you would
need look ahead and see that the variable was not used again and then
decrement the reference count artificially early.
On 1 April 2016 at 11:52, doug moen doug@moens.org wrote:
Abstractly, we could provide a function upd(a,i,x) where a is an array, i
is an index, and x is the new value to store in a[i]. The result is a copy
of a in which a[i] is set to x.
There are various ways to implement this efficiently. There is a field of
study called "functional data structures" which provides solutions. In the
OpenSCAD implementation, values are reference counted. A simple approach is
to check if a has a reference count of 1: If so, we destructively update a,
and return it. If not, we return a copy of a with element i changed to x.
That would be good enough.
On 1 April 2016 at 06:31, nop head nop.head@gmail.com wrote:
How do you update an array element without mutable data?
On 1 April 2016 at 11:19, doug moen doug@moens.org wrote:
Parkinbot wrote:
OpenSCAD is not well enough equipped for implementing iterative
solutions. You have random read access
to vector elements, but a random write forces you to do a full copy of
the
whole structure. You can't break a for loop, have no while loop, can't
even
sum up anything within a for loop - the only way: do it recusively. The
full
curse of functional programming :-)
Thanks for pointing this out. General purpose functional languages have
a better set of control structures, and can efficiently update a single
element of an array. It would be reasonable to add array update to
OpenSCAD. We just need a good syntax. I'll think about this.
OpenSCAD mailing list
Discuss@lists.openscad.org
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
@nop.head: Yes, those are good points. Testing if a variable isn't
referenced again is a standard operation in an optimizing compiler. It
could also be used to optimize away unnecessary reference count
manipulations in function calls. If you are calling f(x), and there are no
further references to x, then you could eliminate an increment of x's
reference count when f is called, and a decrement of x's reference count
after f returns. These optimizations would be easier to implement if we had
a compiler.
On 1 April 2016 at 07:15, nop head nop.head@gmail.com wrote:
I see, so fine to mutate variables as long as no other references exist.
Unless it is a literal or a temporary expression result, won't it always
be bound to at least a local variable?
When calling a function and passing a variable as an argument you would
need look ahead and see that the variable was not used again and then
decrement the reference count artificially early.
On 1 April 2016 at 11:52, doug moen doug@moens.org wrote:
Abstractly, we could provide a function upd(a,i,x) where a is an array, i
is an index, and x is the new value to store in a[i]. The result is a copy
of a in which a[i] is set to x.
There are various ways to implement this efficiently. There is a field of
study called "functional data structures" which provides solutions. In the
OpenSCAD implementation, values are reference counted. A simple approach is
to check if a has a reference count of 1: If so, we destructively update a,
and return it. If not, we return a copy of a with element i changed to x.
That would be good enough.
On 1 April 2016 at 06:31, nop head nop.head@gmail.com wrote:
How do you update an array element without mutable data?
On 1 April 2016 at 11:19, doug moen doug@moens.org wrote:
Parkinbot wrote:
OpenSCAD is not well enough equipped for implementing iterative
solutions. You have random read access
to vector elements, but a random write forces you to do a full copy of
the
whole structure. You can't break a for loop, have no while loop, can't
even
sum up anything within a for loop - the only way: do it recusively.
The full
curse of functional programming :-)
Thanks for pointing this out. General purpose functional languages have
a better set of control structures, and can efficiently update a single
element of an array. It would be reasonable to add array update to
OpenSCAD. We just need a good syntax. I'll think about this.
OpenSCAD mailing list
Discuss@lists.openscad.org
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org
Would it be such an offense to the functional paradigma of OpenSCAD, which I
guess makes a lot of sense for CSG, but not at all for number crunching, to
allow value altering for local variables within function scope?
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16889.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
My vote is in favor. OpenSCAD is a bit more than geometry and topology. In
animation for instance, to produce in an easy way good animated sequences we
need to make state machines and record events.
I know this is not the object of the language, but why not extend it to
facilitate and stimulate creativity. Most of users are not MScs or Phds
searching for the optimal way to produce great geometrical creations. Of
course this aspect is important as it is related with the main goal of the
language. However the popularity of the language in which users have a
fundamental role would be augmented if some new language tools could
facilitate their expression. As an example, some contributors offer great
libraries that most of us can't use because the most common of us don't
understand it and in it's majority, documentation is scarse. OpenSCAD is
becoming more popular and one of it's strengths is the easy way one can
express some creativity at its own level of knowledge.
jpmendes
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16892.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
You don't need state to produce animation sequences. I once designed an
animation system that had to synchronise across lots of adjacent machines
that were connected to their left and right neighbours. From that they
could decide where they were and in the row, how long the row was and the
most left one passed a time variable up stream. From that any arbitrary
animation sequence were run. You simply have to divide down the time
variable to work out what state you are in if your animation has states.
If you want to write in a procedural language why not use the Python
binding, or the C++ binding, or Angel script, etc.
On 1 April 2016 at 20:10, jpmendes jpmendes54@gmail.com wrote:
My vote is in favor. OpenSCAD is a bit more than geometry and topology. In
animation for instance, to produce in an easy way good animated sequences
we
need to make state machines and record events.
I know this is not the object of the language, but why not extend it to
facilitate and stimulate creativity. Most of users are not MScs or Phds
searching for the optimal way to produce great geometrical creations. Of
course this aspect is important as it is related with the main goal of the
language. However the popularity of the language in which users have a
fundamental role would be augmented if some new language tools could
facilitate their expression. As an example, some contributors offer great
libraries that most of us can't use because the most common of us don't
understand it and in it's majority, documentation is scarse. OpenSCAD is
becoming more popular and one of it's strengths is the easy way one can
express some creativity at its own level of knowledge.
jpmendes
--
View this message in context:
http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16892.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
Although strictly not necessary it helps a lot. Specially if we want to keep
track of object mutual interference. Maybe a mathematician can describe a
path for objects interfering with each other and predict the time for the
occurrences, but it is much simpler to record those events. As this is an
off topic let's stop for now. Thanks for the suggestion I will think more
about it.
jpmendes
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16897.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
I think I have a full functional implementation of a O(n^3) ear-cutting
method to triangulate 2D simple polygon. I addopted the data structure
Parkinbot suggested to output the triangles as lists of the indices in the
array of vertices. The code is relatively simple:
// ************* 2D triangulation of a polygon ********************
// triangulate the 2D polygon lp
function triangulate(lp) = triangs(build_data(lp));
function triangs(lsp) = // traiangulate structured 2D polygon
len(lsp) <= 3 ?
[ [ind(lsp[0]), ind(lsp[1]), ind(lsp[2])] ]:
let( ear = any_ear(lsp), l = len(lsp) )
ear== undef? [] : // case S is a non-simple polygon
concat( [ [ ind(lsp[(ear+l-1)%l]), ind(lsp[ear]), ind(lsp[(ear+1)%l])
] ],
triangs(exclude_vertex(ear, lsp)) );
// small data structure to preserve original indexing of points
function build_data(lp) = [ for(i=[0:len(lp)-1]) [ i, lp[i] ] ];
function ind(sp) = sp[0];
function vtx(sp) = sp[1];
function exclude_vertex(k,lsp) =
concat( k<=0? []: [for(i=[0 :k-1]) lsp[i] ] ,
k>=len(lsp)-1? []: [for(i=[k+1:len(lsp)-1]) lsp[i] ] );
// find the index of the best ear according to inner angle criterium
function any_ear(lsp) =
let( l = len(lsp) )
[ for(j=[0:l-1]) if( is_ear(j, lsp) ) j ] [0];
// false if vtx(lsp[k1]) not convex or
// any other point of lsp is inside the triangle
function is_ear(k1,lsp) =
let( l = len(lsp), k0 = (k1+l-1)%l, k2 = (k1+1)%l )
!is_left(vtx(lsp[k1]), vtx(lsp[k0]), vtx(lsp[k2])) // vtx(lsp[k1]) is
convex
&& void( [ for(i=[0:l-1]) // and other vertices are outside
if( (i!=k0) && (i!=k1) && (i!=k2) &&
is_left(vtx(lsp[k0]), vtx(lsp[k1]), vtx(lsp[i])) &&
is_left(vtx(lsp[k1]), vtx(lsp[k2]), vtx(lsp[i])) &&
is_left(vtx(lsp[k2]), vtx(lsp[k0]), vtx(lsp[i]))) 1 ]
) ;
// true if vtx(sp3) is at left or is colinear to the vector
vtx(sp2)-vtx(sp1)
function is_left(a, b, c) =
(b[0] - a[0])(c[1] - a[1]) - (b[1] - a[1])(c[0] - a[0]) <= 0 ;
function void(lst) = len(lst)==0 ;
For almost planar closed polygonals in the 3D space, I do a plane fitting to
the vertices, project them into the plane and define a local 2D coordinate
system before the triangulation. This process is entirely independent of the
triangulation above:
// ******* plane fitting and projection ***************
function project_poly(lp3D) =
let( leigen = least_square_transform_3(lp3D),
pm = leigen[0], // mean point
eigenv = leigen[1], // eigenvectors
projection = [ eigenv[0], eigenv[1] ], // projection transform
lproj2D = [ for(i=[0:len(lp3D)-1]) projection*(lp3D[i]-pm) ] )
is_CCW(lproj2D) ?
lproj2D*[ [-1,0],[0,1] ] :
lproj2D;
// check whether the orientation induced by projection is CCW
function is_CCW(lp) =
let( bd = bounds(lp),
im = index_min(bd) ) // the index of the smaller index in bd
( bd[im] <= bd[(im+1)%4] )
&& ( bd[(im+1)%4] <= bd[(im+2)%4] )
&& ( bd[(im+2)%4] <= bd[(im+3)%4] ) ;
// indices of extreme points of p
function bounds(lp) =
let( indminX = index_min( lp*[ 1, 0]), indmaxX = index_min( lp*[-1,
0]),
indminY = index_min( lp*[ 0, 1]), indmaxY = index_min( lp*[
0,-1]) )
[ indminX, indminY, indmaxX, indmaxY ]; // this order is essential
// the index of the smaller value in the list l
function index_min(l) = let( minv = min(l) )
[ for(i=[0:len(l)-1]) if ( l[i] == minv ) i ][0];
As you can see, the CCW orientation of the vertices is ensured by a well
chosen projection. The data fitting relies on the
least_square_transform_3(). That function returns the three pairs of
eigenvalue/eigenvector of the covariance matrix of the points. As those
functions may be used by other applications, I kept them in a file with
other linear algebra stuffs. This file may be retrived bellow:
linear_algebra.scad
http://forum.openscad.org/file/n16898/linear_algebra.scad
To improve the time complexity of the method for triangulation of the simple
polygon, we need to find an ear in constant time. There is solutions to that
but it requires more sophisticated data structure. I will not do it now, may
be later. Now I will return to the application that demanded the
triangulation.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16898.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Ronaldo, if you are interested, this is my implementation of the subdivision
algorithm.
This time I combined cycle and split to avoid an extra copy.
I works as expected, but not on all generated polygons. I suspect, now is
not the generator, but the algorithm freaking out on rare occasions. Don't
have the time, to find out the reason, without a debugging facility. I
testet the predicates and functions and gave some test code for them.
Maybe four eyes see more than two ...
N = 19;
p = randpolyg(N); // too easy to be perfect
q = triag_div(p);
for (i=[0:len(q)-1]) color([i,i,i]/(N)) polygon(q[i]); translate([0,
0,-1]) polygon(p);
// testing of the subfunctions
echo(is_left([0,0], [1,1], [2,3])); // true
echo(is_left([0,0], [1,1], [3,1])); // false
echo(is_in([[0,0], [1,1], [4,3]], [2,2])); // false
echo(is_in([[0,0], [1,1], [4,3]], [1.1,1])); // true
echo(split([1,2,3,4])); // [[1, 2, 3], [3, 4, 1]]
echo(split([1,2,3,4], 1)); //+cycle [[2, 3, 4], [4, 1, 2]]
echo(split([1,2,3,4,5])); // [[1, 2, 3], [3, 4, 5, 1]]
echo(split([1,2,3,4,5, 6])); // [[1, 2, 3, 4], [4, 5, 6, 1]]
// generate a randomized CW polygon
function randpolyg(N=10) = [for (i=[0:360/N:359.9])
let (rnd = rands(5,10,2)) [rnd[0]*sin(i), rnd[1]*cos(i)]];
// the main algorithm
function triag_div(S, cycle=0) =
len(S)==3? [S]:
let(A = split(S, cycle))
is_division(A) ?
concat(
triag_div(A[0]),
triag_div(A[1]))
:triag_div(S, cycle+1);
// split a polygon in the middle after cycling
function split(S, cycle=0) = let (n=len(S), cut=floor(n/2))
[[for(i=[0:cut]) S[(i+cycle)%n]], [for(i=[cut:n]) S[(i+cycle)%n]]];
// true, if all points of first set outside of second set
// and all points of second set outside of first set
function is_division(A) =
let (res1 = [for(i=[1:len(A[0])-2]) if(is_in(A[1], A[0][i])) 1])
let (res2 = [for(i=[1:len(A[1])-2]) if(is_in(A[0], A[1][i])) 1])
res1==[] && res2==[];
// true if p is in A
function is_in(A, p) = let(n=len(A)) let(res =
[for(i=[0:n-1]) if(is_left(A[i], A[(i+1)%n], p)) 1])
res==[];
// true if c is left of a--->b
function is_left(a, b, c) =
(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-tp16755p16900.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Thank you, Parkinbot. Your code looks nice. I will comment it tomorrow.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16909.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Parkinbot, please try your subdivision algorithm with this polygon
N = 19;
pre_seed = 5528; //-1; // if pre_seed <0 a random seed is generated
seed = pre_seed<0 ? floor(rands(0,10000,1)[0]): pre_seed;
echo(seed=seed);
// generate a randomized CW polygon
function randpolyg(N=10) = let( rnd = rands(5,10,2N+1,seed) )
[for (i=[0:360/N:359.9]) [rnd[2iN/360]sin(i), rnd[2iN/360+1]*cos(i)]];
p = randpolyg(N); // too easy to be perfect
2016-04-02 15:27 GMT-03:00 Ronaldo rcmpersiano@gmail.com:
Thank you, Parkinbot. Your code looks nice. I will comment it tomorrow.
--
View this message in context:
http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16909.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
This is not my area, but maybe you guys interested in tessellation can take
a look at "Wings 3D" sources. Wings 3D is written in Erlang also a
functional language, perhaps some algorithms maybe easily adapted to
OpenSCAD by some experienced developer. I downloaded the sources just for
curiosity and the code seems to be well structured. A specific module on
tessellation exists with some interesting comments. BTW Synwrite has a
lexter for Erlang.
jpmendes
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16917.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Meanwhile I know, what is problem. The code's logic is fine, but the
is_division() predicate is too weak. It can fail on certain non-convex
partitions. I think I'll have to give it new try.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16918.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
I have noticed that. Besides, your main algoritm assumes that, for any
simple polygon, there is a diagonal from some p[i] to p[i+floor(n/2)] for
some i. This is not true in general.
2016-04-02 20:05 GMT-03:00 Parkinbot rudolf@parkinbot.com:
Meanwhile I know, what is problem. The code's logic is fine, but the
is_division() predicate is too weak. It can fail on certain non-convex
partitions. I think I'll have to give it new try.
--
View this message in context:
http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16918.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
You are right. I planned to include a randomized cutting scheme, but as the
algorithm was not traceble for me anymore, I thought it was a bad idea to
put another randomness on top of it.
So my new idea is to check the subdivision line AB against all other lines
PQ built by two subsequent points (with P,Q<>A,B of course). If they
intersect, the polygon APBQ will be convex in a CW or CCW sense. This is
easy to test with is_left()
Being short of time in the next days, I might implement it in Matlab next
weekend. With a debugger one can effectively trace recursions.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16920.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
Ronaldo wrote
I think I have a full functional implementation of a O(n^3) ear-cutting
method to triangulate 2D simple polygon. I addopted the data structure
Parkinbot suggested to output the triangles as lists of the indices in the
array of vertices. The code is relatively simple:
yeah, nice! Much more interative.
Two annotations:
-
As you collect the set of all ears in one cycle, isn't it that you test
the ears already within the set again and again? Avoiding this might bring
you closer to O(n²) :-) -
I would suspect that approximating the best plane is not much safer than
using a suboptimal plane, as long as you don't have a criteria assuring that
the projection to this best (or any other) plane will be homomorph. Indeed
you easily can construct polygons, that don't project homomorph to the best
plane, but to a suboptimal plane. Also Eigenvector calculation takes time.
Projecting the polygon to xy, yz, zx planes and ranking the point triples
defining the largest areas within those projections in 3D, might be a good
and fast heuristics, until you give this criteria.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16923.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
nophead wrote
If you want to write in a procedural language why not use the Python
binding, or the C++ binding, or Angel script, etc.
If you want us all to write in a functional language, why not offer at least
an echo() function primitive which outputs a value to console while
returning it unchanged.
function myfunction(a) = let(x=echo(f(a,i))) // inspect this value
g(x); // calculation some other value from it
function f(a,i) = i;// ... some function
function g(x) = x;// ... some other function
function echo(x) = x; // side effect to console
Also a break-function breaking a for-loop would not be incompatible.
Unbreakable for-loops are in fact foreach-loops and could be offered
additionally.
function myfunction(a) = [for(i=[n-1]) if (mypredicate(a,i))
break(fn(a,i))] ;
function myfunction1(A) = [foreach (a in A) fn(a)];
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16924.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
On 04/03/2016 03:00 PM, Parkinbot wrote:
If you want us all to write in a functional language, why not offer
at least an echo() function primitive which outputs a value to console
while returning it unchanged.
This is already sitting in the pull request queue. Mainly
the backward compatibility needs to be clarified:
Also a break-function breaking a for-loop would not be incompatible.
Unbreakable for-loops are in fact foreach-loops and could be offered
additionally.
That should be possible with the C++ style for loops which are
in the dev version as experimental feature, see
https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/WIP
ciao,
Torsten.
Torsten has a recently active pull request which extends echo() to work
inside functions. The same pull also adds assert().
https://github.com/openscad/openscad/pull/1587
On Sunday, 3 April 2016, Parkinbot rudolf@parkinbot.com wrote:
why not offer at least
an echo() function primitive which outputs a value to console while
returning it unchanged.
function myfunction(a) = let(x=echo(f(a,i))) // inspect this value
g(x); // calculation some other value from it
function f(a,i) = i;// ... some function
function g(x) = x;// ... some other function
function echo(x) = x; // side effect to console
that is good news. Not exactly, what I meant, but I see a progress.
"each" is flattening a list. "foreach" would be a loop type.
the "echo" function seems to do much more, I had ever expected. But I doubt,
it is a good idea to allow it introducing new variables living in outer
scope. An output expression should not be allowed to have side effects on or
be vulnerable by surrounding code, as in most cases it might be used
temporarly for debugging only.
a = 3;
b = 6;
d = 9, // will be superseeded by echo
t6 = echo(c = 2, d = c + 9) abc*d;
echo(t6 = t6);
...
d = 10; // this anywhere in subsequent code, will superseed the output of
echo
t7 = echo(d=t6+a) d; // this one will win and have a side effect on the
first echo output.
--
View this message in context: http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16927.html
Sent from the OpenSCAD mailing list archive at Nabble.com.
On 04/03/2016 04:43 PM, Parkinbot wrote:
the "echo" function seems to do much more, I had ever expected. But I doubt,
it is a good idea to allow it introducing new variables living in outer
scope. An output expression should not be allowed to have side effects on or
be vulnerable by surrounding code, as in most cases it might be used
temporarly for debugging only.
If you use the variables defined in the echo() call, just replace it
with let() and you'll get the same behavior without the output.
ciao,
Torsten.
On 04/03/2016 04:47 PM, Torsten Paul wrote:
On 04/03/2016 04:43 PM, Parkinbot wrote:
the "echo" function seems to do much more, I had ever expected. But I doubt,
it is a good idea to allow it introducing new variables living in outer
scope. An output expression should not be allowed to have side effects on or
be vulnerable by surrounding code, as in most cases it might be used
temporarly for debugging only.
If you use the variables defined in the echo() call, just replace it
with let() and you'll get the same behavior without the output.
Well, that said, there's still the backward compatibility issue to
solve which is mainly that the new stuff should use sequential
assignment where the existing statement level echo() uses parallel
assignment.
So we might end up creating a new function with a different name
that could have the same behavior for both statement and function
level.
In that case we could have something like...
print(name = "x", value = x);
...which might open wishes for the next can of worms: formatting!!1!
ciao,
Torsten.
I agree with Parkinbot. echo(stuff) expr should not introduce local
variable bindings into expr.
The idea behind this particular syntax is that you can take an existing
function definition,
f(x) =
expr;
and you can temporarily introduce an echo for debugging, without
rearranging the code:
f(x) =
echo(stuff)
expr;
Once you are done debugging, you can delete the echo, or you can comment it
out for later use:
f(x) =
// echo(stuff)
expr;
echo() is a debugging mechanism. Commenting out an echo() should not change
the meaning of the code.
Doug.
On 3 April 2016 at 10:43, Parkinbot rudolf@parkinbot.com wrote:
that is good news. Not exactly, what I meant, but I see a progress.
"each" is flattening a list. "foreach" would be a loop type.
the "echo" function seems to do much more, I had ever expected. But I
doubt,
it is a good idea to allow it introducing new variables living in outer
scope. An output expression should not be allowed to have side effects on
or
be vulnerable by surrounding code, as in most cases it might be used
temporarly for debugging only.
a = 3;
b = 6;
d = 9, // will be superseeded by echo
t6 = echo(c = 2, d = c + 9) abc*d;
echo(t6 = t6);
...
d = 10; // this anywhere in subsequent code, will superseed the output
of
echo
t7 = echo(d=t6+a) d; // this one will win and have a side effect on the
first echo output.
--
View this message in context:
http://forum.openscad.org/Simple-polygon-triangulation-tp16755p16927.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
On 04/03/2016 05:36 PM, doug moen wrote:
echo() is a debugging mechanism. Commenting out an echo() should
not change the meaning of the code.
Well, we can't have echo() then. Any ideas for a good name (which
would then be available as both statement and function and likely
also deprecate the existing echo()).
When we at that, we can also change/remove the existing "ECHO:"
prefix.
Also I'm going to strongly suggest that there will be no
accidental HTML handling as it currently happens due to the
console being able to render that.
ciao,
Torsten.