# Surface preserving resampling of labeled images

Grégoire Malandain, Patrick Lemaire

► To cite this version:

Grégoire Malandain, Patrick Lemaire. Surface preserving resampling of labeled images. 2025 IEEE 22nd International Symposium on Biomedical Imaging (ISBI), Apr 2025, Houston, United States. pp.1-4, 10.1109/ISBI60581.2025.10980702. hal-05065550

HAL Id: hal-05065550

https://hal.science/hal-05065550v1

Submitted on 13 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SURFACE PRESERVING RESAMPLING OF LABELED IMAGES

Grégoire Malandain* Patrick Lemaire†

* Université Côte d'Azur, Inria, CNRS, I3S, France

† CRBM, Université de Montpellier, CNRS, France

ABSTRACT

Thanks to the development of microscopy techniques, the acquisi-
tion of temporal series of 3D images is becoming a standard for the
study of evolving phenomena (e.g. developmental biology). In most
cases, live samples are moving/growing during long-term imaging,
therefore it is desirable to compensate for this global 3D motion for
both a more comfortable visualization and analysis. It is implic-
itly assumed that the quantitative properties of the resampled series
are similar (if not equal) to those of the original one. Among these
properties, the surface measurement is quite important since it helps
predicting cell behavior, fate, or to compute symmetry axis. In the
present work, we demonstrated that surface estimation of segmented
objects may not be preserved by the nearest neighbor interpolation,
the usual technique to resampled label images, and proposed a novel
interpolation technique, dedicated for multi-labeled images, that pre-
serves the surface estimation.

***Index Terms***— Surface area estimation, nearest neighbor inter-
polation, labeled image resampling

1. INTRODUCTION

Recent advances in microscopy allow now the acquisition of 3D liv-
ing organs or organisms [1] at a cellular level and with a temporal
frequency that allows to capture cell division, hence to build lin-
eages. Eg, in [2], confocal microscopy was used to image devel-
oping flower meristems; in [3], light-sheet microscopy allowed to
follow the development of ascidian embryos. Subsequent segmen-
tation of such temporal series of 3D images, that yields 3D labeled
images (where each label represent a segmented object), allowed to
extract information about the morphogenesis or the embryogenesis.

Because of the natural growth [2], or because of the un-
avoidable displacement of the imaged object within the viscous
or liquid medium wherein it lies [3], it is mandatory to compute
the transformations that appear between successive acquisitions.
Not only it eases the lineage calculation [3], but it also allows
to compensate for the global motion along the temporal series,
making the visual inspection of the data more efficient. Mor-
phonet [4] offers a web-based visualization tool that permits to
manipulate meshes built from segmented temporal series. The in-
terested reader will find there examples of developing meristems
(eg morphonet.org/dataset/1152) or growing ascidian
embryos (eg morphonet.org/dataset/1095) where the pre-
sented meshes have been computed after motion compensation.

Once every couple of successive acquisitions has been co-
registered, with proven techniques [5] or dedicated ones [6], the
obtained transformations are compounded to provide a transforma-
tions with respect to a reference (eg, the first image of the temporal

series), and all images are resampled in the geometry of this refer-
ence. Acquisitions (ie grey-level images) can be resampled using
(tri-)linear or more complex interpolation schemes (spline, etc)
while the nearest neighbor interpolation scheme seems the most
appropriate for segmented images [7] since it won't create inter-
mediary values between two adjacent objects. By doing so, it is
implicitly assumed that the morphometric characteristics of the seg-
mented objects are preserved by the resampling operation. However,
it appeared that this assumption does not stand for surface measure-
ment, hence the need for a surface-preserving interpolation scheme
for labeled images.

2. SURFACE ESTIMATION IN LABELED IMAGES

Measuring surfaces in 3D images has many application. It may give
access to the microstructure of snow [8]; the leaf area [9] helps un-
derstanding the plant status; cell contact surfaces helps predicting the
geometry of cell division [10], or cell communication [3], or allows
to compute the symmetry axis of embryos [11].

3D digital images are made of voxels. The segmentation opera-
tion groups voxels into individualized objects. The count of voxels
of an object is a good estimate of its volume, while the count of voxel
faces (surfels) is definitively a poor estimate of its surface, hence the
need for dedicated surface estimation methods.

On the one hand, most existing methods are defined for closed
surfaces. Crofton formula [12], which integrates the number of in-
tersections between lines and the surface to be computed, can only
estimated closed surfaces. Integrating the gradient over the sur-
fels defining a discrete surface gives also a good estimate of the
(closed) surfaces [13, 14, 8]. When it comes to grey-level images,
surface may be estimated from the mesh issued from a threshold (eg
marching cubes), it can be extended by binary images with an ad-
ditional smoothing [8]. These methods may require some effort to
be adapted to open surfaces, as contact surfaces between objects in a
multi-labeled image. On the other hand, the use of local configura-
tions (associated with optimized weights) is can trivially be adapted
to open surface estimation. Optimal local area estimations associ-
ated to 3x3x3 voxels configurations are computed in [15, 16], while
2x2x2 voxels configurations, identical to marching cube configura-
tions, are considered in [17, 18, 19]. These articles used different
optimization criteria and thus get different weights (Tab. 1). The
total surface results from the summation of the local areas.

The 2x2x2 neighborhood being symmetric with respect to fore-
ground and background, the estimate does not change if we apply it
to the complementary image and this ensures a quasi-equality of a
contact surface between two adjacent objects. For this last reason,
we will only assess [17, 18, 19] in the following. In [19], optimized
weights are computed for only 5 types of configurations and depend
on two parameters, *r* and *s*: those parameters were set to *r* = 0.5
and *s* = 1.704, and the weights associated to other configurations
<table><caption>**Table 1.** Local area associated to the different configurations for [17, 18, 19]. For [19], some coefficients depends on 2 parameters $r \in [0, 1]$ and $s \in [1.663, 1.745]$; we choose $r = 0.5$ and $s = \frac{1.663+1.745}{2} = 1.704$. Informative configurations are in solid line.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th></tr></thead><tbody><tr><td></td><td>□</td><td>□</td><td>□</td><td>□</td><td>□</td><td>□</td><td>□</td><td>□</td><td>□</td><td>□</td><td>□</td><td>□</td><td>□</td></tr><tr><td>[17]</td><td>0.636</td><td>0.669</td><td>1.272</td><td>1.272</td><td>0.554</td><td>1.305</td><td>1.908</td><td>0.927</td><td>0.422</td><td>1.338</td><td>1.573</td><td>1.190</td><td>2.544</td></tr><tr><td>[18]</td><td>0.376</td><td>0.659</td><td>0.646</td><td>0.588</td><td>0.839</td><td>0.768</td><td>0.813</td><td>0.927</td><td>0.914</td><td>0.785</td><td>0.856</td><td>0.845</td><td>0.874</td></tr><tr><td>[19]</td><td>1.102</td><td>0.679</td><td>0.0</td><td>0.0</td><td>0.75</td><td>0.0</td><td>0.0</td><td>0.960</td><td>0.852</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr></tbody></table>

were set to 0.0 (Tab. 1).

The adaptation to multiple labels is required when there are more than 2 labels in a 2x2x2 configurations. Let exemplify it with the three labels case {$\ell_i, \ell_j, \ell_k$}. To compute surface increments to surfaces $s(\ell_i, \ell_j)$ and $s(\ell_i, \ell_k)$, we first get the local area $ds$ associated to the binary configuration yield by {$\ell_i, \ell_i$} where $\ell_i$ designated the complementary of $\ell_i$ (here $\ell_j$ and $\ell_k$). $ds$ is then proportionally distributed in both $s(\ell_i, \ell_j)$ and $s(\ell_i, \ell_k)$ with coefficients $|\ell_j|/(|\ell_j| + |\ell_k|)$ and $|\ell_k|/(|\ell_j| + |\ell_k|)$, where $|\ell_i|$ denotes the count of label $\ell_i$ in the configuration.

## 2.1. Toy examples

Two toy examples will be used to assess the evaluation of contact surfaces. A cubic image (size of $4L^3$) is built, divided in 8 octants (each has a size of $2L^3$). In the center of the image, a symmetric object of known geometry is drawn, which yield 8 equal contact surfaces for which the ground truth is known. Two objects are considered here (see Fig. 1):

* a sphere of diameter $L + 1$ at the exact center of the image

* a cylinder of height $L$ and diameter $L + 1$.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Toy examples, $L = 20$. Left: sphere; right: cylinder.</figcaption></figure>

## 2.2. Transformations

To assess the robustness of the surface measurement of an object with respect to its orientation, [17, 19] transformed the object (ie apply the transformation to the system of equations defining the object) and then sampled (ie voxelized) it. This method, named *transformation then sampling*, only assess the surface estimation robustness with respect to the object pose. Alternatively, we can consider a first voxelized image of the object (eg Fig. 1) and resample it with a nearest neighbor interpolation. This second method, named *sampling then transformation*, is closer to practical use and also assess the resampling scheme.

The sphere toy has been transformed (with the same transformation) using these two methods in Fig. 2. One can notice that the entire image is filled with the first one, while the second one results in a truncated cube (due to image resampling).

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** The transformed sphere ($L = 20$, rotation of 41 degrees). Right: transformation then sampling; left: sampling then transformation (nearest neighbor interpolation).</figcaption></figure>

## 2.3. Robustness wrt rigid transformation

We consider here random rigid transformations with a rotation direction drawn from the unit sphere, an angle rotation drawn from the interval $[0, 90]$ (in degrees) and a translation in $[-0.5, 0.5]^3$ (in voxels). To assess the robustness, we drawn 1000 random transformations, and grouped the results with respect to the rotation angles in bin of width 5 degrees (Fig. 3).

With the *transformation then sampling* strategy, contact surfaces estimation on the sphere (Fig. 3, upper left) has a constant and small relative error, $\frac{\bar{s}-s}{s}$ where $\bar{s}$ and $s$ are respectively the estimated and the ground truth surfaces, (close to 0% for [17, 18], and close to 0.9% for [19]) whatever the rotation angle, which is expected since the sphere has no preferred orientation. Since sampling depends on the sphere center, displacing it (with the random translation) causes a variation in the surface estimation (depicted by the standard deviations). For the cylinder (Fig. 3, lower left), the three tested methods performed similarly for rotations with an angle larger than 20 degrees. It can be seen that [19] outperformed other methods for smaller angles.

More interesting are the results with the *sampling then transformation* method, where images are resampled with the nearest neighbor interpolation scheme (Fig. 3, middle column). Concerning the sphere, the relative error may reach a maximum about 20% for [17], 10% for [18], and 5% for [19], which represents a huge deterioration of their performances. This deterioration is less severe for the cylinder. The relative error reaches a maximum about 8% for [17], and 5% for both [18] and [19], while the latter has the overall best performance.

We investigate the cause of the performance deterioration. [17] remarked that surfaces are mostly described by 5 configurations (named *basic cases*, depicted in solid line in Tab. 1). For this reason, [19] only optimized the weights associated to these configurations (that they called *informative*) and ignore the others (dashed line in Tab. 1). Fig. 5 demonstrates that the configuration composition drastically changes with the nearest neighbor interpolation with re-
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Relative error in percentage (average ± one standard deviation, bin width = 5 degrees) for contact surface area. Upper row: sphere; lower row: cylinder. From left to right: transformation then sampling; sampling then transformation (nearest neighbor interpolation); sampling then transformation (Gaussian based interpolation, $\sigma = 1.0$). 1000 transformations were drawn.</figcaption></figure>

respect to the *sampling then transformation* case. It is worth noticing that the *other configurations* count may not be negligible. This can also be seen in Fig. 2 (left) where the sphere surface seems to be noisy. One can notice that the counts of these configuration have a similar shape than the relative error of [17]. The large weights associated to them (Tab. 1) may be the main cause for the observed large errors. The same explanation stands, in a less extend, for [18].

### 3. GAUSSIAN BASED INTERPOLATION

Results from Fig. 3 demonstrates that surface estimation is not robust with respect to nearest neighbor interpolation. We propose an other interpolation scheme, which not only is robust to surface estimation, but also allows to smooth surfaces.

**Initialization:** Pick one label $\ell_0$, let $P = SB_0 \circ T_{I \leftarrow I'}$;
Fill $I'$ with $\ell_0$ where $SB_0 \circ T_{I \leftarrow I'}$ is non zero;
Remove $\ell_0$ from $\{\ell_k\}$;
**while** $\{\ell_k\} \neq \emptyset$ **do**
  Pick $\ell_i$ from $\{\ell_k\}$;
  Compute $SB_i = B_i * G_\sigma$;
  **for** Voxels $x$ **do**
    **if** $SB_i \circ T_{I \leftarrow I'}(x) > P(x)$ **then**
      $P(x) \leftarrow SB_0 \circ T_{I \leftarrow I'}(x)$;
      $I'(x) \leftarrow \ell_i$;
    **end**
  **end**
  Remove $\ell_i$ from $\{\ell_k\}$;
**end**

<caption>**Algorithm 1:** Gaussian-based interpolation</caption>

Let $\{\ell_k\}$ be the labels associated to objects in a multi-labeled image $I$ and let $T_{I \leftarrow I'}$ be the transformation that allows to resample $I$ into $I'$ frame, i.e. to compute $I \circ T_{I \leftarrow I'}$. $I$ is represented by $n$ binary images $B_k$, one per label. Each of these binary images is transformed into a kind of probabilistic image $SB_k$, e.g. by smoothing it with a Gaussian (our choice, but any other choice can be made). $SB_k$ is now a grey-valued image that can be interpolated with higher degree scheme than the neighbor neighbor one (we use tri-linear interpolation). In $I'$ the label of a voxel is the one of maximum probability, ie $I'(x) = \arg_i \max SB_i \circ T_{I \leftarrow I'}(x)$. Though concise, the previous equation is not efficient, since it suggests that all $SB_i \circ T_{I \leftarrow I'}$ have to be computed. An iterative update of $I'$

seems to be more adequate (Alg. 1). Choosing the background label as $\ell_0$ (line Init) allows to decrease the computational cost for all labels in our targeted application (embryogenesis, Fig. 4). Indeed, it is sufficient for $B_i$ to be a sub-image of $I$ with sufficient margins to compute $SB_i$ (for Gaussian, this reduces the computational cost of the convolution, line Conv), then the update for $\ell_i$ takes place only in the transformed bounding box of $B_i$ (line Update) and not in the whole image $I'$.

For both the sphere and the cylinder (Fig. 3, right column), we retrieve relative errors comparable to the ones obtained with the *sampling then transformation* method. Fig. 5 (right) demonstrates that the gaussian based interpolation allows also to retrieve the same configuration distribution than the *sampling then transformation* one. For the sphere toy images, the (almost constant) relative error is close to 0.4% for [17, 18], and close to 1.2% for [19]. Interestingly, the gaussian based interpolation allows also to drastically decrease the variation of the surface estimation, as demonstrated by the smallest standard deviations (compare left column of Fig. 3 to the right one).

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Gaussian based interpolation. Upper row: XY-section of a 3D image; lower row: close-up. Left: original data. Middle and right; interpolation with $\sigma = 2.0$ and $\sigma = 10$.</figcaption></figure>

The proposed interpolation scheme has also a smoothing effect on the segmentation as demonstrated by Fig. 4 (where the transformation is set to identity). The drawback is that it tends to shrink the convex parts, thus the choice of $\sigma$ is a compromise between the desired smoothing and the acceptable shrinkage.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Histogram of configuration (Tab. 1) counts in the images used for assessment (average ± one standard deviation, bin width = 5 degrees). Upper row: sphere; lower row: cylinder. From left to right: transformation then sampling; sampling then transformation (nearest neighbor interpolation); sampling then transformation (gaussian based, σ = 1.0). Same transformations as in Fig. 3.</figcaption></figure>

## 4. CONCLUSION

We demonstrated that the surface estimation through local configuration counting is not robust with respect to nearest neighbor interpolation, which is undesirable at least for temporal series where images are resampled for motion compensation. The cause of the poor performances of this estimation is due to the appearance of *non-informative* configurations after interpolation. We introduced then a Gaussian-based interpolation for multi-labeled images that allows not only to preserve the surface estimation after resampling but also to smooth the multi-labeled image.

**Compliance with Ethical Standards:** this is a numerical simulation study for which no ethical approval was required.

**Acknowledgments:** this work was supported by the Cell-whisper (ANR-19-CE13-0020) grant of French National Research Agency.

## 5. REFERENCES

[1] PJ Keller, “Imaging Morphogenesis: Technological Advances and Biological Insights,” *Science*, vol. 340, no. 6137, pp. 1234168+, June 2013.

[2] Y Refahi, A Zardilis, G Michelin, R Wightman, and alii, “A multiscale analysis of early flower development in Arabidopsis provides an integrated view of molecular regulation and growth control,” *Develop. Cell*, vol. 56, no. 4, pp. 540–556, 2021.

[3] L. Guignard et al., “Contact area-dependent cell communication and the morphological invariance of ascidian embryogenesis,” *Science*, vol. 369, no. 6500, pp. 158, July 2020.

[4] B Leggio, J Laussu, A Carlier, C Godin, P Lemaire, and E Faure, “MorphoNet: an interactive online morphological browser to explore complex multi-scale data,” *Nature Communications*, vol. 10, no. 2812, pp. 1–8, 2019.

[5] B Zitová and J Flusser, “Image registration methods: a survey,” *Image Vis Comput*, vol. 21, no. 11, pp. 977–1000, 2003.

[6] C Medioni, F Besse, X Descombes, and G Malandain, “Motion compensation in two-photon microscopy temporal series,” in *ISBI*, Apr. 2015, p. 4.

[7] G Malandain and G Michelin, “Coherent temporal extrapolation of labeled images,” in *International Symposium on Biomedical Imaging (ISBI)*, Apr. 2017, pp. 433–436.

[8] F. Flin, J.-B. Brzoska, D. Coeurjolly, R.A. Pieritz, and alii, “Adaptive estimation of normals and surface area for discrete 3-d objects: application to snow binary data from x-ray tomography,” *IEEE Transactions on Image Processing*, vol. 14, no. 5, pp. 585–596, 2005.

[9] K Itakura and F Hosoi, “Voxel-based leaf area estimation from three-dimensional plant images,” *Journal of Agricultural Meteorology*, vol. 75, no. 4, pp. 211–216, 2019.

[10] J Moukhtar, A Trubuil, K Belcram, D Legland, Z Khadir, A Urbain, JC Palauqui, and P Andrey, “Cell geometry determines symmetric and asymmetric division plane selection in arabidopsis early embryos,” *PLoS Comput Biol*, vol. 15, no. 2, February 2019.

[11] G Malandain and P Lemaire, “Detecting symmetry in ascidian embryo,” in *ISBI 2024 - 21st International Symposium on Biomedical Imaging*, Athenes, Greece, May 2024, p. 4.

[12] YS Liu, J Yi, H Zhang, GQ Zheng, and JC Paul, “Surface area estimation of digitized 3d objects using quasi-monte carlo methods,” *Pattern Recognition*, vol. 43, no. 11, pp. 3900–3909, 2010.

[13] A Lenoir, R Malgouyres, and M Revenu, “Fast computation of the normal vector field of the surface of a 3-d discrete object,” in *DGCI*, 1996, pp. 101–112.

[14] JO Lachaud and A Vialard, “Geometric measures on arbitrary dimensional digital surfaces,” in *DGCI*, 2003, pp. 434–443.

[15] JC Mullikin and PW Verbeek, “Surface area estimation of digitized planes,” *Bioimaging*, vol. 1, pp. 6–16, 1993.

[16] G Windreich, N Kiryati, and G Lohmann, “Voxel-based surface area estimation: from theory to practice,” *Pattern Recognition*, vol. 36, no. 11, pp. 2531–2541, 2003.

[17] J Lindblad, “Surface area estimation of digitized 3d objects using weighted local configurations,” *Image and Vision Computing*, vol. 23, no. 2, pp. 111–122, 2005.

[18] K Schladitz, J Ohser, and W Nagel, “Measuring intrinsic volumes in digital 3d images,” in *DGCI*, 2006, pp. 247–258.

[19] J Ziegel and M Kiderlen, “Estimation of surface area and surface area measure of three-dimensional sets from digitizations,” *Image and Vision Computing*, vol. 28, no. 1, pp. 64–77, 2010.