Fix bug where polygon_normal and plane_from_polygon sometimes give a

wrong-pointing normal.
Rename project_onto_plane to plane_closest_point for consistency with
line_closest_point.  (Also this distinguishes it better from
project_plane.
Other misc doc tweaks
This commit is contained in:
Adrian Mariano 2021-09-11 16:42:40 -04:00
parent 5475c1650c
commit 60dbf8c73e
3 changed files with 66 additions and 76 deletions

View file

@ -155,15 +155,18 @@ function line_normal(p1,p2) =
// of the intersection point along s2. The proportional values run over
// the range of 0 to 1 for each segment, so if it is in this range, then
// the intersection lies on the segment. Otherwise it lies somewhere on
// the extension of the segment. Result is undef for coincident lines.
// the extension of the segment. If lines are parallel or coincident then
// it returns undef.
function _general_line_intersection(s1,s2,eps=EPSILON) =
let(
denominator = det2([s1[0],s2[0]]-[s1[1],s2[1]])
) approx(denominator,0,eps=eps)? [undef,undef,undef] :
)
approx(denominator,0,eps=eps) ? undef :
let(
t = det2([s1[0],s2[0]]-s2) / denominator,
u = det2([s1[0],s1[0]]-[s2[0],s1[1]]) / denominator
) [s1[0]+t*(s1[1]-s1[0]), t, u];
)
[s1[0]+t*(s1[1]-s1[0]), t, u];
@ -215,7 +218,7 @@ function line_intersection(line1, line2, bounded1, bounded2, bounded, eps=EPSILO
assert( is_undef(bounded1) || is_bool(bounded1) || is_bool_list(bounded1,2), "Invalid value for \"bounded1\"")
assert( is_undef(bounded2) || is_bool(bounded2) || is_bool_list(bounded2,2), "Invalid value for \"bounded2\"")
let(isect = _general_line_intersection(line1,line2,eps=eps))
is_undef(isect[0]) ? undef :
is_undef(isect) ? undef :
let(
bounded1 = force_list(first_defined([bounded1,bounded,false]),2),
bounded2 = force_list(first_defined([bounded2,bounded,false]),2),
@ -390,7 +393,8 @@ function plane3pt_indexed(points, i1, i2, i3) =
// plane = plane_from_normal(normal, [pt])
// Topics: Geometry, Planes
// Description:
// Returns a plane defined by a normal vector and a point.
// Returns a plane defined by a normal vector and a point. If you omit `pt` you will get a plane
// passing through the origin.
// Arguments:
// normal = Normal vector to the plane to find.
// pt = Point 3D on the plane to find.
@ -503,10 +507,13 @@ function plane_from_points(points, fast=false, eps=EPSILON) =
function plane_from_polygon(poly, fast=false, eps=EPSILON) =
assert( is_path(poly,dim=3), "Invalid polygon." )
assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." )
len(poly)==3 ? plane3pt(poly[0],poly[1],poly[2]) :
let( triple = sort(noncollinear_triple(poly,error=false)) )
triple==[] ? [] :
let( plane = plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]]))
let(
poly_normal = polygon_normal(poly)
)
is_undef(poly_normal) ? [] :
let(
plane = plane_from_normal(poly_normal, poly[0])
)
fast? plane: points_on_plane(poly, plane, eps=eps)? plane: [];
@ -539,13 +546,13 @@ function plane_offset(plane) =
// Function: projection_on_plane()
// Function: plane_closest_point()
// Usage:
// pts = projection_on_plane(plane, points);
// pts = plane_closest_point(plane, points);
// Topics: Geometry, Planes, Projection
// Description:
// Given a plane definition `[A,B,C,D]`, where `Ax+By+Cz=D`, and a list of 2d or
// 3d points, return the 3D orthogonal projection of the points on the plane.
// 3d points, return the closest 3D orthogonal projection of the points on the plane.
// In other words, for every point given, returns the closest point to it on the plane.
// Arguments:
// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
@ -553,7 +560,7 @@ function plane_offset(plane) =
// Example(FlatSpin,VPD=500,VPT=[2,20,10]):
// points = move([10,20,30], p=yrot(25, p=path3d(circle(d=100, $fn=36))));
// plane = plane_from_normal([1,0,1]);
// proj = projection_on_plane(plane,points);
// proj = plane_closest_point(plane,points);
// color("red") move_copies(points) sphere(d=2,$fn=12);
// color("blue") move_copies(proj) sphere(d=2,$fn=12);
// move(centroid(proj)) {
@ -562,30 +569,15 @@ function plane_offset(plane) =
// %cube([120,150,0.1],center=true);
// }
// }
function projection_on_plane(plane, points) =
function plane_closest_point(plane, points) =
is_vector(points,3) ? plane_closest_point(plane,[points])[0] :
assert( _valid_plane(plane), "Invalid plane." )
assert( is_matrix(points,undef,3), "Invalid list of points or dimension." )
assert( is_matrix(points,undef,3), "Must supply 3D points.")
let(
p = len(points[0])==2
? [for(pi=points) point3d(pi) ]
: points,
plane = normalize_plane(plane),
n = point3d(plane)
)
[for(pi=p) pi - (pi*n - plane[3])*n];
// Function: plane_point_nearest_origin()
// Usage:
// pt = plane_point_nearest_origin(plane);
// Topics: Geometry, Planes, Distance
// Description:
// Returns the point on the plane that is closest to the origin.
// Arguments:
// plane = The `[A,B,C,D]` plane definition where `Ax+By+Cz=D` is the formula of the plane.
function plane_point_nearest_origin(plane) =
let( plane = normalize_plane(plane) )
point3d(plane) * plane[3];
[for(pi=points) pi - (pi*n - plane[3])*n];
// Function: point_plane_distance()
@ -609,7 +601,7 @@ function point_plane_distance(plane, point) =
point3d(plane)* point - plane[3];
// Returns [POINT, U] if line intersects plane at one point.
// Returns [POINT, U] if line intersects plane at one point, where U is zero at line[0] and 1 at line[1]
// Returns [LINE, undef] if the line is on the plane.
// Returns undef if line is parallel to, but not on the given plane.
function _general_plane_line_intersection(plane, line, eps=EPSILON) =
@ -700,10 +692,10 @@ function plane_line_intersection(plane, line, bounded=false, eps=EPSILON) =
function polygon_line_intersection(poly, line, bounded=false, eps=EPSILON) =
assert( is_finite(eps) && eps>=0, "The tolerance should be a positive number." )
assert(is_path(poly,dim=3), "Invalid polygon." )
assert(!is_list(bounded) || len(bounded)==2, "Invalid bound condition(s).")
assert(is_bool(bounded) || is_bool_list(bounded,2), "Invalid bound condition.")
assert(_valid_line(line,dim=3,eps=eps), "Invalid 3D line." )
let(
bounded = is_list(bounded)? bounded : [bounded, bounded],
bounded = force_list(bounded,2),
poly = deduplicate(poly),
indices = noncollinear_triple(poly)
) indices==[] ? undef :
@ -1194,7 +1186,11 @@ function circle_line_intersection(c,r,d,line,bounded=false,eps=EPSILON) =
// test = noncollinear_triple(points);
// Topics: Geometry, Noncollinearity
// Description:
// Finds the indices of three good non-collinear points from the pointlist `points`.
// Finds the indices of three non-collinear points from the pointlist `points`.
// It selects two well separated points to define a line and chooses the third point
// to be the point farthest off the line. The points do not necessarily having the
// same winding direction as the polygon so they cannot be used to determine the
// winding direction or the direction of the normal.
// If all points are collinear returns [] when `error=true` or an error otherwise .
// Arguments:
// points = List of input points.
@ -1590,17 +1586,21 @@ function reverse_polygon(poly) =
// vec = polygon_normal(poly);
// Topics: Geometry, Polygons
// Description:
// Given a 3D planar polygon, returns a unit-length normal vector for the
// clockwise orientation of the polygon. If the polygon points are collinear, returns `undef`.
// It doesn't check for coplanarity.
// Given a 3D simple planar polygon, returns a unit normal vector for the polygon. The vector
// is oriented so that if the normal points towards the viewer, the polygon winds in the clockwise
// direction. If the polygon has zero area, returns `undef`. If the polygon is self-intersecting
// the the result is undefined. It doesn't check for coplanarity.
// Arguments:
// poly = The list of 3D path points for the perimeter of the polygon.
function polygon_normal(poly) =
assert(is_path(poly,dim=3), "Invalid 3D polygon." )
len(poly)==3 ? point3d(plane3pt(poly[0],poly[1],poly[2])) :
let( triple = sort(noncollinear_triple(poly,error=false)) )
triple==[] ? undef :
point3d(plane3pt(poly[triple[0]],poly[triple[1]],poly[triple[2]])) ;
let(
L=len(poly),
area_vec = sum([for(i=idx(poly))
cross(poly[(i+1)%L]-poly[0],
poly[(i+2)%L]-poly[(i+1)%L])])
)
area_vec==0 ? undef : unit(-area_vec);
function _split_polygon_at_x(poly, x) =

View file

@ -26,8 +26,7 @@ test_plane_from_points();
test_plane_from_polygon();
test_plane_normal();
test_plane_offset();
test_projection_on_plane();
test_plane_point_nearest_origin();
test_plane_closest_point();
test_point_plane_distance();
test__general_plane_line_intersection();
@ -149,16 +148,6 @@ module test_plane_intersection(){
*test_plane_intersection();
module test_plane_point_nearest_origin(){
point = rands(-1,1,3)+[2,0,0]; // a non zero vector
plane = [ each point, point*point]; // a plane containing `point`
info = info_str([["point = ",point],["plane = ",plane]]);
assert_approx(plane_point_nearest_origin(plane),point,info);
assert_approx(plane_point_nearest_origin([each point,5]),5*unit(point)/norm(point),info);
}
test_plane_point_nearest_origin();
module test_plane_offset(){
plane = rands(-1,1,4)+[2,0,0,0]; // a valid plane
info = info_str([["plane = ",plane]]);
@ -248,14 +237,14 @@ module test_points_on_plane() {
ang = rands(0,360,1)[0];
normal = rot(a=ang,p=normal0);
plane = [each normal, normal*dir];
prj_pts = projection_on_plane(plane,pts);
prj_pts = plane_closest_point(plane,pts);
info = info_str([["pts = ",pts],["dir = ",dir],["ang = ",ang]]);
assert(points_on_plane(prj_pts,plane),info);
assert(!points_on_plane(concat(pts,[normal-dir]),plane),info);
}
*test_points_on_plane();
module test_projection_on_plane(){
module test_plane_closest_point(){
ang = rands(0,360,1)[0];
dir = rands(-10,10,3);
normal0 = unit([1,2,3]);
@ -265,16 +254,16 @@ module test_projection_on_plane(){
planem = [each normal, normal*dir];
pts = [for(i=[1:10]) rands(-1,1,3)];
info = info_str([["ang = ",ang],["dir = ",dir]]);
assert_approx( projection_on_plane(plane,pts),
projection_on_plane(plane,projection_on_plane(plane,pts)),info);
assert_approx( projection_on_plane(plane,pts),
rot(a=ang,p=projection_on_plane(plane0,rot(a=-ang,p=pts))),info);
assert_approx( move((-normal*dir)*normal,p=projection_on_plane(planem,pts)),
projection_on_plane(plane,pts),info);
assert_approx( move((normal*dir)*normal,p=projection_on_plane(plane,pts)),
projection_on_plane(planem,pts),info);
assert_approx( plane_closest_point(plane,pts),
plane_closest_point(plane,plane_closest_point(plane,pts)),info);
assert_approx( plane_closest_point(plane,pts),
rot(a=ang,p=plane_closest_point(plane0,rot(a=-ang,p=pts))),info);
assert_approx( move((-normal*dir)*normal,p=plane_closest_point(planem,pts)),
plane_closest_point(plane,pts),info);
assert_approx( move((normal*dir)*normal,p=plane_closest_point(plane,pts)),
plane_closest_point(planem,pts),info);
}
*test_projection_on_plane();
*test_plane_closest_point();
module test_line_from_points() {
assert_approx(line_from_points([[1,0],[0,0],[-1,0]]),[[-1,0],[1,0]]);

View file

@ -130,6 +130,7 @@ function is_only_noncolinear_vertex(points, facelist, vertex) =
// face = The face, given as a list of indices into the vertex array `points`.
function triangulate_face(points, face) =
let(
points = path3d(points),
face = deduplicate_indexed(points,face),
count = len(face)
)
@ -158,21 +159,21 @@ function triangulate_face(points, face) =
(clipable_ear)? // There is no point inside the ear.
is_only_noncolinear_vertex(points, face, cv)?
// In the point&line degeneracy clip to somewhere in the middle of the line.
flatten([
triangulate_face(points, select(face, cv, (cv+2)%count)),
triangulate_face(points, select(face, (cv+2)%count, cv))
])
concat(
triangulate_face(points, select(face, cv, (cv+2)%count)),
triangulate_face(points, select(face, (cv+2)%count, cv))
)
:
// Otherwise the ear is safe to clip.
flatten([
[select(face, pv, nv)],
triangulate_face(points, select(face, nv, pv))
])
[
select(face, pv, nv),
each triangulate_face(points, select(face, nv, pv))
]
: // If there is a point inside the ear, make a diagonal and clip along that.
flatten([
concat(
triangulate_face(points, select(face, cv, diagonal_point)),
triangulate_face(points, select(face, diagonal_point, cv))
]);
);
// Function: triangulate_faces()