Let us consider a Gaussian process (GP) model for a two-dimensional (2D) microscopy image, where the intensity at the pixel location x is modeled as y ( x ) = f ( x ) + ϵ, with f ( · ) being a latent Gaussian process (GP) and ϵ ∼ N ( 0 , σ 0 2 ) being a Gaussian noise with variance σ 0 2 [38]. Consider an n 1 × n 2 image Y, where the ( i , j )th term is y ( x i , j ) for the pixel location x i , j = ( x i , 1 , x j , 2 ). The data matrix can be written as: (2.1) Y = F + E , where F follows a matrix normal distribution F ∼ Matrix-Normal μ 1 n 1 × n 2 , σ 2 R 1 , R 2 with μ being a mean parameter, 1 n 1 × n 2 being a n 1 × n 2 matrix of ones, σ 2 being the variance parameter, R 1 and R 2 being correlation matrices for two inputs, respectively, and the vectorized noise matrix follows Vec ( E ) ∼ MN ( 0 , σ 0 2 I N × N ), with N = n 1 n 2 . Additional trend structure can be modeled in the mean function, and the covariance can be generalized to be semi-separable in the model [14].

Vectorize the observations matrix y v = Vec ( Y ). After marginalizing out the latent signal matrix F, the observation vector of the image follows a multivariate normal distribution, (2.2) ( y v ∣ σ 2 , η , R 1 , R 2 ) ∼ MN ( μ 1 N , σ 2 R ˜ ) , where R ˜ = R + η I N , R = R 2 ⊗ R 1 , ⊗ is a Kronecker product and η = σ 0 2 / σ 2 is a nugget parameter. Here the correlation matrices are parameterized kernel functions R 1 ( i , i ′ ) = K 1 ( x i , 1 , x i ′ , 1 ) and R 2 ( j , j ′ ) = K 2 ( x j , 2 , x j ′ , 2 ). Frequently used kernel functions include the power exponential kernel and Matérn kernel [38]. Denote K l ( x i , l , x i ′ , l ) = K l ( d l ) with d l = | x i , l − x i ′ , l |. For instance, the Matérn kernel with roughness parameter 5 / 2 follows [20]: (2.3) K l ( d l ) = 1 + 5 d l γ l + 5 d l 2 3 γ l 2 exp − 5 d l γ l , where γ l is a range parameter to be estimated for l = 1 , 2. A GP with the kernel function in Equation (2.3) is twice mean-squared differentiable [38], and it is used as a default choice of some software packages of GP surrogate models [42, 16]. We use the kernel function in (2.3) for demonstration purposes.

After specifying the model, we compute the likelihood function and predictive distribution for parameter estimation and prediction, respectively. The parameters of the model in (2.2) include the mean, variance, range and nugget parameters { μ , σ 2 , γ , η } with γ = ( γ 1 , γ 2 ) being the range parameters in the kernel function in Equation (2.3). As the number of pixels is large, we employed the maximum likelihood estimator (MLE) to estimate the parameters, as the result is similar to the robust estimation by marginal posterior mode [17]. Given the range and nugget parameters, the MLE of the mean and variance parameters has a closed-form expression: μ ˆ = ( 1 N T R ˜ − 1 1 N ) − 1 1 N T R ˜ − 1 y v , where σ ˆ 2 = S 2 N with S 2 = ( y v − μ ˆ 1 N ) T R ˜ − 1 ( y v − μ ˆ 1 N ). The logarithm of the profile likelihood after plugging the MLE of the mean and variance parameters follows: (2.4) log ( L ( γ , η ) ) = C − 1 2 log ( | R ˜ | ) − N 2 log ( S 2 ) , where C is a constant not related to ( γ , η ).

The range and nugget parameters are then estimated by maximizing the logarithm of the profile likelihood function: (2.5) ( γ ˆ , η ˆ ) = argmax ( γ , η ) log ( L ( γ , η ) ) . Vectorize signal f v = Vec ( F ) from model (2.1). After plugging the MLE of the parameters, the posterior distribution of the image vector follows (2.6) ( f v ∣ y v , μ ˆ , σ ˆ 2 , η ˆ , γ ˆ ) ∼ MN ( f v ∗ , σ ˆ 2 R ∗ ) , where the predictive mean and predictive covariance follow (2.7) f v ∗ = μ ˆ 1 N + R R ˜ − 1 ( y v − μ ˆ 1 N ) , (2.8) R ∗ = R − R R ˜ − 1 R . The predictive mean f v ∗ and predictive variances, the diagonal terms of σ ˆ 2 R ∗ , are required for prediction and uncertainty quantification. However, direct computation of the likelihood function involves inverting the covariance matrix, which takes O ( N 3 ) operations, where the number of pixels N can be at the order of 10 6 − 10 7 for microscopy images. A fast computational way without approximation is introduced in Section 2.2 to solve this computational challenge.

A wide range of approximation approaches of GPs have been developed [48, 28, 43, 7, 23, 13, 19, 24, 57]. For GPs with product covariances on image data, no approximation is required. Denote the eigendecomposition of the subcovariance R l = U l Λ l U l T , where U l and Λ l are matrices of eigenvectors and a diagonal matrix of eigenvalues for l = 1 , 2, respectively. Furthermore, denote λ l = ( λ 1 , l , … , λ n l , l ) T , an n l -vector of the diagonal terms in Λ l , for l = 1 , 2.

The log profile likelihood in Equation (2.4) and the predictive distributions in Equation (2.6) can be written as a function in terms of eigenpairs in Lemma 1 below. Lemma 1 eliminates the need to invert the N × N covariance matrix and reduce the computational complexity to O ( n 1 3 + n 2 3 ), which is significantly smaller than O ( n 1 3 n 2 3 ). The proof of Lemma 1 is provided in the Appendix. Lemma 1.

1.

(Profile likelihood). Denote the transformed likelihood vector y ˜ v = Vec ( Y ˜ ) = Vec ( U 1 T ( Y − μ ˆ 1 n 1 × n 2 ) U 2 ). The logarithm of the profile likelihood in Equation (2.4) can be written as (2.9) log ( L ( γ , η ) ) = C − 1 2 ∑ i = 1 n 1 ∑ j = 1 n 2 log λ i , 1 λ j , 2 + η − N 2 log ∑ i = 1 n 1 ∑ j = 1 n 2 Y ˜ i , j 2 λ i , 1 λ j , 2 + η , where (2.10) μ ˆ = ∑ j = 1 n 2 ∑ i = 1 n 1 u ˜ i , 1 2 u ˜ j , 2 2 λ i , 1 λ j , 2 + η − 1 × ∑ j = 1 n 2 ∑ i = 1 n 1 u ˜ i , 1 Y ˜ i , j , 0 u ˜ j , 2 λ i , 1 λ j , 2 + η , with Y ˜ i , j being the ( i , j )th entry in ( U 1 T ( Y − μ ˆ 1 n 1 × n 2 ) U 2 ), Y ˜ i , j , 0 being the ( i , j )th entry of the n 1 × n 2 matrix U 1 T Y U 2 , u ˜ i , 1 being the ith term of the n 1 vector 1 n 1 T U 1 i = 1 , … , n 1 and u ˜ j , 2 being the jth term of the n 2 vector 1 n 2 T U 2 , for j = 1 , … , n 2 .

Let Y represent the original cell image, which has size N = n 1 × n 2 . We denote y ( x i , j ) as the value of the image at pixel x i , j = ( x i , 1 , x j , 2 ). All possible pixel values are between 0 and 1, as each image in this analysis is in grayscale. Microscopy images of cells can contain more than 10 6 pixels, and some areas of an image can have substantially higher pixel values than other areas. Thus, for a large image Y, it is normally partitioned into a set of cropped images, denoted as { Y k } k = 1 K ˜ , where K ˜ is the total number of cropped images as illustrated in Figure 1. This step can substantially reduce the computational cost, and enable segmentation outcomes more adaptive to the features of the local areas in a large image.

As experimental images often contain substantial noise, for each cropped image, we first use GPs with fast computation introduced in Section 2.1 for denoising the images. The range and nugget parameters of GPs are estimated by maximizing the likelihood in Equation (2.9), and they are held the same for all sub-images as they have similar smoothness properties. The mean and variance parameters are estimated differently for each cropped image to capture the local change of the optical properties. As the sub-images will be turned into a binary image, the discontinuity of the boundary between different sub-images does not impact the image segmentation task.

The predictive mean at the interior pixel of a cell object typically has a larger value than the one in the background of a microscopy image. Inspired by the IsoData algorithm [40] in ImageJ [1], we developed an automated and robust way to determine the threshold value for separating objects from background noise using the predictive mean from GP models.

Suppose each cropped image is n k , 1 × n k , 2 with predictive mean F k ∗ with the ( i , j )th entry being F k ∗ ( x i , j ) and Vec ( F k ∗ ) = f k , v ∗ given in Equation (2.11). As the intensity values of the pixels are between 0 to 1 and the predictive uncertainty is small, the range of F k ∗ ( x i , j ) is also approximately between 0 to 1. As the image has only one type of cells, the change of pixel values F k ∗ ( x i , j ) is large if x i , j is the pixel at the boundary of a cell object, and small elsewhere. Thus we first compute the normalized predictive mean exceeds each threshold value: c k ( α m ) = ∑ i , j 1 F k ∗ ( x i , j ) max ( F k ∗ ) > α m , where { α m } m = 1 M a sequence of equally spaced thresholds from 0 to 1 with the default value M = 100, and max ( F k ∗ ) denotes the maximum value of the predictive mean F k ∗ of the kth cropped image for k = 1 , … , K ˜. In general, the range of { α m } m = 1 M can be chosen to be the range of F k ∗ ( x i , j ) / max ( F k ∗ ). Subsequently, we calculate the difference Δ c k ( α m ) in pixel counts below Δ c k ( α m ) = c k ( α m − 1 ) − c k ( α m ) . As Δ c k ( α m ) may not be smooth over α m , we compute the predictive mean of Δ c k ( α m ), denoted by Δ c k ∗ ( α m ), via a GP model with the default Matérn kernel in Equation (2.3) by the RobustGaSP package [16].

We define α ∗ as the smallest threshold α m after the point of maximum fluctuation, denoted as arg max m Δ c k ∗ ( α m ), where the smoothed differences Δ c k ∗ ( α m ) stabilize within a specified tolerance relative to the variability of Δ c k ( α m ), an example of which is plotted in Figure 2. Formally, let τ k represent the standard deviation of the differences Δ c k ∗ ( α m ). The optimal threshold α k ∗ is given by: (3.1) α k ∗ = min m { α m : Δ c k ∗ ( α m ) − Δ c k ∗ ( α m − 1 ) < 0.05 τ k , m > arg max m Δ c k ∗ ( α m ) } .

Figure 2 illustrates this process for determining the optimal threshold for segmenting cell objects from the background. While the first difference Δ c k ( α m ) measures the immediate change in the number of pixels above each threshold, the second difference more precisely captures the stabilization process of these changes. Initially, at extremely small thresholds, almost all pixels are classified as foreground or cell objects, leading to minimal changes in pixel counts between threshold values, shown as A in Figure 2. However, once the threshold starts to intersect the main background distribution, even small increments can trigger substantial variations, resulting in high second differences. As the threshold increases, the second difference reaches its maximum value at B in Figure 2 and then becomes erratic between B and C in Figure 2. This erratic behavior occurs because, in this transitional range, as the predictive mean values near the decision boundary are very similar, even minor changes in the threshold result in disproportionately large shifts in the number of pixels classified as foreground or cell objects. Beyond this range, the second difference begins to stabilize and falls below a predetermined tolerance level show as E in Figure 2, indicating that further increases yield only negligible changes in segmentation. Any further increase in the threshold results in over-segmentation.

Panel (A) in Figure 3 plots the location of different threshold values in the distribution of pixel intensity values. The optimal threshold selected by the algorithm lies on the right side of the mode. Panel (B) of Figure 3 plots a horizontal plane of the optimal threshold value. A threshold value much smaller than this value misclassifies the background as cell objects, whereas a threshold value much larger than the threshold value misses some small cell objects with low intensity peaks.

After estimating the threshold value α k ∗ , we construct a binary image B k ∗ for each cropped image at each pixel: B k ∗ ( x i , j ) = 1 if F k ∗ ( x i , j ) > α k ∗ max ( F k ∗ ) 0 otherwise for k = 1 , … , K ˜. We complete panel (C) in Figure 1 by generating binary cell masks for all sub-images. To address a rare scenario with a higher-than-expected object count, which only happened once for all sub-images, we implemented a re-thresholding step discussed in Section S3 in the Supplementary Material. Finally, each object group is separated and detected via flood fill [4], as described in Section S4 in the Supplementary Material.

After obtaining the binary image, we apply the watershed algorithm to separate cell objects from EBImage package in R [35], which is particularly useful for segmenting objects connecting to each other. The watershed algorithm creates a topological representation of an image based on a distance map, which is derived from the binary image representing foreground objects (cells). The distance map assigns a value to each foreground pixel based on its Euclidean distance to the nearest background pixel; higher values represent larger distances from the background. These distance values act as heights in a topological landscape, providing a basis for the watershed function to simulate water through. The negative distance map is then calculated by taking the negative of the distance map. Examples of the binary image and negative distance map are provided in Section S5 in the Supplementary Material.

The local minima of the negative distance map generally correspond to the centers of cell objects and are used as the starting point for the “flooding” process in the watershed algorithm. The water fills the cell object from these local minima as the negative distance map is dunked into water, so the cell object acts as a catch basin as the map floods. As water from different starting points meets, watershed lines are formed, which delineate boundaries for distinct binary cell objects. Additionally, whenever the water from a catch basin meets a background pixel, a watershed line is formed [55]. These watershed lines separate and create distinct cell objects, and these cell objects can be assigned unique labels.

An example of the watershed algorithm at distinct water levels is provided in Figure 4 with two cells denoted by d 1 and d 2 . The top row of Figure 4 plots the cell images where the negative heights of a slice over the red line are plotted in the bottom row. Figure 4 (A) shows depths for the two local minima for the 1D example are D ( d 1 ) = − 14.14 and D ( d 2 ) = − 14.21. Each height acts as a water level, and the local minima act as seed points for water to flow through into each valley or catch basin. When the water level rises to touch a local minimum, the water is given the same label as the local minimum. As the water level rises, all unidentified pixels at or below the water level with neighboring labeled pixels are given the same label.

Figure 4 (B) shows the water levels at − 13. Each catch basin is partially filled, as the water level is higher than the two local minima. As the water in each valley is still well separated, all emerged pixels in each valley are either labeled corresponding with d 1 or d 2 , with no neighboring pixels containing a different label.

At water level around − 10.2, shown in Figure 4 (C), the water from the two different basins begins to touch. The water from the two catch basins corresponding to d 1 and d 2 meet at the unlabeled pixel p 0 with a position close to 14, plotted in the x-coordinate. The Euclidean distances between p 0 and the local minima are utilized to determine its assignment, and a watershed line is formed. The Euclidean distance between p 0 and d 2 ( ≈ 5.01 units) is smaller than the distance between p 0 and d 1 ( ≈ 5.62 units), so p 0 is grouped with cell 2. The watershed line will serve as the boundary between the two cells.

At the water level 0, the objects are completely segmented, shown in Figure 4 (D). The watershed operation then halts, as this height is used as a stopping point since all cell objects have negative heights in the negative distance map. Additional details of the watershed algorithm, including definitions and edge cases, can be found in [55].

Once all objects are segmented by watershed, the average foreground object pixel count is calculated. To filter out any noise that may have been included as a foreground object, all objects with pixel counts less than or equal to 15% of the average are removed. If an object is touching the image boundary, the threshold count for removal is reduced to 5% of the average, as some objects could be partially cut off during imaging.

We use microscopy images of human dermal fibroblasts (hDFs) in [32] for testing different approaches. More details on the experimental conditions and techniques can be found in Section S7 in the Supplementary Material. These images contain both nuclei and cytoplasm channels, and evaluation of the proposed GP-based segmentation method was performed on both channels.

To evaluate the accuracy of the truth and segmentation results, we compare the area of overlap between the detected object and the true object relative to the overall union area of the two objects [39]. Let g represent a ground truth mask and p a predicted mask. The Intersection over Union (IoU) metric between g and p is defined as: (5.1) IoU ( g , p ) = | g ∩ p | | g ∪ p | . A higher IoU indicates a better match between the predictive mask and true mask of an object, with the perfect match occurring at an IoU of 1. We test several levels of IoU for each method: α = [ 0.5 , 0.55 , 0.6 , 0.65 , 0.7 , 0.75 , 0.8 ]. When α = 0.6, for instance, it means that whenever the IoU defined in Equation (5.1) of a cell is larger than 0.6 by a certain method, this object will be classified as a true positive for this method. The range of thresholds allows us to assess image segmentation performance under distinct precision levels for defining true positives.

The average precision metric ( AP ( α )) at each IoU threshold α ∈ [ 0 , 1 ] is often used for evaluating the correction of image segmentation [50]. For α ∈ [ 0 , 1 ], a True Positive (TP) means that the ground truth mask g is matched to a predicted mask p with IoU ( g , p ) ≥ α. Otherwise, unmatched predicted masks are counted as False Positives (FP), and unmatched ground truth masks are counted as False Negatives (FN). The AP ( α ) is calculated as follows: (5.2) AP ( α ) = TP ( α ) TP ( α ) + FP ( α ) + FN ( α ) .

We compare ImageJ segmentation and the GP-based unsupervised detection algorithm developed in this work. We followed the conventional processing steps in ImageJ to generate the cell masks [46]. All input images are converted to grayscale. To detect the foreground and background pixels, a threshold was set for each image using the default methods. To reduce noises in ImageJ, a despeckle operation, which applies a 3 × 3 median filter, was then applied to each image. A watershed algorithm was used to ensure better separation of connecting cells. Finally, the cell masks were generated. The cell size is selected to be larger than 30 pixels. Then, the binary image and the labeled mask can be generated in ImageJ. A detailed description of each step in this process for ImageJ can be found in Section S8 in the Supplementary Material.

We first use five images of cell nuclei to test different methods. We choose a wide range of image sizes to test methods for images with different sizes. The smallest dimension of these images is 374 × 250, and the largest image is 962 × 1128. The experimental details of these images are shown in Table S1 in the Supplementary Material.

In panel (A) of Figure 9, we plot the AP ( α ) between the GP-based segmentation and ImageJ segmentation methods by solid lines and dashed lines, respectively, for 5 test images of nuclei at distinct threshold α. For any image and threshold level, the solid line is almost always above the dashed line with the same color, as the GP-based method consistently achieves higher AP than ImageJ at any threshold level α. The performance gap becomes more obvious at intermediate IoU thresholds between 0.6 to 0.7, which means that the GP-based method performs better at recovering the overall structure of the object than ImageJ. The AP score of both methods drops at a high IoU threshold, as the annotated object may not be fully accurate, which can affect defining the truth. Panel (B) in Figure 9 shows the distribution of these AP scores from GPs and ImageJ segmentation methods averaged at each threshold. Similar to the results in panel (A), the GP-based method achieves higher median AP scores across all thresholds. The difference between the GP-based method and ImageJ is more pronounced at a high threshold level, such as α = 0.75 or α = 0.8, where the boxes of two colors nearly have no overlap, indicating that GP-based method better capture the fine-scale boundary details than ImageJ.

The true annotated boundaries, those generated by the GP-based method, and those generated by ImageJ of the fifth test image are plotted in Figure 9 (C)-(E). The GP-based image segmentation method generally produces smoother, more continuous boundaries capturing the actual cell nuclei, even in cases where cells were closely clustered. In contrast, the ImageJ segmentation method sometimes produces fragmented boundaries and erroneously segmented individual nuclei into multiple smaller regions, especially in regions where the variations of pixel intensity are large. Additionally, ImageJ fails to detect multiple cells entirely. Further image refinement may be achieved using additional plug-ins or filters from ImageJ. However, this process would require extensive human intervention for trial and error, which would substantially increase the image processing time.

We compare the performance of the GP-based segmentation method and the ImageJ method for five whole-cell images. The smallest image is 602 × 600, and the largest is 1000 × 1000.

Panel (A) of Figure 10 shows the AP with different thresholds by the two methods for five images of whole cells. The variation of AP by both methods for the five images of whole cells is larger than the five images of the nuclei images. This is expected as segmenting whole cells is much more challenging due to the distinct shapes of the whole-cell images and more objects connecting to each other. Nonetheless, the solid line representing the AP for the GP-based method is again almost always above the dashed line representing the AP of the ImageJ segmentation results with the same color for all thresholds, meaning that the GP method consistently achieves higher AP scores than ImageJ. Both methods show much lower precision at a smaller IoU threshold, as correctly identifying irregular boundary shapes in connecting whole cells is extremely hard to capture. Such boundaries also pose challenges for humans to annotate, and for supervised methods to segment, even when annotated data are available.

Panel (B) of Figure 10 shows the distribution of the AP scores between the two methods. Compared to the segmentation results of the cell nuclei, the improvement by the GP-based segmentation methods is more pronounced for images of whole cells. The variation of AP by the GP-based segmentation is larger for whole-cell images than the ones for the images of nuclei, due to the difficulty in identifying the whole-cell shape either by statistical learning methods or by humans.

The ground truth boundaries, the GP-based method generated boundaries, and the ImageJ generated boundaries of the first text image are given in panels (C)-(E) in Figure 10. The GP-based method generally produces boundaries that better fit the actual cell shapes, preserving the irregular shapes of the cell cytoplasm. Although the GP-based method generally better preserves the true cell shapes, cell objects were also often grouped together. In contrast, ImageJ segmentation often produces fragmented boundaries and splits individual cell shapes into multiple smaller regions.