An Efficient Automatic Midsagittal Plane Extraction in Brain MRI

Abstract

In this paper, a fully automatic and computationally efficient midsagittal plane (MSP) extraction technique in brain magnetic resonance images (MRIs) has been proposed. Automatic detection of MSP in neuroimages can significantly aid in registration of medical images, asymmetric analysis, and alignment or tilt correction (recenter and reorientation) in brain MRIs. The parameters of MSP are estimated in two steps. In the first step, symmetric features and principal component analysis (PCA)-based technique is used to vertically align the bilateral symmetric axis of the brain. In the second step, PCA is used to achieve a set of parallel lines (principal axes) from the selected two-dimensional (2-D) elliptical slices of brain MRIs, followed by a plane fitting using orthogonal regression. The developed algorithm has been tested on 157 real T₁-weighted brain MRI datasets including 14 cases from the patients with brain tumors. The presented algorithm is compared with a state-of-the-art approach based on bilateral symmetry maximization. Experimental results revealed that the proposed algorithm is fast (<1.04 s per MRI volume) and exhibits superior performance in terms of accuracy and precision (a mean z-distance of 0.336 voxels and a mean angle difference of 0.06).

1. Introduction

Segmentation of brain in magnetic resonance images (MRIs) is one of the difficult and crucial steps of clinical diagnostic tools in medical images. The brain is the most complex organ in the human body that can be split into two approximately symmetrical hemispheres using a plane. This plane is known as the midsagittal plane (MSP) [1]. In brain symmetric/asymmetric analysis, automatic MSP extraction that is independent for symmetrical and asymmetrical brain regions is an essential brain segmentation task [2]. Enormous research reflects that the symmetrical structure of the brain deteriorates due to psychological and physical ailments in the brain [3]. Clinical experts use the symmetry of the brain to identify qualitatively asymmetric patterns that signify an ample range of pathologies, such as brain tumors [4,5], brain infections [6], metabolic disorders [7], brain injury [8], and perinatal brain lesions [9]. Similarly, the computer-aided diagnostic and image analysis systems can use the symmetry and asymmetry information as a prior knowledge to embellish the system efficiency in the analysis of altered brain anatomy [10].

Moreover, the detection of MSP is required in registration [11] of medical images as the first step for spatial normalization [12] and anatomical standardization [13] of the brain images. However, legitimate evaluation of symmetric and asymmetric patterns in brain images is possible only when the symmetry axis or the symmetry plane (MSP) is accurately aligned and appropriately oriented within the coordinate system of the MRI scanner [14]. This permits the system to adjust the possible misalignment of brain MRIs. A general phenomenon in brain MRI scanning is that many neuroimaging scanners produce tilted and distorted brain images. The tilt of the head is not always detectable, due to many reasons such as the health conditions, immobility of patients, imprecision of the data calibration systems, and the inexperience of the technicians. Consequently, the slices of the brain MRIs are no more alike within the same orientation, at either the axial or coronal level [15]. Disoriented and misaligned brain MRIs can betray visual inspection and prevalently yield erroneous clinical perception [16]. In summary, assessment of brain MRIs for any anomaly based on cross-referencing of brain hemispheres (left and right), either by a human expert or computer-based software could be affected by false geometrical representation. Consequently, it is essential to correct the tilt and realign the brain MRIs data before further analysis.

Manual misalignment correction is extremely time-consuming and laborious to perform on a huge scale. It also demands an urbane knowledge of brain anatomy. Therefore, it is neither sufficient nor efficient. Alignment or tilt correction of brain MRIs is tantamount to realigning the MSP with the center of the image matrix or image coordinate system [14]. If the MSP is computed precisely, the orientation problem of the MRI volume can be resolved. Thus, the tilt of the head volume can be assessed and adjusted. An ideal MSP can be defined as a virtual geometric plane passing through the interhemispheric fissure (IF) [17], about which the three-dimensional (3-D) anatomical structure of the brain (such as the ventricles, anterior/posterior commissures, corpus callosum, thalamus) exhibits maximum bilateral symmetry [18].

Previously, several approaches that considered the problem of computing the MSP in brain MRIs and other brain image modalities (Computed Tomography (CT), Positron Emission Tomography (PET), Single Photon Emission Computed Tomography (SPECT)) have been published. These approaches can be divided into two distinct groups, varying in their exclusive interpretation of prescribed MSP: (1) shape-based algorithms that identify the location of cerebral IF using features of the head images to estimate MSP; and (2) content-based algorithms that considered MSP as the plane which maximizes the bilateral symmetry of the brain. A comprehensive survey of all the existing MSP extraction methods can be found in a recent review [19].

Shape-based algorithms first segment the longitudinal fissure of the brain MRIs and employ it as a landmark for symmetry analysis and MSP extraction. For instance, Brummer [20] utilized Hough transform for straight line identification on each coronal slice and computed the MSP using interpolation. Guillemaud et al. [21] exploited linear snakes to find the control points on IF lines and estimate MSP plane through these lines using orthogonal regression. Volkau and Nowinski [17,22] and Kuijf et al. [23] proposed simple and accurate methods based on Kullback and Leibler’s (KL) measure. These approaches are computationally efficient and independent of internal asymmetries. However, they became unstable in the presence of strong mass effect near IF or invisibility of IF, which is common in some imaging protocols (CT, PET or SPECT).

Content-based algorithms, also known as the similarity-based methods, maximize some similarity measure between the two halves (hemispheres) of the 3-D head volume. Ardekani et al. [24] proposed an iterative local search-based algorithm that uses the cross-correlation between the voxels of either side of the estimated MSP. This method failed on images having asymmetries due to pathological effects. Liu et al. [18] computed the MSP by extracting the two-dimensional (2-D) symmetry axes on each slice using cross-correlation from an edge image, followed by plane fitting. Another technique based on the similarity between two sides of the head volume using block matching was given by Prima et al. [25]. These methods are computationally intensive due to their iterative nature and optimization scheme. Ruppert et al. [26,27] improved the efficiency of similarity and symmetric-based methods, and developed an algorithm using 3-D Sobel edge operator, downsampling, and a multiscale scheme. Although the algorithm used the sagittal orientation for MSP extraction, it can be applied to other orientations (axial, and coronal) as well. The authors tested the algorithm on limited imaging protocols and it is also sensitive to noise. The MSP extraction technique based on 3-D scale invariant feature transform (SIFT) was formulated by Wu et al. [28]. The authors determined the MSP by parallel 3-D SIFT matching and voting, followed by least median of square (LMS) regression. The paper also compared the results of the algorithm with three other MSP extraction methods [16,27,29]. The authors reported that the algorithm is sensitive to noise, blur, and asymmetry, greater than a certain threshold. Moreover, the parameter setting of the algorithm is somehow complex. A computationally simple and robust MSP extraction algorithm was presented by [30] using curve fitting. The method depends on skull stripping in brain images and the authors reported that the algorithm may fail to identify the MSP correctly if the image slices have a rotation angle of greater than 15° or unsuccessful skull stripping.

Recently, Ferrari et al. [31] devised a new MSP extraction algorithm using a sheetness measure obtained from 3-D phase congruency (PC) responses. The authors reported results on synthetic and real brain MRIs. A comparison study of three MSP extraction algorithms (symmetry-based [27], phase congruency [31], and Hessian-based [32]) is presented in [33]. In spite of the enormous variety of algorithms published on MSP extraction, there is no unanimity among the researchers about the best algorithm, due to the ambiguous longitudinal fissure lines, low-contrast brain images, mass effect, and absence of intensity standardization. Moreover, MSP extraction becomes more difficult and challenging when the brain MRIs having a pathological disorder [18,25].

In this article, we have combined the advantages of both aforementioned techniques (to some extent) and developed a new principal component analysis (PCA) and symmetric feature-based approach to automatically compute and reorient the MSP in T₁-weighted MRIs. In fact, the pathological disorder and variations, such as stroke, brain tumor, bleedings, and brain injury, only alter the local intensities and symmetries of brain MRIs. They do not affect the overall shape topological properties of the 3-D head. Furthermore, when the head volume demonstrates a low signal-to-noise ratio (SNR) and significant artifacts, the segmentation of external surfaces is easier as compared to that of internal structures.

Therefore, by considering all these observations and assuming that the head is an ellipsoid-like 3-D solid object, a PCA-based algorithm is designed for MSP extraction. PCA is a fundamental and prevailing statistical technique also known as Hotelling transform substantially used in digital image processing for data dimension reduction [34], feature pattern recognition [35], quality control [36], data decorrelation [37], data compression [38], and segmentation [39]. It is also acknowledged as a low-level digital image processing tool for tasks such as the orientation assessment and alignment of particular shape objects [40,41]. In this paper, PCA has been used for determining the rotation angle (yaw angle) of the bilateral symmetric axis of the brain. The parameters of MSP (yaw angle, roll angle, and offset) are estimated in two steps. In the first step, a coarse value of yaw angle has been estimated using PCA. The angle value is further refined using a cross-correlation method. After thresholding and elliptical area extraction, PCA is used to achieve a set of parallel lines (principal axes) from the selected 2-D slices of brain MRIs. In the second step, the roll angle and the plane offset (a perpendicular distance of MSP from the origin) have been computed by fitting a plane to these parallel lines using orthogonal regression [42]. Initial slices in brain MRIs, where no or very small brain is present (in size), show ambiguous symmetry features as compared to the slices near the center of the brain. Therefore, selected slices have been used for MSP extraction and automatically discarded the ambiguous slices based on semi-axes (major and minor of the ellipse). Similar to the work by Liu [18] who used a weighted mean due to biasing in mean by the initial slices as compared to the superior slices, the removal of ambiguous slices of brain MRIs makes this technique to perform robustly and efficiently. Finally, an affine transformation has been applied to rotate the 3-D head volume to realign (recenter) within the required coordinate system (scanner coordinate system). The proposed technique is insensitive to pathological asymmetries, acquisition noises, and bias fields.

The rest of the paper is categorized as follows: Section 2 describes the methodology of MSP extraction algorithm. Implementation of the algorithm, the description of datasets used for evaluation, and results are reported in Section 3. Section 4 discusses some limitations of the developed algorithm and concludes the proposed technique.

2. Materials and Methods

2.1. Geometry of MSP

Generally, MRI of the brain consists of 3-D volumetric data in three orientations: axial, coronal, and sagittal orientations. In the proposed algorithm, only the axial orientation was considered as an input to the algorithm. The head coordinate system is defined as the ideal head coordinate system (X_o, Y_o, Z_o) and the imaging coordinate system (X, Y, Z), as described by Liu et al. [18]. The origin of the ideal coordinate is the center of the brain with positive X_o pointing to the right. Anterior and superior directions represent positive Y_o and Z_o axes from the center of the brain, respectively. Mathematically, X_o = 0 is defined to be the MSP with respect to the ideal coordinate system. Practically, the imaging coordinate system (blue) varies from the ideal coordinate system (black) due to translations and three rotations (pitch ω, roll $\phi$, and yaw θ with respect to X_o, Y_o, and Z_o axes, respectively) of the patient head as portrayed in Figure 1. Therefore, the main objective of MSP extraction is to determine the transformation between the two planes, i.e., X_o = 0 and X = 0.

MSP in an image coordinate system can be defined as:

$$aX+bY+cZ+d=0$$

(1)

where parameters $(a,b,c)$ are not all zero and can be scaled by any non-zero scalar, vector $(a,b,c)$ is the normal vector of the MSP, and $d/\sqrt{a^2+b^2+c^2}$ is the perpendicular distance of the plane from the origin.

The intersection of MSP with each axial slice (slice cut perpendicularly to the Z_o) is always a vertical line in the ideal coordinate system. The rth axial slice is represented by a plane equation as:

$$Z=Z_r$$

(2)

The intersection of Equations (1) and (2) is the vertical line (ideally a line of bilateral symmetry of brain) on the rth slice and can be written as:

$$aX+bY+(c{Z}_r+d)=0$$

(3)

This represents a normal equation of the 2-D line in the XY plane. By comparing Equation (3) with the standard normal equation of the line, the orientation θ_r of the line (vertical line) can be computed as:

$${\theta }_r={\tan }^{-1}\left(\frac{b}{a}\right)$$

(4)

Equation (4) exhibits that the orientation of all the 2-D symmetric lines should be the same, irrespective of the position of slice i.e., Z_r. This angle corresponds to the yaw angle of the patient’s head. This angle is estimated using the PCA and cross-correlation technique described in the succeeding paragraph.

Similarly, the perpendicular distance (translation offset) p_r of the line (Equation (3)) from the point (0, 0, Z_r) can be calculated as:

$$p_r=c{Z}_r+d$$

(5)

Equation (5) demonstrates that the offset p_r of the symmetric line on the rth slice (Z = Z_r) is linearly related to slice position as a function of plane parameters c and d. This represents an overdetermined set of linear equations in p_r and Z_r. It can be solved by fitting a plane to a set of parallel lines having the orientation θ_r. We use an orthogonal regression [42] using PCA to fit the plane to these lines in 3-D Euclidian space.

To completely determine the MSP parameter (a, b, c, d), we need to assess the transformation between the two planes X_o = 0 (MSP in the ideal coordinate system) and X = 0 (MSP in the image coordinate system). The derivation of this transformation can be found in [18]. The final expression for MSP X_o = 0 can be written in terms of the imaging coordinate system as:

$$X\cos \phi \cos \theta +Y\cos \phi \sin \theta -Z\sin \phi -\left(n\cdot \Delta \right)=0$$

(6)

where $n={\left[\cos \phi \cos \theta ,\cos \phi \sin \theta ,\sin \theta \right]}^T$ is the unit normal vector of the plane, $(\cdot )$ indicates the standard dot product of the vectors, and $\Delta ={[\Delta X_o,\Delta Y_o,\Delta Z_o]}^T$ is the translation vector.

Dividing by cos φ and abs(φ) ≠ 90°, and comparing the result with Equation (1), we obtain:

$$a=\cos \theta ,b=\sin \theta ,c=-\tan \phi ,d=-\frac{n\cdot \Delta }{\cos \phi }$$

(7)

Therefore, $\phi$ is the roll angle and θ = θ_r is the yaw angle of each axial slice of the head’s imaging coordinate system.

2.2. Estimation of Yaw Angle (θ_r)

The yaw angle is estimated in two stages. A coarse value (θ₁) is estimated in the first stage using PCA after the region of interest (ROI) extraction. Then, the image is vertically aligned by θ₁. The value of the yaw angle is further refined in the second stage and a cross-correlation technique is exploited to measure θ₂. The sum of θ₁ and θ₂ completes the procedure of calculating yaw angle (θ_r).

2.2.1. Region of Interest Extraction

Due to the fact that the human brain is nearly elliptical, PCA is used, as it aligns the data in the direction of maximum variance (data spread). A reference 2-D slice I_o is selected from the 3-D volume of the brain MRIs. The selection of the slice is important. As in fact, the higher slice (before and after the half of the total slices present in the volume), the brain in the 2-D image becomes more elliptical in shape. It can give a good estimate of the angle of brain symmetric axis as compared to other slices [43]. Therefore, a 1/55th slice of the total slices present in the volume of brain MRIs is considered as a reference slice. Next, the image is binarized [44] and noise is removed by a mathematical morphological filtering operation called area opening [45] (pp. 112–114). In this filtering procedure, small connected pixels of which the area (in a number of pixels) is less than a specified threshold (η) are removed. The value of “η” can be varied depending upon the type of brain MRIs and the extent of noise. After several experiments, a value of 100 pixels has been considered a good estimate for η (area threshold) in all the brain MRIs including tumor datasets. Mathematical, area opening is represented as:

$${\mathsf{\gamma }}_{\mathsf{\eta }}={\displaystyle \underset{i}{{\displaystyle \cup }}\left\{{\mathsf{\gamma }}_{{\mathrm{B}}_i}|{{\rm B}}_i\mathrm{is}\mathrm{connected}\mathrm{and}\mathrm{card}({{\rm B}}_i)=\mathsf{\eta }\right\},}i=1,2\dots ,k.$$

(8)

where γ_η is the area opening, η is the value of threshold area of the connected pixels, and card(B) is the number of elements (cardinal number) of B.

Area opening is equivalently defined as the union of all the morphological opening (erosion followed by dilation) with the connected structuring element of which the size is equal to η. Then, a rectangular area is achieved by searching for the first and last nonzero pixels along the rows (top and bottom) and columns (left and right) of the noise-free binary image, as shown in Figure 2.

Some of the ROI pixels inside the rectangular boundary (image I₃) have value 0. To make all the ROI pixels in the image I₃ to one, the complement of noise-free binary image I₂ is multiplied (logical AND) with the rectangular boundary, such that all the pixels inside the boundary are equal to 1 (Figure 3a). Largest connected component (LCC) is chosen (Figure 3b) from the resulted image and added (logical OR) to I₃ (Figure 3c). Lastly, morphological flood-fill operation [45] (p. 208) is used to fill the holes (0-valued pixels) and achieved the binary image I₄ with a single ROI. This image is exploited as an input for PCA (Figure 3d).

2.2.2. Principal Component Analysis

After noise removal and ROI extraction steps, the method uses 2-D coordinates (column and row) of the non-zero pixels of the object’s region as data points and an array X is formed as:

$$X=\left[\begin{array}{cc}{c}_1& r_1\\ c_2& r_2\\ ⋮& ⋮\\ c_n& r_n\end{array}\right]$$

(9)

where r, c, and n represent row, column, and the total number of data points (1-valued pixels) in the ROI, respectively.

The size of X is n × 2, and rows of X are the location coordinate values of each non-zero pixel. Mean of X is a row vector m_x (mean of the elements in each column of X), and can be computed as:

$$m_x=\frac{1}{n}{\displaystyle \sum _{i=1}^n{X}_i}$$

(10)

Similarly, a covariance matrix P_x can be calculated as:

$$P_x=\frac{1}{n-1}{\displaystyle \sum _{i=1}^n(X_i-m_x)(X_i-m_x{)}^T}$$

(11)

The mean subtraction is vital for performing PCA to ensure that the first principal component represents the direction of maximum variance. The covariance matrix P_x is a real and symmetric matrix with a size of 2 × 2. Thus, finding a pair of orthonormal eigenvectors is always possible [46].

Suppose the elements of covariance matrix as:

$$P_x=\left[\begin{array}{cc}f& g\\ g& h\end{array}\right]$$

(12)

where $f=\mathrm{var}(x_{ii})$, $h=\mathrm{var}(x_{jj})$, and $g=\mathrm{cov}(x_{ij})$.

The principal components can be determined by solving an eigenvalue problem as:

$$\left(P_x-\mathsf{\lambda }I\right)e=0$$

(13)

where $\mathsf{\lambda }$, I, and e are the eigenvalue, identity matrix, and eigenvector, respectively.

If Equation (13) is to have a solution other than vector zero, then $\left(P_x-\mathsf{\lambda }I\right)$ must be a nonsingular matrix. Therefore, it leads to a characteristics equation as:

$$\det \left(P_x-\mathsf{\lambda }I\right)e=0$$

(14)

After expansion of Equation (14), it becomes a second-degree equation as:

$${\mathsf{\lambda }}^2-\mathsf{\lambda }\left(f+h\right)+\left(fh-g^2\right)=0$$

(15)

The Equation (15) can be solved using a quadratic formula as:

$${\mathsf{\lambda }}_1,{\mathsf{\lambda }}_2=\frac{tr\left(P_x\right)\pm \sqrt{tr{\left(P_x\right)}^2-4\left|P_x\right|}}{2}$$

(16)

where $tr\left(P_x\right)=f+h$ and $\left|P_x\right|=fh-g^2$.

The corresponding eigenvectors of the eigenvalues $\left({\mathsf{\lambda }}_1,{\mathsf{\lambda }}_2\right)$ can be calculated as:

$$e_j=\frac{1}{\sqrt{g^2+({\mathsf{\lambda }}_j-f)}}\left[\begin{array}{c}g\\ {\mathsf{\lambda }}_j-f\end{array}\right]$$

(17)

where $e_j\left(j=1,2\right)$, and ${\mathsf{\lambda }}_j$ denote the eigenvectors and eigenvalues of $P_x$, respectively.

The eigenvalues represent the length of the semi-axes of the ellipse and the eigenvectors specify the directions of these axes, as illustrated in Figure 4.

The rotation angle can be computed using the eigenvector corresponding to the largest eigenvalue as:

$${\theta }_1={\tan }^{-1}\left(\frac{v_y}{v_x}\right),{\theta }_1\in \left(-{90}^{°},{90}^{°}\right)$$

(18)

where $v_x,v_y$ are the horizontal and vertical components of the eigenvector associated with the largest eigenvalue, respectively.

This procedure provides a rough estimate of brain bilateral symmetric axis orientation θ₁. Next, the slice is realigned with the vertical axis of the image using the angle θ₁. Generally, it was observed that PCA aligned the brain bilateral symmetric axis with the vertical axis of the image with an error of less than 1°, as displayed in Figure 5.

The first row of Figure 5 indicates that bilateral symmetry axis of the brain is successfully aligned with the vertical axis (white vertical line) of the image by PCA, but sometimes as the second row in Figure 5 depicts, the bilateral symmetric axis of the brain is not completely aligned with the vertical axis of the image. Therefore, to ensure the accurate alignment of the bilateral symmetric axis of the brain with the vertical axis of the image, another fine alignment step is applied using a cross-correlation technique to find the angle θ₂ with the vertical axis of the aligned image F (output of the previous step).

2.2.3. Cross-Correlation

In the cross-correlation technique, the aligned image F acquired in the last step is reflected about the current vertical center line to produce a new reflected image F′. The reflected image F′ is rotated by 2θ₂ (θ₂ ∈ [−10°, 10°]) with a 0.5° interval about the center of the image, cross-correlated with the image F, and the maximum correlation score is noted. The value of θ₂ at which the cross-correlation score is maximum will be the required angle of bilateral symmetry with the vertical axis of the image. The cross-correlation is accomplished in the frequency space for greater efficiency. The detail of this method can be found in [18].

If θ₂ is zero, it means that the PCA accurately aligns the brain bilateral symmetric axis with the vertical axis of the image, otherwise θ_r (yaw angle) will be:

$${\theta }_r={\theta }_1+{\theta }_2$$

(19)

where θ₁ is the angle obtained from PCA and θ₂ is the angle yielded from the cross-correlation technique.

After applying the cross-correlation technique, the brain bilateral symmetric axis is completely aligned with the vertical axis of the image, as displayed in Figure 5 (last column). Now, the angle θ_r can be used to estimate the first two parameters of the normal of required MSP plane using Equation (7), i.e., a = cos θ, and b = sin θ.

2.3. Fitting of Plane in Three Dimensions

According to Equation (4), the angle θ_r will be same for all the 2-D axial symmetric lines on each axial slice, irrespective of the position of slice Z_r. Therefore, each image in the 3-D volume of the brain MRIs is rotated by an angle of −θ_r, so that bilateral symmetry axis of the brain becomes vertically oriented in each image. Normally, the translation offset p_r of each symmetric axis can be computed using a cross-correlation of the aligned image with its vertical reflection about the center of the image [18,25]. This technique has two complications. The first one is the existence of outlier in the translation offset due to pathology effects, image artifacts, and symmetry axis ambiguity in the superior slices of brain. The second one is that this technique is also computationally intensive as it takes approximately 10 s for each slice of size (rows, column) 512 × 512 [18].

We have introduced a fast approach to estimating the translation offset independent of these constraints. The technique is based on the straightforward and effective observation that the head in the slices of brain MRIs is shaped like an ellipse (ellipsoid in three dimensions). Moreover, a trapezoid area on the surface of the head, center upon +x and −x (from ear to ear) directions of a height (30–45°) and a width (45–70°), has the most significant geometrical features as described by [16]. Fitting of plane consists of the following three steps:

1. Elliptical Area Extraction

The main objective to extract the elliptical area is to make the offset estimation independent of pathological effects, image artifacts, symmetry axis ambiguity in the higher slices, and computationally efficient. The aligned image is binarized [44] and noise is removed by the same procedure as described in the previous Section 2.2.1. Then, a rectangular area is achieved by searching for the first and last nonzero pixels along the rows (top and bottom) and columns (left and right) of the noise-free binary image (Figure 6d). First and last nonzero columns and rows are denoted by c_1, c_2, r_1, and r₂, respectively. Accordingly, the vertices of the rectangle are labeled as: A(c_1, r₁), B(c_2, r₁), C(c_2, r₂), and D(c_1, r₂), respectively.

Now, the parameters of an ellipse can be determined as:

$$Centerofellipse:O=\left(\frac{c_1+c_2}{2},\frac{r_1+r_2}{2}\right)$$

(20)

$$Semi-majoraxis:a=\left(\frac{r_2-r_1}{2}\right)$$

(21)

$$Semi-minoraxis:b=\left(\frac{c_2-c_1}{2}\right)$$

(22)

Now, ellipse in parametric form can be written as:

$$x=\frac{c_1+c_2}{2}+b\cos \theta$$

(23)

$$y=\frac{r_1+r_2}{2}+a\sin \theta$$

(24)

Now, within I(i, j), which consists of all the pixels inside the elliptical boundary, the pixels are set to “1” based on Equation (25), as shown in Figure 6e,f, respectively.

$$I(i,j)=\{\begin{array}{cc}1& \mathrm{if}(i,j)\mathrm{lies}\mathrm{inside}\mathrm{elliptical}\mathrm{area}\\ 0& \mathrm{otherwise}\end{array}$$

(25)

2. Set of Mid-Parallel Lines Extraction

In the axial orientation, inferior and superior slices of the brain volume have ambiguous symmetry axes as compared to the slices near the center of the brain. Secondly, slices higher in the brain are almost ovals [16]. Therefore, ambiguous symmetry slices are automatically eliminated based on the ratio of the semi-axes of the ellipse, and only those slices of which this ratio is greater than 1.2 are extracted. This significantly improves the accuracy and efficiency of the algorithm. Selected sample slices with their respective elliptical area are illustrated in Figure 7.

The midline (principal axis) on each elliptical area slice is extracted using PCA as described in the preceding Section 2.2.2. An array A_r,s is formed from the location coordinates of three points (top, center, and bottom) on each midline (Figure 6f) as:

$$A_{r,s}={\left[\begin{array}{ccc}{x}_{r,s}& y_{r,s}& z_r.1_s\end{array}\right]}_{rs\times 3}r=1,2,\dots ,N.$$

(26)

where r denotes the rth brain slice, s is the number of points considered on each midline, 1_s is a column vector (s-dimensional) with all its elements equal to one, and N is equal to the total number of slices used (automatically selected).

The first column of A_r,s consists of p_r (offset) values of the symmetric axis and the last column contains the indices of the respective slice. According to Equation (5), this makes an overdetermined set of linear equations in p_r and Z_r as a function of plane parameters c and d, as displayed in Figure 8. It can be solved by fitting a plane in 3-D Euclidean space to a set of parallel lines (midlines) having the orientation θ_r.

3. Fitting of Plane Using Orthogonal Regression

Orthogonal regression is employed using PCA to fit a plane to these midlines. PCA minimizes the orthogonal distances from the data point to the fitting plane (fitting model). In the linear case, it is also known as total least squares [42]. It is appropriate when all the variables are measured with errors. In contrary to the usual regression, where the assumption is that the predictor variables are measured precisely, and only the response variables have the component of error. Singular value decomposition (svd) is used to find the principal components of the PCA. The first step is to center the location data matrix A_r,s. This can be achieved by subtracting each data point of the matrix from its column mean. The resultant matrix is labeled as B. Then, the svd of B can be expressed as:

$$B=US{V}^T$$

(27)

where U is a left-singular vector, S is a diagonal matrix of singular values, and V is a right-singular matrix.

The columns of V will be the required principal components. First two columns (first two principal components) of V define vectors that form a basis for the plane, and the third column (third principal component) is orthogonal to the first two principal components. The coefficients of the third principal component define the normal vector (n) of the MSP having an orientation of θ_r (yaw angle). The third coefficient of the normal vector n can be exploited to estimate the roll angle $\phi$ using Equation (7) as:

$$\phi =-{\tan }^{-1}\left(c\right)$$

(28)

where c is the third coefficient of the normal vector.

Similarly, the parameter d of the MSP can be calculated by taking the scalar product of the normal vector n and the mean of A_r,s.

2.4. Transformation for Tilt Correction

After MSP normal vector computation, the translation matrix and the rotation matrix can be easily determined to correct the tilt (recenter and reorientation) of 3D-volume of brain MRIs. Let R_req = yaw(θ)roll($\phi$)pitch(ω) be the required rotation matrix, which can be written as:

$$R_{req}=\left[\begin{array}{ccc}c\phi c\theta & c\theta s\omega s\phi -c\omega s\theta & c\omega c\theta s\phi +s\omega s\theta \\ s\theta c\phi & c\omega c\theta +s\theta s\omega s\theta & c\omega s\theta s\phi -s\omega c\theta \\ -s\phi & s\omega c\phi & c\omega c\phi \end{array}\right]$$

(29)

where $c\theta \equiv \cos \theta ,s\theta \equiv \sin \theta ,\mathrm{and}\mathrm{so}\mathrm{on}.$

The pitch (ω) angle will be zero for MSP, and yaw (θ) and roll ($\phi$) angles are calculated using Equations (19) and (28), respectively. The translation vector between the two coordinate systems, i.e., the center of the volume and centroid of the image grid, can be written as:

$$\Delta ={\left[\Delta X_o,\Delta Y_o,\Delta Z_o\right]}^T$$

(30)

Trilinear interpolation is used to reslice the head volume after realignment and tilt correction.

3. Results and Discussion

The presented algorithm has been implemented in MATLAB 2018a on a PC with Intel(R) Core (TM) i5-6600 CPU @ 3.30 GHz, 8 GB RAM. Total time for all algorithmic steps is 1.04 s on average. No attempt has been made on optimization of the code. The developed algorithm has been tested on 157 real T₁-weighted brain MRI datasets including 14 cases from the patients with the brain tumors. All the brain MRI datasets are publicly available. The details of the sample dataset images and parameters are given in

Table 1. Datasets used for evaluation.

Datasets	Detail of Images
NFBS [48]	There are 125 T1-weighted MRI scans, 77 females and 48 males in the 21–45 age range (average: 31) with a variety of clinical and subclinical psychiatric symptoms. The size of the individual scan is 256 × 256 × 192 and each voxel size is 1 × 1 × 1 mm3. The first two dimensions in each scan size indicate the individual image size (rows, columns) and the third dimension represents the number of images in the scan.
IBSR [49]	Eighteen volumes of T1-weighted brain MRI from all age groups from juvenile to adult are available online with ground truth. The size of the individual scan is 256 × 256 × 128 and each voxel size is 1.5 × 1.5 × 1.5 mm3. Most of the scans in this database have low-contrast images.
MNI BITE [50]	Real T1-weighted brain MRI of 14 patients with brain tumors (gliomas). We have used scans from Group 2 (pre-operative MRIs) and Group 3 (post-resection MRIs). The size of each scan in Group 2 is 394 × 466 × 378. Group 3 contains scans of different sizes and dimensions.

. The first dataset is from Neurofeedback Skull-stripped (NFBS) repository [48] containing 125 volumes of brain MRI having several clinical and subclinical psychiatric syndromes. There is no ground truth (GT) for MSP available for this database. The second database is from the Internet Brain Segmentation Repository (IBSR) [49], containing 18 volumes of T₁-weighted brain MRI with manually segmented hemispheres of the brain. The third dataset is from Montreal Neurological Institute’s Brain Images of Tumors for Evaluation (MNI BITE) database [50], which consists of 14 patients, 5 women and 9 men with a mean age of 52 years. The dataset includes 4 patients with low-grade gliomas (brain tumors) and 10 patients with high-grade gliomas (brain tumors). The mean tumor volume calculated manually by the experts was 30 cm³.

3.1. Evaluation on Real Datasets

The MSP is extracted from each dataset using the proposed algorithm and some slices perpendicular to the estimated MSP are snipped and displayed in Figure 9. The green line in each image is the intersecting line between the estimated MSP and the corresponding orthogonal slice. The first row (Figure 9a) represents the images from the NFBS database with extracted MSP by the proposed algorithm. The images of the same dataset (NFBS) are synthetically degraded by adding zero-mean Gaussian noise of several levels. Proposed algorithm breaking points can be found by incrementally adding the noise until the algorithm fails to detect the accurate MSP plane. The proposed algorithm successfully estimated the MSP at levels of noise up to SNR = −10.09 decibel (dB). The second row (Figure 9b) indicates representative resulting slices and the estimated MSP from noisy images. Images from the IBSR database are portrayed in Figure 9c. Manual delineation (GT) of brain hemispheres is available only for the IBSR database. The manual delineation boundary is superimposed on the input image with the white pixels by using a morphological gradient and binary skeletonization [51], as shown in Figure 9c. The proposed algorithm successfully extracts the MSP in all the volumes of the IBSR database and no obvious error is detected. The last two rows in Figure 9 contain the MSP-extracted results from the MNI BITE database. Group 2 (Figure 9d) comprises pre-operative MRIs and Group 3 (Figure 9e) includes MRIs acquired at different intervals of time, i.e., before and after surgery. Obvious asymmetries can be seen in the two groups of brain MRIs. The proposed algorithm extracted the symmetry axes from these volumes robustly and accurately.

3.2. Evaluation and Comparison on Synthetic Datasets

The accuracy of the proposed algorithm for extracting MSP is also evaluated by creating a set of 50 symmetrical scans of the real brain MRI from the NFBS database [48]. Each head scan is manually adjusted and perfectly aligned followed by reflecting one half of the head volume about the known MSP to form the other half. The two mirror halves are stitched together to create a symmetrical head scan with known GT MSP. The reason for creating such volumes is that perfectly symmetric head volumes with GT can be manipulated and transformed arbitrarily. In this way, we can avoid typical subjective factors of human visual inspection. Moreover, in reality, no human head scan exhibits perfect digital symmetry [52]. Therefore, GT in real brain MRIs cannot be used directly for MSP algorithm evaluation [18].

The presented algorithm results have been compared with a state-of-the-art MSP extraction method proposed by Ruppert et al. [27]. The algorithm is based on maximization of bilateral symmetry using 3-D Sobel edge operator, thresholding, downsampling, and a multiscale scheme. To improve the quantitative analysis of MSP identification, the authors also introduced a new MSP estimation error metric called average z-distance. The detail of this metric is discussed in the succeeding paragraph.

Metrics: Two metrics are measured to assess how accurately the algorithms detected the MSP. One is the angle difference (in degree) and is defined as the angle between normal vectors of GT MSP and estimated MSP. Mathematically, it can be computed using the inner product as:

$$\alpha ={\cos }^{-1}\left(\frac{<u,v>}{\Vert u\Vert \Vert v\Vert }\right)\times \frac{180}{\mathsf{\pi }}$$

(31)

where α is the angular difference, u is the normal vector of GT MSP, and v is the normal vector of calculated MSP.

When the two planes are parallel, this error metric is not sufficient and may be misleading due to translation between the estimated MSP and GT MSP. This problem is circumvented by Ruppert et al. [27] who proposed a new, simple, and fast metric known as average z-distance or z-score (in voxels), to measure MSP estimation error as a function of the distance between the two planes. This distance can be measured [27,33] by computing z coordinate from the plane equations, for the GT plane and the estimated plane, using each “x” and “y”. It can be written as:

$$zdistance=\frac{{\displaystyle \sum _{(x,y)}\left(\left|z_{coord.}(GT)-z_{coord.}(Est.)\right|\right)}}{dim(x)\times dim(y)}$$

(32)

where dim is the image dimension along x and y axes, and GT, and Est. stands for ground truth and estimated MSPs, respectively.

Both the algorithms (the proposed algorithm and Ruppert et al. algorithm) are evaluated for all the slices in the perfectly symmetrical head volume by determining their average z-distance with respect to each corresponding GT MSP. The mean z-scores obtained for both algorithms on 50 perfectly symmetric datasets are shown in

Table 2. Quantitative results comparison for perfectly symmetric datasets.

	Ruppert et al. Algorithm			Proposed Algorithm
	z Score (Voxels)	Angle Difference (°)	Time (s)	z-Score (Voxels)	Angle Difference (°)	Time (s)
Mean	1.246	0.10	35.02	0.336	0.06	1.04
Std.	2.041	0.22	1.12	0.324	0.21	0.02
Median	0.50	0.00	34.86	0.250	0.00	1.01

.

The plot of average z-distance for the individual scan is illustrated in Figure 10. Ruppert et al. approach was unable to truly estimate the MSP in some scans and showed substantially large values for average z-distances, as indicated by the black square in Figure 10. These scans were not considered when measuring the mean z-score and the mean of angle differences, since they would synthetically intensify the z-score and standard deviation results of Ruppert et al. algorithm.

The values of z-score are almost similar as they reported in their paper, except for some scans having a rotation angle (either yaw or roll) of greater than 5°, as displayed in Figure 11. The situations at which their algorithm relied, i.e., smoothing, Sobel operator followed by thresholding of 5% brightest voxels, and symmetric score, are not satisfied in the presence of noise, image artifacts, and intensity inhomogeneity. This causes the offsets of plane underestimated. From Table 2, it is evident that Ruppert et al. algorithm accurately estimated the orientation of the MSP but failed to correctly calculate the offset of MSP. On the other hand, the proposed algorithm results are more consistent and accurate in orientations as well as in offsets.

To precisely inspect the accuracy of the proposed algorithm, many illustrative slices orthogonal to the estimated MSP are displayed in Figure 11. The lines in different colors (magenta, blue and green) are the intersecting lines between the extracted MSP and the respective orthogonal slice. The first column represents the input images and the second column shows the input images with the GT MSP intersection line (magenta color line). Similarly, third and fourth columns contain the images of Ruppert et al. and proposed algorithm results, respectively. Visual comparison in Figure 11 also reveals that the proposed algorithm outperformed Ruppert et al. algorithm in terms of accuracy, both in orientation and offsets. Ruppert et al. algorithm could not always achieve a rigorous estimate of MSP, particularly when the brain MRIs underwent a considerable transformation (rotation, translation, noise).

3.3. Evaluation and Comparison on Real Datasets

The proposed method has also been compared with Ruppert et al. method on real brain MRIs. Both techniques have been tested on 125 real head scans of the NFBS database [48] and 18 real head volumes of the IBSR database [49]. No GT for MSP is available for NFBS. Therefore, only visually comparison of the MSP extraction results has been reported for both the algorithms. The first row of Figure 12 shows some of the input image slices orthogonal to MSP.

Extracted MSP (green lines) using the proposed algorithm is displayed in the second row of Figure 12. Blue lines in the third row of Figure 12 displays the MSP extraction results of Ruppert et al. algorithm. The same pattern of the results in detecting MSP is shown by Ruppert et al. algorithm in real brain volumes. It estimates the orientation of the MSP more accurately but fails to correctly calculate the offset of MSP. On the other hand, the proposed algorithm detects the MSP more precisely and consistently.

The IBSR database contains the manual delineation (GT) of brain hemispheres. For illustration, the input image from a brain volume is shown in Figure 13a with its manual delineation (mask) of the left (dark) and the right (bright) brain regions (Figure 13b). For comparison purpose, the boundary of the right region is superimposed on the input image with the white pixels using a morphological gradient and the binary skeletonization [51], as shown in Figure 13c. Therefore, this image is used as a GT for evaluation of MSP extraction results. For the precise and accurate result, the intersection line (between the estimated MSP and the corresponding orthogonal slice) should be overlapped or coincided with the white pixels boundary in the GT image.

The first row in Figure 14 consists of all the input images with GT delineation of the right brain. Since the MSP divides the brain into two roughly symmetrical regions, the detected MSP (green line) by the proposed algorithm is almost completely overlapped with the boundary of the GT delineation of the right brain, as shown in the second row of Figure 14. In contrary, the MSP (blue line) estimated by Ruppert et al. algorithm is deteriorated significantly from the real boundary of the GT delineation of the right brain. Note that we did not compare results of brain tumors datasets with Ruppert et al. algorithm because they did not report results on such datasets in their paper.

In short, all the promising results given by the proposed algorithm indicates that the developed technique has the highest accuracy and consistency in extracting the MSP. Finally, the results of automatic MSP detection and tilt correction (recenter and reorientation) in brain MRIs are displayed in Figure 15, where the first row represents the tilted slices of the three distinct brain volumes. After computing the parameters of the MSP and affine transformation, we reoriented and recentered the brain volumes. The corrected volume images are displayed in the second row of Figure 15.

4. Conclusions

In this paper, we have presented a fully automatic and computationally efficient MSP extraction and tilt correction technique in brain MRIs. The proposed method is based on PCA and symmetric features followed by plane fitting using orthogonal regression. Experimental results on 157 real heterogeneous brain MRIs including 14 datasets with brain tumors and comparison with a state-of-the-art method have confirmed that the proposed technique provides consistent performance with the highest accuracy. Moreover, it is 30 times faster than the competitor algorithm and takes only 1.04 s (on average) for all algorithmic steps per MRI volume. It is also robust on pathological brain MRIs having various intensity inhomogeneities, noises and image artifacts.

Some limitation of the proposed algorithm should be taken into consideration. The algorithm can only take T₁-weighted MRI of the brain in axial orientation as an input. Although one can convert from one orientation (sagittal or coronal) to other using existing algorithms, the selection of brain slice in the first step can affect the yaw angle due to the assumption that it should be the same in each axial slice (see Equation (4)). The reason is that in the true anatomical structure of the brain, MSP is not exactly the plane but a curved surface even for a normal brain. Even though the planer estimation is adequate for many applications such as registration and symmetric/asymmetric analysis of brain images, the results obtained from the PCA in estimating the yaw angle can also affect the performance of the algorithm in the presence of high level of noise. When the SNR is less than −10.09 dB, estimates of the yaw angle cannot be trusted. This problem can be circumvented by increasing the angle range in the cross-correlation step. Our future work will include evaluation of the algorithm on brain images of various sources and modalities such as T₂-weighted, Proton Density (PD) weighted, Fluid Attenuated Inversion Recovery (FLAIR), PET, and SPECT.