Medical Image Processing for Fully Integrated Subject Specific Whole Brain Mesh Generation

Currently, anatomically consistent segmentation of vascular trees acquired with magnetic resonance imaging requires the use of multiple image processing steps, which, in turn, depend on manual intervention. In effect, segmentation of vascular trees from medical images is time consuming and error prone due to the tortuous geometry and weak signal in small blood vessels. To overcome errors and accelerate the image processing time, we introduce an automatic image processing pipeline for constructing subject specific computational meshes for entire cerebral vasculature, including segmentation of ancillary structures; the grey and white matter, cerebrospinal fluid space, skull, and scalp. To demonstrate the validity of the new pipeline, we segmented the entire intracranial compartment with special attention of the angioarchitecture from magnetic resonance imaging acquired for two healthy volunteers. The raw images were processed through our pipeline for automatic segmentation and mesh generation. Due to partial volume effect and finite resolution, the computational meshes intersect with each other at respective interfaces. To eliminate anatomically inconsistent overlap, we utilized morphological operations to separate the structures with a physiologically sound gap spaces. The resulting meshes exhibit anatomically correct spatial extent and relative positions without intersections. For validation, we computed critical biometrics of the angioarchitecture, the cortical surfaces, ventricular system, and cerebrospinal fluid (CSF) spaces and compared against literature OPEN ACCESS Technologies 2015, 3 127 values. Volumina and surface areas of the computational mesh were found to be in physiological ranges. In conclusion, we present an automatic image processing pipeline to automate the segmentation of the main intracranial compartments including a subject-specific vascular trees. These computational meshes can be used in 3D immersive visualization for diagnosis, surgery planning with haptics control in virtual reality. Subject-specific computational meshes are also a prerequisite for computer simulations of cerebral hemodynamics and the effects of traumatic brain injury.


Introduction
Medical imaging is widely used for visualizing the interior of the body non-invasively for clinical diagnosis [1][2][3][4]. The two most common imaging techniques are computed tomography (CT) and magnetic resonance imaging (MRI). CT uses ionizing radiation and measures the radiation attenuation which is proportional to the image intensity. Dense structures, such as bones, have high attenuation thus high intensity. Soft structures, such as the brain, have low attenuation and thus low intensity. MRI uses radio frequencies to resonate hydrogen atoms, and detects the difference of signal recovery under influence of designed pulse sequence. It shows contrast between soft tissues. The MR imaging protocol can be customized to adjust image intensity based on desired features, such as flow phenomenon. Both imaging techniques output gray scale intensity maps for distinguishing diseased and normal tissue.
In addition to clinical diagnosis, medical images are also used for constructing meshes for computational modeling in simulation studies, including hemodynamic studies [5][6][7][8], traumatic brain injury [9][10][11], and cerebrospinal fluid flow studies [12][13][14][15][16][17][18]. Simulations with anatomical meshes and computational modeling provide results for experiments that are not easily performed, such as in vivo human experiments. It also provides results within a short amount of time compared to hands on experiments. There also exists mesh-free techniques for simulation [19].
Anatomical meshes are a necessary source for surgical planning with haptics control in virtual reality [20][21][22][23]. By assigning different physical properties to individual mesh structures, one can emulate the impact upon contact to better understand the spatial and physical relationships of the anatomical structures. This allows surgeons to simulate and practice surgeries to better assess and adjust the procedures to reduce the risks involved in invasive surgery.
The meshes also allow better visualization of the anatomy using cutting edge immersive 3D displays [21,22,24]. Current visualization of medical images are confined to 2D display, which shows the image slices sequentially. Segmentation algorithms are used to extract certain features for surface and volumetric rendering. However, the interaction and display is still presented in a planar screen. By combining volume rendering of original medical images with displaying volumetric meshes in an immersive environment, one can explore and interact with the anatomical structures in a global sense to better understand and diagnose the diseased situation, such as intracranial stenosis.
Nevertheless, finite resolution of the medical images lead to uncertainty in identifying the structures for mesh construction. Intensity threshold algorithms may encounter difficulties in identifying voxels around interfaces between structures especially in imaging studies of the brain. Several groups use manual input on top of image processing to better identify the interfaces. Pons proposed Delaunay-based technique for generating watertight surface and volume meshes [25]. The Cerefy atlas [26][27][28] exhibits a detailed whole brain segmentation for one subject using various image processing algorithms with manual modification. Adams [29] uses image processing with manual modification to reconstruct the cerebrospinal fluid (CSF) space. Manual segmentation requires anatomical knowledge with extensive man hours (>20 h) and is highly operator dependent. To date, there exists no fully automatic image segmentation for the whole brain.
In this paper, we propose an automatic image processing pipeline that takes multiple MRI studies and generates computational meshes for major compartments of the brain, composed of grey and white matter, cerebrospinal fluid space, skull and scalp, and arteries and veins. We first describe the imaging protocols that were used to acquire the images. Then, discuss various image processing algorithms that were used to segment the different structures. In results, we demonstrate our proposed pipeline with image data for two volunteers. Finally, we analyze the mesh quality and limitations of the algorithms and discuss possible improvements.

Image Acquisition
Two healthy human subjects with no known cerebral vascular disease were recruited and underwent MR imaging studies on a General Electric 3T MR750 scanner using a 32 channel phased array coil T1 provides contrast for the grey and white matter, skull, and scalp. T2 provided the cerebrospinal fluid space. The MRA and MRV captured major branches of the cerebral arterial and venous systems. This allows us to reconstruct the arterial and venous networks, grey and white matter surface, skull and scalp, and cerebrospinal fluid space for each subject. Automatic rigid coregistration of the MRA, MRV, T1, and T2 using one plus one evolutionary optimizer and Mattes mutual information metrics was performed to compensate for the different resolution settings and motion artifact. The image processing pipeline is illustrated in Figure 1 and discussed in details below.

Figure 1.
Proposed image processing pipeline for automatic reconstruction of subject-specific computational meshes. All medical images were convolved with a Gaussian operator for removal of inhomogeneity. Magnetic resonance angiography and venography is filtered for removal of non-vascular tissues and enhancement of vascular tissues. The filtered image is processed with marching cubes algorithm for capturing Cartesian coordinates of the vessel wall and maximal inscribed sphere is used to determine the centerline and diameter of the vascular network. T1 images are inputted into Freesurfer for segmentation of the grey and white matter as binary image matrices along with an image mask of skull and scalp. The image mask is processed with morphological operations to delineate contiguous surfaces of skull and scalp. CSF is identified using automatic thresholding and subtraction of other identified tissues.

Surface Segmentation
FreeSurfer was used to extract the brain from the T1 MR image set [30][31][32][33]. The brain extraction process removes the skull, scalp and neck from the dataset using a deformable model that is initialized in the center of the brain and then expands to the surface using a functional that weighs an expansion driving internal force and an image derived external force that slows the front at feature boundaries. The gray and white matter is classified by image thresholding and a priori geometric information.

Skull and Scalp Segmentation
Skull and scalp segmentation was performed using a procedure described by Dogdas [34]. The algorithm employs a series of morphological operations and masking using the segmented binary brain image and original T1 image in order to produce binary maps of the scalp and the skull. The algorithm segments the scalp first by using the brain mask from the T1 image with autothresholding for initial decision for distinguishing skull and scalp with a binary map. The binary map was then dilated and eroded to close the cavities. The largest connected region was then determined to be the outer border of the scalp. The skull segmentation was performed by first employing an automatic threshold for the skull for a binary map. This binary map was expanded by finding the union with a dilated brain mask. The largest component of the intersection between the expanded mask and the scalp was selected to be the base of the outer skull, which was dilated and eroded for closing the cavity. The inner skull was determined by first intersecting the original image with the outer skull followed by an automatic thresholding for a binary map. This binary map was then united with a dilated brain map for the base of the inner skull. The base was then eroded and dilated for the inner skull.

Cerebrospinal Fluid Space Segmentation
T2 images were used to delineate the cerebrospinal fluid space through the following steps. We use automatic multilevel threshold with Otsu's method [35,36] to categorize the voxels in four levels. A closing operation was performed to fill holes and maintain connectivity. Finally we isolate the largest connected component as our cerebrospinal fluid space.

Mesh Generation
The binary images for the cerebrospinal fluid space, grey and white matter, skull and scalp, were processed with the marching cubes algorithm [37,38] to find the Cartesian coordinates of the surfaces. The outputs are triangular surface meshes.

Cerebral Vasculature Segmentation
To identify and segment the vessels from the soft tissues, we utilized our in-house Hessian-based vessel filter [39][40][41][42]. Our filter outputs an image where each voxel intensity is proportional to the probability it belongs to the vascular network. At the vessel centerline, the image intensity assumes a local maximum. The centerline of the vascular network is the locus to the maxima in the filtered image, which forms a three dimensional space curve. The filtered image is then used for the vessel mesh generation procedure.

Vessel Mesh Generation
Due to the natural round shape of the vessels, we segmented the vessels based on centerlines and diameter information. The filtered images were processed to create a distance map between the vessel centerlines and the vessel wall, where a zero level set corresponds to vessel walls. The fast marching algorithm was used to detect the network connectivity by solving the Eikonal equation [43]. The connected domain serves as the initial deformable model for the levelset geodesic active contour to compute the distances of each point to the nearest vessel wall. Then, the marching cubes algorithm is applied to track the physical coordinates of the vessel walls and to construct a connected surface mesh of the vascular network. Using the maximal inscribed spheres method [44], the centerline trajectory and the corresponding vessel radius are precisely tracked. The centerline space curve representation and diameters can be used for generation of unstructured volumetric meshes [6].
Additionally, we developed an automatic parametric meshing algorithm to construct structured parametric volumetric meshes for the cerebral vasculature [45]. Our parametric meshing algorithm uses Bezier splines for approximation and smoothing of the vessel centerline. The control points of the splines are then used to define the geometrical separation planes for creating a continuous and smooth volumetric mesh. The parameterized volumetric mesh allows us to define the vessel lumen and vessel walls of the cerebral vasculature for detailed hemodynamic simulation. Figure 2 exhibits the MR imaging studies from both subjects in axial, sagittal and coronal views respectively. The maximum intensity projection is used for the MRA to better delineate the vessels. The imaging protocols inherently brighten grey and white matter in T1, cerebrospinal fluid space in T2, arteries in MRA, and veins in MRV. The scalp exhibits high intensity in both T1 and T2. The skull does not exhibit any signal in all imaging studies. Image registration was performed to compensate for the different resolution and slice settings in each imaging protocol. Figure 3 shows the registered and segmented results with red being arteries, dark blue being veins, dark grey for grey matter, light grey for white matter, light blue for cerebrospinal fluid space, white for skull, and flesh color for scalp. Figure 3 shows that the coregistration successfully aligned the different structures together.

Cortex Surface Segmentation
The T1 and T2 images in Figure 2 show clear delineation between the grey matter, white matter, and cerebrospinal fluid space. Figure 4 shows the segmented grey and white matter for the two subjects. The mesh vertices and faces are shown in Table 1. The meshes are generated with marching cubes along with binary images depicting the respective soft tissues. These images are used as input for skull and scalp segmentation.

Skull and Scalp Segmentation
The brain mask is used to first determine the scalp mesh by removing the brain on the T1 image and thresholding. The skull is determined by removing the scalp and brain on the T1 image and another thresholding. The marching cubes algorithm is utilized to construct the meshes of skull and scalp of the two subjects as shown in Figure 4. The scalp exhibits high intensity outside at the exterior of the head in T1 and T2 images. The skull does not exhibit any signal and is hard to detect in all images. Thus, the morphological operation performed between the brain and scalp is used for approximating the skull. The skull is underestimated due to lack of signal by limitation of the imaging modality. CT can be used to achieve a better recognition of the skull with the tradeoff of exposure to ionizing radiation [46]. As pointed out by Dogdas et al. [34], the morphological operations and automatic thresholding might produce holes and it can be fixed by manually adjusting the threshold value. In our data sets, we did not find holes in our reconstructed meshes.

Cerebrospinal Fluid Space Segmentation
We compute four thresholds using Otsu's method using the T2 image. The first threshold separates the background from other signals. The second threshold contains tissues with low intensity such as gray and white matter. All voxels above the second threshold are categorized as candidates for the cerebrospinal fluid space. The eyes also appears bright in T2 due to the aqueous humour in the eye. To exclude the eye, we look for the largest connected component and produced a binary image for mesh generation shown in Figure 4. Figure 5 exhibits the MRA and MRV along with the filtered images. The MRA has high intensity in the large vessels, but the small vessels have similar intensity to the grey matter. After filtering, all vessels appear brighter which facilitates automatic vessel segmentation. The MRV imaging protocol was designed to remove non-vascular tissues inherently but the small veins are still difficult to detect. After filtering, all veins have high intensity and the filtered image is used for segmentation. The MRA imaging protocol was designed to enhance intensity for high flow areas which also includes signals from the superior sagittal sinus. The contamination from the vein was removed by excluding common voxels between MRA and MRV. After applying our vessel filter, it is shown that the non-vascular tissues were suppressed and small vessels are enhanced. The filter images then served as input for the centerline and diameter extraction. The vessel mesh was reconstructed using the diameter and centerline information shown in Figure 4.  Figure 4 displays the individual meshes reconstructed from the proposed pipeline for the two subjects. Mesh metrics are evaluated in Table 1 and compared with physiological values [47][48][49][50][51][52]. The grey matter volume for both subject falls within the reported range of 710-980 mL. The white matter volume also for both subjects are also within reported values of 260-600 mL. The arterial and venous volume are slightly lower than reported values due to limited image resolution. The CSF volumes for both subjects in the lateral ventricles, third ventricles are fourth ventricles are all comparable to reported values. The scalp surface area for both subjects are approximately 10% higher than the reported value of 600 cm 2 . It demonstrates the ability of our pipeline to process image signals to reconstruct meshes for display and simulation. However, due to limited resolution and partial volume effect, the meshes intersect with each other at the boundaries which is not anatomically accurate. This inaccuracy appears mostly at the deep gyri, in which the small vessels are embedded with the cerebrospinal fluid. To compensate for the hardware limitations, we used the binary images to detect for intersections and excluded common voxels. One can also pursue Pons's method to produce watertight meshes [25].

Parametric Meshing
In addition to constructing cylinder meshes with diameter and centerline information, we also developed an automatic parametric meshing algorithm using cubic Bezier models to generate structured hexahedral meshes as shown in Figure 6. The parametric meshing algorithm takes in the diameter and length information form the vessel segmentation and groups the segments into bifurcation, bifurcation to bifurcation, and bifurcation to terminal. For each group, separation points and planes are computed to define the connection between the different groups. The parametric meshing allows us to control the meshes quality and define regions for vessel wall and vessel lumens for hemodynamic simulations. The surface meshes can be improved be adapting parameterized surfaces such as Bezier surface, B-spline surface, and NURBS [53][54][55][56].

Cerebral Hemodynamics Simulation
The computational meshes can be used for computational fluid dynamics simulations. Detailed hemodynamic simulation can be achieved by using commercial software solving the 3D Navier-Stokes equation as shown in Figure 6. The meshes are loaded into the software and assigned with boundary conditions. Our 3D simulation demonstrates velocity field, pressure gradient, and dye convection with a pulsatile boundary condition. Figure 6 also shows preliminary hemodynamic simulations for blood flow rate and blood pressure using the cylinder mesh and simplified flow principle equations. The computational time for these simulations are less than one minute. Additionally, we also performed dye convection simulation with time step = 0.1 s for 10 s duration to produce artificial dynamic angiography in Figure 6. The computational time for dye convection is dependent on the step size and in this case it took less than five minutes. The simulation results demonstrate the use of our image processing pipeline for rapid simulation to aid clinical diagnosis. In addition to cerebral hemodynamics, we also acquired images for full body CSF reconstruction along with phase contrast MRI for CSF flow measurements as shown in Figure 7. The CSF model allows us to study full body CSF dynamics [14,17]. Combined with the whole brain model, we can investigate the relationship between cerebral hemodynamics and CSF dynamics.  . CSF imaging for the entire central nervous system. T2 images are acquired for the brain, cervical-thoracic region, thoracic-lumbar region, lumbar-sacral region. Image stitching algorithms were utilized to construct a contiguous image of the entire CSF space. The image is segmented to provide a continuous computational mesh of the entire CSF space.

Conclusions
In this paper, we combined several image processing algorithms to construct an image processing pipeline for automatic generation of subject specific computational meshes. The imaging protocols used are provided by the hardware company which does not require additional modification or optimization. Currently, commercial software does not provide specific tools for individual structures, but only generic processing algorithms. Generic processing algorithms are inadequate for segmentation of vessels and smooth extraction of the cortex surface. This limits the accuracy of the reconstruction of the meshes since the algorithms are not specifically designed for different signals and characteristics of the structures. It also requires manual intervention and modification which hinders the reconstruction of subject-specific meshes. Our pipeline is designed to automatically segment individual structures based on their geometry and signal characteristics. By introducing the proposed imaging processing procedures, we hope to reduce the difficulty of reconstructing subject specific computational meshes.
her PhD thesis under the supervision of Andreas Linninger. All authors have read and approved the final manuscript.

Conflicts of Interest
The authors declare no conflict of interest.