Data Acquisition and ProcessingSpecimens

The brains used in this paper were provided by the Brazilian Neurobiodiversity Network Initiative (https://neurobiodiversidade.org), recently established to remediate the lack of neuroanatomical descriptions in non-model animals. One brain was collected from each of three different species: Guiana dolphin (Sotalia guianeneis), Franciscana dolphin (Pontoporia blainvillei) and Steller sea lion (Eumetopias jubatus).

The two dolphins used in this study were collected with authorization from the Chico Mendes Institute for Biodiversity Conservation (ICMBio) and the Brazilian Institute of Environment and Renewable Natural Resources (IBAMA). The research was approved by the Committee on Ethical Animal Use of the Science Center of the Federal University of Rio de Janeiro (process number 01200.001568/2013-87).

The Steller sea lion was hunted as part of pest control at the coast of Shakotan (Hokkaido, Japan; Matsuda, 2021) in accordance with Japanese laws and regulations, a decision unrelated to the current study.

MRI Acquisition

For the two cetaceans, MRI images were generated in a 7 Tesla Scanner (Classic Magneton, Siemens Healthcare, Erlangen, USA) equipped with a 32-channel head coil (Nova Medical, Wilmington, USA) at the "Image Platform in the Autopsy Room" (PISA) facilities of the School of Medicine at the University of São Paulo.

For the Steller Sean Lion, brain scans were generated in a 3 Tesla Scanner (Magnetom Prisma, Siemens Health Care, Erlangen, Germany) with a 16-channel body coil at the facilities of the Medical School of Hokkaido University.

The MRI images were acquired by using the following protocol:

Guiana dolphin: 3D-SPACE Repetition Time (TR) 2000 ms; Echo Time (TE) 109 ms; 0.38 mm isotropic spatial resolution; Field of View (FOV) 192x144; Bandwidth (BW) 610 Hz/px; Flip Angle (FA) 120; Acquisition Time (AT) 1 h 45 min (Avelino-de Souza et al., 2024);

Franciscana dolphin: 3D-SPACE TR 2000 ms; TE 115 ms; 0.27 mm isotropic spatial resolution; FOV 102x102; BW 372 Hz/px; FA 120; AT 2 h 5 min;

Steller sea lion: Axial T2 space TR 4000 ms; TE 410 ms; 0.78 mm isotropic spatial resolution; FOV 237x200; BW 700 Hz/px; FA 120; AT 12 min 24 s.

Processing

After image acquisition, images were used to obtain the contours of anatomically-defined regions of interest (ROI). This so-called segmentation consists of manually tracing the perimeters surrounding such areas (Fig. 1).

All brains were segmented by a single anatomist with the aid of a Wacom Cintiq Pro 13 (Wacom, Saitama, Japan) and the medical image software Osirix MD v.12.0.3 (Pixmeo, Geneva, Switzerland; Rosset et al., 2004). For S. guianeneis, the segmentation was performed every 5 slices, leaving a gap of 1.90 mm between consecutive segmentations. For P. blainvillei, the same spacing of 5 slices was used, leaving a gap of 1.33 mm. Lastly, E. jubatus had all its slices segmented, leaving a gap of 0.78 mm, the same value as the isovoxel length.

Fig. 1

Examples of how the MRI and the manual segmentation of GM, WM and exposed surface ROIs look like for the three species in this study. A Steller sea lion. B Franciscana dolphin; adapted from Avelino-de-Souza (2023). C Guiana dolphin; adapted from Avelino-de Souza et al. (2024)

The manual segmentation data were then rearranged into two groups: full-resolution, when all available slices are used in the reconstruction, and low-resolution, when every other segmentation of the full-resolution is used. The low-resolution case is intended to mimic a scenario where data for cortical reconstruction is sparse.

Each individual segmentation was then saved as a JSON files, which were used as input in the reconstruction that followed. The full process is exemplified in Fig. 2.

Fig. 2

Visual guide of the data preparation for processing on Stitcher. A Example of a segmentation of gray matter, white matter and exposed surface ROIs; adapted from Avelino-de-Souza et al. (2024). B Sagittal view of the same brain; adapted from Avelino-de-Souza et al. (2024). The orange stripe shows the the first slice of segmentation and the red arrow indicates the direction of subsequent slices that must be segmented. C Stack of segmented slices, at full resolution (left), using all available slices; and at low resolution (right), using every other slice

Stitcher

The Stitcher tool is an innovative software package for surface reconstructions, utilizing manual segmentation of contours in slice stacks, and generating a self-consistent lateral interpolating mesh surface as final output. This package (available at https://github.com/hmynssen/Stitcher) is structured as an open-source python library. Upon installation, it will provide 3 classes that can be employed together to 3D-reconstruct the brain. The core of the code implements a 2D search to find a minimum cost path that corresponds to a triangulation between contours in two consecutive slices (Fuchs et al., 1977), therefore creating a stitching pattern that connects all the vertices in both slices while preserving geometric and topological integrity. We use as the cost function the total edge length (therefore generating ’stretchy’ interpolation surfaces, analogous to cling film); but other functions, such as total surface area, can be easily implemented.

The minimum cost path might not be unique but all reconstructed surfaces corresponding to such absolute minimum path will be different triangulations of the same surface: think of a set of four points, two on each slice, forming a square being. Triangulation using either diagonal would be equally valid, but the final surface would not be exactly the same. In this sense, we can call Stitcher a deterministic reconstruction method for always creating equivalent reconstructions.

The first class, Point class, stores 3 number parameters to represent spatial coordinates for each vertex, and contains basic vector operations that simplify the code, such as vector addition, multiplication by scalar and dot product.

The next class is Perimeter, which contains a numpy (Harris et al., 2020) array storing a collection of Points corresponding to the contour exported by the segmentation tool. This sequence of N Points correspond to a closed polygon of \(N-1\) edges and N vertices that encloses one of the three desired regions: GM, WM or exposed surface. From all the methods included in the Perimeter, three can be highlighted: Perimeter addition, that merges two polygons by connecting them in the shortest distance possible; plane self-avoiding correction, that removes edges that cross each other and creates the necessary new edges to preserve the polygon continuity; and orientation flip, that invert the order of the vertices and is used to guarantee that every polygon will be oriented clock-wise.

Lastly, the Surface class is centered around the method build_surface which performs the following algorithmic steps:

Algorithm 1

Class method build_surface logic

However, the minimum cost (5) may create self-intersections, which cannot exist in the real surface of a brain. For this reason, Stitcher imposes an extra condition in the search algorithm:

Algorithm 2

The retraction (6) is performed by removing the last entry from the path and setting it to have an infinite cost. If too many retractions are performed, as set by an end-user criteria, a new starting point is chosen. Usually, the starting point, i.e., the first element of the Cost Matrix, is the closest distance between the above and below Perimeters. Due to the complexity of the brain shape, this starting point sometimes is not the best choice, leading to a mesh that would necessarily have a self-intersection. For that reason, a sorted list of of closest to furthest pairs of Points is computed and used in step 8 if necessary. The end-user may also provide the initial pair of points if so desired.

After step 5 from the first algorithm is done, step 8 is used to create a closed surface. Since the segmentation step only traces perimeter of a polygon, at least the final and last Perimeters must be tiled with a flat triangular mesh to prevent a hole from forming. Additionally, typically in the interior slices one or more splits/bifurcations of the segmentation contours will occur, two contours in one slice will need to be stitched self-consistently to a single contour in the next (much like trousers). Part of the single contour will need to be tiled with a triangular mesh to prevent a hole for forming.

Application of Stitcher for Subcortical and Non-Cortical Structures

Stitcher is capable of building any lateral surface given two polygons in parallel planes (Fig. 3). It can also reconstruct the WM surface, subcortical structures and non-cortical structures. As an example, we reconstruct the ventricle of the Sotalia guianensis, and calculate its area and volume. This multi-lobed but otherwise smooth structure poses no challenges to the algorithm.

Fig. 3

Example of a stack of arbitrary quadrilaterals (left side) and its reconstruction using Stitcher algorithm (right side)

Lastly, by using an atlas, an anatomist could identify structures of very low contrast, which are an even greater challenge for automatic segmentation tools. In this paper, superior and inferior colliculi, hippocampus, amygdala and thalamus from Sotalia guianensis were reconstructed to demonstrate this type of application.

Manual Corrections

The Stitcher algorithm ultimately finds a self-avoiding mesh but that does not assure the absence of artifacts. In fact, one type of artifact can occur in highly gyrified brains that contain many perimeters/segmentation per MRI slice. Called knot artifacts, they can be easily identified visually as in Fig. 4.

Fig. 4

Knot artifact example: the final mesh (top) containing the artifact and the same mesh (bottom) after the manual correction of the input data

These artifacts increase the surface area and do not represent the real topology of the brain. However, the cause is not the stitching process itself, instead it is either a poorly performed segmentation or a failure in the algorithm of merging two Perimeters together, as in Fig. 5. The former can be fixed with a consistent segmentation of the gyri and sulci, while for the latter, new points must be inserted to create a bridge indicating where the connection between the two Perimeters should occur.

Fig. 5

Example of inconsistent segmentation of a human subject gray matter. Two consecutive slices (top and bottom) have the beginning of a deep sulci at opposing sides (indicated by the arrows)

Post-Processing

The final reconstruction from Stitcher can be analyzed and manipulate in any mesh processing software. Here we used Meshlab (Cignoni et al., 2008) with the following pipeline to extract the geometrical value of volume and area of GM, WM and Exposed surfaces:

The first two steps ensure that the surface is orientable and closed, which are necessary respectively for area and volume estimation. But the resulting raw surface will be at this point jagged in a ladder-like manner, which is clearly an artifact of the underlying voxelization. This will distort the area estimation, so for a more reliable result, one needs to smooth away the jagged edges. To that end, we first perform an isotropic remeshing (Hoppe et al., 1993) to re-tile the surface with more regularly-sized triangle edges; this does not alter areas or volumes but ensure that the ensuing smoothing, acting on the mesh vertices, happens homogeneously and isotropically over the entire surface. For the smoothing, we apply first a volume-preserving HC Laplacian flow (Vollmer et al., 1999) that diffusively evens out the positions of neighboring vertices; then a Taubin flow (Taubin, 1995), a modified Laplacian flow that evens out Gaussian curvature.

Visual inspection clearly indicates that the smoothed reconstructed surfaces are much closer to the original cortical surface than the raw ones, with their characteristic ladder-like pattern. The Smooth Surface remains closed and thus can have its geometrical values extracted again. Meshlab provides a python library, pymeshlab (Muntoni & Cignoni, 2021), making the above processing programmable via simple scripting. Here we used the graphical interface for visual inspection and pymeshlab for the all numerical data processing.

Stitcher Usage

Stitcher relies only on sorting the stack of perimeters to start working. This can be done automatically by an adapted shadow casting algorithm that projects the Perimeters contained in one slice into the parallel slice below. All overlapping shadows should be grouped together to be reconstructed by Stitcher.

If this procedure is applied recursively through all pairs of segmentation slices, Stithcer will be able to reconstruct the ROI surface. Thus, Stitcher is able to produce a surface even if not all MRI slices of an acquisition sequence have been manually segmented. However having only a few segmented slices may reduce the anatomical information below the threshold necessary to properly reconstruct the structure. Thus, the number of slices that need to be traced must be determined by the anatomist and the end-user to balance the necessary detail with the required tracing workload.

This pipeline was implemented for the reconstruction of all cortical surfaces, i.e., GM, WM and exposed surfaces, at both full- and low-resolution for comparison purposes. Also, the same pipeline was applied for Guiana dolphin’s hippocampus, amygdala, thalamus and ventricles.

AnalysisDirect Measurements and Reference Values

The reconstruction results will be analyzed visually and numerically. For the GM surface, we extracted the area \(A_t\) (total area) and the volume \(V_t\) (total volume). For the other two, only one of the geometrical measures was needed. For WM, we only used the volume \(V_\), and for the Exposed surface, the area \(A_e\). The cortical thickness T was computed from GM and WM as:

$$\begin \begin T = \frac} . \end \end$$

(1)

Also, we compute the gyrification index (GI) as the ratio between \(A_t\) and \(A_e\).

In the absence of a full surface reconstruction of the brain regions of interest, it is possible to approximate (Ribeiro et al., 2013) their volume V and lateral surface area A using analytical formulas that take into account only the areas and perimeters of the contours in the slice stack:

$$\begin \begin A = \sum _\-S_n)^2 + [h(P_n+P_)/2]^2\}^ \end \end$$

(2)

$$\begin \begin V = \sum _ h[S_n+S_+(S_n\times S_)^]/3 \end \end$$

(3)

where \(S_n\) and \(P_n\) are the area and perimeter of a polygonal contour in slice n and h is the separation between adjacent slices. The Analytical Approximation (AA) process can thus be used in segmentation from consecutive MRI slices to approximate the value of \(A_t\), \(A_e\) and T,. This procedure has been applied to human cortices, and results were shown to be close to the corresponding real values (Ribeiro et al., 2013). As the gyrification increases, the cross-sectional images from MRI may display a series nested perimeters. Consider, for instance, an annular configuration in slice n of two nested contours, so that a region of (say) WM lies entirely inside GM; but within the WM region lies a smaller island of GM. This slice will then have two disjoint and closed contours of WM-GM boundary, each with an area and perimeter. In this case, the WM area \(S_n\) and perimeter \(P_n\) of the annular WM region to be used in Eqs. 2 and 3 must now read

$$\begin \begin S_n = S_ - S_ \end \end$$

(4)

$$\begin \begin P_n = P_ + P_\text \end \end$$

(5)

where the n1 and \(n_2\) indices refer respectively to the outer and inner contour. More generally, if a slice includes multiple disjoint contours, their perimeters should always be added, and their areas summed algebraically with signs given by the parity of the degree of nesting.

This is a complication that usually does not arise often in the human and terrestrial mammal cortices we analyzed, but that clearly has to be taken into account for the more complex gyrification patterns often found in cetaceans.

After accounting for multiple contours, Eqs. 2 and 3 will compute the volume and area of any enclosed contour, with reasonable precision for actual cortical surfaces. But if we need a full mesh reconstruction of said surfaces, we need to apply the Stitcher or some equivalent method. In any case, we can directly compare the values of V and A obtained from the AA with those derived from the Stitcher reconstruction.

Effects of Different Resolutions

By varying the gap between slices, we also vary the amount of information provided to reconstruct the brain. On one extreme, by having only the first and last slices, the mesh would yield a reconstruction resembling a truncated cone. On the other extreme, we would have mesh whose triangles would be so tiny we could consider it a continuous and smooth surface.

The minimum gap value to get proper cortex lies between these extremes in this spectrum of resolutions. To approximate it, we’ll use the condition:

$$\begin \begin \text _} \simeq T \end \end$$

(6)

The rationale behind this are spatial frequency of the image and practical aspects of measurements. Firstly, for any given brain MRI, we expect for structures come in and out from the screen at a certain rate as the slices are looked through. This is the spatial frequency of structures and for the folded GM with thickness T, a sulci will be defined/identified with \(2 \times T\) (to account for GM walls on both sides). In order to obey the Nyquist-Shannon Theorem (Shannon, 1949), we take half of \(2 \times T\), which is just T, to properly represent the image signal.

Lastly, it is important to have an initial guess on the gap before starting the segmentation. Therefore, by making the gap equals the thickness, one may simply use a ruler on any image processing software and directly measure the gap estimation.