Minor corrections

Corrections and simplifications of  polygon_area(), centroid() and point_in_polygon.
This commit is contained in:
RonaldoCMP
2021-03-25 01:57:07 +00:00
parent 4d9bf2d028
commit 80e6cf6316
2 changed files with 29 additions and 24 deletions

View File

@@ -1670,26 +1670,26 @@ function furthest_point(pt, points) =
// area = polygon_area(poly); // area = polygon_area(poly);
// Description: // Description:
// Given a 2D or 3D planar polygon, returns the area of that polygon. // Given a 2D or 3D planar polygon, returns the area of that polygon.
// If the polygon is self-crossing, the results are undefined. For non-planar points the result is undef. // If the polygon is self-crossing, the results are undefined. For non-planar 3D polygon the result is undef.
// When `signed` is true, a signed area is returned; a positive area indicates a counterclockwise polygon. // When `signed` is true, a signed area is returned; a positive area indicates a clockwise polygon.
// Arguments: // Arguments:
// poly = polygon to compute the area of. // poly = polygon to compute the area of.
// signed = if true, a signed area is returned (default: false) // signed = if true, a signed area is returned (default: false)
function polygon_area(poly, signed=false) = function polygon_area(poly, signed=false) =
assert(is_path(poly), "Invalid polygon." ) assert(is_path(poly), "Invalid polygon." )
len(poly)<3 ? 0 : len(poly)<3 ? 0 :
let( cpoly = close_path(simplify_path(poly)) )
len(poly[0])==2 len(poly[0])==2
? sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 ? let( total = sum([for(i=[1:1:len(poly)-2]) cross(poly[i]-poly[0],poly[i+1]-poly[0]) ])/2 )
signed ? total : abs(total)
: let( plane = plane_from_points(poly) ) : let( plane = plane_from_points(poly) )
plane==undef? undef : plane==undef? undef :
let( let(
n = unit(plane_normal(plane)), n = plane_normal(plane),
total = sum([ total = sum([
for(i=[1:1:len(cpoly)-2]) for(i=[1:1:len(poly)-2])
let( let(
v1 = cpoly[i] - cpoly[0], v1 = poly[i] - poly[0],
v2 = cpoly[i+1] - cpoly[0] v2 = poly[i+1] - poly[0]
) )
cross(v1,v2) * n cross(v1,v2) * n
])/2 ])/2
@@ -1853,7 +1853,9 @@ function centroid(poly) =
for(i=[0:len(poly)-1]) for(i=[0:len(poly)-1])
let(segment=select(poly,i,i+1)) let(segment=select(poly,i,i+1))
det2(segment)*sum(segment) det2(segment)*sum(segment)
]) / 6 / polygon_area(poly) ]) / 6 / polygon_area(poly,signed=true)
// polygon_area(concat([[0,0]],segment),signed=true)*sum(segment)
// ]) / 3 / polygon_area(poly,signed=true)
: let( plane = plane_from_points(poly, fast=true) ) : let( plane = plane_from_points(poly, fast=true) )
assert( !is_undef(plane), "The polygon must be planar." ) assert( !is_undef(plane), "The polygon must be planar." )
let( let(
@@ -1879,7 +1881,7 @@ function centroid(poly) =
// the specified 2D polygon using either the Nonzero Winding rule or the Even-Odd rule. // the specified 2D polygon using either the Nonzero Winding rule or the Even-Odd rule.
// See https://en.wikipedia.org/wiki/Nonzero-rule and https://en.wikipedia.org/wiki/Evenodd_rule. // See https://en.wikipedia.org/wiki/Nonzero-rule and https://en.wikipedia.org/wiki/Evenodd_rule.
// The polygon is given as a list of 2D points, not including the repeated end point. // The polygon is given as a list of 2D points, not including the repeated end point.
// Returns -1 if the point is outside the polyon. // Returns -1 if the point is outside the polygon.
// Returns 0 if the point is on the boundary. // Returns 0 if the point is on the boundary.
// Returns 1 if the point lies in the interior. // Returns 1 if the point lies in the interior.
// The polygon does not need to be simple: it can have self-intersections. // The polygon does not need to be simple: it can have self-intersections.
@@ -1890,7 +1892,7 @@ function centroid(poly) =
// poly = The list of 2D path points forming the perimeter of the polygon. // poly = The list of 2D path points forming the perimeter of the polygon.
// nonzero = The rule to use: true for "Nonzero" rule and false for "Even-Odd" (Default: true ) // nonzero = The rule to use: true for "Nonzero" rule and false for "Even-Odd" (Default: true )
// eps = Acceptable variance. Default: `EPSILON` (1e-9) // eps = Acceptable variance. Default: `EPSILON` (1e-9)
function point_in_polygon(point, poly, eps=EPSILON, nonzero=true) = function point_in_polygon(point, poly, nonzero=true, eps=EPSILON) =
// Original algorithms from http://geomalgorithms.com/a03-_inclusion.html // Original algorithms from http://geomalgorithms.com/a03-_inclusion.html
assert( is_vector(point,2) && is_path(poly,dim=2) && len(poly)>2, assert( is_vector(point,2) && is_path(poly,dim=2) && len(poly)>2,
"The point and polygon should be in 2D. The polygon should have more that 2 points." ) "The point and polygon should be in 2D. The polygon should have more that 2 points." )
@@ -1923,12 +1925,12 @@ function point_in_polygon(point, poly, eps=EPSILON, nonzero=true) =
p0 = poly[i]-point, p0 = poly[i]-point,
p1 = poly[(i+1)%n]-point p1 = poly[(i+1)%n]-point
) )
if( ( (p1.y>eps && p0.y<=0) || (p1.y<=0 && p0.y>eps) ) if( ( (p1.y>eps && p0.y<=eps) || (p1.y<=eps && p0.y>eps) )
&& 0 < p0.x - p0.y *(p1.x - p0.x)/(p1.y - p0.y) ) && -eps < p0.x - p0.y *(p1.x - p0.x)/(p1.y - p0.y) )
1 1
] ]
) )
2*(len(cross)%2)-1;; 2*(len(cross)%2)-1;
// Function: polygon_is_clockwise() // Function: polygon_is_clockwise()
@@ -1979,7 +1981,7 @@ function reverse_polygon(poly) =
// n = polygon_normal(poly); // n = polygon_normal(poly);
// Description: // Description:
// Given a 3D planar polygon, returns a unit-length normal vector for the // Given a 3D planar polygon, returns a unit-length normal vector for the
// clockwise orientation of the polygon. // clockwise orientation of the polygon. If the polygon points are collinear, returns `undef`.
function polygon_normal(poly) = function polygon_normal(poly) =
assert(is_path(poly,dim=3), "Invalid 3D polygon." ) assert(is_path(poly,dim=3), "Invalid 3D polygon." )
let( let(
@@ -1989,7 +1991,7 @@ function polygon_normal(poly) =
for (i=[1:1:len(poly)-2]) for (i=[1:1:len(poly)-2])
cross(poly[i+1]-p0, poly[i]-p0) cross(poly[i+1]-p0, poly[i]-p0)
]) ])
) unit(n); ) unit(n,undef);
function _split_polygon_at_x(poly, x) = function _split_polygon_at_x(poly, x) =

View File

@@ -842,7 +842,7 @@ module test_cleanup_path() {
module test_polygon_area() { module test_polygon_area() {
assert(approx(polygon_area([[1,1],[-1,1],[-1,-1],[1,-1]]), 4)); assert(approx(polygon_area([[1,1],[-1,1],[-1,-1],[1,-1]]), 4));
assert(approx(polygon_area(circle(r=50,$fn=1000)), -PI*50*50, eps=0.1)); assert(approx(polygon_area(circle(r=50,$fn=1000),signed=true), -PI*50*50, eps=0.1));
} }
*test_polygon_area(); *test_polygon_area();
@@ -914,7 +914,7 @@ module test_noncollinear_triple() {
module test_centroid() { module test_centroid() {
$fn = 24; $fn = 24;
assert_approx(centroid(circle(d=100)), [0,0]); assert_approx(centroid(circle(d=100)), [0,0]);
assert_approx(centroid(rect([40,60],rounding=10,anchor=LEFT)), [20,0]); assert_approx(centroid(rect([40,60],rounding=10,anchor=LEFT)), [-20,0]);
assert_approx(centroid(rect([40,60],rounding=10,anchor=FWD)), [0,30]); assert_approx(centroid(rect([40,60],rounding=10,anchor=FWD)), [0,30]);
poly = [for(a=[0:90:360]) poly = [for(a=[0:90:360])
move([1,2.5,3.1], rot(p=[cos(a),sin(a),0],from=[0,0,1],to=[1,1,1])) ]; move([1,2.5,3.1], rot(p=[cos(a),sin(a),0],from=[0,0,1],to=[1,1,1])) ];
@@ -943,19 +943,22 @@ module test_point_in_polygon() {
poly2 = [ [-3,-3],[2,-3],[2,1],[-1,1],[-1,-1],[1,-1],[1,2],[-3,2] ]; poly2 = [ [-3,-3],[2,-3],[2,1],[-1,1],[-1,-1],[1,-1],[1,2],[-3,2] ];
assert(point_in_polygon([0,0], poly) == 1); assert(point_in_polygon([0,0], poly) == 1);
assert(point_in_polygon([20,0], poly) == -1); assert(point_in_polygon([20,0], poly) == -1);
assert(point_in_polygon([20,0], poly,EPSILON,nonzero=false) == -1); assert(point_in_polygon([20,0], poly,nonzero=false) == -1);
assert(point_in_polygon([5,5], poly) == 1); assert(point_in_polygon([5,5], poly) == 1);
assert(point_in_polygon([-5,5], poly) == 1); assert(point_in_polygon([-5,5], poly) == 1);
assert(point_in_polygon([-5,-5], poly) == 1); assert(point_in_polygon([-5,-5], poly) == 1);
assert(point_in_polygon([5,-5], poly) == 1); assert(point_in_polygon([5,-5], poly) == 1);
assert(point_in_polygon([5,-5], poly,EPSILON,nonzero=false) == 1); assert(point_in_polygon([5,-5], poly,nonzero=false,eps=EPSILON) == 1);
assert(point_in_polygon([-10,-10], poly) == -1); assert(point_in_polygon([-10,-10], poly) == -1);
assert(point_in_polygon([10,0], poly) == 0); assert(point_in_polygon([10,0], poly) == 0);
assert(point_in_polygon([0,10], poly) == 0); assert(point_in_polygon([0,10], poly) == 0);
assert(point_in_polygon([0,-10], poly) == 0); assert(point_in_polygon([0,-10], poly) == 0);
assert(point_in_polygon([0,-10], poly,EPSILON,nonzero=false) == 0); assert(point_in_polygon([0,-10], poly,nonzero=false) == 0);
assert(point_in_polygon([0,0], poly2,EPSILON,nonzero=true) == 1); assert(point_in_polygon([0,0], poly2,nonzero=true) == 1);
assert(point_in_polygon([0,0], poly2,EPSILON,nonzero=false) == -1); assert(point_in_polygon([0,1], poly2,nonzero=true) == 0);
assert(point_in_polygon([0,1], poly2,nonzero=false) == 0);
assert(point_in_polygon([1,0], poly2,nonzero=false) == 0);
assert(point_in_polygon([0,0], poly2,nonzero=false,eps=EPSILON) == -1);
} }
*test_point_in_polygon(); *test_point_in_polygon();