Coverage for trimesh/path/polygons.py: 85%

1import numpy as np 

2from shapely import ops 

3from shapely.geometry import Polygon 

4 

5from .. import bounds, geometry, graph, grouping 

6from ..constants import log 

7from ..constants import tol_path as tol 

8from ..iteration import reduce_cascade 

9from ..transformations import transform_points 

10from ..typed import ArrayLike, Iterable, NDArray, Number, Seed, float64, int64 

11from ..util import random_generator 

12from .simplify import fit_circle_check 

13from .traversal import resample_path 

14 

15try: 

16 import networkx as nx 

17except BaseException as E: 

18 # create a dummy module which will raise the ImportError 

19 # or other exception only when someone tries to use networkx 

20 from ..exceptions import ExceptionWrapper 

21 

22 nx = ExceptionWrapper(E) 

23try: 

24 from rtree.index import Index 

25except BaseException as E: 

26 # create a dummy module which will raise the ImportError 

27 from ..exceptions import ExceptionWrapper 

28 

29 Index = ExceptionWrapper(E) 

30 

31 

32def enclosure_tree(polygons): 

33 """ 

34 Given a list of shapely polygons with only exteriors, 

35 find which curves represent the exterior shell or root curve 

36 and which represent holes which penetrate the exterior. 

37 

38 This is done with an R-tree for rough overlap detection, 

39 and then exact polygon queries for a final result. 

40 

41 Parameters 

42 ----------- 

43 polygons : (n,) shapely.geometry.Polygon 

44 Polygons which only have exteriors and may overlap 

45 

46 Returns 

47 ----------- 

48 roots : (m,) int 

49 Index of polygons which are root 

50 contains : networkx.DiGraph 

51 Edges indicate a polygon is 

52 contained by another polygon 

53 """ 

54 

55 # nodes are indexes in polygons 

56 contains = nx.DiGraph() 

57 

58 if len(polygons) == 0: 

59 return np.array([], dtype=np.int64), contains 

60 elif len(polygons) == 1: 

61 # `paths_to_polygons` may produce `None` for unrecoverable 

62 # geometry: never emit it as a root, matching the multi-polygon 

63 # code path below where `None` is excluded by the bounds check 

64 if polygons[0] is None: 

65 return np.array([], dtype=np.int64), contains 

66 # add an early exit for only a single polygon 

67 contains.add_node(0) 

68 return np.array([0], dtype=np.int64), contains 

69 

70 # get the bounds for every valid polygon 

71 bounds = { 

72 k: v 

73 for k, v in { 

74 i: getattr(polygon, "bounds", []) for i, polygon in enumerate(polygons) 

75 }.items() 

76 if len(v) == 4 

77 } 

78 

79 # make sure we don't have orphaned polygon 

80 contains.add_nodes_from(bounds.keys()) 

81 

82 if len(bounds) > 0: 

83 # if there are no valid bounds tree creation will fail 

84 # and we won't be calling `tree.intersection` anywhere 

85 # we could return here but having multiple return paths 

86 # seems more dangerous than iterating through an empty graph 

87 tree = Index(zip(bounds.keys(), bounds.values(), [None] * len(bounds))) 

88 

89 # loop through every polygon 

90 for i, b in bounds.items(): 

91 # we first query for bounding box intersections from the R-tree 

92 for j in tree.intersection(b): 

93 # if we are checking a polygon against itself continue 

94 if i == j: 

95 continue 

96 # do a more accurate polygon in polygon test 

97 # for the enclosure tree information 

98 if polygons[i].contains(polygons[j]): 

99 contains.add_edge(i, j) 

100 elif polygons[j].contains(polygons[i]): 

101 contains.add_edge(j, i) 

102 

103 # a root or exterior curve has an even number of parents 

104 # wrap in dict call to avoid networkx view 

105 degree = dict(contains.in_degree()) 

106 # convert keys and values to numpy arrays 

107 indexes = np.array(list(degree.keys())) 

108 degrees = np.array(list(degree.values())) 

109 # roots are curves with an even inward degree (parent count) 

110 roots = indexes[(degrees % 2) == 0] 

111 # if there are multiple nested polygons split the graph 

112 # so the contains logic returns the individual polygons 

113 if len(degrees) > 0 and degrees.max() > 1: 

114 # collect new edges for graph 

115 edges = [] 

116 # order the roots so they are sorted by degree 

117 roots = roots[np.argsort([degree[r] for r in roots])] 

118 # find edges of subgraph for each root and children 

119 for root in roots: 

120 children = indexes[degrees == degree[root] + 1] 

121 edges.extend(contains.subgraph(np.append(children, root)).edges()) 

122 # stack edges into new directed graph 

123 contains = nx.from_edgelist(edges, nx.DiGraph()) 

124 # if roots have no children add them anyway 

125 contains.add_nodes_from(roots) 

126 

127 return roots, contains 

128 

129 

130def edges_to_polygons(edges: NDArray[int64], vertices: NDArray[float64]): 

131 """ 

132 Given an edge list of indices and associated vertices 

133 representing lines, generate a list of polygons. 

134 

135 Parameters 

136 ----------- 

137 edges : (n, 2) 

138 Indexes of vertices which represent lines 

139 vertices : (m, 2) 

140 Vertices in 2D space. 

141 

142 Returns 

143 ---------- 

144 polygons : (p,) shapely.geometry.Polygon 

145 Polygon objects with interiors 

146 """ 

147 

148 assert isinstance(vertices, np.ndarray) 

149 

150 # create closed polygon objects 

151 polygons = [] 

152 # loop through a sequence of ordered traversals 

153 for dfs in graph.traversals(edges, mode="dfs"): 

154 try: 

155 # try to recover polygons before they are more complicated 

156 repaired = repair_invalid(Polygon(vertices[dfs])) 

157 # if it returned a multipolygon extend into a flat list 

158 if hasattr(repaired, "geoms"): 

159 polygons.extend(repaired.geoms) 

160 else: 

161 polygons.append(repaired) 

162 except ValueError: 

163 continue 

164 

165 # if there is only one polygon, just return it 

166 if len(polygons) == 1: 

167 return polygons 

168 

169 # find which polygons contain which other polygons 

170 roots, tree = enclosure_tree(polygons) 

171 

172 # generate polygons with proper interiors 

173 return [ 

174 Polygon( 

175 shell=polygons[root].exterior, 

176 holes=[polygons[i].exterior for i in tree[root].keys()], 

177 ) 

178 for root in roots 

179 ] 

180 

181 

182def polygons_obb(polygons: Iterable[Polygon] | ArrayLike): 

183 """ 

184 Find the OBBs for a list of shapely.geometry.Polygons 

185 """ 

186 rectangles = [None] * len(polygons) 

187 transforms = [None] * len(polygons) 

188 for i, p in enumerate(polygons): 

189 transforms[i], rectangles[i] = polygon_obb(p) 

190 return np.array(transforms), np.array(rectangles) 

191 

192 

193def polygon_obb(polygon: Polygon | NDArray): 

194 """ 

195 Find the oriented bounding box of a Shapely polygon. 

196 

197 The OBB is always aligned with an edge of the convex hull of the polygon. 

198 

199 Parameters 

200 ------------- 

201 polygons : shapely.geometry.Polygon 

202 Input geometry 

203 

204 Returns 

205 ------------- 

206 transform : (3, 3) float 

207 Transformation matrix 

208 which will move input polygon from its original position 

209 to the first quadrant where the AABB is the OBB 

210 extents : (2,) float 

211 Extents of transformed polygon 

212 """ 

213 if hasattr(polygon, "exterior"): 

214 points = np.asanyarray(polygon.exterior.coords) 

215 elif isinstance(polygon, np.ndarray): 

216 points = polygon 

217 else: 

218 raise ValueError("polygon or points must be provided") 

219 

220 transform, extents = bounds.oriented_bounds_2D(points) 

221 

222 if tol.strict: 

223 moved = transform_points(points=points, matrix=transform) 

224 assert np.allclose(-extents / 2.0, moved.min(axis=0)) 

225 assert np.allclose(extents / 2.0, moved.max(axis=0)) 

226 

227 return transform, extents 

228 

229 

230def transform_polygon(polygon, matrix): 

231 """ 

232 Transform a polygon by a a 2D homogeneous transform. 

233 

234 Parameters 

235 ------------- 

236 polygon : shapely.geometry.Polygon 

237 2D polygon to be transformed. 

238 matrix : (3, 3) float 

239 2D homogeneous transformation. 

240 

241 Returns 

242 -------------- 

243 result : shapely.geometry.Polygon 

244 Polygon transformed by matrix. 

245 """ 

246 matrix = np.asanyarray(matrix, dtype=np.float64) 

247 

248 if hasattr(polygon, "geoms"): 

249 result = [transform_polygon(p, t) for p, t in zip(polygon, matrix)] 

250 return result 

251 # transform the outer shell 

252 shell = transform_points(np.array(polygon.exterior.coords), matrix)[:, :2] 

253 # transform the interiors 

254 holes = [ 

255 transform_points(np.array(i.coords), matrix)[:, :2] for i in polygon.interiors 

256 ] 

257 # create a new polygon with the result 

258 result = Polygon(shell=shell, holes=holes) 

259 return result 

260 

261 

262def polygon_bounds(polygon, matrix=None): 

263 """ 

264 Get the transformed axis aligned bounding box of a 

265 shapely Polygon object. 

266 

267 Parameters 

268 ------------ 

269 polygon : shapely.geometry.Polygon 

270 Polygon pre-transform 

271 matrix : (3, 3) float or None. 

272 Homogeneous transform moving polygon in space 

273 

274 Returns 

275 ------------ 

276 bounds : (2, 2) float 

277 Axis aligned bounding box of transformed polygon. 

278 """ 

279 if matrix is not None: 

280 assert matrix.shape == (3, 3) 

281 points = transform_points(points=np.array(polygon.exterior.coords), matrix=matrix) 

282 else: 

283 points = np.array(polygon.exterior.coords) 

284 

285 bounds = np.array([points.min(axis=0), points.max(axis=0)]) 

286 assert bounds.shape == (2, 2) 

287 return bounds 

288 

289 

290def plot(polygon=None, show=True, axes=None, **kwargs): 

291 """ 

292 Plot a shapely polygon using matplotlib. 

293 

294 Parameters 

295 ------------ 

296 polygon : shapely.geometry.Polygon 

297 Polygon to be plotted 

298 show : bool 

299 If True will display immediately 

300 **kwargs 

301 Passed to plt.plot 

302 """ 

303 import matplotlib.pyplot as plt # noqa 

304 

305 def plot_single(single): 

306 axes.plot(*single.exterior.xy, **kwargs) 

307 for interior in single.interiors: 

308 axes.plot(*interior.xy, **kwargs) 

309 

310 # make aspect ratio non-stupid 

311 if axes is None: 

312 axes = plt.axes() 

313 axes.set_aspect("equal", "datalim") 

314 

315 if polygon.__class__.__name__ == "MultiPolygon": 

316 [plot_single(i) for i in polygon.geoms] 

317 elif hasattr(polygon, "__iter__"): 

318 [plot_single(i) for i in polygon] 

319 elif polygon is not None: 

320 plot_single(polygon) 

321 

322 if show: 

323 plt.show() 

324 

325 return axes 

326 

327 

328def resample_boundaries(polygon: Polygon, resolution: float, clip=None): 

329 """ 

330 Return a version of a polygon with boundaries re-sampled 

331 to a specified resolution. 

332 

333 Parameters 

334 ------------- 

335 polygon : shapely.geometry.Polygon 

336 Source geometry 

337 resolution : float 

338 Desired distance between points on boundary 

339 clip : (2,) int 

340 Upper and lower bounds to clip 

341 number of samples to avoid exploding count 

342 

343 Returns 

344 ------------ 

345 kwargs : dict 

346 Keyword args for a Polygon constructor `Polygon(**kwargs)` 

347 """ 

348 

349 def resample_boundary(boundary): 

350 # add a polygon.exterior or polygon.interior to 

351 # the deque after resampling based on our resolution 

352 count = boundary.length / resolution 

353 count = int(np.clip(count, *clip)) 

354 return resample_path(boundary.coords, count=count) 

355 

356 if clip is None: 

357 clip = [8, 200] 

358 # create a sequence of [(n,2)] points 

359 kwargs = {"shell": resample_boundary(polygon.exterior), "holes": []} 

360 for interior in polygon.interiors: 

361 kwargs["holes"].append(resample_boundary(interior)) 

362 

363 return kwargs 

364 

365 

366def stack_boundaries(boundaries): 

367 """ 

368 Stack the boundaries of a polygon into a single 

369 (n, 2) list of vertices. 

370 

371 Parameters 

372 ------------ 

373 boundaries : dict 

374 With keys 'shell', 'holes' 

375 

376 Returns 

377 ------------ 

378 stacked : (n, 2) float 

379 Stacked vertices 

380 """ 

381 if len(boundaries["holes"]) == 0: 

382 return boundaries["shell"] 

383 return np.vstack((boundaries["shell"], np.vstack(boundaries["holes"]))) 

384 

385 

386def medial_axis(polygon: Polygon, resolution: Number | None = None, clip=None): 

387 """ 

388 Given a shapely polygon, find the approximate medial axis 

389 using a voronoi diagram of evenly spaced points on the 

390 boundary of the polygon. 

391 

392 Parameters 

393 ---------- 

394 polygon : shapely.geometry.Polygon 

395 The source geometry 

396 resolution : float 

397 Distance between each sample on the polygon boundary 

398 clip : None, or (2,) int 

399 Clip sample count to min of clip[0] and max of clip[1] 

400 

401 Returns 

402 ---------- 

403 edges : (n, 2) int 

404 Vertex indices representing line segments 

405 on the polygon's medial axis 

406 vertices : (m, 2) float 

407 Vertex positions in space 

408 """ 

409 # a circle will have a single point medial axis 

410 if len(polygon.interiors) == 0: 

411 # what is the approximate scale of the polygon 

412 scale = np.ptp(np.reshape(polygon.bounds, (2, 2)), axis=0).max() 

413 # a (center, radius, error) tuple 

414 fit = fit_circle_check(polygon.exterior.coords, scale=scale) 

415 # is this polygon in fact a circle 

416 if fit is not None: 

417 # return an edge that has the center as the midpoint 

418 epsilon = np.clip(fit["radius"] / 500, 1e-5, np.inf) 

419 vertices = np.array( 

420 [fit["center"] + [0, epsilon], fit["center"] - [0, epsilon]], 

421 dtype=np.float64, 

422 ) 

423 # return a single edge to avoid consumers needing to special case 

424 edges = np.array([[0, 1]], dtype=np.int64) 

425 return edges, vertices 

426 

427 from scipy.spatial import Voronoi 

428 from shapely import vectorized 

429 

430 if resolution is None: 

431 resolution = np.ptp(np.reshape(polygon.bounds, (2, 2)), axis=0).max() / 100 

432 

433 # get evenly spaced points on the polygons boundaries 

434 samples = resample_boundaries(polygon=polygon, resolution=resolution, clip=clip) 

435 # stack the boundary into a (m,2) float array 

436 samples = stack_boundaries(samples) 

437 # create the voronoi diagram on 2D points 

438 voronoi = Voronoi(samples) 

439 # which voronoi vertices are contained inside the polygon 

440 contains = vectorized.contains(polygon, *voronoi.vertices.T) 

441 # ridge vertices of -1 are outside, make sure they are False 

442 contains = np.append(contains, False) 

443 # make sure ridge vertices is numpy array 

444 ridge = np.asanyarray(voronoi.ridge_vertices, dtype=np.int64) 

445 # only take ridges where every vertex is contained 

446 edges = ridge[contains[ridge].all(axis=1)] 

447 

448 # now we need to remove uncontained vertices 

449 contained = np.unique(edges) 

450 mask = np.zeros(len(voronoi.vertices), dtype=np.int64) 

451 mask[contained] = np.arange(len(contained)) 

452 

453 # mask voronoi vertices 

454 vertices = voronoi.vertices[contained] 

455 # re-index edges 

456 edges_final = mask[edges] 

457 

458 if tol.strict: 

459 # make sure we didn't screw up indexes 

460 assert np.ptp(vertices[edges_final] - voronoi.vertices[edges]) < 1e-5 

461 

462 return edges_final, vertices 

463 

464 

465def identifier(polygon: Polygon) -> NDArray[float64]: 

466 """ 

467 Return a vector containing values representative of 

468 a particular polygon. 

469 

470 Parameters 

471 --------- 

472 polygon : shapely.geometry.Polygon 

473 Input geometry 

474 

475 Returns 

476 --------- 

477 identifier : (8,) float 

478 Values which should be unique for this polygon. 

479 """ 

480 result = [ 

481 len(polygon.interiors), 

482 polygon.convex_hull.area, 

483 polygon.convex_hull.length, 

484 polygon.area, 

485 polygon.length, 

486 polygon.exterior.length, 

487 ] 

488 # include the principal second moments of inertia of the polygon 

489 # this is invariant to rotation and translation 

490 _, principal, _, _ = second_moments(polygon, return_centered=True) 

491 result.extend(principal) 

492 

493 return np.array(result, dtype=np.float64) 

494 

495 

496def random_polygon(segments=8, radius=1.0, seed: Seed = None): 

497 """ 

498 Generate a random polygon with a maximum number of sides and approximate radius. 

499 

500 Parameters 

501 --------- 

502 segments : int 

503 The maximum number of sides the random polygon will have 

504 radius : float 

505 The approximate radius of the polygon desired 

506 seed : None or int 

507 Seed for deterministic results, otherwise OS entropy. 

508 

509 Returns 

510 --------- 

511 polygon : shapely.geometry.Polygon 

512 Geometry object with random exterior and no interiors. 

513 """ 

514 random = random_generator(seed) 

515 angles = np.sort(np.cumsum(random.random(segments) * np.pi * 2) % (np.pi * 2)) 

516 radii = random.random(segments) * radius 

517 

518 points = np.column_stack((np.cos(angles), np.sin(angles))) * radii.reshape((-1, 1)) 

519 points = np.vstack((points, points[0])) 

520 polygon = Polygon(points).buffer(0.0) 

521 if hasattr(polygon, "geoms"): 

522 return polygon.geoms[0] 

523 return polygon 

524 

525 

526def polygon_scale(polygon): 

527 """ 

528 For a Polygon object return the diagonal length of the AABB. 

529 

530 Parameters 

531 ------------ 

532 polygon : shapely.geometry.Polygon 

533 Source geometry 

534 

535 Returns 

536 ------------ 

537 scale : float 

538 Length of AABB diagonal 

539 """ 

540 extents = np.ptp(np.reshape(polygon.bounds, (2, 2)), axis=0) 

541 scale = (extents**2).sum() ** 0.5 

542 

543 return scale 

544 

545 

546def paths_to_polygons(paths, scale=None): 

547 """ 

548 Given a sequence of connected points turn them into 

549 valid shapely Polygon objects. 

550 

551 Parameters 

552 ----------- 

553 paths : (n,) sequence 

554 Of (m, 2) float closed paths 

555 scale : float 

556 Approximate scale of drawing for precision 

557 

558 Returns 

559 ----------- 

560 polys : (p,) list 

561 Filled with Polygon or None 

562 

563 """ 

564 polygons = [None] * len(paths) 

565 for i, path in enumerate(paths): 

566 if len(path) < 4: 

567 # since the first and last vertices are identical in 

568 # a closed loop a 4 vertex path is the minimum for 

569 # non-zero area 

570 continue 

571 try: 

572 polygon = Polygon(path) 

573 if polygon.is_valid: 

574 polygons[i] = polygon 

575 else: 

576 polygons[i] = repair_invalid(polygon, scale) 

577 except ValueError: 

578 # raised if a polygon is unrecoverable 

579 continue 

580 except BaseException: 

581 log.error("unrecoverable polygon", exc_info=True) 

582 polygons = np.array(polygons) 

583 

584 return polygons 

585 

586 

587def sample(polygon, count, factor=1.5, max_iter=10, seed: Seed = None): 

588 """ 

589 Use rejection sampling to generate random points inside a 

590 polygon. Note that this function may return fewer or no 

591 points, in particular if the polygon as very little area 

592 compared to the area of the axis-aligned bounding box. 

593 

594 Parameters 

595 ----------- 

596 polygon : shapely.geometry.Polygon 

597 Polygon that will contain points 

598 count : int 

599 Number of points to return 

600 factor : float 

601 How many points to test per loop 

602 max_iter : int 

603 Maximum number of intersection checks is: 

604 > count * factor * max_iter 

605 seed : None or int 

606 Seed for deterministic results, otherwise OS entropy. 

607 

608 Returns 

609 ----------- 

610 hit : (n, 2) float 

611 Random points inside polygon 

612 where n <= count 

613 """ 

614 # do batch point-in-polygon queries 

615 from shapely import vectorized 

616 

617 # TODO : this should probably have some option to 

618 # sample from the *oriented* bounding box which would 

619 # make certain cases much, much more efficient. 

620 

621 # get size of bounding box 

622 bounds = np.reshape(polygon.bounds, (2, 2)) 

623 extents = np.ptp(bounds, axis=0) 

624 

625 # how many points to check per loop iteration 

626 per_loop = int(count * factor) 

627 

628 # start with some rejection sampling 

629 random = random_generator(seed) 

630 points = bounds[0] + extents * random.random((per_loop, 2)) 

631 # do the point in polygon test and append resulting hits 

632 mask = vectorized.contains(polygon, *points.T) 

633 hit = [points[mask]] 

634 hit_count = len(hit[0]) 

635 # if our first non-looping check got enough samples exit 

636 if hit_count >= count: 

637 return hit[0][:count] 

638 

639 # if we have to do iterations loop here slowly 

640 for _ in range(max_iter): 

641 # generate points inside polygons AABB 

642 points = (random.random((per_loop, 2)) * extents) + bounds[0] 

643 # do the point in polygon test and append resulting hits 

644 mask = vectorized.contains(polygon, *points.T) 

645 hit.append(points[mask]) 

646 # keep track of how many points we've collected 

647 hit_count += len(hit[-1]) 

648 # if we have enough points exit the loop 

649 if hit_count > count: 

650 break 

651 

652 # stack the hits into an (n,2) array and truncate 

653 hit = np.vstack(hit)[:count] 

654 

655 return hit 

656 

657 

658def repair_invalid(polygon, scale=None, rtol=0.5): 

659 """ 

660 Given a shapely.geometry.Polygon, attempt to return a 

661 valid version of the polygon through buffering tricks. 

662 

663 Parameters 

664 ----------- 

665 polygon : shapely.geometry.Polygon 

666 Source geometry 

667 rtol : float 

668 How close does a perimeter have to be 

669 scale : float or None 

670 For numerical precision reference 

671 

672 Returns 

673 ---------- 

674 repaired : shapely.geometry.Polygon 

675 Repaired polygon 

676 

677 Raises 

678 ---------- 

679 ValueError 

680 If polygon can't be repaired 

681 """ 

682 if hasattr(polygon, "is_valid") and polygon.is_valid: 

683 return polygon 

684 

685 # basic repair involves buffering the polygon outwards 

686 # this will fix a subset of problems. 

687 basic = polygon.buffer(tol.zero) 

688 # if it returned multiple polygons check the largest 

689 if hasattr(basic, "geoms"): 

690 basic = basic.geoms[np.argmax([i.area for i in basic.geoms])] 

691 

692 # check perimeter of result against original perimeter 

693 if basic.is_valid and np.isclose(basic.length, polygon.length, rtol=rtol): 

694 return basic 

695 

696 if scale is None: 

697 distance = 0.002 * np.ptp(np.reshape(polygon.bounds, (2, 2)), axis=0).mean() 

698 else: 

699 distance = 0.002 * scale 

700 

701 # if there are no interiors, we can work with just the exterior 

702 # ring, which is often more reliable 

703 if len(polygon.interiors) == 0: 

704 # try buffering the exterior of the polygon 

705 # the interior will be offset by -tol.buffer 

706 rings = polygon.exterior.buffer(distance).interiors 

707 if len(rings) == 1: 

708 # reconstruct a single polygon from the interior ring 

709 recon = Polygon(shell=rings[0]).buffer(distance) 

710 # check perimeter of result against original perimeter 

711 if recon.is_valid and np.isclose(recon.length, polygon.length, rtol=rtol): 

712 return recon 

713 

714 # try de-deuplicating the outside ring 

715 points = np.array(polygon.exterior.coords) 

716 # remove any segments shorter than tol.merge 

717 # this is a little risky as if it was discretized more 

718 # finely than 1-e8 it may remove detail 

719 unique = np.append(True, (np.diff(points, axis=0) ** 2).sum(axis=1) ** 0.5 > 1e-8) 

720 # make a new polygon with result 

721 dedupe = Polygon(shell=points[unique]) 

722 # check result 

723 if dedupe.is_valid and np.isclose(dedupe.length, polygon.length, rtol=rtol): 

724 return dedupe 

725 

726 # buffer and unbuffer the whole polygon 

727 buffered = polygon.buffer(distance).buffer(-distance) 

728 # if it returned multiple polygons check the largest 

729 if hasattr(buffered, "geoms"): 

730 areas = np.array([b.area for b in buffered.geoms]) 

731 return buffered.geoms[areas.argmax()] 

732 

733 # check perimeter of result against original perimeter 

734 if buffered.is_valid and np.isclose(buffered.length, polygon.length, rtol=rtol): 

735 log.debug("Recovered invalid polygon through double buffering") 

736 return buffered 

737 

738 raise ValueError("unable to recover polygon!") 

739 

740 

741def projected( 

742 mesh, 

743 normal, 

744 origin=None, 

745 ignore_sign=True, 

746 rpad=1e-5, 

747 apad=None, 

748 tol_dot=1e-10, 

749 precise: bool = False, 

750 precise_eps: float = 1e-10, 

751): 

752 """ 

753 Project a mesh onto a plane and then extract the polygon 

754 that outlines the mesh projection on that plane. 

755 

756 Note that this will ignore back-faces, which is only 

757 relevant if the source mesh isn't watertight. 

758 

759 Also padding: this generates a result by unioning the 

760 polygons of multiple connected regions, which requires 

761 the polygons be padded by a distance so that a polygon 

762 union produces a single coherent result. This distance 

763 is calculated as: `apad + (rpad * scale)` 

764 

765 Parameters 

766 ---------- 

767 mesh : trimesh.Trimesh 

768 Source geometry 

769 normal : (3,) float 

770 Normal to extract flat pattern along 

771 origin : None or (3,) float 

772 Origin of plane to project mesh onto 

773 ignore_sign : bool 

774 Allow a projection from the normal vector in 

775 either direction: this provides a substantial speedup 

776 on watertight meshes where the direction is irrelevant 

777 but if you have a triangle soup and want to discard 

778 backfaces you should set this to False. 

779 rpad : float 

780 Proportion to pad polygons by before unioning 

781 and then de-padding result by to avoid zero-width gaps. 

782 apad : float 

783 Absolute padding to pad polygons by before unioning 

784 and then de-padding result by to avoid zero-width gaps. 

785 tol_dot : float 

786 Tolerance for discarding on-edge triangles. 

787 precise : bool 

788 Use the precise projection computation using shapely. 

789 precise_eps : float 

790 Tolerance for precise triangle checks. 

791 

792 Returns 

793 ---------- 

794 projected : shapely.geometry.Polygon or None 

795 Outline of source mesh 

796 

797 Raises 

798 --------- 

799 ValueError 

800 If max_regions is exceeded 

801 """ 

802 # make sure normal is a unitized copy 

803 normal = np.array(normal, dtype=np.float64) 

804 normal /= np.linalg.norm(normal) 

805 

806 # the projection of each face normal onto facet normal 

807 dot_face = np.dot(normal, mesh.face_normals.T) 

808 if ignore_sign: 

809 # for watertight mesh speed up projection by handling side with less faces 

810 # check if face lies on front or back of normal 

811 front = dot_face > tol_dot 

812 back = dot_face < -tol_dot 

813 # divide the mesh into front facing section and back facing parts 

814 # and discard the faces perpendicular to the axis. 

815 # since we are doing a unary_union later we can use the front *or* 

816 # the back so we use which ever one has fewer triangles 

817 # we want the largest nonzero group 

818 count = np.array([front.sum(), back.sum()]) 

819 if count.min() == 0: 

820 # if one of the sides has zero faces we need the other 

821 pick = count.argmax() 

822 else: 

823 # otherwise use the normal direction with the fewest faces 

824 pick = count.argmin() 

825 # use the picked side 

826 side = [front, back][pick] 

827 else: 

828 # if explicitly asked to care about the sign 

829 # only handle the front side of normal 

830 side = dot_face > tol_dot 

831 

832 # subset the adjacency pairs to ones which have both faces included 

833 # on the side we are currently looking at 

834 adjacency_check = side[mesh.face_adjacency].all(axis=1) 

835 adjacency = mesh.face_adjacency[adjacency_check] 

836 

837 # transform from the mesh frame in 3D to the XY plane 

838 to_2D = geometry.plane_transform(origin=origin, normal=normal) 

839 # transform mesh vertices to 2D and clip the zero Z 

840 vertices_2D = transform_points(mesh.vertices, to_2D)[:, :2] 

841 

842 if precise: 

843 # precise mode just unions triangles as one shapely 

844 # polygon per triangle which historically has been very slow 

845 # but it is more defensible intellectually 

846 faces = mesh.faces[side] 

847 # round the 2D vertices with slightly more precision 

848 # than our final dilate-erode cleanup 

849 digits = int(np.abs(np.log10(precise_eps)) + 2) 

850 rounded = np.round(vertices_2D, digits) 

851 # get the triangles as closed 4-vertex polygons 

852 triangles = rounded[np.column_stack((faces, faces[:, :1]))] 

853 # do a check for exactly-degenerate triangles where any two 

854 # vertices are exactly identical which means the triangle has 

855 # zero area 

856 valid = ~(triangles[:, [0, 0, 2]] == triangles[:, [1, 2, 1]]).all(axis=2).any( 

857 axis=1 

858 ) 

859 # union the valid triangles and then dilate-erode to clean up 

860 # any holes or defects smaller than precise_eps 

861 return ( 

862 ops.unary_union([Polygon(f) for f in triangles[valid]]) 

863 .buffer(precise_eps) 

864 .buffer(-precise_eps) 

865 ) 

866 

867 # a sequence of face indexes that are connected 

868 face_groups = graph.connected_components(adjacency, nodes=np.nonzero(side)[0]) 

869 

870 # reshape edges into shape length of faces for indexing 

871 edges = mesh.edges_sorted.reshape((-1, 6)) 

872 

873 polygons = [] 

874 for faces in face_groups: 

875 # index edges by face then shape back to individual edges 

876 edge = edges[faces].reshape((-1, 2)) 

877 # edges that occur only once are on the boundary 

878 group = grouping.group_rows(edge, require_count=1) 

879 # turn each region into polygons 

880 polygons.extend(edges_to_polygons(edges=edge[group], vertices=vertices_2D)) 

881 

882 padding = 0.0 

883 if apad is not None: 

884 # set padding by absolute value 

885 padding += float(apad)