Development of a cerebral aneurysm segmentation method to prevent sentinel hemorrhage

Image segmentation being the first step is always crucial for brain aneurysm treatment planning; it is also crucial during the procedure. A robust brain aneurysm segmentation has the potential to prevent the blood leakage, also known as sentinel hemorrhage. Here, we present a method combining a multiresolution and a statistical approach in two dimensional domain to segment cerebral aneurysm in which the Contourlet transform (CT) extracts the image features, while the Hidden Markov Random Field with Expectation Maximization (HMRF-EM) segments the image, based on the spatial contextual constraints. The proposed algorithm is tested on Three-Dimensional Rotational Angiography (3DRA) datasets; the average values of segmentation accuracy, DSC, FPR, FNR, specificity, and sensitivity, are found to be 99.72%, 93.52%, 0.07%, 5.23%, 94.77%, and 99.96%, respectively.


Background
The sentinel clot sign is a useful computed tomography finding for the evaluation of probable anatomic sites of hemorrhage.Sentinel leaks produce sudden focal or generalized head pain that may be severe.Sentinel headaches precede the aneurysm rupture.To be further specific, cerebral Aneurysm (CA) is an abnormal inflammation within the brain blood vessels that is commonly detected after the age around forty (Higashida 2003;Nikravanshalmani et al. 2013).The mortality rate of this disease is between 30% and 40%, whereas the morbidity rate ranges between 35% and 60% (Li and Li 2010).Computer Aided Diagnosis (CAD) algorithms have been helpful in reducing the bias, time consumption, increasing the CA detection accuracy, etc. Wong (2005).Medical image segmentation (MIS) algorithms are important in the diagnosis phase.
There have been several methods in image segmentation that include both traditional methods and machinelearning/deep-learning methods.In deep learning (DL)based methods, the neural networks (NN) used for image segmentation, have received much attention due to their end-to-end nature and state-of-the-art performance.Deep neural networks gained more popularity after their success at ImageNet Large Scale Visual Recognition Challenge 2012 (Russakovsky et al. 2015).With the improvement in GPU architectures and the development of dedicated deep learning libraries [e.g., Tensorflow (Abadi et al. 2016) and PyTorch (Paszke et al. 2019)], neural networks have been in use for a wide range of computationally heavy computer vision tasks (e.g., object recognition and segmentation).For image segmentation, the basic schematic diagram of a convolutional neural networks (CNNs) is given in Fig. 1; it comprises a series of convolution and pooling layers to condense the input image into a dense feature representation.The feature representation is utilized to reconstruct the segmentation mask using deconvolution Abbes Amira and Sarada Prasad Dakua contributed equally to this work.

Related work
There is rich literature in conventional MIS that includes semi-automatic and automatic (Dakua 2013;Dakua and Sahambi 2011).Broadly, the semi-automatic approaches include threshold-based (Mitra and Chandra 2013), modelbased (Dakua et al. 2019), graph-based (Dakua andSahambi 2010, 2011) methods that need human intervention (Chen et al. 2014), which is tedious and prone to operator variability; these factors certainly affect the segmentation accuracy.Therefore, automatic approaches were introduced.In Mitra and Chandra (2013), spatial filtering and dual global thresholding are considered, where three parameters (i.e. two threshold values and filter mask size) are user-defined.Yang et al. present a dot enhancement filter that includes thresholding, and a region growing technique.However, it lacks the inconsistency of the sensitivity, where it ranges between 80% and 95% for CAs larger than 5 mm and between 71% and 91% otherwise.Bogunovic et al. (2011) present a Geodesic Active Regions (GAR) based algorithm.However, there is the possibility of merging either a vessel with a CA or two vessels together.This partly happens either because of the insufficiency of the imaging resolution and the low blood flow in the small vessels.In Hentschke et al. (2011), a blobness filter and a k-means technique are adapted.The algorithm results in a false positive rate reaching 20.8%.Jerman et al. (2015), present a blobness enhancement filter combined with Random Forest Decision (RFD) classifier and a grow-cut technique, where it is assumed that a saccular CA consists of more than 15 voxels.Thus, this assumption reduces the chances to detect small CAs.Suniaga et al. (2012), propose a fuzzy logic in a level set framework along with a SVM classifier to segment saccular CAs.In Zhang et al. (2012), a method on Conditional Random Field (CRF) and gentle Adaboost classifier is proposed that is trained on six datasets and tested on two datasets.However, machine Learning (specifically, Deep Learning (DL)) models need extensive data to capture the underlying distribution for generalization in real-world systems.Additionally, the black box nature of the DL models adds uncertainty to their predictions, thus limiting their use as a second opinion in clinical setting (Su et al. 2019).Lu et al. (2016) present a method, where multi-scale filtering is used to enhance the vessels suppressing the noise and a mixture model is built to fit the histogram curve of the preprocessed data.Finally, expectation maximization is used for parameter estimation.Kai et al. (Lawonn et al. 2019) present a graph-cut based method for aneurysm segmentation.Yu et al. (Chen et al. 2020) present a geodesic active contour (GAC) and Euler's elastica model-based method for aneurysm segmentation, where GAC segments the giant aneurysms in high contrast regions and the elastica model estimates the missing boundaries in low contrast regions.Thus, the numerous automatic segmentation algorithms reported so far in the literature have potential limitations that warrant further research on automatic methods.However, the automatic methods are difficult to be controlled, as usually desired by the clinicians during their treatment planning to explore various options.Thus, this paper proposes a semi-automatic algorithm combining Contourlet Transform (CT) and Hidden Markov Random Field with Expectation Maximization (HMRF-EM) to segment CA regardless the CA shape or size.
The rest of this paper is organized as follows: in Sect.3, the foundation and mathematical background for the proposed algorithm are presented.In Sect.4, the proposed algorithm is discussed in detail.In Sect.5, the objective and subjective evaluation of the proposed work are presented along with the dataset description and the environmental setup for the implementation.Section 6 includes discussion and some future work, whereas the Sect.7 concludes the paper.

Preprocessing
Noise is a prime obstacle in image segmentation and it is inherent in medical data.Thus, noise reduction is necessary to improve the segmentation outcome.Traditionally, filtering has been commonly used, however, it results in significant degradation in the input image quality.Stochastic resonance (SR) concept is an emergent one that constructively use the noise present in input data (Dakua et al. 2019).Stochastic resonance occurs when the ratio of signal to noise ration attains the maximum.At this state, the input/output correlation is maximum at a certain noise level.To exhibit SR, a system should possess three basic properties: a non-linearity in terms of threshold, sub-threshold signals like signals with small amplitude and a source of additive noise.This phenomenon occurs frequently in bistable systems or in systems with threshold-like behavior.The general behavior of SR mechanism shows that at lower noise intensities the weak signal is unable to cross the threshold, thus giving a very low SNR.For large noise intensities the output is dominated by the noise, also leading to a low SNR.But, for moderate noise intensities, the noise allows the signal to cross the threshold giving maximum SNR at some optimum additive noise level.Thus, a plot of SNR as a function of noise intensity shows a peak at an optimum noise level (shown in Fig. 2).There have already been many attempts to use SR in different domains such Fourier and spatial domains; however, we have chosen the maximal overlap discrete wavelet transform (MODWT) because of some of its key advantages: (1) MODWT can handle any sample size, (2) the smooth and detail coefficients of MODWT multiresolution analysis are associated with zero phase filters, (3) it is transform invariant, and (4) it produces a more asymptotically efficient wavelet variance estimator than DWT.

Denoising by MODWT
In this methodology, 2D MODWT is applied to the M × N size image I. Applying SR to the approximation and detail coefficients, the stochastically enhanced (tuned) coefficient sets in MODWT domain are obtained as W s (l, p, q) SR and W l 0 , p, q SR .The SR in discrete form is defined as: , where √ D (t) and B sin t represent noise and input, respectively; these are replaced by MODWT sub-band coefficients.The noise term is the factor to produce SR; maximization of SNR occurs at the double well parameter a. Implementation of SR on digital images necessitates the need for solving the stochastic differential equation using Euler-Maruyama's method (Dakua et al. 2016) that gives the iterative discrete equation: where a and e are the bistable parameters, whereas n and Δt represent iteration and sampling time, respectively.Input denotes the sequence of input signal and noise, with the initial condition being x(0) = 0 .The final stochastic simula- tion is obtained after some predefined number of iterations.Given the tuned (enhanced and stabilized) set of wavelet coefficients ( X l 0 , p, q and X s (l, p, q) ), the denoised image I denoised in spatial domain is obtained by inverse maximal overlap discrete wavelet transform (IMODWT) as: The double well parameters a and e are determined from the SNR by differentiating SR with respect to a and equating to zero; in this way, SNR is maximized resulting in a = 2 2 0 for maximum SNR, where 0 is the noise level administered to the input image.The maximum possible value of restoring force ( R = B sin t ) in terms of gradient o f s o m e b i s t a b l e p o t e n t i a l f u n c a 3e .R at this value gives maximum force as 27e .Maximizing the left term (keeping B = 1 ), e < 4a 3 27 .To get the parameter values, we consider (1) 27 ; w and z are weight parameters for a and e.Initially, w is an experimentally chosen constant that later becomes input image standard deviation dependent, while z is a number less than 1 to ensure subthreshold condition of the signal.In this way, the noise in input image is countered and maximum information from the image is achieved.This enhanced image is used subsequently.

Contourlet transform
Multiresolution analysis techniques usually utilize the image features for computer vision.There have been several multiresolution analysis techniques such as Wavelet, Ridgelet, Curvelet, and Contourlet transforms for medical image segmentation (AlZubi et al. 2011;Po and Do 2006;Moayedi et al. 2010).We have preferred the Contourlet Transform (CT) in our proposed approach due to some of its advantages (Do and Vetterli 2003).CT being an ideal 2D transform in the discrete domain has other salient features: multiresolution, localization, critical sampling, anisotropy, and multiple directions for different resolutions.In addition, this transformation provides a sparse representation saving a significant amount of memory by offering simple and fast data processing as it requires O(N) operations for an image with N-pixels (Po and Do 2006).These characteristics capture the geometrical smoothness of the 2D contours.
The Pyramidal Directional Filter Bank (PDFB) (Do and Vetterli 2001) combines the Laplacian Pyramid (LP) and Directional Filter Bank (DFB) to extract the desirable .Subsequently, the lowpass image, a 1 [n] , is passed again through LP to repeat the same process until a certain predefined number of decomposition levels, L j , is reached.Figure 3 illustrates the CT process to decom- pose a 512 × 512 image into two levels, where eight and four directions are applied at each level, respectively.
CT has the adeptness at capturing geometrical smoothness of 2D contours and anisotropy in the discrete domain.In addition, it has a high degree of directionality as it allows to define different directions for different scales, which is not possible in other multiresolution analysis techniques.These advantages help extract the features from images resulting in improved segmentation.
To summarize, the CT takes a 2D image, a 0 [n] , as an input and decomposes it into a lowpass sub-band, a j [n] , and several bandpass directional sub-bands, c L j j,k [n] , which are called as the Contourlet coefficients.This process can be expressed by: where (L)  j,k,n is the LP basis function for scale decomposition composed of the lowpass filter, g k [n] , and the scaling func- tion, j,k (t) ; and (L) j,k,n is the DFB basis function for directional decomposition composed of orthogonal filter, d k [n] , and the directional function, j,n (t) .The parameters j, k, d, and n, used in Eqs. ( 2) and (3), are number of levels/scales, number of directions for each level, dimensionality (in our case, it is 2 since we are working in 2D domain), and a scale parameter along the frequency axis, respectively.The overall decomposition process of the CT is provided in Algorithm 1.
L j ← initialize the number of decomposition levels K ← initialize the number of directions for each The number of the decomposition levels is experimentally selected since it depends on the image size and the amount of the details to be highlighted.For the number of directions, it is recommended to gradually increase them by 2 k . (2)

Hidden Markov random field model
HMRF model is a statistical approach in stochastic domain that provides prior knowledge to simplify the segmentation task.
HMRF model segments the medical images based on the spatial correlation between the neighboring pixels.Some important notions about this model are: • Random field: The random variables are the intensity levels in an image.In HMRF model, two random fields exist: -Hidden random field: X = (x 1 , x 2 , .., x N )‖x i ∈ L, i ∈ S is a random field in a finite state space, L, and indexed by a set, S, with respect to a neighboring system of size, N. The state of this field X is unobservable/hidden and every x i is independent of all other x j .The objective of this assumption is to classify each pixel independently (Marroquin et al. 2003).-Observable random field: Y = {y = (y 1 , y 2 , .., y N )‖y i ∈ D, i ∈ S} is a random field in a finite space, D, and indexed by a set, S, with respect to a neighboring system of size, N.This random field, Y, is observable and can only be defined with respect to X, where y i follows a con- ditional probability distribution given any particular configuration x i = l : p(y i ‖l) = {f (y i ; l ), ∀l ∈ L} , where l is the set of the parameters that are involved.
• Parameters: The set of parameters, l , are unknown.Therefore, a model fitting solution is adopted to estimate them.In our context, the parameters are mainly the mean, , and the standard deviation, .
• Conditional independence: The two random fields, (X, Y), are conditionally independent.• Neighborhood system: It is a way to define the surrounding pixels for a specific pixel (Chen et al. 2011).• Clique: It is a subset of pixels, where every pair of distinct pixels are neighbors.A value is assigned to each clique, c, to define the clique potential V c (x) , where the sum of all of these values results in the energy function, U(x), that we aim to minimize.
To know l , EM algorithm is preferred.The HMRF-EM framework (Zhang et al. 2001) incorporates the EM algorithm in association with HMRF model to estimate the parameters and segment using iterative updates.The framework starts by initializing both: the segmentation and parameters (means, , and standard deviations, ).Then, iteratively, it goes through the Expectation Step (E-Step) and Maximization Step (M-Step) to update the parameters and the initial segmentation until no further development is observed.The E-Step updates the segmentation by assigning to each pixel an estimated class label, x , from a set of labels, L. The assignment is performed based on the MAP criterion, which tries to minimize the posterior energy using the current parameters estimate; during the energy maximization, the conditional posterior probability distribution, P(Y‖X) , gets maximized.Equation ( 5) illustrates the energy calculation: where U(x) is the energy function as illustrated in Eq. ( 4) and U(y‖x) is the likelihood energy illustrated below in Eq. ( 6).
While the M-Step updates the parameters based on ML criterion, which tries to maximize the expected likelihood found in the E-Step.The parameters and are calculated using the Eqs.( 7) and ( 8), respectively.
[θ (0) ,x (0) ] ← initialize the parameters and segmentation EM itr ← initialize the number of EM iterations MAP itr ← initialize the number of MAP iterations i ∈ {1, ..., EM itr} j ∈ {1, ..., MAP itr} x(j) ← Update the segmentation x(j−1) based on the MAP criterion θ (i) ← Update the parameters θ (i−1) based on ML criterion θ (EM itr) , x(MAP itr) This framework works well for small data dimensions; the main advantages are, it: is easy to implement, provides an accurate segmentation, and is less sensitive to noise as compared to other segmentation techniques since it also considers the contextual information (Yazdani et al. 2015).

Modifications to conventional k-means clustering
This is well-known that k-means is a clustering technique maximizing the similarity of intra-class and minimizing the similarity of inter-class; it is computationally fast.In our method, we initialize k automatically based on the image entropy, a statistical measure of randomness that can be used Page 7 of 14 18 to characterize the texture of the gray-scale image (Liang et al. 2012).It is expressed in Eq. ( 9): where X is a vector of all intensities of an image and n is the number of pixels.

Proposed segmentation algorithm
The proposed CA segmentation algorithm is illustrated in Fig. 4. The algorithm starts by feeding the input 2D scans after enhancement.The selection of the Region of Interest (ROI) from the entire cerebral vasculature is done manually.Subsequently, the below steps are performed on each 2D image.
During the first step, CT is applied to extract features from the input image by decomposing it into six pyramidal levels and different directions for each level, where the number of directional decomposition at each pyramidal level (from coarse to fine) are: 2 2 , 2 2 , 4 2 , 4 2 , 8 2 , and 8 2 (Hiremath et al. 2010;Po and Do 2006).As discussed earlier, CT consists of two main filters, LP and DFB.A ladder filter, known as PKVA filter (Phoong et al. 1995), is selected for the first filter.The PKVA filter is effective to localize edge direction as it reduces the inter-direction mutual information (Po and Do 2006).As for the second filter, 9-7 bi-orthogonal Daubechies filter (Cohen and Daubechies 1993) this filter significantly reduces the inter-scale, inter-location, and inter-direction mutual information of the Contourlet (Po and Do 2006).After the decomposition by CT, the lowpass subband image, a j [n] , is used to perform the rest of the steps.
To apply the second step (HMRF-EM) of the segmentation algorithm, two prior steps need to be performed.First, an initial segmentation is generated using the k-means method, which is known with its under-estimation to complement the HMRF-EM framework (AlZubi et al. 2011).In addition, a Canny edge filter (Gonzalez and Wood 2002) is applied to highlight the image edges.After that, the HMRF-EM algorithm iterates between the E-Step and M-Step to enhance the initial segmentation constrained by the Canny edge detector.
The HMRF-EM method starts after getting the initial segmentation, x( 0 After completing these two phases, the reconstruction of all the segmented 2D images is performed to get the final segmented 3D volume of the ROI.The pseudocode for the overall proposed CA segmentation algorithm is presented in Algorithm 3. 5 Datasets and results

Datasets
The aneurysms obtained from the respective 3DRA scanner (Siemens machine) were small ( ≤ 3 mm).Six 3DRA data- sets are provided by Hamad Medical Cooperation (HMC) to validate the proposed segmentation; each dataset consists of 385 2D slices of size 512 × 512 each.In addition, the cor- responding ground truth datasets were obtained from the same hospital that were prepared by three neuro-experts.

Results
For preprocessing/denoisng of the input images, the initial values of Δt and z are taken as 0.007 and 0.000027, respectively.The results are provided in Fig. 5.To determine the quality of the denoised image, we have calculated distribution separation measure (DSM) that estimates the , where E T and O T are the mean of the selected target regions of the denoised and original images, respectively; E B and O B are the mean of the selected background region of the denoised and original image, respectively.Higher the value of DSM, better is the quality.It is observed that the value of DSM is maximum at iteration 200 and then it starts decreasing, therefore, this iteration is considered as the optimal.These denoised data are used subsequently as the input volumetric data.We have now two scenarios: (1) to segment the region of interest, and (2) the entire input volumetric data. in scenario 1, we manually, select the slices to get the ROI and it is, on average, 86 slices and they are continuous, whereas in scenario 2, we feed all the slices as input to the algorithm.Both the quantitative and qualitative results are obtained by comparing the segmented volumes with the respective ground truth volume, as elaborated in Sects.5.2.2 and 5.2.3, respectively.Figures 6, 7, 8, and 9 depict each dataset before and after applying the segmentation.To test, if the results are statistically significant, we set the significance level, = 0.05 and we obtained the p-values smaller than 0.05, indicating that these segmentation results correctly identified the brain aneurysm.

Registration
We have also considered image registration in our work to enable the comparison between the segmented volume and the corresponding ground truth.In this process, one of the images is defined as the target (or the subject), which we wish to apply a transformation, while the other image is defined as the reference (or the source).In our case, the target image is the segmented ROI volume; while the reference image is the ground truth ROI volume.The target image is transformed by means of the selected mapping functions to align it with the reference image (Zitova and Flusser 2003).We have selected an affine transformation; Fig. 10c illustrates the coordinate systems of the segmented volume and the ground truth data before and after registration.

Objective evaluation
Six performance metrics are used to measure the proposed segmentation quantitatively: Dice Similarity Index (DSI), sensitivity, specificity, accuracy, False Positive Ratio (FPR), and False Negative Ratio (FNR).The value of these metrics ranges between 0 and 1 (but we have determined their percentage value).Table 1 provides the definition and the formula of each metric.
Two measures, True Positive (TP) and True Negative(TN), in Table 1 indicate  The values of these performance metrics are presented in Table 2 against each dataset.

Subjective evaluation
Each segmentation has been assessed visually by five observers (neurologists).The score, ranging between 0 and 5, is assigned by each observer, where 5 means that the ground truth volume almost completely matches with segmented one and 0 means that the two volumes do not match.Table 3 presents the observations.

Discussion
Both the qualitative and quantitative results are promising when tested on these 3DRA datasets.In the quantitative evaluation, the average values of accuracy, DSC, FPR, FNR, specificity, and sensitivity, are 99.72% , 93.52% , 0.07% , 5.23% , 94.77% , and 99.96% , respectively; an aver- age of 4.14 over 5 is obtained in the qualitative evaluation.
It may be observed in Tables 2 and 3 that the last dataset has the worst results as compared to the remaining datasets in both the quantitative and qualitative evaluation.This may be due to the fact that the provided ground truth data does not involve the complete brain vessels tree and only a delineated ROI is provided, where some surrounding vessels are excluded.The same may be observed from Fig. 12.
The computational time to segment a CA volume is considerably fast.Table 4 reports the running time of the proposed segmentation algorithm for the ROI in seconds (s) and the whole volume of a subject in minutes (min).
Even though the results are acceptable, the computational time can further be reduced by using Field-Programmable Gate Array (FPGAs) or Graphics Processing Unit (GPUs).
We have also compared the performance of the proposed method with some similar methods that have been published recently.The results are provided in Table 5.The results show that the performance of the proposed method is better than the others.
Additionally, we have compared (in Table 6) the proposed method with some popular DL-based methods, including Voxel-Morph (VM) (Balakrishnan et al. 2019), LT-Net (Wang et al. 2020), and Symmetric Normalization (SyN) (Avants et al. 2008).We have selected these methods because they have been regularly preferred by the research fraternity for comparison purpose.Among these, Voxel-Morph is probably the most famous method in recent years.The results show that the proposed method fairly performs as compared to the DL-based methods although the margin is not very significant.Furthermore,   data to improve the contrast of the input images and make them suitable for segmentation.We have obtained promising results when tested on 3DRA datasets.Since, it is usually difficult to obtain publicly available original and corresponding ground truth data, we had rely on HMC data.The preparation of ground truth data needs a huge effort and thus, we have used only a small number of datasets.In future, we would like to increase the number of datasets and make the segmentation algorithm further robust.Furthermore, since this segmentation algorithm could be used to improve the cerebral aneurysm treatment planning, we also plan to take this into clinical routine after thorough testing and validation.
Acknowledgements Not applicable.we have also compared the computational complexity involved in Table 7.

Conclusions
Sentinel headache is often assumed to portend an increased risk of delayed cerebral ischemia and aneurysm rebleeding.High rates of morbidity and mortality are associated with Sub-Arachnoid Hemorrhage (SAH) that is caused by a ruptured CA.Therefore, it is quite imperative to detect and diagnose CAs at an early stage.In this paper, we have presented a semi-automatic CA segmentation method using Contourlet transform and the hidden Markov random field model.Prior to use the method, we have denoised the input

Fig. 1
Fig. 1 Deep learning-based system for image segmentation

Fig. 2 a
Fig. 2 a SNR vs noise standard deviation curve.b Bistable double potential well

Fig. 3
Fig. 3 Contourlet transform decomposition process of a 512 × 512 image into two levels, where eight and four directions are applied at each level, respectively Fig. 4 Flowchart of the proposed segmentation

Algorithm 3
Proposed algorithm for cerebral aneurysm segmentation.Read the input images of a dataset and store them in V1 Select ROI from V1 and store them in V2 each img ∈ V 2 Apply CT to decompose img and extract the coarsest coefficients Initialize the number of clusters k based on the image entropy Apply k-means on the coarsest coefficients Apply canny edge operator on the coarsest coefficients Apply HMRF-EM algorithm to get the final segmentation Reconstruct the image by applying the ICT Reconstruct the 3D segmented volume, V3

Fig. 5
Fig. 5 a, b Input angiographic images with object of interest being highlighted.Contrast enhancement of angiographic images at: c, d iteration number 200, and e, f iteration number 300

Fig. 6
Fig. 6 Three slices from dataset 1 a before segmentation, b after segmentation, c original ROI, d segmented ROI a correct segmentation, while False Positive (FP) and False Negative (FN) indicate an incorrect segmentation.Figure 11 depicts the meaning of each measure more clearly.

Fig. 9
Fig. 9 Three slices from dataset 4 a before segmentation, b after segmentation, c original ROI, d segmented ROI

Fig. 11
Fig. 11 Four measures used in performance metrics

Fig. 12
Fig. 12 Dataset 4 a Segmented volume from the original DICOM dataset.b Ground truth volume delineated by the experts

Table 4
Time consumption to segment CA using the proposed approach