[question] Finding number of closed potential contours and area enclosed by them

Specific Domains Visualization
, , , ,
9

Here’s how you can get the line coordinates with Contours.jl:

using Contours


len = 50
x = range(-1, 1; length = len)
y = range(-1, 1; length = len)
four_wells(x, y) = x^4 - x^2 + y^4 - y^2
z = [four_wells(x[i], y[j]) for i in 1:len, j in 1:len]

# energy
height = -0.4

cls = contours(x, y, z, [height]) # compute the contours at manually specified levels
level = only(levels(cls)) # extract the only result (since we only asked for one level)

points = [l.vertices for l in lines(level)] # ::Vector{Vector{Tuple{Float64, Float64}}}

Now that we have the lines, we need to check that they’re closed. The nature of Contour.jl’s algorithm makes this simple - we only need to know whether the first and last points are the same.

closed_ring_points = filter(points) do pointvec
    pointvec[begin] == pointvec[end] # you can also make this approximate, if you want
end

Now that you have polygon-like constructs, it’s easy enough to solve for its area. There are multiple solutions, I’m using GeoInterface.jl and GeometryOps.jl here.

import GeoInterface as GI, GeometryOps as GO
# construct geometry from vectors of points
rings = GI.LinearRing.(closed_ring_points)
polygons = [GI.Polygon([ring]) for ring in rings]

areas = GO.area.(polygons) # 4x ~0.1694

# if you want to be efficient / nonallocating
total_area = sum(GO.area, polygons) # 0.6775172177537117

Hope that helps! If you are worried about intersecting polygons, you can also use GeometryOps.jl to fix that:

mp = GI.MultiPolygon(polygons)

fixed_mp = GO.fix(mp; corrections = [GO.UnionIntersectingPolygons()])
2 Likes