1 introduction

This section presents low-star and sub-set filter shapes on image data. This filter shape allows us to use the local minimum or local maximum as the initial time and the saddle point as the death time in the 0-dimensional continuation graph. The following are some necessary libraries:

import numpy as np
import matplotlib.pyplot as plt
import scipy
import PIL
from scipy import ndimage
from persim import plot_diagrams
from ripser import ripser, lower_star_img

The follow-up will define a function for constructing a 0-dimensional low-star filter shape on the image, which constructs a sparse distance matrix, where each pixel in the image is a vertex, and each vertex is connected to eight spatial neighbors, unless is the case on the border. Boundary weights are considered to be the maximum of the two pixel values ​​they connect (hence the "low star").

Although the number of pixels is large, the time complexity of this function is linearly close to the number of pixels, because there is only a finite number of bounds on the pixel scale, and only 0-dimensional homotopy is computed.

2 Gaussian spots

Here is an image with Gaussian blobs: Three negative Gaussians are placed in the image with local minima at -3, -2, and -1:

ts = np.linspace(-1, 1, 100)
x1 = np.exp(-ts**2/(0.1**2))
ts -= 0.4
x2 = np.exp(-ts**2/(0.1**2))
img = -x1[None, :]*x1[:, None] - 2*x1[None, :]*x2[:, None] - 3*x2[None, :]*x2[:, None]
plt.imshow(img)
plt.show()

The output is as follows:

Obtain its 0-dimensional continuation graph:

dgm = lower_star_img(img)

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(img)
plt.colorbar()
plt.title("Test Image")
plt.subplot(122)
plot_diagrams(dgm)
plt.title("0-D Persistence Diagram")
plt.tight_layout()
plt.show()

The output is as follows:

In this example, the three points in the graph correspond to three Gaussian blobs, each initially at its own minimum. Two of them perish around 0, meeting other classes at the saddle point, while the one initial at -3 absorbs other classes and lives forever due to the "elder rule".

3 Cell Biology Images

The schematic image is as follows:

When converting this image to a grayscale image of [0,255], where 0 and 255 correspond to dark and light respectively, we can see that the internal brightness of the cells is very high, and each cell is at a saddle point somewhere on the boundary meet, the saddle point is closer.

In this case, it is assumed that there is a large persistence of local maxima within each cell; that is, we perform superordered sets or ascending star filtering. To make the code work for local maxima but not local minima, one can simply provide the lower star filter function with the negative of the image:

cells_original = (plt.imread("Cells.png") * 255).astype(np.uint8)
cells_grey = np.asarray(PIL.Image.fromarray(cells_original).convert('L'))

plt.subplot(121)
plt.title(cells_original.shape)
plt.imshow(cells_original)
plt.axis('off')
plt.subplot(122)
plt.title(cells_grey.shape)
plt.imshow(cells_grey, cmap='gray')
plt.axis('off')
plt.show()

The output is as follows:

Draw a 0-dimensional continuation graph:

dgm = lower_star_img(-cells_grey)

plt.figure(figsize=(6, 6))
plot_diagrams(dgm, lifetime=True)
plt.show()

The output is as follows:

Now choose a persistence threshold above which a point is considered associated with a cell.

Ripser.py currently does not return representations of 0-dimensional homotopy classes, but can do some workarounds by adding a small amount of uniform noise to each pixel. This makes each pixel have a unique value, so one can simply find the pixel whose value is equal to the initial time of the class you are looking for. Before doing this, local averaging is performed to drive the representation of the maximum closer to the center of the cell:

smoothed = ndimage.uniform_filter(cells_grey.astype(np.float64), size=10)
smoothed += 0.01 * np.random.randn(*smoothed.shape)

plt.figure(figsize=(10, 5))
plt.subplot(121)
im = plt.imshow(cells_grey, cmap='gray')
plt.colorbar(im, fraction=0.03)

plt.subplot(122)
im = plt.imshow(smoothed, cmap='gray')
plt.colorbar(im, fraction=0.03)

plt.tight_layout()
plt.show()

The output is as follows:

Draw the 0D continuation graph again:

dgm = lower_star_img(-smoothed)
plot_diagrams(dgm, lifetime=True)
plt.show()

The output is as follows:

We'll look at a cutoff point and look at all points with a lifetime greater than 70. The following shows each 0D pixel representation highlighted in the original image:

thresh = 70
idxs = np.arange(dgm.shape[0])
idxs = idxs[np.abs(dgm[:, 1] - dgm[:, 0]) > thresh]

plt.figure(figsize=(8, 5))
plt.imshow(cells_original)

X, Y = np.meshgrid(np.arange(smoothed.shape[1]), np.arange(smoothed.shape[0]))
X = X.flatten()
Y = Y.flatten()
for idx in idxs:
    bidx = np.argmin(np.abs(smoothed + dgm[idx, 0]))
    plt.scatter(X[bidx], Y[bidx], 20, 'k')
plt.axis('off')

plt.show()

The output is as follows: