Ref. [1] has made available a large collection of labelled images of BECs containing 0, 1 or multiple dark solitons. The authors of that study trained a convolutional neural network to distinguish these three classes of images. Single dark solitons could be detected with an accuracy of about 90\%. On 2014 MacBook Pro each image took 2.3s to classify using the neural network. Our aim here is to develop a persistent-homology based approach with superior performance (both in accuracy and run-time). Fast and accurate machine learning-based identification of solitons would enable high-throughput experimental study of their dynamics by eliminating the time-consuming and at times error-prone manual identification of the solitons in the experimentally-obtained images.
First we import the python packages needed in the following analysis,
import numpy as np
import matplotlib.pyplot as plt
import scipy
from persim import plot_diagrams
from ripser import ripser, lower_star_img
import os
from timeit import default_timer as timer
Next let's import and inspect one of the sample images provided by Ref. [1].
# read the filenames contained in the image folder "data"
image_filenames = os.listdir('data')
img_test = np.load('data/' + image_filenames[10], allow_pickle=True)
img_data = img_test.item()['data'][:,:,0]
img_label = img_test.item()['label']
#(ny, nx) = img_data.shape # image dimensions, needed later
plt.imshow(img_data, cmap='gray', origin='lower')
plt.title("Test image, label = " + str(img_label))
plt.show()
This trial image contains a clearly-visible intensity dip near the centre of the BEC. We will start by analyzing the simpler case of a 1D cut of the condensate intensity profile.
y = 75
img_cut = img_data[y]
plt.plot(img_cut)
plt.title("Intensity cut along y = " + str(y))
plt.xlabel("x")
plt.show()
In addition to the prominent dip at the dark soliton core (near x = 70), the intensity profile has many other local minima. We can use the sublevel set filtration to reliably identify the most persistent minima and discard those created due to noise-induced local fluctuations of the intensity. We use the example provided in the online documentation of Ref. [2] to define the filtration,
def time_series_filtration(x):
N=x.shape[0]
# Add edges between adjacent points in the time series, with the "distance"
# along the edge equal to the max value of the points it connects
I = np.arange(N-1)
J = np.arange(1, N)
V = np.maximum(x[0:-1], x[1::])
# Add vertex birth times along the diagonal of the distance matrix
I = np.concatenate((I, np.arange(N)))
J = np.concatenate((J, np.arange(N)))
V = np.concatenate((V, x))
#Create the sparse distance matrix
D = scipy.sparse.coo_matrix((V, (I, J)), shape=(N, N)).tocsr()
return ripser(D, maxdim=0, distance_matrix=True)['dgms'][0]
With this filtration, points are added to the simplicial complex whenever their intensity exceeds the filtration value, and only neighbouring points are connected by vertices. The resulting simplicial complex only exhibits 0-dimensional features (cluters). Features are born whenever the filtration value crosses a local minimum of intensity profile, and clusters merge (resulting in death of a feature) whenever the filtration value crosses a local intensity maxima. We compute the persistence diagram of our 1D cut,
rips_result = time_series_filtration(img_cut)
plt.figure(figsize=(10, 5))
plt.subplot(121)
plot_diagrams(rips_result)
plt.title("0-D Persistence Diagram")
plt.subplot(122)
plot_diagrams(rips_result, lifetime=True)
plt.title("0-D Feature Lifetimes")
plt.tight_layout()
plt.show()
There are many features near the diagonal - these have low persistence and are noise-induced. In the upper left there are three highly persistent features, associated with three distinct clusters of low intensity. To see where these clusters are located, we now plot the birth and death times of these features (with lifetimes > 0.1).
thresh = 0.1
lifetimes = rips_result[:,1] - rips_result[:,0]
birth_times = rips_result[lifetimes > thresh,0]
death_times = rips_result[lifetimes > thresh,1]
death_times = death_times[ death_times < np.inf]
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.plot(img_cut)
ax = plt.gca()
ax.set_yticks(birth_times)
ax.set_xticks([])
ax.tick_params(labelleft=False)
plt.grid(linewidth=1, linestyle='--')
plt.title("Feature birth times")
plt.subplot(122)
plt.plot(img_cut)
ax = plt.gca()
ax.set_yticks(death_times)
ax.set_xticks([])
ax.tick_params(labelleft=False)
plt.grid(linewidth=1, linestyle='--')
plt.title("Feature death times")
plt.show()
Thus, the dark soliton can be identified by three persistent clusters of low intensity (to the left of the soliton, at the soliton core, and to the right of the soliton). However, the persistence threshold must be chosen to distinguish dark solitons from noise-induced intensity dips. This might be difficult if the image is highly noisy or has a low contrast.
Let's now repeat the above procedure for each row in the image, highlighting the birth pixels of the persistent clusters.
def find_persistent_minima(image, threshold):
(ny, nx) = image.shape # image dimensions, needed later
Nmin = 20 # estimate of the maximum number of persistent minima per row
minima = np.zeros((ny*Nmin,2)) # array for storing the persistent minima coordinates
count = 0 #counter for number of persistent minima identified
for y in range(0,ny,1):
subimg = image[y]
rips_result = time_series_filtration(subimg)
lifetimes = rips_result[:,1] - rips_result[:,0]
idxs = np.nonzero(lifetimes > threshold)[0]
for i in idxs:
index = np.nonzero(np.abs(subimg-rips_result[i,0])<1e-6)[0][0] #find the pixel position corresponding to the feature birth time
minima[count,:] = [y, index]
count = count + 1
#exclude birth positions at edges of the image (spurious)
minima = minima[(minima[:,0] > 0.1) & (minima[:,0] < nx-0.1) & (minima[:,1] > 0.1) & (minima[:,1] < ny-0.1)]
return minima
min1 = find_persistent_minima(img_data,0.05)
min2 = find_persistent_minima(img_data,0.1)
min3 = find_persistent_minima(img_data,0.15)
plt.figure(figsize=(10, 5))
plt.subplot(221)
plt.imshow(img_data,origin='lower',cmap='gray')
plt.title("Test Image")
ax=plt.subplot(222)
plt.scatter(min1[:,1],min1[:,0],s=1)
plt.title("Persistent minima, threshold = 0.05")
ax.set_aspect(1)
ax=plt.subplot(223)
plt.scatter(min2[:,1],min2[:,0],s=1)
plt.title("Persistent minima, threshold = 0.1")
ax.set_aspect(1)
ax=plt.subplot(224)
plt.scatter(min3[:,1],min3[:,0],s=1)
plt.title("Persistent minima, threshold = 0.15")
ax.set_aspect(1)
plt.tight_layout()
plt.show()
If the persistence threshold is too low we end up catching too many spurious minima. If it is too high we fail to identify the soliton minimum in the low contrast parts of the image. The soliton minimum forms a continous line, whereas noise-induced minima are typically isolated. We can try further suppressing the noise by considering sublevel set filtrations of narrow 2D slices of the image, effectively averaging over rapid noise-induced intensity fluctuations.
def find_persistent_minima_sliced(image, threshold, slice_width):
(ny, nx) = image.shape
Nmin = 20 # estimate of the maximum number of persistent minima per row
minima = np.zeros((ny*Nmin,2)) # array for storing the persistent minima coordinates
count = 0 #counter for number of persistent minima identified
for y in range(0,ny,slice_width):
subimg = image[y:y+slice_width]
rips_result = lower_star_img(subimg)
lifetimes = rips_result[:,1] - rips_result[:,0]
idxs = np.nonzero(lifetimes > threshold)[0]
for i in idxs:
index = np.mod(np.nonzero(np.abs(np.ravel(subimg)-rips_result[i,0])<1e-6)[0][0],nx)
minima[count] = [y+round(slice_width/2),index]
count = count + 1
#exclude birth positions at edges of the image (spurious)
minima = minima[(minima[:,0] > 0.1) & (minima[:,0] < nx-0.1) & (minima[:,1] > 0.1) & (minima[:,1] < ny-0.1)]
return minima
min1 = find_persistent_minima_sliced(img_data,0.05,5)
min2 = find_persistent_minima_sliced(img_data,0.1,5)
min3 = find_persistent_minima_sliced(img_data,0.15,5)
plt.figure(figsize=(10, 5))
plt.subplot(221)
plt.imshow(img_data,origin='lower',cmap='gray')
plt.title("Test Image")
ax=plt.subplot(222)
plt.scatter(min1[:,1],min1[:,0],s=1)
plt.title("Persistent minima, threshold = 0.05")
ax.set_aspect(1)
ax=plt.subplot(223)
plt.scatter(min2[:,1],min2[:,0],s=1)
plt.title("Persistent minima, threshold = 0.1")
ax.set_aspect(1)
ax=plt.subplot(224)
plt.scatter(min3[:,1],min3[:,0],s=1)
plt.title("Persistent minima, threshold = 0.15")
ax.set_aspect(1)
plt.tight_layout()
plt.show()
Using slicing we significantly reduce the number of spurious points appearing at a given persistence threshold.
As a final step, we need to identify the dark soliton as dark line that traverses the entire soliton. We can also do this using persistent homology. First, we combine our point cloud of local minima with points randomly-sampled from the edge of the condensate.
def sample_edge_points(image,Npoints):
(ny, nx) = image.shape
cutoff_intensity = 0.25 # this is the pixel value set for pixels lying outside of the condensate region
edge = np.zeros(((nx+ny)*2,2))
count = 0
for y in range(ny):
im = image[y,:]
x0 = np.argmax(im != cutoff_intensity)
if x0>1:
edge[count,:] = [y,x0-1]
count = count + 1
im2 = im[::-1]
x2 = len(im2) - np.argmax(im2 != cutoff_intensity) -1
if x2 < nx-1:
edge[count,:] = [y,x2]
count = count+1
for x in range(nx):
im = image[:,x]
y0 = np.argmax(im != cutoff_intensity)
if y0>1:
edge[count,:] = [y0-1,x]
count = count + 1
im2 = im[::-1]
y2 = len(im2) - np.argmax(im2 != cutoff_intensity) -1
if y2 < ny-1:
edge[count,:] = [y2,x]
count = count + 1
edge = edge[(edge[:,0] > 0.1) & (edge[:,0] < ny-0.1) & (edge[:,1] > 0.1) & (edge[:,1] < nx-0.1)]
edge = edge[ np.random.randint(0,edge.shape[0],Npoints),:]
return edge
edge = sample_edge_points(img_data,100)
plt.scatter(edge[:,1],edge[:,0])
plt.title("Points on edge of condensate")
plt.show()
minima = find_persistent_minima_sliced(img_data,0.1,6)
edge = sample_edge_points(img_data,200)
point_cloud = np.concatenate((minima, edge))
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(img_data,origin='lower',cmap='gray')
plt.title("Test Image")
plt.subplot(122)
plt.scatter(point_cloud[:,1],point_cloud[:,0],s=1)
plt.title("Persistent minima + edge of condensate")
plt.tight_layout()
plt.show()
The number of dark soliton stripe can be counted by computing the number of persistent one-dimensional features (holes) of this 2D point cloud. If there are no stripes, then the cloud will have a single persistent hole (corresponding to the loop of points along the edge of the condensate). Each dark soliton stripe added to the condensate will add a new one-dimensional feature to the point cloud.
dgms = ripser(point_cloud)['dgms']
plot_diagrams(dgms, lifetime=True)
Here we there are two one-dimensional features (H_1) with a long lifetime. Thus, the single dark soliton was successfully identified. The separation between the H_1 features with long lifetimes and the noise-induced features with short lifetimes can be used as a measure of the uncertainty of the classification. Let us use the average lifetime of the H_0 features (excluding the feature with infinite lifetime) as a persistence threshold for counting the number of dark solitons.
def count_dark_solitons(image,threshold=0.1,slice_size=6,sample_size=100):
minima = find_persistent_minima_sliced(image,threshold,slice_size)
edge = sample_edge_points(img_data,sample_size)
point_cloud = np.concatenate((minima, edge))
dgms = ripser(point_cloud)['dgms']
H0_lifetimes = dgms[0][:,1] - dgms[0][:,0]
H1_lifetimes = dgms[1][:,1] - dgms[1][:,0]
H1_threshold = sum(H0_lifetimes[1:-1]) / H0_lifetimes.shape[0]
return sum(H1_lifetimes > H1_threshold)-1
Let us now repeat the procedure over a subset of the data.
Nsamples = 1000
test_results = np.zeros((Nsamples,2))
start = timer()
for n, i in enumerate(np.random.randint(0,len(image_filenames),Nsamples)):
img_test = np.load('data/' + image_filenames[i], allow_pickle=True)
img_data = img_test.item()['data'][:,:,0]
img_label = img_test.item()['label']
N_solitons = count_dark_solitons(img_data)
if N_solitons < 2:
img_result = N_solitons
else:
img_result = 2 # data uses 2 to label all images which do not have either 0 or 1 soliton
test_results[n,:] = [img_label, img_result]
end = timer()
print("Elapsed time = ", end - start)
print("Time per image = ", (end-start)/Nsamples)
Elapsed time = 84.84474258497357 Time per image = 0.08484474258497357
In comparison, the convolutional neural network approach required ~2.4s per image. So our approach is more than 20 times faster, even without any performance optimization. What about the accuracy?
agreed = test_results[:,0] == test_results[:,1]
disagreed = ~agreed
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.hist(test_results[agreed,0])
plt.title("Distribution of sucessfully classified images")
plt.xlabel("Image label")
plt.ylabel("Number of images")
plt.subplot(122)
plt.hist(test_results[disagreed,1])
plt.title("Distribution of misclassified images")
plt.xlabel("Image label")
plt.ylabel("Number of images")
plt.show()
print("Success probability = ", sum(agreed)/Nsamples)
Success probability = 0.65
The accuracy is already quite good given that we chose the model hyperparameters (threshold for persistent minima, slice sizes, persistence threshold for the H_1 features) using a single good quality image. The next steps will be to optimize these hyper-parameters further, either via direct search or by studying the behaviour of the mis-classified images, and to add some estimation of uncertainty in the algorithm's prediction.
References¶
[1] S. Guo, A. R. Fritsch, C. Greenberg, I. B. Spielman, J.P. Zwolak, Machine-learning enhanced dark soliton detection in Bose-Einstein condensates, arXiv:2101.05404
[2] C. Tralie, N. Saul, R. Bar-On, Ripser.py: A Lean Persistent Homology Libraray for Python, Journal of Open Source Software 3, 925 (2018); https://ripser.scikit-tda.org/en/latest/
No comments:
Post a Comment