A Model for Image Segmentation in Retina

Image segmentation is a challenging and important problem in computer vision and the Berkeley Segmentation Data Set (BSDS) is a standard benchmarking data set for many computer vision image segmentation algorithms [19, 1]. It consists of 500 large ($\sim 400\times 300$ pixels) color images each with multiple ($\sim 5$) human drawn boundary contours (green box in Fig. 2), as well as code provided for standard benchmarking and comparison of algorithms. Since image segmentation is closely related to boundary detection and quantification of boundary detection performance is more straightforward than that of image segmentation, segments in images are often recast as boundaries for benchmarking. Binary boundary pixel locations are compared to human drawn boundaries using the precision, recall, f-measure framework. In this context, "Precision" is the proportion of image pixels hypothesized by a method to belong to segment boundaries that agree with the ground truth. "Recall" is the percentage of ground truth boundary pixels that are found by a method. F-measure is the harmonic mean of Precision and Recall.

In order to leverage the BSDS resource, we must first convert the output of a segmentation model - a phase, spectral or feature activation map (blue box in Fig. 2) - into binary boundaries. Intuitively, a good segmentation of an image has been achieved if the model output map has very similar values within segments and large discontinuities at boundaries. We compute spatial derivatives ($\nicefrac{{\delta}}{{\delta r}}$) in the output map and normalize the values between 0 and 1, allowing us to interpret resulting probabilistic boundary (pb) as the algorithm’s confidence that there is a boundary between segments at a particular image location. We can threshold pb’s at multiple values and compare each resulting binary boundary map (bb) to each human drawn ground-truth boundary map (gT), generating a pixel match set by a logical AND operation. Because human drawn boundaries are not precise down to the pixel, we allow small misalignment between gT and bb pixel including a pair in the match set if they are within $d_{t}$ pixels of one another.

We compared the ability of different phase coupled oscillator models to segment images from the Berkeley Segmentation Dataset (BSDS) [19]. The models differed in the phase couplings. One baseline model contained isotropic couplings, while the couplings in the other models AA, GL, M and TM were the transformations of the adjacency of features described above. We set the parameters $\sigma_{f}$, $\sigma_{d}$ and $\sigma_{\omega}$ to adequate common values and performed for each method a parameter grid search in neighborhood connectivity $R_{M}$ and scale $K_{s}$ to maximize the average $\Delta F_{b}$ across 500 image patches. The oscillator frequency was 60 Hz (typical for retinal gamma oscillations) and we gave the networks 300ms to relax the phases, corresponding to an interval between saccades.

Within a stage of visual processing, in which a set of local visual features is extracted, image segmentation can be viewed as a graph clustering problem [35]. Consider an image and its corresponding neural representation in retina or LGN, in which the activity in individual cells represent the strength of local center-surround features. Image segmentation consists of clustering sets of local image features that share properties and thus likely correspond to larger objects in the image. However, much more efficient than clustering pixel values, is to apply clustering on more sophisticated local image features, e.g., center-surround, edges, multi-scale textures, as done in state-of-the-art segmentation methods [35, 32].

The problem of graph clustering, that is finding “cliques” of strongly connected nodes in a graph, is obviously related to finding pixel sets with similar local features. To recast image segmentation in terms of graph clustering, one first uses kernels to construct an adjacency matrix, in which an element is large whenever two features have similar values and lie nearby each other in the image plane [35, 32]. The segmentation of the image corresponds to finding the communities (subsets of nodes) that are strongly interconnected within the community, and well separated from nodes outside. The goal then is to find non-trivial subsets of nodes that can be separated from one another by cutting through the minimum weight of edges, know as the "mincut" problem. Though a brute force, optimal solution to this problem would be combinatorically intractable, approximate solutions can be found efficiently by leveraging the machinery of "spectral graph theory" [8].

The simplest strategy of graph clustering, referred to as average association (AA), is to analyze the adjacency matrix directly [32]. Eigenvalues of the adjacency quantify the amount of correlated structure and the associated eigenvectors characterize the location of the correlated structure in the image. Other methods of graph clustering utilize transformations of the adjacency matrix, often incorporating the node "degree", $d_{i}=\sum_{j}A_{ij}$, which captures the total weight of connections to each node from all other nodes in the network. One such transformation we considered is the normalized graph Laplacian (GL) or Kirchhoff matrix: $L=D^{-1/2}(D-A)D^{-1/2}$ with diagonal matrix $D_{ij}=\delta_{ij}\sum_{k}{A_{kj}}$, $\delta_{ij}$ the Kronecker symbol. This strategy, combined with more sophisticated image features, forms the basis of a very successful image segmentation algorithm, the “Normalized Cut” [35]. The eigenvectors and associated smallest eigenvalues of the Laplacian matrix find divisions in the input characterized by large feature differences.

A third transformation of the association matrix we considered is modularity (M) [26], which has successfully discovered community structure in social and information networks, outperforming the graph Laplacian in these tasks. The modularity matrix can be written as

$$Q=A-N\quad\text{with}\quad N_{ij}=D_{i}D_{j}\quad\text{and}\quad D_{i}:=\frac{d_{i}}{\sqrt{2m}}$$

(2)

Once an associated matrix representing a graph is constructed, spectral methods have been predominantly used within the graph clustering community to find clusters within because eigenvalues and eigenvectors efficiently find an approximate solution to the combinatorially intractable "mincut" problem. It has been observed on simple networks that the eigenvalue spectrum of an associated matrix resembles the temporal progression of clusters discernible from phases of nodes in a Kuramoto network [2], this time evolution of clusters forming a hierarchical clustering of a network. Given this observation, we compute the time evolution of a phase coupled oscillator network dynamical system as an alternative to eigenvector-based graph clustering methods.

The described graph clustering methods in 2.2 compute the eigenvectors of the associated matrices [8] which, in essence, is assessing anisotropic diffusion in these networks. This process has also been related to the path a random walker would take through the graph where edge weights represent transition probabilities and the distribution of electrical potentials on nodes in a resistor network where an edge weight represents the conductance of a particular resistor [12]. A further parallel has been between eigenvector based methods for graph clustering and the "fundamental mode(s) of a spring-mass system" [35]. To rigorously investigate this last claim, we simulate phase relaxation in a network of Kuramoto coupled oscillators [18] with networks defined by methods described in 2.2.

Here we followed [2] and assessed diffusion properties by relaxing a network of phase-coupled oscillators :

$$\Delta\phi_{i}=\omega_{i}+\sum_{j}{K_{ij}sin(\phi_{i}-\phi_{j})}\text{,}\quad\quad K_{ij}=k_{s}M_{ij}$$

(3)

In the implementation of the model, the overall positive scaling factor $k_{s}$ is critical. If coupling weights are too large, phasers will spin wildly in response to even small phase differences. Conversely, if too small, oscillators will adjust their phase too slowly and the relaxation will not converge in time. Importantly, the phase relaxation was limited to 300ms or 20 periods of the 60Hz signal, which is the average duration of fixation before a saccade brings the eye’s gaze to a new point, refreshing the input and beginning the computation once again. The value for the $k_{s}$ parameter was set for each graph individually based on mathematical considerations in equation 3. A middle value $k_{s}^{\tiny{mid}}$ was chosen so that the phase change of the node with largest degree $D_{k_{max}}$ is limited to $\nicefrac{{\pi}}{{2}}$ radians in one full period of the 60Hz signal when all its neighbors are aligned $\nicefrac{{\pi}}{{2}}$ radians away and exerting maximal pull.

$$k_{s}^{\tiny{mid}}=60Hz\cdot\frac{2\pi}{\nicefrac{{\pi}}{{2}}}\cdot D_{k_{max}}$$

(4)

The final result of the phase relaxation simulation is a phase map with a phase value, $\phi_{i}\in[0,2\pi]$, associated with every node, $i$, in the network and corresponding location $i$ in the image. Spectral methods also yield a value associate with each location, $i$, in the image with $v_{i}\in[-\infty,\infty]$. In order to compare our results to other algorithms using the BSDS resources, we convert these maps to probabilistic boundaries and recast the image segmentation problem as a boundary detection one as discussed in Section 2.1 and illustrated in Fig. 2.

In practice, two meta parameters, $r_{M}$ defining the neighborhood structure of the Adjacency graph and $k_{s}$ defining an overall scaling on the strength of phase interactions in the network, impacted image segmentation performance. They were optimized for each method and results shown are with optimized parameters, shown in Fig. 11. To optimize parameters for each method, we performed segmentation of 500 image patches with four $r_{M}$ values ranging from 1 to 10 and bracketing $k_{s}$ as discussed above and chose the parameter settings with best average performance across all images and across $d_{t}$. Fig. 3 illustrates the procedure for one particular method. It shows average performance across $k_{s}$ values for optimal $r_{M}$ on the left and performance across $r_{M}$ values for optimal $k_{s}$ on the right. Fig. 4 shows the effect of different parameter settings on one example image patch.

An image can be described by a multi-graph, in which pixels or local image features are represented by nodes and each pair of pixels has two different types of edges connecting them. One edge type represents geometric distance in the image plane and the other edge type represents feature differences. The two types of edges are given by adjacency matrices, resulting from the two types of distances and corresponding kernel functions (like a Gaussian kernel), as in Eq. 1. Shi and Malik [35] proposed a way to collapse this multi-graph of an image to an ordinary graph by forming the Hadamard product of the two adjacency matrices. An entry in the resulting single adjacency matrix $A$ represents the two distinct similarities between pixels, geometric proximity and feature similarity by a single number. Specifically, an entry in $A$ can only be large, if both, distance and feature differences are small in the corresponding pair of pixels. In order to find image segments, researchers then used "spectral" graph clustering methods on the matrix $A$ [32, 35].

For some graph clustering methods, such as modularity [26], the collapsing of the multi-graph into an ordinary graph destroys information, which is critical for segmenting images. The modularity matrix consists of the difference of the adjacency matrix and a null model. The null model represents an average adjacency value. In the standard modularity method, Eq. 2, the average is computed from the degrees of the two nodes involved, the row and column sum of the collapsed graph. However, in natural images, the average feature similarity of a pair of pixels is a function of geometric distance [31], see also Fig. 20 in supplemental section B.6. Thus, an appropriate null model for images should also depend on the geometric adjacency matrix.

To address this issue, we devised a novel graph clustering method called topographic modularity (TM) in which the null model takes topographic distance in the image plane into account. Like the standard modularity [26], see Fig. 5a and Eq. 2, an entry of the topographic modularity matrix, $Q^{T}$, is the difference between the entry $A_{ij}$ and the expected connectivity, captured by the null model, $N_{ij}$. Here, the topographic null model $N_{T}$ accounts for distance dependent factors in feature similarity with the $R_{ij}$ term in addition to node degree heterogeneity.

$$Q^{T}_{ij}=A_{ij}-N^{T}_{ij}\quad\text{where}\quad N^{T}_{ij}=D_{i}\cdot D_{j}\cdot R_{ij}$$

(5)

The $R_{ij}$ factor represents the average connectivity between all node pairs that are separated by the same geometric distance as the nodes $i$ and $j$. For a network in space along a 1D line, Fig. 5b, the distance dependent contribution to the null model can be written mathematically as

$$R_{ij}^{(1D)}=(\frac{1}{n-L})\sum_{k=1}^{n-L}{A_{k,k+L}}$$

(6)

where $L$ is the distance separating nodes $i$ and $j$ (i.e., $L=r_{i}-r_{j}$) and $n$ is the total number of nodes (or pixels or features in the image). In 1D, the average connectivity of all nodes separated by a distance L is equal to the mean along the L’th diagonal.

For networks with 2D grid-like geometry, like Adjacency graphs constructed from images, the computation of $R_{ij}^{(2D)}$ is more involved, yet the interpretation is the same. Reshaping a 2D image into a 1D vector so that similarity relationships can be represented in a 2D matrix introduces discontinuities in spatial relationships between entries in the matrix. Weights between nodes separated by a particular distance can be labeled by a mask specific to the dimensions of a particular image. Fig. 5c shows the weights between all neighboring nodes ($L=1$) in the network derived from the 11 x 11 binary image in the inset. For completeness, we show $R_{ij}^{(2D)}$ masks for other pixel separations in supplement section B Fig. 19.

Before comparing the different null models in an image segmentation we compare how well they capture the structure of an adjacency matrix of an image. The null model in Newman’s modularity, by construction, is a "consistent" estimator of node degrees [7], ensuring that $\sum_{j}{N_{ij}}=\sum_{j}{A_{ij}}$ (blue line in Fig. 6 middle). However, it is clearly the wrong null model for natural image Adjacency graphs for two reasons. First, the null model incorrectly contains positive diagonal weights in proportion to $D_{i}^{2}$, although the diagonal elements of the adjacency matrix are zero. Second, it does not capture the distance dependence of the adjacencies, thereby underestimating average adjacency between proximal nodes and overestimating it for distant node pairs. Both problems manifest in the difference between the blue and the dashed lines in Fig. 6 bottom.

While the TM-1D and TM-2D null models are not strictly consistent in node degree or distance dependence, they are nearly so (green and red lines respectively in Fig. 6). Introducing distance-dependent statistics into the TM-1D null model corrects for the spatial "inconsistency", vastly improving estimates of edge-weight over M. TM-2D offers improvement over TM-1D due to further refinement of its null model, see Eq. 6 and surrounding text. Further discussion of null model consistency and bias in the supplemental section LABEL:ModularityConsistency.

Importantly, the difference between an adjacency value and its average in the modularity can become negative. In a Kuramoto net relaxation, these negative weights mediate phase repulsion and introduce targeted phase desynchronization, see Sec. 2.4, at boundaries in an image where gross image statistics change. In contrast, if the modularity value between a node pair is positive, it contributes to phase synchronization. Fig. 7 illustrates image segmentation performance before and after phase relaxation through connections defined by M (in blue), TM-1D (in green) and TM-2D (in red). While M does not significantly change image segmentation performance over Gaussian RF independent sensors, TM-1D does so (p-value $\sim 0.004$) and TM-2D does so even more (p-value $\sim 4\cdot 10^{-7}$). With TM-2D, we see improvement for $\nicefrac{{\sim 460}}{{500}}$ image patches.

Following [15], we investigate the idea whether a phase-coupled network of simple sensors of local image features, similar to those in the retina, could at the same time represent local and contextual image features in its output. Specifically, phase interactions mediated through heterogeneous network edges which are influenced by local features similarities can segment an image, grouping regions within a segment into the same relative phase and introducing phase breaks at segment boundaries. In a biological system, the contextual image information encoded by phase can be represented by the timing of spikes and be multiplexed into spike trains, whose rates represent the local features Fig. 1.

This idea is tested on images provided in the Berkeley Segmentation Data set (Sec. 2.1). For an image patch, we construct a graph based on local features in the image (Sec. 2.2) and segment the image by either computing eigenvectors or by performing anisotropic phase diffusion in a Kuramoto net. Computing a spatial derivative on either eigenvectors or the final phase distributions and normalizing values between 0 and 1 converts the output into probabilistic boundaries, which can be quantitatively compared to assess relative performance of different image segmentation methods.

We ask whether a phase-coupled sensor array can add to an image segmentation that can be done based on the independent sensor measurements alone. Thus, the network computation must outperform two baseline methods. The first method computes normalized spatial gradients on the raw image pixels (magenta, RawPix). In the second method the image pixels are first convolved with Gaussian receptive fields, roughly similar to those measured in retina (cyan, GaussRF). As a third baseline method, we include isotropic diffusion in a network with homogeneous phase couplings between nearest-neighbor nodes (black, ISO).

Probabilistic boundaries (pb) can be interpreted as the algorithm’s confidence that a boundary exists between two segments at a particular location in the image. Fig. 8 shows pb’s resulting from the segmentation of different networks constructed from the same image patch, either by computing eigenvectors and by performing Kuramoto net relaxation. Qualitatively, we observe that eigenvectors seem to focus a spotlight on a region of the image patch while information propagated through the Kuramoto Net covers all parts of the image patch. Regardless of the network method used, boundaries found with the Kuramoto net are crisper and extend further across the image patch than do those found by computing eigenvectors.

To assess whether this trend in image segmentation performance is statistically significant, we calculate Precision, Recall and F-measure across 500 image patches, shown in Fig. 9. Plotting F-measure statistics for network and baseline models segmented by Kuramoto-net and Eigen-methods, we find that segmentation without network computation (magenta and cyan) outperforms the results from the best combination of the top 3 eigenvectors, regardless of the model. We also find that all scatter points lie above the unity line, indicating superior image segmentation performance of anisotropic phase diffusion in a Kuramoto net verses the spectral clustering methods. As a consequence of this observation, we focus in the reminder on the superior methods based on Kuramoto Nets.

We further observe that the features from which networks are constructed influence segmentation performance achieved. This comes as no surprise since state-of-the-art image segmentation algorithms rely on a combination of sophisticated spatially-extended features.

We constrain our investigation to the relatively simple and local stimulus features that retina is supposed to have access to. Specifically, we investigate the difference in segmentation caused by switching between raw pixels and Gaussian receptive fields with different radii. Again, we compare the segmentation performance of networks with phase relaxation to baseline models representing independent sensors, and a model with isotropic diffusion through a homogeneous neigbor connections. We find that Gaussian receptive field features provide better segmentation than raw image pixels both when used as independent sensors and to construct phase interaction networks. Fig. 10 shows the segmentation performance (F-measure and the change in F-measure relative to the independent sensors image pixels baseline model) as a function of pixel match distance tolerance ($d_{t}$).

For small tolerances $d_{t}$ in the F-measure (see section 2.2) the simpler isotropic phase diffusion model was a surprisingly strong competitor, even beating some of the anisotropic networks (black lines in Fig. 10c and d). Isotropic diffusion with optimized parameters provides mild smoothing of image structure, which operates indiscriminately within and across segments. To introduce the effect of smoothing in other models, we introduced Gaussian RF features. The filters corresponded to optical blur and the extended (centers of) receptive fields in retinal ganglion cells. Fig. 10b shows segmentation performance as a function of the width of the Gaussian filter, $\sigma$, and tolerance parameter, $d_{t}$. We find that Gaussian RF features with $\sigma=1$ were beneficial and near optimal across different tolerance values. Interestingly, the size of the optimal Gaussian coincides with the size of retinal ganglion cell receptive fields measured in primate retina [9]. See supporting information A for further discussion.

Fig. 10e and 10f show the improvement in segmentation performance using Gaussian features above using image pixel features. In particular the method TM displayed a significant increase in $\Delta F_{i}$ which became more prominent for larger pixel match distance tolerances $d_{t}$. Among all methods TM was able to improve segmentation performance the most, compared to that achievable with the Gaussian RF independent sensors model.

To assess the overall performance of different models on the diverse input images, each model was run with optimized parameters. Fig. 11 shows image segmentation performance improvement from Gaussian RF independent sensors. Here the models TM-1D, TM-2D and ISO were significantly different from the three other methods that stayed near baseline $\Delta F=0$. ISO stayed below baseline because the input kernels provide near optimal blur and therefore additional isotropic blurring deteriorated the segmentation performance. TM-1D performs well too, but not as well as TM-2D. This is because the null models are increasingly accurate, section 3.1. Shown results are with best matching ground truth. Results hold with average across all human drawn ground truths, though less pronounced.

Segmentation performance via anisotropic phase diffusion in a Kuramoto net depends critically on the structure of the phase couplings. Kuramoto nets using the graph Laplacian, average association or Newman’s modularity as the phase couplings do not improve segmentation performance significantly over the independent sensors Gaussian RF baseline model. Only the Kuramoto model with the topographic modularity as phase couplings increases segmentation information over baseline independent sensors, homogeneous network and competitor heterogeneous network models, as quantified by the F-measure.

The F-measure combines the performance measures Precision and Recall, each with intuitive interpretations described in section 2.1. To analyze the differences between our different models, we separately plot the precision and recall distributions in Fig. 12. Note the position of curves for each network method relative to the independent sensors Gaussian RF baseline model (cyan dashed curve). Focusing first on the F-measure, in panel a, three of the network models (AA in yellow, GL in green, M in blue) did not show significant differences. The ISO model (black) degraded segmentation performance while the TM models (red & magenta) improved relative to the Gaussian RF baseline. In panels b and c, the precision distribution of both TM models shifts significantly to higher values while the recall distribution shifts only slightly to lower values. Thus, the performance improvement of the TM model is mainly caused by increased precision, reflecting superior ability to suppress spurious boundaries, texture or "noise" in the probabilistic boundary maps.

To better understand the computation in the TM-2D model, we visualize changes to Precision and Recall together for individual image patches in Fig. 13.

The TM-2D network relaxation improved segmentation for $\sim 93\%$ of all image patches, in blue, panel a. Clear positive shifts in the precision and F-measure distributions can be observed from the independent sensors Gaussian RF model (dashed cyan) to the phase output from the TM-2D network relaxation (solid blue). No clear trend emerges for the recall distribution with improved images. No clear trend exists for images where the network relaxation decreased performance. For some precision increased, and recall decreased. For others, vice versa.

Finally, to provide some intuition what a $\Delta F$ value means for individual images, some examples are shown in Fig. 14. Compared to the results from other methods, the TM model produces probabilistic boundaries (pb’s) that are often thinner and cleaner.

Further, in Fig. 15, we show samples of image patches with varying image segmentation performance relative to the Gaussian RF independent sensors model.