
Download supplementary movie 1 (‘nmeth.3036sv1.avi’) of Amat et al., 2014, entitled ‘SiMView imaging of Drosophila embryogenesis’.

Construct a function to read the .avi video using the moviepy library with optional rescaling of video frame dimensions through the resize parameter. For more information using moviepy please refer to https://zulko.github.io/moviepy/.
def read_movie(moviefile, resize=1.):
from moviepy.editor import VideoFileClip
from tqdm import tqdm
from skimage.transform import rescale
import numpy as np
vidframes = []
clip = VideoFileClip(moviefile)
count = 0
for frame in tqdm(clip.iter_frames()):
vidframes.append(np.uint8(rescale(frame, 1./resize, preserve_range=True)))
count+=1
vidframes = np.array(vidframes)
return vidframes
moviefile = '../Data/Videos/nmeth.3036sv1.avi'
# read in the movie, downsampling frames
by 4.
movie = read_movie(moviefile, resize=4.)
n_frame, n_rows, n_cols, n_channels =
movie.shape
print('Size
of video: (%d,%d,%d,%d)' %(n_frame,n_rows,n_cols,n_channels))

Set motion parameters as before and extract forward and backward superpixel tracks using dense superpixel tracking.
from MOSES.Optical_Flow_Tracking.superpixel_track import compute_grayscale_vid_superpixel_tracks_FB
# motion extraction
parameters.
optical_flow_params = dict(pyr_scale=0.5, levels=5, winsize=21, iterations=5, poly_n=5, poly_sigma=1.2, flags=0)
# number of superpixels
n_spixels = 1000
# extract forward and
backward tracks.
optflow, meantracks_F, meantracks_B =
compute_grayscale_vid_superpixel_tracks_FB(movie[:,:,:,1], optical_flow_params,
n_spixels, dense=True)

Cluster superpixels according to track similarity to group together similar tracks. This requires the construction of a feature matrix for each superpixel track.

Construct the track feature, f_{ i}for superpixel track i by concatenating the (x,y) coordinates at each time frame up to frame T such that f=[x_1,x_2,…,x_(T1),x_T,y_1,y_2,…,y_(T1),y_T ]. For this example, the final feature matrix, X is a matrix of 2628 x 1104.
import numpy as np
X = np.hstack([meantracks_F[:,:,0], meantracks_F[:,:,1]])

Reduce dimensionality of feature matrix, X using principal components analysis (PCA) to compress features and improve clustering.
from sklearn.decomposition import PCA
pca_model = PCA(n_components = 3, whiten=True, random_state=0)
X_pca = pca_model.fit_transform(X)

Cluster tracks using the PCA reduced features using a Gaussian mixture model (GMM) specifying 10 clusters. The result is given as a vector array, track_labels specifying using unique integer values the group that each superpixel track has been assigned to by the GMM model. Note colors will be different for different runs of the algorithm due to randomness in the setting of the GMM algorithm.
Note: Any other clustering algorithm operating on a matrix can be used such as hierarchical clustering, Kmeans clustering and density based clustering.
from sklearn import mixture
n_clusters = 10
gmm = mixture.GaussianMixture(n_components=n_clusters, covariance_type='full', random_state=0)
gmm.fit(X_pca)
# get the labels
track_labels = gmm.predict(X_pca)

Visualize the clustered superpixel tracks coloring each unique group with a different color (Figure 16).
Note: We confirmed differences in GMM clustering results between different OS and scikitlearn library versions with the same code. The reported figure was obtained for Windows 10, scikit
learn version 0.21.2.
from MOSES.Visualisation_Tools.track_plotting import plot_tracks
import seaborn as sns
import pylab as plt
# Generate colours for
each unique cluster.
cluster_colors = sns.color_palette('hls', n_clusters)
# overlay cluster tracks
and clustered superpixels
fig, ax = plt.subplots(nrows=1,ncols=3, figsize=(15,15))
ax[0].imshow(movie[0]); ax[0].grid('off');
ax[0].axis('off'); ax[0].set_title('Initial Points')
ax[1].imshow(movie[1]); ax[1].grid('off');
ax[1].axis('off'); ax[1].set_title('Final Points')
ax[2].imshow(movie[0]); ax[2].grid('off');
ax[2].axis('off'); ax[2].set_title('Clustered Tracks')
for ii, lab in enumerate(np.unique(track_labels)):
# plot coloured initial
points
ax[0].plot(meantracks_F[track_labels==lab,0,1],
meantracks_F[track_labels==lab,0,0], 'o', color=cluster_colors[ii], alpha=1)
# plot coloured final
points
ax[1].plot(meantracks_F[track_labels==lab,1,1],
meantracks_F[track_labels==lab,1,0], 'o', color=cluster_colors[ii], alpha=1)
# plot coloured tracks
plot_tracks(meantracks_F[track_labels==lab],
ax[2], color=cluster_colors[ii], lw=1.0, alpha=0.7)
plt.show()
Figure 16. Unsupervised clustering of superpixel tracks to group local motion during Drosophila embryogenesis. Initial (A) and final (B) centroid positions of superpixels plotted as a dot and colored by assigned group using Gaussian mixture model for clustering. Corresponding superpixel tracks colored by cluster (C).

Compute the motion saliency map of the forward and backward tracks to find motion sources and sinks (Figure 17).
from MOSES.Motion_Analysis.mesh_statistics_tools import compute_motion_saliency_map
from skimage.exposure import equalize_hist
from skimage.filters import Gaussian
# specify a large threshold
to capture longdistances.
dist_thresh = 20
spixel_size = meantracks[1,0,1]meantracks[1,0,0]
motion_saliency_F,
motion_saliency_spatial_time_F = compute_motion_saliency_map(meantracks, dist_thresh=dist_thresh, shape=movie.shape[1:1], max_frame=None, filt=1, filt_size=spixel_size)
motion_saliency_B,
motion_saliency_spatial_time_B = compute_motion_saliency_map(meantracks_B, dist_thresh=dist_thresh, shape=movie.shape[1:1], max_frame=None, filt=1, filt_size=spixel_size)
# smooth the discrete
looking motion saliency maps.
motion_saliency_F_smooth =
gaussian(motion_saliency_F, spixel_size/2.)
motion_saliency_B_smooth =
gaussian(motion_saliency_B, spixel_size/2.)
# visualise the computed
results.
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(15,15))
ax[0].imshow(movie[0], cmap='gray');
ax[0].grid('off'); ax[0].axis('off')
ax[1].imshow(movie[0], cmap='gray');
ax[1].grid('off'); ax[1].axis('off')
ax[0].set_title('Motion Sinks');
ax[1].set_title('Motion Sources')
ax[0].imshow(equalize_hist(motion_saliency_F_smooth), cmap='coolwarm', alpha=0.5, vmin=0, vmax=1)
ax[1].imshow(equalize_hist(motion_saliency_B_smooth), cmap='coolwarm', alpha=0.5, vmin=0, vmax=1)
plt.show()
Figure 17. Motion saliency map to highlight spatial areas of motion sinks (A) and sources (B) during Drosophila embryogenesis