Hybrid Vessel Extraction Method Based on Tight-Frame and EM Algorithms by Using 2D Dual Tree Complex Wavelet

1Introduction

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.

2Methods

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

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 {ξ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):

(1)

F=18Iid−IidIidIidIidIidIid−IidFξξFηηFηξFξη.

The inverse of F is obtained as follows:

(2)

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.

Fig. 2

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

3Results

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

4Discussion

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

Fig. 4

The number of the wrongly detected pixels.

Example 2.

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.

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

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

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

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

Fig. 9

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.

5Conclusions

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.