Flooding from the dark: sublevel sets

val image = Image(pixels)
val dark = Persistence(image)     // sublevel sets: the darkest pixels first
dark.size                         // 221 bars

Noise makes 221 bars, nearly all of them short. Because we know the noise is at most 0.1, a sensible cut-off is "persists for more than 0.3", three times the noise amplitude. (With point clouds the scale of the data suggested a cut-off; here the noise level of the camera does. dark.significant(), the command-line and MATLAB default of 1% of the range, would keep 170.)

dark.longerThan(0.3)
// [0.000, Infinity)   [0.001, 1.002)   in degree 0
// [0.062, 1.099)      [0.064, 0.900)   in degree 1

Four features survive, and each one is something you can point at in the picture. Starting from the darkest pixels:

• [0, Infinity) in degree 0: the dark background, which exists from the start and is never swallowed.
• [0, 1.00) in degree 0: the dark disc inside the ring. It is a separate piece of darkness, walled in by the ring, until the flood reaches the ring's brightness (about 1.0). Then the ring's own pixels begin to enter, the disc is joined to the background, and its bar ends.
• [0.06, 1.10) in degree 1: a hole. The dark background, once it has formed a connected band around the ring, encircles it, a loop that stays until the ring itself is flooded (at 1.1, the brightest pixel of the ring).
• [0.06, 0.90) in degree 1: the same, around the blob, which is flooded at 0.9.

The death values of the two loops, 1.10 and 0.90, are the brightness of the ring and of the blob: sublevel sets read the brightness of obstacles. Everything else is noise.

Flooding from the bright: superlevel sets

Run the threshold the other way and you track the bright structures instead. image.superlevel is the same image filtered from the bright end; internally the image is negated, so the engine's "smallest first" becomes "brightest first", and the numbers you read back are negative brightnesses (a bar from -1.10 is something born at brightness 1.10):

val bright = Persistence(image.superlevel)
bright.longerThan(0.3)
// [-1.099, Infinity)   [-0.900, -0.063)   in degree 0
// [-1.004, -0.001)                        in degree 1

• [-1.10, Infinity) in degree 0: the ring, a bright piece born at its brightest pixel, never merged into anything.
• [-0.90, -0.06) in degree 0: the blob, born at brightness 0.90 and absorbed into the ring's piece when the threshold falls to the background level, 0.06.
• [-1.00, -0.00) in degree 1: the ring's loop. It is born at brightness 1.00, when the weakest pixel of the ring closes the circle, and dies at 0 when the dark disc in the middle has been flooded as well.

Compare the two views. The bright view says "a ring and a blob, and the ring has a hole"; the dark view says "a background with two obstacles in it, one of which walls in a pocket of darkness". These are the same facts, seen from either side. Which one is easier to read depends on whether the objects you care about are brighter or darker than their surroundings.

A faster engine

Persistence runs a general engine, the same as for point clouds, on all 3,249 cells. Images are so regular that there is a much faster approach exploiting the grid, the fast-cubical engine (see the user guide). Using an engine directly means choosing the coefficient field yourself:

given Double is Field = Field.DoubleApproximated(1e-9)
val fastBars = FastCubicalHomologyEngine[Double]().persistentHomology(CubicalImage.fromGrayscale2D(pixels))
// the same bars as dark.bars

On a 28 × 28 image you will not feel the difference. On a photograph of a few megapixels you will.

Where to go next

Images come from files more often than from code: the input/output page lists the image and grid formats TDA4j reads (Perseus and DIPHA grids, and the image loaders in CubicalImage), and the cubical complexes page covers three-dimensional volumes.

The whole script

import scala.language.experimental.modularity
import org.appliedtopology.tda4j.*

val random = new java.util.Random(5L)
val size = 28
val pixels = Array.tabulate(size, size) { (row, column) =>
  val distanceToRingCentre = math.hypot(row - 12.0, column - 12.0)
  val distanceToBlobCentre = math.hypot(row - 23.0, column - 23.0)
  val signal =
    if distanceToRingCentre >= 6.0 && distanceToRingCentre <= 8.0 then 1.0
    else if distanceToBlobCentre <= 2.5 then 0.8
    else 0.0
  signal + 0.1 * random.nextDouble()
}

val image = Image(pixels)
val dark = Persistence(image)
val bright = Persistence(image.superlevel)

val stream = CubicalImage.fromGrayscale2D(pixels)   // the cubical complex itself
given Double is Field = Field.DoubleApproximated(1e-9)
val fastBars = FastCubicalHomologyEngine[Double]().persistentHomology(stream)

javaaddpath('target/out/jvm/scala-3.9.0/tda4j/tda4j-0.5.0-SNAPSHOT-assembly.jar');   % from the repository root, after sbt assembly
import org.appliedtopology.tda4j.matlab.*;

% The same image as the Scala script, drawn from the same Java random number generator
random = java.util.Random(int64(5));
pixels = zeros(28, 28);
for row = 0:27
    for column = 0:27
        distanceToRingCentre = hypot(row - 12, column - 12);
        distanceToBlobCentre = hypot(row - 23, column - 23);
        if distanceToRingCentre >= 6 && distanceToRingCentre <= 8
            signal = 1.0;
        elseif distanceToBlobCentre <= 2.5
            signal = 0.8;
        else
            signal = 0.0;
        end
        pixels(row + 1, column + 1) = signal + 0.1 * random.nextDouble();
    end
end

% minPersistence 0.3 keeps the bars longer than 0.3, which is the Scala script's 'significant'
dark = TDA4j.computeFromImage(pixels, {'sublevel', 'true', 'minPersistence', '0.3'});
dark.toArray()          % four rows: [0 0 Inf], [0 0.001 1.002], [1 0.062 1.099], [1 0.064 0.900]

% With sublevel false the facade works on the negated image: births and deaths are NEGATED intensities
bright = TDA4j.computeFromImage(pixels, {'sublevel', 'false', 'minPersistence', '0.3'});
bright.toArray()        % three rows: [0 -1.099 Inf], [0 -0.900 -0.063], [1 -1.004 -0.001]

% The fast engine gives the same bars
fast = TDA4j.computeFromImage(pixels, {'sublevel', 'true', 'engine', 'fast-cubical', 'minPersistence', '0.3'});
isequal(sortrows(fast.toArray()), sortrows(dark.toArray()))      % true
size(dark.toArrayUnfiltered(), 1)                                 % 221 bars in the full barcode

The facade takes an image as a matrix and has the engine choice as an option. sublevel false filters the negated image, which is why the bright features have negative births and deaths. MatlabTabsSpec makes these same calls from Scala and checks every number quoted in the MATLAB tab; MATLAB itself is not run by the test suite.