# lCuts: Linear clustering of bacteria using recursive graph cuts

###### Abstract

Bacterial biofilm segmentation poses significant challenges due to lack of
apparent structure, poor imaging resolution, limited contrast between
conterminous cells and high density of cells that overlap. Although there
exist bacterial segmentation algorithms in the existing art, they fail to
delineate cells in dense biofilms, especially in 3D imaging scenarios in which
the cells are growing and subdividing in a complex manner. A graph-based data
clustering method, $\mathcal{L}$Cuts, is presented with the application on
bacterial cell segmentation. By constructing a weighted graph with node
features in locations and principal orientations, the proposed method can
automatically classify and detect differently oriented aggregations of linear
structures (represent by bacteria in the application). The method assists in
the assessment of several facets, such as bacterium tracking, cluster growth,
and mapping of migration patterns of bacterial biofilms. Quantitative and
qualitative measures for 2D data demonstrate the superiority of proposed
method over the state of the art. Preliminary 3D results exhibit reliable
classification of the cells with $97\%$ accuracy.

Index Terms—  Segmentation, bacterial biofilm, clustering, graph cut, point
cloud data

##  1 Introduction

Analyzing cellular behavior of individual bacteria in a biofilm is a key for
biologists and biochemists to understand biofilm growth, in diverse
applications such as electrical power and public health research [1]. Lack of
knowledge in macroscopic biofilm properties (e.g. size, shape, cohesion /
adhesion) that emerge from the behaviors of individual bacteria in different
micro-environment is a major barrier in biofilm studies. To make up for the
deficiency, an advanced image analysis toolkit for segmenting individual cells
is in high demand along with efficient image acquisition methods, such as
using super-resolution technology [2][3] that overcome the diffraction limit
of traditional optical microscopy techniques.

Fig. 1: Performance of $\mathcal{L}$Cuts on 3D point cloud data comparing with
manually grouped ground truth. The counting accuracy is $97\%$ and grouping
accuracy is $90\%$ . Three viewpoints from left to right: 3D view, $xy$-plane
and $yz$-plane.

The segmentation of individual bacterial cells in dense bacterial biofilms is
a challenging problem. One of the major challenges to the state of the art
comes from the presence of inhomogeneous fluorescence intensity within a
single cell or across multiple cells. When using standard level set
segmentation methods [4] and level sets using Legendre polynomials as basis
functions [5], the segmentation fails where the contrast between the cells and
background is weak. The watershed algorithm [6][7] uses the gradient flow to
identify the morphological changes along segment contours, whereas [8] [9]
separate large segments at the concavities. With both approaches, situations
where the intensity of the regions of interest is non-homogeneous often lead
to segmentation errors. Other edge-based parametric methods [10][11]are
insufficient given the subtle and often absent boundaries between cells that
are densely packed in three dimensions.

To achieve 3D cell segmentation, the authors in [12] presented a technique to
track bacteria in a dense mono-layer by way of 3D time-lapse images. This
solution employs an iterative threshold-based approach, which is heavily
dependent on high contrast between the signal and background in the images.
Yan et al. [13] proposed a single cell tracking toolkit based on marker
controlled watershed and threshold techniques. This method allows tracking of
bacterial growth in multi-layered biofilms when florescence intensity is
uniform and void spaces between cells are readily discernable, but struggles
with the detection of individual bacteria when cells are closely packed or
when inter- and intra-cellular fluorescence intensity are not heterogeneous.
Building on the work in [13], Hartmann et al. [14] recently reported a
solution to 3D segmentation in confocal images of biofilms that exploits prior
knowledge of cell size to segment low density biofilms. As this method, like
that of [13], is watershed-based, it suffers from similar drawbacks. In [15],
the authors attempted to solve the problem via constructing single-cell
regions to ensure the gap between neighboring cells in a seeded iterative
active contour approach. The single cell identification performance degrades
in the cases where the contrast between cells and voids in the biofilm is low.

As a solution to overcome the aforementioned limitations, namely the
difficulty in segmenting dense aggregations in large biofilm with non-
homogeneous inter- and intra-cell intensities, a novel approach is proposed in
this paper with two major contributions:

  * •

The bacterial cell segmentation problem is transformed into a data clustering
problem by generating pointillist data that represents the regions of
interest;

  * •

A recursive multi-class linear data clustering algorithm ($\mathcal{L}$Cuts)
is proposed that is capable of finding the linear structures in point cloud
data where cell boundaries may be ambiguous.

Our approach is built on the following insight: Even though the raw image data
does not show distinct boundaries in intensity between densely packed cells,
we are still able to reliably compute local intensity maxima that delineate
the central axis of each cell. Therefore, the proposed $\mathcal{L}$Cuts
algorithm first derives these maximal points and then partitions them based on
the approximate co-linearity of points. Moreover, this local maximum-based
initialization translates seamlessly and robustly into the 3D imaging and 3D
segmentation problem.

##  2 Linear Clustering Algorithm

Numerous algorithms exist in the clustering community that group the data by
finding the similarities between classes. Distance, number of neighbors,
density and predefined probability distribution functions are the major
perspectives for measuring similarities between points in the trending
literature, such as k-means [16], DBSCAN [17], DensityClust [18]. Among those,
density based clustering methods ([17, 18]) detect non-spherical arbitrary
clusters; however, they are still limited in precisely classifying linear
groups as discussed in the comparison in sec 3. The Hough transform [19] is
well known in detecting lines in the space, but the approach is not sufficient
for delineating cells that are intersecting and is also computationally
expensive. Unlike k-means and DensityClust, our approach does not require
manual intervention in order to locate appropriate number of clusters.
Incorporation of structural constraints, such as the distance limit and the
eccentricity of the bacteria into $\mathcal{L}$Cuts obviates the need for a
priori information regarding the number of clusters (in our case, bacteria),
making $\mathcal{L}$Cuts a fully automatic approach.

In the paper, we propose a recursive graph cuts algorithm (see work flow in
Fig. 2) for efficient computation to find the linear groups in the point cloud
data. The algorithm can be primarily divided into three parts: construct the
graph (sec 2.1), compute the bi-partition solution, and recursively re-
partition until the stopping criterion is satisfied (sec 2.2). The bi-
partition solution to separate the nodes (the local maxima) is inspired by
[20]. They addressed the problem ”how to find the precise groups in an image”
as normalized graph partitioning, where an image is partitioned into two
groups, $A$ and $B$, by disconnecting the edges between these two groups.

Fig. 2: Intuitive work flow for the recursive program. Left: an example of bi-
partition decision tree. Right: detailed example for checking the stopping
criterion of component red. Here, $sizeLimit$, $distLimit$ and $eccLimit$ are
parameters based on prior biological information.

###  2.1 Graph construction

Nodes: The nodes (local maxima along the ridgeline of a cell) in the
constructed graph have two features: location ($nodeLoc$) and direction
($nodeDir$). Location is simply the Cartesian position of the node. Direction
of each node is the principal axis direction of the ridgeline computed via
majority voting (see Fig. 3). A ”neighborhood” consists of multi-hop neighbors
that is constructed for voting. In graph theory, a ”hop” between two nodes is
defined as the number of edges that one has to traverse in order to reach from
one node to the other node.

Fig. 3: An illustration of majority voting. (a) A 4-hop ”neighborhood”
example. Each hop-neighbor is found within a specified distance (dashed
circles) to node. (b) The dashed lines connecting target node with all the
other nodes in the neighborhood are possible orientations. (c) Those
orientations that have larger relative angles with respect to the orientations
are excluded from the candidates. (d) The direction to represent the target
node is determined as the average orientation from the candidates.

A $N_{r}\times N_{p}$ accumulator is set up for the majority voting. One
dimension of this accumulator represents the $N_{p}$ possible orientations
($\mathbf{p}$) in $N_{p}$ bins (see Fig. 3b). Another dimension corresponds to
the quantized relative angles ($\phi$) with $N_{r}$ bins, where $\phi$ is
computed from each possible orientation to all the others. Here, $N_{r}$ is
chosen based on the ”hop” number. The accumulator will count the number of
parameter pairs ($\mathbf{p},\phi$) that lie in each bin. Within the first bin
of $\phi$, the orientations with the largest value are selected which give the
candidate directions. These candidates are averaged to yield the major
direction for the target node.

Adjacency matrix: The adjacency matrix reflects the likelihood if two nodes
are in the same group. Suppose there are $N$ nodes in the graph, then the
dimension of the adjacency matrix is $N\times N$. Each attribute in the matrix
represents the connectivity and edge weight between two nodes ($i$, $j$),
which measures the similarity of their features according to:

$w_{ij}=w_{distance}\cdot w_{direction}\cdot w_{intensity}\vspace{-0.2cm}$ (1)

Three similarity measures are involved: the Euclidean distance of locations of
nodes (eq 2), the relative angle between major directions (eq 3) and the
dissimilarity of intensity along the segment connecting two nodes (eq 4).

The first term is straightforward with an additional condition that sets the
weights to be zero when two nodes are farther than a given distance $r$ (set
by maximum cell length).

$w_{distance}=\begin{cases}e^{-D_{ij}^{2}/\sigma_{D}^{2}},&{\text{if
$D_{ij}\leq r$}}\\\ 0,&{\text{otherwise}}\end{cases}\vspace{-0.2cm}$ (2)

where $D_{ij}=||\text{nodeLoc}_{i}-\text{nodeLoc}_{j}||_{2}$ and $\sigma_{D}$
reflects the allowed variance for distance between nodes. The second part
measures the angle difference between two node directions, called relative
angle. Given two node directions, the relative angle ($\theta$) is the cosine
term that varies from 1 to 0 as $\theta$ becomes larger. Then the
corresponding weighting is given by:

$w_{direction}=e^{-(\cos(\theta)-1)^{2}/\sigma_{T}^{2}}\vspace{-0.2cm}$ (3)

By adjusting $\sigma_{T}$ , one can control the variance of relative angles
within each group.

The third term in (1) detects the intensity dissimilarity along the segment
joining two nodes in the image, which is defined as:

$w_{intensity}=\begin{cases}\min I_{i\rightarrow j},&{\text{if $\min
I_{i\rightarrow j}\leq\textit{thresh}$}}\\\
1,&{\text{otherwise}}\end{cases}\vspace{-0.2cm}$ (4)

Here, thresh equals the difference between the midrange (Mid) of all the nodes
and the variance (Var) of the constituent node intensities. In the case that
the nodes have no intensity information, this term can be set as $\mathbf{1}$.
Otherwise, we extract the intensities along the connecting segment from node
$i$ to node $j$ from the image as shown in Fig. 4 and compute the lowest
intensity along the segment and compare to thresh.

Fig. 4: Illustration and motivation of defining intensity on edge weight. a:
Nodes are denoted as red asterisk. b: After extracting the node directions
from the red region in a, it is still hard to separate two groups as the
relative angle (in c) and relative distance (shown in d, the distance is 10)
are close. In this case, we evaluate the intensity along the connection of the
two nodes. The intensity changes are shown in d. The intensity weighting is
then assigned as the lowest intensity lower than $thresh$.

###  2.2 Stopping criterion for recursion

Two stopping conditions are checked after each bi-partition level to decide
the completeness of the recursion.

Criterion 1 - size: The preliminary components that are less than $sizeLimit$
have the potential to be an individual group. This $sizeLimit$ is a user
defined parameter. For the application we discuss in the paper, we used the
prior bio-information of the maximum length of the bacterium to determine the
value.

Criterion 2 - linearity: This criterion is designed for preserving the linear
groups with different size from the potentials (after criterion 1).
Intuitively, if a single component is found (see black nodes in Fig. 2) and it
is less than the maximum size limit as specified, it may not be a finalized
group as linearity remains to be checked. Three aspects are checked to ensure
the linearity: (1) Standard deviation (Std) from nodes in the group to the
least square fitted line; (2) Intensity changes between the nodes within the
group (as explained in Sec 2.1); (3) Eccentricity of the group. This is an
optional condition based on the data type. For linear components, the
eccentricity ($eccLimit$) is closer to 1; while, for circular components, it
is closer to 0.

##  3 Application and Analysis

###  3.1 Experiments on bacterial images

For qualitative and quantitative assessments, $\mathcal{L}$Cuts is tested on
10 two-dimensional point cloud data which are generated from bacterial images
using Airyscan microscopy. From these images, we obtain prior information
regarding the longest bacterium in the dataset (approximately maximum 60
pixels in length and 15 pixels in width, where each pixel is $46nm\times
46nm$). The typical data have 250 to 600 nodes with approximately 20 to 60
cells observed.

Fig. 5: Pipeline for finding nodes from bacterial images. Step 1: Filter the
original image with a Gaussian kernel (a $\rightarrow$ b). Step 2: Enhance the
signals in the image via background subtraction (b $\rightarrow$ c). Step 3:
Find the local maxima (c $\rightarrow$ d). Step 4: Clear the points if they
have no neighbors or overlap with other points and rest are the found nodes
(red asterisks in d).

To build the graph, we generate the point cloud data following the pipeline in
Fig. 5. An experimental result and corresponding node features is shown in
Fig. 6.

Fig. 6: An example performance of $\mathcal{L}$Cuts with constructed graph
features. (a) Nodes are marked in red asterisks. (b) Red lines show the major
direction features for each nodes (blue dots). (c) $\mathcal{L}$Cut clustering
results for the constructed graph.

###  3.2 Qualitative and quantitative comparison

The performance of $\mathcal{L}$Cuts is analyzed qualitatively and
quantitatively by comparing with two current methods used in the bioimaging
community, DensityClust [18] and Single Cell tracking [13]. Based on the
imaging technique and biological cell information, the parameter settings for
$\mathcal{L}$Cuts are $sizeLimit=60\text{ pixels}$, $distLimit=5\text{
pixels}$ (maximum distance between nodes that have neighborhood), and
$eccLimit=0.9$. The parameters are also tuned in the other two algorithms to
achieve optimal performance in each case. In DensityClust, we chose ”Gaussian”
mode for computing densities. Due to the manual input for the selection of
cluster centers, we performed five times for each data and chose the best
performance from all. In Single Cell tracking, the watershed value is the key
to optimizing the algorithm, where a value of one is used. Qualitative
comparison is shown in Fig. 7.

Two measures, grouping accuracy (GAcc) and counting accuracy (CAcc), are
computed for quantitative comparison using $Dice$ = 2TP/(2TP \+ FP +FN), where
$TP$= true positive, $FP$ = false positive, and $FN$ = false negative. GAcc
accounts for the performance of how many nodes are correctly classified in
each group (cell); while CAcc indicates the classification accuracy in terms
of matching the final clusters with individual cells in the image. Here,
individual cell regions are manually labeled as ground truth in the
comparison.

Fig. 7: Qualitative comparison for proposed method (first column) with DensityClust [18] (second column) and Single Cell Tracking [13] (third column). For $\mathcal{L}$Cuts and DensityClust, different groups are marked with different colors and shown on the original image. The results of Single Cell Tracking are shown by overlapping the point cloud data on the segmented image, where different colors represent different single cell groups. |  $\mathcal{L}$Cuts | DensityClust | SCT
---|---|---|---
| GAcc | CAcc | GAcc | CAcc | GAcc | CAcc
Best | 95.9 | 95.1 | 94.6 | 92.3 | 94.4 | 94.1
Worst | 87.8 | 85.2 | 78.3 | 83.7 | 77.8 | 53.5
Avg |  91.6 |  91.2 | 85.9 | 87.2 | 87.7 | 86.5
Table 1: Quantitative comparison of $\mathcal{L}$Cuts with DensityClust [18]
and Single Cell Tracking [13] using Dice scores.

Overall, $\mathcal{L}$Cuts outperforms DensityClust and SCT in GAcc and CAcc
by a margin of at least $4$% on average. There are circumstances that some
cells are misclassified in $\mathcal{L}$Cuts. One cause is the non-linearity
of auto-produced point cloud data, especially when cells are randomly floating
in the three-dimensional space. Another cause is the trade-off between the
tolerance in distance/intensity changes and the continuity of the linear
structure.

$\mathcal{L}$Cuts can be directly applied on three-dimensional data. A
preliminary result is shown in Fig. 1 with a Counting Accuracy of $97\%$. The
point cloud data was generated by biofilm researchers in Gahlmann Lab. They
manually labeled the centers of each bacteria slice by slice from $x$, $y$ and
$z$ directions in Lattice Lightsheet microscopic image. The ground truth was
manually grouped which reflects the actual single bacterium layout in 3D
space.

##  4 Conclusion

We presented $\mathcal{L}$Cuts, a graph-based solution for finding linear
structures in multi-dimensional spaces. $\mathcal{L}$Cuts outperforms the
existing methods in majority of cases. Furthermore, $\mathcal{L}$Cuts enables
automated processing of 2D and 3D images to identify individual bacteria in
biofilms independent of the number of bacteria present. $\mathcal{L}$Cuts
provides quantifiable information in the form of cellular positions,
orientations, and the physical contact points between them. Beyond bacterial
biofilms, $\mathcal{L}$Cuts can be extended to other biological applications
in which boundaries are elusive but ridgelines of objects are accessible.

## References

  * [1] Carey D Nadell, Knut Drescher, and Kevin R Foster,  “Spatial structure, cooperation and competition in biofilms,”  Nature Reviews Microbiology, vol. 14, no. 9, pp. 589, 2016.
  * [2] Andreas Gahlmann and WE Moerner,  “Exploring bacterial cell biology with single-molecule tracking and super-resolution imaging,”  Nature Reviews Microbiology, vol. 12, no. 1, pp. 9, 2014.
  * [3] Steffen J Sahl, Stefan W Hell, and Stefan Jakobs,  “Fluorescence nanoscopy in cell biology,”  Nature reviews Molecular cell biology, vol. 18, no. 11, pp. 685, 2017.
  * [4] Scott T Acton and Nilanjan Ray,  “Biomedical image analysis: Segmentation,”  Synthesis Lectures on Image, Video, and Multimedia Processing, vol. 4, no. 1, pp. 1–108, 2009.
  * [5] Suvadip Mukherjee and Scott T Acton,  “Region based segmentation in presence of intensity inhomogeneity using legendre polynomials,”  IEEE Signal Processing Letters, vol. 22, no. 3, pp. 298–302, 2015\.
  * [6] Luc Vincent and Pierre Soille,  “Watersheds in digital spaces: an efficient algorithm based on immersion simulations,”  IEEE Transactions on Pattern Analysis & Machine Intelligence, , no. 6, pp. 583–598, 1991.
  * [7] Anthony S Wright and Scott T Acton,  “Watershed pyramids for edge detection,”  in IEEE International Conference on Image Processing (ICIP). IEEE, 1997, vol. 2, pp. 578–581.
  * [8] Matthew A Reyer, Eric L McLean, Shriram Chennakesavalu, and Jingyi Fei,  “An automated image analysis method for segmenting fluorescent bacteria in three dimensions,”  Biochemistry, vol. 57, no. 2, pp. 209–215, 2017.
  * [9] Yong He, Hui Gong, Benyi Xiong, Xiaofeng Xu, Anan Li, Tao Jiang, Qingtao Sun, Simin Wang, Qingming Luo, and Shangbin Chen,  “icut: an integrative cut algorithm enables accurate segmentation of touching cells,”  Scientific reports, vol. 5, pp. 12089, 2015.
  * [10] Nilanjan Ray and Scott T Acton,  “Active contours for cell tracking,”  in Image Analysis and Interpretation, 2002. Proceedings. Fifth IEEE Southwest Symposium on. IEEE, 2002, pp. 274–278.
  * [11] A-R Mansouri, Dipti Prasad Mukherjee, and Scott T Acton,  “Constraining active contour evolution via lie groups of transformation,”  IEEE Transactions on Image Processing, vol. 13, no. 6, pp. 853–863, 2004.
  * [12] Sajith Kecheril Sadanandan, Özden Baltekin, Klas EG Magnusson, Alexis Boucharin, Petter Ranefall, Joakim Jaldén, Johan Elf, and Carolina Wählby,  “Segmentation and track-analysis in time-lapse imaging of bacteria,”  IEEE Journal of Selected Topics in Signal Processing, vol. 10, no. 1, pp. 174–184, 2016.
  * [13] Jing Yan, Andrew G Sharo, Howard A Stone, Ned S Wingreen, and Bonnie L Bassler,  “Vibrio cholerae biofilm growth program and architecture revealed by single-cell live imaging,”  Proceedings of the National Academy of Sciences, vol. 113, no. 36, pp. E5337–E5343, 2016.
  * [14] Raimo Hartmann, Praveen K Singh, Philip Pearce, Rachel Mok, Boya Song, Francisco Díaz-Pascual, Jörn Dunkel, and Knut Drescher,  “Emergence of three-dimensional order and structure in growing biofilms,”  Nature Physics, p. 1, 2018.
  * [15] J Wang, R Sarkar, A Aziz, Andrea Vaccari, A Gahlmann, and Scott T Acton,  “Bact-3d: A level set segmentation approach for dense multi-layered 3d bacterial biofilms,”  in 2017 IEEE International Conference on Image Processing (ICIP). IEEE, 2017, pp. 330–334.
  * [16] James MacQueen et al.,  “Some methods for classification and analysis of multivariate observations,”  in Proceedings of the fifth Berkeley symposium on mathematical statistics and probability. Oakland, CA, USA, 1967, vol. 1, pp. 281–297.
  * [17] Martin Ester, Hans-Peter Kriegel, Jörg Sander, Xiaowei Xu, et al.,  “A density-based algorithm for discovering clusters in large spatial databases with noise.,”  in Kdd, 1996, vol. 96, pp. 226–231.
  * [18] Alex Rodriguez and Alessandro Laio,  “Clustering by fast search and find of density peaks,”  Science, vol. 344, no. 6191, pp. 1492–1496, 2014.
  * [19] Priyanka Mukhopadhyay and Bidyut B Chaudhuri,  “A survey of hough transform,”  Pattern Recognition, vol. 48, no. 3, pp. 993–1010, 2015.
  * [20] Jianbo Shi and Jitendra Malik,  “Normalized cuts and image segmentation,”  IEEE Transactions on pattern analysis and machine intelligence, vol. 22, no. 8, pp. 888–905, 2000.