maltezfaria · maltezfaria · Mar 16, 2026 · Mar 16, 2026 · Mar 16, 2026 · Mar 16, 2026
diff --git a/Project.toml b/Project.toml
@@ -7,20 +7,19 @@ authors = ["Luiz M. Faria and Nicolas Lebbe"]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
 MMG_jll = "86086c02-e288-5929-a127-40944b0018b7"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MarchingCubes = "299715c1-40a9-479a-aaf9-4a633d36f717"
-NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 
 [extensions]
 MMGSurfaceExt = ["MMG_jll", "MarchingCubes"]
 MMGVolumeExt = "MMG_jll"
 MakieExt = "Makie"
-ReinitializationExt = ["NearestNeighbors"]
 
 [compat]
 DiffResults = "1.1.0"

diff --git a/docs/make.jl b/docs/make.jl
@@ -16,7 +16,7 @@ const numbered_pages = [
 ]
 
 modules = [LevelSetMethods]
-for extension in [:MakieExt, :MMGSurfaceExt, :MMGVolumeExt, :ReinitializationExt]
+for extension in [:MakieExt, :MMGSurfaceExt, :MMGVolumeExt]
     ext = Base.get_extension(LevelSetMethods, extension)
     if isnothing(ext)
         @warn "extension $extension not loaded"
@@ -40,8 +40,9 @@ makedocs(;
         "interpolation.md",
         "time-integrators.md",
         "boundary-conditions.md",
+        "reinitialization.md",
         "Extensions" =>
-            ["extension-makie.md", "extension-mmg.md", "extension-reinitialization.md"],
+            ["extension-makie.md", "extension-mmg.md"],
         "Examples" => ["example-zalesak.md", "example-shape-optim.md"],
         numbered_pages,
     ),

diff --git a/docs/src/example-mass-conservation.md b/docs/src/example-mass-conservation.md
@@ -30,7 +30,7 @@ w = 0.2
 disk = LevelSetMethods.circle(grid; center, radius)
 rec = LevelSetMethods.rectangle(grid; center = center .- (0, radius), width = (w, h))
 ϕ₀ = setdiff(disk, rec)
-V₀ = volume(ϕ₀)
+V₀ = LevelSetMethods.volume(ϕ₀)
 tf = 2π   # one full revolution
 nframes = 50
 timestamps = range(0, tf; length = nframes + 1)
@@ -48,7 +48,7 @@ function track_volume(integrator)
     V[1] = V₀
     for (i, t) in enumerate(timestamps[2:end])
         integrate!(eq, t)
-        V[i + 1] = volume(current_state(eq))
+        V[i + 1] = LevelSetMethods.volume(current_state(eq))
     end
     return (V .- V₀) ./ V₀
 end

diff --git a/docs/src/extension-makie.md b/docs/src/extension-makie.md
@@ -49,7 +49,7 @@ Once again, you can manually customize the options by calling the `volume` funct
 
 ```@example volume3D
 with_theme(theme) do
-  volume(ϕ; algorithm = :iso, isovalue = 0.5)
+  Makie.volume(ϕ; algorithm = :iso, isovalue = 0.5)
 end
 ```
 

diff --git a/docs/src/extension-reinitialization.md b/docs/src/extension-reinitialization.md
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -113,8 +113,6 @@ optional dependencies:
 
 - **[Makie](https://docs.makie.org)**: Loading a `Makie` backend (like `GLMakie` or `CairoMakie`)
   enables plotting recipes for level sets and equations. See [Makie extension](@ref extension-makie).
-- **[NearestNeighbors](https://github.com/KristofferC/NearestNeighbors.jl)**: Loading
-  `NearestNeighbors` enables high-order Newton-based reinitialization. See [Reinitialization](@ref extension-reinitialization).
 - **[MMG](https://github.com/JuliaBinaryWrappers/MMG_jll.jl.git)**: Loading `MMG_jll` and `MarchingCubes`
   enables exporting level sets as volume or surface meshes. See [MMG extension](@ref extension-mmg).
 

diff --git a/docs/src/interpolation.md b/docs/src/interpolation.md
@@ -21,7 +21,7 @@ bc = ntuple(_ -> (NeumannGradientBC(), NeumannGradientBC()), 2)
 itp = interpolate(ϕ) # cubic interpolation by default (order=3)
 ```
 
-The returned object is a [`PiecewisePolynomialInterpolation`](@ref), which is callable and
+The returned object is a [`PiecewisePolynomialInterpolation`](@ref LevelSetMethods.PiecewisePolynomialInterpolation), which is callable and
 efficient. Once constructed, the interpolant can be used to evaluate the level-set function
 anywhere inside (and even slightly outside, using boundary conditions) the grid:
 

diff --git a/docs/src/reinitialization.md b/docs/src/reinitialization.md
@@ -0,0 +1,118 @@
+# [Reinitialization](@id reinitialization-newton)
+
+Reinitialization transforms a level set function into a signed distance function without
+moving the zero level set. This is useful when the level set has been distorted by
+advection or other terms, since many numerical schemes assume ``|\nabla\phi| \approx 1``.
+
+Two approaches are available:
+
+- **[`EikonalReinitializationTerm`](@ref)**: a PDE-based approach that evolves the level
+  set under the equation ``\phi_t + \text{sign}(\phi)(|\nabla\phi| - 1) = 0`` using the
+  same time-stepping infrastructure as other level-set terms. See the [Reinitialization
+  term](@ref reinitialization) section.
+- **[`NewtonReinitializer`](@ref)** (recommended): a geometry-based approach that samples
+  the interface, builds a KD-tree, and computes the exact signed distance to the interface
+  using Newton's closest-point method. It is applied between time steps and converges in a
+  single pass.
+
+## Usage
+
+Call `reinitialize!` on a `LevelSet` to reinitialize it in place:
+
+```@example reinit
+using LevelSetMethods
+using GLMakie
+
+grid = CartesianGrid((-1, -1), (1, 1), (100, 100))
+sdf = LevelSet(x -> sqrt(x[1]^2 + x[2]^2) - 0.5, grid)
+ϕ = LevelSet(x -> x[1]^2 + x[2]^2 - 0.5^2, grid)
+LevelSetMethods.set_makie_theme!()
+fig = Figure(; size = (800, 300))
+ax1 = Axis(fig[1, 1]; title = "Signed Distance")
+ax2 = Axis(fig[1, 2]; title = "Before Reinitialization", ylabel = "", yticklabelsvisible = false)
+ax3 = Axis(fig[1, 3]; title = "After Reinitialization", ylabel = "", yticklabelsvisible = false)
+
+contour!(ax1, sdf; levels = [0.25, 0, 0.5], labels = true, labelsize = 14)
+contour!(ax2, ϕ; levels = [0.25, 0, 0.5], labels = true, labelsize = 14)
+
+reinitialize!(ϕ, NewtonReinitializer())
+contour!(ax3, ϕ; levels = [0.25, 0, 0.5], labels = true, labelsize = 14)
+fig
+```
+
+You can verify that the reinitialized level set is a signed distance function:
+
+```@example reinit
+max_er = maximum(eachindex(grid)) do i
+  abs(ϕ[i] - sdf[i])
+end
+println("Maximum error after reinitialization: $max_er")
+```
+
+## Automatic reinitialization in `LevelSetEquation`
+
+Pass `reinit` to [`LevelSetEquation`](@ref) to reinitialize automatically every `n` steps:
+
+```julia
+eq = LevelSetEquation(;
+    terms  = (AdvectionTerm(𝐮),),
+    levelset = ϕ,
+    bc     = PeriodicBC(),
+    reinit = 5,          # reinitialize every 5 time steps
+)
+```
+
+The integer shorthand creates a `NewtonReinitializer(; reinit_freq = n)` with default
+settings. For full control over the algorithm parameters, pass a `NewtonReinitializer`
+directly:
+
+```julia
+eq = LevelSetEquation(;
+    terms    = (AdvectionTerm(𝐮),),
+    levelset = ϕ,
+    bc       = PeriodicBC(),
+    reinit   = NewtonReinitializer(; reinit_freq = 5, upsample = 4),
+)
+```
+
+You can also reinitialize the equation's current state manually at any time:
+
+```julia
+reinitialize!(eq)
+```
+
+## `NewtonSDF`: a reusable signed distance function
+
+`LevelSetMethods.NewtonSDF` wraps the same closest-point algorithm in a callable object,
+letting you evaluate the signed distance at arbitrary points without modifying the level
+set in place. This is useful when you need a signed distance function as an ingredient in
+a larger computation (e.g. to measure distances or to build an extension velocity):
+
+```@example reinit
+using StaticArrays
+sdf_obj = LevelSetMethods.NewtonSDF(ϕ; upsample = 8)
+# query the signed distance at arbitrary points
+sdf_obj(SVector(0.0, 0.0))   # distance from origin to the circle
+```
+
+The interface sample points used to build the KD-tree can be retrieved with
+`LevelSetMethods.get_sample_points`:
+
+```@example reinit
+pts = LevelSetMethods.get_sample_points(sdf_obj)
+println("$(length(pts)) interface sample points")
+```
+
+When the underlying level set changes, use `LevelSetMethods.update!` to rebuild the
+KD-tree, reusing the same interpolation order and tolerances:
+
+```@example reinit
+ϕ2 = LevelSet(x -> x[1]^2 + x[2]^2 - 0.3^2, grid)  # smaller circle
+LevelSetMethods.update!(sdf_obj, ϕ2)
+sdf_obj(SVector(0.3, 0.0))   # should be ≈ 0
+```
+
+!!! warning "Thread safety"
+    `NewtonSDF` uses internal mutable buffers in its interpolant and is **not thread-safe**.
+    If you need to evaluate it from multiple threads, use `deepcopy(sdf_obj)` to give
+    each thread its own independent copy.
diff --git a/docs/src/terms.md b/docs/src/terms.md
@@ -253,7 +253,7 @@ fig
 We will now evolve the level-set using the reinitialization term:
 
 ```@example reinitialization-term
-eq = LevelSetEquation(; terms = (ReinitializationTerm(),), levelset = deepcopy(ϕ), bc = PeriodicBC())
+eq = LevelSetEquation(; terms = (EikonalReinitializationTerm(),), levelset = deepcopy(ϕ), bc = PeriodicBC())
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.25, 0.5, 0.75])
     integrate!(eq, t)
@@ -271,10 +271,10 @@ Alternatively, you can use a modified reinitialization term that applies the sig
   \phi_t + \text{sign}(\phi_0) \left( |\nabla \phi| - 1 \right) = 0
 ```
 
-To enable this behavior, simply pass a `LevelSet` object to the `ReinitializationTerm`:
+To enable this behavior, simply pass a `LevelSet` object to the `EikonalReinitializationTerm`:
 
 ```@example reinitialization-term
-eq = LevelSetEquation(; terms = (ReinitializationTerm(ϕ),), levelset = deepcopy(ϕ), bc = PeriodicBC())
+eq = LevelSetEquation(; terms = (EikonalReinitializationTerm(ϕ),), levelset = deepcopy(ϕ), bc = PeriodicBC())
 fig = Figure(; size = (1200, 300))
 for (n,t) in enumerate([0.0, 0.25, 0.5, 0.75])
     integrate!(eq, t)
@@ -285,3 +285,10 @@ fig
 ```
 
 The outcome closely matches that of the previous approach.
+
+!!! tip "Consider Newton reinitialization instead"
+    `EikonalReinitializationTerm` requires many pseudo-time steps to propagate corrections
+    away from the interface and can cause spurious mass loss. For most use cases,
+    [`NewtonReinitializer`](@ref) is a better choice: it samples the interface, builds a
+    KD-tree, and computes the exact signed distance in a single pass. See the
+    [Reinitialization](@ref reinitialization-newton) section for details.