2.1.

High-Quality Three-Dimensional Brain Mesh Generation Pipeline

2.1.1.

Brain segmentation

The brain anatomical modeling pipeline reported in this work starts from a presegmented brain. Here, we want to particularly highlight that brain segmentation is outside the scope of this paper. As noted below, there is an array of dedicated brain segmentation tools, extensively developed and validated over the past several decades. Advanced statistical and template-based segmentation methods have been investigated and rigorously implemented in these tools. It is not in our interests to develop a new segmentation method but to convert these presegmented anatomies into accurate meshes for subsequent computational modeling.

A diagram summarizing the common pathways in segmenting a neuroanatomical scan using popular neuroimaging analysis tools is shown in Fig. 1. In most cases, a tissue probability or multilabel volume is obtained for the white matter (WM), gray matter (GM), and CSF. Some segmentation tools, such as Statistical Parametric Mapping (SPM), allow the use of a matching T2-weighted MRI to improve the CSF segmentation. Utilities such as the FSL brain extraction tool and SPM can provide additional information on the scalp and skull. In this work, we primarily focus on six tissue types—WM, GM, CSF, skull, scalp, and air cavities. Additional classes of tissue (e.g., dura, vessels, fatty tissues, skin, and muscle) are also available in some segmentation outputs, which can be incorporated to our mesh generation pipeline. However, one should be aware that adding additional segmentations may result in increased node numbers and surface complexity, including disconnected surface components.

Fig. 1

Segmentation pathways from anatomical head and brain MRI scans. The common neuroimaging tools/extensions (left) and the corresponding outputs (right, shaded) are listed.

2.1.2.

Segmentation pre-preprocessing

The segmented brain or full-head volume, represented by either a multilabel 3-D mask or a set of probability maps (in floating-point 0 to 1 values), is preprocessed to ensure a layered tissue model—i.e., the WM, GM, CSF, skull, and scalp are incrementally enclosed by or share a common boundary with the later tissue layers; the scalp surface is the outermost layer and the WM is the innermost layer. In previous publications, two adjacent layers must be separated by a nonzero gap.^35,38 In this work, we have overcome this limitation by initially inserting a small gap between successive layers and then applying a postprocessing step to recover the merged boundaries. Moreover, we also consider air cavities, which can be located between two tissue surfaces, for example, inside or outside skull surfaces. In Fig. 2, the workflow to create different brain tissue boundaries is outlined.

Fig. 2

Illustration of the layered tissue model and the segmentation preprocessing workflow. Multiple air cavities are allowed. An arrow represents a thinning ($T_{-}$) or thickening ($T_{+}$) operation between two adjacent regions. Two sample pathways are indicated, shown by black and blue arrows, respectively. The circled numbers indicate the processing order. The gaps inserted between layers can be removed in the postprocessing step to recover shared boundaries (such as the CSF/GM/WM surfaces here).

To avoid intersecting triangles in the generated multilayered surface model, we first insert a small gap between adjacent brain tissue layers (only in the regions where tissues have shared boundaries). This is achieved using either a 3-D image “thickening” or “thinning” operator in the segmented volume. In the case of a thickening operator ($T_{+}$), the outer layer tissue segmentation ($P_{\text{out}}$) is modified as

Eq. (1)

$$T_{+}(P_{\text{out}};P_{\text{in}}):P_{\text{out}}\leftarrow \max [P_{\text{in}}+P_{\text{out}},D_{\epsilon }(P_{\text{in}})],$$

where $P_{\text{in}}$ and $P_{\text{out}}$ represent the inner and outer tissue probabilistic segmentations, respectively; $D$ denotes a “max-filter,” i.e., a volumetric dilation operator defined by replacing each voxel with the maximum value in a cubic neighboring region with a half-edge length of $\epsilon$, i.e.,

Eq. (2)

$$D_{\epsilon }(P):P[i,j,k]\leftarrow \max ({\{P[i+\mathrm{\Delta }i,j+\mathrm{\Delta }j,k+\mathrm{\Delta }k]\}}_{-\epsilon \le \mathrm{\Delta }i,\mathrm{\Delta }j,\mathrm{\Delta }k\le \epsilon }).$$

Similarly, a thinning operator (${\mathrm{T}}_{-}$) is defined as

Eq. (3)

$$T_{-}(P_{\mathrm{in}};P_{\mathrm{out}}):P_{\mathrm{in}}\leftarrow \min [P_{\mathrm{in}},E_{\epsilon }(P_{\text{in}}+P_{\mathrm{out}})],$$

where $E_{\epsilon }(P_{\mathrm{out}})$ is an “erosion” operator of width $\epsilon$, defined by a “min-filter” as

Eq. (4)

$$E_{\epsilon }(P):P[i,j,k]\leftarrow \min ({\{P[i+\mathrm{\Delta }i,j+\mathrm{\Delta }j,k+\mathrm{\Delta }k]\}}_{-\epsilon \le \mathrm{\Delta }i,\mathrm{\Delta }j,\mathrm{\Delta }k\le \epsilon }).$$

The $T_{-}$ operator effectively shrinks the inner layer mask at the intersecting regions with the outer layer. We note here that the above $T_{+}$ and $T_{-}$ operators work for both probabilistic and binary segmentations. We also highlight that the above operators only alter tissue segmentations in the regions where the inner tissue boundaries “merge” or intersect with the outer boundaries. Such regions only account for a small fraction of the brain tissue boundaries generated from realistic data. An optional postprocessing step is applied to “relabel” the elements in the expanded regions to recover the shared tissue boundaries (see below).

2.1.3.

Tissue surface extraction and surface-mesh processing

If the input data are given in the form of a multilabeled volume or tissue probability maps (see Fig. 1), the next step of the mesh generation is to create a triangular surface mesh for each tissue layer. This is achieved using the “$\epsilon$-sample” algorithm using the CGAL Surface Mesh Generation library.³⁹ For each extracted tissue surface, independent mesh quality and density criteria can be defined. In general, a surface mesh extracted from a probabilistic segmentation (grayscale) is smoother than the one derived from a binary segmentation. In addition, we also provide three surface smoothing algorithms, including the Laplacian, Laplacian+HC, and low-pass filters.⁴⁰ Moreover, a surface Boolean operation using on a customized Cork 3-D surface Boolean library⁴¹ is used to avoid surface intersections.

If the segmentation tool directly outputs the surface meshes of GM and WM, such as FreeSurfer, additional surface preprocessing is often required. For example, the FreeSurfer GM/WM surfaces are typically very dense. A surface simplification algorithm using the Lindstrom–Turk algorithm is applied^42,43 to decimate the surface nodes. In another example, if the pial and WM surfaces contain a separate surface for each brain hemisphere, a merge operation is applied to combine them into one closed surface. In addition, the FreeSurfer and FSL pial/WM surfaces do not cover the cerebellum and brainstem regions. To add those two anatomical regions, we first rasterize the pial/WM surfaces and then subtract those from the GW/WM probabilistic segmentations. From the subtracted probability maps, we extract the brainstem and cerebellum WM and GM surfaces and subsequently merge these meshes with the cerebral GM/WM surfaces. A diagram summarizing our surface-based brain mesh generation pipeline is shown in Fig. 3.

Fig. 3

Processing steps for a surface-based mesh generation workflow. The left side shows the steps for processing tissue probability maps and multilabel volumes, and the right side shows additional steps to incorporate precreated pial and WM surfaces. The specific algorithm used in each step is indicated in red, and dashed boxes and arrows indicate optional processing steps.

2.2.

Assessing Impact of Anatomical Models in Light Transport Modeling

The ability to generate high-quality brain mesh models alongside the in-depth understandings to the state-of-the-art voxel and MMC photon transport modeling tools allow us to investigate the impact of brain anatomical models on fNIRS data analysis. Here, we are interested in quantitatively comparing various brain anatomical models in terms of light transport modeling and quantify their differences in the optical parameters essential to fNIRS measurements. The three brain anatomical models that we are evaluating include (1) mesh-based brain models,^6,16 (2) voxel-based brain segmentations,⁴⁵ and (3) simple layered-slab brain models often found in literature.^3,46 We apply our widely adopted and cross-validated MC light transport simulators—MMC^16,17—for mesh- and layered-slab brain simulations, and Monte Carlo eXtreme (MCX)^31,47 for voxel-based brain simulations. Furthermore, a number of publications applied the DA to fNIRS modeling by utilizing an approximated scattering coefficient for CSF,^5,48 as it has very low scattering. It is beneficial to the community to understand the errors caused by such model approximation.

To ensure that the mesh model is “equivalent” to the original segmentation, we calculate the volumetric ratios, denoted as $V_{\mathrm{rel}}$, between the enclosed volume of the tissue boundary over those derived from the corresponding segmentation for a given tissue. A $V_{\mathrm{rel}}$ value close to 1 suggests excellent volume conservation. From all MC simulations, we compute important light transport parameters, such as the average partial pathlengths⁴⁹ in the GM/WM regions (${\mathrm{PPL}}_B$ in millimeters), the fraction of the partial pathlengths in the brain ($R_B$), as well as the optical fluence spatial distributions. In addition, we also compute the percentage fractions of the total energy deposition in the brain regions using both mesh and voxel-based simulations. Such parameter is strongly relevant to photobiomodulation (PBM) applications.⁵⁰

3.1.

High-Quality Tetrahedral Meshes of Human Head and Brain Models

In Fig. 4, a sample full-head mesh model is generated from an SPM segmentation with tissue priors from the Laboratory for Research in Neuroimaging.⁵² The generated mesh contains 181,026 nodes and 1,060,051 tetrahedral elements. A T1-weighted MRI scan of an average head for the 40 to 44 years old age group from the University of South Carolina (USC) Neurodevelopmental MRI database (in the following sections, we will use “USC age group” to refer to an atlas from this database; for example, the atlas used in this example is USC 40 to 44) was used as the input. The segmentation yields five tissue classes that are used in the mesh generation: WM, GM, CSF, skull, and scalp.

Fig. 4

A five-layered full head tetrahedral mesh derived from an atlas head of the USC 40 to 44 atlas. It contains (a) WM, (b) GM, (c) CSF, (d) skull, and (e) scalp layers. (f) A cross-cut view of the tetrahedral mesh is shown.

The mesh density for each tissue surface and volumetric region is fully customizable by setting the following three parameters:

Each surface of the brain tissue layer is extracted using a layer-specific $R_{\max }$ value: $R_{\max }=1.7\mathrm{mm}$ for pial and WM, 2 mm for CSF, 2.5 mm for skull, and 3.5 mm for scalp. The $V_{\max }$ value was defined as $30{\mathrm{mm}}^3$, and $q$ was set to 1.414. The probability threshold for surface extraction was set to 0.5 for each of the tissues. The entire mesh takes 53.47 s to generate using a single CPU thread. The $V_{\mathrm{rel}}$ values computed from the generated mesh layers for WM, GM, CSF, skull, and scalp are 0.9924, 0.9989, 0.9921, 1.0088, and 1.0037, respectively. The excellent match of the volumes is a strong indication that our meshing pipeline preserves the tissue shapes accurately.

The tissue surfaces shown in Fig. 4 are visually smooth with an average Joe–Liu quality metric⁵³ of 0.74, indicating excellent element shape quality. Overall, no degenerated element is found. The element volumes in Fig. 4 are well distributed (not shown) with only a small portion of relatively small elements. The significant improvement in mesh quality is demonstrated in the side-by-side comparison shown in Fig. 5. Here we compare the full head meshes created using the conventional CGAL-based direct meshing approach³² and the results from our surface-based meshing pipeline. As shown in Fig. 5(a)–5(c), there are several notable limitations from the conventional approach, namely (1) the tissue boundaries are not smooth (also evident from Fig. 2 in Ref. 32), (2) it has difficulty processing thin-layered tissue such as CSF [see Fig. 5(b)], and (3) it can produce small isolated “islands” [see Fig. 5(c)] due to the noise present in the volumetric image. In comparison, our surface-based meshing pipeline produces smooth and well-shaped surfaces with correct topological order. We also want to highlight that the CGAL-generated mesh contains over 268,000 nodes, whereas the mesh from our approach has only 158,211 nodes.

Fig. 5

Comparison between (a)–(c) conventional CGAL-based volumetric meshing and (d)–(f) the new surface-based meshing approaches. From left to right, we show sample meshes for (a), (d) GM; (b), (e) CSF; and (c), (f) WM/scalp.

In Fig. 6, we show that users can generate tetrahedral meshes of different densities, by conveniently setting $V_{\max }$ and the sizing-field ($S$) parameters. The set values are $V_{\max }=4$ and $S=2$, shown in Fig. 6(b), and $V_{\max }=2.5$ and $S=1.5$ shown in Fig. 6(c). In the case of Fig. 6(c), we also reduced $R_{\max }$ to 1.4 (in voxel length unit) for WM and GM, 1.7 for CSF, 2.0 for skull, and 2.5 for scalp.

Fig. 6

Demonstration of mesh density control. (a) The mesh contains 181,026 nodes, 1,060,051 elements, with a runtime of 53.47 s. (b) The mesh includes 499,134 nodes, 3,009,706 elements, with runtime 76.03 s, and (c) the mesh has 1,023,739 nodes, 6,220,187 elements, and 135 s runtime. The five layers of brain tissues are, in order from outer to inner: scalp (apricot), skull (light-yellow), CSF (blue), GM (gray), and WM (white).

3.2.

Hybrid Meshing Pipeline Combining Volumetric Segmentations with Tissue Surface Models

In this subsection, we demonstrate our “hybrid” meshing pathway (described in the right-half of Fig. 3). This approach combines tissue surfaces extracted from probabilistic segmentations with the pial and WM surfaces generated by dedicated neuroanatomical analysis tools, such as FreeSurfer and FSL. In the case of FreeSurfer, the raw pial and WM surfaces are very dense. A mesh simplification algorithm is performed to “downsample” the surface to the desired density. In Fig. 7, we characterize the trade-offs between mesh size and the surface error⁵⁴ at various resampling ratios of the dense FreeSurfer-generated pial surface for the USC 30 to 34 atlas.⁴³ The surface error is computed as the absolute value of the Euclidean distance (in millimeters) between each node of the downsampled surface to the closest node on the original surface.⁵⁴ As shown in Fig. 7, at resampling ratio values below 0.1 (i.e., decimating over 90% edges), both gyri and sulci show a surface error above 0.5 mm. When the resampling ratio value is increased to 0.15, the observed error at the gyri becomes minimal. Only a few regions show errors above 1 mm. A resampling ratio of 0.2 is selected in this example, giving a mean surface error below 0.2 mm.

Fig. 7

Box plots of surface errors as a function of resampling ratio (percentage of edges that are preserved) when downsampling a FreeSurfer-generated pial surface. Spatial distributions of the errors are shown as insets.

As we illustrate in Fig. 3, a number of additional step have been taken to process the FreeSurfer pial/WM surfaces. These include the merging of the left- and right-hemisphere surfaces, addition of the CSF ventricles, and the addition of cerebellum and brainstem using SPM segmentations. The final mesh, as shown in Fig. 8, contains 150,999 nodes and 917,212 elements. The mesh generation time was 158.44 s. This meshing time is significantly lower than the reported 3 to 4 h required for creating a similar mesh using mri2mesh.³⁵ While not shown here, this hybrid workflow also accepts probabilistic segmentations produced by FSL or pial/WM surfaces created by BrainSuite and similar combinations.

Fig. 8

Tetrahedral mesh generated from a hybrid meshing pathway combining FreeSurfer surfaces with SPM segmentation outputs for the USC 30 to 34 atlas. The (a) sagittal and (b) coronal views are shown. The tissue layers include scalp (apricot), skull (light-yellow), CSF (blue), GM (gray), and WM (white).

3.3.

Brain Mesh Library Generated from Public Brain Databases

To test the robustness of our meshing workflow described above, we successfully processed many of the publicly available brain segmentation datasets, including the BrainWeb database⁵⁵ and the recently published Neurodevelopmental MRI database.³⁷ For the BrainWeb atlas database, we created the corresponding mesh models directly from the available brain segmentations. For neurodevelopmental atlases, the WM and GM segmentations provided as part of the database were used. However, the CSF and skull segmentations were not directly included by the database because they are generally more difficult to create. For these missing tissues, separate segmentations for CSF and skull were created using SPM. In addition, the scalp surface was extracted from the raw MRI image using an intensity thresholding approach followed by three iterations of Laplacian+HC smoothing.⁴⁰ In Fig. 9, we show nine sample USC atlas brain meshes derived from adult and adolescent scans. In all processed MRI scans, the proposed meshing workflow worked smoothly. The average processing time is less than a minute per mesh when the voxel resolution is $1\times 1\times 1{\mathrm{mm}}^3$ and about 3 min per mesh when the resolution is $0.5\times 0.5\times 0.5{\mathrm{mm}}^3$. It is important to note that the CSF and skull segmentations in Fig. 9 have not been validated and are shown only for illustration purposes.

Fig. 9

Illustrative brain mesh examples (coronal views) produced using the Neurodevelopmental MRI Database, including (a) 16 years, (b) 17.5 years, (c) 25 to 29 years, (d) 30 to 34 years, (e) 35 to 39 years, (f) 40 to 44 years, (g) 50 to 54 years, (h) 60 to 64 years, and (i) 70 to 74 years old. The tissue layers include scalp (apricot), skull (light-yellow), CSF (blue), GM (gray), and WM (white).

3.4.

Comparing Mesh-, Voxel-, and Layer-Based Brain Models in Light Transport Simulations

Next, we demonstrate the impact of different brain anatomical models, particularly between the mesh, voxel, and layered-slab brain representations and highlight their discrepancies in optical parameters estimated from 3-D MC light transport simulations. Here we use our in-house dual-grid mesh-based Monte Carlo (DMMC) simulator¹⁷ for mesh- and layer-based MC simulations and MCX³¹ for voxel-based simulation. An MRI brain atlas (19.5 year group⁵¹) was selected for this comparison, although our methods are readily applicable to other brain models. The SPM segmentation ($166\times 209\times 223$ with $1\times 1\times 1{\mathrm{mm}}^3$ resolution) of the selected atlas and the generated tetrahedron mesh from this segmentation are used for this comparison.

In this case, a tetrahedral mesh with 442,035 nodes and 2,596,064 elements is used for the DMMC simulations. The mesh was created using our aforementioned meshing pipeline with maximum volume size $V_{\max }=30{\mathrm{mm}}^3$ and maximum Delaunay sphere radii $R_{\max }=1.2\mathrm{mm}$ for all tissue layers. In comparison, a much simplified layered-slab brain model is made of slabs of the same five tissue layers with the layer thicknesses calculated based on the mesh model: scalp: 7.25 mm, skull: 4.00 mm, CSF: 2.73 mm, and GM: 3.29 mm. WM tissue fills the remaining space. To minimize boundary effect, the layered-slab brain model has a dimension of $200\times 200\times 50{\mathrm{mm}}^3$.

For the two anatomically derived (voxel and mesh) brain simulations, an inward-pointing pencil beam source is placed at an EEG 10-5 landmark⁵⁶—“C4h”—selected using the “Mesh2EEG” toolbox.⁵⁷ Within the same coronal plane, five 1.5-mm-radius detectors are placed on either side of the source along the scalp, 8.4, 20, 25, 30, and 35 mm from the source [in geodesic distance, see Fig. 10(a)], respectively, determined based on typical fNIRS system settings.^18,58 Similarly, for the simulations with the layered-slab brain model, a pencil beam source pointing down is placed at (99.5, 100, 0) mm. A similar set of detectors are placed on one side of the source due to symmetry [see Fig. 10(b)].

Fig. 10

Comparisons of fluence distributions in an MRI brain atlas (19.5 years) using three different brain models: (a) MC fluence maps using anatomically derived mesh (computed using DMMC) and voxel (computed using MCX) brain representations, and (b) fluence maps computed using the MC and DA in a simple layered-slab brain model. Contour plots, in log-10 scale, are shown along the coronal planes with each brain tissue layer labeled and delineated by black dashed lines. In (a), the “L” and “R” markings (red) indicate the left brain and the right brain, respectively. The comparisons between the mesh and voxel tissue boundaries are shown in the inset of (a).

The 3-D MC photon simulations are performed on all three brain models, where DMMC is used for both the mesh-based and the layered-slab brain models and MCX is used for the voxel-based brain model. The output fluence distributions along the SD plane are compared in Fig. 10(a). From these MC simulations, we also report the average partial pathlengths in the brain regions (${\mathrm{PPL}}_B$), average total pathlengths ($\mathrm{TPL}$) and their percentage ratios $R_B={\mathrm{PPL}}_B/\mathrm{TPL}$ in Table 1. Moreover, we also computed the percentage fraction of the energy deposition in the GM region with respect to the total simulated energy. In addition to MC-based photon modeling, we have also applied the DA, frequently seen in the literature,⁵ to the layered-slab brain model and compared the results with those derived from the MC method in Fig. 10(b). For solving the diffusion equation, we used our in-house diffusion solver, Redbird.⁵⁹ The reduced scattering coefficient of the CSF region is set to $0.3{\mathrm{mm}}^{-1}$ as suggested in Refs. 5 and 48. The Redbird solution matches excellently with that from NIRFAST¹⁵ (not shown).

Table 1

Comparison of key optical parameters derived from MC simulations from an MRI brain atlas (19.5 years): for each detector, we compare the average photon partial pathlengths in the brain region (PPLB), total-pathlengths (TPL), and their percentage ratios (RB) derived from mesh-based (DMMC), voxel-based (MCX), and layered-slab (DMMC) brain representations at various source–detector (SD) separations.

In Fig. 10(a), the fluence contour plots produced by MCX (orange dashed-line) and DMMC (white dashed-line) agree excellently in the vicinity of the source, whereas noticeable discrepancies are observed when moving away from the source. We believe this is a combined result of (1) photon energy deposition variations due to the small disagreement between a terraced tissue boundary and the smooth surface boundary, and (2) the distinct photon reflection behaviors between a voxel- and a mesh-based surface due to the differences in the orientations of surface facets. The effect from the first cause is largely depicted by the deviation between the two solutions in the depth direction near the source. Such difference is particularly prominent near highly curved boundaries or near boundaries with high absorption/scattering contrasts, such as the CSF region beneath detectors #7, #8, and #9 in this plot. The effect from the second cause is highlighted by the worsened discrepancy when moving away from the source along the scalp layer, for example, the scalp region to the left of detector #10. Overall, the second source of error is noticeably prominent than the mismatch resulted from the first cause. This observation is further validated by disabling the refractive index mismatch calculations in both of our simulations (results not shown): the error along the scalp surface was largely removed, but the deviations in the deep-brain regions remain.

In Fig. 10(b), the DA (orange dashed-line) and MC (white dashed-line) produce well-matched fluence contour plots near the source but show significant difference in the regions distal to the source. The difference is particularly noticeable within the CSF and GM regions and above 30-mm SD separation on the scalp surface. We believe it is largely due to the error introduced by the approximated CSF reduced scattering coefficient.⁴⁸

To further quantify the differences caused by different brain representations and their impact on fNIRS brain measurements, in Table 1, we also compare several key photon parameters derived from MC simulations. Here we use the parameters derived from MMC as the reference. This is because MC solutions are typically used as gold standard, and mesh-based shape representations are known to be more accurate than voxelated domains.¹⁶ We observe that simulations using voxel-based and layered-slab brain models tend to overestimate ${\mathrm{PPL}}_B$ compared to the mesh models. For detectors #6 to #10, the voxel-based simulation gives a 21% overestimation at the shortest SD separation (8.4 mm). Such discrepancy is reduced to within $\pm 2\%$ at the largest two separations (30 and 35 mm). The $\mathrm{TPL}$ values are less susceptible to anatomical model accuracy, reporting a percentage difference between 0.1% and 1.8%. As a result, the difference in $R_B$ is largely modulated by that of ${\mathrm{PPL}}_B$, ranging between −1.5% and 21%. However, for detectors #1 to #5 located on a different brain region where the superficial layers are shallower than those under detectors #6 to #10, more pronounced overestimations of ${\mathrm{PPL}}_B$ for all SD separations, ranging from 12% to 23%, are observed, resulting in an $R_B$ percentage difference between 12% and 24%. Similarly, compared to mesh-based model, the layered-slab brain MC simulations report significant overestimation of ${\mathrm{PPL}}_B$ (36% to 166% with the highest difference at the shortest separation) and $\mathrm{TPL}$ (6.3% to 10%), resulting in significant variation in $R_B$: 151% at the shortest SD separation and 78% to 24% for four long separations. Furthermore, we have also computed the percentage fraction of the energy deposition within the GM. This fraction is 1.69% when using mesh-based brain model, and 1.42% when using a voxel-based brain model, resulting in a 16% reduction in brain energy deposition. This result could have some implications to many PBM applications.⁵⁰