Sublevel Image Filtrations

In this example, we will demonstrate sublevel persistent homology on images. We will use the Horizon telescope picture of a black hole. We will use the 768×768 pixel image. Computing the same result with a larger image should not be a problem, but the plots would make this page less responsive.

Start by loading the required packages and the image.

using Ripserer
using Plots
using Images
blackhole = load(joinpath(@__DIR__, "../assets/blackhole768px.jpg"))

plot(blackhole, size=(800, 800))

Convert the image to a matrix of floats.

blackhole_grayscale = Float64.(Gray.(blackhole))

768×768 Array{Float64,2}:
 0.027451  0.027451  0.027451  0.027451  …  0.0313725  0.0313725  0.0313725
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451  …  0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.0235294  0.0235294  0.0235294
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 ⋮                                       ⋱  ⋮                     
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451  …  0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451  …  0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451
 0.027451  0.027451  0.027451  0.027451     0.027451   0.027451   0.027451

And plot it as a heatmap. Note that the image is flipped because matrices start at the top and images start at the bottom.

heatmap(blackhole_grayscale, aspect_ratio=1, size=(800, 800))

Compute the persistent homology.

result_min = ripserer(CubicalFiltration(blackhole_grayscale))

2-element Array{PersistenceDiagram{PersistenceInterval{Float64,Nothing}},1}:
 640-element 0-dimensional PersistenceDiagram
 474-element 1-dimensional PersistenceDiagram

We plot the diagram with the persistence argument. This plots persistence vs birth instead of death vs birth.

plot(result_min, persistence=true)

We notice there is a lot of noise along the diagonal, but four intervals stand out. We remove the noise by supplying the cutoff argument to ripserer. This removes all intervals with persistence strictly smaller than cutoff.

result_min = ripserer(CubicalFiltration(blackhole_grayscale), cutoff=0.1)

2-element Array{PersistenceDiagram{PersistenceInterval{Float64,Nothing}},1}:
 2-element 0-dimensional PersistenceDiagram
 2-element 1-dimensional PersistenceDiagram

plot(result_min, persistence=true)

Now that we know we have a smaller number of intervals, we can compute the representatives.

result_min = ripserer(CubicalFiltration(blackhole_grayscale), cutoff=0.1, representatives=true)

2-element Array{PersistenceDiagram,1}:
 2-element 0-dimensional PersistenceDiagram
 2-element 1-dimensional PersistenceDiagram

We separate the finite and infinite intervals. The finite one corresponds to the local minimum in the middle of the image and the infinite one corresponds to the global minimum.

finite_interval = only(filter(isfinite, result_min[1]))
infinite_interval = only(filter(!isfinite, result_min[1]))

Plot the location of the minima on the image. We have to use the threshold argument to only plot the simplices with diameter lower than or equal to the birth of the interval.

heatmap(blackhole_grayscale, aspect_ratio=1, size=(800, 800))
plot!(finite_interval, blackhole_grayscale,
      color=2,
      markershape=:d,
      markerstrokecolor=2,
      threshold=birth(finite_interval),
      label="local minimum")

plot!(infinite_interval, blackhole_grayscale,
      color=3,
      markershape=:d,
      markerstrokecolor=3,
      threshold=birth(infinite_interval),
      label="global minimum")

Note that there are several points for each minimum because several pixels have the same value.

To plot the region of the local minimum, we skip the threshold and choose a small marker size. We will not plot the global minimum region because it encompasses the whole image.

heatmap(blackhole_grayscale, aspect_ratio=1, size=(800, 800))
plot!(finite_interval, blackhole_grayscale,
      markersize=1,
      color=1,
      markerstrokecolor=1,
      label="local minimum region")