INFORMATICA Informatica 1822-8844 0868-4952 0868-4952 Vilnius University INFOR435 10.15388/20-INFOR435 Research Article Hybrid Vessel Extraction Method Based on Tight-Frame and EM Algorithms by Using 2D Dual Tree Complex Wavelet AbdollahiFaridf.abdollahi@tabrizu.ac.ir1

F. Abdollahi was born in 1988. He is a PhD student of Applied Mathematics at Faculty of Mathematical Sciences, University of Tabriz, Iran. His research interests include image processing and wavelets.

 LakestaniMehrdadlakestani@tabrizu.ac.ir1

M. Lakestani received his BSc in applied mathematics from Department of Mathematics at the University of Tabriz, in 1998, MSc and PhD in applied mathematics in 2000 and 2005, respectively, at Amirkabir University of Technology. He has been a professor of mathematics at the University of Tabriz since 2015. His research interests include wavelets, numerical methods, and image processing.

 RazzaghiMohsenrazzaghi@math.msstate.edu2∗

M. Razzaghi received his BSc degree in mathematics and a PhD degree in applied mathematics both from the University of Sussex in England. Since 1986, he has been with the Department of Mathematics and Statistics at Mississippi State University, where he is currently a professor and the department head. He was the recipient of two Fulbright scholar programmes, one in 2011–2012 and another in 2015–2016, and one Fulbright specialist programme in 2019, in Romania. His current area of research centres on orthogonal functions, optimal control, wavelets, fractional calculus, and their applications in mathematical modelling, and engineering. He has over 190 refereed journal publications in mathematics, mathematical physics, and engineering. One of his papers, coauthored with one of his PhD students, was cited over 700 times.

 Department of Applied Mathematics, Faculty of Mathematical Sciences, University of Tabriz, Tabriz, Iran Department of Mathematics and Statistics, Mississippi State University, MS 39762, USA Corresponding author. 2021 1122021 321122 92019 112020 © 2021 Vilnius University2021 Open access article under the CC BY license.

The vessel extraction is very important for the vascular disease diagnosis and grading of the stenoses and aneurysms in vessels. This aids in brain surgery and making angioplasty. The presence of noise in the MRA image, etc., turns the vessel extraction into a difficult problem. In this paper, we derive a vessel extraction algorithm based on TFA and EMS algorithms. We prove the convergence of the proposed method within a few iterations. Results of applying the presented method on real 2D MRA images demonstrate that our method is very efficient.

 VESSEL extraction MRA images TFA algorithm EMS algorithm 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.

 Methods 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 (ϕi)i∈I in H (Hilbert space) is a tight frame for H if F has the perfect reconstruction property: F∗F=IH , where operators F and F∗ the adjoint of F defined as Kutyniok and Labate (2012). F:H→ℓ2(I),x⟼(⟨x,ϕi⟩)i∈I. So for each vector x∈H x=∑i∈I(⟨x,ϕi⟩)ϕi. There are various tight frame systems. Some of these tight frames are wavelets for Hilbert space L2(R) : {ϕm=ϕ(.−m):m∈Z}∪{ψj,m=ψ(2j.−m):j,m∈Z}, where ϕ∈L2(R) and ψ∈L2(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 L2(R2) 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 ϕ(t)=ϕξ(t)+jϕη(t) and ψ(t)=ψξ(t)+jψη(t) . By taking the real part of each of these six complex wavelets ψ(s)ψ(t) , ψ(s)ψ(t)‾ , ϕ(s)ψ(t) , ψ(s)ϕ(t) , ϕ(s)ψ(t)‾ and ψ(s)ϕ(t)‾ , we obtain real oriented 2D wavelets. Specifically, we obtain the following six wavelets: ψi(s,t)=12(ψ1,i(s,t)−ψ2,i(s,t)),ψi+3(s,t)=12(ψ1,i(s,t)+ψ2,i(s,t)), for i=1,2,3 , where the two separable 2D wavelet bases are defined by: ψ1,1(s,t)=ϕξ(s)ψξ(t),ψ1,1(s,t)=ϕη(s)ψη(t),ψ1,1(s,t)=ψξ(s)ϕξ(t),ψ1,1(s,t)=ψη(s)ϕη(t),ψ1,1(s,t)=ψξ(s)ψξ(t),ψ1,1(s,t)=ψη(s)ψη(t). The support of the spectrum of these real wavelets are oriented at ±15∘ , ±45∘ and ±75∘ (see Fig. 1) (Kingsbury, 2001; Selesnick et al., 2005).

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 {ξ0,ξ1} and {η0,η1} , where ξ0 and η0 denote the low-pass filters, and ξ1 and η1 denote the high-pass filters. More precisely, {ξ0,ξ1} and {η0,η1} satisfy in relations: ϕξ(t)=∑nξ0(n)2ϕξ(2t−m),ψξ(t)=∑nξ1(n)2ϕξ(2t−m),ϕη(t)=∑nη0(n)2ϕη(2t−m),ψη(t)=∑nη1(n)2ϕη(2t−m), where m∈Z . {ξ0,ξ1} and {η0,η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 Fηξ , is defined to be applying filters ηi,i=0,1 along the rows and ξi,i=0,1 along the columns, and define Fξη , Fξξ and Fηη similarly.

Therefore, the 2D dual-tree ℂWT transform is given by Selesnick et al. (2005): F=18Iid−IidIidIidIidIidIid−IidFξξFηηFηξFξη. The inverse of F is obtained as follows: F−1=18Fξξ−1,Fηη−1,Fηξ−1,Fξη−1IidIid−IidIidIidIidIid−Iid, where Iid is the identity matrix. It is easy check to F−1 is exactly the transpose of the F . In this case, F is a tight frame.

Analysis filter bank for the 2D dual-tree discrete complex wavelet transform (Selesnick et al., 2005).

 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 (x1,x2,…,xN) be MRA image, where xj denotes intensity value of j-th pixel such that xj∈R . Assume this image, which consists of two classes, background and vessels, is defined by Gaussian distributions with mean μi (i=1,2) and deviation σi (i=1,2) : Gi(x|μi,σi)=12πσi2exp[−12(x−μiσi)2]. The mixture model is P(x)= ∑i=12ωiGi(x|μi,σi), where ωi (i=1,2) is the weight of each classes i in the mixture model such that ω1+ω2=1 . To estimate the values of μi , σi and ωi (i=1,2) , we use EM algorithm as: ωinew=1N ∑j=1NP(i|xj,μiold,σiold), μinew=∑j=1NP(i|xj,μiold,σiold)xj∑j=1NP(i|xj,μiold,σiold), σinew=∑j=1NP(i|xj,μiold,σiold)(xj−μinew)2∑j=1NP(i|xj,μiold,σiold), where the function P(i|xj,μiold,σiold) is the conditional probability of pixel j (j=1,2,…,N) that belongs to the class i (i=1,2) at the current iteration and is defined as: P(i|xj,μiold,σiold)=ωioldGi(x|μiold,σiold)∑k=12ωkoldGk(x|μkold,σkold). Also, in order to initialize two means μ1 , μ2 , two standard deviations σ1 , σ2 and two weights ω1 , ω2 , in the iteration process, we let μ1=m3 , μ2=2m3 , σ1=σ2=m and ω1=ω2=0.5 , where m=max{x1,x2,…,xN} .

EM algorithm is repeated until the log likelihood log ∏i=12Gi(x|μinew,σinew)−log ∏i=12Gi(x|μiold,σiold)<δ, 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,…,N) belongs to the vessel class if ω2G2(xj|μ2,σ2)⩾ω1G1(xj|μ1,σ1). We construct the EMS mask M as follows: Mj=1,if the pixeljis a vessel pixel,0,otherwise, where Mj denotes the j-th value of mask M.

 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: f(i+1/2)=R(fi), f(i+1)=FTTλ(Ff(i+1/2)),i=1,2,…, where f(i) is an approximate solution at the i-th iteration, R is a problem-dependent operator (Cai et al., 2013). Tλ(.) is the soft thresholding operator given by Tλ(x)=[tλ(x1),tλ(x2),…,tλ(xn)]T, where λ is thresholding parameter, tλ(x)=sgn(x)max{0,|x|−λ} and x=[x1,x2,…,xn]T . The operator F is the tight frame operator defined by (1). Suppose that p=0.08839 , q=0.69588 and r=0.01127 . Define ξ0=[−p,p,q,q,p,−p,r,r]T, ξ1=[−r,r,p,p,−q,q,−p,−p]T, η0=[r,r,−p,p,q,q,p,−p]T, and η1=[−p,−p,q,−q,p,p,r,−r]T. It can be shown that the matrices Fξξ , Fηη , Fξη and Fηξ in Eq. (1) can be found by computing Fξi and Fηi,(i=0,1) as Fξi=ξi(1)ξi(2)ξi(3)…ξi(8)ξi(1)ξi(2)ξi(3)…ξi(8)⋱⋱⋱⋱⋱ξi(1)ξi(2)ξi(3)…ξi(8),Fηi=ηi(1)ηi(2)ηi(3)…ηi(8)ηi(1)ηi(2)ηi(3)…ηi(8)⋱⋱⋱⋱⋱ηi(1)ηi(2)ηi(3)…ηi(8). It can be shown that the frame obtained from Eqs. (14)–(17) has a tight frame property, so F−1=FT . 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 Γ(0) are considered for starting algorithm, where f is the given image and Γ(0) is as follows: Γ(0)={j∈Ω|‖[∇f]j‖⩾ε}, where Ω is the index set of all the pixels in the image and [∇f]j is the discrete gradient of f at the j-th pixel. Below, the i-th iteration is described in details.

First, interval [ai,bi] is estimated as follows. Define the mean pixel value ρ(i)=1|Γ(i)|∑j∈Γ(i)fj(i), where fj(i) is the j-pixel value in the approximation image f(i) . Let Γ−(i)={j∈Γ(i):fj(i)⩽ρ(i)} and Γ+(i)={j∈Γ(i):fj(i)⩾ρ(i)} , then the mean pixel values of the two sets separated by ρ(i) , is computed: ρ−(i)=1|Γ−(i)|∑j∈Γ−(i)fj(i),ρ+(i)=1|Γ+(i)|∑j∈Γ+(i)fj(i). It follows that ai and bi are defined by ai=max{ρ(i)+ρ−(i)2,0},bi=min{ρ(i)+ρ+(i)2,1}. Now threshold image fj(i+1/2) is obtained as follows: The image f(i) is separated into three parts by applying the interval [ai,bi]⊆[0,1] . Those fj(i) ’s that are smaller than ai , threshold to 0, those larger than bi 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 fj(i) between ai and bi , then fj(i+1/2) is a binary image and the algorithm stops. More precisely, define Mi=max{fj(i)|ai⩽fj(i)⩽bi,j∈Γ(i)},mi=min{fj(i)|ai⩽fj(i)⩽bi,j∈Γ(i)}, then threshold image is fj(i+1/2)=0,iffj(i)⩽ai,fj(i)−miMi−mi,ifai⩽fj(i)⩽bi,1,iffj(i)⩾bi. If fj(i+1/2)=0 , then pixel j is in the background, or if fj(i+1/2)=1 , then pixel j is inside the vessel. The remaining pixels are putted into a set of Γ(i+1)={j|0<fj(i+1/2)<1,j∈Ω}.

Next, fj(i+1/2) is denoised on Γ(i+1) by using the tight frame transformation (12) on Γ(i+1) and fj(i+1) is obtained. More precisely, fj(i+1)=fji+1/2,ifj∉Γ(i+1),[FTTλ(Ffi+1/2)]j,otherwise.

TFA algorithm is repeated until Γ(i+1)=∅ , 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. (See Cai <italic>et al.</italic>, <xref ref-type="bibr" rid="j_infor435_ref_005">2013</xref>).

TFA converges to a binary image within a finite number of steps.

 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 fEM=M∘f, where M is the EMS mask which is defined in Eq. (10) and “∘” is the Hadamard product operator. Suppose g=αfEM+βf,0⩽α,β⩽1, 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 Γ(i) (22) in the beginning of the i-th iteration. Then we (a) compute an interval of possible boundary pixel values [ai,bi] (20); (b) use the interval [ai,bi] 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 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: S={f˜|f˜=P(αfEM+βf),0⩽α,β⩽1}. Note that the image segmented by TFA algorithm can be obtained from the presented method by choosing α=0 and β=1 in (25).

Hybrid of EMS and TFA algorithm

1. Input: given MRA image f.
 While happen Eq. (8),
Step 1  Update ωinew , μinew and σinew , 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 fEM .
Step 2 4. Constitute g=αfEM+βf
5. Set g(0)=g and Γ(0) by Eq. (18)
6. While g(i+1) becomes a binary image
 (a). Obtain [ai,bi] 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 Γ(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) .

The presented method converges to a binary image.

Suppose f is a given image. We obtain fEM=M∘f , where M is a binary image that is obtained by Eq. (10). We constitute a new image g=αfEM+βf with linear combination of fEM and f, where 0⩽α , β⩽1 . According to theorem 2, iteratively TFA algorithm converges to a binary image g˜=P(g) , where 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(KNd3) . Moreover, the M-step requires O(KNd2) 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.

 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 λ=0.1 and δ=0.0001 . ε used in (18) is chosen to be ε=0.003 except in Example 1 in which we set ε=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).

The number of pixels in the set Γ(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 Γ(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 Γ(i) at each iteration shows the number of pixels that are unclassified yet in the i-th iteration.

 Discussion

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.

This example is a 512×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 ε=0.02 in (18). We see that in these figures, the EMCTFA method extracts the vessels better than the other mentioned methods.

The number of the wrongly-detected pixels are given in Fig. 4 for noisy simulated images by additive Gaussian noise with mean 0.01 and variances 0.001, 0.0015, 0.002 and 0.0025. We see in Fig. 4 that the EMCTFA method has the least wrongly detected vessel pixels.

The number of the wrongly detected pixels.

This example is a 128×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 α=0.25 and β=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 Γ(i) at each iteration of the presented method is less than that of the TFA algorithm.

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.

This example is a 256×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 α=0.2 and β=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.

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.

This example is one slide of 3D CE-MRA image of the abdominal vascular system whose size is 512×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 α=1 and β=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.

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).

This example is a 512×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 α=1 and β=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.

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).

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 ( α=0 , β=1 ), ( α=0.5 , β=1 ), ( α=0.75 , β=1 ), ( α=1 , β=1 ), ( α=1 , β=.25 ) and ( α=1 , β=0 ), respectively (red lines).

Example 6. Results of vessel segmentation for image (a), obtained with the presented method for (b) α=0 , β=1 ; (c) α=0.5 , β=1 ; (d) α=0.75 , β=1 ; (e) α=1 , β=1 ; (f) α=1 , β=0.25 ; (g) α=1 , β=0 . Blue lines denote boundary of vessel extraction.

 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.

 References Aghazadeh, N., Cigaroudy, L.S. (2014). A multistep segmentation algorithm for vessel extraction in medical imaging. arXiv:1412.8656. Arivazhagan, S., Ganesan, L. (2003). Texture segmentation using wavelet transform. Pattern Recognition Letters, 24(16), 3197–3203. Budak, Ü., Cömert, Z., Çıbuk, M., Şengür, A. (2020). DCCMED-Net: densely connected and concatenated multi Encoder-Decoder CNNs for retinal vessel extraction from fundus images. Medical Hypotheses, 134, 109426. Cai, J.-F., Chan, R., Shen, L., Shen, Z. (2008). Restoration of chopped and nodded images by framelets. SIAM Journal on Scientific Computing, 30(3), 1205–1227. Cai, X., Chan, R., Morigi, S., Sgallari, F. (2013). Vessel segmentation in medical imaging using a tight-frame–based algorithm. SIAM Journal on Imaging Sciences, 6(1), 464–486. Candès, E., Demanet, L., Donoho, D., Ying, L. (2006). Fast discrete curvelet transforms. Multiscal e Modeling & Simulation, 5(3), 861–899. Caselles, V., Catté, F., Coll, T., Dibos, F. (1993). A geometric model for active contours in image processing. Numerische mathematik, 66(1), 1–31. Chan, R.H., Chan, T.F., Shen, L., Shen, Z. (2003). Wavelet algorithms for high-resolution image reconstruction. SIAM Journal on Scientific Computing, 24(4), 1408–1432. Chan, T.F., Vese, L.A. (2001). Active contours without edges. IEEE Transactions on Image Processing, 10(2), 266–277. Cigaroudy, L.S., Aghazadeh, N. (2017). A multiphase segmentation method based on binary segmentation method for Gaussian noisy image. Signal, Image and Video Processing, 11(5), 825–831. Cline, H.E. (2000). Enhanced visualization of weak image sources in the vicinity of dominant sources. Google Patents. US Patent 6,058,218. Daubechies, I., Bates, B.J. (1993). Ten lectures on wavelets. The Journal of the Acoustical Society of America, 93(3), 1671. Dempster, A.P., Laird, N.M., Rubin, D.B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1), 1–22. Do, M.N., Vetterli, M. (2005). The contourlet transform: an efficient directional multiresolution image representation. IEEE Transactions on Image Processing, 14(12), 2091–2106. Dougherty, G. (2011). Medical Image Processing. Techniques and Applications. Springer Science & Business Media. Duffin, R.J., Schaeffer, A.C. (1952). A class of nonharmonic Fourier series. Transactions of the American Mathematical Society, 72(2), 341–366. Fajardo, V.A., Liang, J. (2017). On the EM-Tau algorithm: a new EM-style algorithm with partial E-steps. arXiv:1711.07814. Firoiu, I. (2010). Complex Wavelet Transform. Application to Denoising Universitatea Politechnica, Timisoara. Frangi, A.F., Niessen, W.J., Vincken, K.L., Viergever, M.A. (1998). Multiscale vessel enhancement filtering. In: International Conference on Medical Image Computing and Computer-Assisted InterveNTION, Springer, pp. 130–137. Gerig, G., Kikinis, R., Jolesz, F.A. (1990). Image processing of routine spin-echo MR images to enhance vascular structures: comparison with MR angiography. In: 3D Imaging in Medicine. Springer, pp. 121–132. Guo, K., Labate, D. (2007). Optimally sparse multidimensional representation using shearlets. SIAM Journal on Mathematical Analysis, 39(1), 298–318. Hashemzadeh, M., Azar, B.A. (2019). Retinal blood vessel extraction employing effective image features and combination of supervised and unsupervised machine learning methods. Artificial Intelligence in Medicine, 95, 1–15. Iwasaki, T., Ritman, E.L., Fiksen-Olsen, M.J., Romero, J.C., Knox, F.G. (1985). Renal cortical perfusion-preliminary experience with the dynamic spatial reconstructor (DSR). Annals of Biomedical Engineering, 13(3–4), 259–271. Kakileti, S.T., Venkataramani, K. (2016). Automated blood vessel extraction in two-dimensional breast thermography. In: 2016 IEEE International Conference on Image Processing (ICIP). IEEE, pp. 380–384. Kingsbury, N. (2001). Complex wavelets for shift invariant analysis and filtering of signals. Applied and Computational Harmonic Analysis, 10(3), 234–253. Kirbas, C., Quek, F. (2004). A review of vessel extraction techniques and algorithms. ACM Computing Surveys (CSUR), 36(2), 81–121. Krissian, K., Malandain, G., Ayache, N. (1997). Directional anisotropic diffusion applied to segmentation of vessels in 3D images. In: International Conference on Scale-Space Theories in Computer Vision, Springer, pp. 345–348. Krissian, K., Malandain, G., Ayache, N., Vaillant, R., Trousset, Y. (1998). Model based multiscale detection of 3D vessels. In: Proceedings. Workshop on Biomedical Image Analysis (Cat. No. 98EX162). IEEE, pp. 202–210 Kutyniok, G., Labate, D. (2012). Shearlets: Multiscale Analysis for Multivariate Data. Springer Science & Business Media. Labate, D., Lim, W.-Q., Kutyniok, G., Weiss, G. (2005). Sparse multidimensional representation using shearlets. In: Proceedings of SPIE 5914, Wavelets XI, 59140U, International Society for Optics and Photonics. Lorenz, C., Carlsen, I.-C., Buzug, T.M., Fassnacht, C., Weese, J. (1997). A multi-scale line filter with automatic scale selection based on the Hessian matrix for medical image segmentation. In: International Conference on Scale-Space Theories in Computer Vision. Springer, pp. 152–163. Malladi, R., Sethian, J.A., Vemuri, B.C. (1995). Shape modeling with front propagation: a level set approach. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17(2), 158–175. Mallat, S. (1999). A Wavelet Tour of Signal Processing. Elsevier. Mustafa, W.A.B.W., Yazid, H., Yaacob, S.B., Basah, S.N.B. (2014). Blood vessel extraction using morphological operation for diabetic retinopathy. In: 2014 IEEE Region 10 Symposium. IEEE, pp. 208–212. Navid, M., Hamidpour, S.S.F., Khajeh-Khalili, F., Alidoosti, M. (2020). A novel method to infrared thermal images vessel extraction based on fractal dimension. Infrared Physics & Technology, 107, 103297. Prinet, V., Monga, O., Ge, C., Xie, S., Ma, S. (1996). Thin network extraction in 3D images: application to medical angiograms. In: Proceedings of 13th International Conference on Pattern Recognition, Vol. 3. IEEE, pp. 386–390. Ron, A., Shen, Z. (1997). Affine systems inL2 (Rd): the analysis of the analysis operator. Journal of Functional Analysis, 148(2), 408–447. Rothwell, P.M., Eliasziw, M., Gutnikov, S., Fox, A.J., Taylor, D.W., Mayberg, M., Warlow, C.P., Barnett, H., Collaboration, C.E.T., et al.(2003). Analysis of pooled data from the randomised controlled trials of endarterectomy for symptomatic carotid stenosis. The Lancet, 361(9352), 107–116. Saha, P.K., Udupa, J.K., Odhner, D. (2000). Scale-based fuzzy connected image segmentation: theory, algorithms, and validation. Computer Vision and Image Understanding, 77(2), 145–174. Sato, Y., Nakajima, S., Shiraga, N., Atsumi, H., Yoshida, S., Koller, T., Gerig, G., Kikinis, R. (1998). Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images. Medical Image Analysis, 2(2), 143–168. Selesnick, I.W., Baraniuk, R.G., Kingsbury, N.C. (2005). The dual-tree complex wavelet transform. IEEE Signal Processing Magazine, 22(6), 123–151. Suri, J.S., Laxminarayan, S. (2003). Angiography and Plaque Imaging: Advanced Segmentation Techniques. CRC Press. Udupa, J.K., Samarasekera, S. (1996). Fuzzy connectedness and object definition: theory, algorithms, and applications in image segmentation. Graphical Models and Image Processing, 58(3), 246–261. Udupa, J.K., Odhner, D., Tian, J., Holland, G., Axel, L. (1997). Automatic clutter-free volume rendering for MR angiography using fuzzy connectedness. In: Medical Imaging 1997: Image Processing, Vol. 3034. International Society for Optics and Photonics. pp. 114–119. Unser, M. (1995). Texture classification and segmentation using wavelet frames. IEEE Transactions on Image Processing, 4(11), 1549–1560. Wells, W.M., Grimson, W.E.L., Kikinis, R., Jolesz, F.A. (1996). Adaptive segmentation of MRI data. IEEE Transactions on Medical Imaging, 15(4), 429–442. Wilson, D.L., Noble, J.A. (1997). Segmentation of cerebral vessels and aneurysms from MR angiography data. In: Biennial International Conference on Information Processing in Medical Imaging. Springer, pp. 423–428. Wilson, D.L., Noble, J.A. (1999). An adaptive segmentation algorithm for time-of-flight MRA data. IEEE Transactions on Medical Imaging, 18(10), 938–945.