Content-based Propagation of User Markings for
Interactive Segmentation of Patterned Images

Benefits of user input with real-time feedback have been recognised in image segmentation. A comprehensive summary of interactive approaches can be found in boykov2014computer. Here, we review some important advances to place our method in the existing framework, and to explain how our method differs from the current trends in interactive segmentation.

Early interactive techniques for segmentation of highly complex images include intelligent scissors mortensen1995intelligent or live wire falcao1998user, where the user cuts out an object by placing markers along its boundary. These algorithms are computationally cheap but require a lot of user effort to obtain a segmentation. Less user input is required when using interactive graph cuts boykov2001interactive; boykov2006graph, which often give very impressive results with only a few seeds provided by the user. In the GrabCut method rother2004grabcut the user provides a bounding rectangle, often leading to very precise foreground-background separation. Extensions of GrabCut include shape priors price2010geodesic and and improvement to graph cut energy representation tang2013grabcut. An alternative to combinatorial graph-based solutions is the use of a curve to represent the segmentation boundaries. Such interactive active contours often minimise an energy functional in a variational framework unger2008tvseg; santner2011interactive.

Common to the described methods is the focus on segmenting relatively large foreground objects, which justifies using regularisation on the length or the curvature of the segmentation boundary. In some applications, it is, however, not possible to use a strong regulariser. For example, when segmenting the bee eyes shown in Fig. 1, regularisation could remove or merge small regions.

The need for segmenting a number of small objects is often seen in areas like microscopy for life science or materials science. The appearance of such images can vary significantly, with texture as well as intensity carrying information that is useful for obtaining the desired segmentation. A specialist would use such clues to distinguish amongst structures, but automating the segmentation task typically requires highly sophisticated and problem-adapted methods. While there are situations that justify the development of a specialised method, in many cases a reasonable result with modest interactive effort would be preferred.

When segmenting small image structures, e.g. cells, a well-suited approach is classification of pixels. This is the basis for the ilastik segmentation tool berg2019ilastik; sommer2011ilastik, which employs a random forest classifier breiman2001random trained on image features including colour, edges, orientation and texture. The features are computed from the image before starting the interactive labelling of image structures, while parameters of the random forest classifier are learned from the manual labelling. When a user updates the labels to improve the segmentation, the parameters of the classifier need to be re-learned, which is computationally costly and causes a noticeable delay in the feedback. Another specialised tool for segmentation of microscopy images is the trainable Weka segmentation arganda2017trainable (a part of the Fiji schindelin2012fiji distribution of ImageJ) which utilises a data mining and machine learning toolkit for solving pixel classification problems. A user can choose from a variety of image features and interactively re-train the classifier.

Frameworks using neural networks are increasingly popular in pixel classification, and often yield impressive results lecun2015deep. A neural network operates on features extracted locally from the image. This input is fed through a series of multidimensional linear functions, with a non-linear activation between them, ending up in a probabilistic output. The weights of the linear functions need to be trained by optimising the performance on the usually large set of prelabelled data. This provides extreme flexibility to the method and, provided adequate training, neural networks may solve pixel classification problems as accurately as specialists. However, neural networks are dependent on large training sets and require computationally costly training, which makes them less convenient for the task of segmenting a small set of images.

Our approach shares some similarities with neural networks. We also feed the input through linear functions with non-linear steps in between. However, we use the features extracted from the image to construct the linear functions in a preprocessing step. The functions are then kept fixed, while they operate on the interactively provided user input, resulting in a probabilistic output. Due to the fixed linear functions, our method is not as adaptable as neural networks. For example, our approach is less fit for semantic segmentation of photographs. Nevertheless, we achieve excellent results when segmenting patterned images, without requiring a large set of labelled data and without performing a costly optimisation during interactive update.

The foundation of our method is a linear operator encoding image content using image–dictionary relationship. Similar approach, without the interactive update, has been used for evolving deformable models dahl2014dictionary; dahl2015dictionary; dahl2017probabilistic, quantifying composite materials emerson2018statistical; emerson2018quantifying and measuring retinal microvasculature engberg2020automated. In this work, we use the image–dictionary relationship to propagate the brush strokes provided by the users.

The aim of preprocessing is to find the structures in the image without considering the user-provided labels. In the framework of semi-supervised learning, a cluster assumption states that, if points are in the same cluster, they are likely to be of the same class – which does not imply that each class forms a single cluster chapelle2009semi. For our purpose, we assume that image features tend to form discrete clusters and that image features in the same cluster are more likely to share a class. However, we do not assume that each class is represented by only one cluster, so we will need many more clusters than classes. Therefore, we create a multitude of clusters to capture the variety of features present in the image.

In Sec. 3 we will explain the implementation details and some more advanced ways of accomplishing clustering. In this section, we outline the basic approach, which operates on intensity patches. For this case, only two parameters are required: the number of clusters $K$ and the size of the patches $M$. The number of clusters should be large, measured in hundreds or thousands, and is roughly reflecting the variability in the image. The size of the patches should reflect the scale of the distinctive image features and could, for example, be 9 pixels. For simplicity, we always assume that the size of the image patches $M$ is odd and patches are centred around the central pixel.

We extract patches of size $M$-by-$M$ from the image $I$, treat each patch as a vector containing the pixel intensities and group those vectors into $K$ clusters, e.g. using $k$-means clustering based on Euclidean distance. The resulting collection of cluster centres represents the content of the image. As these basic elements are inferred by grouping features from the image, we call the collection of $K$ cluster centres an intensity dictionary, and each of its elements (each cluster centre) is denoted dictionary patch. Every image pixel $(x,y)$ in the centre of an $M$-by-$M$ image patch is, by means of clustering, uniquely assigned to one cluster. We represent this using an assignment image $A$. For boundary pixels we define $A(x,y)=0$.

According to the cluster assumption, image patches assigned to the same dictionary patch are more likely to belong to the same class. Unique for our method is that we use this assumption on a pixel level, and not on a patch level. That is, if two image patches are assigned to the same dictionary patch, their corresponding pixels (i.e. the pixels at the same position in the patch) are more likely to belong to the same class. In other words, for every dictionary patch, there is a certain (unknown) classification of its individual pixels, which all assigned patches are likely to share.

To exploit this assumption, we define a binary relation between corresponding pixels assigned to the same dictionary pixel. For example, a central pixel of an image patch assigned to a certain dictionary patch relates to central pixels of all other patches assigned to the same dictionary patch. Likewise, the pixel directly above the central pixel relates to corresponding pixels in other patches, and a similar relation extends to all positions in a patch. This results in $M^{2}K$ cliques of pixels, one for every pixel in the intensity dictionary. Due to the overlap between image patches, every non-boundary pixel belongs to $M^{2}$ different cliques.

The central part of our method is concerned with transforming a user-provided partial labelling to pixelwise probabilities. The transformation matrix we use has a very simple decomposition, which makes our method efficient and allows for immediate feedback to the user. The construction of the transformation matrix is therefore fundamental for our method. However, describing how this matrix is constructed provides little intuition about our method, so we start by motivating our approach.

As covered previously, the assignment image $A$, obtained in an unsupervised manner, contains information on clusters of structures in the image $I$. At the same time, image $I$ is accompanied by the user-provided partial labelling $L$. To combine the two sources of information, we create a dictionary of labels to accompany our intensity dictionary. For each dictionary patch $k\in\{1,\dots,K\}$ we use $A$ to identify the locations of all image patches assigned to it. At those locations in the image grid we extract corresponding patches but from the labelling image $L$. For the set of related labelling patches we compute a pixelwise average for every layer. As a result, every $M$-by-$M$ dictionary patch now has a corresponding $M$-by-$M$ labelling representation consisting of $C$ layers.

When the image is fully labelled, the label image $L$ sums to one in every pixel, as only one out of $C$ classes has a label of 1. Consequently, the labelling representation of every dictionary patch also sums to one in every pixel. However, due to the pixelwise averaging, the values of this representation are not binary, they instead encode the normalised frequency of a dictionary pixel being labelled as belonging to class $c$ in the current labelling image. For this reason, we think of this labelling representation as of pixelwise probabilities of belonging to class $c$, and we call them dictionary probabilities.

Dictionary probabilities can now be pasted back into an $X$-by-$Y$ image grid, again using the location information from $A$, and again averaged in every pixel. This results in an $X$-by-$Y$ probability image $P$ consisting of $C$ layers, where $P$ is a diffused version of $L$. In other words, we use the self-similarity information encoded by $A$ to propagate the user-provided markings from $L$ onto the rest of the image.

In light of this motivation, now we turn to explain the construction of the transformation matrices used for efficient computation of dictionary probabilities and image probabilities. Fundamental for this transformation is the relation between the $X$-by-$Y$ image grid and the $M$-by-$M$-by-$K$ dictionary grid. This relation will be encoded using an $n$-by-$m$ biadjacency matrix $\mathbf{B}$, where $n=XY$ and $m=M^{2}K$. For this purpose, we need a linear (single) index for the pixels in the image and the pixels in the dictionary grid.

The linear index of an image pixel $(x,y)$ is

As for the dictionary grid, we use $(0,0,k)$ for the central pixel of the $k$-th dictionary element, and coordinates of other pixels in the patch are defined in terms of within-patch displacements $\Delta x$ and $\Delta y$, both from $\{-s,\dots,0,\dots,s\}$ with $s=(M-1)/2$. A dictionary pixel at coordinates $(\Delta x,\Delta y,k)$ has a linear index

Each assignment of an image patch centered around $(x,y)$ to a $k$-th dictionary patch centered around $(0,0,k)$ induces a relation between the $M^{2}$ image pixels and the $M^{2}$ dictionary pixels, see Fig. 3. Using $\sim$ for denoting a relation between image pixels and dictionary pixels gives

Since image patches are overlapping, every non-boundary image pixel relates to $M^{2}$ dictionary pixels. Image pixels in a boundary relate to less than $M^{2}$ dictionary pixels, and the four corner pixels relate to only one dictionary pixel. In total there are $(X-2s)(Y-2s)M^{2}$ relations between the image pixels and the dictionary pixels.

We represent the relations between $n$ image pixels and $m$ dictionary pixels using an $n$-by-$m$ biadjacency matrix $\mathbf{B}$, with elements

where $i$ and $j$ are linear indices of an image pixel and a dictionary pixel. The algorithm for constructing $\mathbf{B}$ is summarised in Alg. 1.

The biadjacency matrix $\mathbf{B}$ defines the linear mapping used to propagate the information from the image to the dictionary and vice versa. Consider a quantity defined on the image grid (e.g. user-provided markings indicating pixels which belong to class 1) arranged into a length $n$ vector $\mathbf{v}$ such that the $i$-th element contains the value of the $i$-th image pixel. Propagating these values to the dictionary is carried out by calculating a length $m$ vector

where $\mathbf{1}$ denotes a column vector of ones, while $\mathrm{diag}(\cdot)$ denotes a diagonal matrix with the diagonal defined by the argument. The $j$-th element of $\mathbf{d}$ contains the value of the $j$-th dictionary pixel computed by averaging the values of the related image pixels. The summation is accomplished by multiplying with $\mathbf{B}^{\mathsf{T}}$ while the diagonal matrix accomplishes the division with the total number of related pixels.

For this reason we define the $m$-by-$n$ transformation matrix for mapping from the image to the dictionary as

Similarly, mapping from the dictionary to the image is given by the $n$-by-$m$ matrix

Those two transformation matrices are fundamental for our method. The propagation of user-provided markings (as described in the motivational paragraphs) is computed as

where $\mathbf{L}$ is the user-provided labelling $L$ arranged in a $n$-by-$C$ matrix, while the resulting $n$-by-$C$ matrix $\mathbf{P}$ needs to be arranged back into a layered image $P$.

When equipping our method with the user-provided interactive update, we run into choices with regards to: i) how we treat unlabelled pixels, ii) the number of applied diffusion steps, and the way of treating intermediate results between the steps, and iii) the possibility of changing the number of segmentation classes. After testing many types of interactive updates, we kept three main versions. In all versions the number of classes $C$ is chosen during initialisation and kept fixed during the update.

How we handle pixels that have not been labelled by the user is also common to all versions. Such pixels are initially assigned equal probability of belonging to each class. As a result, before the user places the first label, all probabilities are equal and no segmentation is possible.

The user starts the interaction by choosing a pencil corresponding to one of the $C$ classes applies markings to some pixels. The partial labelling information is immediately transformed to the probability image and shown to the user as an image segmentation, with every pixel placed in the class with the highest probability. After the first pencil stroke, only one class will have values larger than $\frac{1}{C}$ in the label image $L$, and the same applies for the probability image $P$ computed using (8). Thus, at first, many pixels will belong to the first marked class and no pixels will be assigned to the classes that have not been marked yet. As the user adds markings for the other classes, those will appear in probability image $P$.

Thanks to the real-time feedback, the user can quickly improve the result by placing markings in misclassified regions (the regions that have been incorrectly classified). With many unlabelled pixels in $L$, the image $P$ will typically have many values that only differ slightly from $\frac{1}{C}$. Those small deviations carry the information needed for inferring the class of the unlabelled pixels.

As for the number of applied diffusion steps, we use either one or two. When using two diffusion steps, instead of continuing to diffuse the (already diffused) probability image, we can apply additional non-linear operations between the two diffusions. Very good results are obtained if we apply binarisation of the labels between the two diffusion steps. For binarisation, we identify the class of the highest probability for each pixel, and apply $\{0,1\}$ labelling. If there are pixels with no clear probability maximum, we let them retain their unresolved labels. Consequently, for the second iteration of the diffusion, many pixels act as labelled, and this improves the quality of the result.

The options for the two-step diffusion and binarisation are implemented in our segmentation tool, such that the user can quickly switch between the variants of the method and decide which one yields the best results for the data at hand. Likewise, the user can quickly determine whether the quality of the results is sufficient or additional markings should be placed.

The user can inspect the output of the classification displayed as a final segmentation based on the resulting probability image. Alternatively, there is an option for inspecting the $C$ probability images, which often gives a better insight into the quality of the result.

Our approach allows for various postprocessing options, which may be grouped into two postprocessing strategies. One strategy involves processing the probability image to obtain the segmentation or detection of interesting features from the probability image. These operations are application-driven and examples are illustrated in Sec. 4.

The second strategy involves reusing the information stored in the dictionary and the associated dictionary probabilities. The linear transformation (8), which is core to our method, first transforms the user-provided markings from $L$ to the dictionary space (using matrix $\mathbf{T}_{1}$) and then back to the image space (using matrix $\mathbf{T}_{2}$). Consider only the first product

This is an $m$-by-$C$ matrix containing the pixelwise probabilities of the dictionary pixels (i.e. the dictionary probabilities) which can be useful for processing a previously unseen image similar to $I$.

Processing a new image $\hat{I}$ requires extracting all $M$-by-$M$ patches for every pixel of $\hat{I}$ and assigning those patches to the existing dictionary, i.e. the dictionary created using patches from $I$. Just like before, this assignment defines an image-to-dictionary and we can compute the two associated transformation matrices. Here we are interested in the dictionary-to-image transformation $\hat{\mathbf{T}}_{2}$. To compute the probability image corresponding to the unlabelled image $\hat{I}$ we therefore need to compute

and rearrange the result into $\hat{P}$.

This way of using our method fits into the framework of supervised learning. The original image $I$ and the computed labelling $L$ can in this context be seen as a (labelled) training set (ignoring the fact that the labelling is computed in a semi-supervised way). Our method is then capable of producing the probability image $\hat{P}$ for the new, unlabelled image $\hat{I}$. The approach will work as long as the initial clustering captures the features present in $\hat{I}$, which holds for similar images.