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

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.

`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

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