

This post was updated on .
Hi
Anyone know an easy way to get the area of a polygon? Since a polygon is essentially made up of triangles I presume this is not too hard. Anyone got a solution? I don't want to waste valuable time reinventing the wheel.
Thanks
Ex.


It is absolutely NOT NECESSARY to triangulate in order to compute area. You can walk the perimeter and compute the area of the polygon. Google “Newman’s Scheme”, or directly apply either Green’s Theorem or Stokes’s Theorem, depending on your training.
That was the short, enigmatic answer (Google is your friend). The long answer is:
a) pick a point in the plane  any point. Call it O. [Hint: [0,0] works best] b) for each edge in the polygon, construct the triangle O,v[i], v[i+1 mod n] where n is the number of vertices and the numbering begins with 0 and the edge goes from v[i] to v[i+1 mod n] c) compute the SIGNED area of this triangle and accumulate the sum of the SIGNED areas d) the total SIGNED area will give the area of the polygon and the SIGN will tell you if the polygon is correctly oriented.
A quick and dirty first implementation might stop here. To make it faster:
e) for extra credit, rearrange the terms to eliminate O [that’s why I suggested that [0,0] works best]. You will arrive at a computation that looks something like:
AREA = 0.5 * SUM(v[i] CROSS v[i+1 mod n])
where: the representation of the point v[i] is interpreted as a the representation of a vector from O to v[i] (no, points and vectors are not the same, but SOME computations can pretend that they are) vertices are numbered 0..n1
I have a beautiful proofbydiagram of this, but it does not fit in ASCII format.
If you want to do Volumes, substitute triangular faces of the polyhedron for the edges of the polygon. Use the scalar triple product instead of cross product, and divide by 6 instead of 2.
Proof left as an exercise for the reader.
Green and Stokes showed that you can compute properties of the area by walking around the perimeter. Their followers proved that one follows from the other, and vice versa. Newman did the algebra and optimized the computation in the context of computer graphics.
Newman’s scheme dates to AT LEAST 1965. Green and Stokes were much earlier.
 Kenneth Sloan Vision is the art of seeing what is invisible to others.
Well, first you triangulate the polygon using the code posted to the recent "Simple Polygon Triangulation" thread, then you add up the area of each triangle?
_______________________________________________ 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


If the polygon is in a 2D plane, then the cross product produces a scalar. (At least, that’s the lie they teach to engineers; you may prefer to think of it as a vector in the z direction). If the polygon lies in 3space (but, one hopes, is planar), then, yes  the cross product produces a vector. The crossproduct of two 3D vectors is a vector which has a length twice the area of the triangle defined by the two input vectors. The direction of the cross product is the surface normal.
The accumulated sums of all the crossproducts produced by creating triangles from O and the successive vertices of the polygon is thus a vector with a length twice the area of the polygon (if the polygon is planar) and a direction perpendicular to the plane.
For extra credit: what is the interpretation of the length and direction of the computed vector when the vertices of the “polygon” do not lie in a plane?
Answer (I’m not cruel…): the direction is perpendicular to a “bestfit” plane through the wrinkled “polygon”, and the length is twice the area of the projection of the “polygon” onto that bestfit plane.
So…if you have wrinkled polygons and you want the area of the wrinkled surface…then you *do* have to triangulate. Question: in this case, is the answer unique?  Kenneth Sloan Vision is the art of seeing what is invisible to others.
Nice solution! If I undertand it well the AREA evaluated by the expression in e) is a vector. And the scalar area will be its norm.
_______________________________________________ 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


It should be noted that all this is true for simple polygons.
function area(p) = abs( sum_list( [ for(i=[0:len(p)1]) cross2D( p[i], p[(i+1)%len(p)] ) ] ) )/2;
function sum_list(v, i=0, sum=0) = i==len(v) ? sum : sum_list(v, i+1, sum+v[i] );
function cross2D(p, q) = p[0]*q[1]  p[1]*q[0];
p = 10*[ [1,0], [2,0], [2,1], [1,1] ];
q = 10*[ [1,20], [2,21], [2,20], [1,21] ];
polygon(p);
polygon(q);
echo(area_p=area(p));
echo(area_q=area(q));


Thanks, Kenneth. I translated your algorithm into this:
function area(vs) = let (n = len(vs)) abs(0.5 * sum([for (i=[0:n1]) cross(vs[i], vs[(i+1)%n])]));
function sum(v,i=0) = i < len(v) ? v[i] + sum(v,i+1) : 0;
I used abs() to convert the signed area into an absolute value, since the polygon() module doesn't care about the winding order of the vertices.
This turns out to be the same as Marius's solution, except that it returns an absolute value, and also, the division by 2 happens at the end, instead of once for each vertex.
I ran some tests, and it mostly works as expected. Degenerate polygons work (area==0).
Selfintersecting polygons do not work, which is inconsistent with the polygon() primitive. So the area won't match what polygon() draws on the screen in this case.
Here's that last test case:
//self intersecting quad weird = [ [0,0], [0,1], [1,0], [1,1], ]; echo(area(weird)); // 0 polygon(weird);
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


doug.moen wrote
Selfintersecting polygons do not work, which is inconsistent with the
polygon() primitive. So the area won't match what polygon() draws on the
screen in this case.
To calculate area or even to fill a selfintersecting polygons is ambiguous. The drawing of a selfintersecting polygon done by primitive polygon() is just one of the alternatives. Consider an extrusion of the polygon:
p = [ [0,0], [7,0], [7,7], [3,7], [3,3], [10,3], [10,10], [0,10] ];
polygon(p);
which has a hole in it. If you process that extrusion with slic3r and Simplify3D you will get rather different results. Each one choose its own interpretation of the hole of polygon(). Slic3r does preserve the hole, Simplify3D fill it. And nobody can argue that one of them is wrong.


Hi
Thanks everyone for your answers, I shall play
Ex


Ronaldo gave an example of an "ambiguous" model:  Consider an extrusion of the polygon:
> p = [ [0,0], [7,0], [7,7], [3,7], [3,3], [10,3], [10,10], [0,10] ]; > polygon(p);
which has a hole in it. If you process that extrusion with slic3r and Simplify3D you will get rather different results. Each one choose its own interpretation of the hole of polygon(). Slic3r does preserve the hole, Simplify3D fill it.

We don't really like to tolerate ambiguity in interchange formats for 3D objects. Ideally, everybody should interpret a given file the same way, or report an error that the file is invalid.
Unfortunately, STL doesn't have a well defined interpretation. So what OpenSCAD tries to do is give a warning message when exporting an ambiguous STL file. In your example, there is no warning, and that is a bug. So I reported it: https://github.com/openscad/openscad/issues/1621
In the case of AMF, the standard explicitly forbids self intersection. When you try to export this model to AMF, you get a fatal error message on the console. So that works.
In the case of 3MF, the standard explicitly supports self intersection, and dictates that the Nonzero Winding Number Rule (aka Positive Fill Rule) should be used to disambiguate. We don't support 3MF export yet, but if we did, then this model would not be ambiguous.
Doug.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


Kenneth Sloan wrote
Answer (I’m not cruel…): the direction is perpendicular to a “bestfit” plane through the wrinkled “polygon”, and the length is twice the area of the projection of the “polygon” onto that bestfit plane.
So…if you have wrinkled polygons and you want the area of the wrinkled surface…then you *do* have to triangulate. Question: in this case, is the answer unique?
Is the "bestfit" you mentioned the least square one? Do you suggest any simple criteria to evaluate how much a polygonal in 3D is nonplanar?


Ronaldo wrote
Kenneth Sloan wrote
Answer (I’m not cruel…): the direction is perpendicular to a “bestfit” plane through the wrinkled “polygon”, and the length is twice the area of the projection of the “polygon” onto that bestfit plane.
So…if you have wrinkled polygons and you want the area of the wrinkled surface…then you *do* have to triangulate. Question: in this case, is the answer unique?
Is the "bestfit" you mentioned the least square one? Do you suggest any simple criteria to evaluate how much a polygonal in 3D is nonplanar?
Well, I found this one in case you use eigenvectors/eigenvalues method for least square fitting: the ratio of the minimum eigenvalue (of the covariance matrix) and the maximum one gives a criteria to decide. A zero ratio means the polygonal is planar, a low ratio means it is "almost planar". The ratio grows quadratically with the relative distance from the points to the fitting plane and it is rather independent of the number of points.

