152	M.A. Dabbah et al.

The dual-model consists of a 2D Gabor wavelet (foreground model) and a Gaussian envelope (background model), which are applied to the original CCM images. The detection relies on estimating the correct local and dominant orientation of the nerve fibers. We evaluate our dual-model in comparison with feature detectors described in Sect. 7.2 that are well established for linear and more general image features. In addition to the evaluation of the nerve fiber detection responses, we have also evaluated the clinical utility of the method by a comparison with manual analysis.

7.4.1 Foreground and Background Adaptive Models

For this purpose, the foreground model MF is an even-symmetric and real-valued Gabor [9,38] wavelet and the background model MB is a two-dimensional Gaussian envelope, Fig. 7.4.

2π

1

x2

γ 2y2

MF(xθ , yθ ) = cos

xθ + ϕ

exp −

θ

+

θ

(7.1)

λ

2

σx2

σy2

1

x2

γ 2y2

MB(xθ , yθ ) = α exp

θ

+

θ

(7.2)

−2

σx2

σy2

xθ = x cos θ + ysin θ

(7.3)

yθ = −x sin θ + ycosθ

(7.4)

The x and y axes of the dual-model coordinate frame xθ and yθ are defined by a rotation of θ , which is the dominant orientation of the nerve fibers in a particular region within the image (see Sect. 7.4.2). This dual-model is used to generate the positive response RP = MF + MB and the negative response RN = MF − MB that are applied to the original CCM image and can be represented as in (7.5) and (7.6), respectively.

2π

1

x2

γ 2y2

RP(xθ , yθ ) =

cos

xθ + ϕ

+ α

exp

−

θ

+

θ

(7.5)

λ

2

σx2

σy2

2π

1

x2

γ 2y2

RN(xθ , yθ ) =

cos

xθ + ϕ

− α

exp

−

θ

+

θ

(7.6)

λ

2

σx2

σy2

The equations of RP and RN assume that the Gaussian envelope of both responses are identical, that is, they have the same variances σ 2(x, y) and the same aspect ratio γ . The magnitude of the Gaussian envelope α defines the threshold in which a

7 Detecting and Analyzing Linear Structures in Biomedical Images: A Case Study...

153

Fig. 7.4 Foreground and background models for the nerve fibers. (a) the two-dimensional Gabor wavelet at a particular orientation and frequency. It represents the foreground model of the nerve fibers. (b) the Fourier transforms of (a). (c) the two-dimensional Gaussian envelope that represents the background model and (d) its Fourier transform

nerve fiber can be distinguished from the background image. The value of α can be set empirically to control sensitivity and accuracy of detection. The wavelength λ defines the frequency band of the information to be detected in the CCM image, and is related to the width of the nerve fibers (see Fig. 7.2). Its value might be computed for a sub-region within the image that has significant variability of nerve fiber width. However, for simplicity, λ is chosen to be a global estimate of the entire image based on empirical results.

7.4.2 Local Orientation and Parameter Estimation

In CCM images, the nerve fibers flow in locally constant orientations. In addition, there is a global orientation that dominates the general flow. The orientation field describes the coarse structure of nerve fibers in the CCM images and has been proven to be of a fundamental importance in many image analysis applications [39, 40]. Using the least mean square algorithm [41], the local orientation θ (i, j) of the block centered at pixel (i, j) (7.9), is computed using the following equations [39].

154									M.A. Dabbah et al.
ω			ω
Vx(i, j) = ∑iu+=i2	ω	∑vj=+ j2			ω	∂x2(u, v) − ∂y2(u, v)			(7.7)
−	2			−	2
Vy(i, j) = ∑i+	ω	ω ∑j+			ω	ω 2∂x(u, v)∂y(u, v)
	2				2				(7.8)
u=i−		2	v= j−			2
θ (i, j) = π /2 +			1	tan−1			Vy(i, j)		(7.9)
			2
						Vx(i, j)

The gradients ∂x(u, v) and ∂y(u, v) are computed at each pixel (u, v) and may vary from the simple Sobel operator to the more complex Canny operator depending on the computational requirements. ω is the width of the block centered at pixel (i, j). The orientation field is then smoothed by convolving the x and y vector field components in (7.7) and (7.8), respectively, with a low-pass Gaussian filter. This

smoothed orientation field is calculated by (7.14), where Φx(i, j) the smoothed continuous x and y vector field components.

Φx(i, j) = cos(2θ (i, j))

Φy(i, j) = sin (2θ (i, j))

and Φy(i, j) are

(7.10)

(7.11)

According to the original algorithm [41], the low-pass 2-dimensional Gaussian filter G is applied on the block level ω of the orientation field computed earlier in (7.9). The filter has a unit integral and a kernel size of ωΦ × ωΦ . However, since the orientation in CCM images varies at a slow rate, the low-pass filter is applied globally to further reduce errors at near-nerve fiber and nonnerve fiber regions. The estimated orientation is not always correct, hence, the low-pass filter tries to rectify the error given that the orientation in the local neighborhood varies slowly;

ωΦ		ωΦ
	ωΦ ∑ 2
Φx(i, j) = ∑ 2			ωΦ
u=−	2	v=−	2
ωΦ		ωΦ
	ωΦ ∑ 2
Φy(i, j) = ∑ 2			ωΦ
u=−	2	v=−	2

O(i, j) = 12 tan−1

G(u, v)Φx(i − u, j − v)			(7.12)
G(u, v)Φy(i − u, j − v)			(7.13)
	Φy(i, j)		(7.14)
	Φx(i, j)

The least square estimate produces a stable smooth orientation field in the region of the nerve fibers. However, when applied on the background of the image, that is, between fibers, the estimate is dominated by noise due to the lack of structure and uniform direction, which is expected and understandable. Figure 7.5 shows a CCM image and its orientation field estimate.

7 Detecting and Analyzing Linear Structures in Biomedical Images: A Case Study...

155

Fig. 7.5 An illustration of the orientation field (right) of the original CCM image (left). The orientations on the nerve fibers and their surrounding are similar and follow the predominant orientation in the image, while orientations everywhere else (background) are random and noisy

7.4.3 Separation of Nerve Fiber and Background Responses

The models are applied on the image pixel-wise. During this operation, they are adjusted to suit the local neighborhood characteristics of the reference pixel at f (i, j) by modifying their parameters of the foreground and background separately in (7.5) and (7.6). The dot products of the models and the reference pixel’s neighborhood ((7.15) and (7.16)) are then combined to generate the final enhanced value of this particular reference pixel g(i, j) (7.17).

Γp(i, j) = fω (i, j), RP		(7.15)
Γn(i, j) = fω (i, j), RN		(7.16)
g(i, j) =	Γ p(i, j)	(7.17)
	1 + exp{−2kΓn(i, j)}

The neighborhood area of the reference pixel is defined by the width ω . The transition from foreground to background at a particular pixel g(i, j) occurs at Γn = 0. The sharpness of this transition is controlled by k: larger k results in sharper transition. This in turn enhances the nerve fibers that are oriented in the dominant direction, and decreases noisy structures that are oriented differently by increasing the contrast between the foreground and the noisy background, whilst effectively reducing noise around the nerve fiber structure as shown in Fig. 7.6.