<figure><img src="image_1.png" /></figure>

HAL
open science

# Prediction of Parametric Surfaces for Multi-Object Segmentation in 3D Biological Imaging

Quentin Rapilly, Anaïs Badoual, Pierre Maindron, Guenaelle Bouet, Charles Kervrann

Kervrann

► To cite this version:

Quentin Rapilly, Anaïs Badoual, Pierre Maindron, Guenaelle Bouet, Charles Kervrann. Prediction of Parametric Surfaces for Multi-Object Segmentation in 3D Biological Imaging. SSVM 2025 - 10th International Conference on Scale Space and Variational Methods in Computer Vision, May 2025, Totnes, United Kingdom. pp.255-268, 10.1007/978-3-031-92366-1_20. http://hal-04978619v2

HAL Id: hal-04978619

https://hal.science/hal-04978619v2

Submitted on 16 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Prediction of Parametric Surfaces for Multi-Object Segmentation in 3D Biological Imaging

Quentin RAPILLY¹, Anaïs BADOUAL¹, Pierre MAINDRON²,
Guenaelle BOUET², and Charles KERVRANN¹

¹Centre Inria de l'Université de Rennes, Equipe SAIRPICO, Unité Chimie et Biologie de la Cellule, U1143 INSERM, UMR3666 CNRS, Institut Curie, Université PSL

²Mines Saint Etienne, Université Jean Monnet Saint-Étienne, INSERM, SAINBIOSE U1059, F-42023, SAINT-ÉTIENNE, FRANCE

## Abstract

Multi-object segmentation algorithms are of great interest in a very large range of fields. Deep learning brought major improvements in terms of processing speed or prediction accuracy. Nevertheless, some traditional methods such as active surfaces have features that conventional deep learning methods cannot provide, especially representing the object in a continuous geometrical way and encoding prior information on the shapes to segment. Those features are of particular interest in biology to efficiently segment noisy and poorly resolved data, and then understand the interactions between segmented cells. We introduce NAGINI-3D (N-Active shapes for seGmentING 3D biological Images), a new hybrid segmentation method dedicated to multi-object segmentation of 3D images that combines the efficiency of deep learning and the powerful representation of active surfaces. We evaluate our method on real and synthetic 3D datasets of fluorescence microscopy.

**Keywords**— Multi-object segmentation, Parametric surfaces, Convolutional Neural Networks.

This version of the paper is not the peer-reviewed one. To access it, consult SSVM 2025 conference proceedings at https://doi.org/10.1007/978-3-031-92366-1_20.

## 1 Introduction

The study of biological phenomena requires to capture a lot of observations and produces many microscopy images. That big amount of data can be tedious for experts to analyse. It has led to the development of automatic methods for various tasks, such as counting, classifying or segmenting objects.

More than three decades ago, active contours in 2D, also called "snakes", then generalised to active surfaces in 3D, were widely used for segmentation tasks in biomedical imaging. In recent years, deep learning brought new powerful segmentation methods able to overcome difficulties with snakes, in particular the management of multiple objects. They are also faster and more accurate but cannot represent objects in a continuous way. In this paper, we propose a method that combines the benefits of both paradigms by using a Convolutional Neural Network (CNN) to predict, for each object in an image, the optimal parameters of a surface representing its boundaries.

Snakes are known to be very popular for various segmentation tasks in medical imaging [1, 2, 3]. Initially proposed for 2D image segmentation, parametric snakes consist in a contour parameterized by control points that are roughly initialized to approximate the boundaries of the object of interest. Their positions are then iteratively optimized by minimizing an energy functional, to delineate the object. Their main advantage is that they can represent the object to segment with few parameters, in a continuous and smooth way, that can be used afterwards for morphological and geometrical analysis. This method has been generalized to 3D [4, 5] with active surfaces. However, the need of an initialization and the time consuming optimization made them quite obsolete for many use-case despite their interesting properties. Moreover, most of them are dedicated to single object segmentation and are limited for multi-object problems.

More recently, the revolution of deep learning led to powerful generic methods dedicated to segmentation tasks that produce pixel-wise segmentation masks and deal with multi-object easily [6, 7]. Recent foundation models [8] can segment objects unseen in
their training dataset, as they learn the concept of boundaries between objects. However, these models often require a huge quantity of annotated data to be trained. If microscopy imaging produces big amount of data, few labels are available. It makes generic models less efficient than deep learning methods dedicated to biological images that incorporate geometrical features: StarDist represents objects with star-convex polygon [9] while Cellpose and more recently Omnipose [10, 11] predict gradients used to construct the cells by solving the heat diffusion equation. These methods, dedicated to 2D images, are widely used by biologists because of their accurate multi-object detection and segmentation results. Concerning the segmentation of 3D images, there is no consensus within the community. Some methods are direct 3D extensions of 2D counterparts such as StarDist [12] and Omnipose which provides a reconstruction algorithm with 3D gradient flows [11]. These methods are often criticized because 3D annotated datasets are rare and the trained algorithms often lack of generalization. Other approaches overcome this limitation by applying the 2D algorithm on each slice of a 3D volume and combine the 2D results in the final step [10, 13, 14]. These methods can suffer from spatial incoherence along the third dimension and require a lot of hyper-parameters tuning. Moreover, they are often less accurate than the 3D methods if 3D annotations are available.

Among hybrid methods, SplineDist [15] proposed an approach for 2D images using a CNN to predict parametric contours. Inspired from StarDist 2D, SplineDist offers a way to segment objects in a continuous manner and is not restricted to star-convex shapes. A variant of StarDist 3D, using spherical harmonics instead of star-convex polyhedron to represent object continuously, is presented in [16]. Nevertheless, spherical harmonics still limit the approach to star-convex objects. To our knowledge, the idea of SplineDist has never been extended to 3D before because of several technical barriers explained below. We present a first approach to overcome them.

In this paper, we propose a segmentation method that is not limited to star-convex shapes, inspired from both SplineDist 2D and StarDist 3D, that 1) performs multi-object segmentation on 3D images in a fast manner, taking benefits of the CNN architectures; 2) represents objects in a smooth and continuous way using parametric surfaces; 3) approximates any object homeomorphic to a sphere. In our experiments, we show that our method achieves competitive results with the state-of-the-art algorithms. We validate our approach on a synthetic dataset depicting non star-convex objects and compare our method on a dataset depicting nuclei in fluorescence microscopy images. Our code is available on GitHub at https://github.com/QuentinRapilly/NAGINI-3D/.

## 2 Parametric Spline Surfaces for Object Representation

We represent each object with a parametric surface:

$$ \forall (u, v) \in [0, 1]^2, \quad s(u, v) = \sum_{k=-\infty}^{\infty} \sum_{l=-\infty}^{\infty} c_{k,l} \varphi_1(M_1 u - k) \varphi_2(M_2 v - l), \quad (1) $$

where {$c_{k,l}$}$_{k \in \mathbb{Z}}$ ∈ ($\mathbb{R}^3$)$^{Z \times Z}$ are the 3D control points describing the shape. In (1), $M_1$ and $M_2$ are two parameters determining the complexity of the surface. The higher they are, the more accurately the surface can represent local variations. Moreover, $\varphi_1$ and $\varphi_2$ are the basis functions, whose choice allows to enforce various properties to the surface (smoothness, reproduction properties, ...). Here, we use the parameterization as described in [5] for surfaces with a spherical topology. The functions $\varphi_1 = \phi_{M_1,per}$ and $\varphi_2 = \phi_{M_2}$ allows to reproduce ellipsoids using (1), where the function $\phi_{M,per}(t) = \sum_{n=-\infty}^{\infty} \phi_M(t - Mn)$ is the M periodic version of $\phi_M$, an exponential B-Spline defined as:

$$ \phi_M : t \to \frac{1}{1 - \cos \frac{2\pi}{M}} \left[ \mathbb{1}_{[0, \frac{1}{2}]}(|t|) \left( \cos \frac{2\pi|t|}{M} \cos \frac{\pi}{M} - \cos \frac{2\pi}{M} \right) + \mathbb{1}_{[\frac{1}{2}, \frac{3}{2}]}(|t|) \left( \sin \frac{\pi(\frac{3}{2} - |t|)}{M} \right)^2 \right]. $$

Using $\varphi_1$ periodicity and $\varphi_2$ compact support, we rewrite (1) in the form:

$$ \forall (u, v) \in [0, 1]^2, \quad s(u, v) = \sum_{k=0}^{M_1-1} \sum_{l=-1}^{M_2+1} c_{k,l} \varphi_1(M_1 u - k) \varphi_2(M_2 v - l). \quad (2) $$

Note that, for surfaces homeomorphic to a sphere, the variables u and v can be viewed as the longitude and the latitude parameters of the surface.
In (2), to guarantee that the surface is closed and smooth, the control points have to satisfy the following constraints [5]:

$$c_{k,-1} = c_{k,1} + \frac{c_{M_1}[k]T_{1,N} + s_{M_1}[k]T_{2,N}}{M_2\varphi'_2(1)}, \quad (3a)$$

$$c_{k,0} = \frac{c_N}{\varphi_2(0)} - \frac{\varphi_2(1)(c_{k,-1} + c_{k,1})}{\varphi_2(0)}, \quad (3b)$$

$$c_{k,M_2+1} = c_{k,M_2-1} - \frac{c_{M_1}[k]T_{1,S} + s_{M_1}[k]T_{2,S}}{M_2\varphi'_2(1)}, \quad (3c)$$

$$c_{k,M_2} = \frac{c_S}{\varphi_2(0)} - \frac{\varphi_2(1)(c_{k,M_2-1} + c_{k,M_2+1})}{\varphi_2(0)}, \quad (3d)$$

with $c_{M_1}[n] = \Gamma_{M_1} \cos \frac{2\pi n}{M_1}$, $s_{M_1}[n] = \Gamma_{M_1} \sin \frac{2\pi n}{M_1}$ and $\Gamma_{M_1} = \frac{2(1-\cos \frac{2\pi}{M_1})}{\cos \frac{\pi}{M_1}-\cos \frac{3\pi}{M_1}}$. Here, $c_S$, $c_N$ denote the south and north poles of the surface respectively, and the vectors $T_{1,N}$, $T_{2,N} \in \mathbb{R}^3$, ($T_{1,S}$, $T_{2,S} \in \mathbb{R}^3$, respectively) characterize the tangent plane to the north pole (south, respectively).

Finally, the set of $I = M_1(M_2 - 1) + 6$ parameters $\{f_i\}_{1 \le i \le I} \in (\mathbb{R}^3)^I$ provides a way to represent any object homeomorphic to a sphere by a closed and smooth parametric surface:

$$\{f_i\}_{1 \le i \le I} = \{c_N, c_S, T_{1,N}, T_{2,N}, T_{1,S}, T_{2,S}, c_{0,1}, \dots, c_{M_1-1, M_2-1}\}. \quad (4)$$

It contains a subpart of the control points $\{c_{k,l}\}_{\substack{0 \le k \le M_1-1 \\ 1 \le l \le M_2-1}}$. The missing ones are constructed using (3). An example of parameterized surface is shown on Fig 1a.

# 3 A CNN for the Prediction of the Surface Parameters

Our idea is to train a CNN to predict for each object in an image, an optimal set of parameters $\{f_i\}_{1 \le i \le I}$ that describes its boundaries in a single forward pass.

The CNN architecture we consider is a U-Net [7]. The input of the network is a 3D image of shape ($C_{in}$, $H$, $W$, $D$), where $H$, $W$ and $D$ denote the height, width and depth of the image, respectively, and $C_{in}$ is the number of channels. A channel in microscopy often represents the intensity of the signal representing a dedicated molecule. The network predicts a volume of shape ($C_{out}$, $H$, $W$, $D$). For each voxel of the input image $x = (x, y, z) \in \Omega$ with $\Omega = \{1, \dots, H\} \times \{1, \dots, W\} \times \{1, \dots, D\}$, the network outputs $C_{out} = 3I + 1$ channels. The 3I first channels correspond to the 3 spherical coordinates of the $I$ estimated parameters $\{\hat{f}_{x,i}\}$ of a surface delineating the object to which $x$ belongs. The spherical coordinates are predicted in a frame centered on $x$. The last channel is the probability $\hat{p}(x)$ of $x \in \Omega$ to be the barycenter of an object. This channel is used in the object detection step during the inference. It allows to extract a single location in each object. At each location, we use the predicted parameters to construct the surface representing the corresponding object. More details are given in Section 3.2.1.

## 3.1 Training Description

The training pipeline (see Fig 2) consists in comparing the predicted channels to the ground-truth (GT). To compare the predicted probability map with a GT one, we use the Binary Cross Entropy (BCE). However, evaluating how well a predicted surface describes the object of interest is not a trivial task. One cannot just compare the predicted control points with "optimal ones", as very different sets of parameters can describe nearly identical surfaces. Instead of comparing the parameters, we compare the surface generated with (2) to the GT surface computed from the boundaries of a binary mask.

The idea, inspired from [15] (in 2D), is to sample the two surfaces and compare the resulting points clouds. However, the translation from 2D contours to surfaces is still difficult for the two following reasons: 1) Sampling the predicted surface cannot be performed by evenly sampling the variables $u$ and $v$ definition domains, as it would result in a concentration of points in the vicinity of the poles (Fig 1b); 2) The points sampled on a contour are ordered (by following it) which leads to a very simple formulation for the training loss [17]. On surfaces, the sampled points are unordered.

### 3.1.1 Surface Sampling Procedure

We aim at sampling points, evenly distributed, on a surface described by (2), using parameters predicted by the CNN. Evenly partitioning the definition domains of variables
<figure><img src="image_3.png" /><figcaption>Figure 1: (a) Parametric surface (green arrows: normal at poles $T_{1,N}, T_{2,N}, T_{1,S}, T_{2,S}$; green dots: control points {$c_{k,l}$}, see (4)). (b) Sampling method evenly dividing $u$ and $v$ domains. (c) Sampling method using the Fibonacci lattice.</figcaption></figure>

u and v: $(u_p)_{0 \le p \le P-1}$ and $(v_q)_{0 \le q \le Q-1}$ with $u_p = \frac{p}{P}$ and $v_q = \frac{q}{Q}$ would result in oversampled regions around the poles and undersampled regions around the equator of the surface [18]. A solution consists in generating a Fibonacci lattice containing $P = 2N + 1$ points ($N \in \mathbb{N}$), giving the latitudes and longitudes $(u_i, v_i)_{-N \le i \le N}$ to use to sample the surface: $\forall i \in \{-N, ..., N\}, u_i = \text{frac}(i \cdot G)$ and $v_i = \frac{1}{\pi} \cdot \arcsin(\frac{2 \cdot i}{P}) + 0.5$, where $\text{frac}: x \to x - \lfloor x \rfloor$ denotes the fractional part function and $G = \frac{\sqrt{5-1}}{2}$. The two samplings methods are shown on Fig 1b and 1c. By substituting these values in (2) we obtain P points $\hat{S} = (\hat{S}_i)_{-N \le i \le N}$ evenly sampling the surface:

$$ \forall i \in \{-N, ..., N\}, \quad \hat{S}_i = \sum_{k=0}^{M_1-1} \sum_{l=-1}^{M_2+1} c_{k,l} \varphi_1(M_1 u_i - k) \varphi_2(M_2 v_i - l). \quad (5) $$

For fixed values of $P, M_1$ and $M_2$, we pre-compute and store the terms $w_{i,k,l}^{P,M_1,M_2} = \varphi_1(M_1 u_i - k)\varphi_2(M_2 v_i - l)$, $\forall (i, k, l) \in \{-N, ..., N\} \times \{0, ..., M_1 - 1\} \times \{-1, ..., M_2 + 1\}$, which are identical for all surfaces. This leads to an efficient way to sample any surface by just computing a weighted sum of the control points:

$$ \forall i \in \{-N, ..., N\}, \quad \hat{S}_i = \sum_{k=0}^{M_1-1} \sum_{l=-1}^{M_2+1} c_{k,l} w_{i,k,l}^{P,M_1,M_2}. \quad (6) $$

This formulation makes the computation of multiple sampling in the mean time very fast, as it can be parallelized on GPU architectures.

To sample a GT surface, we first compute the boundary of the segmentation mask (using a dilation deformation). Then we apply a "Farthest Point Sampling" (FPS) algorithm [19], that samples a surface evenly. It consists in choosing one initial point on the surface then iteratively adding to the sampling the farthest point on the surface (using geodesic distance) to all the others already added. The process is stopped when the desired number of points is reached. That leads to a set of P points $S = (S_i)_{1 \le i \le P}$ evenly distributed on the GT surface. The sampling of each GT surface is computed once and stored.

### 3.1.2 Surface matching

The next issue to deal with is to perform matching between the samplings of the GT and predicted surfaces. According to [17], the ordering of the points on sampled contours and cyclic permutations allows to recover the best matching. Nevertheless, in 3D, the points are unordered on the surface and there is no such trick to speed up the computation of a distance to our knowledge. To overcome this difficulty, the Sliced Wasserstein Distance (SW) [20], has been investigated to compare points clouds [21]. SW is more representative than the chamfer distance and faster to compute than the Wasserstein Distance (WD). In our case, SW of order $p$ between $\hat{S}$ and $S$ is defined as:

$$ SW_p(\hat{S}, S) = \int_{\theta \in \mathcal{U}} W_p(\pi_{\theta}(\hat{S}), \pi_{\theta}(S)) d\theta, \quad (7) $$

where $\mathcal{U}$ is the unit sphere, $\pi_{\theta}(S) = \{\theta^T S_i, S_i \in S\}$ is the set of projections of each element of $S$ on the 1D subspace of $\mathbb{R}^3$ directed by $\theta$. In our framework, $W_p(X, Y)$ is the 1D WD between two discrete one-dimensional distributions. It consists in sorting the two sets in increasing order and comparing the elements two by two: $W_p(X, Y) = (\frac{1}{P} \sum_{i=1}^{P} |x_{\sigma_X(i)} - y_{\sigma_Y(i)}|^p)^{1/p}$, where $\sigma_X$ and $\sigma_Y$ are the permutations operators sorting X and Y in increasing order.

According to [21], the integral in (7) can be approximated by a Monte-Carlo procedure that randomly picks $K \ge 100$ directions $\theta$. Finally, we can rewrite the distance:

$$ SW_p(\hat{S}, S) = \frac{1}{K} \sum_{\theta \in \mathcal{U}_K} W_p(\pi_{\theta}(\hat{S}), \pi_{\theta}(S)), \quad (8) $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Training loop of the network.</figcaption></figure>

where $\cup_K$ contains $K$ random elements, uniformly sampled in $\cup$. In the following, we use $p = 1$ and note $SW$ the $SW_1$ of order 1.

### 3.1.3 Computation of the GT Probability Map

For each voxel $x \in \Omega$, we define the GT probability $p(x)$ as follow: $p(x) = 0$ if $x$ belongs to the background. If $x$ is inside an object of index $j$ ($x \in \Omega_j$), we define $d_j(x) = \|b_j - x\|_2$ the distance to the barycenter $b_j$ of the object $\Omega_j$. We note $D(x) = \max_{y \in \Omega_j} d_j(y) - d_j(x)$ and then $p(x) = \frac{D(x)}{\max_{y \in \Omega_j} D(y)}$ its normalized version. With this definition, $p(x)$ is positive, increases as $x$ get closer to the barycenter of an object of interest and is maximum equal to 1 if $\exists j, x = b_j$. For particular shapes, when the barycenter does not belong to the interior of the object, it is maximum for the point in this object that is the closest to the barycenter.

### 3.1.4 Training Loss

We define the training loss for each voxel $x$ of an image $\Omega$ as:

$$L(p(x), \hat{p}(x), S_x, \hat{S}_x, \{\hat{f}_{x,i}\}) = BCE(p(x), \hat{p}(x)) + \lambda_s L_s(p(x), S_x, \hat{S}_x, \{\hat{f}_{x,i}\}), \quad (9)$$

where the parameter $\lambda_s$ weights the importance given to the surface loss $L_s$. The surface loss compares the two samplings $\hat{S}_x$ and $S_x$ using $SW$, if $x$ belongs to an object ($x \in \Omega_{fg}$). For the background voxels ($x \in \Omega_{bg}$) it minimizes the norm of the predicted parameters. This acts as a regularization term to penalize the values of the parameters to be too large (leading to unrealistic shapes). The surface loss is defined as:

$$L_s(p(x), S_x, \hat{S}_x, \{\hat{f}_{x,i}\}) = p(x) \cdot \mathbb{1}_{\Omega_{fg}}(x) \cdot SW(\hat{S}_x, S_x) + \lambda_{reg} \cdot \mathbb{1}_{\Omega_{bg}}(x) \cdot \sum_i \|\hat{f}_{x,i}\|_2, \quad (10)$$

where $\mathbb{1}$ denotes the indicator function. The SW is weighted by the GT probability to give more importance to the voxels that are more likely to be a voxel of interest. The parameter $\lambda_{reg}$ controls the importance of the regularization term.

## 3.2 Inference Procedure for Segmentation

### 3.2.1 Surface Prediction

During the inference step, the goal is to predict one surface for each object on a given input image $\Omega$. To determine at which positions we need to construct parametric surfaces we extract the local maxima ($x_q$)$_{1 \le q \le Q}$ of $\hat{p}$ having a probability greater than a chosen threshold $\gamma > 0$. These voxels are supposed to represent the barycenters (or the closest voxel to them) in the objects of interest. We generate the surface parameterized by $\{\hat{f}_{x_q,i}\}_{1 \le i \le I}$ for each selected location $x_q$. Finally, we process to a Non Maximum Suppression (NMS) that consists in evaluating the Intersection-over-Union (IoU) between each pair of surface and removing potential duplicated surfaces representing the same object. Our inference pipeline is illustrated in Fig 3.

### 3.2.2 Refinement of Predicted Surfaces

We observed experimentally that the network is much easier to train with small values of $M_1$ and $M_2$ in (2), mainly because it avoids unrealistic self-intersecting shapes. In a refinement step, we decided to add an optional post-processing algorithm: we up-sample
<figure><img src="image_5.png" /><figcaption>Figure 3: Inference pipeline applied to generate multi-object segmentation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Accuracy curves $Acc(\sigma)$ depending on the IoU threshold $\sigma$ for each method.</figcaption></figure>

the number of control points of the predicted shapes and then apply a classical snake optimization process with an energy based on image gradients just before the NMS step. We can then train a network for small values of $M_1$ and $M_2$, and still segment more complex shapes (with local variations or non-convex parts) that require a higher number of control points to be represented. Even with this step, segmenting an image that contains hundreds of objects takes only a few seconds.

## 3.3 Settings of Model Hyper Parameters

The U-Net implementation we used is the very versatile one "General U-Net" described in [22]. It has 3 layers. Each layer is a double 3 × 3 convolution, with batch-normalisation and ReLU activation function. The number of channels of the first layer is 32. At each down-sampling step, the number of channels is multiplied by a factor 2 and each dimension of the image is divided by a factor 2. The parameters used are $\lambda_s = 1$ and $\lambda_{reg} = 10^{-3}$. We trained our model over 6000 epochs, with a batch size of 1, using the AdamW optimizer with a learning rate of $10^{-4}$. We applied data augmentation on the images using random flips along each axis. The number of sampled points is set to $P = 151$. We chose small values of the surface parameters $M_1 = 4$ and $M_2 = 2$ ($I = 10$) as explained in Section 3.2.1. For the experiments, we present the results with the up-sampling of control points $M_1 \leftarrow 2M_1$ and $M_2 \leftarrow 2M_2$ ($I = 30$) and the snake refinement step, denoted as **Our+Opt**. The results without refinement are denoted as **Our**. The algorithms are trained on a Nvidia A100 GPU.

The probability threshold $\gamma$ used during the local maxima extraction and the IoU threshold used for the NMS step are selected automatically to maximize $\int_0^1 Acc(\sigma)d\sigma$ (see in Section 4), on the validation set, at the end of the training.

# 4 Experiments

In this section, we present the results obtained on a dataset of real microscopy images depicting nuclei and a dataset of synthetic images: **BONES** is a dataset of 28 images of fluorescence microscopy representing nuclei of bone cells cultivated on furrows of various angles to evaluate the influence of support on cell development. We split the set into 18 for training, 3 for validation and 7 for testing. Each image (80 – 320 × 180 × 310 voxels) (the first dimension depends of the furrow angle) contains approximately 240 objects. The GT masks were segmented automatically and curated by an expert. **CAPS** is a synthetic dataset we created, representing spherical caps (portion of a sphere cut off by a plane). The geometry of these objects could be assimilated to those of erythrocytes (red blood cells). We generated caps of various size and thickness randomly cut off to represent star-convex and non star-convex objects. The proportion of non-star-convex objects is driven by the ranges of possible positions of the cut-off plane (see Fig 6a) and
<table><caption>Table 1: Detection and segmentation scores of tested methods. The recall and precision are both computed with TP(0), FP(0) and FN(0) (i.e. σ = 0) to assess the quality of the detection. The mean IoU is given among the TP to assess the quality of the segmentation for the objects correctly detected.</caption><thead><tr><th></th><th>Star3D</th><th>Omni3D</th><th>Sti+CP</th><th>Sti+Star</th><th>Useg</th><th>CP2.5D</th><th>Our</th><th>Our+Opt</th></tr></thead><tbody><tr><td><strong>BONES</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Recall</td><td>0.965</td><td><strong>0.984</strong></td><td>0.930</td><td>0.976</td><td>0.888</td><td>0.955</td><td>0.947</td><td>0.939</td></tr><tr><td>Precision</td><td>0.882</td><td>0.356</td><td>0.755</td><td>0.593</td><td>0.325</td><td>0.809</td><td><strong>0.884</strong></td><td><strong>0.884</strong></td></tr><tr><td>Mean IoU</td><td>0.793</td><td>0.717</td><td>0.511</td><td>0.555</td><td>0.313</td><td>0.584</td><td>0.606</td><td><strong>0.796</strong></td></tr><tr><td><strong>CAPS</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Recall</td><td><strong>1.0</strong></td><td><strong>1.0</strong></td><td>0.979</td><td>0.983</td><td>0.974</td><td>0.986</td><td><strong>1.0</strong></td><td>0.999</td></tr><tr><td>Precision</td><td>0.958</td><td>0.162</td><td>0.905</td><td>0.729</td><td>0.996</td><td>0.994</td><td><strong>1.0</strong></td><td><strong>1.0</strong></td></tr><tr><td>Mean IoU</td><td>0.779</td><td>0.688</td><td>0.764</td><td>0.692</td><td>0.592</td><td>0.776</td><td>0.678</td><td><strong>0.858</strong></td></tr></tbody></table>

of possible thicknesses. The thicker the sphere is and the smaller part of the sphere the cut-off plane keeps, the more likely the corresponding cap will be star-convex. We added blur and Poisson-Gaussian noise to simulate fluorescence microscopy images. We generated 10 training images, 3 for validation, 5 for testing. Each image (200 × 200 × 200 voxels) contains 220 objects.

We compared our method with **direct 3D methods**: StarDist 3D denoted **Star3D** with 96 radii as in [12], and Omnipose 3D [11] denoted **Omni3D**. We trained each method on the dataset considered prior to evaluation. We also compared with **aggregation algorithms** using pre-trained 2D models: U-Segment3D [14] denoted **Useg** provides an implementation using pre-trained 2D Cellpose models, Cellstitch [13] which can easily be linked after any 2D segmentation models, we applied it with StarDist 2D (**Sti+Star**) and Cellpose 2D (**Sti+CP**). For StarDist we used the pre-trained model dedicated to fluorescence. For Cellpose we tried the nuclei and cytoplasm models and kept the best for each use case. We also evaluated the aggregation method built into Cellpose, that have been improved in version 3.1.0, denoted **CP2.5D**. For each method we tuned the hyper-parameters to get the best segmentation results as possible.

We assessed the tested algorithms using the method depicted in [12]. We find matching between the predicted and GT cells by solving an assignment problem maximizing the sum of the IoU of all matches. A GT cell with no corresponding prediction is considered a false-negative (FN) and a predicted cell with no corresponding GT is considered a false-positive (FP). A predicted cell having a corresponding GT is considered a true-positive (TP) prediction and the corresponding IoU is computed. For each IoU threshold $\sigma \in [0, 1]$, TP($\sigma$) is the number of TP with an IoU greater or equal than $\sigma$. FN($\sigma$) (FP($\sigma$), respectively) is the number of GTs (predicted cells, respectively) with no matching prediction (GT, respectively) or in a matching with an IoU smaller than $\sigma$. The evaluation metric that we consider is the accuracy for a given $\sigma$: $Acc(\sigma) = \frac{TP(\sigma)}{TP(\sigma)+FN(\sigma)+FP(\sigma)}$. Fig 4 shows the accuracy curve for each method and dataset. In Table 1, we reported several scores assessing the quality of detection and segmentation.

On BONES, we note that the 3D methods outperform the aggregation algorithms in term of both detection and segmentation. The only exception is **Omni3D**, that predicts very accurate masks but tends to predict residual groups of pixels not merged with other labels or nuclei split into multiple labels, resulting in many FP as shown by its precision score. The best methods for this dataset are **Star3D** and **Our+Opt** that have very similar performances in detection and segmentation. On CAPS, the difference between 3D methods and aggregation algorithms is less noticeable: **Useg** gets detection scores similar to the best 3D methods; **Sti+CP** perform really well in term of segmentation accuracy; **CP2.5D** achieves detection and segmentation scores comparable to 3D methods. Unlike **Star3D**, whose performances drop to that of the best aggregation algorithms, **Our+Opt** still outperforms all the other methods in terms of detection and segmentation. This difference can be explained by the presence of non-star-convex objects in this dataset: **Star3D** cannot represents such objects and split them into several parts (Fig 5), leading to FP predictions and under-segmented objects. In our framework, surfaces are not limited to star-convex shapes and can represent all the objects of the set, allowing **Our+Opt** to get almost perfect detection scores and the best mean IoU by a noticeable margin.

To demonstrate that our method performs equivalently on star-convex and non-star-convex objects, we evaluated it on images with a higher proportion of non-star-convex objects, driven by the range of possible positions of the cut-off plane when creating the caps (see Fig 6a). We thus, generated two other test sets of 5 images, one with significantly more non-star-convex objects but still an important part of star-convex objects (denoted "medium prop. set") and one with almost no star-convex object (denoted "strong prop. set"). On Fig 6, we see that there is a drop in performance for almost every methods
<figure><img src="image_7.png" /><figcaption>Figure 5: Example of **Star3D** splitting non star-convex objects in two. **Our+Opt** predicts a single object. The first row shows the 3D volumes and the second row 2D slices.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: (a) Cap (in red) generated by the cut-off plane (dashed line) randomly positioned in the chosen blue range; (b) and (c) Evolution of accuracy according to the proportion of non-star-convex objects in testing images.</figcaption></figure>

as the proportion of non-star-convex objects increases, especially for **Star3D**. However, **Our**, **CP2.5D** and **Our+Opt** show very little loss of performance.

# 5 Conclusion

We presented an approach to train a CNN to predict parametric surfaces allowing to perform multi-object segmentation on 3D images in a fast way. Our modeling represents any closed shape continuously. We demonstrated that it outperforms aggregation algorithms based on 2D models and is competitive against other direct 3D algorithms, especially on images depicting complicated shapes such as non star-convex objects.

## Author contributions

Q.R. designed the method, implemented the code and conducted the experiments. P.M. and G.B. provided the fluorescence images and the corresponding labels. Q.R. generated the synthetic images and labels. Q.R., A.B and C.K. wrote the article. A.B. and C.K. supervised the research.

## Acknowledgment

Many thanks to the BioDac project (EPSRC grant EP/Y008715/1) for access to computing resources.

# References

[1] M. Kass, A. Witkin, and D. Terzopoulos, "Snakes: Active contour models," *Int J. Comp. Vision*, vol. 1, no. 4, pp. 321–331, 1988.

[2] A. Blake and M. Isard, *Active Contours: The Application of Techniques From Graphics, Vision, Control Theory and Statistics To Visual Tracking of Shapes In Motion*. Springer Science & Business Media, 2012.

[3] A. Badoual, D. Schmitter, V. Uhlmann, and M. Unser, "Multiresolution Subdivision Snakes," *IEEE Trans. Image Process.*, vol. 26, no. 3, pp. 1188–1201, 2016.

[4] A. Dufour, R. Thibeaux, E. Labruyere, N. Guillen, and J.-C. Olivo-Marin, "3-D active meshes: fast discrete deformable models for cell tracking in 3-D time-lapse microscopy," *IEEE Trans. Image Process.*, vol. 20, no. 7, pp. 1925–1937, 2010.
[5] R. Delgado-Gonzalo, N. Chenouard, and M. Unser, "Spline-based deforming ellipsoids for interactive 3D bioimage segmentation," *IEEE Trans. Image Process.*, vol. 22, no. 10, pp. 3926-3940, 2013.

[6] K. He, G. Gkioxari, P. Dollar, and R. Girshick, "Mask R-CNN," in *IEEE Int. Conf. Comp. Vision (ICCV)*, Oct 2017.

[7] O. Ronneberger, P. Fischer, and T. Brox, "U-Net: Convolutional Networks for Biomedical Image Segmentation," in *Med. Image Comp. Comp-Assist. Interv. (MICCAI)*, pp. 234-241, 2015.

[8] A. Kirillov, E. Mintun, N. Ravi, H. Mao, and et al, "Segment Anything," in *IEEE Int. Conf. Comp. Vision (ICCV)*, pp. 4015-4026, October 2023.

[9] U. Schmidt, M. Weigert, C. Broaddus, and G. Myers, "Cell detection with star-convex polygons," in *Med. Image Comp. Comp-Assist. Interv. - (MICCAI)*, pp. 265-273, 2018.

[10] C. Stringer, T. Wang, M. Michaelos, and M. Pachitariu, "Cellpose: a generalist algorithm for cellular segmentation," *Nat. Meth.*, vol. 18, no. 1, pp. 100-106, 2021.

[11] K. J. Cutler, C. Stringer, T. W. Lo, L. Rappez, N. Stroustrup, S. Brook Peterson, P. A. Wiggins, and J. D. Mougous, "Omnipose: a high-precision morphology-independent solution for bacterial cell segmentation," *Nat. Meth.*, vol. 19, no. 11, pp. 1438-1448, 2022.

[12] M. Weigert, U. Schmidt, R. Haase, K. Sugawara, and G. Myers, "Star-convex polyhedra for 3D object detection and segmentation in microscopy," in *IEEE Winter Conf. on App. of Comp. Vision (WACV)*, pp. 3655-3662, 2020.

[13] Y. Liu, Y. Jin, E. Azizi, and A. J. Blumberg, "Cellstitch: 3D cellular anisotropic image segmentation via optimal transport," *BMC Bioinformatics*, vol. 24, no. 1, p. 480, 2023.

[14] F. Y. Zhou, C. Yapp, Z. Shang, S. Daetwyler, Z. Marin, M. T. Islam, B. Nanes, E. Jenkins, G. M. Gihana, B.-J. Chang, et al., "A general algorithm for consensus 3D cell segmentation from 2D segmented stacks," *BioRxiv*, 2024.

[15] S. Mandal and V. Uhlmann, "Splinedist: Automated cell segmentation with spline curves," in *IEEE Int. Symp. Biomedical Imaging (ISBI)*, pp. 1082-1086, 2021.

[16] D. Eschweiler, M. Rethwisch, S. Koppers, and J. Stegmaier, "Spherical harmonics for shape-constrained 3D cell segmentation," in *IEEE Int. Symp. Biomedical Imaging (ISBI)*, pp. 792-796, 2021.

[17] S. Mandal and V. Uhlmann, "A learning-based formulation of parametric curve fitting for bioimage analysis," in *Num. Math. and Adv. App. (ENUMATH): Eur. Conf., Egmond aan Zee, The Netherlands*, pp. 1031-1038, 2021.

[18] Á. González, "Measurement of areas on a sphere using fibonacci and latitude-longitude lattices," *Math. Geosciences*, vol. 42, pp. 49-64, 2010.

[19] G. Peyré and L. D. Cohen, *Geodesic Methods for Shape and Surface Processing*, pp. 29-56. Dordrecht: Springer Netherlands, 2009.

[20] N. Bonneel, J. Rabin, G. Peyré, and H. Pfister, "Sliced and radon wasserstein barycenters of measures," *J. Math. Imaging Vision*, vol. 51, pp. 22-45, 2015.

[21] T. Nguyen, Q.-H. Pham, T. Le, T. Pham, N. Ho, and B.-S. Hua, "Point-set distances for learning representations of 3D point clouds," in *IEEE Int. Conf. Comp. Vision (ICCV)*, pp. 10478-10487, 2021.

[22] E. Meunier and P. Bouthemy, "Unsupervised space-time network for temporally-consistent segmentation of multiple motions," in *IEEE/CVF Conf. Comp. Vision Pattern Recogn.*, pp. 22139-22148, 2023.