Tissue segmentation‐based electron density mapping for MR‐only radiotherapy treatment planning of brain using conventional T1‐weighted MR images

Abstract Purpose Magnetic resonance imaging (MRI) is the primary modality for targeting brain tumors in radiotherapy treatment planning (RTP). MRI is not directly used for dose calculation since image voxel intensities of MRI are not associated with EDs of tissues as those of computed tomography (CT). The purpose of the present study is to develop and evaluate a tissue segmentation‐based method to generate a synthetic‐CT (sCT) by mapping EDs to corresponding tissues using only T1‐weighted MR images for MR‐only RTP. Methods Air regions were contoured in several slices. Then, air, bone, brain, cerebrospinal fluid (CSF), and other soft tissues were automatically segmented with an in‐house algorithm based on edge detection and anatomical information and relative intensity distribution. The intensities of voxels in each segmented tissue were mapped into their CT number range to generate a sCT. Twenty‐five stereotactic radiosurgery and stereotactic ablative radiotherapy patients’ T1‐weighted MRI and coregistered CT images from two centers were retrospectively evaluated. The CT was used as ground truth. Distances between bone contours of the external skull of sCT and CT were measured. The mean error (ME) and mean absolute error (MAE) of electron density represented by standardized CT number was calculated in HU. Results The average distance between the contour of the external skull in sCT and the contour in coregistered CT is 1.0 ± 0.2 mm (mean ± 1SD). The ME and MAE differences for air, soft tissue and whole body voxels within external body contours are −4 HU/24 HU, 2 HU/26 HU, and −2 HU/125 HU, respectively. Conclusions A MR‐sCT generation technique was developed based on tissue segmentation and voxel‐based tissue ED mapping. The generated sCT is comparable to real CT in terms of anatomical position of tissues and similarity to the ED assignment. This method provides a feasible method to generate sCT for MR‐only radiotherapy treatment planning.


| INTRODUCTION
Magnetic resonance imaging (MRI) is the modality of choice for defining brain tumor volume in precision radiotherapy techniques such as stereotactic radiosurgery (SRS) and stereotactic ablative radiotherapy (SABR) because MRI provides high soft-tissue contrast, functional information, and high resolution, which is superior to what can be provided by a planning CT. However, MRI alone is not sufficient in radiotherapy treatment planning (RTP) due to its lack of electron density (ED) information in its image voxel intensity values for radiation dose calculation. Conventionally, the radiation dose has to be calculated using CT images which are registered with MR images. This image registration process will introduce errors in defined tumor volumes. 1 Although these errors are usually small due to the rigid structure of the skull, image registration errors can be significant in some cases when MR and CT scans have differences in positions or setup. 2 This may lead to geometrical miss of target volumes as the target volumes are usually small and increase in dose to nearby critical organs. This, in turn, would compromise the effectiveness of treatment and the patient's quality of life. MR-only treatment planning will eliminate image registration error, minimize patient setup error, and reduce unnecessary radiation to patients from multiple CT scans.
One of the major benefits of MR-only RTP is improved workflow and reduced burden on the patient in terms of multiple visits/scans. Also, MR-only RTP is what will be needed in the future for new technologies such as the MR-Linac.
One solution to the problem is to derive ED information from MR images to generate synthetic CT (sCT) images. It is not a trivial voxel intensity standardization algorithms such as those based on tissue segmentation were developed. 3 In literature, atlas-based deformable image registration, MR bone imaging-based tissue classification or segmentation, and machine learning-based convolutional neural network (CNN) are the main approaches used to generate MR-based sCT. [3][4][5][6][7][8][9][10][11][12][13][14][15][16] Atlas-based approaches usually applied deformable registration algorithms to a pair of population based coregistered MR/CT images from one or multiple atlases to match a new patient's MR images.
The voxels of the patient's MR images were assigned ED values by warping electron densities from the deformed CT images. 3,5,6 However, the uncertainty of image registration increases with the variation of a patient's anatomy and pathology with respect to the atlases. Some patients were excluded from these studies due to significant anatomic or pathologic variations. The deformable registration error can be as much as 3 mm reported by a study evaluating some deformable registration methods. 4 Several approaches utilized a specialized MR sequence such as ultra-short echo time (UTE) to differentiate air and bone. 7-11 These sequences typically have a low signal to noise ratio and tissue contrast, bone, and other soft tissues such as fat and muscle cannot be clearly differentiated. Therefore, typical bone scan sequences in these studies included multiple UTE, T1, and T2 sequences. Some postscan image processing methods were developed to process multiple MR image series obtained from the bone scan sequences to classify image voxels as bone, air, and other soft tissues. Then, either bulk or voxel-based ED assignment was used to generate sCTs. [9][10][11][12][13] It is a costly approach. Errors due to MR artifacts such as chemical shift and motion in any of these MR sequences can be inherited by the generated sCT.
Deep learning-based CNN techniques used a trained convolutional kernel to predict the ED of each voxel of a patient's MR images to generate sCT. 15,16 The kernel consists of many layers of small image filters associated with many parameters trained using coregistered MR and CT image datasets. These methods used only one set of T1-weighted MR images and did not need deformable registration. High-resolution MR images were required since small sizes of filters were used in the kernel to avoid blurring. Therefore, these methods require computers with high computing power such as one with GPUs and large random access memory (RAM) for training kernels and generating sCTs. Large number of high-quality training data are essential for this approach as any error introduced into kernel can produce inherited errors into generated sCTs.
In our previous study, 17 a semiautomated method was developed which can segment cortical bone from T1-weighted MR and generate "CT-like" MR for radiation treatment image verification by identifying air region (sinuses and air way) through several slices of manual contours. The generated MR-based digitally reconstructed radiographs (DRR) have been verified for geometric accuracy of segmented bone. However, to generate sCT for dose calculation, electron densities of soft tissues have to be accurately assigned, especially for those soft tissues such as brain that may have MR intensities overlapping with those of trabecular bone and fat. Contrast enhanced T1weighted MR images clinically used have high contrast (high intensity value) where a tumor is located resulting in contrast-related artifacts. Cerebrospinal fluid usually has low MR intensity values which may partially overlap with those of bones.
In MR neuroimaging, methods have been developed that can automatically segment brain tissues using conventional MR sequences such as T1-weighted images. 18,19 These methods can segment brain tissue without using atlas and deformable registration.
Moreover, these methods are not very sensitive to MR voxel intensity variation due to anatomical variation or scan parameters since the methods do not rely on the thresholding of standardized voxel intensity values of MR images obtained from a set of training data.
The MR voxel intensity range of brain for an individual patient was found by searching representative voxels of brain tissue with certain anatomic or textural properties in the patient's MR images.
In this study, a voxel-based ED assignment method has been developed by segmenting tissues: air, brain, CSF, cortical bone, trabecular bone, eyes, muscles and fat from a patient's MR images. A sCT was generated by assigning voxel-based electron density values to corresponding types of tissue's voxels.

2.A | Tissue segmentation
This process was used to separate tissues with similar intensities in T1-weighted MR images and a large difference in electron density such as cortical bone, trabecular bone, brain, CSF, eyes, muscle and fat. Tissues were automatically segmented by utilizing information such as high soft-tissue contrast, edge detection, anatomical location, shapes, and statistical features of the tissue such as mean (μ), standard deviation (σ), and uniformity ( Fig. 1).
The uniformity of a tissue U tissue within a defined window such as 1 cm × 1 cm is calculated as where voxel p i ∈ voxels within the sample window of a tissue.

2.A.1 | Air mask generation
As described in our previous study, 17 separating air from cortical bone was achieved by manually contouring air regions enclosing sinuses and airway in less than 12 slices. The entire air regions in the head then were obtained by interpolating these contoured slices to slices in between. This is the only manual step in this process.
Afterward, an in-house algorithm developed with IDL8.7 (ITT Visual Information Solutions) calculated the statistic mean (μ Air ) and standard deviation (σAir) of MR intensity values of the air, then generated an air mask by subtracting soft tissues enclosed in the air regions. Subsequently, it will generate other tissue segmentation as follows.

2.A.2 | Automated brain and CSF segmentation
Determination of statistical intensity range of brain First, search a voxel with the largest intensity uniformity calculated using eq. (1) within a defined 1cm x 1cm window centered on the voxel in the middle coronal slice. The mean intensity of voxels within the window centered on the voxel is the reference intensity of brain (wmIntensity). The statistic brain tissue intensity range is from 0:53 Â wmIntensity; 1:35 Â wmIntensity ½ . 15 Edge detection MR images were processed by thresholding using intensity range found in Section "Determination of statistical intensity range of brain" to separate the skull from brain tissue. To enhance edge detection, all voxel intensities out of this range were set to 0. Then, Canny edge detection algorithm 20 was applied to each slice in axial and sagittal directions. The largest external contours of these edges within the skull for each slice were identified for brain segmentation.
The broken edges were enhanced by searching the external path of the edges and closing any gaps more than one voxel distance between two closest voxels of the edge by adding missing voxels into the edges along the shortest path between these edge voxels.
All the edge voxels were set to maximum intensity values to prevent leakage when using region growing algorithm.

Initial brain segmentation
Firstly, brain mainly consisting of white matter with high uniformity was segmented by applying region growing algorithm from the reference voxel of white matter in three dimensions with intensity range 0:7 Â wmIntensity; x007E; 1:2 Â wmIntensity ½ .

Determine statistical intensity range of CSF
Cerebrospinal fluid reference voxels were searched within an initial segmented brain in a window size of 5 mm × 5 mm × 5 mm with a maximum uniformity and intensity less than 0.53 times of wmIntensity. The mean (μ CSF ) and standard deviation (σ CSF ) of the voxels within the window were calculated.

Skull Mask
Due to surgery on the skull, segmented brain volume may spread into the spongy bone of the skull or soft tissue outside of the skull with a region growing algorithm.
Reference voxels of the skull were obtained by thresholding with an intensity range 0; μ AIR þ 3 Â σ AIR ½ from a region within 1 cm from external body contour. A skull mask was generated by region growing algorithm from the skull reference voxels with a range 0; μ CSF þ 2 Â σ CSF ½ in regions external to the segmented brain and CSF, and the edges of the brain identified in Section "Edge detection". All the soft tissues voxels external to the skull relative to that brainin distance were added into the skull mask.

CSF segmentation
CSF voxels were segmented using region growing algorithm starting from CSF reference voxels which were found in Section "Determine statistical intensity range of CSF," in brain, around brain, and spinal cord with a range from μ CSF À 2 Â σ CSF ; μ CSF þ 2 Â σ CSF ½ constrained by the skull mask.

Final brain segmentation
In this final step, the brain voxels (e.g., gray matter) within groves among initial segmented brain in Section "Initial brain segmentation" were segmented when these voxels were connected with the segmented brain tissue with intensity range found in Section "Determination of statistical intensity range of brain" and did not fall within the skull mask.

2.A.3 | Automated eye segmentation
Eyes have well-defined boundaries in T1-weighted MR images and their MR intensities are lower than those of brain. Left and right eyes are located at the anterior hemisphere relative to the brain and are symmetrical to the middle line of brain. Therefore, the eye segmentation was designed using these intensity and anatomical information. First, a brain slice with the maximum anterior posterior width was identified. In an image box which is 5 cm superior and 7 cm inferior to that brain slice at the anterior hemisphere of the head, connected regions with voxel intensity values less than 0.7 times of typical brain intensity wmIntensity were separated. From these regions, any region which was connected to the cross sagittal midline of brain, or had a too small volume (less than half of a typical eye size: 4 cm in diameter), or too large volume (more than two times of a typical eye's size) or not a spherical shape (defined as when the maximum difference of diameters in three dimensions was more than 1 cm), were excluded from eye regions. Regions which were satisfied all above criteria were identified as left or right eye.

2.A.4 | Automated bone segmentation
As described in our previous study, 17  Trabecular bones such as those in spinal bone may not be entirely enclosed by a cortical bone mask due to "thinness" of cortical bone, thickness of image slices, and partial volume effect. These trabecular bone volumes were recovered by locally searching voxels in a small neighborhood window near cortical bone voxels with intensities falling within the interval μ spongy À 2 Â σ spongy ; μ spongy þ 2 Â σ spongy Â Ã .

2.A.5 | Other soft tissue segmentation
Other soft tissues, excluding brain, CSF, bone, and air cavities mainly consist of fat and muscle. Fat usually has high intensity value in T1weighted MR images, but relatively low electron density compared to that of muscle. The minimum (min soft ), mean (μ soft ), and standard deviation (σ soft ) of soft tissues were calculated excluding bone and air. To cope with intensity variation due to the nonuniform magnetic field across the volume, μ soft and σ soft were calculated for each axial slice by including all soft tissue voxels from two adjacent slices to the slice.
Muscle intensity was assumed to fall within the an intensity range min soft ; μ soft ½ . Fat was segmented with an intensity range . Contrast artifacts were obtained by thresholding voxels with an intensity range larger than μ soft þ 2 Â σ soft .  Table 1.

2.C | Validation materials
To validate the accuracy of the generated sCT images in the head,

2.D | Evaluation methods
The generated sCT was compared with coregistered CT to evaluate its geometric similarity and accuracy of ED assignment by using the following criteria:

2.D.1 | Distances between bone contours of sCT and CT
On three axial, coronal, and sagittal slices with the largest width of external bone contours, for each voxel on the external bone contours of sCT, automatically searching for a voxel on the external bone contour of CT with a shortest distance. The shortest distance is regarded as one measurement of distance between voxels on bone contours of sCT and CT as shown in Fig. 3. The average of the shortest distances of all voxels on the evaluated external bone contours are recorded as the average distance difference of bone contours between sCT and CT.

2.D.2 | Dice similarity coefficient of bone and brain
The bone volumes were generated using a thresholding 200 HU respectively for both CT and sCT. The brain segmentation was evaluatedusing manual segmentation by a radiation oncologist for five patients in MR images. The Dice Similarity Coefficient (DSC) is calculated as.
The mean error (ME) of electron density represented by standardized CT number -HU and the mean absolute error (MAE) of electron density error for a given segmented tissue voxels (N) was calculated using following equations.
The volumes of soft tissues came from MR tissue segmentations of sCTs within the external body contours of corresponding CTs. The external contours were obtained with a thresholding of -400 HU.
Due to significant differences in air cavities because of different patient setups between MR and CT scans, the air volume considered only air voxels within the overlap of sCT and CT (less than −800 HU). ME and MAE of whole body voxels were calculated from the voxels within the overlap of external head contours of sCT and CT. The HU differences were calculated in corresponding voxels in sCT and CT.
3 | RESULTS Figure 2 shows the sCT and CT in three views. sCT is visually similar to the CT. The main differences appear close to the edges of brain and CSF. The soft-tissue contrast was preserved in sCT. Brain, CSF, brain ventricle, and brain stem can be clearly seen.  patient near the table top due to different shape of table tops (curved and flat) used in acquiring MR and CT scans. Significant differences can also be seen at mastoid sinuses which consist of small air bubbles separated by thin lamina bones.
The dice similarity coefficients for bone and brain (including CSF) are 0.72 (σ = 0.04) and 0.86(σ = 0.02), respectively. The ME and MAE for difference in soft tissues are shown in Table 2. The geometrical accuracy of generated sCT was evaluated by calculating average distances between external contours of the skulls of sCT and CT. The average difference was about1 mm. The largest difference was in the longitudinal direction. This may be due to the slice thickness of MR images which was 2.5 mm or 1.5 mm and was thicker than that 1.0 mm of CT. Different scan position was another factor. As shown in Fig. 3, the mismatching of skull contours of sCT and CT was mainly due to small differences in the positional angles between the neck and the base of the skull between MR and CT scans which cannot be corrected with rigid image registration.
As shown in Fig. 4, large ED errors were present in mastoid sinuses which are complex structures consisting of small air bubbles with thin bony separators. They were not always distinguishable from surrounding skull bone. Large ED assignment errors in these regions were also observed in other atlas, UTE, and CNN-based techniques. 12,[14][15][16] Due to the small volume of mastoid sinuses and the average effect of mixture of bone and air components, the impact of the ED errors on the overall accuracy of sCT is minimal.
Errors also occurred near oral cavities because of different patient setup and at the neck where it is close to the edges of field of view (FOV) because of MR artifacts.
In our method, voxel-based ED mapping was performed by assigning statistical ranges of CT numbers to voxels of individual tissues to minimize ED errors for each type of tissue so that the errors of overall ED assignment was minimized. Studies show that sCT generated using bulk ED assignment or voxel-based ED assignment all achieved high overall dosimetric accuracy (<5%) compared with CT as ground truth. [22][23][24][25][26][27] The large dosimetric errors were found where targets and OARs were located near air cavities or bones when bulk ED assignment was used. 10,26 The ED of cortical bone can be as high as 2000 HU, and so the ED errors can be more than 1000 HU when 800 HU was used for bulk ED assignment of bone. Therefore, this method used voxel-based approach to map EDs for cortical bone and trabecular bone. It achieved the average MAE of bone 244 HU which is significantly better than 422 HU reported in Ref. [9]. which used UTE sequence and bulk density assignment. It is worse than 130 HU of a study 7 using multiple atlases. However  Some patients even have streak dental artifacts on the CT.
One main advantage of our method compared to other methods is that it is robust with respect to patient-related variations due to differences in anatomy, pathological conditions, and aging process.
We Another main advantage of our method compared to other methods is that it is robust to tissue intensity variation due to differences in scan protocols, scan parameters, and/or scanner type. Our method utilized prior knowledge of each type of tissues such as its location, relative position, tissue contrast, edges, shape and size to automatically search for the representative tissue voxels within each type of tissue region to learn the intensity variation of different tissues for each individual patient. Therefore, our method is not susceptible to scan-related intensity variation acquired with heterogeneous scan protocols and parameters, and mitigates the impact of MR intensity variation due to nonuniformity of magnetic field and other image artifacts. We have tested the robustness of our method using validation MR images from two cancer programs which were acquired with different MR scan protocols, scan parameters, and scanners. Other atlas and deep learning-based methods reduce the impact of MR scan-related variation by strictly using same scan protocol, scan parameters, and high resolution to acquire and/or select "good" training data. The approach  can be extended to automatically generate contours of organs of risk in the head which is an advantage for MR-only RTP. The work flow with our method is quite simple to implement and more economic than the current approaches.
In future work, sCTs generated by this method will be evaluated dosimetrically using SRS and SABR patient data from multiple centers. T1-weighted MR images of these patients are acquired using different parameters, scanners, and protocols. The feasibility of using sCT to substitute CT for patient setup verification by registering with cone beam CT will be evaluated.

| CONCLUSION
A new MR electron density mapping technique was developed based on tissue segmentation information from only one set of clinical T1weighted images with contrast. The generated synthetic-CT is comparable to that of CT in terms of anatomical position of tissues and similarity of ED assignment. This method is a practical method for MR-only radiotherapy. It applies to a variety of patient anatomy and is robust to MR intensity variation.

ACKNOWLEDGMENTS
This work was supported by a project grant from Medbuy Corporationa not-for-profit organization.

CONF LICT OF I NTEREST
No conflicts of interest.