Improved interpolations

Project.toml

@@ -1,10 +1,9 @@
 name = "TopoPlots"
 uuid = "2bdbdf9c-dbd8-403f-947b-1a4e0dd41a7a"
 authors = ["Benedikt Ehinger", "Simon Danisch", "Beacon Biosignals"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
-Delaunay = "07eb4e4e-0c6d-46ef-bc4e-83d5e5d860a9"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -17,7 +16,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 julia = "1"
-Delaunay = "1.2"
 Dierckx = "0.5"
 GeometryBasics = "0.4"
 Makie = "0.17.8"

docs/src/functions.md

@@ -4,7 +4,3 @@
 TopoPlots.enclosing_geometry
 TopoPlots.labels2positions
 ```
-
-```@docs
-TopoPlots.pad_boundary!
-```

docs/src/general.md

@@ -13,7 +13,6 @@ TopoPlots.topoplot
 The recipe supports different interpolation methods, namely:
 
 ```@docs
-TopoPlots.DelaunayMesh
 TopoPlots.ClaughTochter
 TopoPlots.SplineInterpolator
 TopoPlots.NullInterpolator
@@ -32,7 +31,7 @@ using TopoPlots, CairoMakie
 data, positions = TopoPlots.example_data()
 
 f = Figure(resolution=(1000, 1000))
-interpolators = [DelaunayMesh(), ClaughTochter(), SplineInterpolator()]
+interpolators = [ClaughTochter(), SplineInterpolator()]
 data_slice = data[:, 360, 1]
 
 for (i, interpolation) in enumerate(interpolators)
@@ -47,27 +46,6 @@ end
 f
 ```
 
-## Interactive exploration
-
-`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))
-s = Slider(f[:, 1], range=1:size(data, 2), startvalue=351)
-data_obs = map(s.value) do idx
-    data[:, idx, 1]
-end
-TopoPlots.topoplot(
-    f[2, 1],
-    data_obs, positions,
-    interpolation=DelaunayMesh(),
-    labels = string.(1:length(positions)),
-    colorrange=(-1, 1),
-    colormap=:viridis,
-    axis=(title="delaunay mesh",aspect=DataAspect(),))
-f
-```
-
 ## Different geometry
 
 The bounding geometry pads the input data with more points in the form of the geometry.

src/TopoPlots.jl

@@ -2,7 +2,6 @@ module TopoPlots
 
 using Makie
 using SciPy
-using Delaunay
 using Dierckx
 using LinearAlgebra
 using Statistics
@@ -28,6 +27,6 @@ include("core-recipe.jl")
 include("eeg.jl")
 
 # Interpolators
-export ClaughTochter, SplineInterpolator, DelaunayMesh
+export ClaughTochter, SplineInterpolator
 
 end

src/core-recipe.jl

@@ -64,17 +64,37 @@ macro plot_or_defaults(var, defaults)
     return :(plot_or_defaults($(esc(var)), $(esc(defaults)), $(QuoteNode(var))))
 end
 
+get_bounding_rect(rect::Rect) = rect
+
+function get_bounding_rect(circ::Circle)
+    xmin, ymin = minimum(circ)
+    xmax, ymax = maximum(circ)
+    Rect2f(xmin, ymin, xmax - xmin, ymax - ymin)
+end
+
+function get_bounded(rect, positions, data, padding, pad_value)
+    bb_xy = decompose(Point2f, rect)
+    bb_v = fill(pad_value, 4)
+    padding ≥ 0 || error("cannot deal with negative padding")
+    allunique(positions) || error("positions must be unique")
+    if padding == 0
+        i = findall(∈(positions), bb_xy)
+        deleteat!(bb_xy, i)
+        deleteat!(bb_v, i)
+    end
+    append!(bb_xy, positions)
+    append!(bb_v, data)
+    return bb_xy, bb_v
+end
+
 function Makie.plot!(p::TopoPlot)
     npositions = Observable(0; ignore_equal_values=true)
-    geometry = lift(enclosing_geometry, p.bounding_geometry, p.positions, p.padding; ignore_equal_values=true)
-    p.geometry = geometry # store geometry in plot object, so others can access it
     # positions changes with with data together since it gets into convert_arguments
     positions = lift(identity, p.positions; ignore_equal_values=true)
-    padded_position = lift(positions, geometry, p.resolution; ignore_equal_values=true) do positions, geometry, resolution
-        points_padded = append!(copy(positions), decompose(Point2f, geometry))
-        npositions[] = length(points_padded)
-        return points_padded
-    end
+    geometry = lift(enclosing_geometry, p.bounding_geometry, positions, p.padding; ignore_equal_values=true)
+    p.geometry = geometry # store geometry in plot object, so others can access it
+    bounding_box = lift(get_bounding_rect, geometry)
+    bounded = lift(get_bounded, bounding_box, positions, p.data, p.padding, p.pad_value)
 
     xg = Observable(LinRange(0f0, 1f0, p.resolution[][1]); ignore_equal_values=true)
     yg = Observable(LinRange(0f0, 1f0, p.resolution[][2]); ignore_equal_values=true)
@@ -88,22 +108,22 @@ function Makie.plot!(p::TopoPlot)
     end
     notify(p.resolution) # trigger above (we really need `update=true` for onany)
 
-    padded_data = lift(pad_data, p.data, npositions, p.pad_value)
-
-    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, padded_position)
-        mesh!(p, m, color=padded_data, colorrange=p.colorrange, colormap=p.colormap, shading=false)
-    else
-        data = lift(p.interpolation, xg, yg, padded_position, padded_data) do interpolation, xg, yg, points, data
-            return interpolation(xg, yg, points, data)
-        end
-        heatmap!(p, xg, yg, data, colormap=p.colormap, colorrange=p.colorrange, interpolate=true)
-        contours = to_value(p.contours)
-        attributes = @plot_or_defaults contours Attributes(color=(:black, 0.5), linestyle=:dot, levels=6)
-        if !isnothing(attributes)
-            contour!(p, xg, yg, data; attributes...)
+    data = lift(p.interpolation, xg, yg, bounded, geometry) do interpolation, xg, yg, (points, data), geometry
+        z = interpolation(xg, yg, points, data)
+        if geometry isa Circle
+            c = geometry.center
+            r = geometry.r
+            out = [norm(Point2f(x, y) - c) > r for x in xg, y in yg]
+            z[out] .= NaN
         end
+        return z
+    end
+
+    heatmap!(p, xg, yg, data, colormap=p.colormap, colorrange=p.colorrange, interpolate=true)
+    contours = to_value(p.contours)
+    attributes = @plot_or_defaults contours Attributes(color=(:black, 0.5), linestyle=:dot, levels=6)
+    if !isnothing(attributes)
+        contour!(p, xg, yg, data; attributes...)
     end
     label_scatter = to_value(p.label_scatter)
     attributes = @plot_or_defaults label_scatter Attributes(markersize=p.markersize, color=p.data, colormap=p.colormap, colorrange=p.colorrange, strokecolor=:black, strokewidth=1)
@@ -139,29 +159,3 @@ function enclosing_geometry(::Type{Rect}, positions, enlarge=0.0)
     mini = minimum(rect) .- ((padded_w .- w) ./ 2)
     return Rect2f(mini, padded_w)
 end
-
-"""
-    pad_boundary(::Type{Geometry}, positions, enlarge=0.2) where Geometry
-
-Adds new points to positions, adding the boundary from enclosing all positions with `Geometry`.
-See [`TopoPlots.enclosing_geometry`](@ref) for more details about the boundary.
-"""
-function pad_boundary!(::Type{Geometry}, positions, enlarge=0.2) where Geometry
-    c = enclosing_geometry(Geometry, positions, enlarge)
-    return append!(positions, decompose(Point2f, c))
-end
-
-function pad_data(data::AbstractVector, positions::AbstractVector, value::Number)
-    pad_data(data, length(positions), value)
-end
-
-function pad_data(data::AbstractVector, npositions::Integer, value::Number)
-    ndata = length(data)
-    if npositions == ndata
-        return data
-    elseif npositions < ndata
-        error("To pad the data for new positions, we need more positions than data points")
-    else
-        vcat(data, fill(value, npositions - ndata))
-    end
-end

src/interpolators.jl

@@ -58,28 +58,28 @@ function (sp::SplineInterpolator)(
     return evalgrid(spl, xrange, yrange)'
 end
 
-
+# currently, Delaunay doesn't work at all. The reason is that we can't interpolate using this mesh so we can't find the values of the data at the novel coordinates of the perimeter (geometry).
 # TODO how to properly integrade delauny with the interpolation interface,
 # if the actualy interpolation happens inside the plotting framework (or even on the GPU for (W)GLMakie).
-
-"""
-    DelaunayMesh()
-
-Creates a delaunay triangulation of the points and linearly interpolates between the vertices of the triangle.
-Really fast interpolation that happens on the GPU (for GLMakie), so optimal for exploring larger timeseries.
-
-!!! warning
-    `DelaunayMesh` won't allow you to add a contour plot to the topoplot.
-"""
-struct DelaunayMesh <: Interpolator
-end
-
-(::DelaunayMesh)(positions::AbstractVector{<: Point{2}}) = delaunay_mesh(positions)
-
-function delaunay_mesh(positions::AbstractVector{<: Point{2}})
-    m = delaunay(convert(Matrix{Float64}, hcat(first.(positions), last.(positions))))
-    return GeometryBasics.Mesh(Makie.to_vertices(m.points), Makie.to_triangles(m.simplices))
-end
+#
+# """
+#     DelaunayMesh()
+#
+# Creates a delaunay triangulation of the points and linearly interpolates between the vertices of the triangle.
+# Really fast interpolation that happens on the GPU (for GLMakie), so optimal for exploring larger timeseries.
+#
+# !!! warning
+#     `DelaunayMesh` won't allow you to add a contour plot to the topoplot.
+# """
+# struct DelaunayMesh <: Interpolator
+# end
+#
+# (::DelaunayMesh)(positions::AbstractVector{<: Point{2}}) = delaunay_mesh(positions)
+#
+# function delaunay_mesh(positions::AbstractVector{<: Point{2}})
+#     m = delaunay(convert(Matrix{Float64}, hcat(first.(positions), last.(positions))))
+#     return GeometryBasics.Mesh(Makie.to_vertices(m.points), Makie.to_triangles(m.simplices))
+# end
 
 
 #=

test/runtests.jl

@@ -7,7 +7,7 @@ data, positions = TopoPlots.example_data()
 
 begin
     f = Figure(resolution=(1000, 1000))
-    interpolators = [DelaunayMesh(), ClaughTochter(), SplineInterpolator()]
+    interpolators = [ClaughTochter(), SplineInterpolator()]
 
     s = Slider(f[:, 1], range=1:size(data, 2), startvalue=351)
     data_obs = map(s.value) do idx
@@ -24,30 +24,57 @@ begin
     @test_figure("all-interpolations", f)
 end
 
-
 begin # empty eeg topoplot
     f, ax, pl = TopoPlots.eeg_topoplot(1:length(TopoPlots.CHANNELS_10_20),TopoPlots.CHANNELS_10_20; interpolation=TopoPlots.NullInterpolator(),)
     @test_figure("nullInterpolator", f)
 end
 
-begin
-    f = Figure(resolution=(1000, 1000))
-    s = Slider(f[:, 1], range=1:size(data, 2), startvalue=351)
-    data_obs = map(s.value) do idx
-        data[:, idx, 1]
-    end
-    TopoPlots.topoplot(
-        f[2, 1],
-        data_obs, positions,
-        interpolation=DelaunayMesh(),
-        labels = string.(1:length(positions)),
-        colorrange=(-1, 1),
-        colormap=[:red, :blue],
-        axis=(title="delaunay mesh", aspect=DataAspect(),))
-    display(f)
-    @test_figure("delaunay-with-slider", f)
+@testset "peaks" begin
+    # 4 coordinates with one peak
+    positions = Point2f[(-1, 0), (0, -1), (1, 0), (0, 1)]
+    i = 1
+    peak_xy = positions[i]
+    data = zeros(length(positions))
+    data[i] = 1
+    fig = topoplot(data, positions)
+    # tighten the limits so that the limits of the axis and the data will match
+    tightlimits!(fig.axis)
+
+    # retrieve the interpolated data
+    m = fig.plot.plots[].color[]
+    # get the limits of the axes and data
+    rect = fig.axis.targetlimits[]
+    minx, miny = minimum(rect)
+    maxx, maxy = maximum(rect)
+    # recreate the coordinates of the data
+    x = range(minx, maxx, length=size(m, 1))
+    y = range(miny, maxy, length=size(m, 2))
+    xys = Point2f.(x, y')
+
+    # find the highest point
+    _, i = findmax(x -> isnan(x) ? -Inf : x, m)
+    xy = xys[i]
+    @test isapprox(xy, peak_xy; atol=0.02)
 end
 
+# begin
+#     f = Figure(resolution=(1000, 1000))
+#     s = Slider(f[:, 1], range=1:size(data, 2), startvalue=351)
+#     data_obs = map(s.value) do idx
+#         data[:, idx, 1]
+#     end
+#     TopoPlots.topoplot(
+#         f[2, 1],
+#         data_obs, positions,
+#         interpolation=DelaunayMesh(),
+#         labels = string.(1:length(positions)),
+#         colorrange=(-1, 1),
+#         colormap=[:red, :blue],
+#         axis=(title="delaunay mesh", aspect=DataAspect(),))
+#     display(f)
+#     @test_figure("delaunay-with-slider", f)
+# end
+
 begin
     f, ax, pl = TopoPlots.topoplot(
         data[:, 340, 1], positions,