isoband.jl

using RCall
using AbstractPlotting
using AbstractPlotting.MakieLayout
using GLMakie
using EarCut
using AbstractPlotting.ColorSchemes
using AbstractPlotting.GeometryBasics
using PolygonOps


# R"""install.packages("isoband")"""
R"""library("isoband")"""

xs = LinRange(0, 10, 100)
ys = LinRange(0, 10, 100)
zs = [sin(x) * cos(y) for x in xs, y in ys]


function isobands(xs, ys, zs, low, high)
    result = first(values(rcopy(R"""isobands($xs, $ys, $zs, $low, $high)""")))
    points = Point2f0.(result[:x], result[:y])
    ids = result[:id]
    (points = points, ids = ids)
end

points, ids = isobands(xs, ys, zs, -0.8, 0.8)
unique(ids)


function group_polys(points, ids)

    polys = [points[ids .== i] for i in unique(ids)]
    npolys = length(polys)

    polys_lastdouble = [push!(p, first(p)) for p in polys]

    # this matrix stores whether poly i is contained in j
    # because the marching squares algorithm won't give us any
    # intersecting or overlapping polys, it should be enough to
    # check if a single point is contained, saving some computation time
    containment_matrix = [
        p1 != p2 &&
        inpolygon(first(p1), p2) == 1
        for p1 in polys_lastdouble, p2 in polys_lastdouble]

    unclassified_polyindices = collect(1:size(containment_matrix, 1))
    # @show unclassified_polyindices

    # each group has first an outer polygon, and then its holes
    groups = Vector{eltype(polys)}[]
    # a dict that maps index in `polys` to index in `groups` for outer polys
    outerindex_groupdict = Dict{Int, Int}()

    # all polys have to be classified
    while !isempty(unclassified_polyindices)
        to_keep = ones(Bool, length(unclassified_polyindices))

        # println("containment matrix")
        # display(containment_matrix)

        # go over unclassifieds and find outer polygons in the remaining containment matrix
        for (ii, i) in enumerate(unclassified_polyindices)
            # an outer polygon is not inside any other polygon of the matrix
            if sum(containment_matrix[ii, :]) == 0
                # an outer polygon
                # println(i, " is an outer polygon")
                push!(groups, [polys_lastdouble[i]])
                outerindex_groupdict[i] = length(groups)
                # delete this poly from further rounds
                to_keep[ii] = false
            end
        end

        # go over unclassifieds and find hole polygons
        for (ii, i) in enumerate(unclassified_polyindices)
            # the hole polygons can only be in one polygon from the current group
            # if they are in more than one, they are "inner outer" or inner hole polys
            # and will be handled in one of the following passes
            if sum(containment_matrix[ii, :]) == 1
                outerpolyindex_of_unclassified = findfirst(containment_matrix[ii, :])
                outerpolyindex = unclassified_polyindices[outerpolyindex_of_unclassified]
                # a hole
                # println(i, " is an inner polygon of ", outerpolyindex)
                groupindex = outerindex_groupdict[outerpolyindex]
                push!(groups[groupindex], polys_lastdouble[i])
                # delete this poly from further rounds
                to_keep[ii] = false
            end
        end
    
        unclassified_polyindices = unclassified_polyindices[to_keep]
        containment_matrix = containment_matrix[to_keep, to_keep]
    end
    groups
end




function contourf!(ax, xs, ys, zs, nlevels)

    mi, ma = extrema(zs)
    levels = collect(LinRange(mi, ma, nlevels+1))

    scene = Scene()

    for (lo, hi, i) in zip(levels[1:end-1], levels[2:end], 1:nlevels)
        points, ids = isobands(xs, ys, zs, lo, hi)

        groups = group_polys(points, ids)

        for group in groups
            faces = EarCut.triangulate(group)
            mesh!(ax, vcat(group...), faces, color = ColorSchemes.get(ColorSchemes.viridis, (i - 1) / (nlevels - 1)), shading = false)
        end
    end

    scene
end

##
using CairoMakie
CairoMakie.activate!(type = "png")

##
using RDatasets
RDatasets.datasets("datasets")

volcano = dataset("datasets", "volcano")
volcano = Array(volcano)

volcano_nan = Float64.(volcano)
volcano_nan[rand(length(volcano_nan)) .< 0.01] .= NaN

##
scene, layout = layoutscene()
ax = layout[1, 1] = LAxis(scene)

contourf!(ax, 1:size(volcano, 2), 1:size(volcano, 1), volcano, 10)
tightlimits!(ax)
ax.aspect = 1

scene