Improving model-based functional near-infrared spectroscopy analysis using mesh-based anatomical and light-transport models

Abstract. Significance: Functional near-infrared spectroscopy (fNIRS) has become an important research tool in studying human brains. Accurate quantification of brain activities via fNIRS relies upon solving computational models that simulate the transport of photons through complex anatomy. Aim: We aim to highlight the importance of accurate anatomical modeling in the context of fNIRS and propose a robust method for creating high-quality brain/full-head tetrahedral mesh models for neuroimaging analysis. Approach: We have developed a surface-based brain meshing pipeline that can produce significantly better brain mesh models, compared to conventional meshing techniques. It can convert segmented volumetric brain scans into multilayered surfaces and tetrahedral mesh models, with typical processing times of only a few minutes and broad utilities, such as in Monte Carlo or finite-element-based photon simulations for fNIRS studies. Results: A variety of high-quality brain mesh models have been successfully generated by processing publicly available brain atlases. In addition, we compare three brain anatomical models—the voxel-based brain segmentation, tetrahedral brain mesh, and layered-slab brain model—and demonstrate noticeable discrepancies in brain partial pathlengths when using approximated brain anatomies, ranging between −1.5% to 23% with the voxelated brain and 36% to 166% with the layered-slab brain. Conclusion: The generation and utility of high-quality brain meshes can lead to more accurate brain quantification in fNIRS studies. Our open-source meshing toolboxes “Brain2Mesh” and “Iso2Mesh” are freely available at http://mcx.space/brain2mesh.


Introduction
Functional near-infrared spectroscopy (fNIRS) has played an increasingly important role in functional neuroimaging. 1 Using light in the red and near-infrared range, the hemodynamic response of the brain is probed through careful placement of sources and detectors on the scalp surface at multiple wavelengths. Relative changes in oxygenated (HbO 2 ) and deoxygenated hemoglobin (HbR) concentrations, as a result of neural activities, lead to variations in light intensities at the detectors that are used to infer the locations of these activities. The accuracy of this inference depends greatly not only upon an accurate representation of the complex human brain anatomy but also on the surrounding tissues that affect the migrations of photons from the sources to the detectors. Using the anatomical scans of the patient's head or a resembling atlas, Monte Carlo (MC) simulations are often used in conjunction with tissue optical properties to approximate the photon pathlengths that are used to reconstruct the changes in HbR and HbO 2 . While much simplified brain models, such as planar 2,3 or spherical layers, 4 as well as approximated photon propagation models, such as the diffusion approximation (DA), 5 have been widely utilized by the research community, their limitations are recognized by a number of studies. 3,6 In addition, modeling brain anatomy accurately also plays important roles in other quantitative neuroimaging modalities, such as electroencephalography (EEG) 7 and magnetoencephalography. 8 Whereas a voxelated brain representation has been dominantly used for acquiring and storing a three-dimensional (3-D) neuroanatomical volume, the terraced boundary shape in a voxelated space has difficulty in representing a smooth and curved boundary that typically delineates human tissues, resulting in a loss of accuracy, especially when modeling complex cortical surfaces with limited resolution. In addition, the uniform grid structure of the voxel space also demands a large number of cells in order to store brain anatomy without losing spatial details. This can cause prohibitive memory allocation and runtime in applications where solving sophisticated numerical models is necessary. Another approach-octree-uses nested voxel refinement near curved boundaries. This improves memory efficiency significantly but still suffers from terraced mesh boundaries. [9][10][11] Mesh-based brain/head models made of triangular surfaces or tetrahedral elements have advantages in both improved boundary accuracy and high flexibility compared to voxelated domains. Mesh-based models are not only the most common choice in computer graphics and 3-D visualization of brain structures, but they are also the primary format for finite element analysis (FEA) and image reconstructions in many neuroimaging studies. 12,13 In fNIRS, tetrahedral meshes have been reported in several studies to model light propagation and recover brain hemodynamic activation using the finite-element 14,15 or mesh-based Monte Carlo (MMC) method. 16,17 Despite its importance for quantitative fNIRS analysis, the available mesh-based brain models remain limited 6,18,19 in part due to the difficulties in generating accurate brain tetrahedral meshes.
The importance of creating high-quality brain mesh models is not limited to fNIRS. EEG relies on surface and volumetric head/brain meshes to quantitatively estimate the brain cortical activities. 20 The effects of transcranial magnetic stimulation and transcranial direct current stimulation (tDCS) can also be simulated on realistic mesh models to evaluate brain damages 13,21 or by measuring tDCS effects on major brain disorders. 22 In addition, FEA of brain tissue deformation using mesh models can assist neurosurgeons in the study of traumatic brain injuries or surgical planning. 23 On the other hand, while it has been generally agreed that mesh-based light transport models are more accurate than those using a voxel-based domain, there has not been a systematic study to investigate such difference and its impact on fNIRS. The shape differences between a voxelated and a mesh-based boundary not only influence how photons get absorbed by tissue, thus altering fluence distributions, but also greatly impact photon reflection/transmission characteristics near a tissue/air boundary. Such error could also be amplified with the presence of lowscattering/low-absorption tissue, such as cerebrospinal fluid (CSF) in the brain, and result in inaccurate estimations in fNIRS. To quantify the impact of brain anatomical models in light modeling using the state-of-the-art voxel and MMC simulators could be greatly beneficial for the community to design more efficient study protocols and instruments.
Despite the broad awareness of anatomical model differences, fNIRS studies utilizing mesh brain models are quite limited, largely due to the challenges to create brain meshes and lack of publicly available meshing tools. A large portion of these studies rely on previously created meshes 24 or using general-purpose tetrahedral mesh generators that are not optimized for meshing the brain. For example, a voxel conforming mesh generation approach, 25,26 the marching cubes algorithm 27,28 and the "Cleaver" software 29 can achieve good surface accuracy but at the cost of highly dense elements near the boundaries due to the octree-like refinement. A general-purpose 3-D mesh generation pipeline was proposed alongside with an open-source meshing software BioMesh3D. 30 This approach makes use of physics-based optimizations to obtain a high-quality multimaterial feature-preserving tetrahedral mesh models but at the expense of lengthy runtime, on the order of several hours. In another Delaunay-based meshing pipeline provided in the Computational Geometry Algorithms Library (CGAL), also supported in "Iso2Mesh" 31 in 2009 and NIRFAST in 2013, 32 the tetrahedral mesh is generated from a random point set that is iteratively refined. 33,34 This procedure is relatively fast, parallelizable, and robust. However, nonsmooth boundaries are often observed between tissue regions (example shown below). Commercially available tools, such as Mimics (Materialise, Leuven, Belgium) and ScanIP (Simpleware, Exeter, UK), offer integrated interfaces for image segmentation and mesh generation but often require manual mesh editing. A streamlined high-quality mesh-processing pipeline from neuroanatomical scans remains a challenge.
Meshing tools specifically optimized for brain mesh generation are very limited. In 2010, we reported a surface-based brain meshing approach. 16 In 2013, a similar approach was reported, 35 in which an automated meshing pipeline "mri2mesh" was reported, incorporating FreeSurfer surface models and the scalp, skull, and CSF segmentations from FMRIB Software Library (FSL) with the assistance of a surface decoupling step. 36 Although mri2mesh yielded smooth boundaries and high-quality tetrahedral elements, the reported meshing times are on the order of 3 to 4 h, in addition to the time for segmentation. Moreover, mri2mesh can only process FreeSurfer and FSL outputs, limiting its integration with other segmentation tools.
In this work, we address two challenges in model-based neuroimaging analysis. First, we report a fully automated, surface-based brain/full head 3-D meshing pipeline "Brain2Mesh"a specialized wrapper built upon our widely adopted mesh generation toolbox, Iso2Mesh, dedicated toward high-quality brain mesh generation. A major difference separating this work from the more conventional CGAL-based volumetric meshing approach 32,33 is the use of surfacebased meshing workflow. This allows it to produce brain mesh models with significantly higher quality. It is also much more flexible, processing data conveniently from multilabel (discrete) or probabilistic (continuous) segmentations and surface models of the brain. Second, this work quantitatively demonstrates that by utilizing an accurate mesh-based brain representation, one can potentially improve the accuracy in fNIRS data analysis. We analyze the errors in both fluence and photon partial pathlengths in the MC simulation outputs comparing between a layeredslab, a voxel-based, and a mesh-based brain model.
In the remainder of this paper, we first detail the brain mesh generation pipeline that we have developed to create high-quality brain mesh models. In Sec. 3, we show a variety of examples of the brain/full head meshes created using different types of neuroanatomical data inputs, including an open-source brain mesh libraries created based on the recently published Neurodevelopmental Magnetic Resonance Imaging (MRI) Database. 37 In addition, we also perform MMC simulation using the produced mesh model and compare the results with those from voxelated and layered-slab brain models. We highlight the modeling errors of using a voxelated model in both fluence and partial pathlength in comparison with a mesh-based brain.

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.

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.
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 out ) is modified as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 8 2 T þ ðP out ; P in Þ∶P out ← max½P in þ P out ; D ε ðP in Þ; (1) 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.
where P in and P 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 ε, i.e., E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 4 4 9 Similarly, a thinning operator (T − ) is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 4 0 5 where E ε ðP out Þ is an "erosion" operator of width ε, defined by a "min-filter" as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 3 6 2 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).

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 "ε-sample" algorithm using the CGAL Surface Mesh Generation library. 39 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 lowpass filters. 40 Moreover, a surface Boolean operation using on a customized Cork 3-D surface Boolean library 41 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 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).
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.

Volumetric mesh generation and postprocessing
With the above-derived combined multilayer head surface model, a full head tetrahedral mesh can be finally generated using a constrained Delaunay tetrahedralization algorithm, achieved using an open-source meshing utility TetGen. 44 The output mesh quality and density are fully controlled by a set of user-defined meshing criteria. An un-normalized radius-edge ratio (q) lower bound (see below) can be specified to control the overall quality of the tetrahedral elements. In addition, one can set an upper bound for the tetrahedral element volume globally for the entire mesh or for a particular tissue label. Spatially varying mesh density is also supported via user-defined "sizing field." 44 After tessellation, each enclosed region in a multilayered brain surface model is filled by tetrahedral elements and is assigned with a unique region label to distinguish different tissue types.
If the small gaps inserted by the aforementioned thickening and thinning operations are not desired, an optional "relabeling" step is performed to recover the originally merged tissue boundaries. To do this, we use the centroids of the tetrahedral mesh elements in each unique region to determine the innermost surface that encloses this region volumes, based on which we can retag these elements using the correct tissue labels.

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, 45 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 rel , between the enclosed volume of the tissue boundary over those derived from the corresponding segmentation for a given tissue. A V 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 49 in the GM/WM regions (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. 50

Results and Discussion
In the below sections, we first showcase the robustness and flexibility of our aforementioned brain meshing pipeline by processing a wide range of complex brain anatomical scans. Various meshing pathways of our meshing pipeline are validated by using publicly available brain atlas datasets, including the Neurodevelopmental MRI database. 37,51 In addition, we use a sample fullhead mesh generated from the MRI database and report their differences in key optical parameters by performing 3-D mesh-, voxel-and layered-domain MC transport simulations at a range of source-detector (SD) separations. All computational times were benchmarked on an Intel i7-8700K processor using a single thread.

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. 52 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.
The mesh density for each tissue surface and volumetric region is fully customizable by setting the following three parameters: • R max : The maximum radius of the Delaunay sphere 39 (in voxel unit) that bounds each of the triangles in a given surface mesh. For example, setting R max ¼ 2 requires that the circumscribed sphere of each surface triangle in the generated mesh must have a radius less than 2 (voxel length). A smaller R max value is generally needed when meshing objects with sharp features.
• V max : The maximum tetrahedral element volume 44 (in cubic voxel-length unit). For example, setting V max to 4 means no tetrahedral element in the generated mesh can exceed 4 cubic voxels in volume.
• q: The lower bound of the radius-to-edge ratio 44 (measuring mesh quality), defined by q ¼ minðfR i ∕L i g i¼1;2;: : : Þ, where R i is the radius of the circumsphere of the i'th tetrahedron and L i is the shortest edge length of that tetrahedron.
Each surface of the brain tissue layer is extracted using a layer-specific R max value: R max ¼ 1.7 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 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 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 53 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 32 and the results from our surface-based meshing pipeline. As shown in 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.
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.

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 54 at various  resampling ratios of the dense FreeSurfer-generated pial surface for the USC 30 to 34 atlas. 43 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. 54 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. 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. 35 While not shown here, this hybrid workflow also accepts probabilistic segmentations produced by FSL or pial/WM surfaces created by BrainSuite and similar combinations.

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 55 and the recently published Neurodevelopmental MRI database. 37 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. 40 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 × 1 × 1 mm 3 and about 3 min per mesh when the resolution is 0.5 × 0.5 × 0.5 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.

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 dualgrid mesh-based Monte Carlo (DMMC) simulator 17 for mesh-and layer-based MC simulations and MCX 31 for voxel-based simulation. An MRI brain atlas (19.5 year group 51 ) was selected for this comparison, although our methods are readily applicable to other brain models. The SPM segmentation (166 × 209 × 223 with 1 × 1 × 1 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 mm 3 and maximum Delaunay sphere radii R max ¼ 1.2 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 × 200 × 50 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 56 -"C4h"-selected using the "Mesh2EEG" toolbox. 57 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  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).
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)]. 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 (PPL B ), average total pathlengths (TPL) and their percentage ratios 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, 5 to the layeredslab 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. 59 The reduced scattering coefficient of the CSF region is set to 0.3 mm −1 as suggested in Refs. 5 and 48. The Redbird solution matches excellently with that from NIRFAST 15 (not shown).
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 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 (PPL B ), total-pathlengths (TPL), and their percentage ratios (R B ) derived from meshbased (DMMC), voxel-based (MCX), and layered-slab (DMMC) brain representations at various source-detector (SD) separations.

SD (mm)
Mesh-based brain model (DMMC) Voxel-based brain model (MCX) Layered-slab brain model (DMMC) 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 wellmatched 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. 48 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. 16 We observe that simulations using voxel-based and layered-slab brain models tend to overestimate 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 AE2% at the largest two separations (30 and 35 mm). The 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 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 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 PPL B (36% to 166% with the highest difference at the shortest separation) and 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. 50

Conclusion
In this work, we address the increasing needs for accurate and high-quality brain/head anatomical models that arise in fNIRS and many other neuroimaging modalities for brain function quantification, image reconstruction, multiphysics modeling, and visualization. Combined with the advance in light transport simulators, 16,17,47 our proposed brain mesh generation pipeline enables fNIRS research community to utilize more accurate anatomical representations of the human brain to improve quantification accuracy and make atlas-based as well as subject-specific fNIRS analysis more feasible. This also gives us an opportunity to systematically investigate how neuroanatomical models-ranging from the simple layered-slab brain model to voxel-based and mesh-based models-impact the estimations of optical parameters that are essential to fNIRS imaging. Specifically, we first described a fast and robust brain mesh generation algorithm and demonstrated that our MATLAB-based open-source toolbox, Brain2Mesh, can produce high-quality brain and full-head tetrahedral meshes from multilabel or probabilistic segmentations with full automation. The abilities to create tissue boundaries from grayscale probabilistic maps and incorporate detailed surface models from FreeSurfer/FSL ensure smoothness and high accuracy in representing brain tissue boundaries. The output meshes generally exhibit excellent shape quality without needing to generate excessive number of small elements, such as in many existing mesh generators. For most of the included examples, the processing time ranges between 1 and a few minutes using only a single CPU thread. This is dramatically faster than most previously published brain meshing tools. 30,35,60 Moreover, the entire meshing pipeline was developed based on our open-source meshing toolbox, Iso2Mesh, and other open-source meshing utilities such as CGAL, TetGen, and Cork. This ensures excellent accessibility of this tool to the community. In addition to developing this brain mesh generation toolkit, we have also produced a set of high-quality brain atlas mesh models, including the widely used BrainWeb, Colin27, and MNI atlases. We believe these ready-to-use brain/full head models will be valuable resources for the broad neuroimaging community.
Another important aspect of this study is that we demonstrate how tissue boundary representations, especially layered-, voxel-and mesh-based anatomical models, could impact light transport simulations in fNIRS data analysis. While the modeling error caused by voxelization in MC simulations has been previously reported, 61 we believe this is the first time such discrepancy has been quantified, particularly in the context of brain imaging, enabled by our unique access to high-quality brain meshes and highly accurate MMC simulation tools. We believe such findings could provide guidance for advancing fNIRS toward improved accuracy and broad utility. Our open-source meshing software and the brain mesh library are freely available at Ref. 62.

Disclosures
No conflicts of interest, financial or otherwise, are declared by the authors.