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, float64, int64 

11from .simplify import fit_circle_check 

12from .traversal import resample_path 

13 

14try: 

15 import networkx as nx 

16except BaseException as E: 

17 # create a dummy module which will raise the ImportError 

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

19 from ..exceptions import ExceptionWrapper 

20 

21 nx = ExceptionWrapper(E) 

22try: 

23 from rtree.index import Index 

24except BaseException as E: 

25 # create a dummy module which will raise the ImportError 

26 from ..exceptions import ExceptionWrapper 

27 

28 Index = ExceptionWrapper(E) 

29 

30 

31def enclosure_tree(polygons): 

32 """ 

33 Given a list of shapely polygons with only exteriors, 

34 find which curves represent the exterior shell or root curve 

35 and which represent holes which penetrate the exterior. 

36 

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

38 and then exact polygon queries for a final result. 

39 

40 Parameters 

41 ----------- 

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

43 Polygons which only have exteriors and may overlap 

44 

45 Returns 

46 ----------- 

47 roots : (m,) int 

48 Index of polygons which are root 

49 contains : networkx.DiGraph 

50 Edges indicate a polygon is 

51 contained by another polygon 

52 """ 

53 

54 # nodes are indexes in polygons 

55 contains = nx.DiGraph() 

56 

57 if len(polygons) == 0: 

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

59 elif len(polygons) == 1: 

60 # add an early exit for only a single polygon 

61 contains.add_node(0) 

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

63 

64 # get the bounds for every valid polygon 

65 bounds = { 

66 k: v 

67 for k, v in { 

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

69 }.items() 

70 if len(v) == 4 

71 } 

72 

73 # make sure we don't have orphaned polygon 

74 contains.add_nodes_from(bounds.keys()) 

75 

76 if len(bounds) > 0: 

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

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

79 # we could return here but having multiple return paths 

80 # seems more dangerous than iterating through an empty graph 

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

82 

83 # loop through every polygon 

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

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

86 for j in tree.intersection(b): 

87 # if we are checking a polygon against itself continue 

88 if i == j: 

89 continue 

90 # do a more accurate polygon in polygon test 

91 # for the enclosure tree information 

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

93 contains.add_edge(i, j) 

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

95 contains.add_edge(j, i) 

96 

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

98 # wrap in dict call to avoid networkx view 

99 degree = dict(contains.in_degree()) 

100 # convert keys and values to numpy arrays 

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

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

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

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

105 # if there are multiple nested polygons split the graph 

106 # so the contains logic returns the individual polygons 

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

108 # collect new edges for graph 

109 edges = [] 

110 # order the roots so they are sorted by degree 

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

112 # find edges of subgraph for each root and children 

113 for root in roots: 

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

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

116 # stack edges into new directed graph 

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

118 # if roots have no children add them anyway 

119 contains.add_nodes_from(roots) 

120 

121 return roots, contains 

122 

123 

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

125 """ 

126 Given an edge list of indices and associated vertices 

127 representing lines, generate a list of polygons. 

128 

129 Parameters 

130 ----------- 

131 edges : (n, 2) 

132 Indexes of vertices which represent lines 

133 vertices : (m, 2) 

134 Vertices in 2D space. 

135 

136 Returns 

137 ---------- 

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

139 Polygon objects with interiors 

140 """ 

141 

142 assert isinstance(vertices, np.ndarray) 

143 

144 # create closed polygon objects 

145 polygons = [] 

146 # loop through a sequence of ordered traversals 

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

148 try: 

149 # try to recover polygons before they are more complicated 

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

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

152 if hasattr(repaired, "geoms"): 

153 polygons.extend(repaired.geoms) 

154 else: 

155 polygons.append(repaired) 

156 except ValueError: 

157 continue 

158 

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

160 if len(polygons) == 1: 

161 return polygons 

162 

163 # find which polygons contain which other polygons 

164 roots, tree = enclosure_tree(polygons) 

165 

166 # generate polygons with proper interiors 

167 return [ 

168 Polygon( 

169 shell=polygons[root].exterior, 

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

171 ) 

172 for root in roots 

173 ] 

174 

175 

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

177 """ 

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

179 """ 

180 rectangles = [None] * len(polygons) 

181 transforms = [None] * len(polygons) 

182 for i, p in enumerate(polygons): 

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

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

185 

186 

187def polygon_obb(polygon: Polygon | NDArray): 

188 """ 

189 Find the oriented bounding box of a Shapely polygon. 

190 

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

192 

193 Parameters 

194 ------------- 

195 polygons : shapely.geometry.Polygon 

196 Input geometry 

197 

198 Returns 

199 ------------- 

200 transform : (3, 3) float 

201 Transformation matrix 

202 which will move input polygon from its original position 

203 to the first quadrant where the AABB is the OBB 

204 extents : (2,) float 

205 Extents of transformed polygon 

206 """ 

207 if hasattr(polygon, "exterior"): 

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

209 elif isinstance(polygon, np.ndarray): 

210 points = polygon 

211 else: 

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

213 

214 transform, extents = bounds.oriented_bounds_2D(points) 

215 

216 if tol.strict: 

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

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

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

220 

221 return transform, extents 

222 

223 

224def transform_polygon(polygon, matrix): 

225 """ 

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

227 

228 Parameters 

229 ------------- 

230 polygon : shapely.geometry.Polygon 

231 2D polygon to be transformed. 

232 matrix : (3, 3) float 

233 2D homogeneous transformation. 

234 

235 Returns 

236 -------------- 

237 result : shapely.geometry.Polygon 

238 Polygon transformed by matrix. 

239 """ 

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

241 

242 if hasattr(polygon, "geoms"): 

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

244 return result 

245 # transform the outer shell 

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

247 # transform the interiors 

248 holes = [ 

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

250 ] 

251 # create a new polygon with the result 

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

253 return result 

254 

255 

256def polygon_bounds(polygon, matrix=None): 

257 """ 

258 Get the transformed axis aligned bounding box of a 

259 shapely Polygon object. 

260 

261 Parameters 

262 ------------ 

263 polygon : shapely.geometry.Polygon 

264 Polygon pre-transform 

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

266 Homogeneous transform moving polygon in space 

267 

268 Returns 

269 ------------ 

270 bounds : (2, 2) float 

271 Axis aligned bounding box of transformed polygon. 

272 """ 

273 if matrix is not None: 

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

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

276 else: 

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

278 

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

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

281 return bounds 

282 

283 

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

285 """ 

286 Plot a shapely polygon using matplotlib. 

287 

288 Parameters 

289 ------------ 

290 polygon : shapely.geometry.Polygon 

291 Polygon to be plotted 

292 show : bool 

293 If True will display immediately 

294 **kwargs 

295 Passed to plt.plot 

296 """ 

297 import matplotlib.pyplot as plt # noqa 

298 

299 def plot_single(single): 

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

301 for interior in single.interiors: 

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

303 

304 # make aspect ratio non-stupid 

305 if axes is None: 

306 axes = plt.axes() 

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

308 

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

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

311 elif hasattr(polygon, "__iter__"): 

312 [plot_single(i) for i in polygon] 

313 elif polygon is not None: 

314 plot_single(polygon) 

315 

316 if show: 

317 plt.show() 

318 

319 return axes 

320 

321 

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

323 """ 

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

325 to a specified resolution. 

326 

327 Parameters 

328 ------------- 

329 polygon : shapely.geometry.Polygon 

330 Source geometry 

331 resolution : float 

332 Desired distance between points on boundary 

333 clip : (2,) int 

334 Upper and lower bounds to clip 

335 number of samples to avoid exploding count 

336 

337 Returns 

338 ------------ 

339 kwargs : dict 

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

341 """ 

342 

343 def resample_boundary(boundary): 

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

345 # the deque after resampling based on our resolution 

346 count = boundary.length / resolution 

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

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

349 

350 if clip is None: 

351 clip = [8, 200] 

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

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

354 for interior in polygon.interiors: 

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

356 

357 return kwargs 

358 

359 

360def stack_boundaries(boundaries): 

361 """ 

362 Stack the boundaries of a polygon into a single 

363 (n, 2) list of vertices. 

364 

365 Parameters 

366 ------------ 

367 boundaries : dict 

368 With keys 'shell', 'holes' 

369 

370 Returns 

371 ------------ 

372 stacked : (n, 2) float 

373 Stacked vertices 

374 """ 

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

376 return boundaries["shell"] 

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

378 

379 

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

381 """ 

382 Given a shapely polygon, find the approximate medial axis 

383 using a voronoi diagram of evenly spaced points on the 

384 boundary of the polygon. 

385 

386 Parameters 

387 ---------- 

388 polygon : shapely.geometry.Polygon 

389 The source geometry 

390 resolution : float 

391 Distance between each sample on the polygon boundary 

392 clip : None, or (2,) int 

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

394 

395 Returns 

396 ---------- 

397 edges : (n, 2) int 

398 Vertex indices representing line segments 

399 on the polygon's medial axis 

400 vertices : (m, 2) float 

401 Vertex positions in space 

402 """ 

403 # a circle will have a single point medial axis 

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

405 # what is the approximate scale of the polygon 

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

407 # a (center, radius, error) tuple 

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

409 # is this polygon in fact a circle 

410 if fit is not None: 

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

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

413 vertices = np.array( 

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

415 dtype=np.float64, 

416 ) 

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

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

419 return edges, vertices 

420 

421 from scipy.spatial import Voronoi 

422 from shapely import vectorized 

423 

424 if resolution is None: 

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

426 

427 # get evenly spaced points on the polygons boundaries 

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

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

430 samples = stack_boundaries(samples) 

431 # create the voronoi diagram on 2D points 

432 voronoi = Voronoi(samples) 

433 # which voronoi vertices are contained inside the polygon 

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

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

436 contains = np.append(contains, False) 

437 # make sure ridge vertices is numpy array 

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

439 # only take ridges where every vertex is contained 

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

441 

442 # now we need to remove uncontained vertices 

443 contained = np.unique(edges) 

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

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

446 

447 # mask voronoi vertices 

448 vertices = voronoi.vertices[contained] 

449 # re-index edges 

450 edges_final = mask[edges] 

451 

452 if tol.strict: 

453 # make sure we didn't screw up indexes 

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

455 

456 return edges_final, vertices 

457 

458 

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

460 """ 

461 Return a vector containing values representative of 

462 a particular polygon. 

463 

464 Parameters 

465 --------- 

466 polygon : shapely.geometry.Polygon 

467 Input geometry 

468 

469 Returns 

470 --------- 

471 identifier : (8,) float 

472 Values which should be unique for this polygon. 

473 """ 

474 result = [ 

475 len(polygon.interiors), 

476 polygon.convex_hull.area, 

477 polygon.convex_hull.length, 

478 polygon.area, 

479 polygon.length, 

480 polygon.exterior.length, 

481 ] 

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

483 # this is invariant to rotation and translation 

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

485 result.extend(principal) 

486 

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

488 

489 

490def random_polygon(segments=8, radius=1.0): 

491 """ 

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

493 

494 Parameters 

495 --------- 

496 segments : int 

497 The maximum number of sides the random polygon will have 

498 radius : float 

499 The approximate radius of the polygon desired 

500 

501 Returns 

502 --------- 

503 polygon : shapely.geometry.Polygon 

504 Geometry object with random exterior and no interiors. 

505 """ 

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

507 radii = np.random.random(segments) * radius 

508 

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

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

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

512 if hasattr(polygon, "geoms"): 

513 return polygon.geoms[0] 

514 return polygon 

515 

516 

517def polygon_scale(polygon): 

518 """ 

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

520 

521 Parameters 

522 ------------ 

523 polygon : shapely.geometry.Polygon 

524 Source geometry 

525 

526 Returns 

527 ------------ 

528 scale : float 

529 Length of AABB diagonal 

530 """ 

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

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

533 

534 return scale 

535 

536 

537def paths_to_polygons(paths, scale=None): 

538 """ 

539 Given a sequence of connected points turn them into 

540 valid shapely Polygon objects. 

541 

542 Parameters 

543 ----------- 

544 paths : (n,) sequence 

545 Of (m, 2) float closed paths 

546 scale : float 

547 Approximate scale of drawing for precision 

548 

549 Returns 

550 ----------- 

551 polys : (p,) list 

552 Filled with Polygon or None 

553 

554 """ 

555 polygons = [None] * len(paths) 

556 for i, path in enumerate(paths): 

557 if len(path) < 4: 

558 # since the first and last vertices are identical in 

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

560 # non-zero area 

561 continue 

562 try: 

563 polygon = Polygon(path) 

564 if polygon.is_valid: 

565 polygons[i] = polygon 

566 else: 

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

568 except ValueError: 

569 # raised if a polygon is unrecoverable 

570 continue 

571 except BaseException: 

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

573 polygons = np.array(polygons) 

574 

575 return polygons 

576 

577 

578def sample(polygon, count, factor=1.5, max_iter=10): 

579 """ 

580 Use rejection sampling to generate random points inside a 

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

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

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

584 

585 Parameters 

586 ----------- 

587 polygon : shapely.geometry.Polygon 

588 Polygon that will contain points 

589 count : int 

590 Number of points to return 

591 factor : float 

592 How many points to test per loop 

593 max_iter : int 

594 Maximum number of intersection checks is: 

595 > count * factor * max_iter 

596 

597 Returns 

598 ----------- 

599 hit : (n, 2) float 

600 Random points inside polygon 

601 where n <= count 

602 """ 

603 # do batch point-in-polygon queries 

604 from shapely import vectorized 

605 

606 # TODO : this should probably have some option to 

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

608 # make certain cases much, much more efficient. 

609 

610 # get size of bounding box 

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

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

613 

614 # how many points to check per loop iteration 

615 per_loop = int(count * factor) 

616 

617 # start with some rejection sampling 

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

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

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

621 hit = [points[mask]] 

622 hit_count = len(hit[0]) 

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

624 if hit_count >= count: 

625 return hit[0][:count] 

626 

627 # if we have to do iterations loop here slowly 

628 for _ in range(max_iter): 

629 # generate points inside polygons AABB 

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

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

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

633 hit.append(points[mask]) 

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

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

636 # if we have enough points exit the loop 

637 if hit_count > count: 

638 break 

639 

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

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

642 

643 return hit 

644 

645 

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

647 """ 

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

649 valid version of the polygon through buffering tricks. 

650 

651 Parameters 

652 ----------- 

653 polygon : shapely.geometry.Polygon 

654 Source geometry 

655 rtol : float 

656 How close does a perimeter have to be 

657 scale : float or None 

658 For numerical precision reference 

659 

660 Returns 

661 ---------- 

662 repaired : shapely.geometry.Polygon 

663 Repaired polygon 

664 

665 Raises 

666 ---------- 

667 ValueError 

668 If polygon can't be repaired 

669 """ 

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

671 return polygon 

672 

673 # basic repair involves buffering the polygon outwards 

674 # this will fix a subset of problems. 

675 basic = polygon.buffer(tol.zero) 

676 # if it returned multiple polygons check the largest 

677 if hasattr(basic, "geoms"): 

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

679 

680 # check perimeter of result against original perimeter 

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

682 return basic 

683 

684 if scale is None: 

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

686 else: 

687 distance = 0.002 * scale 

688 

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

690 # ring, which is often more reliable 

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

692 # try buffering the exterior of the polygon 

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

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

695 if len(rings) == 1: 

696 # reconstruct a single polygon from the interior ring 

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

698 # check perimeter of result against original perimeter 

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

700 return recon 

701 

702 # try de-deuplicating the outside ring 

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

704 # remove any segments shorter than tol.merge 

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

706 # finely than 1-e8 it may remove detail 

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

708 # make a new polygon with result 

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

710 # check result 

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

712 return dedupe 

713 

714 # buffer and unbuffer the whole polygon 

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

716 # if it returned multiple polygons check the largest 

717 if hasattr(buffered, "geoms"): 

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

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

720 

721 # check perimeter of result against original perimeter 

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

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

724 return buffered 

725 

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

727 

728 

729def projected( 

730 mesh, 

731 normal, 

732 origin=None, 

733 ignore_sign=True, 

734 rpad=1e-5, 

735 apad=None, 

736 tol_dot=1e-10, 

737 precise: bool = False, 

738 precise_eps: float = 1e-10, 

739): 

740 """ 

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

742 that outlines the mesh projection on that plane. 

743 

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

745 relevant if the source mesh isn't watertight. 

746 

747 Also padding: this generates a result by unioning the 

748 polygons of multiple connected regions, which requires 

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

750 union produces a single coherent result. This distance 

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

752 

753 Parameters 

754 ---------- 

755 mesh : trimesh.Trimesh 

756 Source geometry 

757 normal : (3,) float 

758 Normal to extract flat pattern along 

759 origin : None or (3,) float 

760 Origin of plane to project mesh onto 

761 ignore_sign : bool 

762 Allow a projection from the normal vector in 

763 either direction: this provides a substantial speedup 

764 on watertight meshes where the direction is irrelevant 

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

766 backfaces you should set this to False. 

767 rpad : float 

768 Proportion to pad polygons by before unioning 

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

770 apad : float 

771 Absolute padding to pad polygons by before unioning 

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

773 tol_dot : float 

774 Tolerance for discarding on-edge triangles. 

775 precise : bool 

776 Use the precise projection computation using shapely. 

777 precise_eps : float 

778 Tolerance for precise triangle checks. 

779 

780 Returns 

781 ---------- 

782 projected : shapely.geometry.Polygon or None 

783 Outline of source mesh 

784 

785 Raises 

786 --------- 

787 ValueError 

788 If max_regions is exceeded 

789 """ 

790 # make sure normal is a unitized copy 

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

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

793 

794 # the projection of each face normal onto facet normal 

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

796 if ignore_sign: 

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

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

799 front = dot_face > tol_dot 

800 back = dot_face < -tol_dot 

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

802 # and discard the faces perpendicular to the axis. 

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

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

805 # we want the largest nonzero group 

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

807 if count.min() == 0: 

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

809 pick = count.argmax() 

810 else: 

811 # otherwise use the normal direction with the fewest faces 

812 pick = count.argmin() 

813 # use the picked side 

814 side = [front, back][pick] 

815 else: 

816 # if explicitly asked to care about the sign 

817 # only handle the front side of normal 

818 side = dot_face > tol_dot 

819 

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

821 # on the side we are currently looking at 

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

823 adjacency = mesh.face_adjacency[adjacency_check] 

824 

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

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

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

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

829 

830 if precise: 

831 # precise mode just unions triangles as one shapely 

832 # polygon per triangle which historically has been very slow 

833 # but it is more defensible intellectually 

834 faces = mesh.faces[side] 

835 # round the 2D vertices with slightly more precision 

836 # than our final dilate-erode cleanup 

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

838 rounded = np.round(vertices_2D, digits) 

839 # get the triangles as closed 4-vertex polygons 

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

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

842 # vertices are exactly identical which means the triangle has 

843 # zero area 

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

845 axis=1 

846 ) 

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

848 # any holes or defects smaller than precise_eps 

849 return ( 

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

851 .buffer(precise_eps) 

852 .buffer(-precise_eps) 

853 ) 

854 

855 # a sequence of face indexes that are connected 

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

857 

858 # reshape edges into shape length of faces for indexing 

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

860 

861 polygons = [] 

862 for faces in face_groups: 

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

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

865 # edges that occur only once are on the boundary 

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

867 # turn each region into polygons 

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

869 

870 padding = 0.0 

871 if apad is not None: 

872 # set padding by absolute value 

873 padding += float(apad) 

874 if rpad is not None: 

875 # get the 2D scale as the longest side of the AABB 

876 scale = np.ptp(vertices_2D, axis=0).max() 

877 # apply the scale-relative padding 

878 padding += float(rpad) * scale 

879 

880 # if there is only one region we don't need to run a union 

881 elif len(polygons) == 1: 

882 return polygons[0] 

883 elif len(polygons) == 0: 

884 return None