Implement interpolation methods from NaturalNeighbours.jl (#46)

* Explicitly mention how to do regridding using the TopoPlots API * Bump version + add NaturalNeighbours compat * Implement NaturalNeighbours.jl interpolation as `NaturalNeighboursMethod` similar to how `ScatteredInterpolationMethod` currently exists. * Document NaturalNeighboursMethod * Fix Makie warning about `shading=false` being invalid Replaced by `NoShading` * Describe the method for scattered and nn interpolation in the docs * Associate `lift`s with the plot in the topoplot recipe This makes for cleaner GC/windups when deleting a plot - the lift observer functions are known to the plot object, and can be removed immediately instead of relying on the GC to detect that those variables are no longer in use. * Implement gridding * Add an "interpolator reference image" page to the docs * Add NaturalNeighbours.jl to the docs Project * Add a brief summary of the advantage of natural neighbour interpolation * Warn against using DelaunayMesh in CairoMakie * Fix error in doc build by sampling without replacement * Fix typos / omissions * finally * fix more doc errors * ugh * Add compat for Makie v0.21 * Address code review * Restore the option for the derivative calculation to also be parallelized
MakieOrg · Jun 14, 2024 · 7f81e8d · 7f81e8d
1 parent 84122ef
commit 7f81e8d
Show file tree

Hide file tree

Showing 9 changed files with 235 additions and 24 deletions.
diff --git a/Project.toml b/Project.toml
@@ -11,6 +11,7 @@ GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+NaturalNeighbours = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 ScatteredInterpolation = "3f865c0f-6dca-5f4d-999b-29fe1e7e3c92"
@@ -24,6 +25,7 @@ GeometryBasics = "0.4"
 InteractiveUtils = "1"
 LinearAlgebra = "1"
 Makie = "0.17.8, 0.18, 0.19, 0.20, 0.21"
+NaturalNeighbours = "1"
 Parameters = "0.12"
 PrecompileTools = "1"
 ScatteredInterpolation = "0.3.6"

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -2,6 +2,7 @@
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NaturalNeighbours = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 ScatteredInterpolation = "3f865c0f-6dca-5f4d-999b-29fe1e7e3c92"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TopoPlots = "2bdbdf9c-dbd8-403f-947b-1a4e0dd41a7a"

diff --git a/docs/make.jl b/docs/make.jl
@@ -16,7 +16,8 @@ makedocs(;
         "Home" => "index.md",
         "General TopoPlots" => "general.md",
         "EEG" => "eeg.md",
-        "Function reference" => "functions.md"
+        "Function reference" => "functions.md",
+        "Interpolator reference images" => "interpolator_reference.md"
     ],
 )
 

diff --git a/docs/src/general.md b/docs/src/general.md
@@ -9,34 +9,47 @@ TopoPlots.topoplot
 
 ## Interpolation
 
+TopoPlots provides access to interpolators from several different Julia packages through its [`TopoPlots.Interpolator`](@ref) interface.
+
+They can be accessed via plotting, or directly by calling the instantiated interpolator object as is shown below, namely with the arguments `(::Interpolator)(xrange::LinRange, yrange::LinRange, positions::AbstractVector{<: Point{2}}, data::AbstractVector{<:Number})`.  This is similar to using things like Matlab's `regrid` function.  You can find more details in the [Interpolation](@ref) section.
+
 The recipe supports different interpolation methods, namely:
 
 ```@docs
 TopoPlots.DelaunayMesh
 TopoPlots.CloughTocher
 TopoPlots.SplineInterpolator
 TopoPlots.ScatteredInterpolationMethod
+TopoPlots.NaturalNeighboursMethod
 TopoPlots.NullInterpolator
 ```
-One can define your own interpolation by subtyping:
+
+You can define your own interpolation by subtyping:
 
 ```@docs
 TopoPlots.Interpolator
 ```
 
+and making your interpolator `SomeInterpolator` callable with the signature 
+```julia        
+(::SomeInterpolator)(xrange::LinRange, yrange::LinRange, positions::AbstractVector{<: Point{2}}, data::AbstractVector{<:Number}; mask=nothing)
+```
+
 The different interpolation schemes look quite different:
 
 ```@example 1
-using TopoPlots, CairoMakie, ScatteredInterpolation
+using TopoPlots, CairoMakie, ScatteredInterpolation, NaturalNeighbours
 
 data, positions = TopoPlots.example_data()
 
-f = Figure(resolution=(1000, 1000))
+f = Figure(resolution=(1000, 1250))
 
 interpolators = [
     DelaunayMesh() CloughTocher();
     SplineInterpolator() NullInterpolator();
-    ScatteredInterpolationMethod(ThinPlate()) ScatteredInterpolationMethod(Shepard(3))]
+    ScatteredInterpolationMethod(ThinPlate()) ScatteredInterpolationMethod(Shepard(3));
+    NaturalNeighboursMethod(Sibson(1)) NaturalNeighboursMethod(Triangle());
+    ]
 
 data_slice = data[:, 360, 1]
 
@@ -61,6 +74,9 @@ for idx in CartesianIndices(interpolators)
         axis=(type=Axis, title="$(typeof(interpolation))()",aspect=DataAspect(),))
 
    ax.title = ("$(typeof(interpolation))() - $(round(t, digits=2))s")
+   if interpolation isa Union{NaturalNeighboursMethod, ScatteredInterpolationMethod}
+       ax.title = "$(typeof(interpolation))($(typeof(interpolation.method))) - $(round(t, digits=2))s"
+   end
 end
 f
 ```
@@ -100,7 +116,7 @@ f
 `DelaunayMesh` is best suited for interactive data exploration, which can be done quite easily with Makie's native UI and observable framework:
 
 ```@example 1
-f = Figure(resolution=(1000, 1000))
+f = Figure(resolution=(1000, 1250))
 s = Slider(f[:, 1], range=1:size(data, 2), startvalue=351)
 data_obs = map(s.value) do idx
     data[:, idx, 1]

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -21,3 +21,5 @@ f
 Find more documentation for `topoplot` in [Recipe for General TopoPlots](@ref).
 
 It also contains some more convenience methods for EEG data, which is explained in [EEG Topoplots](@ref).
+
+You can also use TopoPlots' interpolators as a simple interface to regrid irregular data.  See [Interpolation](@ref) for more details.
diff --git a/docs/src/interpolator_reference.md b/docs/src/interpolator_reference.md
@@ -0,0 +1,104 @@
+# Interpolator reference
+
+This file contains reference figures showing the output of each interpolator available in TopoPlots, as well as timings for them. 
+
+It is a more comprehensive version of the plot in [Interpolation](@ref).
+
+### Example data
+
+```@example 1
+
+using TopoPlots, CairoMakie, ScatteredInterpolation, NaturalNeighbours
+
+data, positions = TopoPlots.example_data()
+
+f = Figure(size=(1000, 1500))
+
+interpolators = [
+    SplineInterpolator() NullInterpolator() DelaunayMesh();
+    CloughTocher() ScatteredInterpolationMethod(ThinPlate()) ScatteredInterpolationMethod(Shepard(3));
+    ScatteredInterpolationMethod(Multiquadratic()) ScatteredInterpolationMethod(InverseMultiquadratic()) ScatteredInterpolationMethod(Gaussian());
+    NaturalNeighboursMethod(Hiyoshi(2)) NaturalNeighboursMethod(Sibson()) NaturalNeighboursMethod(Laplace());
+    NaturalNeighboursMethod(Farin()) NaturalNeighboursMethod(Sibson(1)) NaturalNeighboursMethod(Nearest());
+    ]
+
+data_slice = data[:, 360, 1]
+
+for idx in CartesianIndices(interpolators)
+    interpolation = interpolators[idx]
+
+    # precompile to get accurate measurements
+    TopoPlots.topoplot(
+        data_slice, positions;
+        contours=true, interpolation=interpolation,
+        labels = string.(1:length(positions)), colorrange=(-1, 1),
+        label_scatter=(markersize=10,),
+        axis=(type=Axis, title="...", aspect=DataAspect(),))
+
+    # measure time, to give an idea of what speed to expect from the different interpolators
+    t = @elapsed ax, pl = TopoPlots.topoplot(
+        f[Tuple(idx)...], data_slice, positions;
+        contours=true,
+        interpolation=interpolation,
+        labels = string.(1:length(positions)), colorrange=(-1, 1),
+        label_scatter=(markersize=10,),
+        axis=(type=Axis, title="$(typeof(interpolation))()",aspect=DataAspect(),))
+
+   ax.title = ("$(typeof(interpolation))() - $(round(t, digits=2))s")
+   if interpolation isa Union{NaturalNeighboursMethod, ScatteredInterpolationMethod}
+       ax.title = "$(typeof(interpolation))() - $(round(t, digits=2))s"
+       ax.subtitle = string(typeof(interpolation.method))
+   end
+end
+f
+```
+
+### Randomly sampled function
+
+```@example 1
+using TopoPlots, CairoMakie, ScatteredInterpolation, NaturalNeighbours
+using TopoPlots.Makie.Random: randsubseq
+
+data = Makie.peaks(100)
+sampling_points = randsubseq(CartesianIndices(data), 0.011)
+data_slice = data[sampling_points]
+positions = Point2f.(Tuple.(sampling_points))
+
+interpolators = [
+    SplineInterpolator(; smoothing = 7) NullInterpolator() DelaunayMesh();
+    CloughTocher() ScatteredInterpolationMethod(ThinPlate()) ScatteredInterpolationMethod(Shepard(3));
+    ScatteredInterpolationMethod(Multiquadratic()) ScatteredInterpolationMethod(InverseMultiquadratic()) ScatteredInterpolationMethod(Gaussian());
+    NaturalNeighboursMethod(Hiyoshi(2)) NaturalNeighboursMethod(Sibson()) NaturalNeighboursMethod(Laplace());
+    NaturalNeighboursMethod(Farin()) NaturalNeighboursMethod(Sibson(1)) NaturalNeighboursMethod(Nearest());
+    ]
+
+f = Figure(; size = (1000, 1500))
+
+for idx in CartesianIndices(interpolators)
+    interpolation = interpolators[idx]
+
+    # precompile to get accurate measurements
+    TopoPlots.topoplot(
+        data_slice, positions;
+        contours=true, interpolation=interpolation,
+        labels = string.(1:length(positions)), colorrange=(-1, 1),
+        label_scatter=(markersize=10,),
+        axis=(type=Axis, title="...", aspect=DataAspect(),))
+
+    # measure time, to give an idea of what speed to expect from the different interpolators
+    t = @elapsed ax, pl = TopoPlots.topoplot(
+        f[Tuple(idx)...], data_slice, positions;
+        contours=true,
+        interpolation=interpolation,
+        labels = string.(1:length(positions)), colorrange=(-1, 1),
+        label_scatter=(markersize=10,),
+        axis=(type=Axis, title="$(typeof(interpolation))()",aspect=DataAspect(),))
+
+   ax.title = ("$(typeof(interpolation))() - $(round(t, digits=2))s")
+   if interpolation isa Union{NaturalNeighboursMethod, ScatteredInterpolationMethod}
+       ax.title = "$(typeof(interpolation))() - $(round(t, digits=2))s"
+       ax.subtitle = string(typeof(interpolation.method))
+   end
+end
+f
+```
diff --git a/src/TopoPlots.jl b/src/TopoPlots.jl
@@ -7,12 +7,13 @@ using GeometryBasics
 using GeometryBasics: origin, radius
 using Parameters
 using InteractiveUtils # needed for subtypes
+using PrecompileTools
+# Import the various interpolation packages
 using Dierckx
 using ScatteredInterpolation
-using PrecompileTools
-
+using NaturalNeighbours # voronoi tessellation based methods
 using Delaunator # DelaunayMesh
-using CloughTocher2DInterpolation  # pure julia implementation
+using CloughTocher2DInterpolation  # pure julia implementation of the algorithm used by scipy etc.
 
 assetpath(files...) = normpath(joinpath(dirname(@__DIR__), "assets", files...))
 
@@ -50,7 +51,7 @@ include("core-recipe.jl")
 include("eeg.jl")
 
 # Interpolators
-export CloughTocher, SplineInterpolator, DelaunayMesh, NullInterpolator, ScatteredInterpolationMethod
+export CloughTocher, SplineInterpolator, DelaunayMesh, NullInterpolator, ScatteredInterpolationMethod, NaturalNeighboursMethod
 @deprecate ClaughTochter(args...; kwargs...) CloughTocher(args...; kwargs...) true
 # Extrapolators
 export GeomExtrapolation, NullExtrapolation

diff --git a/src/core-recipe.jl b/src/core-recipe.jl
@@ -72,13 +72,13 @@ function Makie.plot!(p::TopoPlot)
     npositions = Obs(0)
 
     # positions changes with with data together since it gets into convert_arguments
-    positions = lift(identity, p.positions; ignore_equal_values=true)
-    geometry = lift(enclosing_geometry, p.bounding_geometry, positions, p.enlarge; ignore_equal_values=true)
+    positions = lift(identity, p, p.positions; ignore_equal_values=true)
+    geometry = lift(enclosing_geometry, p, p.bounding_geometry, positions, p.enlarge; ignore_equal_values=true)
 
     xg = Obs(LinRange(0f0, 1f0, p.interp_resolution[][1]))
     yg = Obs(LinRange(0f0, 1f0, p.interp_resolution[][2]))
 
-    f = onany(geometry, p.interp_resolution) do geometry, interp_resolution
+    f = onany(p, geometry, p.interp_resolution) do geometry, interp_resolution
         (xmin, ymin), (xmax, ymax) = extrema(geometry)
         xg[] = LinRange(xmin, xmax, interp_resolution[1])
         yg[] = LinRange(ymin, ymax, interp_resolution[2])
@@ -89,11 +89,11 @@ function Makie.plot!(p::TopoPlot)
 
     p.geometry = geometry # store geometry in plot object, so others can access it
 
-    padded_pos_data_bb = lift(p.extrapolation, p.positions, p.data) do extrapolation, positions, data
+    padded_pos_data_bb = lift(p, p.extrapolation, p.positions, p.data) do extrapolation, positions, data
         return extrapolation(positions, data)
     end
 
-    colorrange = lift(p.data, p.colorrange) do data, crange
+    colorrange = lift(p, p.data, p.colorrange) do data, crange
         if crange isa Makie.Automatic
             return Makie.extrema_nan(data)
         else
@@ -103,15 +103,15 @@ function Makie.plot!(p::TopoPlot)
 
     if p.interpolation[] isa DelaunayMesh
         # TODO, delaunay works very differently from the other interpolators, so we can't switch interactively between them
-        m = lift(delaunay_mesh, p.positions)
-        mesh!(p, m, color=p.data, colorrange=colorrange, colormap=p.colormap, shading=false)
+        m = lift(delaunay_mesh, p, p.positions)
+        mesh!(p, m, color=p.data, colorrange=colorrange, colormap=p.colormap, shading=NoShading)
     else
-        mask = lift(xg,yg,geometry) do xg,yg,geometry
+        mask = lift(p, xg, yg, geometry) do xg,yg,geometry
             pts = Point2f.(xg' .* ones(length(yg)), ones(length(xg))' .* yg)
-            return in.(pts,Ref(geometry))
+            return in.(pts, Ref(geometry))
         end
 
-        data = lift(p.interpolation, xg, yg, padded_pos_data_bb,mask) do interpolation, xg, yg, (points, data, _, _),mask
+        data = lift(p, p.interpolation, xg, yg, padded_pos_data_bb,mask) do interpolation, xg, yg, (points, data, _, _),mask
             z = interpolation(xg, yg, points, data;mask=mask)
 #            z[mask] .= NaN
             return z