We leveraged five public datasets that are available online in the paleo.esrf.eu database, a public resource for palaeontological volumes scanned at the European Synchrotron Radiation Facility (ESRF).

The datasets used in this study represent a range of fossil types and scanning conditions, each presenting its unique challenges to segmentation. The Burrow dataset, from the Early Triassic Karoo Basin in South Africa, contains complete skeletons of a Thrinaxodon and a Broomistega and was scanned at 45.5 µm resolution by Fernandez et al., (2013). This dataset benefits from a relatively clear contrast between the fossils and the surrounding matrix, allowing for high segmentation accuracy. In contrast, the Middle Triassic Nothosaurus skull from the Muschelkalk of the Netherlands scanned at 12.82 µm resolution (Voeten et al. 2018) presents lower contrast. This makes it more difficult to distinguish between fossilised bone and matrix. The Nothosaurus stack was also post-processed into a recoded version of the three-dimensional data based on local texture complexity to reveal low-contrast features (Voeten et al. 2018). The end-Cretaceous paddlefish Parapsephurus willybemisi, from the Tanis deposit in North Dakota, was scanned at 45.53 µm and 13.67 µm resolution (During et al. 2022 and Hilton et al. 2023). These datasets highlight the role of scan resolution in segmentation accuracy. The higher-resolution ∼13 µm dataset allows for the segmentation of more detailed phenomena, while the lower-resolution ∼45 µm dataset does not, as key anatomical features become more difficult to discern. Finally, an Early-Cretaceous squamate embryo inside an egg from Thailand (Phu Phok) was scanned at 5.06 µm resolution by Fernandez et al. (2015). This dataset presents a unique challenge due to the fine structures, compounded by low contrast, leading to a significant challenge in segmentation for both humans and machine learning-based methods. For additional details about the subsamples of these datasets used for our evaluation, see Table 1.

The ml4paleo web application comprises four components; one main application server, and three task runners that operate on an on-disk queue.

Application server : The web application server is written in Flask (Grinberg 2018). The server can use an on-disk serialisation format for storing dataset metadata, or it can be configured to run a standalone database daemon. By default, our application uses a standalone JSON file as a metadata store to reduce the software dependencies of running a minimal deploy.

Conversion queue : Researchers commonly receive concatenated image stacks as output from the synchrotron facility. Such image stacks tend to be too large and inefficient to access for common machine learning workflows. Thus, the first memory-intensive stage of processing is to convert this inefficient data format into a chunked data format. The user can easily upload their dataset in most common formats, including image stacks and some common chunked formats. Chunking high-dimensional data involves re-organising a large dataset of 2D images into smaller, more compact 3D (volumetric) or 4D (volumetric time series) pieces, known as chunks. These chunks can then be treated as independent small volumes or manipulated in spatial sequence. This conversion enables the ml4paleo tools to operate on sub-volume cut-outs without loading complete slices into memory. It is this memory efficiency that enables our infrastructure to work equally well on low-performance and high-grade computers alike.

As a core design principle for this software package, we aimed to select default algorithms that impose minimal restrictions on user hardware, ensuring accessibility for a wide range of systems. In other words, our goal was to accommodate the lowest common denominator in terms of hardware capabilities. This required our segmentation algorithms to be compatible with consumer-grade hardware while operating at the necessary spatial scales.

Two primary decisions emerged from this requirement: first, the dimensionality of the segmentation process, and second, the selection of machine learning models, based on both computational demands and dimensionality considerations.

While 2D stacked segmentation, where individual slices are segmented and later reconstituted into 3D, can deliver acceptable results in many cases, it is inferior to true 3D segmentation. The latter benefits from the inclusion of additional spatial context within the 3D data stack. However, memory limitations significantly reduce the efficiency of 3D segmentation, as its computational complexity scales with N³ compared to N² for 2D segmentation. Given that X-ray tomography (XRT) data is isotropic, with uniform voxel sizes across all dimensions, 2D models have proven sufficient for our purposes thus far. Nevertheless, our codebase includes 3D models, allowing users to switch between 2D and 3D approaches as needed (see Data & Code Availability).

As the default segmentation algorithm, we utilise a simple Random Forest (RF) model trained on heuristic-based pixel features (Pedregosa et al. 2011). This is implemented using the "skimage.feature.multiscale_basic_features" module with default parameters, all of which are adjustable through the administrator settings of the web application. The feature extractor relies on basic image properties such as texture, intensity, and edge detection, progressively applying Gaussian filters to estimate these features on ”downsampled” versions of the images.

For training efficiency, we optionally provide a parameter to subsample pixels, using every nth pixel for training rather than all annotated pixels. By default, 50% of the available image/segmentation pairs are used during the training process.

Our codebase also includes a convolutional residual U-Net (Ronneberger et al. 2015), which requires a GPU for efficient training, though it otherwise operates with modest computational requirements. These tools are implemented through the "ml4paleo.segmentation.Segmenter3D" interface. Both this and the "ml4paleo.segmentation.Segmenter2D" interface are designed to support easy integration of more advanced or specialised segmentation algorithms by third-party contributors.

In order to process arbitrarily large volumes of 3D imagery, it is necessary to decrease the size of the volume that is fed into the pipeline. In volumetric data processing, this is referred to as “chunking” the volume into subvolumes with contiguous areas nearby in memory. Several data standards have been built in a variety of scientific domains to meet this need (Miles et al. 2020; Matelsky et al. 2021; Hider et al. 2022). For most purposes, interchangeable data standards make the distinctions between these chunked formats insignificant, and for this project, we have used the standard “zarr” format (Miles et al. 2020). This chunking technique makes it possible to independently segment subvolumes of the dataset, which can then be combined in a segmentation fusion step later. In their corresponding Python Application Programming Interfaces (APIs), these “datastores” tend to meet the standard numpy-like slicing interface in order to be maximally interoperable with the machine-learning ecosystem (Harris et al. 2020). This also means that data can be lazily loaded from disk so that the majority of the volume can remain compressed and at rest while subvolumes are accessed in memory.

The fusion step can be quite nuanced: in some cases, a single pixel may be justifiably assigned to more than one segmentation type, and therefore it is generally advised to store segmentation as a multi-channel output. While this enables higher-quality segmentation, it comes at the cost of greatly increased storage space, as each segmentation mask must be stored independently. Because our intention was to minimise the computational and storage costs for new projects, we opted here to store single-channel segmentation by default, though this is configurable by the administrator.

We provide a simple web application to generate human expert ground truth to train the segmentation pipeline. This tool (Fig. 2) provides a resizable brush for short annotation tasks of small 2D slices. The user can adjust the contrast and scroll several slices up and down in the stack to provide spatial context during annotation. The tasks are designed to be very short in duration to avoid attentional errors and to make exhaustive annotation feasible on resource-constrained hardware, following design principles from Bishop et al. (2021).

Figure 2. 

The annotation tool on an excerpt slice of the Burrow dataset (screenshot). The user is given a resizable brush to manually segment fossils from the matrix in a small 2D area. Here, shown in red, the ml4paleo automated segmentation has already been run on this slice, so the user only needs to proofread the machine’s work, by adding with a brush or removing with an eraser, rather than annotate de novo.

Pairs of human-annotated segmentation masks and the original corresponding imagery are then fed to the selected machine learning model for training. By default, ml4paleo uses a random forest with local pixel-neighbourhood features (Pedregosa et al. 2011), though the segmentation algorithm in use can be replaced by a simple U-net neural network (Ronneberger et al. 2015) or other machine learning models that implement our “ml4paleo.segmentation.Segmenter3D” interface. A checkpoint of the trained model is then saved for reuse under a unique identifier that contains both the timestamp of training and the specific training samples used. This means that as more data is generated, a library of models can be accumulated to enable analysis and understanding of the relationship between training data volume and segmentation quality.

Segmentation queue : Once a machine learning model has been trained for image segmentation, the dataset is queued for segmentation. In this stage, each worker loads a subvolume of the imagery, performs a dense segmentation, and saves a new segmentation mask volume to disk.

This new layer is identified by the unique model identifier that was generated during the model training; thus, segmentation layers can be uniquely associated with the model that generated them, for provenance and reproducibility.

Mesh queue : As a postprocessing stage after segmentation, a user may opt to download their dataset as a mesh in OBJ or STL formats for rendering in 3D software. This mesh process also operates on one chunk at a time, after which the sub-meshes are stitched into one large mesh.

Visualisation in three dimensions : To visualise the scan as well as the annotations and segmentations, we utilise “Neuroglancer” (Fig. 3), a high-performance, WebGL-based visualisation framework developed by the Google Connectomics Team, designed for efficiently viewing and analysing large-scale volumetric neuroscience datasets (Charles et al. (2020). A deployment of Neuroglancer can be packaged with deployments of the ml4paleo software suite, or any public Neuroglancer application can be used and pointed at ml4paleo data volumes through our volumetric API.

Figure 3. 

Browser-based visualisation. This stage in the workflow uses Neuroglancer to render imagery and annotations. This figure shows a user visualising the Parapsephurus willybemisi 13.67 µm resolution scan (During et al. 2022) in interactive 3D in a modern web browser. This tool is included with the ml4paleo web suite.

Once the user has produced at least one saved model — requiring a minimum of one annotated slice — they can optionally allow the machine to produce an initial hypothesis segmentation before each subsequent annotation task. This has the potential to greatly reduce the labour cost of individual annotation tasks, since the machine may produce a nearly complete segmentation mask, depending on the image properties of the underlying volume. Furthermore, the available eraser tool enables users to refine segmented slices, significantly reducing the effort compared to re-segmenting from scratch. Through this machine guidance, the user will become aware of the qualitative performance characteristics of the segmentation model and can monitor its improvement throughout the annotation process.

In our approach, for the datasets tested in this study, a total of 20 slices (from patches of full slides) were initially annotated by the same expert palaeontology segmentation author (MD) for each dataset using the annotation web application. However, due to the specific challenges posed by the Phu Phok and Nothosaur datasets, additional annotations were required to improve to the required level of performance. These properties and improvements are reported below.

Fossil embryos, such as those in the Phu Phok dataset, are typically poorly ossified (Fernandez et al. 2015), resulting in significant contrast variations between bones, soft tissues, and matrix. This contrast variability makes it very challenging for the human eye to differentiate structures. The same difficulty is encountered by automated segmentation systems, which struggle with this inconsistent presentation of biologically similar structures. Likewise, the Nothosaur dataset contained a crack running through the fossil’s skull, which created areas of phase contrast that obscured the underlying bone texture. Cracks like these are common in fossil scans and pose difficulties for both human and AI segmentation efforts.

Due to these challenges, after the initial 20 slices were annotated for both datasets, an additional 10 slices were segmented and the model was retrained. A final set of 5 slices was then annotated, resulting in a total of 35 annotated slices for the Nothosaur and Phu Phok Egg datasets to achieve satisfactory segmentation results. In our experiments, switching the default segmentation tool from the random forest to a more sophisticated (neural) model resolved many of these failures, at the cost of additional execution time and GPU requirements.

After the user is satisfied with the model performance, they can deploy the segmentation model on the complete data volume. Below, we provide quantitative evaluation of the ml4paleo web application default segmentation, showing performance metrics for leave-one-out cross-validated pairs of human-annotated segmentation and imagery. A visualisation of the ground truth versus the prediction is provided in Fig. 4.

Figure 4. 

Performance on the burrow dataset. We picture one 2D slice of the original volume (top left), as well as human (top right) and machine annotated segmentation masks (bottom right). On the bottom left, we render the residual errors.

Because ml4paleo datasets are stored in a chunked volumetric data format, it is easy to parallelise an operation across the dataset in small sub-volume increments. Using dataset slicing code borrowed from the electron microscopy neuroscience community (Matelsky et al. 2021), we run the latest model checkpoint from the training stage on each sub-volume of data. When processed in memory, the sub-volumes do not have to match the default size of the chunked format subvolumes on disk (128 × 128 × 128). Our informal testing suggests that using a larger segmentation sub-volume size of 256 × 256 × 256 can significantly improve throughput without compromising segmentation quality. However, this depends heavily on the available computing resources: memory-constrained environments may benefit from smaller segmentation chunk sizes, while environments with limited CPU cores may perform better with larger chunk sizes.

It is beneficial to have a small overlap of the processed chunks so that the computed pixel-wise features on the edges of each chunk accurately reflect their spatial context (i.e. in the “middle” of the volume, not on an edge). However, we discovered in our tests that this was not always necessary, especially in more easy-to-segment datasets. Samples with a clear brightness contrast between matrix and fossil tended not to need additional context, and fossils distinguished from the matrix only by texture tended to need a broader margin of additional context from neighbouring chunks. A further evaluation and benchmark of this dataset property is left for future work.

To better understand the segmentation performance, we compared it qualitatively to prior segmentations by the original researchers (Fernandez et al. 2013, 2015; Voeten et al. 2018; During et al. 2022) and quantitatively with our own cross-validation evaluation suite. We report the following:

Precision measures the proportion of correctly identified fossil pixels relative to all pixels predicted as fossil by the model. High precision indicates that the model made fewer false positive predictions.

Recall assesses the model’s ability to identify all actual fossil pixels within the dataset. A higher recall means that the model successfully identified more true fossil pixels.

F1-score (Dice coefficient) is the harmonic mean of precision and recall, providing a balanced measure of segmentation quality. It is particularly useful when precision and recall values are unbalanced.

Accuracy measures the proportion of correctly classified pixels (both fossil and background), a high accuracy therefore reflects a low systematic error.

Jaccard index, also called intersection over union, measures similarity between the ground truth and the prediction.

The segmentation model demonstrated strong performance, with high precision and recall, suggesting a robust ability to correctly identify and segment relevant features. The accuracy was also high (0.982), while the Jaccard Index, a measure of similarity between predicted and true masks, was 0.759, indicating a solid overlap between machine predictions and human annotations (Fig. 6).

Figure 6. 

Burrow segmentation samples. The high contrast of the fossil versus the matrix was well-suited for the simple multiscale pixel features used by our default model.

Paddlefish ∼13 µm Dataset (Fig. 7): For this high-resolution scan, performance was moderate, with precision (0.583) and recall (0.492), resulting in an F1-score of 0.488. The segmentation shows a reasonable but not perfect overlap between predicted and true segmentation masks.

Figure 7. 

Paddlefish ∼13 µm segmentation samples. While the segmentation tends to include the right gross morphology, the model notably failed to segment the internal volume of bones correctly. This interestingly has the effect of producing visually “correct” meshes, since the volume of the segmentation is irrelevant to the surface of the generated OBJ or STL meshes.

Paddlefish ∼45 µm Dataset (Fig. 8): Performance on this lower-resolution scan was weaker overall, with precision (0.718) and recall (0.191), resulting in an F1-score of 0.185. However, the model achieved high accuracy (0.979), suggesting it effectively avoided false positives, though at the cost of lower sensitivity.

Figure 8. 

Paddlefish ∼45 µm segmentation samples. The model tended to confuse bone for other high-brightness objects in the scene, likely due to underparametrisation and an inability to distinguish bone from other bright sources using the necessary 3D context.

The model’s performance declined significantly on this dataset, achieving precision (0.602) and recall (0.496). The low recall suggests that the model struggled to capture true positives. Despite this, accuracy remained reasonably high, likely due to the background dominating the images (Fig. 9).

Figure 9. 

Nothosaur segmentation samples. Two samples of the Nothosaur dataset, illustrating that the texture of the imagery (left panels) was beyond the capabilities of a simple random forest segmentation model.

This dataset posed the greatest challenge, with extremely low precision which was improved by annotating more slices (n = 20 0.232, n = 35 0.401) and recall (n = 20 0.06, n = 35 0.02). Visually, the segmentation quality by the simple random forest was poor overall, reflecting the difficulty encountered by human annotators alike (Fig. 10).

Figure 10. 

Phu Phok Egg segmentation samples. This challenging dataset illustrates a major failure of the default random forest segmentation model with the multiscale basic pixelwise features in 2D. The high phase contrast in the scan likely caused the model to excessively interpret areas of strong phase contrast as bone, leading it to incorrectly identify other high-contrast edges, such as the boundary of the sample, rather than the low-contrast bones.