Skip to content

reach_for_the_spheres

ReachForTheSpheresState

An object to keep state during the iterations of the Reach for the Spheres method. Meant to be used in conjunction with reach_for_the_spheres_iteration.

Parameters:

Name Type Description Default
V (n,dim) numpy double array

Matrix of mesh vertices of the initial mesh

required
F (m,dim) numpy int array

Matrix of triangle indices into V of the initlal mesh

required
U (k,dim) numpy double array, optional (default: None)

Matrix of SDF sample positions. The sdf will be sampled at these points

None
S (k,dim) nummpy double array, optional (default: None)

Matrix of SDF samples at the sample positions in U. If this is not provided, it will be computed using sdf.

None
sdf

The signed distance function to be reconstructed. This function must take in a (ks,dim) matrix of sample points and return a (ks,) array of SDF values at these points.

None
V_active (na,dim) numpy double array, optional (default: None)

Matrix of mesh vertices active during iteration

None
F_active (ma,dim) numpy int array, optional (default: None)

Matrix of triangle indices into V_active of the mesh active during iteration

None
V_inactive (na,dim) numpy double array, optional (default: None)

Matrix of mesh vertices inactive during iteration

None
F_inactive (ma,dim) numpy int array, optional (default: None)

Matrix of triangle indices into V_inactive of the mesh inactive during iteration

None
rng numpy random generator, optional (default: None)

numpy rng object used to create randomness during iterations

None
h float (default None)

The method's initial target mesh length for the reconstructed mesh. This will change during iteration, set min_h if you want to control the minimum edge length overall.

None
min_h float (default None)

The method's minimal target edge length for the reconstructed mesh.

None
best_performance float (default None)

Remembers the method's best performance

None
convergence_counter int (default None)

The method's convergence counter, used to track for how long progress has not been made.

None
best_avg_error float (default None)

Remembers the method's average error.

None
resample_counter int (default None)

Tracks how often the SDF has been resampled.

None
V_last_converged (nl,dim) numpy double array, optional (default: None)

Matrix of mesh vertices for the last known converged mesh.

None
F_last_converged (ml,dim) numpy int array, optional (default: None)

Matrix of triangle indices into V_last_converged for the last known converged mesh.

None
U_batch (kb,dim) numpy double array, optional (default: None)

Matrix of SDF sample positions as used during the last batching operation.

None
S_batch (kb,dim) nummpy double array, optional (default: None)

Matrix of SDF samples at the sample positions in U as used during the last batching operation.

None
Notes

If you create this object yourself, you should supply V, F, U, S, sdf. Supply all other parameters only as explicitly needed.

Examples:

import gpytoolbox as gpy
import numpy as npy

# Get a signed distance function
V,F = gpy.read_mesh("my_mesh.obj")

# Create an SDF for the mesh
j = 20
sdf = lambda x: gpy.signed_distance(x, V, F)[0]
gx, gy, gz = np.meshgrid(np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1))
U = np.vstack((gx.flatten(), gy.flatten(), gz.flatten())).T

# Create ReachForTheSpheresState to use in iterative method later
V0, F0 = gpy.icosphere(2)
state = ReachForTheSpheresState(V=V0, F=F0, sdf=sdf, U=U)

# Run one iteration of Reach for the Spheres
converged = gpy.reach_for_the_spheres_iteration(state)

# Reconstruct triangle mesh
Vr,Fr = state.V,state.F

#The reconstructed mesh after one iteration is now Vr,Fr
Source code in src/gpytoolbox/reach_for_the_spheres.py
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
class ReachForTheSpheresState:
    """An object to keep state during the iterations of the Reach for the
    Spheres method.
    Meant to be used in conjunction with `reach_for_the_spheres_iteration`.

    Parameters
    ----------
    V : (n,dim) numpy double array
        Matrix of mesh vertices of the initial mesh
    F : (m,dim) numpy int array
        Matrix of triangle indices into V of the initlal mesh
    U : (k,dim) numpy double array, optional (default: None)
        Matrix of SDF sample positions.
        The sdf will be sampled at these points
    S : (k,dim) nummpy double array, optional (default: None)
        Matrix of SDF samples at the sample positions in U.
        If this is not provided, it will be computed using sdf.
    sdf: lambda function that takes in a (ks,dim) numpy double array, optional (default: None)
        The signed distance function to be reconstructed.
        This function must take in a (ks,dim) matrix of sample points and
        return a (ks,) array of SDF values at these points.
    V_active : (na,dim) numpy double array, optional (default: None)
        Matrix of mesh vertices active during iteration
    F_active : (ma,dim) numpy int array, optional (default: None)
        Matrix of triangle indices into V_active of the mesh active during iteration
    V_inactive : (na,dim) numpy double array, optional (default: None)
        Matrix of mesh vertices inactive during iteration
    F_inactive : (ma,dim) numpy int array, optional (default: None)
        Matrix of triangle indices into V_inactive of the mesh inactive during iteration
    rng : numpy random generator, optional (default: None)
        numpy rng object used to create randomness during iterations
    h : float (default None)
        The method's initial target mesh length for the reconstructed mesh.
        This will change during iteration, set min_h if you want to control the
        minimum edge length overall.
    min_h : float (default None)
        The method's minimal target edge length for the reconstructed mesh.
    best_performance : float (default None)
        Remembers the method's best performance
    convergence_counter : int (default None)
        The method's convergence counter, used to track for how long progress
        has not been made.
    best_avg_error : float (default None)
        Remembers the method's average error.
    resample_counter : int (default None)
        Tracks how often the SDF has been resampled.
    V_last_converged : (nl,dim) numpy double array, optional (default: None)
        Matrix of mesh vertices for the last known converged mesh.
    F_last_converged : (ml,dim) numpy int array, optional (default: None)
        Matrix of triangle indices into V_last_converged for the last known
        converged mesh.
    U_batch : (kb,dim) numpy double array, optional (default: None)
        Matrix of SDF sample positions as used during the last batching
        operation.
    S_batch : (kb,dim) nummpy double array, optional (default: None)
        Matrix of SDF samples at the sample positions in U as used during the
        last batching operation.

    Notes
    --------
    If you create this object yourself, you should supply V, F, U, S, sdf.
    Supply all other parameters only as explicitly needed.

    Examples
    --------
    ```python
    import gpytoolbox as gpy
    import numpy as npy

    # Get a signed distance function
    V,F = gpy.read_mesh("my_mesh.obj")

    # Create an SDF for the mesh
    j = 20
    sdf = lambda x: gpy.signed_distance(x, V, F)[0]
    gx, gy, gz = np.meshgrid(np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1))
    U = np.vstack((gx.flatten(), gy.flatten(), gz.flatten())).T

    # Create ReachForTheSpheresState to use in iterative method later
    V0, F0 = gpy.icosphere(2)
    state = ReachForTheSpheresState(V=V0, F=F0, sdf=sdf, U=U)

    # Run one iteration of Reach for the Spheres
    converged = gpy.reach_for_the_spheres_iteration(state)

    # Reconstruct triangle mesh
    Vr,Fr = state.V,state.F

    #The reconstructed mesh after one iteration is now Vr,Fr
    ```
    """

    def __init__(self,
        V, F,
        U=None,
        S=None,
        sdf=None,
        V_active=None, F_active=None,
        V_inactive=None, F_inactive=None,
        rng=None,
        h=None, min_h=None, its=None,
        best_performance=None,
        convergence_counter=None,
        best_avg_error=None,
        resample_counter=None,
        V_last_converged=None, F_last_converged=None,
        U_batch=None, S_batch=None
        ):

        self.V=V
        self.F=F
        self.U=U
        self.S=S
        self.sdf=sdf
        self.V_active = V_active
        self.F_active = F_active
        self.V_inactive = V_inactive
        self.F_inactive = F_inactive
        self.rng = rng
        self.h = h
        self.min_h = min_h
        self.its = its
        self.best_performance = best_performance
        self.convergence_counter = convergence_counter
        self.best_avg_error = best_avg_error
        self.resample_counter = resample_counter
        self.V_last_converged = V_last_converged
        self.F_last_converged = F_last_converged
        self.U_batch = U_batch
        self.S_batch = S_batch

reach_for_the_spheres(U, sdf, V, F, S=None, return_U=False, verbose=False, max_iter=None, tol=None, h=None, min_h=None, linesearch=None, min_t=None, max_t=None, dt=None, inside_outside_test=None, resample=None, resample_samples=None, feature_detection=None, output_sensitive=None, remesh_iterations=None, batch_size=None, fix_boundary=None, clamp=None, pseudosdf_interior=None)

Creates a mesh from a signed distance function (SDF) using the "Reach for the Spheres" method of S. Sellán, C. Batty, and O. Stein [2023].

This method takes in an sdf, sample points (do not need to be on a grid), and an initial mesh. It then flows this initial mesh into a reconstructed mesh that fulfills the signed distance function.

This method works in dimension 2 (dim==2), where it reconstructs a polyline, and in dimension 3 (dim==3), where it reconstructs a triangle mesh.

Parameters:

Name Type Description Default
U (k,dim) numpy double array

Matrix of SDF sample positions. The sdf will be sampled at these points

required
sdf

The signed distance function to be reconstructed. This function must take in a (ks,dim) matrix of sample points and return a (ks,) array of SDF values at these points.

required
V (n,dim) numpy double array

Matrix of mesh vertices of the initial mesh

required
F (m,dim) numpy int array

Matrix of triangle indices into V of the initlal mesh

required
S (k,dim) nummpy double array, optional (default: None)

Matrix of SDF samples at the sample positions in U. If this is not provided, it will be computed using sdf.

None
return_U bool, optional (default: False)

Whether to return the matrix of SDF sample positions along with the reconstructed mesh.

False
verbose bool (default false)

Whether to print method statistics during operation.

False
max_iter int (default None)

The maximum number of iterations to perform for the method. If not supplied, a sensible default is used.

None
tol float (default None)

The method's tolerance for the sphere tangency test. If not supplied, a sensible default is used.

None
h float (default None)

The method's initial target mesh length for the reconstructed mesh. This will change during iteration, set min_h if you want to control the minimum edge length overall. If not supplied, a sensible default is used.

None
min_h float (default None)

The method's minimal target edge length for the reconstructed mesh. If not supplied, a sensible default is used.

None
linesearch bool (default None)

Whether to use a linesearch-like heuristic for the method's timestep. If not supplied, linesearch is used.

None
min_t float (default None)

The method's minimum timestep. If not supplied, a sensible default is used.

None
max_t float (default None)

The method's minimum timestep. If not supplied, a sensible default is used.

None
dt float (default None)

The method's default timestep. If not supplied, a sensible default is used.

None
inside_outside_test bool (default None)

Whether to use inside-outside testing when projecting points to be tangent to the sphere. Turn this off if your distance function is unsigned. If not supplied, inside-outside test is used

None
resample int (default None)

How often to resample the SDF after convergence to extract more information. If not supplied, resampling is not performed.

None
resample_samples int (default None)

How many samples to use when resampling. If not supplied, a sensible default is used.

None
feature_detection string (default None)

Which feature detection mode to use. If not supplied, aggressive feature detection is used.

None
output_sensitive bool (default None)

Whether to use output-sensitive remeshing. If not supplied, remeshing is output-sensitive.

None
remesh_iterations int (default None)

How many iterations of the remesher to run each step. If not supplied, a sensible default is used.

None
batch_size int (default None)

For large amounts of sample points, the method is sped up using sample point batching. This parameter specifies how many samples to take for each batch. Set it to 0 to disable batching. If not supplied, a sensible default is used.

None
fix_boundary bool (default None)

Whether to fix the boundary of the mesh during iteration. If not supplied, the boundary is not fixed.

None
clamp float (default None)

If sdf is a clamped SDF, the clamp value to use. np.inf for no clamping. If not supplied, there is no clamping.

None
pseudosdf_interior bool (default None)

If enabled, treats every negative SDF value as a bound on the signed distance, as opposed to an exact signed distance, for use in SDFs resulting from CSG unions as described by Marschner et al. "Constructive Solid Geometry on Neural Signed Distance Fields" [2023]. If not supplied, this feature is disabled.

None

Returns:

Name Type Description
Vr (nr,dim) numpy double array

Matrix of mesh vertices of the reconstructed triangle mesh

Fr (mr,dim) numpy int array

Matrix of triangle indices into Vr of the reconstructed mesh

Ur (kr,dim) numpy double array, if requested

Matrix of SDF sample positions. This can be different from the supplied Ur if the method is set to resample.

See Also

marching_squares, marching_cubes

Notes

This method has a number of limitations that are described in the paper. E.g., the method will only work for SDFs that describe surfaces with the same topology as the initial surface.

Examples:

import gpytoolbox as gpy
import numpy as npy

# Get a signed distance function
V,F = gpy.read_mesh("my_mesh.obj")

# Create an SDF for the mesh
j = 20
sdf = lambda x: gpy.signed_distance(x, V, F)[0]
gx, gy, gz = np.meshgrid(np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1))
U = np.vstack((gx.flatten(), gy.flatten(), gz.flatten())).T

# Choose an initial surface for reach_for_the_spheres
V0, F0 = gpy.icosphere(2)

# Reconstruct triangle mesh
Vr,Fr = gpy.reach_for_the_spheres(U, sdf, V0, F0)

#The reconstructed triangle mesh is now Vr,Fr.
Source code in src/gpytoolbox/reach_for_the_spheres.py
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
def reach_for_the_spheres(U, sdf, V, F, S=None,
    return_U=False,
    verbose=False,
    max_iter=None, tol=None, h=None, min_h=None,
    linesearch=None, min_t=None, max_t=None,
    dt=None,
    inside_outside_test=None,
    resample=None,
    resample_samples=None,
    feature_detection=None,
    output_sensitive=None,
    remesh_iterations=None,
    batch_size=None,
    fix_boundary=None,
    clamp=None, pseudosdf_interior=None):
    """Creates a mesh from a signed distance function (SDF) using the
    "Reach for the Spheres" method of S. Sellán,  C. Batty, and O. Stein [2023].

    This method takes in an sdf, sample points (do not need to be on a grid),
    and an initial mesh.
    It then flows this initial mesh into a reconstructed mesh that fulfills
    the signed distance function.

    This method works in dimension 2 (dim==2), where it reconstructs a polyline,
    and in dimension 3 (dim==3), where it reconstructs a triangle mesh.

    Parameters
    ----------
    U : (k,dim) numpy double array
        Matrix of SDF sample positions.
        The sdf will be sampled at these points
    sdf: lambda function that takes in a (ks,dim) numpy double array
        The signed distance function to be reconstructed.
        This function must take in a (ks,dim) matrix of sample points and
        return a (ks,) array of SDF values at these points.
    V : (n,dim) numpy double array
        Matrix of mesh vertices of the initial mesh
    F : (m,dim) numpy int array
        Matrix of triangle indices into V of the initlal mesh
    S : (k,dim) nummpy double array, optional (default: None)
        Matrix of SDF samples at the sample positions in U.
        If this is not provided, it will be computed using sdf.
    return_U : bool, optional (default: False)
        Whether to return the matrix of SDF sample positions along with the
        reconstructed mesh.
    verbose : bool (default false)
        Whether to print method statistics during operation.
    max_iter : int (default None)
        The maximum number of iterations to perform for the method.
        If not supplied, a sensible default is used.
    tol : float (default None)
        The method's tolerance for the sphere tangency test.
        If not supplied, a sensible default is used.
    h : float (default None)
        The method's initial target mesh length for the reconstructed mesh.
        This will change during iteration, set min_h if you want to control the
        minimum edge length overall.
        If not supplied, a sensible default is used.
    min_h : float (default None)
        The method's minimal target edge length for the reconstructed mesh.
        If not supplied, a sensible default is used.
    linesearch : bool (default None)
        Whether to use a linesearch-like heuristic for the method's timestep.
        If not supplied, linesearch is used.
    min_t : float (default None)
        The method's minimum timestep.
        If not supplied, a sensible default is used.
    max_t : float (default None)
        The method's minimum timestep.
        If not supplied, a sensible default is used.
    dt : float (default None)
        The method's default timestep.
        If not supplied, a sensible default is used.
    inside_outside_test : bool (default None)
        Whether to use inside-outside testing when projecting points to be
        tangent to the sphere.
        Turn this off if your distance function is *unsigned*.
        If not supplied, inside-outside test is used
    resample : int (default None)
        How often to resample the SDF after convergence to extract more
        information.
        If not supplied, resampling is not performed.
    resample_samples : int (default None)
        How many samples to use when resampling.
        If not supplied, a sensible default is used.
    feature_detection : string (default None)
        Which feature detection mode to use.
        If not supplied, aggressive feature detection is used.
    output_sensitive : bool (default None)
        Whether to use output-sensitive remeshing.
        If not supplied, remeshing is output-sensitive.
    remesh_iterations : int (default None)
        How many iterations of the remesher to run each step.
        If not supplied, a sensible default is used.
    batch_size : int (default None)
        For large amounts of sample points, the method is sped up using sample
        point batching.
        This parameter specifies how many samples to take for each batch.
        Set it to 0 to disable batching.
        If not supplied, a sensible default is used.
    fix_boundary : bool (default None)
        Whether to fix the boundary of the mesh during iteration.
        If not supplied, the boundary is not fixed.
    clamp : float (default None)
        If sdf is a clamped SDF, the clamp value to use.
        np.inf for no clamping.
        If not supplied, there is no clamping.
    pseudosdf_interior : bool (default None)
        If enabled, treats every negative SDF value as a bound on the signed
        distance, as opposed to an exact signed distance, for use in SDFs
        resulting from CSG unions as described by Marschner et al.
        "Constructive Solid Geometry on Neural Signed Distance Fields" [2023].
        If not supplied, this feature is disabled.

    Returns
    -------
    Vr : (nr,dim) numpy double array
        Matrix of mesh vertices of the reconstructed triangle mesh
    Fr : (mr,dim) numpy int array
        Matrix of triangle indices into Vr of the reconstructed mesh
    Ur : (kr,dim) numpy double array, if requested
        Matrix of SDF sample positions.
        This can be different from the supplied Ur if the method is set to
        resample.

    See Also
    --------
    marching_squares, marching_cubes

    Notes
    --------
    This method has a number of limitations that are described in the paper.
    E.g., the method will only work for SDFs that describe surfaces with the
    same topology as the initial surface.

    Examples
    --------
    ```python
    import gpytoolbox as gpy
    import numpy as npy

    # Get a signed distance function
    V,F = gpy.read_mesh("my_mesh.obj")

    # Create an SDF for the mesh
    j = 20
    sdf = lambda x: gpy.signed_distance(x, V, F)[0]
    gx, gy, gz = np.meshgrid(np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1))
    U = np.vstack((gx.flatten(), gy.flatten(), gz.flatten())).T

    # Choose an initial surface for reach_for_the_spheres
    V0, F0 = gpy.icosphere(2)

    # Reconstruct triangle mesh
    Vr,Fr = gpy.reach_for_the_spheres(U, sdf, V0, F0)

    #The reconstructed triangle mesh is now Vr,Fr.
    ```
    """


    state = ReachForTheSpheresState(V=V, F=F, sdf=sdf, U=U, S=S,
        h=h, min_h=min_h)
    converged = False

    while not converged:
        converged = reach_for_the_spheres_iteration(state,
            max_iter=max_iter, tol=tol,
            linesearch=linesearch, min_t=min_t, max_t=max_t,
            dt=dt,
            inside_outside_test=inside_outside_test,
            resample=resample,
            resample_samples=resample_samples,
            feature_detection=feature_detection,
            output_sensitive=output_sensitive,
            remesh_iterations=remesh_iterations,
            batch_size=batch_size,
            verbose=verbose,
            fix_boundary=fix_boundary,
            clamp=clamp, pseudosdf_interior=pseudosdf_interior)

    if return_U:
        return state.V, state.F, state.U
    else:
        return state.V, state.F

reach_for_the_spheres_iteration(state, max_iter=None, tol=None, linesearch=None, min_t=None, max_t=None, dt=None, inside_outside_test=None, resample=None, resample_samples=None, feature_detection=None, output_sensitive=None, remesh_iterations=None, batch_size=None, fix_boundary=None, clamp=None, pseudosdf_interior=None, verbose=False)

Performs one iteration of the "Reach for the Spheres" method of S. Sellán, C. Batty, and O. Stein [2023]. This method is used to create a mesh from a signed distance function (SDF).

This method takes in the current state of the method in the form of a ReachForTheSpheresState object, and stores its results as well as any temporary information needed in that state object.

This method works in dimension 2 (dim==2), where it reconstructs a polyline, and in dimension 3 (dim==3), where it reconstructs a triangle mesh.

Parameters:

Name Type Description Default
state

Stores all needed information about the current state of the method.

required
max_iter int (default None)

The maximum number of iterations to perform for the method. If not supplied, a sensible default is used.

None
tol float (default None)

The method's tolerance for the sphere tangency test. If not supplied, a sensible default is used.

None
linesearch bool (default None)

Whether to use a linesearch-like heuristic for the method's timestep. If not supplied, linesearch is used.

None
min_t float (default None)

The method's minimum timestep. If not supplied, a sensible default is used.

None
max_t float (default None)

The method's minimum timestep. If not supplied, a sensible default is used.

None
dt float (default None)

The method's default timestep. If not supplied, a sensible default is used.

None
inside_outside_test bool (default None)

Whether to use inside-outside testing when projecting points to be tangent to the sphere. Turn this off if your distance function is unsigned. If not supplied, inside-outside test is used

None
resample int (default None)

How often to resample the SDF after convergence to extract more information. If not supplied, resampling is not performed.

None
resample_samples int (default None)

How many samples to use when resampling. If not supplied, a sensible default is used.

None
feature_detection string (default None)

Which feature detection mode to use. If not supplied, aggressive feature detection is used.

None
output_sensitive bool (default None)

Whether to use output-sensitive remeshing. If not supplied, remeshing is output-sensitive.

None
remesh_iterations int (default None)

How many iterations of the remesher to run each step. If not supplied, a sensible default is used.

None
batch_size int (default None)

For large amounts of sample points, the method is sped up using sample point batching. This parameter specifies how many samples to take for each batch. Set it to 0 to disable batching. If not supplied, a sensible default is used.

None
fix_boundary bool (default None)

Whether to fix the boundary of the mesh during iteration. If not supplied, the boundary is not fixed.

None
clamp float (default None)

If sdf is a clamped SDF, the clamp value to use. np.inf for no clamping. If not supplied, there is no clamping.

None
pseudosdf_interior bool (default None)

If enabled, treats every negative SDF value as a bound on the signed distance, as opposed to an exact signed distance, for use in SDFs resulting from CSG unions as described by Marschner et al. "Constructive Solid Geometry on Neural Signed Distance Fields" [2023]. If not supplied, this feature is disabled.

None
verbose bool (default false)

Whether to print method statistics during operation.

False

Returns:

Name Type Description
converged bool

Whether the method has converged after this iteration or not.

See Also

reach_for_the_spheres, marching_squares, marching_cubes

Notes

This method has a number of limitations that are described in the paper. E.g., the method will only work for SDFs that describe surfaces with the same topology as the initial surface.

Examples:

import gpytoolbox as gpy
import numpy as npy

# Get a signed distance function
V,F = gpy.read_mesh("my_mesh.obj")

# Create an SDF for the mesh
j = 20
sdf = lambda x: gpy.signed_distance(x, V, F)[0]
gx, gy, gz = np.meshgrid(np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1))
U = np.vstack((gx.flatten(), gy.flatten(), gz.flatten())).T

# Create ReachForTheSpheresState to use in iterative method later
V0, F0 = gpy.icosphere(2)
state = ReachForTheSpheresState(V=V0, F=F0, sdf=sdf, U=U)

# Run one iteration of Reach for the Spheres
converged = gpy.reach_for_the_spheres_iteration(state)

# Reconstruct triangle mesh
Vr,Fr = state.V,state.F

#The reconstructed mesh after one iteration is now Vr,Fr
Source code in src/gpytoolbox/reach_for_the_spheres.py
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
def reach_for_the_spheres_iteration(state,
    max_iter=None, tol=None,
    linesearch=None, min_t=None, max_t=None,
    dt=None,
    inside_outside_test=None,
    resample=None,
    resample_samples=None,
    feature_detection=None,
    output_sensitive=None,
    remesh_iterations=None,
    batch_size=None,
    fix_boundary=None,
    clamp=None, pseudosdf_interior=None,
    verbose=False):
    """Performs one iteration of the "Reach for the Spheres" method of
    S. Sellán,  C. Batty, and O. Stein [2023].
    This method is used to create a mesh from a signed distance function (SDF).

    This method takes in the current state of the method in the form of a 
    ReachForTheSpheresState object, and stores its results as well as any
    temporary information needed in that state object.

    This method works in dimension 2 (dim==2), where it reconstructs a polyline,
    and in dimension 3 (dim==3), where it reconstructs a triangle mesh.

    Parameters
    ----------
    state: ReachForTheSpheresState object
        Stores all needed information about the current state of the method.
    max_iter : int (default None)
        The maximum number of iterations to perform for the method.
        If not supplied, a sensible default is used.
    tol : float (default None)
        The method's tolerance for the sphere tangency test.
        If not supplied, a sensible default is used.
    linesearch : bool (default None)
        Whether to use a linesearch-like heuristic for the method's timestep.
        If not supplied, linesearch is used.
    min_t : float (default None)
        The method's minimum timestep.
        If not supplied, a sensible default is used.
    max_t : float (default None)
        The method's minimum timestep.
        If not supplied, a sensible default is used.
    dt : float (default None)
        The method's default timestep.
        If not supplied, a sensible default is used.
    inside_outside_test : bool (default None)
        Whether to use inside-outside testing when projecting points to be
        tangent to the sphere.
        Turn this off if your distance function is *unsigned*.
        If not supplied, inside-outside test is used
    resample : int (default None)
        How often to resample the SDF after convergence to extract more
        information.
        If not supplied, resampling is not performed.
    resample_samples : int (default None)
        How many samples to use when resampling.
        If not supplied, a sensible default is used.
    feature_detection : string (default None)
        Which feature detection mode to use.
        If not supplied, aggressive feature detection is used.
    output_sensitive : bool (default None)
        Whether to use output-sensitive remeshing.
        If not supplied, remeshing is output-sensitive.
    remesh_iterations : int (default None)
        How many iterations of the remesher to run each step.
        If not supplied, a sensible default is used.
    batch_size : int (default None)
        For large amounts of sample points, the method is sped up using sample
        point batching.
        This parameter specifies how many samples to take for each batch.
        Set it to 0 to disable batching.
        If not supplied, a sensible default is used.
    fix_boundary : bool (default None)
        Whether to fix the boundary of the mesh during iteration.
        If not supplied, the boundary is not fixed.
    clamp : float (default None)
        If sdf is a clamped SDF, the clamp value to use.
        np.inf for no clamping.
        If not supplied, there is no clamping.
    pseudosdf_interior : bool (default None)
        If enabled, treats every negative SDF value as a bound on the signed
        distance, as opposed to an exact signed distance, for use in SDFs
        resulting from CSG unions as described by Marschner et al.
        "Constructive Solid Geometry on Neural Signed Distance Fields" [2023].
        If not supplied, this feature is disabled.
    verbose : bool (default false)
        Whether to print method statistics during operation.

    Returns
    -------
    converged : bool
        Whether the method has converged after this iteration or not.

    See Also
    --------
    reach_for_the_spheres, marching_squares, marching_cubes

    Notes
    --------
    This method has a number of limitations that are described in the paper.
    E.g., the method will only work for SDFs that describe surfaces with the
    same topology as the initial surface.

    Examples
    --------
    ```python
    import gpytoolbox as gpy
    import numpy as npy

    # Get a signed distance function
    V,F = gpy.read_mesh("my_mesh.obj")

    # Create an SDF for the mesh
    j = 20
    sdf = lambda x: gpy.signed_distance(x, V, F)[0]
    gx, gy, gz = np.meshgrid(np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1), np.linspace(-1.0, 1.0, j+1))
    U = np.vstack((gx.flatten(), gy.flatten(), gz.flatten())).T

    # Create ReachForTheSpheresState to use in iterative method later
    V0, F0 = gpy.icosphere(2)
    state = ReachForTheSpheresState(V=V0, F=F0, sdf=sdf, U=U)

    # Run one iteration of Reach for the Spheres
    converged = gpy.reach_for_the_spheres_iteration(state)

    # Reconstruct triangle mesh
    Vr,Fr = state.V,state.F

    #The reconstructed mesh after one iteration is now Vr,Fr
    ```
    """


    assert isinstance(state, ReachForTheSpheresState), "State must be a ReachForTheSpheresState"
    dim = state.V.shape[1]
    assert dim==state.F.shape[1]
    if state.U is not None:
        assert dim==state.U.shape[1]

    #Set default parameters.
    #Assign different 2D and 3D default parameters
    default_params = {
        #For each dimension
        2: {'max_iter':10000, 'tol':5e-3, 'h':0.1,
        'linesearch':True, 'min_t':1e-6, 'max_t':50.,
        'dt':10.,
        'inside_outside_test':True,
        'resample':0, 'resample_samples':2*int(np.ceil(np.sqrt(state.V.shape[0] if state.U is None else state.U.shape[0]))),
        'return_U':False,
        'feature_detection':'aggressive', 'output_sensitive':True,
        'remesh_iterations':1,
        'batch_size':20000,
        'fix_boundary':False,
        'clamp':np.Inf, 'pseudosdf_interior':False},
        3: {'max_iter':20000, 'tol':1e-2, 'h':0.2,
        'linesearch':True, 'min_t':1e-6, 'max_t':50.,
        'dt':10.,
        'inside_outside_test':True,
        'resample':0, 'resample_samples':2*int(np.ceil(np.cbrt(state.V.shape[0] if state.U is None else state.U.shape[0]))),
        'feature_detection':'aggressive', 'output_sensitive':True,
        'remesh_iterations':1,
        'visualize':False,
        'batch_size':20000,
        'fix_boundary':False,
        'clamp':np.Inf, 'pseudosdf_interior':False}
    }
    if max_iter is None:
        max_iter = default_params[dim]['max_iter']
    if tol is None:
        tol = default_params[dim]['tol']
    if linesearch is None:
        linesearch = default_params[dim]['linesearch']
    if min_t is None:
        min_t = default_params[dim]['min_t']
    if max_t is None:
        max_t = default_params[dim]['max_t']
    if dt is None:
        dt = default_params[dim]['dt']
    if inside_outside_test is None:
        inside_outside_test = default_params[dim]['inside_outside_test']
    if resample is None:
        resample = default_params[dim]['resample']
    if resample_samples is None:
        resample_samples = default_params[dim]['resample_samples']
    if feature_detection is None:
        feature_detection = default_params[dim]['feature_detection']
    if output_sensitive is None:
        output_sensitive = default_params[dim]['output_sensitive']
    if remesh_iterations is None:
        remesh_iterations = default_params[dim]['remesh_iterations']
    if batch_size is None:
        batch_size = default_params[dim]['batch_size']
    if fix_boundary is None:
        fix_boundary = default_params[dim]['fix_boundary']
    if clamp is None:
        clamp = default_params[dim]['clamp']
    if pseudosdf_interior is None:
        pseudosdf_interior = default_params[dim]['pseudosdf_interior']

    if state.h is None:
        state.h = default_params[dim]['h']
    if state.U is None:
        assert state.sdf is not None, "If you do not provide U, you must provide an sdf function to sample."
        state.U = _sample_sdf(state.sdf, state.V, state.F)
    if state.S is None:
        assert state.sdf is not None, "If you do not provide S, you must provide an sdf function to sample."
        state.S = state.sdf(state.U)
    if state.min_h is None:
        # use a kdtree and get the average distance between two samples
        # import cKDTree
        tree = cKDTree(state.U)
        dists, _ = tree.query(state.U, k=2)
        state.min_h = 2.*np.mean(dists[:,1])
        state.min_h = np.clip(state.min_h, 0.001, 0.1)

    nu_0 = state.U.shape[0]
    nu_max = 2*nu_0

    if state.rng is None:
        state.rng = np.random.default_rng(68)
    if state.its is None:
        state.its = 0
    if state.best_performance is None:
        state.best_performance = np.Inf
    if state.convergence_counter is None:
        state.convergence_counter = 0
    if state.best_avg_error is None:
        state.best_avg_error = np.Inf
    # if state.use_features is None:
    #     state.use_features = False
    if state.V_last_converged is None:
        state.V_last_converged = state.V.copy()
    if state.F_last_converged is None:
        state.F_last_converged = state.F.copy()
    if state.resample_counter is None:
        state.resample_counter = 0
    if state.U_batch is None:
        state.U_batch = state.U.copy()
    if state.S_batch is None:
        state.S_batch = state.S.copy()

    remeshing = True

    #Do we prematurely abort? Increment iteration.
    if state.its>=max_iter:
        return True

    #Algorithm
    if batch_size>0 and batch_size<state.U.shape[0]:
        inds = state.rng.choice(state.U.shape[0], batch_size, replace=False)
        state.U_batch = state.U[inds,:]
        state.S_batch = state.S[inds]
        # include all inside points
        state.U_batch = np.concatenate((state.U_batch, state.U[state.S<=0,:]), axis=0)
        state.S_batch = np.concatenate((state.S_batch, state.S[state.S<=0]), axis=0)
    d2, I, b = squared_distance(state.U_batch, state.V, state.F,
        use_cpp=True, use_aabb=True)
    d = np.sqrt(d2)
    g = np.abs(state.S_batch)-d

    #closest point on edge to u
    pe = np.sum(state.V[state.F[I,:],:]*b[...,None], axis=1)
    pemU = pe-state.U_batch
    if inside_outside_test:
        s = -np.sign(np.sum(pemU*_normals(state.V,state.F)[I,:], axis=-1))
        hit_sides = np.any(b<1e-3, axis=-1) #too close to vertex
        s[hit_sides] = np.sign(state.S_batch)[hit_sides]
    else:
        s = np.sign(state.S_batch)
    #closest signed point on sphere to pe
    ps = state.U_batch+pemU*((s*state.S_batch/np.maximum(0.5*tol,d))[:,None])
    valid_U = np.abs(g) < tol
    invalid_U = np.logical_not(valid_U)
    F_invalid = np.unique(I[invalid_U])
    if feature_detection=='aggressive':
        V_invalid = np.unique(state.F[F_invalid,:].ravel())
        feature_vertices = np.setdiff1d(np.arange(state.V.shape[0]), V_invalid)
    else:
        F_valid = np.setdiff1d(np.arange(state.F.shape[0]), F_invalid)
        feature_vertices = np.unique(state.F[F_valid,:])

    #Build matrices
    wu = np.ones(state.U_batch.shape[0])
    clamped_g = np.where((np.abs(state.S_batch)==clamp)*(g<0.))
    if pseudosdf_interior:
        clamped_g = np.where((state.S_batch>0)*(g<0.))
    wu[clamped_g] = 0.0
    A = sp.sparse.csc_matrix(((wu[:,None]*b).ravel(),
        (np.tile(np.arange(state.U_batch.shape[0])[:,None],
            (1,state.F.shape[1])).ravel(),
            state.F[I,:].ravel())),
        (state.U_batch.shape[0],state.V.shape[0]))
    c = wu[:,None]*ps
    M = _massmatrix(state.V, state.F)
    rho = 1.0*np.ones(state.U_batch.shape[0]) / A.shape[0]
    R = sp.sparse.spdiags([rho], 0, rho.shape[0], rho.shape[0])

    #Compute timestep t
    if linesearch:
        # Backtracking linesearch
        # Nocedal & Wright Algorithm 3.1
        n_c = 0.01
        n_p = -A.transpose()*R*(A*state.V-c)
        t = -(np.sum((A*state.V-c)*(R*(A*n_p))) + n_c*np.sum(n_p**2)) \
        / np.sum((A*n_p)*(R*(A*n_p)))
        t = np.nan_to_num(t, nan=0., posinf=0., neginf=0.)
        t = min(max_t, max(t,min_t))
    else:
        t = dt

    #Minimize
    Q = M + t*(A.transpose()*R*A)
    b = M*state.V + t*A.transpose()*(R*c)
    if fix_boundary:
        bd = boundary_vertices(state.F)
        precomp = fixed_dof_solve_precompute(Q, k=bd)
        state.V = precomp.solve(b=b, y=state.V[bd,:])
    else:
        state.V = sp.sparse.linalg.spsolve(Q,b)

    # catching flow singularities so we fail gracefully

    there_are_non_manifold_edges = False
    if dim==3:
        there_are_non_manifold_edges = len(non_manifold_edges(state.F))>0
    elif dim==2:
        he_nm = np.sort(state.F, axis=1)
        # print(he)
        he_u_nm = np.unique(he_nm, axis=0, return_counts=True)
        # print(he)
        ne_nm = he_u_nm[0][he_u_nm[1]>2]
        there_are_non_manifold_edges = len(ne_nm)>0

    if np.any((np.isnan(state.V))) or np.any(doublearea(state.V, state.F)==0) or there_are_non_manifold_edges : 

        if verbose:
            print("we found a flow singularity. Returning the last converged solution.")
        state.V = state.V_last_converged.copy()
        state.F = state.F_last_converged.copy()
        return True

    #Convergence determination
    state.avg_error = np.linalg.norm(A*state.V-c) / A.shape[0]
    if verbose:
        print("Iteration:", state.its, "Counter:",state.convergence_counter, "h:",state.h, "Average error:", state.avg_error, "Best avg error:",state.best_avg_error, "Max error:", np.max(np.linalg.norm(A*state.V-c,axis=1)))
    if state.avg_error+1e-3*tol >= state.best_avg_error:
        state.convergence_counter = state.convergence_counter + 1
    else:
        state.convergence_counter = 0
        state.best_avg_error = state.avg_error
    if state.convergence_counter > 10:
        # if h==min_h:
        #     remeshing = False
        if state.h>state.min_h:
            state.V_last_converged = state.V.copy()
            state.F_last_converged = state.F.copy()
            state.best_avg_error = np.Inf
            state.convergence_counter = 0
        state.h = np.maximum(state.h/2,state.min_h)
    if state.convergence_counter > 100 or F_invalid.shape[0] == 0:
        if state.resample_counter<resample:
            state.U = _sample_sdf(state.sdf, state.V, state.F, state.U,
                new_n=resample_samples,
                trial_n=int(50*resample_samples), max_n=nu_max,
                remove_samples=True, keep_these_samples=np.arange(nu_0),
                rng=state.rng)
            state.S = state.sdf(state.U)
            state.U_batch = state.U.copy()
            state.S_batch = state.S.copy()
            state.resample_counter += 1
            state.best_performance = np.Inf
            state.convergence_counter = 0
            state.best_avg_error = np.Inf
            if verbose:
                print(f"Resampled, I now have {state.U.shape[0]} sample points.")
        else:
            # We have converged.
            return True

    #Remeshing
    if remeshing:
        if (output_sensitive and F_invalid.shape[0] > 0):
            # we find the invalid faces
            F_invalid = np.unique(I[invalid_U])
            # We compute the face adjacency matrix
            TT = _face_adjacency(state.F)
            # We find the set of invalid faces and their neighbors
            F_invalid_neighbors = np.unique(TT[F_invalid,:].ravel())
            # also add the invalid faces
            # F_invalid_neighbors = np.unique(np.hstack((F_invalid_neighbors,F_invalid)))
            F_invalid_neighbors = _one_ring(state.F,F_invalid)
            # do another round of one ring
            F_invalid_neighbors = _one_ring(state.F,F_invalid_neighbors)
            # We find the set of invalid vertices
            # F_active = F[F_invalid_neighbors,:]
            state.V_active, state.F_active, _, _ = remove_unreferenced(state.V,
                state.F[F_invalid_neighbors,:], return_maps=True)

            if (state.F_active.shape[0] < state.F.shape[0] and
                state.F_active.shape[0] > 0):
                # set of inactive faces
                state.F_inactive = np.setdiff1d(
                    np.arange(state.F.shape[0]), F_invalid_neighbors)
                # set of inactive vertices
                state.V_inactive, state.F_inactive, _, _ = remove_unreferenced(
                    state.V, state.F[state.F_inactive,:],return_maps=True)
                # Remesh only the active part
                state.V_active, state.F_active = _remesh(
                    state.V_active, state.F_active, i=remesh_iterations,
                    h=state.h, project=True)
                # We merge the active and inactive parts
                state.V, state.F = _merge_meshes(state.V_active, state.F_active, state.V_inactive, state.F_inactive)

            else:
                state.V, state.F = _remesh(state.V, state.F,
                    i=remesh_iterations, h=state.h, project=True)
        else:
            state.V, state.F = _remesh(state.V, state.F,
                i=remesh_iterations, h=state.h, project = True,
                feature=feature_vertices)
            state.V_active = None
            state.F_active = None
            state.V_inactive = None
            state.F_inactive = None
    else:
        state.V_active = None
        state.F_active = None
        state.V_inactive = None
        state.F_inactive = None

    #Have we converged?
    state.its = state.its+1
    if state.its>=max_iter:
        return True
    else:
        return False