8 Material and Methods
2.1 Materials
The MRI scans are obtained from Integrated Brain Research Unit (IBRU) of Taipei Veterans General Hospital. The MR images were acquired on a 1.5 Tesla GE MR scanner, using 3D-FSPGR pulse sequence (TR = 8.548 ms, TE = 1.836 ms, FOV = 260 × 260 × 1.5, TI = 400 ms, NEX = 1, flip angle = 15◦, bandwidth = 15.63 kHz, matrix size = 256 × 256 × 124, voxel size = 1.02 × 1.02 × 1.5); and diffusion weighted images (TR = 17000 ms, TE
= 70.2 ms, FOV = 260 × 260 × 2.2, NEX = 6, matrix size = 128 × 128 × 70, voxel size = 2.03 × 2.03 × 2.2) consisted of 13 gradient directions with b = 900 s/mm2, and one b = 0 s/mm2 image. The subject group includes 20 males and 44 females in total 64 normal subjects. Age range is from 21 to 62.
2.2 Methods
The purpose of our study is to create a template which can provide a stereotactic space for a subject group. Additionally, the subject images can be registered to this template space with little deformation. Hence, the goal is to minimize the variance between template and subject images. An adapted template construction and an accurate registration are needed.
We modified the procedure of MRI template construction proposed by Lee [15] to construct our MRI and DTI templates. The new strategy of registration will be described in Section 2.3.
2.2.1 Template construction
Figure 2.1, 2.2 and 2.3 illustrate our procedure for MR/DT template construction.
The first step is to align MRI to the first volume (baseline image) of DWIs for each intra-subject by using rigid registration (Figure 2.1). It is to ensure that MRI and DWIs are in the same stereotactic space. The second step is the pre-processing for MRI and DWIs as described in Section 2.2.2.
To construct a customized template, an initial reference template R0is needed. Accord-ing to Lee [15], we selected a representative subject as an initial template. The
representa-2.2 Methods 9
Figure 2.1: Intra-subject alignment
tive subject had the smallest deformation magnitude kφk during registering this subject to other subjects. To find the representative subject, all possible pairs for registering a source image i to a target image j were considered. The equation can be expressed as follows:
R0 = arg min
i {kφi,j(x)k | ∀j 6= i, x ∈ brain coordinate} (2.1) Note that, we only registered MR images to find the initial template, and directly as-signed the same representative subject of DT image as the initial template to reduce the calculus.
Following, the initial template was set as the temporary template (See Figure 2.3). each source image (of MRI) was registered to the temporary template of MRI. Each affine trans-formation will produce an affine matrix, which including four parameters: shearing, scal-ing, rotation, and translation. Then, we averaged the inverse of these four parameters from all affine matrices. Next, the average parameters were composed into an affine matrix. To avoid the temporary template will be blurred during iterations, the affine matrix was ap-plied on the representative subject instead of the temporary template. Finally, the affine result was set to the temporary template for next iteration. Repeated above step until the temporary template converges into a stable state (More detail about stopping criterion will describe in Section 2.2.4). As a consequence, an affine matrix and a correspondingly affine
10 Material and Methods
Figure 2.2: Brief flowchart of the proposed methods for MR/DT template construction
2.2 Methods 11
Figure 2.3: Detailed flowchart of the proposed methods for MR/DT template construction.
Color blue and purple show the construction of MRI and DTI respectively.
12 Material and Methods
template of MRI were generated.
We utilized the affine matrix of MRI to co-register the initial template of DTI to the same template space. Subsequently, the features were extracted from preprocessed images and affine templates of both MRI and DTI. More detail about feature extraction and the new similarity function for non-rigid registration will describe in Section 2.2.3 and Section 2.2.4.
Iteratively affine transformation supplied the temporary template to iteratively non-rigid registration. Each source image was registered to the temporary template by two step:
affine transformation and non-rigid registration. Affine transformation was a preliminary registration and it provided an affine matrix to non-rigid registration. Non-rigid registration generated a forward and an inverse deformation fields. In the similar way, the temporary template was updated by averaging the inverse deformation fields from each transforma-tion. Repeated these step until the template is convergent (More detail about stopping cri-terion will describe in Section 2.2.4). As a result, the deformation field and corresponding template (of MRI), called representative template, were created.
All subject images (of MRI) were registered to the representative template (of MRI).
Finally, we averaged all registered subject images to create the average template (of MRI).
The construction of DTI template is similar to MRI. The representative subject of DTI was co-registered to the same template space by using the deformation field of MRI. To maintain the property of DTI, the interpolation of DTI was computed in Log-Euclidean framework [2] (Section 2.2.6). Moreover, Figure 2.4 shows that after transformation, the position of voxels have been changed, but the direction have not. In order to maintain the consistency of the anatomical structure, tensor reorientation (Section 2.2.5) was applied.
Then, the representative template of DTI was completed.
All subject images (of DTI) were co-registered to the representative template (of DTI) by using the deformation field of MRI. The average template of DTI was constructed by averaging the registered subject images as well. Nevertheless, we computed the average operators in Log-Euclidean space [2] (Section 2.2.6) instead of Euclidean space to avoid tensor swelling effect.
To summarize, we constructed four template: representative template and average tem-plate of both MRI and DTI. These temtem-plate are in the same stereotactic space, and they can
2.2 Methods 13
Figure 2.4: This diagram shows the tensor reorientation. Left : Four diffusive directions are located on the original stereotactic space. Mid : Four diffusive directions are located on the new stereotactic space after transformation. Right: Four diffusive directions are located on the new stereotactic space after transformation and tensor reorientation.
represent the subject group.
2.2.2 Preprocessing procedure
Preprocessing of MRI includes brain extraction and intensity inhomogeneity correction.
We used Brain Extraction Tool (BET) version 2.1 [13, 23] and Bias Field Corrector (BFC) [22] to realize them separately.
Brain extraction is to remove non-brain tissues, which will affect brain tissues registra-tion. Besides, brain extraction benefits to increase the accuracy of intensity inhomogeneity correction. Magnetic field inhomogeneity in MRI causes intensity non-uniformity. As-sume that a received image can be separate to two parts, which are an original image and a bias field. Bias field correction predicts a bias field and removes it to obtain the original image. The tool called BET has the additional value on intensity inhomogeneity correc-tion. Figure 2.5 indicates the non-uniform intensity was corrected after brain extraccorrec-tion.
The non-uniform intensity will induce the error at image registration. For this reason, we applied BFC again for a further correction.
DWIs preprocessing includes DTI estimation and brain extraction as shown in Figure 2.6.
Fillard et al., [8] proposed a maximum likelihood strategy with Rician noise model to estimate diffusion tensor. In order to further reduce the influence of the noise, they
14 Material and Methods
Figure 2.5: This figure demonstrates each step of MRI preprocessing from left to right.
One slice of two distinct MRI images were presented in top row and bottom row. Left : Original MRI image. Middle : The result of brain extraction. Right : The result of inten-sity inhomogeneity correction. The red circle was used to show the non-uniform inteninten-sity within the same tissue.
Figure 2.6: This figure demonstrates each step of DTI preprocessing from left to right. One slice of the first volume of DWIs and DTI images were presented. Left : Original DWIs image. Middle : Estimated DTI image. Right : The result of brain extraction.
2.2 Methods 15
combined estimation and regularization that results in a maximum a posteriori estimation.
Moreover, the log-Euclidean metrics [2] (Section 2.2.6) was used to cancel the swelling effect in regularization. They developed a tool, called Medical Image Navigation and Re-search Tool by INRIA (MedINRIA), was used in our study to estimate the diffusion tensor.
It can overcome the low SNR in clinical DWIs and ensure the positive definiteness of all tensors which are estimated from noise-sensitive DWIs.
Most brain extraction methods, such as BET, are designed for MRI. Little or no studies have ever tried to extract brain based on DWIs/DTI. For an accuracy result, we extracted the brain mask from the MR image and applied the mask on DTI instead of extracting brain region directly from DWIs or DTI.
2.2.3 Feature extraction
The features extracted from both MR and DT images. The first feature is the intensity of T1-weighted image. The second one is the fractional anisotropy (FA). They are commonly used in the non-rigid registration of MRI and DTI, respectively. Figure 2.7 shows that these two features have significantly distinct information at most regions. Hence, we hope the T1+FA feature can provide complementary information for each other at non-rigid registration.
Figure 2.7: This figure displays three different views of MRI (top row) and FA (bottom row) brain image. From left to right are coronal view, sagittal view and horizontal view.
16 Material and Methods
2.2.4 Image registration
Affine transformation
Affine transformation is a preliminary registration. It aligns the global shape from source image to target image. Here, FMRIB’s Linear Image Registration Tool (FLIRT) [11, 14] was used.
Non-rigid registration
Lin [16] proposed a non-rigid registration algorithm for MRI images. This algorithm supports two important properties for a custom template construction: diffeomorphic and symmetric (or inverse-consistency). To put it differently, the algorithm guarantees the map-ping function with smooth, invertible, and one-to-one relationship. In addition, if a trans-formation is generated by registering a source image to a target image, the unique inverse of the transformation will exactly register the target image to the source image.
Averaging the inverse deformation fields is an essential step in our iteratively non-rigid registration. In the previous construction proposed by Lee [15], the inverse of a transfor-mation is a critical issue because the non-rigid registration without the inverse-consistency property. The invertible can decrease the error may induce by estimating location and interpolation. Hence, the algorithm proposed by Lin [16] was employed in this work. Fur-thermore, we modified the similarity criterion to efficiently utilize the features extracted from DTI and MRI as mentioned. The new similarity evaluation function between target image Itand source image Iscan be expressed as:
S(It, Is, Φ) = SCR(aintensityt , aintensitys , Φ) + SCR(aF At , aF As , Φ) (2.2) The similarity between target image and source image of the two features is assessed by correlation ratio SCR [21]. More detail about correlation ratio describe in Section 2.2.4.
Because we want to improve the accuracy for non-rigid registration of both MRI and DTI, the influence of MRI and DTI is set as the same. For inverse-consistency, the equation can be rewritten as shown below.
2.2 Methods 17
S(It, Is, Φ) = S(Is, It, Φ−1)
= SCR(aintensityt , aintensitys , Φ) + SCR(aF At , aF As , Φ) + SCR(aintensitys , aintensityt , Φ−1) + SCR(aF As , aF At , Φ−1)
(2.3)
The correlation ratio SCR in Equation 2.3 had been normalized to the range of zero to one. Large value indicates high image similarity. Hence, we maximized Equation 2.3 to obtain a better alignment.
Correlation ratio
Figure 2.8: Two subjects with different intensity contrast
The correlation ratio is an efficiency and accuracy measurement for image registration.
It is robust to different intensity contrast (See Figure 2.8), intensity inhomogeneity and noise [14, 21]. The correlation ratio on image similarity can be calculated by [17] (Figure 2.9): The range of features (intensity or FA) in target image was divided into NBbins. Voxels in the same bin was gathered into a set Xj, which contains Nj voxels. We collected their corresponding voxels in source image at set Xj, and calculated the variance from features of these voxels. Let N be the number of voxels in the overlapping region Ω between source image and target image. If the variance ratio of each set Xj to total volume Ω is small, the source image and target image are well aligned.
18 Material and Methods
Figure 2.9: Correlation ratio
Stopping criterion
We used an iterative strategy to create a stable template in both affine and non-rigid registration. A stopping criterion proposed by Lee [15] is used. Since the registration will not steady, if the difference between a new temporary template Ri+1 and last temporary template Ri is small enough, it is called a stable state. For this reason, the slope IDP (Equation 2.6) of intensity differences was compared to a threshold, α (Equation 2.7), for stopping iterative registration. The intensity difference IDi between Ri and Ri+1 at each voxel in brain volume Ω can be written as Equation 2.5.
IDi =X
v∈Ω
(Ri(v) − Ri+1(v))/Ω (2.5)
IDPi+1 = (IDi− IDi+1)/IDi (2.6)
If (IDPj < α) stop (2.7)
Outlier removing
Since the aim of our study is to create a stereotactic template for a subject group, the template should not be biased to any subject. Thus, if a subject image is registered to the
2.2 Methods 19
template with a larger deformation than other subjects, it would be considered as an outlier and would be removed from the template construction. We applied the strategy of outlier removing proposed by Lee [15] and replenished the outlier removing for DT image. First, the magnitude of shearing and scaling were examined when iteratively affine transforma-tion. Second, the magnitude of deform-only (not includes the deformation derived from affine transformation) along x, y, and z axes were inspected when iteratively non-rigid registration. Third, If a non-positive eigenvalue decomposed from a tensor of a subject, the subject would be set as an outlier. As a consequence, only the subjects passed these examinations can be constructed the average template.
2.2.5 DTI reorientation
Diffusion tensor contains direction of water diffusion in tissues. After affine transfor-mation, we only change position of voxels in DTI, but not direction. Likewise, non-rigid registering affects the orientation of the tissue structure. Therefore, we intend to maintain consistency of the anatomical structure by rotating each diffusion tensor according to the image transformation. Alexander et al., [1] developed a preservation of principal direction (PPD) method to estimate local rotation matrices from affine transformation or higher order transformations. The aim of PPD method is to preserve the principal direction and a plane spanned by principal direction and the second direction of a tensor after transformation.
The PPD algorithm is:
• Step 1: Find the first eigenvector e1 and the second eigenvector e2 of a diffusion tensor.
• Step 2: Compute the unit vectors n1 and n2from transformed e1 and e2, which were applied by transformation matrix, F , respectively.
• Step 3: Rotate e1 to n1with rotation matrix, R1, and rotate e2with R1 as well.
• Step 4: Find a projection, P (n2), of n2onto a plane. The plane is spanned by n1and n2, and it is perpendicular to n1. Besides, R1e2(rotate e2with R1) already lies in this plane.
20 Material and Methods
• Step 5: Calculate the rotation matrix, R2, which maps R1e2 to P (n2).
• Step 6: Set the local rotation matrix R = R2R1 and reorient a tensor D by: D0 = RDRT.
The transformation matrix F of higher order transformations, for instance a non-rigid transformation, is consists of an identity matrix I and the Jacobian matrix of displacement field Φ by: F = I + JΦ.
2.2.6 Log-Euclidean metrics
Constructing average template of DTI is calculated in Euclidean space [2]. Log-Euclidean space can avoid the defects of Log-Euclidean calculus, such as tensor swelling effect and non-positive eigenvalues. Moreover, in this space, tensors can be looked upon as vec-tors. In other words, operations of vector can be used directly on tensors in Log-Euclidean space. Changing a tensor into Log-Euclidean space is as same as to calculate a matrix log-arithm of a tensor. Additionally, we can efficiently obtain the matrix loglog-arithm of a tensor by only three steps [2]:
Where M is a diagonal matrix with eigenvalues and R is a rotation matrix with eigen-vectors. They were factored from spectral decomposition of a tensor D. In contrast, chang-ing a tensor back into Euclidean space, that is, calculatchang-ing the matrix exponential of a tensor also can be obtained efficiently by using exponential substituted for logarithm at step 2. Moreover, to reduce the computations, the 3 × 3 symmetric matrix, log(D), can be represented by a 6D vectors as follows:
log(D) '−→
2.2 Methods 21
Let N be the number of registered subject images and x be a voxel of tensor D, we obtained to construct our Average Template of DTI.
Besides, interpolation of tensors is also calculated in Log-Euclidean space. The prop-erties of tensors can be adequately maintained in Log-Euclidean space compared to in Euclidean space [2]. The interpolation with Log-Euclidean framework can be expressed as
D(x) = exp (
N
X
i=1
wilog(Di(x))) (2.10)
, where wi is the trilinear weights.
2.2.7 Evaluation methods
Overlap index (OVL)
Basser and Pajevic [4] developed a measurement of tensor overlap based on eigenvalue-eigenvector pairs. If two different tensors are perfect match, the value of OVL is close to one. Conversely, the value is close to zero. The average overlap of a volume R, such as the region of white matter or brain region, was usually employed to evaluate the accuracy of DTI registration [10, 18]. The OVL and averaged OVL are given by:
OV L =
Sanchez Castro et al., [5] utilized diffusion tensor images to evaluate MRI registration.
They co-registered DTI with the transformation of MRI, and used the standard deviation,
22 Material and Methods
called error, to check the consistency of tensor after registration. They also performed operations on tensor in Log-Euclidean framework. The error can be computed as followed:
DLOG= exp( 1
The error also can be applied to a specific volume by average each voxel in these vol-ume. In our study, the error is computed for both MRI and DTI. In fact, DLOGis exactly our Average Template of DTI. Similarly, registered subject images of MRI are averaged on Euclidean framework is the Average Template of MRI. In consequence, the error can rep-resent the discrepancy between our Average Template and registered subject images, and assists in evaluating the registration method we proposed. The MRI error can be expressed as follows: