**4.2.1 Retinal vasculature extraction using canny edge detector**

The Canny operator [20] is less likely than the others to be "fooled" by noise, and more likely to detect true weak edges. The Canny operator preserves most edges among all other edge detectors. Therefore, the Canny edge detector is employed in this research to extract the retinal vasculature edges. There are two criteria used in the Canny Operator to locate the rapidly changed intensity pixels. They are:

1. Pixels where the 1st derivative of the intensity is larger in magnitude than a certain threshold;

Image 1 Image 2 Image n

Image Segmentation

Control Point Matching

Objective Function Maximization; Iteratively Adjusted Control Point Pairs

Final Fused Image

**4.2.1 Retinal vasculature extraction using canny edge detector** 

Generation Tools:

The Canny operator [20] is less likely than the others to be "fooled" by noise, and more likely to detect true weak edges. The Canny operator preserves most edges among all other edge detectors. Therefore, the Canny edge detector is employed in this research to extract the retinal vasculature edges. There are two criteria used in the Canny Operator to locate the

1. Pixels where the 1st derivative of the intensity is larger in magnitude than a certain

Tools:

Control Point Detection Tools: Adaptive Exploratory

Tools:

Tools:

Edge Detector

Algorithm

Shape Similarity

Heuristic Optimization Algorithm

Mathematical Transformation

Modal

Criteria

Edge Detection

Input: **… …** 

Image Segmentation:

Feature Extraction:

Feature Matching:

Optimization:

Output - Fusion:

threshold;

Diagram 1. Image fusion procedure.

rapidly changed intensity pixels. They are:


$$\frac{d^2\left(G\times I\right)}{dn^2} = \frac{d\left(\left(\frac{dG}{dn}\right)\times I\right)}{dn} \tag{1}$$

where, *n* is the direction of the gradient of the image; *G* is the edge signal; *I* is the image intensity. The zero-crossings of Canny's method correspond to the first directionalderivative's maxima and minima in the direction of the gradient. Edges will be identified as the maxima in magnitude. Each pixel's edge gradient is computed and compared with the gradients of its neighbors along the gradient direction. If the central pixel is smaller, mark the current edge's intensity as 0; if largest among all neighbors, keep the original intensity. Based on the nine-pixel neighborhood, the normal to the edge direction has two *ux* and *uy*. In order to estimate the gradient on the discrete sampling, two pixels closest to *u* are selected. A plane can be identified by the gradient magnitudes of three pixels. By using this plane, the gradient magnitude and the intensity at each pixel on the line can be locally estimated. The gradient magnitude at *Px+1, y+1* and *P x-1, y-1* (Figure 18) can be calculated as:

$$P\_{\mathbf{x}^{y}=1,\mathbf{y}^{y}=1}. \qquad \qquad \mathbf{G}(P\_{\mathbf{x}+1,\mathbf{y}+1}) = \frac{\mu}{\mu\_{\mathbf{y}}} \mathbf{G}(\mathbf{x}+1,\mathbf{y}+1) + \frac{\mu\_{\mathbf{y}}-\mu\_{\mathbf{x}}}{\mu\_{\mathbf{y}}} \mathbf{G}(\mathbf{x},\mathbf{y}+1) \tag{2}$$

$$P\_{x=1,y=1} \colon \qquad \mathbf{G}(P\_{x=1,y=1}) = \frac{\mu}{u\_y} \mathbf{G}(x-1, y-1) + \frac{u\_y - u\_x}{u\_y} \mathbf{G}(x, y-1) \tag{3}$$

If the gradient at *Px,y* is greater than both of 1, 1 ( ) *G Px y* + + and 1, 1 ( ) *G Px y* − − , *Px,y* will be identified as a maximum.

Fig. 18. Canny Edge Detection – Localization of Maxima [34].

In order to make the localization of magnitude maxima accurate, Canny defined a filter by optimizing a performance index which enhances real positive and real negative. The filter is used to minimize the probability of non-detected edge points and false detection.

$$SNR = \frac{\left| \int^{w} G(-\mathbf{x}) f(\mathbf{x}) d\mathbf{x} \right|}{n\_0 \sqrt{\int^{w} f^2(\mathbf{x}) d\mathbf{x}}} \tag{4}$$

where, SNR stands for Signal-to-Noise Ratio; *f* is the filter; denominator is the RMSE response to noise *n(x)*. The identification of the real edge localization is defined as:

$$\text{Localization} = \frac{1}{\sqrt{\mathbb{E}[\mathbf{x}\_0^2]}} = \frac{\left| \int\_{-w}^{w} G'(\mathbf{x}) f'(\mathbf{x}) d\mathbf{x} \right|}{n\_0 \sqrt{\int\_{-w}^{w} f'^2(\mathbf{x}) d\mathbf{x}}} \tag{5}$$

Two adaptive thresholds are used in Canny's method. They are high threshold and low threshold. The high threshold is used to find the start point of strong edges. Any points that meet the high threshold will be selected as the edge point. These start points are growing into different directions as long as there is no edge strength falling below the low threshold. Figure 19 is the 3D shaded surface plots of the original retinal angiogram image. The *X-Y*  axis corresponds to the original image size. The height *Z* axis is a single-valued function defined over a geometrically rectangular grid. *Z* specifies the color data as well as surface height, so color is proportional to surface height with range of [0, 1] of each pixel on the image. All the retinal salient features are preserved in the Canny edges. Those salient features are the retinal vessel bifurcations, from which the control points will be selected using the Adaptive Exploratory Algorithm. Figure 20 and 21 show the retinal vessel edges detected by the Canny operator.

Fig. 19. 3D shaded surface plot of the angiogram image.

## **4.2.2 Control point detection**

462 Biomedical Science, Engineering and Technology

0

*w w*

−

= ∫

response to noise *n(x)*. The identification of the real edge localization is defined as:

*SNR*

*Localization*

detected by the Canny operator.

**Z** 

**Y** 

Fig. 19. 3D shaded surface plot of the angiogram image.

*w*

2

*n f x dx*

*w*

∫

−

−

where, SNR stands for Signal-to-Noise Ratio; *f* is the filter; denominator is the RMSE

2

Two adaptive thresholds are used in Canny's method. They are high threshold and low threshold. The high threshold is used to find the start point of strong edges. Any points that meet the high threshold will be selected as the edge point. These start points are growing into different directions as long as there is no edge strength falling below the low threshold. Figure 19 is the 3D shaded surface plots of the original retinal angiogram image. The *X-Y*  axis corresponds to the original image size. The height *Z* axis is a single-valued function defined over a geometrically rectangular grid. *Z* specifies the color data as well as surface height, so color is proportional to surface height with range of [0, 1] of each pixel on the image. All the retinal salient features are preserved in the Canny edges. Those salient features are the retinal vessel bifurcations, from which the control points will be selected using the Adaptive Exploratory Algorithm. Figure 20 and 21 show the retinal vessel edges

1

= =

*E x*

( )

'

*w*

∫

−

*w*

−

*n f x dx*

( ) '( )

*G x f x dx*

0 '2 0

*w*

−

*w*

∫

[ ] ( )

(4)

(5)

**X** 

( )()

*G x f x dx*

A good-guess of the initial control point selection ensures fused image generated at an efficient computational time. Bad control point selection will significantly increase the computation cost, or even cause the image fusion fail. Vessels or some particular abnormalities make images not necessarily matching the retina structures. Even when structure and function correspond, the abnormality still happens sometimes if inconsistence exists between structural and functional changes. Further more, angiogram images usually have higher resolution and are rich in information, whereas fundus images have lower resolution and are indeed abstract with some details or even missing some small vessels. Practically, those situations are unavoidable and will create difficulties in extracting the control points because the delineation of the vein boundaries may not be precise. In this study, control points are detected using the adaptive exploratory algorithm (Figure 20 and 21) [21].

Fig. 20. Angiogram grayscale reference image's control point selection.

Fig. 21. Fundus true color input image's control point selection.

### **4.2.3 Heuristic optimization algorithm**

An optimization procedure is required to adjust the initial good-guess control points in order to achieve the optimal result. The process can be formulated as a heuristic problem of optimizing an objective function that maximizes the Mutual-Pixel-Count between the reference and input images. The algorithm finds the optimal solution by refining the transformation parameters in an ordered way. By maximizing the objective function, one image's vessels are supposed to be well overlaid onto those of the other image (Figure 23). Mutual-Pixel-Count measures the optic nerve head vasculature overlapping for corresponding pixels in both images. It is assumed that the retinal vessels are represented by 0 (black pixel) and background is represented by 1 (white pixels) in the binary 2D map. When the vasculature pixel's transformed *(u, v)* coordinates on the input image correspond to the vasculature pixel's coordinates on the reference image, the MPC is incremented by 1 (Figure 22). MPC is assumed be maximized when the image pair is perfectly geometrically aligned by the transformation. After pre-processing, the binary images of the reference and input images are obtained, i.e. *Iref* and *Iinput*. Only black pixels from both images contribute to MPC. The ideal case is that all zero pixels of the input image are mapped onto zero pixels of the reference image. The problem can be mathematically formulated as the maximization of the following objective function:

$$f\_{mpc}(\mathbf{x}, y, \mu, \upsilon) = \sum\_{\substack{\mathbf{l}, \upsilon \in \mathcal{R} \cup \mathcal{I} \\ \mathsf{l}\_{\mathsf{r}\mathbf{f}} \{\mathsf{x}, \mathsf{y}\} = 1}}^{\mathsf{u}, \upsilon \in \mathcal{R} \cup \mathcal{I}} I\_{input}(T\_{\mathbf{x}}(\mathsf{u}, \upsilon), T\_{\mathbf{y}}(\mathsf{u}, \upsilon)) \tag{6}$$

where *fmpc* denotes the value of the Mutual-Pixel-Count. *Tx* and *Ty* are the transformations for u and v coordinates of the input image. The ROI (Region-of-Interest) is the vasculature region where the MPC is calculated on.

Coordinates' adjustment is iteratively implemented until one of the following convergence criteria is reached: (1). Predefined maximum number of loops is reached; or (2). the updated *fMPC* is smaller than ε, i.e.

$$\left| \mathbf{f}\_{\text{MPC}^{n+1}} \left( \mathbf{x}, \mathbf{y}, \mathbf{u}, \mathbf{v} \right) \ast \mathbf{f}\_{\text{MPC}^n} \left( \mathbf{x}, \mathbf{y}, \mathbf{u}, \mathbf{v} \right) \right| < \varepsilon \tag{7}$$

where εis a very small non-negative threshold.

Fig. 22. fMPC increasing during the iteration (Y-axis is fMPC; X-axis is the loop count).

An optimization procedure is required to adjust the initial good-guess control points in order to achieve the optimal result. The process can be formulated as a heuristic problem of optimizing an objective function that maximizes the Mutual-Pixel-Count between the reference and input images. The algorithm finds the optimal solution by refining the transformation parameters in an ordered way. By maximizing the objective function, one image's vessels are supposed to be well overlaid onto those of the other image (Figure 23). Mutual-Pixel-Count measures the optic nerve head vasculature overlapping for corresponding pixels in both images. It is assumed that the retinal vessels are represented by 0 (black pixel) and background is represented by 1 (white pixels) in the binary 2D map. When the vasculature pixel's transformed *(u, v)* coordinates on the input image correspond to the vasculature pixel's coordinates on the reference image, the MPC is incremented by 1 (Figure 22). MPC is assumed be maximized when the image pair is perfectly geometrically aligned by the transformation. After pre-processing, the binary images of the reference and input images are obtained, i.e. *Iref* and *Iinput*. Only black pixels from both images contribute to MPC. The ideal case is that all zero pixels of the input image are mapped onto zero pixels of the reference image. The problem can be mathematically formulated as the maximization of

,

*u v ROI mpc input x y I x y and I u v f x y u v I T uv T uv* ∈

*ref input*

(,)1 (,)1

= =

where *fmpc* denotes the value of the Mutual-Pixel-Count. *Tx* and *Ty* are the transformations for u and v coordinates of the input image. The ROI (Region-of-Interest) is the vasculature

Coordinates' adjustment is iteratively implemented until one of the following convergence criteria is reached: (1). Predefined maximum number of loops is reached; or (2). the updated

> f x,y,u,v - f x,y,u,v MPCn+1 ( )( ) MPC<sup>n</sup> <

Fig. 22. fMPC increasing during the iteration (Y-axis is fMPC; X-axis is the loop count).

(,,,) ( ( , ), ( , ))

<sup>=</sup> ∑ (6)

ε

(7)

**4.2.3 Heuristic optimization algorithm** 

the following objective function:

region where the MPC is calculated on.

ε, i.e.

is a very small non-negative threshold.

*fMPC* is smaller than

where ε

(e)

Fig. 23. Fused image improvement during the iteration. From (a)-(e) fMPC = 5144, 7396, 7484, 7681, 7732.

### **4.2.4 Affine transformation model**

A mathematical model is the tool for transforming the target images and fusing them into one single volume. Affine model is a basic geometrical transformation in image processing and is defined as Eq. 8 and 9. The DOF of the affine model is 6 because it has six parameters, i.e. *a1*, *a2, b1, b2, a3*, and *a4*.

$$
\begin{pmatrix} \mathcal{U} \\ \mathcal{V} \\ \mathcal{V} \end{pmatrix} = \begin{pmatrix} a\_1 & a\_2 & b\_1 \\ a\_3 & a\_4 & b\_2 \\ a\_2 & a\_3 & b\_3 \end{pmatrix} \begin{pmatrix} \mathcal{x} \\ \mathcal{y} \\ \mathcal{z} \end{pmatrix} \\
\begin{array}{ll} \mathcal{U}(\mathcal{x}, \mathcal{y}) = a\_1\mathcal{x} + a\_2\mathcal{y} + b\_1 \\ \mathcal{V}(\mathcal{x}, \mathcal{y}) = a\_3\mathcal{x} + a\_4\mathcal{y} + b\_2 \\ \end{array} \\ \tag{8}
$$

$$
\begin{pmatrix} V \\ 1 \end{pmatrix} = \begin{vmatrix} a\_3 & a\_4 & b\_2 \\ 0 & 0 & 1 \end{vmatrix} \begin{pmatrix} y \\ 1 \end{pmatrix} \quad V(x, y) = a\_3 x + a\_4 y + b\_2 \tag{9}
$$

Affine model's advantage lies in that it can measure lost information such as skew, translation, rotation, shearing and scaling that maps finite points to finite points and parallel lines to parallel lines (Figure 24 and Table 1). Its drawback lies in the strict requirement that at least six pairs of control points are needed [19].

Fig. 24. Transformation Graphs.


Table 1. Transformation Descriptions.
