Ischemic stroke segmentation in multi-sequence MRI by symmetry determined superpixel based hierarchical clustering

Journal Pre-proof Ischemic stroke segmentation in multi-sequence MRI by symmetry determined superpixel based hierarchical clustering Anusha Vupputuri, Stephen Ashwal, Bryan Tsao, Nirmlaya Ghosh PII:

S0010-4825(19)30395-6

DOI:

https://doi.org/10.1016/j.compbiomed.2019.103536

Reference:

CBM 103536

To appear in:

Computers in Biology and Medicine

Received Date: 17 July 2019 Revised Date:

7 November 2019

Accepted Date: 7 November 2019

Please cite this article as: A. Vupputuri, S. Ashwal, B. Tsao, N. Ghosh, Ischemic stroke segmentation in multi-sequence MRI by symmetry determined superpixel based hierarchical clustering, Computers in Biology and Medicine (2019), doi: https://doi.org/10.1016/j.compbiomed.2019.103536. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Ischemic Stroke Segmentation in Multi-sequence MRI by Symmetry determined Superpixel based Hierarchical Clustering Anusha Vupputuria,, Stephen Ashwalb , Bryan Tsaoc , Nirmlaya Ghosha a Department

of Electrical Engineering,Indian Institute of Technology, Kharagpur, India-721302 b Department of Pediatrics, Loma Linda University, Loma Linda, CA-92354, USA c Department of Neurology, Loma Linda University, Loma Linda, CA 92354, USA

Abstract Automated estimation of ischemic stroke evolution across different brain anatomical regions has immense potential to revolutionize stroke treatment. Multisequence Magnetic Resonance Imaging (MRI) techniques provide information to characterize abnormal tissues based on their anatomy and physical properties. Asymmetry of the right and left hemispheres of the brain is an important cue for abnormality estimation but using it alone is susceptible to occasional error due to self-asymmetry of the brain. A precise estimate of the symmetry axis is therefore essential for accurate asymmetry identification, which holds the key to the proposed method. The proposed symmetry determined superpixel based hierarchical clustering (SSHC) method initially estimates the lesion from interhemispheric asymmetry. This asymmetry further determines the thresholding parameter for hierarchically clustering the superpixels leading to an automated and accurate lesion delineation. A multi-sequence MRI based pipeline also combines the estimations from individual sequences. SSHC is evaluated on different sequences of the Loma Linda University (LLU) dataset with 26 patients and the Ischemic Stroke Lesion Segmentation (ISLES’15) dataset with 28 patients. SSHC eliminates the need for manual determination of threshold for combining

Email addresses: [email protected] (Anusha Vupputuri), [email protected] (Stephen Ashwal), [email protected] ( Bryan Tsao), [email protected] (Nirmlaya Ghosh)

Preprint submitted to Journal of LATEX Templates

November 11, 2019

the superpixel clusters and is more reliable as it derives the information from the quick estimation of asymmetry. SSHC outperforms the state-of-the-art resulting in a high Dice similarity score of 0.704±0.27 and a recall of 0.85±0.01 which are 6% and 35% respectively higher than the challenge winning method. SSHC thus demonstrates a promising potential in the automated detection of (sub-)acute adult ischemic stroke. Keywords: Axial-symmetry, Hierarchical clustering, Ischemic stroke, Multi-sequence MRI, Superpixel

1. Introduction The majority of acute arterial strokes are ischemic (∼85%) in nature [1]. Ischemia results from inadequate blood supply to the brain tissue which in turn hinders the oxygen supply causing hypoxia. Detection, segmentation and quantitative assessment of ischemic stroke lesions using Magnetic Resonance Imaging (MRI) provide valuable information to analyze the inherent neurological pathologies, identify and monitor disease progression, and plan treatments as well as predict outcomes [1]. However, delineation of ischemic stroke injury is challenging because of the ambiguity and complexity of lesion shapes and the difficulty in obtaining the groundtruth for comparision, which requires high clinical expertise and anatomical knowledge. This study focuses on a new automated technique, Symmetry determined Superpixel based Hierarchical Clustering (SSHC), to segment the ischemic stroke lesion from inter-hemispheric asymmetries. We propose that SSHC can appraise the volume and chronology of both acute and sub-acute ischemic stroke and has the future potential to estimate any salvageable tissue present around necrotic core. Diffusion Weighted Imaging (DWI), which is the common MRI sequence used in stroke protocols, is observed to be efficient and reliable because of its sensitivity to ischemic tissue [1]. At the same time DWI is less sensitive to lesion variations in the sub-acute phase following stroke. Other MR sequences include Fluid Attenuated Inversion Recovery (FLAIR) and T2-Weighted Imag-

2

ing (T2-WI) which are considered for correlating injury severity with the final outcome [2]. These sequences often exhibit significant variations in occurrence, lesion shape, volume and signal intensity and often fail to differentiate between the pre- and post- ischemic lesions which hinder accurate assessment of acute stroke treatment. FLAIR lesion contours interfere with anatomical structures (and even calcification or bleeding) and display higher intensity which can be eliminated only by manual correction even after automated segmentation [3]. Therefore, relying on a single MRI sequence for estimating diverse ischemic lesions might be insufficient for reliable treatment. Significance of multi-sequence MRI data: MRI multi-sequence lesion segmentation provides the advantage of combining diverse information inherent to different imaging modalities. DWI, which relies on diffusion of water assesses the morphology of ischemic tissue in acute focal regions, characterizing cytotoxic edema. These injury regions appear hyperintense on DWI whereas in the sub-acute and chronic phases during the first 4 days of the ischemic event they appear normal to hypointense [4], [5]. Chronic stroke lesions which appear hyperintense on FLAIR and T2 sequences can also appear hyperintense on DWI because of the T2 shine through phenomenon [6]. ADC maps differentiate the falsely hyperintense chronic lesions from acute ischemia lesions, which also have an accentuated hyperintensity on both DWI and the ADC maps. This makes ADC maps important for the correction of T2 shine through effect in DWI. The major difficulties in performing multi-sequence imaging are getting the same patient imaged for multiple sequence scans around the same time. In addition, there is a requirement of accurate physical alignment, through registration of multiple scans in order to compare and fuse the detected lesion areas from different sequences. Scant research is performed on the area of multi-sequence ischemic lesion segmentation [7] and most of the methods concentrate on major abnormalities like ischemic stroke detection, tumor segmentation [8] and multiple sclerosis lesion delineation [9] from a single MRI sequence and validate them using small datasets. Challenge: Delineation of ischemic stroke injury is challenging because 3

of the complexity of lesion shapes and obtaining a ground truth is difficult as it requires high clinical expertise and anatomical knowledge. Multiple MR sequences have multi-spectral information of lesion which has to be fused for the final infarct assesessment. Approach: Brain is mostly symmetric along the axial plane while the stroke often occurs asymmetrically i.e., in only one hemisphere of the brain or generally asymmetric even for bilateral stroke. This study focuses on a new automated technique, SSHC, to segment the ischemic stroke lesion from inter-hemispheric asymmetries. After injury is identified and approximated from its asymmetry, the accuracy of the detection is further improved by preserving the injury boundary using superpixel based clustering. Impact and Contributions: The main impact and contributions of the proposed SSHC method are described as follows: • The proposed method is completely automated, eliminating the need for manual interventions such as seed selection and threshold determination for clustering or classification. It is framed on the superpixel based clustering of injury, leveraged on an asymmetry cue which aids automated delineation of ischemic stroke. • Validation of the SSHC on two different datasets with mutisequence MRI data, and with varying stroke severity levels outperforms the state-of-theart. The task of segmentation is simplified with the use of symmetry which eliminates the need of high-end computations and training. • The evaluation also proposes that diffusion sequences are reliable for acute stroke detection and efficient information can be obtained from the FLAIR sequence as a follow up scan during the sub-acute stage. Rest of the paper is organized as follows. The importance of brain symmetry and the existing methodologies for ischemic lesion delineation are summarized in Section 2. The proposed SSHC technique is elaborated in Section 3 and the comparative results of the SSHC with various state-of-the-art techniques on two

4

challenging datasets, are presented in Section 4. Finally, Section 5 concludes the paper.

2. Related Work 2.1. Brain symmetry and its clinical relevance The human brain exhibits extensive bilateral symmetry based on morphological conditions [10]. The regions of either hemispheres exhibit identical physiological properties like the diffusion, perfusion and similar intensity ranges in a particular MRI sequence. During medical examination of the MRI scans, radiologists intuitively tend to compare both right and left hemispheres of the brain to differentiate pathologies. Although the human brain displays high symmetry, it is not perfectly symmetric. For instance it is reported that due to functional differences, the left hemisphere is greater in size than the right hemisphere [11], mostly due to right-handedness in majority of the subjects. We believe that automated and quantitative asymmetric correlation between the right and left hemispheres of the brain can be a key heuristic to identify pathologies [12]. To evaluate this, the most essential requirement is to determine the symmetry axis or mid-sagittal plane (MSP). Determining MSP also helps to align different radiological scans. Ideally, MSP presents maximum mirror symmetry [11] and there are shape-based and content-based methods in the literature for determining MSP. Inter-hemispheric cerebral fissure is used as a shape based prior to determine MSP but this method is vulnerable to invisibility of the fissure due to other masses [13]. Minovic et al. proposed a 3D symmetry detection method which involves geometric object fitting assuming that the plane of symmetry is perpendicular to the orthogonal X-Y axis [14]. This method results in inaccurate symmetry planes when the dataset is incomplete or deformed (due to acute or chronic injury) and the object fitting is improper. Hachaj et al. presented an automated lesion detection algorithm using perfusion asymmetry from Computed Tomography (CT) images [15]. Most of the parameters used for asymmetry detection are empirically determined and

5

are constrained to that particular dataset where the scanner dependent variations and rotational errors are not considered. Rajini et al. proposed an injured and non-injured hemispheric classification of CT images based on symmetry [16]. The multi-feature vectors from either hemispheres train the classifier to differentiate between injury and non-injury slices. This method resulted in high classification accuracy but the emperical parameters like the fixed brain region diameter, and pre-determined values to determine the asymmetry map limit its adaptability to other datasets. Sun et al. proposed a Scale Invariant Feature Transform (SIFT)-based detection method from perfectly symmetric objects [17] and is adapted by Ghosh et al. for ischemic lesion segmentation from neonatal stroke data of T2-WI [18]. These methods did not take the self-asymmetries and rotational variance in normal brain into account and also require partially hand crafted parametric values like sequence based intensity ranges and area thresholds. In [19], the symmetry axis is determined in animal T2-MRI slices by finding the averaged Inter-Hemispheric Fissure (IHF) to initiate the starting point of the symmetry axis. The symmetry axis is projected onto the local minimum constrained by a manually determined region. Though this method depicts satisfactory symmetry axis detection in the animal data, it might result in multiple minima points due to more number of sulci in the clinical data compared to the animal data. Peter et al. used symmetry to semi automatically rigid register multi-sequences in CT to a pre-determined perfectly symmetric slice [20]. Later the contralateral features are classified to detect the lesion. The asymmetric region detected with this method can be erroneous as the brain shape is highly varying and can be deformed due to its dependency on the reference slice used for the rigid registration. In our previous work, a mid sagittal axis (MSA) detection joining the falx point and the centroid of the MRI volume is proposed which has a neccesity of prominent falx point to be identified without which the later steps of the algorithm could fail [21]. Recently an 8-layer deep symmetry based 3D Convolutional Neural Network (CNN) is proposed by Wang et al. for segmenting stroke lesion which utilizes the concept of feeding symmetric voxels to the network and simultaneously identifying the 6

necrotic tissue [22]. Though this method seemed promising, the methodology for detecting symmetric voxels from both the hemispheres is not mentioned and the proposed network is evaluated only on severe stroke injury and only for the DWI sequence. 2.2. State-of-the-art ischemic lesion segmentation methods Numerous automated and interactive state-of-the-art segmentation algorithms have been proposed for ischemia delineation which include hierarchical thresholding in both animal and clinical studies [23],[24],[25], automated random forest based methods [26]-[27], interactive random walker based methods [28],[23], and deep neural networks based on CNN and stacked autoencoders [29],[30],[31],[32],[33]. Discrete curvelet transform features derived from MRI are fed to the SVM classifier to distinguish injured and non-injured slices [34]. U.Acharya et al. proposed an interesting technique to segment ischemic stroke from only DWI sequence using multiple entropy features derived from higherorder spectra [35]. These features are further fed to an SVM classifier to classify the input slices into three different stroke intensities. A. Subudhi et al. combined Delaunay triangulation for grouping similar pixels along with fractional order Darwinian particle swarm optimization to automatically segment the ischemic lesion. They further classified the lesion into three different stroke severity levels using lesion derived Gray Level Co-occurance Matrix (GLCM) features [36]. The segmentation performance of all these methods hinges on the image features and classification models. Since deep learning techniques are capable of learning high-level and task-adaptive features from training data, they have been adopted for stroke injury segmentation [30]. However, most of the existing segmentation methods based on deep learning yield segmentation results with appearance and spatial inconsistency. Super pixel based methods are being widely used in computer vision tasks and segmentation [37]-[38]. They contribute efficiently to segmentation tasks as they have a major advantage of computational efficiency which can be attributed to similar characteristics of intensity and location features of groups of pixels, instead of single pixels. This 7

property already gives the advantage of working on small groups of clustered pixels and reduces the chance of wrongly labeling individual pixels in the final segmentation. Simple Linear Iterative Clustering (SLIC) method of determining superpixels is fast and has the flexibility of changing the number of superpixels depending on the application and size of the segmentated area [39]. 3D Superpixels have been adopted for segmenting prostate region from MRI using graph-based energy minimization and shape features [40], and quantification and classification of breast cancer cells from histopathological images [41]. Recently SLIC based superpixels and region growing based segmentation of glottal area from video laryngoscopy images is proposed [42]. Multiple methods of superpixel determination like recursive updating of superpixels for more adaptability [37] and K-means based SLIC determination of superpixel have been proposed.

3. Proposed approach for segmenting ischemic lesion using symmetrydetermined superpixel based hierarchical clustering We highlight two novel propositions in the SSHC. • Mid-sagittal axis is determined from the axial MR volumes based on the IHF using multiple Hough lines. This relates to the ’symmetry’ part in SSHC. • The MRI voxels are grouped into superpixels using SLIC. The superpixels are combined hierarchically into a dendogram tree. The predetermined asymmetry information from the inter-hemispheric comparison is used to automatically split the dendogram at an optimal level, segregating ischemic lesion and normal brain matter. This relates to the ’superpixel based hierarchical clustering’ part in SSHC. These contributions make SSHC completely automated, based on a simple but significant symmetry based intuition while keeping the method robust, fast, reliable and efficient by segmenting challenging mild lesions also with high accuracy. Each contribution is further detailed in the following subsections 3.1 and 3.2. 8

The graphical abstract of the SSHC method, as shown in Figure 1, comprises of two components in the pipeline. At first, the brain MSA is determined by eliminating all the effects of rotation. After dividing the brain into two hemispheres they are compared for the initial estimate of asymmetry. Second, the MR images are grouped into superpixels based on the intensity and proximity criterion and then an agglomerative clustering is performed. The asymmetry-determined approximate lesion helps determine the range of intensities or average intensity of the lesion for that particular sequence. This intensity is mathematically related to the threshold of superpixel based clustering to determine the final injury volume. 3.1. Pre-processing and symmetry axis detection The Mid Sagittal Plane (MSP) is the best approximation of the symmetric plane of axial brain and is often considered as the plane orthogonal to the axis of the brain in the axial view. MSP seperates the brain into two hemispheres but tilt or disorientation is common but often not tractable. Reorienting the MR volumes and readjusting the MSP before asymmetrical injury detection would lead to more meaningful data assessment for cross comparison of the two hemispheres. MSP is a volumetric combination of MSAs over individual axial or coronal MRI slices. A 2D based MSA for symmetry axis detection is proposed here because, the common drawback of 3D methods (considering the entire volume) is the high computational cost. To determine the 3D MSA, the optimization algorithm would search over the entire space to estimate the extent of hemispheric mismatch. As a preprocessing step, a Hough transform based multiple line estimation is also proposed in this study to detect the IHF and check for its disorientation. This is followed by the rotational correction for all MRI slices if needed and the experimental validation for the existing patient dataset. The orientation of MR image is estimated from the Mid Axial Slice (MAS) and also the pre- and postconsecutive slices for validation. IHF is best determined from the MAS. So, we detected the edges using the Sobel operator [43]. Multiple disconnected straight 9

Figure 1: Graphical abstract – This flowchart describes the pipeline of the proposed algorithm starting with the determination of symmetry axis for the MR volume and estimation of asymmetric region information which is fed to the hierarchical clustering model along with the Simple Linear Iterative Clustering (SLIC) super pixels to agglomeratively cluster and classify the superpixels as lesion or Normal Appearing Brain Matter (NABM). Acronyms- SSD: sum of squared distances.

10

lines using Hough transform are detected from the edge image as shown in Figure 1. Instead of relying on a single Hough line, the proposed method considers multiple Hough lines and gives a reliable approximate of the disorientation as shown in Figure 2. Each K th Hough line is identified by the upper and lower end points Pku (x, y) and Pkl (x, y). To determine the initial estimate of the symmetry axis, linear regression is performed on the data pairs (xi , yi ) of collective end points as in (1).   (xi , yi ) ∈ Pu1 (x, y) , Pl1 (x, y) , . . . , Puk (x, y) , Plk (x, y) k

(1)

u

∈ ∪ (Pi (x, y), Pli (x, y)) i=1

Linear regression to find the approximate estimates of the coefficients m, ˆ n ˆ of the best-fit line is computed by minimizing the Least Square Error (LSE) given by O(x,y) in (2)-(3). y = m + nx

O (x, y) =

k X

εˆ2i =

i=1

k X

(2)

(yi − m − nxi )2

(3)

i=1

Finally, approximates of the coefficients of best fitting MSA is given by (4)-(5). Pk n ˆ= 1 where x ¯= N

¯)(yi − i=1 (xi − x Pk 2 ¯) i=1 (xi − x k X i=1

y ¯)

,

k 1 X xi and y ¯= yi N i=1

m ˆ =y ¯−n ˆx ¯

(4)

(5)

Orientation of the IHF affects the estimation of MSA due to inaccurate points as determined by disoriented Hough lines. The accuracy of initial estimate of the MSA is improved by initial rotational correction of the MRI volumes. For this correction the initial MSA is estimated from the disoriented or rotated image using (1)-(5) and it’s orientation is corrected with reference to the vertical 900 axis using the angular orientation difference (Orientdiff ) given by (6)-(7). l PU sym (x, y) and Psym (x, y) are the finite end-points of the symmetry axis of the

11

Figure 2: Refined MSA detection – (a) Multiple Hough lines in green detected on the edge detected brain slice. (Red and yellow crosses show the start and end point of the Hough lines) (b) Binary window mask to refine the Hough lines so that MSA detection is unbiased by vertical edges detecting further away from the MSA. The red dot denotes the volumetric centroid and L is the length of the minor axis of the brain area. (c) Refined Hough points used for Least square error minimised MSA detection. (d) Accurate MSA determination for dividing the brain into left and right hemispheres.(Acronyms: MSA: Mid saggital axis)

12

rotated slice. The remaining slices of the entire MRI volume are also corrected using the (Orientdiff ) in (6). Further, on the rotational corrected volume the multiple Hough line based MSA estimation is performed using the same procedure described above. After rotational correction of the MRI volume, the MSA is very much in alignment with the orthogonal axis. Therefore, to estimate the final and accurate MSA by regression based line fitting, only the Hough points lying in a longitudinal window of width β = ±6% from the brain center are considered as shown in Figure 2. This β value is determined by visually validating the alignment of the MSA on the MAS of both the datasets.

Orientdiff = 90 − Orientsym axis ( −1

Orientsym axis = tan

l PU sym (y) − Psym (y) U Psym (x) − Plsym (x)

(6) ) (7)

3.2. Hemisphere based comparison for injury detection The symmetry axis is essential to divide the MR slice into two symmetric hemispheres and the asymmetry of the hemispheres can be marked as injury. To determine the injured part, patch-wise comparison of both the hemispheres is performed using the Sum of Squared Distances (SSD) between the intensities of the left and right hemispheric patches (pr , pl ) using (8). The centroids of these patches are located symmetrically equi-distant from the MSA. SSD [pr , pl ] =

size pX size pX i=1

[pr (i, j) − pl (i, j)]2

(8)

j=1

A simple and quick Otsu thresholding method [44] is followed on the difference image to identify the major asymmetric region. The median intensity of this identified asymmetric region is asymed . Small outliers in this asymmetric region are rejected using morphological area based connectivity and outlier rejection. But the decision based only on SSD fails to indicate the hemisphere which is actually injured because it only detects asymmetric regions from the difference of both hemispheres. To avoid this ambiguity and also to determine 13

the hemisphere which is injured, a sequence based information criterion is taken into account. The median intensity of the asymmetric region is calculated separately from either hemisphere. For instance, the lesion contrast in DWI and FLAIR is higher compared to the Normal Appearing Brain Matter (NABM). Therefore the median intensity of the lesion Region of Interest (ROI) in the injured hemisphere is higher than the median intensity of its symmetric counterpart in the other hemisphere. Assuming the lesion is in the right hemisphere, subtracting the median of the asymmetric region of the left hemisphere from the injured part of the right hemisphere will lead to a positive value because the injury results in higher intensity. The same differencing approach will result in a negative value if the lesion is in the left hemisphere. The positive or negative value of the median intensity difference can thus determine the injured hemisphere. For a sequence like Apparent Diffusion Coefficient (ADC) in acute stroke where the injury is darker (hypo-intensity) than NABM, the difference will be positive if injury is on the left hemisphere and vice-versa. In this case the median intensity of left asymmetric region is separated from the right hemisphere’s asymmetric region. Injury intensity contrast based prior information from the sequence is thus combined with the symmetry cue to determine the injured hemisphere. 3.3. Symmetry determined superpixel based injury segmentation Asymmetry based lesion estimation alone has the disadvantage of overestimating or under-estimating the injury due to self-asymmetry of normal brain. Therefore, superpixel based clustering of the MR slices utilizing the intensity information of the asymmetry detected lesion is proposed to differentiate between the lesion and NABM. Algorithm 1 describes the sequential procedure for hierarchically clustering the superpixels and using the asymed determined earlier to differentiate injury and NABM. Each slice of size M × N is divided into K superpixel clusters with cluster center placed at an interval of p Sin = (M × N ) /K . Based on the intensity (l) difference criterion Din cri and spatial proximity criterion Dsp cri given in (9)-(11), the pixels surrounding 14

the corresponding cluster center are grouped into superpixels. Parameter α determines the relative weightage of these criteria and higher the value, the more regular (circular in size) are the superpixels. Din cri = |li − lj |

Dsp cri =

q

2

(9)

2

(xi − xj ) + (yi − yj )

Dres = Din cri +

α Dsp cri Sin

(10)

(11)

By experimentation, K =1000 is selected (Figure 6 plots the accuracy Vs. no. of superpixels graph) for accurately assessing severe to mild lesions. After superpixel estimation, the goal is to differentiate superpixels of lesion from NABM. To further divide the superpixels into multi-class regions a hierarchical clustering technique is adopted which is driven by the asymmetry based intensity information obtained earlier. Based on the dissimilarity measure of the superpixel with its neighboring superpixels, a dendogram tree is formed. The dendogram exhibits the connectivity between every two closest superpixel clusters by combining them in a bottom-up approach. It is further split at an arbitrary level to combine superpixel clusters into R arbitrary regions. The median of the asymmetry region: asymed decides the splitting threshold Gth to classify the superpixels into lesion and NABM, and thus simultaneously preserving the superpixel boundaries. 3.4. Multi-sequence data fusion model Combining multiple estimates from multi-sequence MRI data for ischemic lesion assessment and delineation is essential for accurate analysis. Ischemic lesion can be estimated from multiple sequences of ADC and DWI (using three orthogonal planes of p, r, and s), as in the LLU dataset. This eliminates certain irregularities like T2-shine through which are mostly evident in MRI scans involving ADC and also improving the overall accuracy of segmentation. In the

15

Algorithm 1: Pseudo-code for superpixel(SP) based clustering determined by axial-asymmetry Input: MRI volume to be segmented, Number of SPs= (K ) Input: Median of Asymmetric region intensity (asymed ) Output: Each voxel labelled as lesion or NABM 1

Initialize superpixel cluster centers CK = [IK , XK , YK ] at intervals Sin

2

Iteratively cluster pixels based on Dres in 11 to form final SPs

3

for K ← 1 to N do

4 5

SPmean (K)← cluster mean of intensities of the k th SPs end /* Hierarchical superpixel clustering tree

*/

/* For K SPs, lowest level (level low) of the hierarchical cluster tree starts with (K(K − 1))/2 leaf nodes and ends at a highest level (level high) with 2 nodes 6 7

*/

for cluster level ← level low to level high do for each SP pair (a,b) in a cluster level do Similarity index of SPs pair (SPa , SPb ) ← Euclidean

8

distance(SPmean a , SPmean b ) 9 10

end Cluster pairs of superpixels into new clusters based on highest similarity index

11

Update new cluster nodes

12

Increment cluster level

13

end

14

Split Hierarchical superpixel clustering tree into R arbitrary regions based on single linkage of clusters.

15

Classify the R Hierarchically Clustered Superpixels (HCS) as injury or NABM based on asymed obtained from asymmetry of hemispheres

16 17

for j ← 1 to R do if mean(HCSj )>asymed then label HCS as injury

18 19

end

20

else label HCS as NABM

21 22 23

end end

16

proposed multi-sequence fusion method, MR volumes of individual sequences are skull-stripped using binerization followed by morphological cleaning. SSHC is applied to each of the multi-sequence MR volumes of the same patient and the individual lesion estimates are obtained as shown in Figure 3. The final ischemic injury is decided by the intersection of DWIp, DWIr and DWIs injury (hyperintensity) followed by union with ADC injury (hypointensity). This includes the entire lesion ROI reflected in either of the DWI scans and then only that part of the lesion is considered which is common to both DWI and ADC, eliminating any effect of T2-shine through in DWI. Similarly for the ISLES0 15 dataset ischemic injury is decided by union of T1-WI, T2-WI and DWI injury followed by intersection with FLAIR lesion.

4. Results and Discussion The proposed SSHC method is experimentally validated on multi-sequence MRI of two widely varying datasets. The consistency in performance of the algorithm is further evaluated in terms of qualitative measures. Individual analysis of the algorithm on multiple sequences with highly varying properties is essential for evaluation of the significance of that particular sequence in ischemia detection. 4.1. Dataset and ground-truth generation The proposed symmetry-determined SSHC method is tested on two sets of multisequence MRI data whose details are mentioned in Table 1. 4.1.1. Loma Linda University (LLU) Dataset The first clinical dataset (IRB-approved) is obtained from the Departments of Neurology and Radiology, Loma Linda University (LLU), Loma Linda, California, USA. For this research, 26 adult stroke patients (>18 years of age) admitted to the medical center are selected based on the following inclusion criteria: (1) hyperacute ischemic stroke (admitted <12 hours from stroke onset), (2) visually detectable injury in DWI, and (3) high National Institute 17

Figure 3: Multisequence estimation –The proposed algorithm is applied only on brain mask resulting in individual lesion estimations from ADC, DWIp, DWIr and DWIs. Then the lesion ROIs are fused by logical operations (OR’ed between DWI lesion ROIs and AND’ed with ADC lesion ROI) to produce final lesion estimation – where estimated result (green) and ground-truth (red) are depicted. (Acronyms: ROI: Region of interest; ADC: Apparent diffusion coefficient; DWI: Diffusion weighted imaging)

18

Table 1: Details of the ISLES and LLU Datasets

Dataset

ISLES- SISS

Number of

Multi-Sequences

Volume

Stroke

Subjects

available

dimensions

severity

230 × 230 × 154

Sub-acute

192 × 192 × 35

Acute

28

FLAIR, DWI, T1-WI, T2-WI

LLU

26

DWIp, DWIr, DWIs, ADC

of Health (NIH) stroke scale score (>10). MRI data selection of the patients and its analysis is performed retrospectively after proper de-identification. It is well-established that manual delineation is considered the ’gold standard’ for comparison in any lesion quantification research. For robustness, manually demarked groundtruth of the lesions are generated by fusing volumetric intersection (logical AND-ing) across multi-timepoint delineations by the same observer (to reduce intra-observer variability). Further, volumetric union (logical ORing) of detections by multiple observers reduces the inter-observer variability [45]. The clinical dataset of 26 patients from LLU has multispectral data of 4 sequences: ADC, along with DWI in the three planes of DWIp, DWIr, and DWIs. 4.1.2. ISLES’15 Dataset The second dataset is the 28 patient sub-acute ischemic lesion segmentation (SISS) challenge dataset of ISLES’15 which is held as a part of medical image computing and computer assisted intervention (MICCAI) conference. This dataset consists of four sequences: DWI, FLAIR, T1-WI, and T2-WI [46]. Lesions are considered as sub-acute infarct if the pathologic signal of that MRI volume is identified simultaneously in FLAIR and DWI images. Infarct lesions with signal variation due to ischemic transformation are only included. The ground truth for this dataset is determined by one radiologist manually delineating the injury from a single sequence at a single time point [47]. 19

(a) GCAT

(b) ISLES

Figure 4: Ischemic lesion delineation from multiple modalities of GCAT and ISLES’15 datasets are shown here and performance of the algorithm is compared across various injury classes.

20

Automatically estimated result (green) and ground truth (red) are superimposed for comparison. Note that the ground truth for ISLES is generated from FLAIR data only. (All the results are generated using K =1000 superpixels)

4.2. Experimental Results The proposed SSHC method has been evaluated qualitatively and quantitatively on both the LLU dataset and ISLES’15 dataset by comparing automaticallyidentified ischemic lesions and manually identified ground-truth. The injury is estimated using parametric measures of precision (12) which estimates what amount of the lesion predicted is actually lesion; recall (13)or sensitivity estimates how well the positive lesion voxels are identified by the SSHC algorithm with respect to the groundtruth; Dice Similarity Score (DSC )(14), which measures the extent of overlap between the lesion detected and grondtruth. Accuracy (15) measures how well both the lesion and NABM voxels are predicted by SSHC. Ground truth labels (grnd lbl), predicted labels (pred lbl) and brain map (Br map ; entire brain region excluding the background and skull) are used for the qualitative comparison of segmentation by the following parametersP recision =

Accuracy = 1 −

(12)

(pred lbl ∩ grnd lbl) grnd lbl

(13)

2(pred lbl ∩ grnd lbl) pred lbl ∪ grnd lbl

(14)

Recall =

DSC =

(pred lbl ∩ grnd lbl) pred lbl

(pred lbl ∪ grnd lbl) − (pred lbl ∩ grnd lbl) Br map

(15)

To explicitly assess the performance of the lesion detection pipeline, three levels of injury severity are decided based on the percentage of total injury volumes of the patient with respect to the entire brain ROI volume. The injury volumes ranging above 30% of brain map are categorized as severe injury, and those volumes ranging in between 15-30% are categorized as moderate injury and less than 15% are considered as mild cases [45]. The algorithm is implemented and all the simulations are performed on a MATLAB environment using 16GB RAM. Performance of the algorithm on both injury severity as well as multiple modalities is depicted in Figure 4 for both ISLES’15 dataset and LLU 21

Table 2: Performance measures of the three injury severity classes of ISLES Dataset

ISLES’15

Severe

Moderate

Mild

Average

Precision

0.713±0.18

0.645±0.21

0.674±0.02

0.677±0.14

Recall

0.871±0.00

0.887±0.01

0.817±0.00

0.8•5±0.01

Dice score

0.79±0.17

0.673±0.39

0.648±0.34

0.704±0.27

Accuracy

0.835±0.02

0.941±0.02

0.865±0.02

0.880±0.02

datasets. It also demonstrates the efficiency of the algorithm in both datasets by differentiating either hyperintensities or hypointensities based on a particular MRI sequence. 4.2.1. Performance evaluation the proposed SSHC method on ISLES’15 and the LLU datasets We observed that the SSHC method provides comparatively better performance for the LLU dataset in the three injury severity classes compared to the ISLES’15 dataset. Table 2 shows the severity classes of the ISLES’15, Table 3 shows the the LLU dataset, and in Table 4 the proposed SSHC method is compared with the top 4 methods [48], [49], [50], [51] of the ISLES’15 challenge along with the state-of-the-art methods in [52], [53], [54], [55] and [56]. It is to be noted that the comparison techniques are re-implemented using the code available from few of the authors for the parametric measures and the remaining values are retained from the respective published results available. Compared to all the recent methods published in recent years SSHC outperformed all the methods in terms of DSC for the SISS dataset yielding a value of 0.704 which is 6% higher than the 3D CNN method which was the challenge winner for ISLES-SISS. The sensitivity value of 0.85 is 35% higher than the best method. Achieving this improvement in performance without any training phase unlike machine learning and deep learning algorithms increases its ease of deployment in the scanner software in the future.

22

Table 3: Performance measures of the three injury severity classes of LLU Dataset

LLU

Severe

Moderate

Mild

Average

dataset Precision

0.924±0.06

0.728±0.05

0.729±0.10

0.794±0.07

Recall

0.983±0.01

0.995±0.00

0.99±0.00

0.993±0.00

Dice score

0.823±0.04

0.754±0.04

0.80±0.08

0.793±0.05

Accuracy

0.979±0.01

0.989±0.00

0.99±0.00

0.989±0.00

Table 4: Qualitative parametric comparison of the proposed SSHC method with the state-ofthe-art (best result in bold)

Algorithm

Method

DSC

Precision

Recall

Konstan

3D CNN+

0.66±0.24

0.683 ±0.24

0.63 ± 0.25

et al. [48] [29]

CRF

Chaolu

Fuzzy

0.55±0.30

0.65 ± 0.26

0.69 ± 0.23

Feng et al.[49]

C-means

Hanna

Random

0.47±0.32

0.53 ± 0.22

0.68 ± 0.25

Halme et al.[50]

Forest

Syed

Local

0.43±0.27

0.51 ± 0.25

0.79 ± 0.15

Reza et al. [51]

Gradient

Zhang et al.[52]

3-D Densenets

0.58±-

0.6±-

0.68±-

S.B.Melingi et al.[53]

Sample

0.62±-

0.952±-

0.65±-

0.56 ± 0.26

0.594 ± 0.24

0.422 ± 0.05

weighted RF Subbanna N.K.

Customised

et.al [54]

MRF

A.Gautam et al. [55]

RF

0.67±-

0.70±-

0.71±-

R.Karthik et al. [56]

FCN

0.70±-

-

-

Proposed method

SSHC

0.704±0.27

0.677±0.14

0.85±0.01

23

4.3. Discussion As the initial estimation of brain injury is dependent on the asymmetry of the right and left hemispheres from axial brain images, handling disoriented brain volumes is a major challenge, though not yet addressed by most of the symmetry detection methods [17]. In addition, most of the methods reviewed in the literature of [24] required manual determination of seed points to guide the algorithm to initiate lesion estimation. The SSHC method accurately and automatically determines the MSA by realigning disoriented volumes. Symmetry-based intensity determination of injury eliminates the need of manual determination of threshold unlike most intensity threshold methods in the literature [17],[18],[28] and is robust to changes in modalities. Further, superpixels are efficient in preserving the boundary and in accurately estimating the injury. This combined method of using symmetry estimates to determine superpixel clustering is completely automated. SLIC superpixels proposed by Achanta et al. [39] are adopted in this work for lesion estimation as they have a major advantage of adherence to contrast based boundary differentiation which is essential for ischemic lesion segmentation. After dividing the MRI slice into superpixels a relevant clustering or thresholding is necessary to differentiate superpixels belonging to injury and NABM. Depending on the MRI sequence chosen for segmentation at least a rough estimate of the lesion intensities and parameters (like minimum area threshold) are required in the prior art [25]. Whereas, the current study proposes utilizing symmetry based information to determine intensity threshold to cluster and differentiate superpixels for final injury estimation. Hierarchical clustering being a cluster analysis method, combines regions or pixels into a hierarchy of clusters based on relevant criteria [57]-[58]. It also requires manual determination of the number of clusters into which the resulting clusters are to be segregated. Taking the asymmetry determined information into account this manual intervention is also eliminated and the proposed methodology gives accurate results with lesser computational expense. One major issue to be addressed in ischemic lesion detection by the automated algorithms is erroneous identification of white matter hyperintensities(WMH) as ischemic lesion [59]. 24

WMH in most of the data are observed to be symmterically existing near the ventricles. As SSHC bases its cue for lesion identification on asymmetry, it automatically excludes the WMH injury from interfering the ischemic lesion. 4.3.1. Validation of SSHC on ISLES and LLU datasets Selection of the number of superpixels (K): All experiments are evaluated with 1000 superpixels per MR slice. Increase in number of superpixels might increase the performance of few measures like recall of lesion segmentation but it will also increase the computational complexity. To maintain the tradeoff between accuracy and computation the optimum number of superpixels is chosen to be 1000 which gave satisfactory segmentation results across most of the data. Though higher number of superpixels may improve the recall but higher number of superpixels require greater execution time across the volume of a single sequence. The evaluation presented in Figure 5 compares the superpixel map, segmentation map and execution time for varying superpixel count. To justify the selection of the number of superpixels K, all the performance measures are evaluated on the maximum area slices from all patients in both the datasets and the average values are compared in the Figure 6. Though true negative rate or recall is slightly better for K =1500 superpixels, considering the overall performance for other measures, K =1000 is chosen as the optimal number of superpixels as shown in Figure 6. 4.3.2. Multi-sequence data validation The ground truth for the ISLES data is obtained by manual delineation on only FLAIR sequence. This estimation differs greatly from the lesion reflected on other MRI sequences like DWI and T1-WI. Because all MRI sequences in this dataset are previously co-aligned among each other, the groundtruth (detected from a single sequence) remains the same for all MRI sequences, while automated estimation from a particular MRI sequence generally varies significantly. In the severe class of DWI data (Figure 4) it can be observed that the automated detection or SSHC does not match with the groundtruth because it

25

Figure 5: Selection of the optimum number of superpixels with a tradeoff between performance and execution time. Top row shows the variation in superpixel clustering with varying number of superpixels. Bottom row shows the segmentation result for the varying superpixel count.

Figure 6: Selection of the optimum number of superpixels with comparison of performance measures.

26

Figure 7: Detection of multifocal (distributed) injury from different data and modalities of both ISLES and LLU datasets. These results show the performance of the algorithm in estimating distributed injury in these example MR slices. Automatically estimated result (green) and ground-truth (red) are superimposed for comparison.

is delineated from FLAIR and not from DWI. 4.3.3. Time complexity analysis of SSHC compared to state-of-the-art The time complexity of an algortihm is dependent on multiple factors like processor configuration, RAM, usage of GPU in case of high end machine learning and deep learning algorithms. Considering these variabilities the time complexity of SSHC is compared with the existing state-of-the-art methods already mentioned in Table 4. The common base methods from all the available methods in Table 4 are chosen and included in the Table 5. It can be observed that the advantages of deep learning methods come with high computational expense of training time and manually tuning the hyperparameters by the expert. Though machine learning methods like Random Forest (RF) and Markov Random Field 27

Table 5: Time complexity comparison of the proposed SSHC method with the state-of-the-art

Algorithm

3D-CNN

Sequences

Training time

Inference time

RAM

included

(in hours)

/volume (in s)

in (GB)

∼6.3

247.34

NA

FLAIR, T1, T2, DWI FLAIR, T1,

GPU

NVIDIA TitanX-12GB

None

20.44 × 103

NA

No GPU used

FLAIR

None

45.59 × 103

NA

No GPU used

Random forests (RF)

FLAIR

∼0.817

12.78 × 103

NA

No GPU used

Fully connected

FLAIR, T1,

∼6.5

150.92

32

Neural Net (FCN)

T2, DWI

Fuzzy C-means

T2, DWI Markov Random Field (MRF)

Proposed SSHC

FLAIR, T1,

NVIDIA Quadro P4000 GPU-8 GB

None

187.57

16

T2, DWI

(MRF) are reliable, they too require few hours of training time and higher inference times compared to SSHC which can be enhanced by utilizing a GPU. But these methods also rely on feature selection and parametric choices like depth of tree, number of trees and atlas generation which can highly vary the end-results. Though FCN has lowest inference time of 151 s, this is a result of usage of GPU and higher RAM. Therefore SSHC is a reliable method utilizing optimized resources and simultaneously providing accurate results. 4.3.4. Comparison between symmetry only method and SSHC Only symmetry based injury detection is highly dependent on the detection of accurate mid-sagittal axis (MSA) and interhemispheric difference is estimated accurately only after rotational correction. Further, the choice of the suitable patch-size for interhemispheric comparison also has its importance in estimating the differences. Symmetry based interhemispheric differencing is good at estimating the approximate area of lesion but fails when lesion has to be ac28

No GPU used

curately determined. The differences in estimation between the symmetry only method with rotational correction and SSHC method are evident from Figure 8. The inaccuracies in highly varying lesions in the 4 slices is presented where symmetry only method fails to detect the scattered lesions. Moreover, there are cases where symmetry only method tends to over-estimate and under estimate the lesion compared to the groundtruth, whereas SSHC accurately delineates stroke in all the patients considered for evaluation. As SSHC method includes the benefits of both symmetry and superpixels, the combination of intensity and local textural information results in efficient boundary preserving leading to highly accurate results. From Figure 8 it can also be observed that only symmetry based injury detection can result in non-injured regions from asymmetric nature of brain which cannot be excluded even after post processing techniques like region connectivity and area based thresholding. SLIC superpixel based detection requires a manual threshold for further clustering the superpixels into lesion and NABM. SLIC superpixel algorithm considers proximity as well as intensity difference while grouping the pixels into superpixels, so it is likely for the algorithm to group the pixels of nearest intensities with minute differences into the same superpixel. For example, in a sequence where the lesion is near to the ventricles has nearly same intensity as the white matter hyperintensities near the ventricles might be grouped together as injury because of spatial proximity. Though it can be overcome by increasing the number of superpixels, this would increase the computational time. But the proposed algorithm eliminates this ambiguity by automatically specifying the appropriate range for clustering the superpixels. The SSHC pipeline for stroke segmentation performs better on the LLU dataset than the ISLES dataset because for the LLU dataset the manual delineation is done by multiple persons and at multiple time points and the result is compared with the resultant ground-truth by minimizing the inter-observer and intra-observer variabilities. The algorithm results in high sensitivity and DSC for severe and mild classes of both the datasets (Tables 2,3), and even performs better in the moderate cases as compared to the stateof-the-art (Table 4). Improvement in detection using proposed SSHC compared 29

Figure 8: Differences between the symmetry only method (1st row) and SSHC method (2nd row). Automatically estimated result (green) and ground-truth (red) are superimposed for comparison.

with both the symmetry based injury detection alone and the superpixel based injury detection alone on both ISLES dataset and the LLU dataset are shown in Figures 9 and 10. 4.3.5. Contribution assessment of each sequence in the multi-sequence data The proposed algorithm fuses multi-sequence detection to provide accurate final injury estimation. Figure 4 also suggests that FLAIR data is more efficient for estimating stroke lesion. The ISLES dataset as mentioned earlier is based on patients with sub-acute stroke where traditionally DWI is more reliable, and using FLAIR for groundtruth generation might be argued. This may explain why a significant variation is observed between automated estimations of multiple sequences of ISLES data. SSHC depicts consistent efficiency in estimating severe to mild injuries. Severe injury estimation though is important, cannot contribute much to further treatment planning (due to less recoverability), whereas automated delineation in moderate and mild injury can help the neuro-

30

Figure 9: Comparison of improvement in detection using proposed symmetry determined superpixel based method with only symmetry based injury detection and only superpixel based injury detection on ISLES dataset.

Figure 10: Comparison of improvement in detection using proposed symmetry determined superpixel based method with only symmetry based injury detection and only superpixel based injury detection on LLU dataset.

31

radiologist to assess and prevent further necrosis by proper treatment planning. This analysis also determines which MRI sequences are essential for determining ischemic injury and which can be eliminated to reduce the imaging as well as computational expenses. To further substantiate the above analysis Region of convergence (ROC) plots for all the sequences of ISLES and LLU dataset, are plotted in Figure 11 and Figure 12 respectively. ROC plots are used here to assess the performance of the SSHC method. Sensitivity versus the fall out rate in ROC measures how well the SSHC method classifies lesion and NABM voxels. In addition, the area under the ROC curve (AUC) is another effective parameter for estimating the performance. Higher the AUC (in the range of 0.5-1) indicates better classification of voxels into the correct classes. The varying parameter used for the ROC plots here is the asymed which is assessed from the inter-hemispheric comparison. It is already mentioned in the methodology that the asymed is used to assess the automated superpixel hierarchical clustering threshold. In order to validate the same, asymed was varied from (0-1) keeping the superpixel count fixed at K=1000. The ROC curve is plotted for each type of sequence present in the dataset and it is validated that highest true positive rate and lowest false-positive rate are approximately obtained for the asymed for that particular sequence. To identify the extent of significance and contribution of each sequence (FLAIR, T2, T1c and DWI in sub-acute stroke; DWIp, DWIr, DWIs and ADC in acute stroke) in lesion estimation, the AUC for each sequence in the SISS and LLU ROC plots are estimated. To plot the ROC, only the hierarchical clustering threshold which is the asymed is varied keeping every other parameter constant for the entire dataset. The AUC of FLAIR in ISLES curve plot is noticeably higher than the other sequences except for T2-WI which is also comparable to FLAIR. Similary, the AUC for LLU dataset is almost the same for the diffusion sequences but ADC performs slightly lower than the others. FLAIR is found to be better for sub-acute stage of stroke (as in ISLES dataset), while DWI/ADC is better for acute stroke (as in LLU dataset). This algorithm has proved its efficiency by performing better than the state-of-the-art and also 32

Figure 11: ROC plot for FLAIR, T2-WI, T1-WI, DWI sequences of ISLES dataset.The AUC of FLAIR is noticeably higher than other sequences.

estimates distributed injury with significant reliability as shown in Figure 7. This figure shows multi-foci injury estimation from multiple modalities of both the datasets considered for experimentation. 4.4. Auxiliary Analysis In the proposed SSHC method MSA is detected using multi-line Hough transform with iterative rotational correction. The asymmetric ROI is obtained by comparing the left and the right hemispheres. We compare the proposed MSA detection and the resulting asymmetric map for lesion segmentation, with two other methods described in [16], [23], [34]. Rotational correction is essential for comparing both the brain hemispheres and SSHC automatically aligns the rotated slices iteratively while determining the MSA. In the approach adopted by Rajini et al. the mid-line is detected using 3 control points joined by a Beizer curve [16]. These control points are obtained from the symmetric map by comparing 48 pixel intensities of either hemispheres. In the first method (Method 1) for comparision we adopted the MSA detection by Rajini et al.

33

Figure 12: ROC plot for DWIp, DWIr, DWIs and ADC sequences of LLU dataset. The AUC of ADC is noticeably lower than other sequences. But almost all the sequences show similar variations.

[16]. Tang et al. presented an automated mid-line detection method which uses a row of slice intensities for cross correlation with the reversed row to identify maximally uncorrelated points as MSA points [23]. Karthik et al. used a curvelet transform to classify injured and non-injured slices from CT scans [34]. Though this method does not include MSA determination using symmetry directly, the wedges of scale 4 and 5 of the curvelet transform of an image carry essential textural and intensity information of the ROI. This information can be used to identify the asymmetry map by comparing the hemispheres. In the second method for comparison (Method 2) we combined the MSA detection from [23] and obtained the asymmetric map from [34] by comparing the entropy difference of curvelet wedges of the scales 4 and 5 of both the hemispheres. The results of SSHC in terms of symmetry axis detection (MSA), asymmetric map of inter-hemispheric comparison and final injury estimation are presented in Figures 13 and 14 on the two slices with varying levels of MSA rotation. It is to be noted that Method 1 and Method 2 do not include an automated rotational correction like the iterative multi-Hough MSA detection in SSHC [16], [34], [23].

34

Table 6: Symmetry based performance comparison of the proposed SSHC method with two other methods by Rajini et al. [16] and Karthik + Tang et al. [34],[23] (Highest values are marked in bold)

Symmetry

DSC

Precision

Recall

0.45±0.29

0.56±0.21

0.65±0.26

0.39±0.20

0.41±0.19

0.45±0.22

0.704±0.27

0.68±0.14

0.85±0.01

Detection Method Method 1 (Rajini et al.) [16] Method 2 (Karthik + Tang et al.)[34],[23] Proposed SSHC

Therefore after the MSA detection in Methods 1 and 2, we realigned the MSA’s in both these methods by adopting the rotational correction part of SSHC. MSA in the proposed method is accurate in both cases owing to the combined estimation from multiple Hough lines and rotational correction. Whereas, the MSA in Method 1 is slightly deviated as it relies only on 3 control points for estimation. MSA in Method 2 follows the midline for case 1 (less rotation/inclination, see Figure 13), as the slice has lesser inclination and is nearly bilaterally symmetric. But in case 2 (high rotation/inclination, see Figure. 14) the detected MSA deviates from the actual due to self-asymmetric correlation discrepancy. There is no direct measure for parametric evaluation of inter-hemispheric asymmetry after identifying the MSA. Hence, to compare the three methods, final lesion segmentation as the result of asymmetry is considered. The same parametric measures used for evaluation of SSHC with the state-of-the-art are reconsidered for measuring the effect of accuracy of MSA and asymmetric maps. The parametric values as a result of this comparison are presented in Table 6. It can be considered that the inaccuracies in determining the MSA can greatly effect the resulting asymmetry. The lower DSC and precision values of Methods 1 and 2 can be attributed to the misaligned overlap of the hemispheres due to inaccurate MSA detection.

35

Figure 13: Symmetry comparison of the proposed with 2 other methods for a slice with Low MSA rotation: Column a) Symmetry axis (MSA) detection; Column b) Asymmetric map obtained by inter-hemispheric comparison; Column c) Asymmetric ROI detection estimated from the asymmetric map; Column d) Final lesion of SSHC, Method 1 (Rajini et al. [16]) and Method 2 (Karthik + Tang et al.[34],[23]) marked in red, compared with manual delineation marked in green.

36

Figure 14: Symmetry comparison of the proposed with 2 other methods for a slice with high MSA rotation: Column a) Symmetry axis (MSA) detection; Column b) Asymmetric map obtained by inter-hemispheric comparison; Column c) Asymmetric ROI detection estimated from the asymmetric map; Column d) Final lesion of SSHC, Method 1 (Rajini et al. [16]) and Method 2 (Karthik + Tang et al. [34],[23]) marked in red, compared with manual delineation marked in green.

37

4.5. Advantages of SSHC over the state-of-the-art The main advantages of the proposed SSHC method compared to the aforementioned methods and all others mentioned in the literature are as described: • SSHC is based on a simple yet effective proposition of inter-hemispheric asymmetry based intensity differencing for initial lesion estimation. This proposition yields results which are very similar to the gold standard manual ground truth because even radiologists intuitively delineate the injury by comparing the intensity variations within the hemispheres and also the region surrounding the lesion. • Patch based asymmetry is obtained rather than pixel to pixel comparison to determine the MSA from optimum slice. The MSA remains constant across the volume and yields results as efficient as 3D asymmetry, simultaneously making the method computationally efficient. • To calculate asymmetry determination of MSA is of utmost importance. Instead of using single Hough line for mid-sagittal axis detection, multiple Hough lines are combined which makes it least affected by the mid-line shift. The improved symmetry axis determination also includes rotational correction making it more reliable and least prone to erroneous MSA. • To further strengthen the accuracy of asymmetry determined lesion along with intensity differencing, local voxel information is combined with the asymmetric initial lesion estimate. Superpixels combine the local textural, intensity and neighbouring voxel information of the MR images, and these superpixels are hierarchically clustered into a dendogram tree. The predetermined asymmetry intensity segregates the dendogram tree nodes into lesion and NABM making SSHC completely automated without the use of hand crafted features and also eliminating the need of training. • Strong validation of SSHC is performed by evaluating it on both acute and sub-acute stroke datasets. To ensure a fair comparison of the method

38

it is validated on a bench-marked ISLES dataset for sub-acute stroke and another challenging LLU dataset for acute stroke. Parametric evaluation of SSHC is done on 54 patients with 4 sequences each and total of 20,888 slices yielding results at par and often better than machine learning and deep learning methods.

5. Conclusion Our proposed SSHC technique provides a novel segmentation method that is more accurate than the existing baseline segmentation methods presented in the ISLES’15 challenge. This automated algorithm is simple yet robust, fast and accurate, with minimal manual input. It eliminates the need of manual ad-hoc selection of threshold values both for clustering and for classifying superpixels – as stroke or NABM, by inclusion of asymmetry based information. Aggregation of asymmetry information and superpixel clustering makes it a novel approach with encouraging results in cases where stroke lesions are spread not only in one hemisphere but in the whole brain as it can also determine the hemisphere of injury by intensity comparison. It can be concluded that FLAIR and T2-WI are better for detection of sub-acute stroke (ISLES data) while DWI/ADC is better for acute stroke (LLU data). This proposed method of assessing ischemic lesion volume can be further extended and used as an approach for estimating salvageable tissue in stroke.

Acknowledgment The authors would like to thank Dr. Elia Haddad, Dr. Karen Tong, Dr. Sheri Harder, Dr. J. Paul Jacobson, Amy Plaia, Jonathan Burden and Paul S. from Departments of Radiology, Neurology and Pediatrics from the School of Medicine, Loma Linda University, for providing the clinical and imaging data, clinical assistance with generating ground truth and funding to carry out the battery of manual detection techniques.

39

References [1] N. Hoggard, I. D. Wilkinson, P. D. Griffiths, The imaging of ischaemic stroke, Clinical Radiology 56 (3) (2001) 171–183. [2] M. A. Jacobs, P. Mitsias, H. Soltanian-Zadeh, S. Santhakumar, A. Ghanei, R. Hammond, D. J. Peck, M. Chopp, S. Patel, Multiparametric MRI tissue characterization in clinical stroke with correlation to clinical outcome: part 2, Stroke 32 (4) (2001) 950–957. [3] M. Artzi, O. Aizenstein, T. Jonas-Kimchi, V. Myers, H. Hallevi, D. B. Bashat, Flair lesion segmentation: application in patients with brain tumors and acute ischemic stroke, European Journal of Radiology 82 (9) (2013) 1512–1518. [4] C. J. Ledezma, J. B. Fiebach, M. Wintermark, Modern imaging of the infarct core and the ischemic penumbra in acute stroke patients: Ct versus mri, Expert Review of Cardiovascular Therapy 7 (4) (2009) 395–403. [5] J. G. Merino, S. Warach, Imaging of acute stroke, Nature Reviews Neurology 6 (10) (2010) 560. [6] A. J. Yoo, B. Pulli, R. G. Gonzalez, Imaging-based treatment selection for intravenous and intra-arterial stroke therapies: a comprehensive review, Expert Review of Cardiovascular Therapy 9 (7) (2011) 857–876. [7] J. Weinman, G. Bissias, J. Horowitz, E. Riseman, A. Hanson, Nonlinear diffusion scale-space and fast marching level sets for segmentation of mr imagery and volume estimation of stroke lesions, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, 2003, pp. 496–504. [8] M. Prastawa, E. Bullitt, S. Ho, G. Gerig, A brain tumor segmentation framework based on outlier detection, Medical Image Analysis 8 (3) (2004) 275–283.

40

[9] K. Van Leemput, F. Maes, D. Vandermeulen, A. Colchester, P. Suetens, et al., Automated segmentation of multiple sclerosis lesions by model outlier detection, IEEE Transactions on Medical Imaging 20 (8) (2001) 677–688. [10] F. Maes, K. Van Leemput, L. E. DeLisi, D. Vandermeulen, P. Suetens, Quantification of cerebral grey and white matter asymmetry from MRI, in: International Conference on Medical Image Computing and ComputerAssisted Intervention, Springer, 1999, pp. 348–357. [11] J. F. Iaccino, Left brain-right brain differences: Inquiries, evidence, and new approaches, Psychology Press, 2014. [12] W. G. Webster, Territoriality and the evolution of brain asymmetry, Annals of the New York Academy of Sciences 299 (1) (1977) 213–221. [13] M. E. Brummer, Hough transform detection of the longitudinal fissure in tomographic head images, IEEE Transactions on Medical Imaging 10 (1) (1991) 74–81. [14] P. Minovic, S. Ishikawa, K. Kato, Symmetry identification of a 3-D object represented by octree, IEEE Transactions on Pattern Analysis and Machine Intelligence 15 (5) (1993) 507–514. [15] T. Hachaj, M. R. Ogiela, Automatic detection and lesion description in cerebral blood flow and cerebral blood volume perfusion maps, Journal of Signal Processing Systems 61 (3) (2010) 317–328. [16] N. H. Rajini, R. Bhavani, Computer aided detection of ischemic stroke using segmentation and texture features, Measurement 46 (6) (2013) 1865–1874. [17] Y. Sun, B. Bhanu, Reflection symmetry-integrated image segmentation, IEEE Transactions on Pattern Analysis and Machine Intelligence 34 (9) (2012) 1827–1841. [18] N. Ghosh, Y. Sun, B. Bhanu, S. Ashwal, A. Obenaus, Automated detection of brain abnormalities in neonatal hypoxia ischemic injury from MR images, Medical Image Analysis 18 (7) (2014) 1059–1069. 41

[19] E. T. Esfahani, D. W. McBride, S. B. Shafiei, A. Obenaus, A real-time analysis of traumatic brain injury from t2 weighted magnetic resonance images using a symmetry-based algorithm, in: Video Bioinformatics, Springer, 2015, pp. 99–117. [20] R. Peter, P. Korfiatis, D. Blezek, A. O. Beitia, I. Stepan-Buksakowska, D. Horinek, K. D. Flemming, B. J. Erickson, A quantitative symmetrybased analysis of hyperacute ischemic stroke lesions in noncontrast computed tomography, Medical Physics 44 (1) (2017) 192–199. [21] A. Vupputuri, S. Dighade, P. Prasanth, N. Ghosh, Symmetry determined superpixels for efficient lesion segmentation of ischemic stroke from MRI, in: 2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), IEEE, 2018, pp. 742–745. [22] Y. Wang, A. K. Katsaggelos, X. Wang, T. B. Parrish, A deep symmetry convnet for stroke lesion segmentation, in: 2016 IEEE International Conference on Image Processing (ICIP), IEEE, 2016, pp. 111–115. [23] F.-h. Tang, D. K. Ng, D. H. Chow, An image feature approach for computer-aided detection of ischemic stroke, Computers in Biology and Medicine 41 (7) (2011) 529–536. [24] I. Rekik, S. Allassonni`ere, T. K. Carpenter, J. M. Wardlaw, Medical image analysis methods in MR/CT-imaged acute-subacute ischemic stroke lesion: segmentation, prediction and insights into dynamic evolution simulation models. a critical appraisal, NeuroImage: Clinical 1 (1) (2012) 164–178. [25] N. Ghosh, Y. Sun, C. Turenius, B. Bhanu, A. Obenaus, S. Ashwal, Computational analysis: a bridge to translational stroke treatment, in: Translational Stroke Research, Springer, 2012, pp. 881–909. [26] J. Mitra, P. Bourgeat, J. Fripp, S. Ghose, S. Rose, O. Salvado, A. Connelly, B. Campbell, S. Palmer, G. Sharma, et al., Lesion segmentation from mul-

42

timodal MRI using random forest following ischemic stroke, NeuroImage 98 (2014) 324–335. [27] O. Maier, M. Wilms, J. von der Gablentz, U. M. Kr¨amer, T. F. M¨ unte, H. Handels, Extra tree forests for sub-acute ischemic stroke lesion segmentation in MR sequences, Journal of Neuroscience Methods 240 (2015) 89–100. [28] Y. Zheng, D. Ai, P. Zhang, Y. Gao, L. Xia, S. Du, X. Sang, J. Yang, Feature learning based random walk for liver segmentation, PloS One 11 (11) (2016) e0164098. [29] K. Kamnitsas, C. Ledig, V. F. Newcombe, J. P. Simpson, A. D. Kane, D. K. Menon, D. Rueckert, B. Glocker, Efficient multi-scale 3D CNN with fully connected CRF for accurate brain lesion segmentation, Medical Image Analysis 36 (2017) 61–78. [30] L. Chen, P. Bentley, D. Rueckert, Fully automatic acute ischemic lesion segmentation in DWI using convolutional neural networks, NeuroImage: Clinical 15 (2017) 633–643. ¨ Yıldırım, U. R. Acharya, Application of deep [31] M. Talo, U. B. Baloglu, O. transfer learning for automated brain abnormality classification using MR images, Cognitive Systems Research 54 (2019) 176–188. [32] S. Pereira, A. Pinto, V. Alves, C. A. Silva, Brain tumor segmentation using convolutional neural networks in MRI images, IEEE Transactions on Medical Imaging 35 (5) (2016) 1240–1251. [33] G. Praveen, A. Agrawal, P. Sundaram, S. Sardesai, Ischemic stroke lesion segmentation using stacked sparse autoencoder, Computers in Biology and Medicine 99 (2018) 38–52. [34] R. Karthik, R. Menaka, A multi-scale approach for detection of ischemic stroke from brain mr images using discrete curvelet transformation, Measurement 100 (2017) 223–232. 43

[35] U. R. Acharya, K. M. Meiburger, O. Faust, J. E. W. Koh, S. L. Oh, E. J. Ciaccio, A. Subudhi, V. Jahmunah, S. Sabut, Automatic detection of ischemic stroke using higher order spectra features in brain MRI images, Cognitive Systems Research (2019). [36] A. Subudhi, U. R. Acharya, M. Dash, S. Jena, S. Sabut, Automated approach for detection of ischemic stroke using Delaunay Triangulation in brain MRI images, Computers in Biology and Medicine 103 (2018) 116– 129. [37] Y. Zhang, X. Li, X. Gao, C. Zhang, A simple algorithm of superpixel segmentation with boundary constraint, IEEE Transactions on Circuits and Systems for Video Technology 27 (7) (2017) 1502–1514. [38] X. Ren, J. Malik, Learning a classification model for segmentation, in: Computer Vision(ICCV),2003 IEEE International Conference on Computer Vision, IEEE, 2003, p. 10. [39] R. Achanta, A. Shaji, K. Smith, A. Lucchi, P. Fua, S. S¨ usstrunk, SLIC superpixels compared to state-of-the-art superpixel methods, IEEE Transactions on Pattern Analysis and Machine Intelligence 34 (11) (2012) 2274– 2282. [40] Z. Tian, L. Liu, Z. Zhang, B. Fei, Superpixel-based segmentation for 3D prostate MR images, IEEE Transactions on Medical Imaging 35 (3) (2016) 791–801. [41] G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, C. I. S´anchez, A survey on deep learning in medical image analysis, Medical Image Analysis 42 (2017) 60–88. [42] H. I. Turkmen, A. Albayrak, M. E. Karsligil, I. Kocak, Superpixel-based segmentation of glottal area from videolaryngoscopy images, Journal of Electronic Imaging 26 (6) (2017) 061608. 44

[43] R. C. Gonzalez, R. E. Woods, Digital Image Processing, 2nd Edition, Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 2001. [44] N. Otsu, A threshold selection method from gray-level histograms, IEEE Transactions on Systems, Man, and Cybernetics 9 (1) (1979) 62–66. [45] A. Vupputuri, S. Ashwal, B. Tsao, E. Haddad, N. Ghosh, MRI based objective ischemic core-penumbra quantification in adult clinical stroke, in: 2017 39th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), IEEE, 2017, pp. 3012–3015. [46] O. Maier, B. H. Menze, J. von der Gablentz, L. H¨ani, M. P. Heinrich, M. Liebrand, S. Winzeck, A. Basit, P. Bentley, L. Chen, et al., ISLES 2015A public evaluation benchmark for ischemic stroke lesion segmentation from multispectral MRI, Medical Image Analysis 35 (2017) 250–269. [47] M. I. challenge, ISLES homepage, http://www.isles-challenge.org/ ISLES2015/, [Online; accessed 9-July-2019] (2015). [48] K. Kamnitsas, L. Chen, C. Ledig, D. Rueckert, B. Glocker, MultiScale 3D Convolutional Neural Networks for Lesion Segmentation in Brain MRI, http://www.isles-challenge.org/ISLES2015/articles/ kamnk1.pdf, [Online; accessed July 7, 2019] (2015). [49] C. Feng, D. Zhao, M. Huang, Segmentation of Stroke Lesions in Multi-spectral MR Images Using Bias Correction Embedded FCM and Three Phase Level Set, http://www.isles-challenge.org/ISLES2015/ articles/fengc1.pdf, [Online; accessed July 7, 2019] (2015). [50] H. Halme, A. Korvenoja, E. Salli, ISLES (SISS) challenge 2015: segmentation of stroke lesions using spatial normalization, Random Forest classifcation and contextual clustering, http://www.isles-challenge.org/ ISLES2015/articles/halmh1.pdf, [Online; accessed July 7, 2019] (2015). [51] S. M. S. Reza, L. Pei, K. M. Iftekharuddin, Ischemic Stroke Lesion Segmentation Using Local Gradient and Texture Features, http://www. 45

isles-challenge.org/ISLES2015/articles/rezas1.pdf, [Online; accessed July 7, 2019] (2015). [52] R. Zhang, L. Zhao, W. Lou, J. M. Abrigo, V. C. Mok, W. C. Chu, D. Wang, L. Shi, Automatic segmentation of acute ischemic stroke from dwi using 3d fully convolutional densenets, IEEE Transactions on Medical Imaging 37 (9) (2018) 2149–2160. [53] M. S. Babu, V. Vijayalakshmi, An effective approach for sub-acute ischemic stroke lesion segmentation by adopting meta-heuristics feature selection technique along with hybrid naive bayes and sample-weighted random forest classification, Sensing and Imaging 20 (1) (2019) 7. [54] N. K. Subbanna, D. Rajasheka, B. Cheng, G. Thomalla, J. Fiehler, T. Arbel, N. D. Forkert, Stroke lesion segmentation in FLAIR MRI datasets using customized Markov random fields, Frontiers in Neurology 10 (2019) 541. [55] A. Gautam, B. Raman, Segmentation of ischemic stroke lesion from 3d mr images using random forest, Multimedia Tools and Applications 78 (6) (2019) 6559–6579. [56] R. Karthik, U. Gupta, A. Jha, R. Rajalakshmi, R. Menaka, A deep supervised approach for ischemic lesion segmentation from multimodal mri using fully convolutional network, Applied Soft Computing 84 (2019) 105685. [57] P. Berkhin, A survey of clustering data mining techniques, in: Grouping Multidimensional Data, Springer, 2006, pp. 25–71. [58] P. P. Mohanta, D. P. Mukherjee, S. T. Acton, Agglomerative clustering for image segmentation, in: Object Recognition Supported by User Interaction for Service Robots, Vol. 1, IEEE, 2002, pp. 664–667. [59] R. Ortiz-Ram´ on, M. d. C. V. Hern´andez, V. Gonz´alez-Castro, S. Makin, P. A. Armitage, B. S. Aribisala, M. E. Bastin, I. J. Deary, J. M. Wardlaw,

46

D. Moratal, Identification of the presence of ischaemic stroke lesions by means of texture analysis on brain magnetic resonance images, Computerized Medical Imaging and Graphics 74 (2019) 12–24.

47

Highlights of Symmetry determined Superpixel based Hierarchical Clustering (SSHC) • • • •

Automated delineation of ischemic injury in acute and sub-acute stroke. Superpixel based clustering of injury leveraged on asymmetry cue. Validation on two different datasets with mutisequence MRI data. Eliminates the need of high-end computations and training.

Conflict of Interest statement

None declared