From c5da41ed8ca653a15bf7ff7bf6de4b0b25d4c5b0 Mon Sep 17 00:00:00 2001 From: RonaldoCMP <rcmpersiano@gmail.com> Date: Mon, 1 Nov 2021 04:42:02 +0000 Subject: [PATCH] Correct bugs in polygon_triangulation and _cleave_connected_region( --- geometry.scad | 213 +++++++++++++++++++++++++++------------ tests/test_all.scad | 32 ++++++ tests/test_geometry.scad | 8 +- vectors.scad | 9 +- vnf.scad | 125 ++++++++++++++++++++++- 5 files changed, 307 insertions(+), 80 deletions(-) create mode 100644 tests/test_all.scad diff --git a/geometry.scad b/geometry.scad index 37326e1..c09604b 100644 --- a/geometry.scad +++ b/geometry.scad @@ -38,14 +38,10 @@ function _is_point_on_line(point, line, bounded=false, eps=EPSILON) = t = v0*v1/(v1*v1), bounded = force_list(bounded,2) ) - abs(cross(v0,v1))<eps*norm(v1) + abs(cross(v0,v1))<=eps*norm(v1) && (!bounded[0] || t>=-eps) && (!bounded[1] || t<1+eps) ; -function xis_point_on_line(point, line, bounded=false, eps=EPSILON) = - assert( is_finite(eps) && (eps>=0), "The tolerance should be a non-negative value." ) - point_line_distance(point, line, bounded)<eps; - ///Internal - distance from point `d` to the line passing through the origin with unit direction n ///_dist2line works for any dimension @@ -61,7 +57,67 @@ function _valid_line(line,dim,eps=EPSILON) = function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps); -/// Internal Function: point_left_of_line2d() +/// Internal Function: _is_at_left() +/// Usage: +/// pt = point_left_of_line2d(point, line); +/// Topics: Geometry, Points, Lines +/// Description: +/// Return true iff a 2d point is on or at left of the line defined by `line`. +/// Arguments: +/// pt = The 2d point to check position of. +/// line = Array of two 2d points forming the line segment to test against. +/// eps = Tolerance in the geometrical tests. +function _is_at_left(pt,line,eps=EPSILON) = _tri_class([pt,line[0],line[1]],eps) <= 0; + + +/// Internal Function: _degenerate_tri() +/// Usage: +/// degen = _degenerate_tri(triangle); +/// Topics: Geometry, Triangles +/// Description: +/// Return true for a specific kind of degeneracy: any two triangle vertices are equal +/// Arguments: +/// tri = A list of three 2d points +/// eps = Tolerance in the geometrical tests. +function _degenerate_tri(tri,eps) = + max(norm(tri[0]-tri[1]), norm(tri[1]-tri[2]), norm(tri[2]-tri[0])) < eps ; + + +/// Internal Function: _tri_class() +/// Usage: +/// class = _tri_class(triangle); +/// Topics: Geometry, Triangles +/// Description: +/// Return 1 if the triangle `tri` is CW. +/// Return 0 if the triangle `tri` has colinear vertices. +/// Return -1 if the triangle `tri` is CCW. +/// Arguments: +/// tri = A list of the three 2d vertices of a triangle. +/// eps = Tolerance in the geometrical tests. +function _tri_class(tri, eps=EPSILON) = + let( crx = cross(tri[1]-tri[2],tri[0]-tri[2]) ) + abs( crx ) <= eps*norm(tri[1]-tri[2])*norm(tri[0]-tri[2]) ? 0 : sign( crx ); + + +/// Internal Function: _pt_in_tri() +/// Usage: +/// class = _pt_in_tri(point, tri); +/// Topics: Geometry, Points, Triangles +/// Description: +/// Return 1 if point is inside the triangle interion. +/// Return =0 if point is on the triangle border. +/// Return -1 if point is outside the triangle. +/// Arguments: +/// point = The point to check position of. +/// tri = A list of the three 2d vertices of a triangle. +/// eps = Tolerance in the geometrical tests. +function _pt_in_tri(point, tri, eps=EPSILON) = + min( _tri_class([tri[0],tri[1],point],eps), + _tri_class([tri[1],tri[2],point],eps), + _tri_class([tri[2],tri[0],point],eps) ); + + +/// Internal Function: _point_left_of_line2d() /// Usage: /// pt = point_left_of_line2d(point, line); /// Topics: Geometry, Points, Lines @@ -72,10 +128,11 @@ function _valid_plane(p, eps=EPSILON) = is_vector(p,4) && ! approx(norm(p),0,eps /// Arguments: /// point = The point to check position of. /// line = Array of two points forming the line segment to test against. -function _point_left_of_line2d(point, line) = +function _point_left_of_line2d(point, line, eps=EPSILON) = assert( is_vector(point,2) && is_vector(line*point, 2), "Improper input." ) - cross(line[0]-point, line[1]-line[0]); - +// cross(line[0]-point, line[1]-line[0]); + _tri_class([point,line[1],line[0]],eps); + // Function: is_collinear() // Usage: @@ -1648,18 +1705,19 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // . // The function produce correct triangulations for some non-twisted non-simple polygons. // A polygon is non-twisted iff it is simple or there is a partition of it in -// simple polygons with the same winding. These polygons may have "touching" vertices +// simple polygons with the same winding such that the intersection of any two partitions is +// made of full edges of both partitions. These polygons may have "touching" vertices // (two vertices having the same coordinates, but distinct adjacencies) and "contact" edges // (edges whose vertex pairs have the same pairwise coordinates but are in reversed order) but has -// no self-crossing. See examples bellow. If all polygon edges are contact edges, -// it returns an empty list for 2d polygons and issues an error for 3d polygons. +// no self-crossing. See examples bellow. If all polygon edges are contact edges (polygons with +// zero area), it returns an empty list for 2d polygons and issues an error for 3d polygons. // . -// Self-crossing polygons have no consistent winding and usually produce an error but -// when an error is not issued the outputs are not correct triangulations. The function +// Twisted polygons have no consistent winding and when input to this function usually produce +// an error but when an error is not issued the outputs are not correct triangulations. The function // can work for 3d non-planar polygons if they are close enough to planar but may otherwise // issue an error for this case. // Arguments: -// poly = Array of vertices for the polygon. +// poly = Array of the polygon vertices. // ind = A list indexing the vertices of the polygon in `poly`. // eps = A maximum tolerance in geometrical tests. Default: EPSILON // Example(2D,NoAxes): @@ -1686,7 +1744,7 @@ function point_in_polygon(point, poly, nonzero=false, eps=EPSILON) = // Example(2D,NoAxes): a polygon with "contact" edges and no holes // poly = [ [0,0], [10,0], [10,10], [0,10], [0,0], [3,3], [7,3], // [7,7], [7,3], [3,3] ]; -// tris = polygon_triangulate(poly); // see from the top +// tris = polygon_triangulate(poly); // see from above // color("lightblue") for(tri=tris) polygon(select(poly,tri)); // color("blue") up(1) for(tri=tris) { stroke(select(poly,tri),.15,closed=true); } // color("magenta") up(2) stroke(poly,.25,closed=true); @@ -1703,97 +1761,118 @@ function polygon_triangulate(poly, ind, eps=EPSILON) = assert(is_undef(ind) || (is_vector(ind) && min(ind)>=0 && max(ind)<len(poly) ), "Improper or out of bounds list of indices") + (! is_undef(ind) ) && len(ind) == 0 ? [] : let( ind = is_undef(ind) ? count(len(poly)) : ind ) len(ind) == 3 - ? _is_degenerate([poly[ind[0]], poly[ind[1]], poly[ind[2]]], eps) ? [] : + ? _degenerate_tri([poly[ind[0]], poly[ind[1]], poly[ind[2]]], eps) ? [] : // non zero area assert( norm(scalar_vec3(cross(poly[ind[1]]-poly[ind[0]], poly[ind[2]]-poly[ind[0]]))) > 2*eps, "The polygon vertices are collinear.") [ind] : len(poly[ind[0]]) == 3 - ? // represents the polygon projection on its plane as a 2d polygon + ? // find a representation of the polygon as a 2d polygon by projecting it on its own plane let( ind = deduplicate_indexed(poly, ind, eps) ) len(ind)<3 ? [] : let( pts = select(poly,ind), - nrm = polygon_normal(pts) + nrm = -polygon_normal(pts) ) assert( nrm!=undef, - "The polygon has self-intersections or its vertices are collinear or non coplanar.") + "The polygon has self-intersections or zero area or its vertices are collinear or non coplanar.") let( imax = max_index([for(p=pts) norm(p-pts[0]) ]), v1 = unit( pts[imax] - pts[0] ), v2 = cross(v1,nrm), - prpts = pts*transpose([v1,v2]) + prpts = pts*transpose([v1,v2]) // the 2d projection of pts on the polygon plane ) [for(tri=_triangulate(prpts, count(len(ind)), eps)) select(ind,tri) ] - : let( cw = is_polygon_clockwise(select(poly, ind)) ) - cw - ? [for(tri=_triangulate( poly, reverse(ind), eps )) reverse(tri) ] - : _triangulate( poly, ind, eps ); + : is_polygon_clockwise(select(poly, ind)) + ? _triangulate( poly, ind, eps ) + : [for(tri=_triangulate( poly, reverse(ind), eps )) reverse(tri) ]; -function _triangulate(poly, ind, eps=EPSILON, tris=[]) = +// poly is supposed to be a 2d cw polygon +// implements a modified version of ear cut method for non-twisted polygons +// the polygons accepted by this function are (tecnically) the ones whose interior +// is homeomoph to the interior of a simple polygon +function _triangulate(poly, ind, eps=EPSILON, tris=[]) = len(ind)==3 - ? _is_degenerate(select(poly,ind),eps) - ? tris // last 3 pts perform a degenerate triangle, ignore it + ? _degenerate_tri(select(poly,ind),eps) + ? tris // if last 3 pts perform a degenerate triangle, ignore it : concat(tris,[ind]) // otherwise, include it : let( ear = _get_ear(poly,ind,eps) ) +/* +let( x= [if(is_undef(ear)) echo(ind=ind) 0] ) +is_undef(ear) ? tris : +*/ assert( ear!=undef, - "The polygon has self-intersections or its vertices are collinear or non coplanar.") - is_list(ear) // degenerate ear - ? _triangulate(poly, select(ind,ear[0]+2, ear[0]), eps, tris) // discard it + "The polygon has twists or all its vertices are collinear or non coplanar.") + is_list(ear) // is it a degenerate ear ? + ? len(ind) <= 4 ? tris : + _triangulate(poly, select(ind,ear[0]+3, ear[0]), eps, tris) // discard it : let( ear_tri = select(ind,ear,ear+2), - indr = select(ind,ear+2, ear) // remaining point indices + indr = select(ind,ear+2, ear) // indices of the remaining path ) - _triangulate(poly, indr, eps, concat(tris,[ear_tri])); + _triangulate(poly, indr, eps, concat(tris,[ear_tri])); // a returned ear will be: -// 1. a CCW (non-degenerate) triangle, made of subsequent vertices, without other -// points inside except possibly at its vertices +// 1. a CW non-reflex triangle, made of subsequent poly vertices, without any other +// poly points inside except possibly at its own vertices // 2. or a degenerate triangle where two vertices are coincident // the returned ear is specified by the index of `ind` of its first vertex -function _get_ear(poly, ind, eps, _i=0) = - _i>=len(ind) ? undef : // poly has no ears +function _get_ear(poly, ind, eps, _i=0) = + let( lind = len(ind) ) + lind==3 ? 0 : let( // the _i-th ear candidate p0 = poly[ind[_i]], - p1 = poly[ind[(_i+1)%len(ind)]], - p2 = poly[ind[(_i+2)%len(ind)]] + p1 = poly[ind[(_i+1)%lind]], + p2 = poly[ind[(_i+2)%lind]] ) - // degenerate triangles are returned codified - _is_degenerate([p0,p1,p2],eps) ? [_i] : - // if it is not a convex vertex, check the next one - _is_cw2(p0,p1,p2,eps) ? _get_ear(poly,ind,eps, _i=_i+1) : - let( // vertex p1 is convex - // check if the triangle contains any other point - // except possibly its own vertices - to_tst = select(ind,_i+3, _i-1), - q = [(p0-p2).y, (p2-p0).x], // orthogonal to ray [p0,p2] pointing right - r = [(p2-p1).y, (p1-p2).x], // orthogonal to ray [p2,p1] pointing right - s = [(p1-p0).y, (p0-p1).x], // orthogonal to ray [p1,p0] pointing right - inside = [for(p=select(poly,to_tst)) // for vertices other than p0, p1 and p2 - if( (p-p0)*q<=0 && (p-p2)*r<=0 && (p-p1)*s<=0 // p is on the triangle - && norm(p-p0)>eps // but not on any vertex of it - && norm(p-p1)>eps - && norm(p-p2)>eps ) - p ] + // if vertex p1 is a convex candidate to be an ear, + // check if the triangle [p0,p1,p2] contains any other point + // except possibly p0 and p2 + // exclude the ear candidate central vertex p1 from the verts to check + _tri_class([p0,p1,p2],eps) > 0 + && _none_inside(select(ind,_i+2, _i),poly,p0,p1,p2,eps) ? _i : // found an ear + // otherwise check the next ear candidate + _i<lind-1 ? _get_ear(poly, ind, eps, _i=_i+1) : + // poly has no ears, look for wiskers + let( + wiskers = [for(j=idx(ind)) if(norm(poly[ind[j]]-poly[ind[(j+2)%lind]])<eps) j ] ) - inside==[] ? _i : // found an ear - // check the next ear candidate - _get_ear(poly, ind, eps, _i=_i+1); - - -// true for some specific kinds of degeneracy -function _is_degenerate(tri,eps) = - norm(tri[0]-tri[1])<eps || norm(tri[1]-tri[2])<eps || norm(tri[2]-tri[0])<eps ; - - -function _is_cw2(a,b,c,eps=EPSILON) = cross(a-c,b-c)<eps*norm(a-c)*norm(b-c); + wiskers==[] ? undef : [wiskers[0]]; + + +// returns false ASA it finds some reflex vertex of poly[idxs[.]] +// inside the triangle different from p0 and p2 +// note: to simplify the expressions it is assumed that the input polygon has no twists +function _none_inside(idxs,poly,p0,p1,p2,eps,i=0) = + i>=len(idxs) ? true : + let( + vert = poly[idxs[i]], + prev_vert = poly[select(idxs,i-1)], + next_vert = poly[select(idxs,i+1)] + ) + // check if vert prevent [p0,p1,p2] to be an ear + // this conditions might have a simpler expression + _tri_class([prev_vert, vert, next_vert],eps) <= 0 // reflex condition + && ( // vert is a cw reflex poly vertex inside the triangle [p0,p1,p2] + ( _tri_class([p0,p1,vert],eps)>0 && + _tri_class([p1,p2,vert],eps)>0 && + _tri_class([p2,p0,vert],eps)>=0 ) + // or it is equal to p1 and some of its adjacent edges cross the open segment (p0,p2) + || ( norm(vert-p1) < eps + && ( _is_at_left(p0,[prev_vert,p1],eps) + && _is_at_left(p2,[p1,next_vert],eps) ) + ) + ) + ? false + : _none_inside(idxs,poly,p0,p1,p2,eps,i=i+1); // Function: is_polygon_clockwise() diff --git a/tests/test_all.scad b/tests/test_all.scad new file mode 100644 index 0000000..e51434a --- /dev/null +++ b/tests/test_all.scad @@ -0,0 +1,32 @@ +include <test_affine.scad> +include <test_attachments.scad> +include <test_comparisons.scad> +include <test_coords.scad> +include <test_cubetruss.scad> +include <test_distributors.scad> +include <test_drawing.scad> +include <test_edges.scad> +include <test_fnliterals.scad> +include <test_geometry.scad> +include <test_hull.scad> +include <test_linalg.scad> +include <test_linear_bearings.scad> +include <test_lists.scad> +include <test_math.scad> +include <test_mutators.scad> +include <test_paths.scad> +include <test_quaternions.scad> +include <test_regions.scad> +include <test_rounding.scad> +include <test_screw_drive.scad> +include <test_shapes2d.scad> +include <test_shapes3d.scad> +include <test_skin.scad> +include <test_strings.scad> +include <test_structs.scad> +include <test_transforms.scad> +include <test_trigonometry.scad> +include <test_utility.scad> +include <test_vectors.scad> +include <test_version.scad> +include <test_vnf.scad> diff --git a/tests/test_geometry.scad b/tests/test_geometry.scad index 7eff2bc..1b649d1 100644 --- a/tests/test_geometry.scad +++ b/tests/test_geometry.scad @@ -84,14 +84,14 @@ module test_polygon_triangulate() { poly1 = [ [-10,0,-10], [10,0,10], [0,10,0], [-10,0,-10], [-4,4,-4], [4,4,4], [0,2,0], [-4,4,-4] ]; poly2 = [ [0,0], [5,5], [-5,5], [0,0], [-5,-5], [5,-5] ]; poly3 = [ [0,0], [10,0], [10,10], [10,13], [10,10], [0,10], [0,0], [3,3], [7,3], [7,7], [7,3], [3,3] ]; - tris0 = sort(polygon_triangulate(poly0)); + tris0 = (polygon_triangulate(poly0)); assert(approx(tris0, [[0, 1, 2]])); tris1 = (polygon_triangulate(poly1)); - assert(approx(tris1,( [[2, 3, 4], [6, 7, 0], [2, 4, 5], [6, 0, 1], [1, 2, 5], [5, 6, 1]]))); + assert(approx(tris1,( [[2, 3, 4], [6, 7, 0], [2, 4, 5], [6, 0, 1], [1, 2, 5], [5, 6, 1]]))); tris2 = (polygon_triangulate(poly2)); - assert(approx(tris2,([[0, 1, 2], [3, 4, 5]]))); + assert(approx(tris2,( [[3, 4, 5], [1, 2, 3]]))); tris3 = (polygon_triangulate(poly3)); - assert(approx(tris3,( [[5, 6, 7], [7, 8, 9], [10, 11, 0], [5, 7, 9], [10, 0, 1], [4, 5, 9], [9, 10, 1], [1, 4, 9]]))); + assert(approx(tris3,( [[5, 6, 7], [11, 0, 1], [5, 7, 8], [10, 11, 1], [5, 8, 9], [10, 1, 2], [4, 5, 9], [9, 10, 2]]))); } module test__normalize_plane(){ diff --git a/vectors.scad b/vectors.scad index 7343941..687f0d4 100644 --- a/vectors.scad +++ b/vectors.scad @@ -491,12 +491,11 @@ function _bt_tree(points, ind, leafsize=25) = bounds = pointlist_bounds(select(points,ind)), coord = max_index(bounds[1]-bounds[0]), projc = [for(i=ind) points[i][coord] ], - pmc = mean(projc), - pivot = min_index([for(p=projc) abs(p-pmc)]), + meanpr = mean(projc), + pivot = min_index([for(p=projc) abs(p-meanpr)]), radius = max([for(i=ind) norm(points[ind[pivot]]-points[i]) ]), - median = median(projc), - Lind = [for(i=idx(ind)) if(projc[i]<=median && i!=pivot) ind[i] ], - Rind = [for(i=idx(ind)) if(projc[i] >median && i!=pivot) ind[i] ] + Lind = [for(i=idx(ind)) if(projc[i]<=meanpr && i!=pivot) ind[i] ], + Rind = [for(i=idx(ind)) if(projc[i] >meanpr && i!=pivot) ind[i] ] ) [ ind[pivot], radius, _bt_tree(points, Lind, leafsize), _bt_tree(points, Rind, leafsize) ]; diff --git a/vnf.scad b/vnf.scad index dafad7f..b2325c1 100644 --- a/vnf.scad +++ b/vnf.scad @@ -318,14 +318,13 @@ function vnf_merge(vnfs, cleanup=false, eps=EPSILON) = cleanup? _vnf_cleanup(verts,faces,eps) : [verts,faces]; - function _vnf_cleanup(verts,faces,eps) = let( dedup = vector_search(verts,eps,verts), // collect vertex duplicates map = [for(i=idx(verts)) min(dedup[i]) ], // remap duplic vertices offset = cumsum([for(i=idx(verts)) map[i]==i ? 0 : 1 ]), // remaping face vertex offsets map2 = list(idx(verts))-offset, // map old vertex indices to new indices - nverts = [for(i=idx(verts)) if(map[i]==i) verts[i] ], // eliminates all unreferenced vertices + nverts = [for(i=idx(verts)) if(map[i]==i) verts[i] ], // this doesn't eliminate unreferenced vertices nfaces = [ for(face=faces) let( @@ -388,7 +387,7 @@ function _join_paths_at_vertices(path1,path2,v1,v2) = // Given a region that is connected and has its outer border in region[0], // produces a polygon with the same points that has overlapping connected paths // to join internal holes to the outer border. Output is a single path. -function _cleave_connected_region(region) = +function _old_cleave_connected_region(region) = len(region)==0? [] : len(region)<=1? clockwise_polygon(region[0]) : let( @@ -411,8 +410,126 @@ function _cleave_connected_region(region) = ] ) assert(len(orgn)<len(region)) - _cleave_connected_region(orgn); + _old_cleave_connected_region(orgn); +/// Internal Function: _cleave_connected_region(region, eps) +/// Description: +/// Given a region that is connected and has its outer border in region[0], +/// produces a polygon with the same points that has overlapping connected paths +/// to join internal holes to the outer border. Output is a single path. +/// It expect that region[0] be a simple closed CW path and that each hole, +/// region[i] for i>0, be a simple closed CCW path. +/// The paths are also supposed to be disjoint except for common vertices and +/// common edges but no crossing. +/// This function implements an extension of the algorithm discussed in: +/// https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf +function _cleave_connected_region(region, eps=EPSILON) = + len(region)==1 ? region[0] : + let( + outer = deduplicate(clockwise_polygon(region[0])), // + holes = [for(i=[1:1:len(region)-1]) // possibly unneeded + let(poly=region[i]) // + deduplicate( ccw_polygon(poly) ) ], // + extridx = [for(li=holes) max_index(column(li,0)) ], + // the right extreme vertex for each hole sorted by decreasing x values + extremes = sort( [for(i=idx(holes)) [ i, extridx[i], -holes[i][extridx[i]].x] ], idx=2 ) + ) + _polyHoles(outer, holes, extremes, eps, 0); + + +// connect the hole paths one at a time to the outer path. +// 'extremes' is the list of the right extreme vertex of each hole sorted by decreasing abscissas +// see _cleave_connected_region(region, eps) +function _polyHoles(outer, holes, extremes, eps=EPSILON, n=0) = + let( + extr = extremes[n], // + hole = holes[extr[0]], // hole path to bridge to the outer path + ipt = extr[1], // index of the hole point with maximum abscissa + brdg = _bridge(hole[ipt], outer, eps) // the index of a point in outer to bridge hole[ipt] to + ) + assert(brdg!=undef, "Error: check input polygon restrictions") + let( + l = len(outer), + lh = len(hole), + // the new outer polygon bridging the hole to the old outer + npoly = + approx(outer[brdg], hole[ipt], eps) + ? [ for(i=[brdg: 1: brdg+l]) outer[i%l] , + for(i=[ipt+1:1: ipt+lh-1]) hole[i%lh] ] + : [ for(i=[brdg: 1: brdg+l]) outer[i%l] , + for(i=[ipt:1: ipt+lh]) hole[i%lh] ] + ) + n==len(holes)-1 ? npoly : + _polyHoles(npoly, holes, extremes, eps, n+1); + +// find a point in outer to be connected to pt in the interior of outer +// by a segment that not cross or touch any non adjacente edge of outer. +// return the index of a vertex in the outer path where the bridge should end +// see _polyHoles(outer, holes, extremes, eps) +function _bridge(pt, outer,eps) = + // find the intersection of a ray from pt to the right + // with the boundary of the outer cycle + let( + l = len(outer), + crxs = + [for( i=idx(outer) ) + let( edge = select(outer,i,i+1) ) + // consider just descending outer edges at right of pt crossing ordinate pt.y + if( (edge[0].y> pt.y) + && (edge[1].y<=pt.y) + && ( norm(edge[1]-pt)<eps // accepts touching vertices + || _tri_class([pt, edge[0], edge[1]], eps)>0 ) ) + [ i, + // the point of edge with ordinate pt.y + abs(pt.y-edge[1].y)<eps ? edge[1] : + let( u = (pt-edge[1]).y / (edge[0]-edge[1]).y ) + (1-u)*edge[1] + u*edge[0] + ] + ] + ) + assert(crxs!=[], "Error: check input polygon restrictions") + let( + // the intersection point nearest to pt + minX = min([for(p=crxs) p[1].x]), + crxcand = [for(crx=crxs) if(crx[1].x < minX+eps) crx ], + nearest = min_index([for(crx=crxcand) outer[crx[0]].y]), + proj = crxcand[nearest], + vert0 = outer[proj[0]], // the two vertices of the nearest crossing edge + vert1 = outer[(proj[0]+1)%l], + isect = proj[1] // the intersection point + ) + // if pt touches the middle of an outer edge -> error + assert( ! approx(pt,isect,eps) || approx(pt,vert0,eps) || approx(pt,vert1,eps), + "There is a forbidden self_intersection" ) + norm(pt-vert0) < eps ? proj[0] : // if pt touches an outer vertex, return its index + norm(pt-vert1) < eps ? (proj[0]+1)%l : + let( + // the edge [vert0, vert1] necessarily satisfies vert0.y > vert1.y + // indices of candidates to an outer bridge point + cand = + (vert0.x > pt.x) + ? [ proj[0], + // select reflex vertices inside of the triangle [pt, vert0, isect] + for(i=idx(outer)) + if( _tri_class(select(outer,i-1,i+1),eps) <= 0 + && _pt_in_tri(outer[i], [pt, vert0, isect], eps)>=0 ) + i + ] + : [ (proj[0]+1)%l, + // select reflex vertices inside of the triangle [pt, isect, vert1] + for(i=idx(outer)) + if( _tri_class(select(outer,i-1,i+1),eps) <= 0 + && _pt_in_tri(outer[i], [pt, isect, vert1], eps)>=0 ) + i + ], + // choose the candidate outer[i] such that the line [pt, outer[i]] has minimum slope + // among those with minimum slope choose the nearest to pt + slopes = [for(i=cand) 1-abs(outer[i].x-pt.x)/norm(outer[i]-pt) ], + min_slp = min(slopes), + cand2 = [for(i=idx(cand)) if(slopes[i]<=min_slp+eps) cand[i] ], + nearest = min_index([for(i=cand2) norm(pt-outer[i]) ]) + ) + cand2[nearest]; // Function: vnf_from_region()