# Extracting multiple surfaces from 3D microscopy images in complex biological tissues with the Zellige software tool

Céline Trébeau, Jacques Boutet de Monvel, Gizem Altay, Jean-Yves Tinevez, Raphaël Etournay

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

## To cite this version:

Céline Trébeau, Jacques Boutet de Monvel, Gizem Altay, Jean-Yves Tinevez, Raphaël Etournay. Extracting multiple surfaces from 3D microscopy images in complex biological tissues with the Zellige software tool. BMC Biology, 2022, 20 (1), pp.183. 10.1186/s12915-022-01378-0. pasteur-03983119v2

HAL Id: pasteur-03983119

https://pasteur.hal.science/pasteur-03983119v2

Submitted on 10 Feb 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
SOFTWARE

Open Access

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

Check for
updates

# Extracting multiple surfaces from 3D microscopy images in complex biological tissues with the Zellige software tool

Céline Trébeau<fnref n="1" />, Jacques Boutet de Monvel<fnref n="1" />, Gizem Altay<fnref n="1" />, Jean-Yves Tinevez<fnref n="2" />* and Raphaël Etournay<fnref n="1" />*

## Abstract

**Background:** Efficient tools allowing the extraction of 2D surfaces from 3D-microscopy data are essential for studies aiming to decipher the complex cellular choreography through which epithelium morphogenesis takes place during development. Most existing methods allow for the extraction of a single and smooth manifold of sufficiently high signal intensity and contrast, and usually fail when the surface of interest has a rough topography or when its localization is hampered by other surrounding structures of higher contrast. Multiple surface segmentation entails laborious manual annotations of the various surfaces separately.

**Results:** As automating this task is critical in studies involving tissue-tissue or tissue-matrix interaction, we developed the Zellige software, which allows the extraction of a non-prescribed number of surfaces of varying inclination, contrast, and texture from a 3D image. The tool requires the adjustment of a small set of control parameters, for which we provide an intuitive interface implemented as a Fiji plugin.

**Conclusions:** As a proof of principle of the versatility of Zellige, we demonstrate its performance and robustness on synthetic images and on four different types of biological samples, covering a wide range of biological contexts.

**Keywords:** Image analysis, Morphology, 3D imaging, Tissue imaging, Image segmentation, Surface extraction, Fiji plugin

## Background

The interplay between gene regulatory networks and physical forces in driving collective cell behaviors is key to tissue morphogenesis during development and to tissue homeostasis throughout life. Recent quantitative studies of epithelial morphogenesis have begun to unravel the basic cellular and physical principles of tissue development, by providing the tools to integrate multiple scales of tissue dynamics [1–4]. These tools are instrumental to quantify how cell shape changes, cell divisions,

cell rearrangements, and cell extrusions contribute to tissue remodeling, and to establish data-driven computational models of tissue morphogenesis.

Quantitative analysis of an epithelium starts with the extraction of its apical surface from 3D-microscopy images (z-stacks of xy-optical sections) encompassing the volume immediately surrounding the epithelium. To this end, the abovementioned studies relied on specific preparations of the specimen, allowing to expose the entire epithelial surface labeled with junctional fluorescent markers to reveal the network formed by epithelial cell-cell contacts. Once the epithelial surface has been extracted, automated cell segmentation and cell contour tracking tools can be used to follow the dynamics of every cell within the epithelium. However, extracting the epithelial surface remains a difficult task because this

<footnote>*Correspondence: tinevez@pasteur.fr; raphael.etournay@pasteur.fr</footnote>

<footnote n="1">¹ Institut Pasteur, Université Paris Cité, Inserm, Institut de l'Audition, F-75012 Paris, France</footnote>

<footnote n="2">² Institut Pasteur, Université Paris Cité, Image Analysis Hub, F-75015 Paris, France</footnote>
surface is usually not flat (it is best modelled as a curved surface, or 2D submanifold embedded in 3D space), and is often surrounded by other biological structures such as cell layers, acellular membranes, extracellular matrix, and vesicles that hamper its visualization and reconstruction. Even more difficult situations occur when the epithelium surface cannot be readily exposed such as in epithelial organoids in which the epithelial surface faces the lumen inside the cyst. More generally, there is a need for efficient tools allowing one to discriminate a surface of interest from surrounding structures of high intensity and contrast. This is required to study how epithelial cell dynamics account for tissue shape changes in complex specimens, in a wide range of contexts including developmental biology, host-pathogen interaction, and tissue regeneration studies.

Another experimental limitation of previous studies is that some of the structures surrounding the epithelium exert external physical constraints that are known to critically affect epithelial morphogenesis by directing cellular dynamics and signalling pathways [1, 5–7]. To understand the physical forces controlling tissue morphogenesis, it is thus essential to also characterize how the dynamics of these extra-epithelial surfaces relate to that of the epithelium (see [8] for review). This calls for the development of dedicated tools allowing the automated extraction of information from several surfaces of interest in a given sample, since the sheer volume of the data challenges any attempt at a manual analysis [9].

Several surface extraction tools have been developed, some of which are available as open access software, such as PreMosa [10], FastSME [11], and LocalZProjector [12]. These tools focus on the extraction of a single, near-horizontal epithelial layer, which is assumed to (i) be sufficiently smooth, (ii) show enough contrast against surrounding background signals, and (iii) cover the entire image field-of-view. Specifically, it is assumed that the fluorescent marker used to label the epithelial cell network should provide the highest contrast in the image and allow to select it out from autofluorescent extracellular structures such as the cuticle in flies, or other acellular membranes in mammalian epithelia. The surface is then localized using heuristic algorithms based on the detection of the voxels of maximum contrast and/or brightness. However, applying these tools on more complex biological images with several epithelia of weaker contrast often leads to incorrect localization of the surface of interest, and its blending with the nearby unwanted biological structures.

MinCostZ on the other hand is the only available open-source tool that allows the extraction of up to two surfaces from a 3D stack, and imposes explicit continuity constraints on the reconstructed surfaces.

MinCostZ surface extraction relies on a previously developed formulation of the problem as a graph-cut optimization [13]. It is implemented as an ImageJ plugin [14], taking as control parameters, the number of surfaces to extract, the maximum slope, and the range of distances allowed between the surfaces, as well as some user-defined cost function that should reflect the characteristics of the surfaces in term of signal intensity, contrast, and texture. Despite its interest, this approach remains computationally costly and difficult to apply in practice due to the non-trivial choice of the cost function and the need to know beforehand the relative positions of the surfaces to be extracted.

Alternatively, one can segment the surfaces of interest by using supervised machine learning tools such as the software solutions Weka [15] or Ilastik [16], as proposed in the ImSAnE surface reconstruction framework [17]. A deep learning approach, using a network of the U-net type to segment the voxels belonging to a single surface of interest, has also recently been reported [18]. While promising as they can provide state of the art segmentations of epithelial surfaces in difficult imaging conditions, machine learning approaches require the prior manual annotation of a sufficiently large set of surfaces to generate suitable training sets. This process can be very time consuming, often necessitating several rounds of trials and errors to obtain satisfactory results, without guarantees to be generic, i.e., to generalize to a wide range of datasets. So far, no solution to the multiple surface extraction problem has been proposed, which is satisfactory both in terms of genericity and ease of use.

However, such a tool is highly desirable for modern biology studies. Indeed, tissue organization in the context of developmental biology emerges from the interaction of several neighboring structures [8] through the interplay of molecular signals [19, 20], as well as electrical [21, 22], hydraulic [23, 24], and mechanical contact interactions [25, 26]. The importance of such interaction is exemplified in embryonic explanted tissue cultures that develop abnormally when separated from their neighboring structures [27]. Similarly, in the context of tissue engineering, stem cell-derived aggregates harbor various types of tissues surrounding the genuine organoid, and these tissues presumably influence organoid shape, fate, and differentiation (see [28] for review). The ability to simultaneously study the dynamics of neighboring structures together with the structure of interest is therefore essential for an integrated understanding of tissue development, and for any attempt to harness tissue self-organization in vitro.

Here, we introduce Zellige, a tool based on a novel constructive approach that allows the automatic extraction of a non-prescribed number of surfaces from a 3D
image. Unlike most other tools, Zellige does not rely on
the localization of voxels of highest contrast or bright-
ness. Rather, it performs the surface reconstruction in
two steps: it first classifies voxels most likely belonging
to surface structures by a combination of maxima detec-
tion along the z-axis and the application of two bright-
ness and contrast classifiers; it then extracts each surface
present in the image volume with an assembly algorithm
implementing explicit continuity constraints. This two-
step approach allows the reconstruction of multiple sur-
faces of different contrast characteristics to be selectively
reconstructed.

To do this, the user is only required to adjust a small set of intuitive control parameters, a task largely facilitated by a user-friendly interface implemented as a plugin for the open-source Fiji platform [29]. We tested the performance and robustness of Zellige for multiple surface extraction by applying it to synthetic and 3D microscopy images containing multiple surfaces of interest of widely varying texture and contrast. Among the four biological samples we considered, three previously required laborious manual extraction of the epithelial surface (fly pupa histoblast, mouse cochlea, and inner ear organoid) therefore greatly challenging the study of the cellular mechanisms involved in the development of these large organs [9]. We also analyzed a bronchial epithelium specimen that highlights the use of Zellige in a recent study to address how SARS-CoV-2 infection interferes with the epithelial function [30]. These experiments demonstrate the ability of the approach to extract several surfaces of potentially very low contrast, selectively from other highly contrasted and complex structures, with a single set of reconstruction parameters. A sensitivity analysis also reveals a high robustness of Zellige against small variations of these parameters. This will make it a tool of choice in terms of versatility and ease of use for the investigation of biological surfaces.

Implementation

Zellige was devised with the goal of achieving accurate
segmentation of multiple biological surfaces from 3D
confocal images. Unlike other existing surface extraction
tools, it makes no assumption on the number of surfaces
to be extracted and does not require the surfaces of inter-
est to be the structures of highest contrast in the image.
Zellige is written in Java, relying on the ImgLib2 library
[31] and is distributed as a Fiji plugin, with a graphical
user interface (GUI) designed to allow users to quickly
find a good set of extraction parameters for a given
image.

Zellig extracts each surface present in the image in
the form of a height map (or z-map), that is, a mapping
z=h(x,y), which associates to each point (x,y) over which

the surface projects, the z-coordinate of the unique voxel
(x,y,z) belonging to the surface. Each extracted height
map is then used to produce a 2D projection of the 3D
stack restricted to a small subvolume (of user-selected
width) centered around the corresponding surface. To
achieve a robust extraction, Zellige proceeds in two algo-
rithmic steps, which are only outlined below (see Addi-
tional file 1: Supplementary Note 1 for implementation
details) [32].

In the first step, or *surface voxel selection* step, a segmentation is applied to the 3D image to select voxels that likely belong to a surface of interest (Fig. 1 and Additional file 2: Figure S1A). These putative surface voxels are detected as local maxima of image intensity along the z-axis, after applying a 3D Gaussian filter to the image to reduce noise artifacts, followed by two independent binary classifiers, one based on voxel contrast and the other one on voxel intensity (see Additional file 1: Supplementary Note 1 for details). Five adjustable parameters control the selection step: two threshold parameters (*T*<sub>A</sub> and *T*<sub>otsu</sub>) control the strength of the binary classifiers applied on contrast and intensity, respectively, and three parameters (*S*<sub>min</sub>, *σ*<sub>xy</sub>, and *σ*<sub>z</sub>) control clean-up operations applied at the end of the classification (removal of small isolated spots, and local averaging along the *xy* plane and the *z*-axis, respectively).

In the second step, or *surface assembly step*, an iterative
algorithm is used to extract the height maps of each of the
surfaces present in the image (Fig. 1 and Additional file 2:
Figure S1A). The assembly starts by grouping together
contiguous chains of putative surface voxels within each
orthogonal *xz* section (alternatively, within each *yz* section)
of the 3D image, in order to form a set of building blocks
referred to as *orthogonal surface elements* (OSEs). These
building blocks are then used to assemble the surfaces, in
a process analogous to jigsaw puzzles, where OSEs adja-
cent to the surface boundary are added if they match this
boundary, and rejected otherwise, according to matching
criteria that impose explicit continuity constraints to the
surface under construction. (See Additional file 1: Supple-
mentary Note 1 for details on these matching criteria and
on the construction of the OSEs.) The surface assembly
proceeds until no matching OSE can be found (Additional
file 2: Figure S1). In order to increase the robustness of the
surface assembly step, Zellige applies it in two rounds, con-
structing the OSEs within *yz* sections only and performing
the assembly along the *x*-axis during the first round, and
constructing OSEs within the *xz* sections and proceeding
along the *y*-axis during the second round (see Additional
file 1: Supplementary Note 1). Each round is controlled by
3 adjustable parameters: a threshold parameter (0 ≤ *T*<sub>OSE</sub>
≤1) sets a minimum size for the building blocks that can
be used as seeds to initiate the assembly of a surface; and
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Flowchart of Zellige's algorithmic steps. Surface voxel selection (step 1), surface assembly in the form of a height map (step 2), and subsequent projection localized to the height map, are schematically depicted in the case of a 3D image containing 4 surfaces of interest</figcaption></figure>

two other parameters ($0 \le R \le 50$ and $0 \le C \le 1$) set the matching constraints used to accept or reject the addition of OSEs to a surface. The assembly step is thus controlled overall by 6 parameters, i.e., two groups of 3 parameters ($T_{OSE1}$, $R_1$, $C_1$) and ($T_{OSE2}$, $R_2$, $C_2$) controlling the first and second assembly rounds, respectively.

Finally, the height maps of each of the reconstructed surfaces are used to obtain a corresponding 2D projection (Fig. 1). In practice, a maximum projection restricted to a subvolume of width $\delta z$ ($\delta z$ being a user-defined parameter) centered around the surface of interest is performed.

# Results and discussion

## Proof of concept of multiple surface extraction on a synthetic image

Figure 2 shows the results produced by Zellige on a phantom 3D image [33] containing three distinct synthetic

surfaces generated as described in Additional file 1: Supplementary Note 2 [34]. The three surfaces are extracted with little errors. We assessed the quality of the reconstruction by comparing each of the height maps produced by Zellige to the corresponding ground truth (GT) height map, which is exactly known in this case (Fig. 2A–C, and Additional file 1: Supplementary Note 3) [12, 35]. For the three surfaces, the absolute value $|h(x,y) - h_{GT}(x,y)|$ of the difference between the reconstructed and the GT height maps is $\le 1$ (in pixel units) over >99% of the GT pixels (Fig. 2D–E), with a root mean square error value (RMSE; defined over the pixels common to the reconstructed and the GT height maps, see Eq. 1 in Additional file 1: Supplementary Note 3) of $\le 0.6$ in pixel units, showing that the surface localization is highly accurate. In addition, the coverage, which measures the proportion of the reconstructed surface relative to the GT, is near 100% for the three surfaces. To achieve these results, the
control parameters of the two steps of the surface extrac-
tion were adjusted manually to some adequate reference
values (see Additional file 1: Supplementary Table S1)
using the Zellige Fiji interface. Only the parameters con-
trolling the voxel classification step (amplitude and Otsu
threshold parameters T<sub>A</sub> and T<sub>otsu</sub>, minimal island size
S<sub>min</sub>, and smoothing parameters σ<sub>xy</sub> and σ<sub>z</sub>) did actually
require a modest adjustment. The parameters of the sur-
face assembly step (parameters T<sub>OSE1</sub>, R<sub>1</sub>, C<sub>1</sub> and T<sub>OSE2</sub>,
R<sub>2</sub>, C<sub>2</sub> of the 1st and 2nd construction rounds, respec-
tively) were set to their default reference values (see
Additional file 3: Figure S2 and Additional file 1: Supple-
mentary Note 4) and did not need to be adjusted.

Thus, using a single set of control parameters, Zellige
can extract multiple surfaces of various shapes and tex-
tures with very little error, without requiring the user to
provide information about their number or relative posi-
tion, nor about their shape or texture characteristics.

Performance of Zellige on biological samples

Example 1: extracting multiple surfaces from an image
of a pupal fly specimen

Over the past few decades, the *Drosophila* model has
been invaluable to decipher the molecular and cellular
mechanisms underlying organ embryogenesis [36, 37].
Epithelium morphogenesis studies not only revealed
the importance of mechanical stresses (including stress
boundary conditions) and planar polarity signalling on
cell dynamics to generate tissues of reproducible sizes
and shapes, it also highlighted the importance of extra-
cellular matrix attachments in constraining the tis-
sue stresses that guide patterning [1, 38]. At the pupal
stage, the fly undergoes dramatic remodeling of its larval
organs into adult organs. Large-scale tissue flows initiate
at a timing that coincides with molting, when the epithe-
lium contracts away from the overlaying cuticular sac, a
protective acellular membrane that imposes mechanical
constraints to the tissue.

Figure 3 shows the results of applying Zellige on a 3D
image of a *Drosophila* pupa acquired with a spinning
disk confocal microscope [39]. The sample expresses
Ecadherin-GFP, a fluorescent marker of cell-cell junc-
tions, and encompasses a portion of the pupa’s abdomen
and a small portion of its wing. Four surfaces of interest

can be identified, with varying signal intensities, noise
levels, and features (Fig. 3A,B). The abdomen is formed
of an epithelium (surface S2) overlaid by a cuticle (surface S1). Lying just beneath these two surfaces, one can
observe globular structures showing in some places a
higher intensity than the signal coming from the surfaces.
The wing also consists of an epithelium of low-intensity
signal (surface S4), and an overlying cuticle (surface S3).
These two surfaces are relatively flat, except for surface
S3 which is very steep near one of its edges.

Figure 3C shows a 3D graphical representation of the height maps reconstructed by Zellige (green) and those reconstructed by an expert biologist (blue), taken as ground truth (GT). While these height maps clearly show greater roughness than those of the synthetic surfaces presented earlier, they could again be obtained with a single set of control parameters that were adjusted manually with the Zellige interface (see Additional file 1: Supplementary Table S1). We observe an excellent match between the four reconstructed and corresponding GT height maps, despite the rather complex topography of surfaces S1 and S2 (with slopes reaching up to 45°), and the near-vertical inclination of surface S3 at its boundary. Yet, small deviations may be seen in the regions of highest slope of the surfaces. Some of these deviations are likely attributable to uncertainties in the definition of the GT height maps, whose accuracy depends on the expert.

Figure 3D shows the differences between the recon-
structed and GT height maps, plotted as color-coded
error maps. These differences are <2 (in pixel units)
for most pixels, while some regions of higher error
values can be seen locally in surfaces S1 and S2, and
at the boundary of surface S3. Note that for surfaces
S2 and S4, which contain junctional epithelial meshes
composed of larger and smaller cells, respectively, the
GT height map encompass not only the mesh but also
the interior of the cells, where no junctional signal is
detected. The distance calculated inside the cells is thus
more subjected to intensity fluctuations, especially for
surface S2. Nonetheless, the RMSEs of surfaces S1,
S2, and S4 are less than 1, showing that on average the
reconstructed height maps match the corresponding
ground truth with subpixel accuracy (i.e., RMSE < 1).
The higher RMSE (1.25) of surface S3 is largely due to

(See figure on next page.)

**Fig. 2** Multiple surface extraction on a synthetic 3D image. **A,B** The image contains 3 phantom surfaces (S1, S2, S3) of different shapes (sinusoidal, flat, and paraboloidal, respectively), and different textures (surface S1 has constant intensity, while surfaces S2 and S3 are supported by Voronoi meshes of different cell sizes). **C** 3D representations of the height maps extracted by Zellige (in green) and of the ground truth (GT, in blue) height maps of surfaces S1, S2, and S3. **D** Error maps displaying the distance along the z-axis between the reconstructed and GT height maps for surfaces S1, S2, and S3. **E** Projections of the 3D image localized to the different surfaces S1–S3 (maximum intensity projections over a subvolume of a width δz=1 pixel above or below the corresponding height maps). Upper and lower panels show the projections based on the GT and the reconstructed height maps, respectively
<figure><img src="image_5.png" /><figcaption>**Fig. 2** (See legend on previous page.)</figcaption></figure>
the region of steep slope at the edge of this surface (yellow region on the error map for this surface, Fig. 3D). The coverage of the reconstructed height maps is excellent (≥ 96%) for surfaces S1 and S2, and slightly lower, but still very good (≥ 93%) for the smaller surfaces S3 and S4. Figure 3E shows the 2D projections of the 3D image obtained for each of the reconstructed surfaces and for the corresponding ground truths. The inaccuracies visible on the error maps (see Fig. 3D) do not significantly impact these projections, which appear very similar to the projections obtained with the corresponding ground truths. Thus, while the biological sample contains significant noise and shows a much more variable contrast (especially with the presence of high-intensity globular structures underneath surfaces S1 and S2), Zellige makes it possible to segment these surfaces selectively, with a quality of segmentation comparable to that obtained by manual expert segmentation.

This possibility brings several perspectives that are not offered by single surface extraction algorithms. First, it opens the possibility to systematically study the tissue axial movements (along *z*) relative to the cuticle during molting, allowing for example to gain insights into the early tissue contraction of the wing hinge that acts as a mechanical inducer over the wing blade [8]. Second, Zellige makes it possible to automatically extract structures such as the abdomen epithelium, which is usually segmented manually [9], due to the difficulty to separate the large larval cells from the cuticle mesh and from other globular structures (such as fat bodies or macrophages) present underneath the epithelium. All these structures become intertwined when using other extracting tools. In this context, Zellige opens new opportunities to study collective cell behavior during epithelial morphogenesis in vivo, and to integrate in the analysis the surrounding surface-like structures involved in the mechanics of the system.

**Example 2: extracting a thin cochlear epithelium surface from a multilayer dataset**

As the first model in which planar cell polarity signaling was shown to be conserved in vertebrates [40], the mammalian auditory organ, the cochlea, is arguably our most valuable model to study epithelial patterning and morphogenesis beyond the fly and zebrafish [41, 42]. Cochlear morphogenesis involves complex and tightly controlled patterning processes during which the cochlear sensory epithelium extends and develops its characteristic coiled snail shape, while adopting a striking cellular mosaic organization, with graded changes of morphogenetic parameters along the cochlea [43, 44]. These morphogenetic processes are well recapitulated in organotypic cultures, on the condition that the mesenchyme that underlies the epithelium be preserved. The cultures are then amenable to live imaging [42, 45], pharmacological [46], and genetic manipulations.

Figure 4A shows a confocal swept-field microscope acquisition of an embryonic mouse cochlea [47]. The sample contains only one surface of interest, the cochlear epithelium, but this surface lays on top of a thick tangled mesh of non-epithelial cells originating from the mesenchyme. The whole biological tissue is stained for filamentous actin (F-actin) using phalloidin. The epithelium surface presents a non-uniform signal included in a small *z*-range (6 ≤ *z* ≤ 10), and a mesh of very heterogeneous size. Between sections *z* = 10 and *z* = 14, one can observe the basolateral region of the epithelial cells, also stained for F-actin. The particularity of this sample is that the mesenchyme presents an intense and contrasted signal over a wide range of *z*-values (14 ≤ *z* ≤ 43). This makes it challenging to extract the surface of the epithelium, which is characterized by low intensity and low contrast.

Figure 4B shows a 3D representation of the height map reconstructed by Zellige and the corresponding GT height map (again reconstructed manually by an expert).

(See figure on next page.)

**Fig. 3** Fly specimen. A,B Volume rendering (A) and orthogonal sections (B) of a 3D image of fly embryo taken around 24 h after puparium formation, covering a portion of the abdomen (showing histoblast cells and larval cells), and a portion of the developing wing. Scale bar 50 µm. Four surfaces of interest may be identified in the dataset (of dimensions 1200 × 1200 × 51 voxels): surfaces S1 and S2 are relatively close to one another and located within overlapping *z*-ranges (8 ≤ *z* ≤ 50 and 20 ≤ *z* ≤ 50, respectively). Surfaces S3 and S4 (located in the *z*-ranges 42 ≤ *z* ≤ 50 and 9 ≤ *z* ≤ 50, respectively) are relatively far from each other and can nearly be separated by a plane. C 3D representations of the height maps extracted by Zellige (in green) and of the ground truth height maps (GT, in blue) of surfaces S1–S4. The reconstructed height maps of all surfaces S1–S4 cover >93% of the area of the corresponding GT (cf. Additional file 3: Figure S2 and Additional file 1: Supplementary Table S1). To reduce the staircase artifacts (more or less visible depending on the surface) due to the digitization of the GT and reconstructed height maps, all height maps were smoothed with a 2D Gaussian filter with a standard radius of 5 pixels (cf. Additional file 1: Supplementary Note 1). D Error maps (color-coded distance along the *z*-axis between the reconstructed and the GT height maps) plotted for each of the reconstructed surfaces. The large majority of pixels on the reconstructed height maps (98%, 96%, 91%, and 99% for surfaces S1 to S4, respectively) display errors of <2 pixels. The height maps of surfaces S1, S2, S4 show subpixel accuracy on average (RMSE < 1), while that of surface S3 is slightly less accurate (RMSE = 1.25). E Projections of the 3D image localized to the different surfaces S1–S4 (in this and all subsequent figures, these are maximum intensity projections over a subvolume of width δz=±1 pixel above or below the corresponding height maps). Upper and lower panels show the projections based on the GT and the reconstructed height maps, respectively. Scale bar 50 µm
<figure><img src="image_6.png" /><figcaption>**Fig. 3** (See legend on previous page.)</figcaption></figure>
<figure><img src="image_7.png" /></figure>

On the corresponding error map (Fig. 4C), most (~83%) pixels of the reconstructed height map show subpixel accuracy (with distances < 1 to the corresponding pixels of the GT height map). The errors are greater in regions where the cell size is larger, as well as in the area where the signal intensity is very low. However, they remain smaller than 2 for > 95% of the pixels. This result is consistent with the low value of the RMSE (1.1). The surface is also reconstructed with an excellent coverage (> 99%, Additional file 1: Supplementary Table S1).

As the sample contains a single epithelial surface of interest, we compared the performance of Zellige with three other software that can extract only a single surface (Fig. 4D). The projections of PreMosa, FastSME, and LocalZProjector completely miss the epithelium, despite our attempts to find an optimal set of parameters using a multi-dimensional parameter screening (as done previously in Herbert et al. [12]). Only regions of high contrast corresponding to the mesenchyme are projected. In contrast, Zellige generates a projection very close to the

ground truth. This demonstrates the efficiency of Zellige to selectively extract a low contrast surface, despite the presence of several structures of higher contrast. Indeed, Zellige detects every structure as a possible surface seed without any assumption on its contrast, and only extends this seed into a surface if enough spatial continuity is found in the surrounding signal. This feature allows to separate individual surfaces from other structures spatially, which should greatly facilitate the analysis of live imaging experiments.

**Example 3: extracting a single bronchial epithelial surface rendered abnormally rough by SARS-CoV-2**

Recently, we used Zellige to extract the surface of a primary culture of bronchial epithelial cells following infection by the SARS-CoV-2 virus [30]. The infection causes the surface of the epithelium to become abnormally rough due to cell damages as seen from discontinuities within the cell layer. The sample we chose from this study is a 3D confocal image of the epithelium
<figure><img src="image_8.png" /></figure>

responding to SARS-CoV-2 infection [48] (Fig. 5A). The surface of interest in this image corresponds to the layer of epithelial cells stained for the tight junction protein Zona Occludens-1 (ZO-1). The surface roughness causes the network of junctions to extend over the height of the z-stack, with a signal of varying intensity (Fig. 5A). In addition, the junctional network remains non-plar-nar even at the level of a single cell, hence violating the smoothness condition commonly assumed to hold in the context of epithelial surface extraction. We also observe the presence of nearby punctiform structures of high contrast that are mainly located outside of the epithelium surface. This sample therefore provides an example of a surface with a complex landscape, interspersed with a constellation of signals which may interfere with the sur-face segmentation.

The 3D representation of the reconstructed and corre-
sponding GT height maps (Fig. 5B) makes it possible to
appreciate the roughness of the surface of interest. The
two height maps overlap quite satisfactorily. As shown on
the error map (Fig. 5C), a large majority (71.1%) of pix-
els of the reconstructed height map show errors smaller
than 1 pixel (~ 96% of them showing errors smaller than

2 pixels). The error is however larger in regions where the
cell size is larger, as well as in areas where the ZO-1 signal
intensity is very low, preventing a complete reconstruc-
tion of the junctions. Nevertheless, the overall RMSE
remains small, with a value of 0.81 (Fig. 5). The coverage
of the reconstructed surface is also excellent (98% of the
GT height map), despite the abovementioned discontinu-
ities. Figure 5D shows the comparison of the projection
generated by Zellige by manually adjusting parameters to
those produced by PreMosa, FastSME, and LocalZPro-
jector after a search of an optimal set of parameters (see
Herbert et al. [12]). The projection generated by Pre-
Mosa misses many junctions of the epithelial network,
but it removes quite well the punctiform signal originat-
ing from other optical sections. FastSME performs better
than PreMosa in reconstructing the junctions, but it pro-
duces a projection where the punctiform signal remains
strong. In contrast, Zellige and LocalZProjector manage
to both reconstruct the surface well and to filter out the
punctiform signal quite effectively. This result demon-
strates the efficiency of Zellige to extract a surface with
complex topography by excluding intense and contrasted
spurious signals away from the epithelium surface.
**Example 4: extracting the apical and basal layers of a dome-shaped epithelium (developing inner ear organoid)**

Organoids are stem cell-derived and self-organizing 3D tissue structures that can mimic certain organ structures. They have emerged as promising in vitro models for developmental biology research, as well as biomedical translational research applications. Here we take the example of mouse stem cell-derived inner ear organoids that form vesicular structures composed of an epithelium harboring sensory cells. These organoids are part of a cellular aggregate that also contains other tissues such as the mesenchyme [49] adjacent to the organoids. The epithelial cells of the forming inner ear organoids acquire a basal-apical polarity, with their apical side facing the lumen of the organoid, and their basal side facing outwards. The apical junctional network of the epithelium is difficult to visualize in microscopy images as it is seen from below, through the basal layer. Another difficulty is the spherical geometry of the vesicle system, which makes the epithelial surface of interest difficult to extract in regions of high inclination relative to the focal plane.

Figure 6 shows the result of applying Zellige on a 3D confocal microscopy image of half of a developing inner ear organoid at 14 days of culture, a stage at which markers characteristic of the mouse otic vesicle can be detected [50]. The sample was fixed and stained for F-actin to visualize all cellular structures including the epithelium. Two surfaces of interest can be identified, namely the basal side of the epithelium and the apical junctional network (Fig. 6A,B). Both surfaces are mesh-like structures of high inclination, high signal intensity, and high contrast. The vesicle lumen also contains cell debris of high intensity and contrast that are not part of any surface of interest.

Figure 6C shows a 3D graphical representation of the height maps reconstructed by Zellige and those reconstructed by an expert biologist, taken as ground truth (GT). Due to their dome-shaped topography, the manual

segmentation of these surfaces was rather laborious and is more likely prone to errors in the regions of high inclination. Figure 6D shows that the distance between the two height maps is <2 for the large majority (> 92%) of pixels, while larger error values occur locally in regions of near-vertical inclination of the surfaces. The RMSEs of both apical and basal surfaces are close to 1, showing that on average the reconstructed height maps match the corresponding GT with about 1 pixel of accuracy. The coverage of the reconstructed height map is nearly 100% for the basal surface, and ~88% for the apical surface. Figure 6E shows the 2D projections of the 3D image obtained for each of the reconstructed surfaces, as compared to the projected GT height maps. Thus, for the extraction for these highly inclined surfaces, Zellige produces height maps of a quality comparable to those obtained from manual expert segmentation.

In this example, Zellige could be combined with a 2D cell tracking framework such as TissueMiner [3] to perform cell shape and topology analysis. Note that, using complementary tools such as DeProj [12], geometric distortions introduced in the projected surfaces by the epithelium inclination could be corrected for (Fig. 6F). This approach could provide a means to quantitatively address how an inner ear organoid epithelium patterns at the cellular and organoid scales, while quantifying the epithelial thickness changes due to cellular intercalation or cell shape changes in the depth of the epithelium. This would also permit to better characterize the variability of inner ear organoids within a given aggregate, and it could allow one to explore how the organoid interacts with surrounding tissues and how these interactions influence the differentiation of their constituent sensory cells.

**A sensitivity analysis reveals the robustness of Zellige in extracting surfaces from biological images**

To evaluate further the quality and robustness of the segmentation obtained by Zellige, we carried out a sensitivity

(See figure on next page.)

**Fig. 6** Organoid specimen. A,B Volume rendering (A) and orthogonal sections (B) of a confocal 3D image of a (half of) inner ear organoid, which has been fixed and stained with phalloidin to reveal F-actin. The dataset (of dimensions 520 × 465 × 35 voxels) includes two dome-shaped epithelial surfaces of interest, forming the apical (inward) and basal (outward) sides of the organoid. C 3D representations of the height map extracted by Zellige (in green) and the GT height map (in blue), of the epithelium surface. D Color-coded error maps of the reconstructed height maps for the apical (left) and basal (right) epithelial surfaces of the organoid. The surfaces of interest are reconstructed with an error of < 2 pixels over a large majority (96% and 93% for the apical and basal surfaces, respectively) of pixels, as well as on average (RMSE ~ 0.8 and 1.1 for the apical and the basal surfaces, respectively). E Projections localized to the GT height maps of the epithelium surface (panels on the left), and the height maps extracted by Zellige (panels on the right). F Using the DeProj tool (Herbert et al. [12]) to generate cell morphology measurements corrected for the projection factor associated with the local surface slope. Cells of the surface S2, which represents the epithelial apical surface, were segmented using the TissueAnalyzer tool (Etournay et al. [3]). The height map obtained from Zellige and the segmentation mask were used to calculate geometrical parameters of the surface and its constituent cells. From left to right: the local slope, the local mean curvature, the projection error on cell area due to the local slope, and the corrected cell area in a 3D representation (axis ratio 1:1:3), for the reconstructed surface S2. Scale bar 100 µm
A

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

B

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

C

GT

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

Zellige

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

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

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

D

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

E

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

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

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

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

F

Slope of the plane (°)

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

Local mean curvature (10⁻³ µm⁻¹)

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

Error on cell area (%)

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

Cell apical area (µm²)

<figure><img src="image_23.png" /><figcaption>**Fig. 6** (See legend on previous page.)</figcaption></figure>

analysis of the reconstruction on each of the samples tested. This analysis consisted in varying one control parameter at a time (Additional files 3, 4, 5, 6 and 7: Figure S2-S6, Fig. 7, Additional file 1: Supplementary Note 4), while keeping the other parameters fixed at a nominal value (Additional file 1: Supplementary Table S1). The RMSE and coverage of each of the reconstructed surfaces

were evaluated and plotted as a function of the value of the parameter that was varied.

Additional file 4: Figure S3 shows the results of the sensitivity analysis carried out on the image of Fig. 3 (Example 1, pupal fly specimen), when varying the parameters controlling step 1, i.e., the selection of putative surface voxels (parameters $T_A$, $T_{otsu}$, $S_{min}$, $\sigma_{xy}$, and $\sigma_z$). As can be
<figure><img src="image_24.png" /><figcaption>**Fig. 7** Summary of the sensitivity analysis. The intervals indicated in gray for each parameter and each of the images tested correspond to the parameter values for which the reconstruction satisfies high-quality criteria defined by RMSE ≤ 1.5 and coverage ≥ 85%. Black marks indicate the reference value obtained by manual adjustment for each image (cf. Additional file 1: Supplementary Note 2). **A** Parameters of the surface selection step. **B** Parameters of the surface assembly step</figcaption></figure>

seen, the variations of the two classification threshold parameters $T_A$, $T_{otsu}$, and of the minimum island size $S_{min}$ in their respective intervals do not substantially modify the RMSE and the coverage of the reconstructed surfaces, whose values remain roughly constant for $T_A \le 8$, $T_{otsu} \le 12$, and over the entire $S_{min}$ interval (Additional file 4: Figure S3A). Within these intervals, the RMSEs of all the reconstructed surfaces remain $\le 1.5$, while the coverage values are >95% for surfaces S1 and S2, and >85% for surfaces S3 and S4. The surfaces S1 and S3 show lower signal intensity and lower contrast than surfaces S2 and S3, making them more difficult to extract. Surface S4 has the lowest contrast and fails to be reconstructed if the classification threshold values are too stringent (namely for $T_{otsu} > 12$, or $T_A > 8$). Nevertheless, the intervals of stability of $T_A$ and $T_{otsu}$ (that is, the intervals of values over which a high-quality extraction of all surfaces is obtained) remain relatively wide (cf. Fig. 7). The smoothing parameters $\sigma_{xy}$ and $\sigma_z$ also have some effect on the reconstruction of the surfaces. When $\sigma_{xy}$ is less than 3, the RMSE is higher for the mesh-like epithelial surfaces S2 and S4 (formed by the junctional network of the epithelium). A minimal smoothing along the axial direction is also important to ensure that the reconstructed surfaces are not too fragmented, preventing their complete reconstruction. Yet, $\sigma_z$ should not be chosen too large either,

to avoid merging nearby surfaces along the z-axis. In this case, the closely positioned surfaces S1 and S2 are well-separated if setting $\sigma_z$ close to 1, but they become merged when $\sigma_z > 2$. In general, the surface construction parameters have little effect on RMSE and coverage (Additional file 4: Figure S3B). The sensitivity analysis on this challenging specimen shows a good robustness of Zellige to extract the four surfaces with a single set of parameters, each of which can be chosen in a reasonably wide interval considering the other fixed.

The results of the sensitivity analyses performed with the other biological image stacks (examples 2, 3, and 4 described above) are shown in Additional files 5-7: Figures S4-S6 (Additional file 1: Supplementary Note 4). These results are summarized in Fig. 7, which show the stability intervals over which the extracted height maps satisfy the criteria RMSE $\le 1.5$ and coverage $\ge 85\%$, for which the reconstruction can be considered of high quality. Overall, the stability intervals for the two classification threshold parameters ($T_A$ and $T_{otsu}$) are narrower for specimens containing a surface of low signal intensity and low contrast relatively to other structure of higher signal, as in examples 2 and 3. The graphical user interface of Zellige allows the user to adjust $T_A$ and $T_{otsu}$ interactively, making it intuitive to search for reasonable values. We found that $2 \le \sigma_{xy} \le 3$ and $\sigma_z = 1$
generally give high-quality results for all tested specimens (Fig. 7A). We therefore expect only little adjustment to be required by the user on the smoothing parameters from their default permissive values (set to $σ_{xy} = 1$ and $σ_z = 1$). Values of $σ_z$ that are too large may lead to the merging of surfaces with a nearby structure of high contrast (surface or else), as it happens for the epithelium surface in the cochlea specimen, which merges with the underlying mesenchyme signal when $σ_z$ is greater than 2 (Fig. 4 and Additional file 5: Figure S4). The effect is even more pronounced for closely positioned surfaces (Figs. 2 and 6 and Additional files 3, 7: Figures S2, S6).

Regarding the control parameters of step 2 (the assembly step), for the four examples presented except for example 1, the quality of the reconstructed surfaces is stable and high over the main part of their respective intervals of variation, deteriorating occasionally only when extreme values are used for these parameters (see Additional files 3, 4, 5, 6 and 7: Figures S2-S6, Fig. 7B and Additional file 1: Supplementary Note 4). Example 1 poses particularly stringent constraints on the control parameters of the reconstruction due to the requirement of reconstructing the 4 surfaces of different contrast and texture present in this sample using the same set of parameters.

Finally, the computation times for running Zellige on a given dataset ranged between a few seconds and a minute on a standard PC computer (see Additional file 8: Figure S7 and Additional file 1: Supplementary Note 5), except in a few exceptional cases corresponding to extreme values of the control parameters (not shown). As a safeguard, a stopping criterion could be implemented so as to exit the run (declaring the current parameter values invalid) if the surface assembly computation exceeds a user-prescribed duration.

Overall, the sensitivity analysis indicates that the surface extraction performed by Zellige is robust to variations of the control parameters of step 1 (surface voxel selection step). In general, the reconstruction is more sensitive to the amplitude threshold parameter ($T_A$), and the Otsu threshold parameter ($T_{otsu}$) should be kept sufficiently low for samples containing surfaces of intensity close to the background. However, in some cases such as in our example 3, the opposite is true. Thus, the two threshold parameters play somewhat complementary roles, and the possibility to adjust them independently is useful in practice to be able to cover as many cases as possible. A smoothing along *xy* appears necessary to correctly reconstruct the surfaces supported by a junctional mesh. Not surprisingly, best results are obtained when the radius of the Gaussian filter used for this (parameter $σ_{xy}$) is adapted to the mesh (or cell) size. Likewise,

a smoothing along *z* is beneficial, but the extent of this smoothing (parameter $σ_z$) should not be too large to avoid causing the fusion of nearby surfaces. With a few exceptions, the values of the RMSE and coverage show little sensitivity to the values of the parameters controlling step 2 (surface assembly step), at least once putative surface voxels have been properly selected. In the presence of several surfaces of potentially very different sizes, the parameter controlling the fraction of OSE sizes allowed for OSE seeds (parameter $T_{OSE1}$) should be relatively large (≥0.5 or greater, i.e., allowing more than 50% of the largest OSE sizes for seeds) to allow the extraction of a surface of small size (for example, to extract the surface S4 of example 1, which covers less than 20% of the *xy*-field of view, $T_{OSE1}$ must be larger than 0.6). Extreme values (close to 1) for the connectivity rates (C1 and C2) are too stringent and lead to a drop in the coverage of the reconstructed surfaces. To sum up, we see that the most critical parameters for a satisfactory extraction of the different surfaces are those controlling step 1. In most cases, the parameters controlling step 2 do not need to be adjusted and can be fixed to their default reference values.

## Conclusion

We have developed Zellige, a new tool to extract multiple surfaces from 3D fluorescence microscopy images. Zellige automatically finds surfaces by first identifying putative voxels that are likely to belong to a biological surface, and second by assembling a surface through connection of adjacent voxels satisfying natural proximity constraints. By using Zellige on synthetic epithelium images, we have shown that it accurately reconstructs a surface with excellent performances in terms of both the distance to the ground truth height map and the surface coverage (Fig. 2). Zellige can deal with complex images containing multiple surfaces, with computation times not exceeding a few tens of seconds on a standard computer. Importantly, the user is not required to specify the number of surfaces to be extracted. In the *Drosophila* specimen (Fig. 3), the software readily extracts the 4 surfaces of interest that could be identified. Since Zellige detects putative surface voxels in the first step by combining local and global thresholds, it can discriminate between a surface of low intensity and some nearby surface or other structures of higher contrast, as is the case in the mouse cochlear embryo (Fig. 4). With this difficult dataset, we could also confirm Zellige's robustness against very low signal-to-noise levels. The constructive approach of surface region growing used by Zellige in its second step enables it to circumvent the surface smoothness requirement that is classically assumed by other surface
extraction tools. For instance, it could reconstruct the highly irregular surface of a bronchial tissue infected by SARS-CoV-2 (Fig. 5).

The robustness and flexibility of Zellige come at a price, namely, the requirement to specify 12 parameters when running the surface extraction. However, the sensitivity analysis we performed shows that adjusting only 4 of these parameters is enough in practice to handle a wide range of image types. These parameters correspond to intuitive notions (e.g., thresholding and smoothing), which makes Zellige particularly easy to use. The Fiji interface that we implemented to perform this adjustment should make Zellige even more user-friendly and effective for biological applications.

To our knowledge, Zellige is the only open-source tool that can extract an unspecified number of epithelial surfaces from a 3D volume, possibly larger than two. This unique feature is especially useful in complex images that could be processed only by specialized tools before. For instance, Zellige can extract surfaces with projections on the *xy* plane that completely overlap, such as the basal and apical epithelia in the organoid image of Fig. 6. Previously, such an image could be processed only by tools that relied on segmenting a mesh around the object surface, such as MorphoGraphX or ImSAnE [17, 51].

The flexibility and robustness of Zellige should allow to considerably relax the constraints that were previously imposed on the sample preparation and the image acquisition steps by the subsequent analysis. Indeed, Zellige can accommodate any number of surfaces in the acquired volume, overlapping or not, and of different contrast features. Zellige also showed excellent robustness against image noise. This should make it particularly useful in imaging contexts that are not easily amenable to automated analysis, such as intravital imaging. Finally, it is worth noting that Zellige is a generalist and modular method. With some adaptation of the surface voxel selection step, it could be used with imaging modalities beyond the scope of this article, for instance, in extracting the irregular and noisy surfaces of biological objects imaged with 3D electron microscopy images.

# Methods

## Human bronchial epithelium imaging

The data used in Fig. 5 were taken from the recent study [30] to which we refer for the preparation and imaging of human bronchial epithelium cultures. In brief, MucilAirTM were purchased from Epithelix (Saint-Julien-en-Genevois, France) and cultured for at least 4 weeks to reconstruct a differentiated human bronchial epithelium in vitro and stained as previously described. Images of the cultures were acquired using an inverted Zeiss LSM

710 confocal microscope controlled by the ZEN pro 2.3 software. Z-stack images of whole-mount samples were acquired with a Zeiss Plan Apochromat ×63 oil immersion lens (NA=1.4). The image used here was published in Robinot et al. [30] under the CC-BY-4 license.

## Drosophila imaging

Flies were raised at 25°C under standard conditions. Pupae were collected for imaging as described previously [52]. Ecad::GFP flies [53] were used for live imaging as previously described [1]. In brief, images were acquired with a spinning disk microscope from Gataca Systems driven by the MetaMorph software. The system is equipped with an inverted Nikon TI2E stand, a motorized XYZ stage, and a Nikon Plan Apo ×60 oil immersion (NA=1.4) lens and with a Prime95B camera.

## Cochlea imaging

The inner ears from wild-type (C57BL/6) mice were rapidly dissected from temporal bones at embryonic stages E14.5 in HEPES-buffered (10 mM, pH 7.4) Hanks' balanced salt solution and fixed in 4% paraformaldehyde, 1 h at room temperature. Specimens were permeabilized and stained for phalloidin-Atto 565 (Sigma) as previously described [44]. Fluorescence images were obtained with a swept-field confocal microscope (Opterra II) from Brucker. This system is equipped with a Nikon Plan Fluor ×60 oil immersion lens (NA=1.4).

## Inner ear organoid imaging

ESCs derived from blastocyst-stage embryos of R1 mice (mESCs) (ATCC, SCRC-1036) were maintained in feeder-free culture on 0.1% w/v gelatin (Sigma) coated substrates using LIF-2i medium as established previously [54]. The organoids were generated following the previously published protocol [49, 54]. Aggregates were harvested at day 14 and fixed in 4% v/v PFA (Electron Microscopy Sciences) overnight at 4 °C. After blocking (PBS; 10% v/v normal goat serum; 0.1% v/v Triton X-100), the aggregates were stained for phalloidin Atto 565 (1:1000) (Sigma) overnight at RT on a shaker and washed three times with PBS containing 0.1% v/v Triton X-100 for 1 h each at RT. Prior to imaging, the aggregates were incubated in a modified version of ScaleS solution containing 4 M Urea (Sigma), 40% w/v D-Sorbitol (Sigma), and 0.1% v/v Triton X-100, for 3–5 days to clarify the tissue. Finally, the aggregates were whole-mounted using the ScaleS solution and imaged using a confocal laser scanning microscope (A1R HD25, Nikon) equipped with a Nikon ×25 silicon oil immersion lens (NA=1.05).
# Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s12915-022-01378-0.

**Additional file 1: Supplementary Notes 1-5. Note 1.** Overview of Zellige implementation. **Note 2.** Generation of phantom 3D images. **Note 3.** Comparing the reconstructed and ground truth height maps. **Note 4.** Sensitivity Analysis. **Note 5.** Computational/processing time of Zellige. **Table S1.** Zellige parameters.

**Additional file 2: Figure S1.** Zellige implementation. (A) The two main algorithmic steps of Zellige. Upper part: Surface voxel selection step (step 1). Lower part: Surface assembly step (step 2). (B) Determination of the amplitude $A(\rho)$ of some local maximum $\rho$ along the z-axis, relative to its two closest local minima $p_1$ and $p_2$. (C,D) Connectivity rules used to connect putative surface voxels together for the construction of orthogonal surface elements (OSEs). These rules are based on a 6-connectivity relationship defined within each orthogonal (xz or yz) section as shown in (C). This relationship is extended to allow the presence of single straight gaps, while being constrained to forbid the occurrence of forking points. The allowed local neighborhood configurations along an OSE are illustrated in (D) in the case of a construction within xz sections. The OSEs are formally defined as the connected components of the graph $G_{OSE}$ defined by these rules within each orthogonal section. (E) Compatibility rules used in the surface assembly step, illustrated here in the case where OSEs have been constructed within xz sections, and assembly proceeds along the y-axis. To validate the addition of a new OSE $\sigma$ constructed within section y+1 (in blue) to the surface S under construction, whose intersection with section y is shown (surface line $l_y$ in green), two quantities are computed: the overlap $R(S,\sigma)$ is defined as the number of pixels of $\sigma$ that share the x coordinate of some pixel of $l_y$. The connectivity $C(S,\sigma)$ is defined as the fraction of the overlapping pixels of $\sigma$ that are 6-connected, in their respective yz section, to the corresponding point of $l_y$ (according to the 6-connectivity relationship shown in C, middle panel). In the depicted example, there are 2 OSEs (in green) with, for both of them, $R(S,\sigma)=3$ and $C(S,\sigma)=3/3=1$. The OSE will be added to S if $R(S,\sigma) \ge R_0$ and $C(S,\sigma) \ge C_0$, where $R_0$ and $C_0$ are tunable thresholds that set the stringency of the matching condition.

**Additional file 3: Figure S2.** Sensitivity analysis of Zellige on the phantom image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line.

**Additional file 4: Figure S3.** Sensitivity analysis of Zellige on the pupal fly image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line. (C) RMSE and coverage results obtained for different values of the $\sigma_{xy}$ blur taken within and outside the interval defined by the parameter sweep in (A, "$\sigma_{xy}$" panel). (D) Visualization of the projections on the extracted surfaces obtained for the different $\sigma_{xy}$ blur values shown in (C).

**Additional file 5: Figure S4.** Sensitivity analysis of Zellige on the cochlear epithelium image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line.

**Additional file 6: Figure S5.** Sensitivity analysis of Zellige on the primary epithelium culture image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line.

**Additional file 7: Figure S6.** Sensitivity analysis of Zellige on the inner ear organoid image. (A) Surface voxel selection parameters. (B) Surface assembly parameters. Reference values are indicated by the dashed line.

**Additional file 8: Figure S7.** Computational time analysis of Zellige. (A) Computation times (step 1 in blue, step 2 in brown) for processing the various images tested on a PC notebook computer with (processor Intel Core i9 2,4 GHz with 32 Gb of Ram). Apart for the case of a highly rough surface (culture specimen) the computation time is largely dominated by the surface voxel selection step (step 1). (B) The same computation times are re-plotted as a function of image size N (number of voxels). Note the linear growth of the computational time of step 1 as a function of N, while that of step 2 shows much slower growth.

## Acknowledgements

We thank Maia Brünstein of the Hearing Institute Bioimaging Core Facility - C2RT/C2RA and Florian Rückerl of the Photonic Bioimaging - C2RT/C2RA, for sharing their expertise on light microscopy. We thank Romain Levayer for sharing files and lab space with us. We thank Gilles Trébeau (https://fr.linkedin.com/in/gilles-trebeau-427601212) for helpful discussions regarding the Java frontend deployment of Zellige. We thank Lisa Chakrabarti, Rémy Robinot, and Vincent Michel for sharing data.

## Authors' contributions

CT designed the program and wrote the Java code. JYT supervised the Fiji plugin development and reviewed the Java code. JBM wrote the MATLAB scripts and generated all synthetic images. RE and GA performed some of the experiments presented here. CT, JMB, JYT, and RE jointly wrote the article. All authors participated in regular group discussions to develop the ideas presented here, which are the result of multiple cycles of tests on synthetic and biological data. The authors read and approved the final manuscript.

## Authors' information

Twitter handles: @jytinevez (Jean-Yves Tinevez) Twitter handles are transferred to the "Authors' information" section., @institutpasteur (the Institut Pasteur).

## Funding

This work was supported by Fondation pour l'Audition (FPA-IDA-STARTING-GRANT); Institut Pasteur (PTR#272); DIM ELICIT's grant from Région Ile-de-France (QP-IDF-DIM-ELICIT-2019); The French Agence Nationale de la Recherche (ANR-21-CE13-0038-013-01 and LabEx LIFESENSES ANR-10-LABX-65).

## Availability of data and materials

The project homepage below contains the source code, installation instructions, and documentation.

### Zellige

* Project name: Zellige.

* Project homepage: https://gitlab.pasteur.fr/ida-public/zellige-core

* URL for the Fiji plugin (Update > Manage update sites): https://sites.imagej.net/Zellige/

* Gitlab Branch: master

* Operating systems: Platform independent.

* Programming language: Java.

* Compiled in Java8

* Other requirements: Runs from Fiji [29].

* License: BSD 2

* Any restrictions to use by non-academics: License needed.

### Scripts to create Phantoms

* Project name: Phantoms.

* Project homepage / Archived version: https://doi.org/10.5281/zenodo.6414596

* Operating systems: Platform independent.

* Programming language: MATLAB.

* License: BSD 2

* Any restrictions to use by non-academics: License needed.

### Data sets

Image 3D stacks, ground truth height map, and height maps produced by Zellige are available on Zenodo under the CC-BY license: https://zenodo.org/communities/zellige/ [33, 39, 47, 48, 50].

## Declarations

### Ethics approval and consent to participate

RE holds a designer certificate of animal experimentation (level 2), allowing to perform experimental work on animals in strict accordance with the European directive 2010/60/EU, and French regulations. The Ethics Committee of the Institut Pasteur (Comité d'Ethique en Experimentation Animale - CETEA) has approved this study with the project identifier dha170006. This approval is based on careful compliance to the 3Rs principle in the care and use of animals (Annex IV - 2010/60/EU).

### Consent for publication

Not applicable.
## Competing interests

The authors declare that they have no competing interests.

Received: 5 May 2022 Accepted: 22 July 2022

Published online: 23 August 2022

## References

1. Etournay R, Popović M, Merkel M, Nandi A, Blasse C, Aigouy B, et al. Interplay of cell dynamics and epithelial tension during morphogenesis of the drosophila pupal wing. eLife. 2015;4:e07090. https://doi.org/10.7554/eLife.07090.

2. Guirao B, Rigaud SU, Bosveld F, Bailles A, López-Gay J, Ishihara S, et al. Unified quantitative characterization of epithelial tissue development. eLife. 2015;4:e08519. https://doi.org/10.7554/eLife.08519.

3. Etournay R, Merkel M, Popović M, Brandl H, Dye NA, Aigouy B, et al. TissueMiner: a multiscale analysis toolkit to quantify how cellular processes create tissue dynamics. eLife. 2016;5. https://doi.org/10.7554/eLife.14334.

4. Saadaoui M, Rocancourt D, Roussel J, Corson F, Gros J. A tensile ring drives tissue flows to shape the gastrulating amniote embryo. Science. 2020;367:453–8. https://doi.org/10.1126/science.aaw1965.

5. Butler LC, Blanchard GB, Kabla AJ, Lawrence NJ, Welchman DP, Mahadevan L, et al. Cell shape changes indicate a role for extrinsic tensile forces in drosophila germ-band extension. Nat Cell Biol. 2009;11:859–64. https://doi.org/10.1038/ncb1894.

6. Collinet C, Rauzi M, Lenne P-F, Lecuit T. Local and tissue-scale forces drive oriented junction growth during tissue extension. Nat Cell Biol. 2015;17:1247–58. https://doi.org/10.1038/ncb3226.

7. Kong D, Wolf F, Großhans J. Forces directing germ-band extension in drosophila embryos. Mech Dev. 2017;144:11–22. https://doi.org/10.1016/j.mod.2016.12.001.

8. Villedieu A, Bosveld F, Bellaïche Y. Mechanical induction and competence in epithelial morphogenesis. Curr Opin Genet Dev. 2020;63:36–44. https://doi.org/10.1016/j.gde.2020.03.003.

9. Prat-Rojo C, Pouille P-A, Buceta J, Martin-Blanco E. Mechanical coordination is sufficient to promote tissue replacement during metamorphosis in drosophila. EMBO J. 2020;39:e103594. https://doi.org/10.15252/embj.2019103594.

10. Blasse C, Saalfeld S, Etournay R, Sagner A, Eaton S, Myers EW. PreMosa: extracting 2D surfaces from 3D microscopy mosaics. Bioinforma Oxf Engl. 2017. https://doi.org/10.1093/bioinformatics/btx195.

11. Basu S, Rexhepaj E, Spassky N, Genovesio A, Paulsen RR, Shihavuddin A. FastSME: faster and smoother manifold extraction from 3D stack. In: 2018 IEEE/CVF conference on computer vision and pattern recognition workshops (CVPRW); 2018. p. 2362–23628. https://doi.org/10.1109/CVPRW.2018.00305.

12. Herbert S, Valon L, Mancini L, Dray N, Caldarelli P, Gros J, et al. LocalZProjector and DeProj: a toolbox for local 2D projection and accurate morphometrics of large 3D microscopy images. BMC Biol. 2021;19:136. https://doi.org/10.1186/s12915-021-01037-w.

13. Li K, Wu X, Chen DZ, Sonka M. Optimal surface segmentation in volumetric images-a graph-theoretic approach. IEEE Trans Pattern Anal Mach Intell. 2006;28:119–34. https://doi.org/10.1109/TPAMI.2006.19.

14. Lombardot B. Min cost Z surface projection Fiji. Java. 2017; https://imagej.net/Minimum_Cost_Zsurface_Projection.

15. Hall M, Frank E, Holmes G, Pfahringer B, Reutemann P, Witten IH. The WEKA data mining software: an update. ACM SIGKDD Explor Newsl. 2009;11:10–8. https://doi.org/10.1145/1656274.1656278.

16. Sommer C, Straehle C, Kothe U, Hamprecht FA. Ilastik: interactive learning and segmentation toolkit. In: 2011 IEEE international symposium on biomedical imaging: from Nano to macro; 2011. p. 230–3. https://doi.org/10.1109/ISBI.2011.5872394.

17. Heemskerk I, Streichan SJ. Tissue cartography: compressing bio-image data by dimensional reduction. Nat Methods. 2015;12:1139–42. https://doi.org/10.1038/nmeth.3648.

18. Haertter D, Wang X, Fogerson S, Ramkumar N, Crawford J, Poss K, et al. DeepProjection: rapid and structure-specific projections of tissue sheets embedded in 3D microscopy stacks using deep learning; 2021. https://doi.org/10.1101/2021.11.17.468809.

19. Green JBA, Sharpe J. Positional information and reaction-diffusion: two big ideas in developmental biology combine. Development. 2015;142:1203–11. https://doi.org/10.1242/dev.114991.

20. Miller CJ, Davidson LA. The interplay between cell signalling and mechanics in developmental processes. Nat Rev Genet. 2013;14:733–44. https://doi.org/10.1038/nrg3513.

21. Levin M, Martyniuk CJ. The bioelectric code: an ancient computational medium for dynamic control of growth and form. Biosystems. 2018;164:76–93. https://doi.org/10.1016/j.biosystems.2017.08.009.

22. Duclut C, Prost J, Jülicher F. Hydraulic and electric control of cell spheroids. Proc Natl Acad Sci. 2021;118. https://doi.org/10.1073/pnas.2021972118.

23. Ruiz-Herrero T, Alessandri K, Gurchenkov BV, Nassoy P, Mahadevan L. Organ size control via hydraulically gated oscillations. Development. 2017;144:4422–7. https://doi.org/10.1242/dev.153056.

24. Dumortier JG, Le Verge-Serandour M, Tortorelli AF, Mielke A, de Plater L, Turlier H, et al. Hydraulic fracturing and active coarsening position the lumen of the mouse blastocyst. Science. 2019;365:465–8. https://doi.org/10.1126/science.aaw7709.

25. Mammoto T, Mammoto A, Ingber DE. Mechanobiology and developmental control. Annu Rev Cell Dev Biol. 2013;29:27–61. https://doi.org/10.1146/annurev-cellbio-101512-122340.

26. Ladoux B, Mêge R-M. Mechanobiology of collective cell behaviours. Nat Rev Mol Cell Biol. 2017;18:743–57. https://doi.org/10.1038/nrm.2017.98.

27. Närhi K. Embryonic explant culture: studying effects of regulatory molecules on gene expression in craniofacial tissues. In: Seymour GJ, Cullinan MP, Heng NCK, editors. Oral biology: molecular techniques and applications. New York: Springer; 2017. p. 367–80. https://doi.org/10.1007/978-14939-6685-1_21.

28. Zeevaert K, Elsafi Mabrouk MH, Wagner W, Goetzke R. Cell mechanics in embryoid bodies. Cells. 2020;9:2270. https://doi.org/10.3390/cells9102270.

29. Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, et al. Fiji: an open-source platform for biological-image analysis. Nat Methods. 2012;9:676–82. https://doi.org/10.1038/nmeth.2019.

30. Robinot R, Hubert M, de Melo GD, Lazarini F, Bruel T, Smith N, et al. SARS-CoV-2 infection induces the dedifferentiation of multiciliated cells and impairs mucociliary clearance. Nat Commun. 2021;12:4354. https://doi.org/10.1038/s41467-021-24521-x.

31. Pietzsch T, Preibisch S, Tomančák P, Saalfeld S. ImgLib2—generic image processing in Java. Bioinformatics. 2012;28:3009–11. https://doi.org/10.1093/bioinformatics/bts543.

32. Otsu N. A threshold selection method from gray-level histograms. IEEE Trans Syst Man Cybern. 1979;9:62–6. https://doi.org/10.1109/TSMC.1979.4310076.

33. Boutet de Monvel J, Trébeau C, Altay G, Tinevez J-Y, Etournay R. Zellige example dataset: synthetic image dataset. 2022. https://doi.org/10.5281/zenodo.6376584.

34. Boutet de Monvel J. Scripts to create phantom images containing surfaces. Zenodo. 2022. https://doi.org/10.5281/zenodo.6414596.

35. de Chaumont F, Dallongeville S, Chenouard N, Hervé N, Pop S, Provoost T, et al. Icy: an open bioimage informatics platform for extended reproducible research. Nat Methods. 2012;9:690–6. https://doi.org/10.1038/nmeth.2075.

36. Jülicher F, Eaton S. Emergence of tissue shape changes from collective cell behaviours. Semin Cell Dev Biol. 2017;67:103–12. https://doi.org/10.1016/j.semcdb.2017.04.004.

37. Morata G, Lawrence P. An exciting period of drosophila developmental biology: of imaginal discs, clones, compartments, parasegments and homeotic genes. Dev Biol. 2022;484:12–21. https://doi.org/10.1016/j.ydbio.2022.01.008.

38. Ray RP, Matamoro-Vidal A, Ribeiro PS, Tapon N, Houle D, Salazar-Ciudad I, et al. Patterned anchorage to the apical extracellular matrix defines tissue shape in the developing appendages of drosophila. Dev Cell. 2015;34:310–22. https://doi.org/10.1016/j.decelv.2015.06.019.

39. Trébeau C, Boutet de Monvel J, Altay G, Tinevez J-Y, Etournay R. Zellige example dataset: drosophila pupal wing and abdomen; 2022. https://doi.org/10.5281/zenodo.6376566.

40. Zallen JA. Planar polarity and tissue morphogenesis. Cell. 2007;129:1051–63. https://doi.org/10.1016/j.cell.2007.05.050.
41. Montcouquiol M, Rachel RA, Lanford PJ, Copeland NG, Jenkins NA, Kelley MW. Identification of Vangl2 and Scrb1 as planar polarity genes in mammals. Nature. 2003;423:173–7. https://doi.org/10.1038/nature01618.

42. Cohen R, Amir-Zilberstein L, Hersch M, Woland S, Loza O, Taiber S, et al. Mechanical forces drive ordered patterning of hair cells in the mammalian inner ear. Nat Commun. 2020;11:5137. https://doi.org/10.1038/s41467-020-18894-8.

43. McKenzie E, Krupin A, Kelley MW. Cellular growth and rearrangement during the development of the mammalian organ of Corti. Dev Dyn. 2004;229:802–12. https://doi.org/10.1002/dvdy.10500.

44. Etournay R, Lepelletier L, Boutet de Monvel J, Michel V, Cayet N, Leibovici M, et al. Cochlear outer hair cells undergo an apical circumference remodeling constrained by the hair bundle shape. Dev Camb Engl. 2010;137:1373–83. https://doi.org/10.1242/dev.045138.

45. Driver EC, Northrop A, Kelley MW. Cell migration, intercalation and growth regulate mammalian cochlear extension. Dev Camb Engl. 2017;144:3766–76. https://doi.org/10.1242/dev.151761.

46. Qian D, Jones C, Rzadzinska A, Mark S, Zhang X, Steel KP, et al. Wnt5a functions in planar cell polarity regulation in mice. Dev Biol. 2007;306:121–33. https://doi.org/10.1016/j.ydbio.2007.03.011.

47. Trébeau C, Boutet de Monvel J, Altay G, Tinevez J-Y, Etournay R. Zellige example dataset: mouse embryonic cochlea stained for F-actin; 2022. https://doi.org/10.5281/zenodo.6376542.

48. Robinot R, Chakrabarti L, Michel V, Trébeau C, Boutet de Monvel J, Altay G, et al. Zellige example dataset: primary culture of human bronchial cells infected by SARS-CoV-2; 2022. https://doi.org/10.5281/zenodo.6376594.

49. Koehler KR, Mikosz AM, Molosh AI, Patel D, Hashino E. Generation of inner ear sensory epithelia from pluripotent stem cells in 3D culture. Nature. 2013;500:217–21. https://doi.org/10.1038/nature12298.

50. Altay G, Trébeau C, Boutet de Monvel J, Tinevez J-Y, Etournay R. Zellige example dataset: inner ear organoid epithelium; 2022. https://doi.org/10.5281/zenodo.6376582.

51. Barbier de Reuille P, Routier-Kierzkowska A-L, Kierzkowski D, Bassel GW, Schüpbach T, Tauriello G, et al. MorphoGraphX: a platform for quantifying morphogenesis in 4D. eLife. 2015;4:e05864. https://doi.org/10.7554/eLife.05864.

52. Classen A-K, Aigouy B, Giangrande A, Eaton S. Imaging drosophila Pupal wing morphogenesis. In: Dahmann C, editor. Drosophila: methods and protocols. Totowa: Humana Press; 2008. p. 265–75. https://doi.org/10.1007/978-1-59745-583-1_16.

53. Huang J, Zhou W, Dong W, Watson AM, Hong Y. Directed, efficient, and versatile modifications of the drosophila genome by genomic engineering. Proc Natl Acad Sci. 2009;106:8284–9. https://doi.org/10.1073/pnas.0900641106.

54. Nie J, Koehler KR, Hashino E. Directed differentiation of mouse embryonic stem cells into inner ear sensory epithelia in 3D culture. In: Organ regeneration. New York: Humana Press; 2017. p. 67–83. https://doi.org/10.1007/978-1-4939-6949-4_6.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

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