Polygonize is double polygonizing in some places

[asinghvi17]

Not sure exactly where, still have to narrow this down. But take a look at this example and the video, note that there are a few areas that have two polygons over them...

using DimensionalData

import GeometryOps as GO, GeoInterface as GI

using CairoMakie, GeoInterfaceMakie


# generate a dataset of moving waves

xvals = X(LinRange(-9, 9, 49))
yvals = Y(LinRange(-9, 9, 49))

tvals = Ti(LinRange(0, 6pi, 100))

waves = @d sin.(hypot.(xvals, yvals) .+ tvals)


# animate

ti_obs = Observable(1)
waves_obs = @lift waves[Ti($ti_obs)]
f, a, p = heatmap(waves_obs)

record(f, "wave_tracking.mp4", 1:length(val(tvals)), framerate=24) do i
    ti_obs[] = i
end

# now, track the wave

# let's say that any value above 0.8 is a part of a "wave"
# then, we can polygonize each wave

polygons_per_time_idx = map(eachslice(waves, dims=Ti)) do wave
    xs = val(val(dims(wave, X)))
    ys = val(val(dims(wave, Y)))
    return GO.polygonize(xs, ys, wave.data .> 0.8).geom # this returns a GI.MultiPolygon, which we decompose into a list of polygons
end

# let's animate this too

ti_obs = Observable(1)
polygons_obs = @lift polygons_per_time_idx[Ti($ti_obs)]

f, a, p = poly(polygons_obs; axis = (; aspect = DataAspect()))

record(f, "wave_tracking_polygons2.mp4", 1:length(val(tvals)), framerate=24) do i
    ti_obs[] = i
end

wave_tracking.mp4

wave_tracking_polygons2.mp4

cc @rafaqz

I'll make this an MWE later.

[rafaqz]

Probably making a wrong turn somewhere that I didn't notice because it rasterizes flat again. It's better than missing pixels I guess!

[asinghvi17]

huh, this is actually kind of a Makie issue? Take a look at this...it's a timeseries of intersecting polygons across timesteps

wave_tracking_cc1.mp4

If you take the wave cube at Ti=1 and polygonize it, then one of the polygons has 3 rings, which causes the doubling up. Maybe one of those rings needs to be rejected if it's contained within the other interior ring?

[asinghvi17]

For the record, here's how you convert these polygons into a graph (though it's a bit unrelated to the issue):

n_polygons = sum(length, polygons_per_time_idx)

using Graphs, SimpleWeightedGraphs

graph = SimpleWeightedDiGraph(n_polygons)

# Step through each time index, and add an edge for polygons which intersect each other
current_polygon_idx = 1
for ti in 1:(length(polygons_per_time_idx)-1)
    current_polys = polygons_per_time_idx[Ti=ti]
    next_polys = polygons_per_time_idx[Ti=ti+1]

    for (i, current_poly) in enumerate(current_polys)
        for (j, next_poly) in enumerate(next_polys)
            if GO.intersects(current_poly, next_poly)
                add_edge!(
                    graph, 
                    current_polygon_idx,  # the index of the current polygon
                    current_polygon_idx + length(current_polys) - i + j,  # the index of the next polygon
                    GO.area(GO.intersection(current_poly, next_poly; target = GI.PolygonTrait))
                )
            end
        end
        current_polygon_idx += 1
    end
end

# We can isolate individual wave systems by finding the connected components of the graph
ccs = Graphs.connected_components(graph)

# Let's plot one such wave system
flattened_polygons = reduce(vcat, polygons_per_time_idx)

f, a, p = poly(flattened_polygons[ccs[1][1]]; axis = (; aspect = DataAspect()))
record(f, "wave_tracking_cc1.mp4", ccs[1], framerate=5) do i
    p.args[1][] = flattened_polygons[i]
end


using GraphMakie

f, a, p = graphplot(graph; layout = NetworkLayout.Shell(; nlist = Graphs.connected_components(graph)))

[rafaqz]

Actually looking at the video better it seems like a nested rings/holes problem. Maybe we don't handle that properly.