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	Matrix of SDF sample positions. The sdf will be sampled at these points	None
S	(k,dim) nummpy double array, optional (default	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	Matrix of mesh vertices active during iteration	None
F_active	(ma,dim) numpy int array, optional (default	Matrix of triangle indices into V_active of the mesh active during iteration	None
V_inactive	(na,dim) numpy double array, optional (default	Matrix of mesh vertices inactive during iteration	None
F_inactive	(ma,dim) numpy int array, optional (default	Matrix of triangle indices into V_inactive of the mesh inactive during iteration	None
rng	numpy random generator, optional (default	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	Matrix of mesh vertices for the last known converged mesh.	None
F_last_converged	(ml,dim) numpy int array, optional (default	Matrix of triangle indices into V_last_converged for the last known converged mesh.	None
U_batch	(kb,dim) numpy double array, optional (default	Matrix of SDF sample positions as used during the last batching operation.	None
S_batch	(kb,dim) nummpy double array, optional (default	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

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	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	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			required
Vr	(nr,dim) numpy double array	Matrix of mesh vertices of the reconstructed triangle mesh	required
Fr	(mr,dim) numpy int array	Matrix of triangle indices into Vr of the reconstructed mesh	required
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.	required

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

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			required
converged	bool	Whether the method has converged after this iteration or not.	required

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

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=V[bd,:])
    else:
        state.V = sp.sparse.linalg.spsolve(Q,b)

    # catching flow singularities so we fail gracefully
    if np.any((np.isnan(state.V))):
        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 = np.vstack((state.V_active, state.V_inactive))
                state.F = np.vstack((state.F_active,
                    state.F_inactive + state.V_active.shape[0]))
                # We remove the duplicate vertices
                state.V,_,_,state.F = remove_duplicate_vertices(
                    state.V,faces=state.F,
                    epsilon=np.sqrt(np.finfo(state.V.dtype).eps))
            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