Implement interpolation methods from NaturalNeighbours.jl

Project.toml

@@ -1,7 +1,7 @@
 name = "TopoPlots"
 uuid = "2bdbdf9c-dbd8-403f-947b-1a4e0dd41a7a"
 authors = ["Benedikt Ehinger", "Simon Danisch", "Beacon Biosignals, Inc."]
-version = "0.1.7"
+version = "0.1.8"
 
 [deps]
 CloughTocher2DInterpolation = "b70b374f-000b-463f-88dc-37030f004bd0"
@@ -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"
@@ -23,7 +24,8 @@ Dierckx = "0.5"
 GeometryBasics = "0.4"
 InteractiveUtils = "1"
 LinearAlgebra = "1"
-Makie = "0.17.8, 0.18, 0.19, 0.20"
+Makie = "0.17.8, 0.18, 0.19, 0.20, 0.21"
+NaturalNeighbours = "1"
 Parameters = "0.12"
 PrecompileTools = "1"
 ScatteredInterpolation = "0.3.6"

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"

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"
     ],
 )
 

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]

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.

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
+```

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

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