1 Introduction
Vascular diseases are one of the major reasons of deaths and disability for human health in the world (Rothwell et al., 2003; Suri and Laxminarayan, 2003). Visualization of vessels is a fundamental part of the early detection and diagnosis of such vascular disease (Dougherty, 2011). This aids surgeons, radiologists and oncology specialists in the diagnosis of abnormalities and surgical planning. Vessels serve as landmarks or road maps, before and during the surgery, and also help decision-making in the operating room in real time and postoperative monitoring. Analysis of vessels is very challenging due to complexity, variety of shape, branches, densities, small diameters and dynamic range of intensity vessels. There exist several imaging modalities consisting of magnetic resonance angiography (MRA), digital subtraction angiography, positron emission tomography and computed tomography angiography (Suri and Laxminarayan, 2003). Images generated by current imaging modalities are often unsatisfactory because of the presence of noise, artifacts, low intensity and the complex structure of vessels. Hence, there is a need to the accurate vessel extraction algorithms to overcome the limitations. Automatic or semi-automatic vessel extraction aids the clinician in making an accurate diagnosis and grading of the stenoses and aneurysms in vessels (Suri and Laxminarayan, 2003).
The purpose of vessel extraction of MRA images is to segment the image into parts of the vessel and the background. In fact, a vessel extraction method transforms given images into a binary image of zero and one. Then the boundaries of the vessels are pixels that take the value between zero and one.
The researchers are trying to use computer vision techniques to do vessel extraction, and have recently been much more interested in using the automatic or semi-automatic vessel extraction from the MRA data set. Extracting vessels from the medical imaging modalities has existed for more than 45 years, but computer-assisted extraction has begun in the past 25 years (Suri and Laxminarayan, 2003). A vessel extraction on angiography began in 1985, when the digital subtraction angiography began (Gerig et al., 1990; Iwasaki et al., 1985). Many approaches exist for vessel extraction. Cline (2000) introduces a mathematical morphology-based approach on the nonlinear mathematical operators. The Fuzzy method is used by the Fuzzy connectivity-based technique for extracting the vessels from MRA images (Saha et al., 2000; Udupa et al., 1997; Udupa and Samarasekera, 1996). Prinet et al. (1996) use a geometric differential to do vessel segmentation. In this approach, MRA images are treated as hyper surfaces. Centrelines of the vessel are obtained by linking the crest points, which are the extreme of curvature on the hyper surface. Multiscale filtering has been suggested for medical images segmentation by convolving the image with Gaussian filters (Frangi et al., 1998; Krissian et al., 1998; Lorenz et al., 1997; Sato et al., 1998). The directional anisotropic diffusion method has been suggested by Krissian et al. (1997) for vessel extraction, which uses an anisotropic diffusion to reduce noise without removing small vessels. Caselles et al. (1993) and Malladi et al. (1995) use propagating interfaces under a curvature dependent speed function to model anatomical shapes. Kirbas and Quek (2004) provide a further review on vessel segmentation.
Recently, papers have appeared that analyse vessel extraction for various problems. Kakileti and Venkataramani (2016) present an automated algorithm for detection of blood vessels in 2D-thermographic images for breast cancer screening. Navid et al. (2020) introduce a novel method to infrared thermal images vessel extraction based on fractal dimension. The retinal vessel segmentation has become an attractive subject. In Budak et al. (2020), a densely connected and concatenated multi encoder-decoder is proposed for segmentation of retinal vessels in colour fundus images. An effective image features a combination of supervised and unsupervised machine learning methods that are used for retinal blood vessel extraction (Hashemzadeh and Azar, 2019). This method first extracts the thick and clear vessels in an unsupervised manner, and then, it extracts the thin vessels in a supervised way. The proposed methods in Mustafa et al. (2014) utilized morphological operation for Diabetic Retinopathy.
Kirbas and Quek (2004) provide a further review on vessel segmentation. In the following, we concentrate on two approaches that relate to the presented method: the EMS and TFA algorithms.
Wells et al. (1996) introduce statistical method EMS for segmenting a data set to arbitrary classes. This method proposes a mixture model whose parameters can be estimated by using a modified expectation-maximization (EM) algorithm (Dempster et al., 1977). In addition to the above methods, new methods based on the concept of a tight-frame have been introduced in Arivazhagan and Ganesan (2003), Unser (1995) to image segmentation. The tight-frame method is a very useful tool for many different image processing applications. Recently, Cai et al. (2013) proposed TFA algorithm based on tight-frame for vessel extraction. The TFA algorithm iteratively purifies a area that surrounds the possible boundary of the vessels, or in other words, iteratively updates an interval of potential boundary pixels. In each iteration, they use tight-frame transformation to denoise and smooth the possible boundary. This algorithm automatically can segment twisted, convoluted and occluded structures (Cai et al., 2013). But the obtained results demonstrate the better performance of the proposed method compared to TFA algorithm in MRA images, meaning that vessel intensities are close to the intensity of the background.
There are various tight-frame systems. Some of these tight frames are shearlets (Guo and Labate, 2007; Labate et al., 2005), framelets (Ron and Shen, 1997), contourlets (Do and Vetterli, 2005) and curvelets (Candès et al., 2006), etc. Cigaroudy and Aghazadeh (2017), Aghazadeh and Cigaroudy (2014) propose an iterative procedure for tubular structure segmentation of 2D images based on tight frame of curvelet transform. The dual-tree complex wavelet transform (ℂWT) was introduced by Kingsbury and their colleagues (Kingsbury, 2001; Selesnick et al., 2005). It has additional properties: it is nearly shift invariant and high directionally selective. The 2D dual-tree complex wavelet transform is nonseparable, but is based on the separable filter bank (Kingsbury, 2001; Selesnick et al., 2005).
In this paper, we derive a vessel extraction algorithm that uses TFA and EMS algorithms. We use the output of the EMS algorithm to construct a new image in which vessel pixels are brightened and noise pixels are darkened. Then we use 2D dual tree complex wavelet tight-frame for denoising of the new image by determining transform matrix. In the wavelet shrinkage procedure, the nonlinear soft thresholding transform is used. We set parameters for the brightness increasing. The presented method can segment complexity structures; it can follow the branching of vessels, from thinner to larger structures; it can remove more artifacts. Also, the presented method extracts well vessels where their intensity is closer to the background. Moreover, we prove that the presented method converges to a binary image. For more comparison, we use B-spline and complex wavelet tight frames. Comparison of methods in Cai et al. (2013), Wilson and Noble, (1997, 1999) on real 2D MRA images show that our EMCTFA method gives more accurate vessel extraction. EMS, STFA and CTFA algorithms exhibit many artifacts which are well removed by the presented EMCTFA algorithm, and our method needs few iterations, unlike TFA algorithm. Numerical experiments demonstrate that when the presented method is used, just after two iterations more than 90% of the pixels are segmented.
2 Methods
2.1 2-D Dual-Tree ℂWT
The concept of a frame, originally defined by Duffin and Schaeffer (1952) and later revived by Daubechies and Bates (1993), guarantees stability while allowing non-unique decompositions. In fact, a sequence ${({\phi _{i}})_{i\in I}}$ in $\mathcal{H}$ (Hilbert space) is a tight frame for $\mathcal{H}$ if $\mathcal{F}$ has the perfect reconstruction property: ${\mathcal{F}^{\ast }}\mathcal{F}={I_{\mathcal{H}}}$, where operators $\mathcal{F}$ and ${\mathcal{F}^{\ast }}$ the adjoint of F defined as Kutyniok and Labate (2012).
\[ \mathcal{F}:\mathcal{H}\to {\ell ^{2}}(I),\hspace{2em}x\longmapsto {\big(\langle x,{\phi _{i}}\rangle \big)_{i\in I}}.\]
So for each vector $x\in \mathcal{H}$
There are various tight frame systems. Some of these tight frames are wavelets for Hilbert space ${L^{2}}(\mathbb{R})$:
\[ \big\{{\phi _{m}}=\phi (.-m):m\in \mathbb{Z}\big\}\cup \big\{{\psi _{j,m}}=\psi \big({2^{j}}.-m\big):j,m\in \mathbb{Z}\big\},\]
where $\phi \in {L^{2}}(\mathbb{R})$ and $\psi \in {L^{2}}(\mathbb{R})$ are called scale and wavelet functions, respectively. The translation of the scaling function takes care of the low-frequency region and the wavelet terms of the high-frequency region (Kutyniok and Labate, 2012). In order to apply the wavelet transform to 2D signals (images, for instance), we need to use the DWTs extension to two dimensions, namely the 2D DWT. Wavelet basis of Hilbert space ${L^{2}}({\mathbb{R}^{2}})$ is obtained from products of scaling and wavelet functions, being associated to a one-dimensional wavelet (Firoiu, 2010; Mallat, 1999). By applying the 2D separable wavelet to an image, we can only observe the vertical details, the horizontal details and the diagonal details. This means that the separable real DWT has poor directional selectivity. The complex wavelet transforms (ℂWT) were introduced in an attempt to overcome this limitation and other limitations of the 2D separable wavelet (Kingsbury, 2001).Let complex wavelet $\phi (t)={\phi _{\xi }}(t)+j{\phi _{\eta }}(t)$ and $\psi (t)={\psi _{\xi }}(t)+j{\psi _{\eta }}(t)$. By taking the real part of each of these six complex wavelets $\psi (s)\psi (t)$, $\psi (s)\overline{\psi (t)}$, $\phi (s)\psi (t)$, $\psi (s)\phi (t)$, $\phi (s)\overline{\psi (t)}$ and $\psi (s)\overline{\phi (t)}$, we obtain real oriented 2D wavelets. Specifically, we obtain the following six wavelets:
\[\begin{aligned}{}& {\psi _{i}}(s,t)=\frac{1}{\sqrt{2}}\big({\psi _{1,i}}(s,t)-{\psi _{2,i}}(s,t)\big),\\ {} & {\psi _{i+3}}(s,t)=\frac{1}{\sqrt{2}}\big({\psi _{1,i}}(s,t)+{\psi _{2,i}}(s,t)\big),\end{aligned}\]
for $i=1,2,3$, where the two separable 2D wavelet bases are defined by:
\[\begin{aligned}{}& {\psi _{1,1}}(s,t)={\phi _{\xi }}(s){\psi _{\xi }}(t),\hspace{2em}{\psi _{1,1}}(s,t)={\phi _{\eta }}(s){\psi _{\eta }}(t),\\ {} & {\psi _{1,1}}(s,t)={\psi _{\xi }}(s){\phi _{\xi }}(t),\hspace{2em}{\psi _{1,1}}(s,t)={\psi _{\eta }}(s){\phi _{\eta }}(t),\\ {} & {\psi _{1,1}}(s,t)={\psi _{\xi }}(s){\psi _{\xi }}(t),\hspace{2em}{\psi _{1,1}}(s,t)={\psi _{\eta }}(s){\psi _{\eta }}(t).\end{aligned}\]
The support of the spectrum of these real wavelets are oriented at $\pm {15^{\circ }}$, $\pm {45^{\circ }}$ and $\pm {75^{\circ }}$ (see Fig. 1) (Kingsbury, 2001; Selesnick et al., 2005).Fig. 1
The real oriented 2D dual (ℂWT) transform. (a) Explains each of the wavelet in the space domain. (b) Explains the support of the spectrum of the wavelets in the 2D frequency plane (Selesnick et al., 2005).
Let the two different sets of orthonormal filters $\{{\xi _{0}},{\xi _{1}}\}$ and $\{{\eta _{0}},{\eta _{1}}\}$, where ${\xi _{0}}$ and ${\eta _{0}}$ denote the low-pass filters, and ${\xi _{1}}$ and ${\eta _{1}}$ denote the high-pass filters. More precisely, $\{{\xi _{0}},{\xi _{1}}\}$ and $\{{\eta _{0}},{\eta _{1}}\}$ satisfy in relations:
\[\begin{aligned}{}& {\phi _{\xi }}(t)=\sum \limits_{n}{\xi _{0}}(n)\sqrt{2}{\phi _{\xi }}(2t-m),\hspace{1em}{\psi _{\xi }}(t)=\sum \limits_{n}{\xi _{1}}(n)\sqrt{2}{\phi _{\xi }}(2t-m),\\ {} & {\phi _{\eta }}(t)=\sum \limits_{n}{\eta _{0}}(n)\sqrt{2}{\phi _{\eta }}(2t-m),\hspace{1em}{\psi _{\eta }}(t)=\sum \limits_{n}{\eta _{1}}(n)\sqrt{2}{\phi _{\eta }}(2t-m),\end{aligned}\]
where $m\in \mathbb{Z}$. $\{{\xi _{0}},{\xi _{1}}\}$ and $\{{\eta _{0}},{\eta _{1}}\}$ make perfect reconstruction filter banks (Selesnick et al., 2005). The decomposition process performed by the filter bank has been shown in Fig. 2. The 2D separable wavelet transform, denoted by ${\mathcal{F}_{\eta \xi }}$, is defined to be applying filters ${\eta _{i}},i=0,1$ along the rows and ${\xi _{i}},i=0,1$ along the columns, and define ${\mathcal{F}_{\xi \eta }}$, ${\mathcal{F}_{\xi \xi }}$ and ${\mathcal{F}_{\eta \eta }}$ similarly.Therefore, the 2D dual-tree ℂWT transform is given by Selesnick et al. (2005):
The inverse of $\mathcal{F}$ is obtained as follows:
where ${I_{id}}$ is the identity matrix. It is easy check to ${\mathcal{F}^{-1}}$ is exactly the transpose of the $\mathcal{F}$. In this case, $\mathcal{F}$ is a tight frame.
(1)
\[ \mathcal{F}=\frac{1}{\sqrt{8}}\left[\begin{array}{c@{\hskip4.0pt}c@{\hskip4.0pt}c@{\hskip4.0pt}c}{I_{id}}& -{I_{id}}\\ {} {I_{id}}& {I_{id}}\\ {} & & {I_{id}}& {I_{id}}\\ {} & & {I_{id}}& -{I_{id}}\end{array}\right]\left[\begin{array}{c}{\mathcal{F}_{\xi \xi }}\\ {} {\mathcal{F}_{\eta \eta }}\\ {} {\mathcal{F}_{\eta \xi }}\\ {} {\mathcal{F}_{\xi \eta }}\end{array}\right].\](2)
\[ {\mathcal{F}^{-1}}=\frac{1}{\sqrt{8}}\left[{\mathcal{F}_{\xi \xi }^{-1}},{\mathcal{F}_{\eta \eta }^{-1}},{\mathcal{F}_{\eta \xi }^{-1}},{\mathcal{F}_{\xi \eta }^{-1}}\right]\left[\begin{array}{c@{\hskip4.0pt}c@{\hskip4.0pt}c@{\hskip4.0pt}c}{I_{id}}& {I_{id}}\\ {} -{I_{id}}& {I_{id}}\\ {} & & {I_{id}}& {I_{id}}\\ {} & & {I_{id}}& -{I_{id}}\end{array}\right],\]Fig. 2
Analysis filter bank for the 2D dual-tree discrete complex wavelet transform (Selesnick et al., 2005).
2.2 Expectation Maximization Segmentation (EMS) Method
Wells et al. (1996) propose an iterative segmentation method based on the expectation-maximization (EM) algorithm. This algorithm is used for data clustering into K regions by using the EM algorithm (Dempster et al., 1977) to determine parameters of a mixture of K Gaussian in the regions. The algorithm for vessel extraction is as follows:
Let $({x_{1}},{x_{2}},\dots ,{x_{N}})$ be MRA image, where ${x_{j}}$ denotes intensity value of j-th pixel such that ${x_{j}}\in \mathbb{R}$. Assume this image, which consists of two classes, background and vessels, is defined by Gaussian distributions with mean ${\mu _{i}}$ $(i=1,2)$ and deviation ${\sigma _{i}}$ $(i=1,2)$:
where the function $P(i|{x_{j}},{\mu _{i}^{\mathit{old}}},{\sigma _{i}^{\mathit{old}}})$ is the conditional probability of pixel j $(j=1,2,\dots ,N)$ that belongs to the class i $(i=1,2)$ at the current iteration and is defined as:
Also, in order to initialize two means ${\mu _{1}}$, ${\mu _{2}}$, two standard deviations ${\sigma _{1}}$, ${\sigma _{2}}$ and two weights ${\omega _{1}}$, ${\omega _{2}}$, in the iteration process, we let ${\mu _{1}}=\frac{m}{3}$, ${\mu _{2}}=\frac{2m}{3}$, ${\sigma _{1}}={\sigma _{2}}=m$ and ${\omega _{1}}={\omega _{2}}=0.5$, where $m=\max \{{x_{1}},{x_{2}},\dots ,{x_{N}}\}$.
\[ {G_{i}}(x|{\mu _{i}},{\sigma _{i}})=\frac{1}{\sqrt{2\pi {\sigma _{i}^{2}}}}\exp \bigg[-\frac{1}{2}{\bigg(\frac{x-{\mu _{i}}}{{\sigma _{i}}}\bigg)^{2}}\bigg].\]
The mixture model is
where ${\omega _{i}}$ $(i=1,2)$ is the weight of each classes i in the mixture model such that ${\omega _{1}}+{\omega _{2}}=1$. To estimate the values of ${\mu _{i}}$, ${\sigma _{i}}$ and ${\omega _{i}}$ $(i=1,2)$, we use EM algorithm as: (4)
\[\begin{aligned}{}& {\omega _{i}^{\mathit{new}}}=\frac{1}{N}{\sum \limits_{j=1}^{N}}P\big(i\big|{x_{j}},{\mu _{i}^{\mathit{old}}},{\sigma _{i}^{\mathit{old}}}\big),\end{aligned}\](5)
\[\begin{aligned}{}& {\mu _{i}^{\mathit{new}}}=\frac{{\textstyle\textstyle\sum _{j=1}^{N}}P(i|{x_{j}},{\mu _{i}^{\mathit{old}}},{\sigma _{i}^{\mathit{old}}}){x_{j}}}{{\textstyle\textstyle\sum _{j=1}^{N}}P(i|{x_{j}},{\mu _{i}^{\mathit{old}}},{\sigma _{i}^{\mathit{old}}})},\end{aligned}\](6)
\[\begin{aligned}{}& {\sigma _{i}^{\mathit{new}}}=\frac{{\textstyle\textstyle\sum _{j=1}^{N}}P(i|{x_{j}},{\mu _{i}^{\mathit{old}}},{\sigma _{i}^{\mathit{old}}}){({x_{j}}-{\mu _{i}^{\mathit{new}}})^{2}}}{{\textstyle\textstyle\sum _{j=1}^{N}}P(i|{x_{j}},{\mu _{i}^{\mathit{old}}},{\sigma _{i}^{\mathit{old}}})},\end{aligned}\](7)
\[ P\big(i\big|{x_{j}},{\mu _{i}^{\mathit{old}}},{\sigma _{i}^{\mathit{old}}}\big)=\frac{{\omega _{i}^{\mathit{old}}}{G_{i}}(x|{\mu _{i}^{\mathit{old}}},{\sigma _{i}^{\mathit{old}}})}{{\textstyle\textstyle\sum _{k=1}^{2}}{\omega _{k}^{\mathit{old}}}{G_{k}}(x|{\mu _{k}^{\mathit{old}}},{\sigma _{k}^{\mathit{old}}})}.\]EM algorithm is repeated until the log likelihood
where δ is an accuracy parameter.
It is shown in Dempster et al. (1977) that in many cases the EM algorithm enjoys pleasant convergence properties, namely that iterations will never worsen the value of the objective function.
Now the pixel j $(j=1,2,\dots ,N)$ belongs to the vessel class if
We construct the EMS mask M as follows:
where ${M_{j}}$ denotes the j-th value of mask M.
(9)
\[ {\omega _{2}}{G_{2}}({x_{j}}|{\mu _{2}},{\sigma _{2}})\geqslant {\omega _{1}}{G_{1}}({x_{j}}|{\mu _{1}},{\sigma _{1}}).\](10)
\[ {M_{j}}=\left\{\begin{array}{l@{\hskip4.0pt}l}1,\hspace{1em}& \text{if the pixel}\hspace{2.5pt}j\hspace{2.5pt}\text{is a vessel pixel},\\ {} 0,\hspace{1em}& \text{otherwise},\end{array}\right.\]2.3 Tight Frame Method for Vessel Extraction Based on 2-D Dual-Tree ℂWT
The tight- rame algorithm (TFA) used in Cai et al. (2008), Chan et al. (2003). This algorithm can be presented in the following form:
where ${f^{(i)}}$ is an approximate solution at the i-th iteration, $\mathcal{R}$ is a problem-dependent operator (Cai et al., 2013). ${\mathcal{T}_{\lambda }}(.)$ is the soft thresholding operator given by
where λ is thresholding parameter, ${t_{\lambda }}(x)=\operatorname{sgn}(x)\max \{0,|x|-\lambda \}$ and ${\mathbf{x}=[{x_{1}},{x_{2}},\dots ,{x_{n}}]^{T}}$. The operator $\mathcal{F}$ is the tight frame operator defined by (1). Suppose that $p=0.08839$, $q=0.69588$ and $r=0.01127$. Define
and
It can be shown that the matrices ${F_{\xi \xi }}$, ${F_{\eta \eta }}$, ${F_{\xi \eta }}$ and ${F_{\eta \xi }}$ in Eq. (1) can be found by computing ${\mathcal{F}_{{\xi _{i}}}}$ and ${\mathcal{F}_{{\eta _{i}}}},(i=0,1)$ as
(13)
\[ {\mathcal{T}_{\lambda }}(\mathbf{x})={\big[{t_{\lambda }}({x_{1}}),{t_{\lambda }}({x_{2}}),\dots ,{t_{\lambda }}({x_{n}})\big]^{T}},\]
\[\begin{aligned}{}& {\mathcal{F}_{{\xi _{i}}}}=\left[\begin{array}{l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l}{\xi _{i}}(1)\hspace{1em}& {\xi _{i}}(2)\hspace{1em}& {\xi _{i}}(3)\hspace{1em}& \dots \hspace{1em}& \hspace{1em}& {\xi _{i}}(8)\hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}\\ {} \hspace{1em}& \hspace{1em}& {\xi _{i}}(1)\hspace{1em}& {\xi _{i}}(2)\hspace{1em}& {\xi _{i}}(3)\hspace{1em}& \dots \hspace{1em}& \hspace{1em}& {\xi _{i}}(8)\hspace{1em}& \hspace{1em}& \hspace{1em}\\ {} \hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}& \ddots \hspace{1em}& \ddots \hspace{1em}& \ddots \hspace{1em}& \ddots \hspace{1em}& \hspace{1em}& \ddots \hspace{1em}\\ {} \hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}& {\xi _{i}}(1)\hspace{1em}& {\xi _{i}}(2)\hspace{1em}& {\xi _{i}}(3)\hspace{1em}& \dots \hspace{1em}& \hspace{1em}& {\xi _{i}}(8)\end{array}\right],\\ {} & {\mathcal{F}_{{\eta _{i}}}}=\left[\begin{array}{l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l@{\hskip4.0pt}l}{\eta _{i}}(1)\hspace{1em}& {\eta _{i}}(2)\hspace{1em}& {\eta _{i}}(3)\hspace{1em}& \dots \hspace{1em}& \hspace{1em}& {\eta _{i}}(8)\hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}\\ {} \hspace{1em}& \hspace{1em}& {\eta _{i}}(1)\hspace{1em}& {\eta _{i}}(2)\hspace{1em}& {\eta _{i}}(3)\hspace{1em}& \dots \hspace{1em}& \hspace{1em}& {\eta _{i}}(8)\hspace{1em}& \hspace{1em}& \hspace{1em}\\ {} \hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}& \ddots \hspace{1em}& \ddots \hspace{1em}& \ddots \hspace{1em}& \ddots \hspace{1em}& \hspace{1em}& \ddots \hspace{1em}\\ {} \hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}& \hspace{1em}& {\eta _{i}}(1)\hspace{1em}& {\eta _{i}}(2)\hspace{1em}& {\eta _{i}}(3)\hspace{1em}& \dots \hspace{1em}& \hspace{1em}& {\eta _{i}}(8)\end{array}\right].\end{aligned}\]
It can be shown that the frame obtained from Eqs. (14)–(17) has a tight frame property, so ${F^{-1}}={F^{T}}$. It also performs denoising on the image. There is no loss of generality in assuming that the intensity of the given image f is in $[0,1]$.The initial approximation image ${f^{(0)}}=f$ and boundary pixels ${\Gamma ^{(0)}}$ are considered for starting algorithm, where f is the given image and ${\Gamma ^{(0)}}$ is as follows:
where Ω is the index set of all the pixels in the image and ${[\nabla f]_{j}}$ is the discrete gradient of f at the j-th pixel. Below, the i-th iteration is described in details.
(18)
\[ {\Gamma ^{(0)}}=\big\{j\in \Omega \big|\hspace{0.1667em}\big\| {[\nabla f]_{j}}\big\| \geqslant \varepsilon \big\},\]First, interval $[{a_{i}},{b_{i}}]$ is estimated as follows. Define the mean pixel value
where ${f_{j}^{(i)}}$ is the j-pixel value in the approximation image ${f^{(i)}}$. Let ${\Gamma _{-}^{(i)}}=\{j\in {\Gamma ^{(i)}}:{f_{j}^{(i)}}\leqslant {\rho ^{(i)}}\}$ and ${\Gamma _{+}^{(i)}}=\{j\in {\Gamma ^{(i)}}:{f_{j}^{(i)}}\geqslant {\rho ^{(i)}}\}$, then the mean pixel values of the two sets separated by ${\rho ^{(i)}}$, is computed:
Now threshold image ${f_{j}^{(i+1/2)}}$ is obtained as follows: The image ${f^{(i)}}$ is separated into three parts by applying the interval $[{a_{i}},{b_{i}}]\subseteq [0,1]$. Those ${f_{j}^{(i)}}$’s that are smaller than ${a_{i}}$, threshold to 0, those larger than ${b_{i}}$ to 1, and those between, are stretched between 0 and 1 by applying a simple linear contrast stretch as follows, and if there isn’t any ${f_{j}^{(i)}}$ between ${a_{i}}$ and ${b_{i}}$, then ${f_{j}^{(i+1/2)}}$ is a binary image and the algorithm stops. More precisely, define
If ${f_{j}^{(i+1/2)}}=0$, then pixel j is in the background, or if ${f_{j}^{(i+1/2)}}=1$, then pixel j is inside the vessel. The remaining pixels are putted into a set of
(19)
\[ {\rho ^{(i)}}=\frac{1}{|{\Gamma ^{(i)}}|}\sum \limits_{j\in {\Gamma ^{(i)}}}{f_{j}^{(i)}},\]
\[ {\rho _{-}^{(i)}}=\frac{1}{|{\Gamma _{-}^{(i)}}|}\sum \limits_{j\in {\Gamma _{-}^{(i)}}}{f_{j}^{(i)}},\hspace{2em}{\rho _{+}^{(i)}}=\frac{1}{|{\Gamma _{+}^{(i)}}|}\sum \limits_{j\in {\Gamma _{+}^{(i)}}}{f_{j}^{(i)}}.\]
It follows that ${a_{i}}$ and ${b_{i}}$ are defined by
(20)
\[ {a_{i}}=\max \bigg\{\frac{{\rho ^{(i)}}+{\rho _{-}^{(i)}}}{2},0\bigg\}\hspace{0.1667em},\hspace{2em}{b_{i}}=\min \bigg\{\frac{{\rho ^{(i)}}+{\rho _{+}^{(i)}}}{2},1\bigg\}.\]
\[\begin{aligned}{}& {M_{i}}=\max \big\{{f_{j}^{(i)}}\big|\hspace{0.1667em}{a_{i}}\leqslant {f_{j}^{(i)}}\leqslant {b_{i}},\hspace{0.1667em}j\in {\Gamma ^{(i)}}\big\},\\ {} & {m_{i}}=\min \big\{{f_{j}^{(i)}}\big|\hspace{0.1667em}{a_{i}}\leqslant {f_{j}^{(i)}}\leqslant {b_{i}},\hspace{0.1667em}j\in {\Gamma ^{(i)}}\big\},\end{aligned}\]
then threshold image is
(21)
\[ {f_{j}^{(i+1/2)}}=\left\{\begin{array}{l@{\hskip4.0pt}l}0,\hspace{1em}& \text{if}\hspace{2.5pt}{f_{j}^{(i)}}\leqslant {a_{i}},\\ {} \frac{{f_{j}^{(i)}}-{m_{i}}}{{M_{i}}-{m_{i}}},\hspace{1em}& \text{if}\hspace{2.5pt}{a_{i}}\leqslant {f_{j}^{(i)}}\leqslant {b_{i}},\\ {} 1,\hspace{1em}& \text{if}\hspace{2.5pt}{f_{j}^{(i)}}\geqslant {b_{i}}.\end{array}\right.\]Next, ${f_{j}^{(i+1/2)}}$ is denoised on ${\Gamma ^{(i+1)}}$ by using the tight frame transformation (12) on ${\Gamma ^{(i+1)}}$ and ${f_{j}^{(i+1)}}$ is obtained. More precisely,
TFA algorithm is repeated until ${\Gamma ^{(i+1)}}=\varnothing $, or equivalently, all the pixels of ${f^{(i+1/2)}}$ are either of value 0 or 1. More details about the TFA algorithm for segmentation can be reviewed in Cai et al. (2013).
The following theorem shows the convergence of TFA algorithm.
2.4 Hybrid of EMS and TFA Methods
In our studies, as we observed in numerical examples, despite of all the accuracy in the vessel extraction by EMS algorithm, it has artifacts in its segmentation and doesn’t show the boundaries of the vessels smoothly, while the TFA algorithm exhibits the boundaries smoothly and reduces noise. On the other hand, TFA algorithm is wrongly segmented in the regions of image with intensity close to the background, while EMS algorithm has worked well. However, in those situations where EMS algorithm has unsatisfactory results, TFA algorithm has worked well and vice versa. Thus, the main idea of our method is to use the combination of EMS and TFA algorithms for vessel extraction.
Let f be the MRA image and
where M is the EMS mask which is defined in Eq. (10) and “∘” is the Hadamard product operator. Suppose
where parameters of α and β are constant scalars. By choosing the appropriate parameters α and β we can get better segmentation results than the TFA and EMS algorithms.
Now, we apply TFA algorithm on image g. Assume, we have new image ${g^{(i)}}$ and a set of all possible boundary pixels ${\Gamma ^{(i)}}$ (22) in the beginning of the i-th iteration. Then we (a) compute an interval of possible boundary pixel values $[{a_{i}},{b_{i}}]$ (20); (b) use the interval $[{a_{i}},{b_{i}}]$ to separate the image into background pixels, possible boundary pixels and pixels in the vessels (21); (c) we use the tight frame transformation to denoise on the possible boundary pixels to get a approximation image ${g^{(i+1)}}$ (23). We stop when the image becomes binary. We have summarized the hybrid method in Algorithm 1.
If the operator $\mathcal{P}$ maps a given image to image segmented by TFA algorithm then the presented method provides a collection of segmented images for different values of α and β for the given image f as:
\[ \mathcal{S}=\big\{\tilde{f}\big|\hspace{0.1667em}\tilde{f}=\mathcal{P}(\alpha {f_{\textit{EM}}}+\beta f),\hspace{0.1667em}0\leqslant \alpha ,\beta \leqslant 1\big\}.\]
Note that the image segmented by TFA algorithm can be obtained from the presented method by choosing $\alpha =0$ and $\beta =1$ in (25).| 1. Input: given MRA image f. | |
| While happen Eq. (8), | |
| Step 1 | Update ${\omega _{i}^{\mathit{new}}}$, ${\mu _{i}^{\mathit{new}}}$ and ${\sigma _{i}^{\mathit{new}}}$, for $(i=1,2)$ |
| by Eqs. (4)–(6), respectively. | |
| 2. Compute EMS mask by Eq. (10). | |
| 3. Do EMS mask on f and set in ${f_{EM}}$. | |
| Step 2 | 4. Constitute $g=\alpha {f_{EM}}+\beta f$ |
| 5. Set ${g^{(0)}}=g$ and ${\Gamma ^{(0)}}$ by Eq. (18) | |
| 6. While ${g^{(i+1)}}$ becomes a binary image | |
| (a). Obtain $[{a_{i}},{b_{i}}]$ by Eq. (20). | |
| (b). Obtain ${g^{(i+1/2)}}$ by Eq. (21). | |
| Step 3 | (c). Stop if ${g^{(i+1/2)}}$ is a binary image. |
| (d). Obtain ${\Gamma ^{(i+1)}}$ by Eq. (22). | |
| (e). Update ${g^{(i+1/2)}}$ to ${g^{(i+1)}}$ by Eq. (23). | |
| 7. Output: binary image ${g^{(i+1)}}$. |
Algorithm 1
Hybrid of EMS and TFA algorithm
Proof.
Suppose f is a given image. We obtain ${f_{EM}}=M\circ f$, where M is a binary image that is obtained by Eq. (10). We constitute a new image $g=\alpha {f_{\mathit{EM}}}+\beta f$ with linear combination of ${f_{\mathit{EM}}}$ and f, where $0\leqslant \alpha $, $\beta \leqslant 1$. According to theorem 2, iteratively TFA algorithm converges to a binary image $\tilde{g}=\mathcal{P}(g)$, where $\tilde{g}$ is the result of the presented method. □
Finally, let us estimate the computation cost of our method for a given image with n pixels. The complexity of TFA algorithm is $O(n)$ per iteration (Cai et al., 2013). It bears mentioning that in each iteration we need to evaluate K Gaussian densities for N points in the E-step, and that this scales as $O(KN{d^{3}})$. Moreover, the M-step requires $O(KN{d^{2}})$ work in order to update the estimates of the Gaussian parameters, which d is dimension of data (Fajardo and Liang, 2017). In this paper, dimension of data is $d=1$, $K=2$. Thus the complexity of EM algorithm for our data is $O(n)$ per iteration. Therefore, the complexity of the presented method is $O(n)$ per iteration.
3 Results
In this section, we test the presented method on five different images that include the simulated, carotid, kidney, abdominal and circle of Willis inverted MIP of vascular systems. The thresholding parameters λ and the accuracy parameter δ used in (12) and (8) are chosen to be $\lambda =0.1$ and $\delta =0.0001$. ε used in (18) is chosen to be $\varepsilon =0.003$ except in Example 1 in which we set $\varepsilon =0.02$. Weights α and β are chosen manually. We show the results for the tight frame of 2D dual-tree complex wavelet in TFA algorithm (CTFA) and the presented method (EMCTFA). The number of wavelet levels is 4. Also, we give the results of our method by the cubic B-spline wavelets as a tight frame (EMSTFA). We compare the results of the presented method with EMS algorithm (Wilson and Noble, 1997, 1999), Chan-Vese active contour model (Chan and Vese, 2001) (Chan-Vese), B-spline wavelet tight frame algorithm (STFA) and the dual-tree complex wavelet tight frame algorithm (CTFA).
Table 1
The number of pixels in the set ${\Gamma ^{(i)}}$ at each iteration.
| Examples | Method | $i=1$ | $i=2$ | $i=3$ | $i=4$ | $i=5$ | $i=6$ | $i=7$ | $i=8$ |
| Annuluses | TFAS | 29529 | 4411 | 1010 | 238 | 54 | 11 | 0 | – |
| CTFA | 29415 | 3121 | 694 | 184 | 50 | 8 | 0 | – | |
| EMSTFA | 22468 | 4089 | 974 | 244 | 65 | 10 | 0 | – | |
| EMCTFA | 22443 | 3007 | 664 | 167 | 45 | 11 | 0 | – | |
| Kidney | STFA | 9985 | 2092 | 472 | 110 | 31 | 7 | 0 | – |
| CTFA | 9985 | 2051 | 473 | 134 | 36 | 10 | 0 | – | |
| EMSTFA | 7510 | 1837 | 512 | 130 | 44 | 9 | 0 | – | |
| EMCTFA | 4704 | 1201 | 305 | 87 | 28 | 5 | 0 | – | |
| Carotid | STFA | 1799 | 341 | 84 | 0 | – | – | – | – |
| CTFA | 1799 | 372 | 94 | 25 | 6 | 0 | – | – | |
| EMSTFA | 1419 | 219 | 69 | 0 | – | – | – | – | |
| EMCTFA | 1419 | 294 | 73 | 11 | 0 | – | – | – | |
| Abdonomial | STFA | 30871 | 7031 | 1709 | 576 | 169 | 37 | 11 | 0 |
| CTFA | 30871 | 7517 | 1964 | 581 | 188 | 62 | 17 | 0 | |
| EMSTFA | 3194 | 738 | 257 | 67 | 18 | 2 | 0 | – | |
| EMCTFA | 4146 | 1128 | 318 | 87 | 24 | 5 | 0 | – | |
| Wills | STFA | 33595 | 8670 | 2176 | 553 | 137 | 34 | 4 | 0 |
| CTFA | 33595 | 8738 | 2393 | 697 | 221 | 70 | 20 | 3 | |
| EMSTFA | 2176 | 518 | 163 | 39 | 6 | 0 | – | – | |
| EMCTFA | 2176 | 636 | 171 | 56 | 12 | 0 | – | – |
The cardinality of ${\Gamma ^{(i)}}$ at each iteration and the number of iterations given in Table 1, show the convergence speed of STFA, CTFA, EMSTFA, and EMCTFA methods for the presented examples. In fact, the cardinality of ${\Gamma ^{(i)}}$ at each iteration shows the number of pixels that are unclassified yet in the i-th iteration.
4 Discussion
Fig. 3
Example 1. (a) Simulated image; (b) Ground truth image; (c) Noisy image; (d)–(i) Results of methods EMS, Chan-Vese, STFA, CTFA, EMSTFA and EMCTFA, respectively; (j)–(n) Differences between the results of mentioned methods (d)–(i) and the ground truth image (b). Numbers of wrongly-detected pixels are shown in braces.
The results of six examples show that EMCTFA method is better than other mentioned methods. In this section, we investigate the results of six examples.
Example 1.
This example is a $512\times 512$ simulated image of vascular system (Fig. 3(a)). The ground truth (Fig. 3(b)), is constituted of four straight stripes and four annuluses with varying widths.
Similar to Blood flowing in the vessels, the centre line of our simulated image has higher intensity and is less towards the boundary. The image shown in Fig. 3(a) is the result of a simulation of the vascular system. We obtain the noisy image in Fig. 3(c) by applying Gaussian noise with mean 0.01 and variance 0.001 to the simulated image. Figures 3(d)–(i) show the resulting images for the EMS, Chan-Vese, STFA, CTFA, EMSTFA and EMCTFA (the presented method) methods. Also Figs. 3(j)–(n), show their differences by the ground truth image (Fig. 3(b)). In the presented method, we used $\varepsilon =0.02$ in (18). We see that in these figures, the EMCTFA method extracts the vessels better than the other mentioned methods.
Example 2.
This example is a $128\times 128$ MRA image (Fig. 5(a)). Figures 5(b)–(g) show the results of vessel extraction by EMS, STFA, Chan-Vese, CTFA, EMSTFA and EMCTFA methods (the presented method). We used $\alpha =0.25$ and $\beta =1$ in the presented method to segment this MRA image.
For this image, EMS and STFA algorithms do not obtain good results, since they have created many artifacts near the boundary of the vessels, while EMCTFA method has acceptable results and performed well. The results by Chan-Vese are not satisfactory since Chan-Vese cannot detect a large part of the vessels. To compare EMS and CTFA and EMCTFA methods more closely, we enlarge the top part of Figs. 5(c), (d) and (f) and depict them in Figs. 5(g)–(i) that show the presented method removes some artifacts close to the boundary.
For this MRA image, Table 1 shows that CTFA and EMCTFA methods converged with the same speed (within 5 iterations) but the cardinality of ${\Gamma ^{(i)}}$ at each iteration of the presented method is less than that of the TFA algorithm.
Fig. 5
Example 2. Carotid vascular system extraction. (a) Given image; (b)–(g) Results by EMS, STFA, Chan-Vese, CTFA and EMSTFA algorithms, respectively. (f) Results of the presented method; (h)–(j) are the zoomed-in top parts of (c), (d) and (f), respectively. Green lines denote boundary of vessel extraction.
Example 3.
This example is a $256\times 256$ MRA image (Fig. 6(a)). Figures 6(b)–(g) show the results of vessel extraction by EMS, STFA, Chan-Vese, CTFA and the presented methods, respectively. In the presented method, we used $\alpha =0.2$ and $\beta =1$ to segment this MRA image.
Figure 6 shows that EMS, STFA and EMSTFA methods give unsatisfactory results since they have obviously wrongly detected pixels in their result of segmentation and do not detect smoothness of boundaries well. Also, Chan-Vese is unable to recover the small occlusions along the coherence direction. EMCTFA method has the same result as CTFA algorithm and detected smoothness of boundaries and also reconstructed structures which present small occlusions along the coherence direction. As Table 1 shows, the presented method converges faster than other mentioned methods.
Fig. 6
Example 3. Kidney vascular system extraction. (a) Given image; (b)–(f) Results by EMS, STFA, Chan-Vese, CTFA, EMSTFA methods, respectively; (g) Results of the presented method. Green lines denote boundary of vessel extraction.
Example 4.
This example is one slide of 3D CE-MRA image of the abdominal vascular system whose size is $512\times 512$ (Fig. 7(a)) from http://www.mr-tip.com. Figures 7(b)–(g) show the results of vessel extraction by EMS, STFA, Chan-Vese, CTFA and the presented methods, respectively. In our method, we used $\alpha =1$ and $\beta =0.15$ to segment this MRA image.
Obviously, the segmented vessels by Chan-Vese give unsatisfactory results. Although STFA and CTFA algorithms detect some thinner vessels, it is unable to separate regions of the background whose intensity is close to the vessel intensity. For this reason those methods have unsatisfactory results in this case. To compare EMS algorithm and the presented method closely, we enlarge the rectangular boxes in Figs. 7(e)–(g). They exhibit that EMS algorithm has many artifacts which are well reduced by our method. Table 1 shows that EMCTFA method converges in 6 iterations. Also, it shows that the presented method converges faster than CTFA method. In the presented method, at the first iteration, the number of pixels that are not classified is 4146 pixels, while in the CTFA method it is 30871 pixels.
Fig. 7
Example 4. CE-MRA image of abdominal vascular system extraction. (a) Given image; (b)–(e) Results by EMS, STFA, Chan-Vese, CTFA and EMSTFA methods, respectively; (g) Results of our EMCTFA method; (h)–(j) are the zoomed-in red rectangular parts of (a), (b) and (g).
Example 5.
This example is a $512\times 512$ TOF-MRA Circle of willis inverted MIP of carotid vascular system (Fig. 8(a)) from http://www.mr-tip.com. Figures 8(b)–(g) show the results of vessel extraction by EMS, STFA, Chan-Vese, CTFA and the presented methods, respectively. In the presented method we used $\alpha =1$ and $\beta =0.25$ to segment this MRA image.
The extracted vessels by Chan-Vese, STFA and CTFA methods give unsatisfactory results. These methods are unable to recover the vessels of some regions, for example, see the upper right side of Figs. 8(a), (b), (d) and (f) (for better viewing, we enlarged this part in Figs. 8(e)–(h)). EMS method has artifacts near vessel boundaries and other regions, and the presented method removes most of them. This example also shows the ability of the presented method to segment regions of the background so that their intensity is close to the vessel and does not detect them as a vessel inaccurately. As shown in Table 1, in this example, EMCTFA method converges by less iteration than CTFA algorithm too.
Fig. 8
Example 5. TOF-MRA Circle of Willis Inverted MIP of carotid vascular system extraction. (a) Given image; (b)–(f) Results by EMS, Chan-Vese, STFA, CTFA and EMSTFA methods, respectively; (g) Results of the presented method; (h)–(k) are the zoomed-in parts of (a), (b), (d) and (e).
Example 6.
In this example, we applied our method on kidney MRA image for some values of α and β. Figure 9(a) shows the original image. Figures 9(b)–(g) show results of applying the presented method with values of ($\alpha =0$, $\beta =1$), ($\alpha =0.5$, $\beta =1$), ($\alpha =0.75$, $\beta =1$), ($\alpha =1$, $\beta =1$), ($\alpha =1$, $\beta =.25$) and ($\alpha =1$, $\beta =0$), respectively (red lines).
Fig. 9
Example 6. Results of vessel segmentation for image (a), obtained with the presented method for (b) $\alpha =0$, $\beta =1$; (c) $\alpha =0.5$, $\beta =1$; (d) $\alpha =0.75$, $\beta =1$; (e) $\alpha =1$, $\beta =1$; (f) $\alpha =1$, $\beta =0.25$; (g) $\alpha =1$, $\beta =0$. Blue lines denote boundary of vessel extraction.
5 Conclusions
This paper presents an efficient vessel extraction method based on tight frame and EMS algorithms. The presented method produces sets of image segmentation for a given image. By choosing the appropriate parameters, we can obtain a well-defined segmental image that has the advantages of both the TFA and EMS algorithms relatively. The test of algorithm on real MRA images demonstrates the ability of the presented method in vessel extraction.
Our method has an advantage of fast implementation, gives very accurate vessel extraction, extracts vessel boundary smoothly and avoids artifacts. We have proved convergence in the presented method to a binary image.
