Coverage for trimesh/points.py: 77%

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

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 axis0 = [normal[0], normal[2], -normal[1]]

140 axis1 = np.cross(normal, axis0)

141 else:

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

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

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

145 axis0 = np.cross(start, normal)

146 axis1 = np.cross(axis0, normal)

147 vectors = points - origin

148 # calculate the angles of the points on the axis

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

150 # return the points sorted by angle

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

152

153

154def project_to_plane(

155 points,

156 plane_normal,

157 plane_origin,

158 transform=None,

159 return_transform=False,

160 return_planar=True,

161):

162 """

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

164

165 Parameters

166 -----------

167 points : (n, 3) float

168 Points in space.

169 plane_normal : (3,) float

170 Unit normal vector of plane

171 plane_origin : (3,)

172 Origin point of plane

173 transform : None or (4, 4) float

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

175 return_transform : bool

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

177 return_planar : bool

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

179 """

180

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

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

183

184 if transform is None:

185 transform = plane_transform(plane_origin, plane_normal)

186

187 transformed = transformations.transform_points(points, transform)

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

189

190 if return_transform:

191 polygon_to_3D = np.linalg.inv(transform)

192 return transformed, polygon_to_3D

193 return transformed

194

195

196def remove_close(points, radius):

197 """

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

199 points where no point is closer than radius.

200

201 Parameters

202 ------------

203 points : (n, dimension) float

204 Points in space

205 radius : float

206 Minimum radius between result points

207

208 Returns

209 ------------

210 culled : (m, dimension) float

211 Points in space

212 mask : (n,) bool

213 Which points from the original points were returned

214 """

215 from scipy.spatial import cKDTree

216

217 tree = cKDTree(points)

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

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

220

221 # how often each vertex index appears in a pair

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

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

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

225

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

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

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

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

230

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

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

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

234

235 # mask the vertices by index

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

237 mask[highest] = False

238

239 if tol.strict:

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

241 test = cKDTree(points[mask])

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

243

244 return points[mask], mask

245

246

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

248 """

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

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

251

252 Parameters

253 ----------

254 points: (n, d) float

255 Points in space

256 k : int

257 Number of centroids to compute

258 **kwargs : dict

259 Passed directly to scipy.cluster.vq.kmeans

260

261 Returns

262 ----------

263 centroids : (k, d) float

264 Points in some space

265 labels: (n) int

266 Indexes for which points belong to which centroid

267 """

268 from scipy.cluster.vq import kmeans

269 from scipy.spatial import cKDTree

270

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

272 points_std = points.std(axis=0)

273 points_std[points_std < tol.zero] = 1

274 whitened = points / points_std

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

276 centroids = centroids_whitened * points_std

277

278 # find which centroid each point is closest to

279 tree = cKDTree(centroids)

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

281

282 return centroids, labels

283

284

285def tsp(points, start=0):

286 """

287 Find an ordering of points where each is visited and

288 the next point is the closest in euclidean distance,

289 and if there are multiple points with equal distance

290 go to an arbitrary one.

291

292 Assumes every point is visitable from every other point,

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

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

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

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

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

298

299 Parameters

300 ---------------

301 points : (n, dimension) float

302 ND points in space

303 start : int

304 The index of points we should start at

305

306 Returns

307 ---------------

308 traversal : (n,) int

309 Ordered traversal visiting every point

310 distances : (n - 1,) float

311 The euclidean distance between points in traversal

312 """

313 # points should be float

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

315

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

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

318

319 # start should be an index

320 start = int(start)

321

322 # a mask of unvisited points by index

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

324 unvisited[start] = False

325

326 # traversal of points by index

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

328 traversal[0] = start

329 # list of distances

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

331 # a mask of indexes in order

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

333

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

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

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

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

338

339 # loop through all points

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

341 # which point are we currently on

342 current = points[traversal[i]]

343

344 # do NlogN distance query

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

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

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

348

349 # minimum distance index

350 min_index = dist.argmin()

351 # successor is closest unvisited point

352 successor = index_mask[unvisited][min_index]

353 # update the mask

354 unvisited[successor] = False

355 # store the index to the traversal

356 traversal[i + 1] = successor

357 # store the distance

358 distances[i] = dist[min_index]

359

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

361 distances **= 0.5

362

363 return traversal, distances

364

365

366def plot_points(points, show=True):

367 """

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

369

370 Parameters

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

372 points : (n, 3) float

373 Points in space

374 show : bool

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

376 """

377 # TODO : should this just use SceneViewer?

378 import matplotlib.pyplot as plt # noqa

379 from mpl_toolkits.mplot3d import Axes3D # noqa

380

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

382

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

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

385

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

387 fig = plt.figure()

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

389 ax.scatter(*points.T)

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

391 plt.scatter(*points.T)

392 else:

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

394

395 if show:

396 plt.show()

397

398

399class PointCloud(Geometry3D):

400 """

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

402 in a scene.

403 """

404

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

406 """

407 Load an array of points into a PointCloud object.

408

409 Parameters

410 -------------

411 vertices : (n, 3) float

412 Points in space

413 colors : (n, 4) uint8 or None

414 RGBA colors for each point

415 metadata : dict or None

416 Metadata about points

417 """

418 self._data = caching.DataStore()

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

420 self.metadata = {}

421 if metadata is not None:

422 self.metadata.update(metadata)

423

424 # load vertices

425 self.vertices = vertices

426

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

428 colors = kwargs["vertex_colors"]

429

430 # save visual data to vertex color object

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

432

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

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

435

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

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

438

439 @property

440 def shape(self):

441 """

442 Get the shape of the pointcloud

443

444 Returns

445 ----------

446 shape : (2,) int

447 Shape of vertex array

448 """

449 return self.vertices.shape

450

451 @property

452 def is_empty(self):

453 """

454 Are there any vertices defined or not.

455

456 Returns

457 ----------

458 empty : bool

459 True if no vertices defined

460 """

461 return len(self.vertices) == 0

462

463 def copy(self):

464 """

465 Safely get a copy of the current point cloud.

466

467 Copied objects will have emptied caches to avoid memory

468 issues and so may be slow on initial operations until

469 caches are regenerated.

470

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

472

473 Returns

474 ---------

475 copied : trimesh.PointCloud

476 Copy of current point cloud

477 """

478 copied = PointCloud(vertices=None)

479

480 # copy vertex and face data

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

482

483 # copy visual data

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

485

486 # get metadata

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

488

489 # make sure cache is set from here

490 copied._cache.clear()

491

492 return copied

493

494 def hash(self):

495 """

496 Get a hash of the current vertices.

497

498 Returns

499 ----------

500 hash : str

501 Hash of self.vertices

502 """

503 return self._data.__hash__()

504

505 @property

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

507 """

508 Return a simple array representing this PointCloud

509 that can be used to identify identical arrays.

510

511 Returns

512 ----------

513 identifier : (9,)

514 A flat array of data representing the cloud.

515 """

516 return self.moment_inertia.ravel()

517

518 @property

519 def identifier_hash(self) -> str:

520 """

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

522 to detect duplicates.

523 """

524 return sha256(

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

526 ).hexdigest()

527

528 def merge_vertices(self):

529 """

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

531 """

532 # run unique rows

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

534

535 # apply unique mask to vertices

536 self.vertices = self.vertices[unique]

537

538 # apply unique mask to colors

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

540 self.colors = self.colors[unique]

541

542 def apply_transform(self, transform):

543 """

544 Apply a homogeneous transformation to the PointCloud

545 object in- place.

546

547 Parameters

548 --------------

549 transform : (4, 4) float

550 Homogeneous transformation to apply to PointCloud

551 """

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

553 return self

554

555 @property

556 def bounds(self):

557 """

558 The axis aligned bounds of the PointCloud

559

560 Returns

561 ------------

562 bounds : (2, 3) float

563 Minimum, Maximum verteex

564 """

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

566

567 @property

568 def extents(self):

569 """

570 The size of the axis aligned bounds

571

572 Returns

573 ------------

574 extents : (3,) float

575 Edge length of axis aligned bounding box

576 """

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

578

579 @property

580 def centroid(self):

581 """

582 The mean vertex position

583

584 Returns

585 ------------

586 centroid : (3,) float

587 Mean vertex position

588 """

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

590

591 @caching.cache_decorator

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

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

594

595 @property

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

597 """

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

599

600 Returns

601 -----------

602 weights : (n,)

603 A per-vertex weight.

604 """

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

606 if current is None:

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

608 self._data["weights"] = ones

609 return ones

610

611 return current

612

613 @weights.setter

614 def weights(self, values: ArrayLike):

615 """

616 Assign a weight to each point for later computation.

617

618 Parameters

619 -----------

620 values : (n,)

621 Weights for each vertex.

622 """

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

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

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

626 self._data["weights"] = values

627

628 @property

629 def vertices(self):

630 """

631 Vertices of the PointCloud

632

633 Returns

634 ------------

635 vertices : (n, 3) float

636 Points in the PointCloud

637 """

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

639

640 @vertices.setter

641 def vertices(self, values):

642 """

643 Assign vertex values to the point cloud.

644

645 Parameters

646 --------------

647 values : (n, 3) float

648 Points in space

649 """

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

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

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

653

654 @property

655 def colors(self):

656 """

657 Stored per- point color

658

659 Returns

660 ----------

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

662 Per- point RGBA color

663 """

664 return self.visual.vertex_colors

665

666 @colors.setter

667 def colors(self, data):

668 self.visual.vertex_colors = data

669

670 @caching.cache_decorator

671 def kdtree(self):

672 """

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

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

675

676 Returns

677 ---------

678 tree : scipy.spatial.cKDTree

679 Contains mesh.vertices

680 """

681

682 from scipy.spatial import cKDTree

683

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

685 return tree

686

687 @caching.cache_decorator

688 def convex_hull(self):

689 """

690 A convex hull of every point.

691

692 Returns

693 -------------

694 convex_hull : trimesh.Trimesh

695 A watertight mesh of the hull of the points

696 """

697 from . import convex

698

699 return convex.convex_hull(self.vertices)

700

701 def scene(self):

702 """

703 A scene containing just the PointCloud

704

705 Returns

706 ----------

707 scene : trimesh.Scene

708 Scene object containing this PointCloud

709 """

710 from .scene.scene import Scene

711

712 return Scene(self)

713

714 def show(self, **kwargs):

715 """

716 Open a viewer window displaying the current PointCloud

717 """

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

719

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

721 """

722 Export the current pointcloud to a file object.

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

724 Supported formats are xyz

725 Parameters

726 ------------

727 file_obj: open writeable file object

728 str, file name where to save the pointcloud

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

730 file_type: str

731 Which file type to export as.

732 If file name is passed this is not required

733 """

734 from .exchange.export import export_mesh

735

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

737

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

739 """

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

741 Parameters

742 ------------

743 input_points : (n, 3) float

744 Input query points

745 kwargs : dict

746 Arguments for proximity.query_from_points

747 result : proximity.NearestQueryResult

748 Result of the query.

749 """

750 from .proximity import query_from_points

751

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

753

754 def __add__(self, other):

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

756 colors = None

757 else:

758 # preserve colors

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

760 other_colors = (

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

762 if len(other.colors) == 0

763 else other.colors

764 )

765 self_colors = (

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

767 if len(self.colors) == 0

768 else self.colors

769 )

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

771 return PointCloud(

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

773 )