Coverage for trimesh/nsphere.py: 82%
55 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 04:40 +0000
1"""
2nsphere.py
3--------------
5Functions for fitting and minimizing nspheres:
6circles, spheres, hyperspheres, etc.
7"""
9import numpy as np
11from . import convex, util
12from .constants import log, tol
14try:
15 # scipy is a soft dependency
16 from scipy import spatial
17 from scipy.optimize import leastsq
18except BaseException as E:
19 # raise the exception when someone tries to use it
20 from . import exceptions
22 leastsq = exceptions.ExceptionWrapper(E)
23 spatial = exceptions.ExceptionWrapper(E)
26def minimum_nsphere(obj):
27 """
28 Compute the minimum n- sphere for a mesh or a set of points.
30 Uses the fact that the minimum n- sphere will be centered at one of
31 the vertices of the furthest site voronoi diagram, which is n*log(n)
32 but should be pretty fast due to using the scipy/qhull implementations
33 of convex hulls and voronoi diagrams.
35 Parameters
36 ----------
37 obj : (n, d) float or trimesh.Trimesh
38 Points or mesh to find minimum bounding nsphere
40 Returns
41 ----------
42 center : (d,) float
43 Center of fitted n- sphere
44 radius : float
45 Radius of fitted n-sphere
46 """
48 # reduce the input points or mesh to the vertices of the convex hull
49 # since we are computing the furthest site voronoi diagram this reduces
50 # the input complexity substantially and returns the same value
51 points = convex.hull_points(obj)
53 # we are scaling the mesh to a unit cube
54 # this used to pass qhull_options 'QbB' to Voronoi however this had a bug somewhere
55 # to avoid this we scale to a unit cube ourselves inside this function
56 points_origin = points.min(axis=0)
57 points_scale = np.ptp(points, axis=0).min()
58 points = (points - points_origin) / points_scale
60 # if all of the points are on an n-sphere already the voronoi
61 # method will fail so we check a least squares fit before
62 # bothering to compute the voronoi diagram
63 fit_C, fit_R, fit_E = fit_nsphere(points)
64 # return fit radius and center to global scale
65 fit_R = (((points - fit_C) ** 2).sum(axis=1).max() ** 0.5) * points_scale
66 fit_C = (fit_C * points_scale) + points_origin
68 if fit_E < 1e-6:
69 # points were on an n-sphere so just return fit
70 return fit_C, fit_R
72 # calculate a furthest site voronoi diagram
73 # this will fail if the points are ALL on the surface of
74 # the n-sphere but hopefully the least squares check caught those cases
75 # , qhull_options='QbB Pp')
76 voronoi = spatial.Voronoi(points, furthest_site=True)
78 # find the maximum radius^2 point for each of the voronoi vertices
79 # this is worst case quite expensive but we have taken
80 # convex hull to reduce n for this operation
81 # we are doing comparisons on the radius squared then rooting once
82 try:
83 # cdist is massivly faster than looping or tiling methods
84 # although it does create an extremely large intermediate array
85 # first, get an order of magnitude memory size estimate
86 # a float64 would be 8 bytes per entry plus overhead
87 memory_estimate = len(voronoi.vertices) * len(points) * 9
88 # if this `cdist` would require more than 4gb of memory fall
89 # back to the more memory-efficent but much slower loop
90 if memory_estimate > 4e9:
91 raise MemoryError
92 radii_2 = spatial.distance.cdist(
93 voronoi.vertices, points, metric="sqeuclidean"
94 ).max(axis=1)
95 except MemoryError:
96 # log the MemoryError
97 log.warning("MemoryError: falling back to slower check!")
98 # fall back to a potentially very slow list comprehension
99 radii_2 = np.array(
100 [((points - v) ** 2).sum(axis=1).max() for v in voronoi.vertices]
101 )
103 # we want the smallest sphere so take the min of the radii
104 radii_idx = radii_2.argmin()
106 # return voronoi radius and center to global scale
107 radius_v = np.sqrt(radii_2[radii_idx]) * points_scale
108 center_v = (voronoi.vertices[radii_idx] * points_scale) + points_origin
110 if radius_v > fit_R:
111 return fit_C, fit_R
113 return center_v, radius_v
116def fit_nsphere(points, prior=None):
117 """
118 Fit an n-sphere to a set of points using least squares.
120 Parameters
121 ------------
122 points : (n, d) float
123 Points in space
124 prior : (d,) float
125 Best guess for center of nsphere
127 Returns
128 ---------
129 center : (d,) float
130 Location of center
131 radius : float
132 Mean radius across circle
133 error : float
134 Peak to peak value of deviation from mean radius
135 """
136 # make sure points are numpy array
137 points = np.asanyarray(points, dtype=np.float64)
138 # create ones so we can dot instead of using slower sum
139 ones = np.ones(points.shape[1])
141 def residuals(center):
142 # do the axis sum with a dot
143 # this gets called a LOT so worth optimizing
144 radii_sq = np.dot((points - center) ** 2, ones)
145 # residuals are difference between mean
146 # use our sum mean vs .mean() as it is slightly faster
147 return radii_sq - (radii_sq.sum() / len(radii_sq))
149 if prior is None:
150 guess = points.mean(axis=0)
151 else:
152 guess = np.asanyarray(prior)
154 center_result, return_code = leastsq(residuals, guess, xtol=1e-8)
156 if return_code not in [1, 2, 3, 4]:
157 raise ValueError("Least square fit failed!")
159 radii = util.row_norm(points - center_result)
160 radius = radii.mean()
161 error = np.ptp(radii)
162 return center_result, radius, error
165def is_nsphere(points):
166 """
167 Check if a list of points is an nsphere.
169 Parameters
170 -----------
171 points : (n, dimension) float
172 Points in space
174 Returns
175 -----------
176 check : bool
177 True if input points are on an nsphere
178 """
179 _center, _radius, error = fit_nsphere(points)
180 check = error < tol.merge
181 return check