Coverage for trimesh/points.py: 79%

1"""

2points.py

3-------------

5Functions dealing with (n, d) points.

6"""

8import copy

9from hashlib import sha256

11import numpy as np

12from numpy import float64

14from . import caching, grouping, transformations, util

15from .constants import tol

16from .geometry import plane_transform

17from .inertia import points_inertia

18from .parent import Geometry3D

19from .typed import ArrayLike, NDArray, NDArray2D

20from .visual.color import VertexColor

23def point_plane_distance(points, plane_normal, plane_origin=None):

24 """

25 The minimum perpendicular distance of a point to a plane.

27 Parameters

28 -----------

29 points : (n, 3) float

30 Points in space

31 plane_normal : (3,) float

32 Unit normal vector

33 plane_origin : (3,) float

34 Plane origin in space

36 Returns

37 ------------

38 distances : (n,) float

39 Distance from point to plane

40 """

41 points = np.asanyarray(points, dtype=float64)

42 if plane_origin is None:

43 w = points

44 else:

45 w = points - plane_origin

46 distances = np.dot(plane_normal, w.T) / np.linalg.norm(plane_normal)

47 return distances

50def major_axis(points):

51 """

52 Returns an approximate vector representing the major

53 axis of the passed points.

55 Parameters

56 -------------

57 points : (n, dimension) float

58 Points in space

60 Returns

61 -------------

62 axis : (dimension,) float

63 Vector along approximate major axis

64 """

65 _U, S, V = np.linalg.svd(points)

66 axis = util.unitize(np.dot(S, V))

67 return axis

70def plane_fit(points):

71 """

72 Fit a plane to points using SVD.

74 Parameters

75 ---------

76 points : (n, 3) float or (p, n, 3,) float

77 3D points in space

78 Second option allows to simultaneously compute

79 p centroids and normals

81 Returns

82 ---------

83 C : (3,) float or (p, 3,) float

84 Point on the plane

85 N : (3,) float or (p, 3,) float

86 Unit normal vector of plane

87 """

88 # make sure input is numpy array

89 points = np.asanyarray(points, dtype=float64)

90 assert points.ndim == 2 or points.ndim == 3

91 # with only one point set, np.dot is faster

92 if points.ndim == 2:

93 # make the plane origin the mean of the points

94 C = points.mean(axis=0)

95 # points offset by the plane origin

96 x = points - C[None, :]

97 # create a (3, 3) matrix

98 M = np.dot(x.T, x)

99 else:

100 # make the plane origin the mean of the points

101 C = points.mean(axis=1)

102 # points offset by the plane origin

103 x = points - C[:, None, :]

104 # create a (p, 3, 3) matrix

105 M = x.swapaxes(1, 2) @ x

106 # run SVD

107 N = np.linalg.svd(M)[0][..., -1]

108 # return the centroid(s) and normal(s)

109 return C, N

110

111

112def radial_sort(points, origin, normal, start=None):

113 """

114 Sorts a set of points radially (by angle) around an

115 axis specified by origin and normal vector.

116

117 Parameters

118 --------------

119 points : (n, 3) float

120 Points in space

121 origin : (3,) float

122 Origin to sort around

123 normal : (3,) float

124 Vector to sort around

125 start : (3,) float

126 Vector to specify start position in counter-clockwise

127 order viewing in direction of normal, MUST not be

128 parallel with normal

129

130 Returns

131 --------------

132 ordered : (n, 3) float

133 Same as input points but reordered

134 """

135

136 # create two axis perpendicular to each other and

137 # the normal and project the points onto them

138 if start is None:

139 axis1, axis0, _ = util.generate_basis(normal)

140 # Preserve the starting angle for a positive Z normal.

141 axis0 = -axis0

142 else:

143 normal, start = util.unitize([normal, start])

144 if np.abs(1 - np.abs(np.dot(normal, start))) < tol.zero:

145 raise ValueError("start must not parallel with normal")

146 axis0 = np.cross(start, normal)

147 axis1 = np.cross(axis0, normal)

148 vectors = points - origin

149 # calculate the angles of the points on the axis

150 angles = np.arctan2(np.dot(vectors, axis0), np.dot(vectors, axis1))

151 # return the points sorted by angle

152 return points[angles.argsort()[::-1]]

153

154

155def project_to_plane(

156 points,

157 plane_normal,

158 plane_origin,

159 transform=None,

160 return_transform=False,

161 return_planar=True,

162):

163 """

164 Project (n, 3) points onto a plane.

165

166 Parameters

167 -----------

168 points : (n, 3) float

169 Points in space.

170 plane_normal : (3,) float

171 Unit normal vector of plane

172 plane_origin : (3,)

173 Origin point of plane

174 transform : None or (4, 4) float

175 Homogeneous transform, if specified, normal+origin are overridden

176 return_transform : bool

177 Returns the (4, 4) matrix used or not

178 return_planar : bool

179 Return (n, 2) points rather than (n, 3) points

180 """

181

182 if np.all(np.abs(plane_normal) < tol.zero):

183 raise NameError("Normal must be nonzero!")

184

185 if transform is None:

186 transform = plane_transform(plane_origin, plane_normal)

187

188 transformed = transformations.transform_points(points, transform)

189 transformed = transformed[:, 0 : (3 - int(return_planar))]

190

191 if return_transform:

192 polygon_to_3D = np.linalg.inv(transform)

193 return transformed, polygon_to_3D

194 return transformed

195

196

197def remove_close(points, radius):

198 """

199 Given an (n, m) array of points return a subset of

200 points where no point is closer than radius.

201

202 Parameters

203 ------------

204 points : (n, dimension) float

205 Points in space

206 radius : float

207 Minimum radius between result points

208

209 Returns

210 ------------

211 culled : (m, dimension) float

212 Points in space

213 mask : (n,) bool

214 Which points from the original points were returned

215 """

216 from scipy.spatial import cKDTree

217

218 tree = cKDTree(points)

219 # get the index of every pair of points closer than our radius

220 pairs = tree.query_pairs(radius, output_type="ndarray")

221

222 # how often each vertex index appears in a pair

223 # this is essentially a cheaply computed "vertex degree"

224 # in the graph that we could construct for connected points

225 count = np.bincount(pairs.ravel(), minlength=len(points))

226

227 # for every pair we know we have to remove one of them

228 # which of the two options we pick can have a large impact

229 # on how much over-culling we end up doing

230 column = count[pairs].argmax(axis=1)

231

232 # take the value in each row with the highest degree

233 # there is probably better numpy slicing you could do here

234 highest = pairs.ravel()[column + 2 * np.arange(len(column))]

235

236 # mask the vertices by index

237 mask = np.ones(len(points), dtype=bool)

238 mask[highest] = False

239

240 if tol.strict:

241 # verify we actually did what we said we'd do

242 test = cKDTree(points[mask])

243 assert len(test.query_pairs(radius)) == 0

244

245 return points[mask], mask

246

247

248def k_means(points, k, **kwargs):

249 """

250 Find k centroids that attempt to minimize the k- means problem:

251 https://en.wikipedia.org/wiki/Metric_k-center

252

253 Parameters

254 ----------

255 points: (n, d) float

256 Points in space

257 k : int

258 Number of centroids to compute

259 **kwargs : dict

260 Passed directly to scipy.cluster.vq.kmeans

261

262 Returns

263 ----------

264 centroids : (k, d) float

265 Points in some space

266 labels: (n) int

267 Indexes for which points belong to which centroid

268 """

269 from scipy.cluster.vq import kmeans

270 from scipy.spatial import cKDTree

271

272 points = np.asanyarray(points, dtype=float64)

273 points_std = points.std(axis=0)

274 points_std[points_std < tol.zero] = 1

275 whitened = points / points_std

276 centroids_whitened, _distortion = kmeans(whitened, k, **kwargs)

277 centroids = centroids_whitened * points_std

278

279 # find which centroid each point is closest to

280 tree = cKDTree(centroids)

281 labels = tree.query(points, k=1)[1]

282

283 return centroids, labels

284

285

286def tsp(points, start=0):

287 """

288 Find an ordering of points where each is visited and

289 the next point is the closest in euclidean distance,

290 and if there are multiple points with equal distance

291 go to an arbitrary one.

292

293 Assumes every point is visitable from every other point,

294 i.e. the travelling salesman problem on a fully connected

295 graph. It is not a MINIMUM traversal; rather it is a

296 "not totally goofy traversal, quickly." On random points

297 this traversal is often ~20x shorter than random ordering,

298 and executes on 1000 points in around 29ms on a 2014 i7.

299

300 Parameters

301 ---------------

302 points : (n, dimension) float

303 ND points in space

304 start : int

305 The index of points we should start at

306

307 Returns

308 ---------------

309 traversal : (n,) int

310 Ordered traversal visiting every point

311 distances : (n - 1,) float

312 The euclidean distance between points in traversal

313 """

314 # points should be float

315 points = np.asanyarray(points, dtype=float64)

316

317 if len(points.shape) != 2:

318 raise ValueError("points must be (n, dimension)!")

319

320 # start should be an index

321 start = int(start)

322

323 # a mask of unvisited points by index

324 unvisited = np.ones(len(points), dtype=bool)

325 unvisited[start] = False

326

327 # traversal of points by index

328 traversal = np.zeros(len(points), dtype=np.int64) - 1

329 traversal[0] = start

330 # list of distances

331 distances = np.zeros(len(points) - 1, dtype=float64)

332 # a mask of indexes in order

333 index_mask = np.arange(len(points), dtype=np.int64)

334

335 # in the loop we want to call distances.sum(axis=1)

336 # a lot and it's actually kind of slow for "reasons"

337 # dot products with ones is equivalent and ~2x faster

338 sum_ones = np.ones(points.shape[1])

339

340 # loop through all points

341 for i in range(len(points) - 1):

342 # which point are we currently on

343 current = points[traversal[i]]

344

345 # do NlogN distance query

346 # use dot instead of .sum(axis=1) or np.linalg.norm

347 # as it is faster, also don't square root here

348 dist = np.dot((points[unvisited] - current) ** 2, sum_ones)

349

350 # minimum distance index

351 min_index = dist.argmin()

352 # successor is closest unvisited point

353 successor = index_mask[unvisited][min_index]

354 # update the mask

355 unvisited[successor] = False

356 # store the index to the traversal

357 traversal[i + 1] = successor

358 # store the distance

359 distances[i] = dist[min_index]

360

361 # we were comparing distance^2 so take square root

362 distances **= 0.5

363

364 return traversal, distances

365

366

367def plot_points(points, show=True):

368 """

369 Plot an (n, 3) list of points using matplotlib

370

371 Parameters

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

373 points : (n, 3) float

374 Points in space

375 show : bool

376 If False, will not show until plt.show() is called

377 """

378 # TODO : should this just use SceneViewer?

379 import matplotlib.pyplot as plt # noqa

380 from mpl_toolkits.mplot3d import Axes3D # noqa

381

382 points = np.asanyarray(points, dtype=float64)

383

384 if len(points.shape) != 2:

385 raise ValueError("Points must be (n, 2|3)!")

386

387 if points.shape[1] == 3:

388 fig = plt.figure()

389 ax = fig.add_subplot(111, projection="3d")

390 ax.scatter(*points.T)

391 elif points.shape[1] == 2:

392 plt.scatter(*points.T)

393 else:

394 raise ValueError(f"points not 2D/3D: {points.shape}")

395

396 if show:

397 plt.show()

398

399

400class PointCloud(Geometry3D):

401 """

402 Hold 3D points in an object which can be visualized

403 in a scene.

404 """

405

406 def __init__(self, vertices, colors=None, metadata=None, **kwargs):

407 """

408 Load an array of points into a PointCloud object.

409

410 Parameters

411 -------------

412 vertices : (n, 3) float

413 Points in space

414 colors : (n, 4) uint8 or None

415 RGBA colors for each point

416 metadata : dict or None

417 Metadata about points

418 """

419 self._data = caching.DataStore()

420 self._cache = caching.Cache(self._data.__hash__)

421 self.metadata = {}

422 if metadata is not None:

423 self.metadata.update(metadata)

424

425 # load vertices

426 self.vertices = vertices

427

428 if "vertex_colors" in kwargs and colors is None:

429 colors = kwargs["vertex_colors"]

430

431 # save visual data to vertex color object

432 self.visual = VertexColor(colors=colors, obj=self)

433

434 def __setitem__(self, *args, **kwargs):

435 return self.vertices.__setitem__(*args, **kwargs)

436

437 def __getitem__(self, *args, **kwargs):

438 return self.vertices.__getitem__(*args, **kwargs)

439

440 @property

441 def shape(self):

442 """

443 Get the shape of the pointcloud

444

445 Returns

446 ----------

447 shape : (2,) int

448 Shape of vertex array

449 """

450 return self.vertices.shape

451

452 @property

453 def is_empty(self):

454 """

455 Are there any vertices defined or not.

456

457 Returns

458 ----------

459 empty : bool

460 True if no vertices defined

461 """

462 return len(self.vertices) == 0

463

464 def copy(self):

465 """

466 Safely get a copy of the current point cloud.

467

468 Copied objects will have emptied caches to avoid memory

469 issues and so may be slow on initial operations until

470 caches are regenerated.

471

472 Current object will *not* have its cache cleared.

473

474 Returns

475 ---------

476 copied : trimesh.PointCloud

477 Copy of current point cloud

478 """

479 copied = PointCloud(vertices=None)

480

481 # copy vertex and face data

482 copied._data.data = copy.deepcopy(self._data.data)

483

484 # copy visual data

485 copied.visual = copy.deepcopy(self.visual)

486

487 # get metadata

488 copied.metadata = copy.deepcopy(self.metadata)

489

490 # make sure cache is set from here

491 copied._cache.clear()

492

493 return copied

494

495 def hash(self):

496 """

497 Get a hash of the current vertices.

498

499 Returns

500 ----------

501 hash : str

502 Hash of self.vertices

503 """

504 return self._data.__hash__()

505

506 @property

507 def identifier(self) -> NDArray[float64]:

508 """

509 Return a simple array representing this PointCloud

510 that can be used to identify identical arrays.

511

512 Returns

513 ----------

514 identifier : (9,)

515 A flat array of data representing the cloud.

516 """

517 return self.moment_inertia.ravel()

518

519 @property

520 def identifier_hash(self) -> str:

521 """

522 A hash of the PointCloud's identifier that can be used

523 to detect duplicates.

524 """

525 return sha256(

526 (self.identifier * 1e5).round().astype(np.int64).tobytes()

527 ).hexdigest()

528

529 def merge_vertices(self):

530 """

531 Merge vertices closer than tol.merge (default: 1e-8)

532 """

533 # run unique rows

534 unique, inverse = grouping.unique_rows(self.vertices)

535

536 # apply unique mask to vertices

537 self.vertices = self.vertices[unique]

538

539 # apply unique mask to colors

540 if self.colors is not None and len(self.colors) == len(inverse):

541 self.colors = self.colors[unique]

542

543 def apply_transform(self, transform):

544 """

545 Apply a homogeneous transformation to the PointCloud

546 object in- place.

547

548 Parameters

549 --------------

550 transform : (4, 4) float

551 Homogeneous transformation to apply to PointCloud

552 """

553 self.vertices = transformations.transform_points(self.vertices, matrix=transform)

554 return self

555

556 @property

557 def bounds(self):

558 """

559 The axis aligned bounds of the PointCloud

560

561 Returns

562 ------------

563 bounds : (2, 3) float

564 Minimum, Maximum verteex

565 """

566 return np.array([self.vertices.min(axis=0), self.vertices.max(axis=0)])

567

568 @property

569 def extents(self):

570 """

571 The size of the axis aligned bounds

572

573 Returns

574 ------------

575 extents : (3,) float

576 Edge length of axis aligned bounding box

577 """

578 return np.ptp(self.bounds, axis=0)

579

580 @property

581 def centroid(self):

582 """

583 The mean vertex position

584

585 Returns

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

587 centroid : (3,) float

588 Mean vertex position

589 """

590 return self.vertices.mean(axis=0)

591

592 @caching.cache_decorator

593 def moment_inertia(self) -> NDArray[float64]:

594 return points_inertia(points=self.vertices, weights=self.weights)

595

596 @property

597 def weights(self) -> NDArray[float64]:

598 """

599 If each point has a specific weight assigned to it.

600

601 Returns

602 -----------

603 weights : (n,)

604 A per-vertex weight.

605 """

606 current = self._data.get("weights")

607 if current is None:

608 ones = np.ones(len(self.vertices), dtype=np.float64)

609 self._data["weights"] = ones

610 return ones

611

612 return current

613

614 @weights.setter

615 def weights(self, values: ArrayLike):

616 """

617 Assign a weight to each point for later computation.

618

619 Parameters

620 -----------

621 values : (n,)

622 Weights for each vertex.

623 """

624 values = np.asanyarray(values, dtype=np.float64)

625 if values.shape != (self.shape[0],):

626 raise ValueError("Weights must match vertices!")

627 self._data["weights"] = values

628

629 @property

630 def vertices(self):

631 """

632 Vertices of the PointCloud

633

634 Returns

635 ------------

636 vertices : (n, 3) float

637 Points in the PointCloud

638 """

639 return self._data.get("vertices", np.zeros(shape=(0, 3), dtype=float64))

640

641 @vertices.setter

642 def vertices(self, values):

643 """

644 Assign vertex values to the point cloud.

645

646 Parameters

647 --------------

648 values : (n, 3) float

649 Points in space

650 """

651 if values is None or len(values) == 0:

652 return self._data.data.pop("vertices", None)

653 self._data["vertices"] = np.asanyarray(values, order="C", dtype=float64)

654

655 @property

656 def colors(self):

657 """

658 Stored per- point color

659

660 Returns

661 ----------

662 colors : (len(self.vertices), 4) np.uint8

663 Per- point RGBA color

664 """

665 return self.visual.vertex_colors

666

667 @colors.setter

668 def colors(self, data):

669 self.visual.vertex_colors = data

670

671 @caching.cache_decorator

672 def kdtree(self):

673 """

674 Return a scipy.spatial.cKDTree of the vertices of the mesh.

675 Not cached as this lead to observed memory issues and segfaults.

676

677 Returns

678 ---------

679 tree : scipy.spatial.cKDTree

680 Contains mesh.vertices

681 """

682

683 from scipy.spatial import cKDTree

684

685 tree = cKDTree(self.vertices.view(np.ndarray))

686 return tree

687

688 @caching.cache_decorator

689 def convex_hull(self):

690 """

691 A convex hull of every point.

692

693 Returns

694 -------------

695 convex_hull : trimesh.Trimesh

696 A watertight mesh of the hull of the points

697 """

698 from . import convex

699

700 return convex.convex_hull(self.vertices)

701

702 def scene(self):

703 """

704 A scene containing just the PointCloud

705

706 Returns

707 ----------

708 scene : trimesh.Scene

709 Scene object containing this PointCloud

710 """

711 from .scene.scene import Scene

712

713 return Scene(self)

714

715 def show(self, **kwargs):

716 """

717 Open a viewer window displaying the current PointCloud

718 """

719 self.scene().show(**kwargs)

720

721 def export(self, file_obj=None, file_type=None, **kwargs):

722 """

723 Export the current pointcloud to a file object.

724 If file_obj is a filename, file will be written there.

725 Supported formats are xyz

726 Parameters

727 ------------

728 file_obj: open writeable file object

729 str, file name where to save the pointcloud

730 None, if you would like this function to return the export blob

731 file_type: str

732 Which file type to export as.

733 If file name is passed this is not required

734 """

735 from .exchange.export import export_mesh

736

737 return export_mesh(self, file_obj=file_obj, file_type=file_type, **kwargs)

738

739 def query(self, input_points, **kwargs):

740 """

741 Find the the closest points and associated attributes from this PointCloud.

742 Parameters

743 ------------

744 input_points : (n, 3) float

745 Input query points

746 kwargs : dict

747 Arguments for proximity.query_from_points

748 result : proximity.NearestQueryResult

749 Result of the query.

750 """

751 from .proximity import query_from_points

752

753 return query_from_points(self.vertices, input_points, self.kdtree, **kwargs)

754

755 def __add__(self, other):

756 if len(other.colors) == len(self.colors) == 0:

757 colors = None

758 else:

759 # preserve colors

760 # if one point cloud has no color property use black

761 other_colors = (

762 [[0, 0, 0, 255]] * len(other.vertices)

763 if len(other.colors) == 0

764 else other.colors

765 )

766 self_colors = (

767 [[0, 0, 0, 255]] * len(self.vertices)

768 if len(self.colors) == 0

769 else self.colors

770 )

771 colors = np.vstack((self_colors, other_colors))

772 return PointCloud(

773 vertices=np.vstack((self.vertices, other.vertices)), colors=colors

774 )

775

776 def estimate_normals(self, cluster_size=30) -> NDArray2D[np.float64]:

777 """

778 Estimate normals of PointCloud from nearest neighbor points

779

780 Parameters

781 ------------

782 cluster_size : float

783 Number of nearest neighbors for each point to compute normals from

784 """

785 _, indices = self.kdtree.query(self.vertices, k=cluster_size)

786 points = self.vertices[indices]

787 centroids = np.mean(points, axis=1, keepdims=True)

788 C = points - centroids

789 covariances = np.matmul(C.transpose(0, 2, 1), C) / (cluster_size - 1)

790 _, v = np.linalg.eigh(covariances)

791 return v[:, :, 0]