# Analysis and Implementation of Kidney Stone Detection by Reaction Diffusion Level Set Segmentation Using Xilinx System Generator on FPGA

Authors: Kalannagari Viswanath; Ramalingam Gunasundari
Journal: VLSI Design
 (2015)

---

## Abstract

Ultrasound imaging is one of the available imaging techniques used for diagnosis of kidney abnormalities, which may be like change in shape and position and swelling of limb; there are also other Kidney abnormalities such as formation of stones, cysts, blockage of urine, congenital anomalies, and cancerous cells. During surgical processes it is vital to recognize the true and precise location of kidney stone. The detection of kidney stones using ultrasound imaging is a highly challenging task as they are of low contrast and contain speckle noise. This challenge is overcome by employing suitable image processing techniques. The ultrasound image is first preprocessed to get rid of speckle noise using the image restoration process. The restored image is smoothened using Gabor filter and the subsequent image is enhanced by histogram equalization. The preprocessed image is achieved with level set segmentation to detect the stone region. Segmentation process is employed twice for getting better results; first to segment kidney portion and then to segment the stone portion, respectively. In this work, the level set segmentation uses two terms, namely, momentum and resilient propagation (R
prop) to detect the stone portion. After segmentation, the extracted region of the kidney stone is given to Symlets, Biorthogonal (bio3.7, bio3.9, and bio4.4), and Daubechies lifting scheme wavelet subbands to extract energy levels. These energy levels provide evidence about presence of stone, by comparing them with that of the normal energy levels. They are trained by multilayer perceptron (MLP) and back propagation (BP) ANN to classify and its type of stone with an accuracy of 98.8%. The prosed work is designed and real time is implemented on both Filed Programmable Gate Array Vertex-2Pro FPGA using Xilinx System Generator (XSG) Verilog and Matlab 2012a.


---

## Body

## 1. Introduction

Kidney stone disease is one of the major life threatening ailments persisting worldwide. The stone diseases remain unnoticed in the initial stage, which in turn damages the kidney as they develop. A majority of people are affected by kidney failure due to diabetes mellitus, hypertension, glomerulonephritis, and so forth. Since kidney malfunctioning can be menacing, diagnosis of the problem in the initial stages is advisable. Ultrasound (US) image is one of the currently available methods with noninvasive low cost and widely used imaging techniques for analyzing kidney diseases [1]. Shock wave lithotripsy (SWL), percutaneous nephrolithotomy (PCNL), and relative super saturation (RSS) are the available practices to test urine. The Robertson Risk Factor Algorithms (RRFA) are open and are used for laparoscopic surgery; these algorithms are assigned for exceptional [2] special cases. Hyaluronan is a large (>106 Da) linear glycosaminoglycan composed of repeating units of glucuronic acid (GlcUA) and N-acetyl glucosamine (GlcNAc) disaccharides [3]. It has a significant role in a number of processes that can eventually lead to renal stone disease, including urine concentration, uric acid, salt form crystal, crystallization inhibition, crystal retention, magnesium ammonium phosphate, and amino acid.Rahman and Uddin have proposed diminution of speckle noise and segmentation from US image. It not only detects problem in the kidney region but also provides image quality enhancement [1]. Hafizah proposed kidney US images and divides them into four dissimilar categories: normal, bacterial infection, cystic disease, and kidney stones, using gray level cooccurrence matrix (GLCM). This categorization helps doctors to identify the abnormalities in kidney [4]. Rathi and Palani have proposed a Hierarchical Self-Organizing Map (HSOM) for brain tumours using the segmentation technique and wavelets packets. Accuracy of the results was found to be correct up to 97% [5]. Norihiro Koizumi has proposed high intensity focused ultrasound (HIFU) technique for terminating tumours and stones [6, 7]. Viswanath and Gunasundari propose content descriptive multiple stones detection using level set segmentation, wavelets processing for identification of kidney stone, and artificial neural network (ANN) for classification. The results yielded a maximum accuracy of 98.66% [8]. The MLP-BP ANN is found to perform better in terms of accuracy 92% with a speed of 0.44 sec and it is found to be very sensitive [9, 10]. The noninvasive combination of renal using pulsed cavitation US therapy proposed that shock wave lithotripsy (ESWL) has become a customary for the treatment of calculi located in the kidney and ureter [11]. Tamilselvi and Thangaraj have proposed seeded region growing based on segmentation and classification of kidney images with stone sizes using CAD system [12]. Bagley et al. estimate location of urinary stones with unenhanced computed tomography (CT) using half-radiation (low) dose compared with the standard dose. Out of the 50 patients examined, 35 patients were found to have a single stone while the rest had multiple stones [13]. The solution for local minima and segmentation problem was proposed by Thord Andersson, Gunnar Lathen, with modified gradient search and level set segmentation technique [2]. Templates based technique was proposed by Emmanouil Skounakis for 3D detection of kidneys and their pathology in real time. Its accuracy was found to be 97.2% and abnormalities in kidneys had an accuracy of 96.1% [7]. Gabor function is used for achievement of optimal sharpening and smoothening of 2D image in both time and frequency resolutions [14]. Chen et al. have proposed the finite element method based 3D tumor growth prediction using longitudinal kidney tumor images [15]. Using the linear elastic theory Owen et al. proposed pressure finding in fluid for calculating the depth of shock wave scattering by kidney stone in water [16]. The pH value based prediction of stone formation epidemiologically has been proposed by Kok [17]. Datar has proposed the segmentation of the desired portion using initial seed selection, growing, and region merging without any edge detection [18]. Multilayer perceptron and back propagation implementation on FPGA and ASIC design were carried out by Raj and Pinjare [19].This research paper ensues as follows. In Section2 problem statement is demarcated, Section 3 defines the proposed method, in Section 4 image segmentation to locate the kidney stone, in Section 5 optimized energy calculation of segmented portion is discussed, in Section 6 wavelets based energy extraction, Section 7 explains the artificial neural networks classifiers used, in Section 8 experiments results are discussed, and in the last section conclusion of the paper with scope for future work is given.


## 2. Problem Statement

The kidney malfunctioning can be life intimidating. Hence early detection of kidney stone is essential. Precise identification of kidney stone is vital in order to ensure surgical operations success. The ultrasound images of kidney comprise speckle noise and are of low contrast which makes the identification of kidney abnormalities a difficult task. As a result, the doctors may find identification of small stones and the type is difficult and challenging for identify the small kidney stones and their type appropriately. To address this issue, a reaction diffusion level set segmentation is proposed to identify location of the stone; it is implemented in real time on Vertex-2Pro FPGA with Verilog HDL using Xilinx System Generator blocks from Matlab 2012a which is compatible with xilinx13.4 ISE and lifting scheme wavelets subbands are employed for extraction of the energy levels of the stone. The results are analyzed using MLP-BP ANN algorithms for classification and its type of stone [20].


## 3. Methodology

Figure1 shows the overall block diagram of the proposed method. It consists of the following blocks via kidney image database, image preprocessing, image segmentation, wavelet processing, and ANN classification.Figure 1
Proposed block diagram for kidney stone detection.


### 3.1. Kidney Image Database

Kidney image database consists of nearly 500 US kidney images collected from different individuals of various hospitals. It consists of both normal and abnormal images stored in the database. One of the images is selected from the database and subjected to stone detection process.


### 3.2. Image Preprocessing

The aim of preprocessing is to improve the acquired low contrast ultrasound image with speckle noise. It suppresses the undesired distortions and enhances certain image features significant for further processing and stone detection. Without preprocessing, the US image quality may not be good for analyzing. For surgical operations, it is essential to identify the location of kidney stone accurately. Preprocessing helps to overcome this issue of low contrast and speckle noise reduction. Figure2 shows the steps involved in preprocessing of US image, which are as follows:(1)
image restoration,
(2)
smoothing and sharpening,
(3)
contrast enhancement.Figure 2
Preprocessing of kidney image.


#### 3.2.1. Image Restoration

Image restoration is meant to mitigate the degradation of the US image. Degradation may be due to motion blur, noise, and camera misfocus. The main purpose of image restoration is to reduce the degradations that are caused during acquisition of US scanning. In this system, level set function is used for proper orientation. Using plane curve motion, curve smoothers, shrinks are eventually removed [1]. Thus Merriman and Sethian proposed evolution between m
a
x
(
k
,
0
) and m
i
n
(
k
,
0
): (1)
f
x
=
max
⁡
k
,
0
,
if
a
x
,
y
<
G
x
,
y
min
⁡
k
,
0
,
otherwise,where a
(
x
,
y
) is average intensity small neighborhood and G
(
x
,
y
) is median in the same neighborhood.


#### 3.2.2. Smoothing and Sharpening

The restored image is enriched with optimal resolution in both spatial and frequency domains using Gabor filter. This filter acts as a band pass filter with local spatial frequency distribution [6]. Image smoothing and removal of noise is performed using convolution operator. The standard deviation of the Gaussian function can be modified to tune the degree of smoothening and its hardware results are shown in Figure 3.Figure 3
Hardware results of smoothening filter.


#### 3.2.3. Contrast Enhancement

Histogram equalization is employed for improvement of the low contrast US image and achievement of the uniform intensity. This approach can be used on the image as a whole or to a part of an image. In this system, contrast enhancement of the images is executed by transforming the image intensity values, such that the histogram of the output image approximately matches a specified histogram and its results are shown in Figure4. The input and output signal are of same data type.Figure 4
Hardware result of contrast enhancement.


## 3.1. Kidney Image Database

Kidney image database consists of nearly 500 US kidney images collected from different individuals of various hospitals. It consists of both normal and abnormal images stored in the database. One of the images is selected from the database and subjected to stone detection process.


## 3.2. Image Preprocessing

The aim of preprocessing is to improve the acquired low contrast ultrasound image with speckle noise. It suppresses the undesired distortions and enhances certain image features significant for further processing and stone detection. Without preprocessing, the US image quality may not be good for analyzing. For surgical operations, it is essential to identify the location of kidney stone accurately. Preprocessing helps to overcome this issue of low contrast and speckle noise reduction. Figure2 shows the steps involved in preprocessing of US image, which are as follows:(1)
image restoration,
(2)
smoothing and sharpening,
(3)
contrast enhancement.Figure 2
Preprocessing of kidney image.


### 3.2.1. Image Restoration

Image restoration is meant to mitigate the degradation of the US image. Degradation may be due to motion blur, noise, and camera misfocus. The main purpose of image restoration is to reduce the degradations that are caused during acquisition of US scanning. In this system, level set function is used for proper orientation. Using plane curve motion, curve smoothers, shrinks are eventually removed [1]. Thus Merriman and Sethian proposed evolution between m
a
x
(
k
,
0
) and m
i
n
(
k
,
0
): (1)
f
x
=
max
⁡
k
,
0
,
if
a
x
,
y
<
G
x
,
y
min
⁡
k
,
0
,
otherwise,where a
(
x
,
y
) is average intensity small neighborhood and G
(
x
,
y
) is median in the same neighborhood.


### 3.2.2. Smoothing and Sharpening

The restored image is enriched with optimal resolution in both spatial and frequency domains using Gabor filter. This filter acts as a band pass filter with local spatial frequency distribution [6]. Image smoothing and removal of noise is performed using convolution operator. The standard deviation of the Gaussian function can be modified to tune the degree of smoothening and its hardware results are shown in Figure 3.Figure 3
Hardware results of smoothening filter.


### 3.2.3. Contrast Enhancement

Histogram equalization is employed for improvement of the low contrast US image and achievement of the uniform intensity. This approach can be used on the image as a whole or to a part of an image. In this system, contrast enhancement of the images is executed by transforming the image intensity values, such that the histogram of the output image approximately matches a specified histogram and its results are shown in Figure4. The input and output signal are of same data type.Figure 4
Hardware result of contrast enhancement.


## 3.2.1. Image Restoration

Image restoration is meant to mitigate the degradation of the US image. Degradation may be due to motion blur, noise, and camera misfocus. The main purpose of image restoration is to reduce the degradations that are caused during acquisition of US scanning. In this system, level set function is used for proper orientation. Using plane curve motion, curve smoothers, shrinks are eventually removed [1]. Thus Merriman and Sethian proposed evolution between m
a
x
(
k
,
0
) and m
i
n
(
k
,
0
): (1)
f
x
=
max
⁡
k
,
0
,
if
a
x
,
y
<
G
x
,
y
min
⁡
k
,
0
,
otherwise,where a
(
x
,
y
) is average intensity small neighborhood and G
(
x
,
y
) is median in the same neighborhood.


## 3.2.2. Smoothing and Sharpening

The restored image is enriched with optimal resolution in both spatial and frequency domains using Gabor filter. This filter acts as a band pass filter with local spatial frequency distribution [6]. Image smoothing and removal of noise is performed using convolution operator. The standard deviation of the Gaussian function can be modified to tune the degree of smoothening and its hardware results are shown in Figure 3.Figure 3
Hardware results of smoothening filter.


## 3.2.3. Contrast Enhancement

Histogram equalization is employed for improvement of the low contrast US image and achievement of the uniform intensity. This approach can be used on the image as a whole or to a part of an image. In this system, contrast enhancement of the images is executed by transforming the image intensity values, such that the histogram of the output image approximately matches a specified histogram and its results are shown in Figure4. The input and output signal are of same data type.Figure 4
Hardware result of contrast enhancement.


## 4. Image Segmentation

In the segmentation process, five level set methods in all are discussed, all implemented and compared.


### 4.1. Conventional LSS

In the conventional level set method, consider a closed parameterized planar curve or surface, denoted byC
(
y
,
t
)
:
[
0,1
]
x
R
+
→
R
n, where n
=
2 is for planar curve and n
=
3 is for surface, and t is the artificial time generated by the movement of the initial curve or surface C
0
(
y
) in its inward normal direction N
~. The curve or surface evolution equation is as follows [21]:(2)
C
y
,
t
=
0
=
C
0
y
,
C
t
=
F
N
¯
,where F is the force function.In the above equation intrinsic drawback of interactively solving (2) lies in its difficulty to handle topological changes of the moving front, such as splitting and merging [21]. This problem can be eliminated by using the level set method (LSM) by modifying the above equation (2) by taking the derivative with respect to time t on the both sides, yielding the following equation:(3)
ϕ
t
+
∇
ϕ
·
C
t
=
ϕ
t
+
Δ
ϕ
·
F
N
¯
=
0
,
ϕ
x
,
t
=
0
=
ϕ
0
x
,where gradient operator ∇
·
=
∂
·
/
∂
x
1
,
∂
·
/
∂
x
2
,
…
,
∂
(
·
)
/
∂
x
n and ϕ
0
(
x
) is the initial LSF C
0
(
y
)
=
x
∣
ϕ
0
(
x
)
=
0.But (3) fails for too flat or too steep near the zero level set to address this issue re-initialization is introduced.


### 4.2. Reinitialization LSS

But, during evolution, the level set function (LSF) fails for too flat or too steep near the zero level set, causing serious numerical errors. Therefore, a procedure called reinitialization is periodically employed to reshape it to be a signed distance function (SDF). In reinitialization the distance signed function isϕ
(
x
)
=
1
±
d
i
s
t
2
(
x
) where d
i
s
t
(
·
) is a distance function and ± denotes the sign inside and outside the contour [21]. But it has many problems, such as expensive computational cost, blocking the emerging of new contours [21], failures when the LSF deviates much from an SDF, and inconsistency between theory and implementation using XSG shown in Figure 5. Therefore, some formulation has been proposed to regularize the variational LSF to eliminate the reinitialization and computational cost. The following reinitialization equation is given by(4)
ϕ
t
+
S
ϕ
0
∇
ϕ
=
0
,where ϕ
0
=
ϕ
0
/
ϕ
0
2
+
Δ
x
2, ϕ
0 is the initial LSF, and Δ
x is the spatial step.Figure 5
Hardware implementation of signed distance function using XSG.


### 4.3. Distance Regularized Level Set Evolution (DRLSE)

Li et al., proposed a signed distance penalizing energy functional is given by(5)
P
ϕ
1
2
∫
Ω
∇
ϕ
-
1
2
d
x
.Equation (5) measures the closeness between an LSF (ϕ) and an SDF in the domain Ω
⊂
R
n, n
=
2 or 3. By calculus of variation [21], the gradient flow of P
(
ϕ
) is obtained as (6)
ϕ
t
=
-
P
ϕ
ϕ
=
div
⁡
r
1
ϕ
∇
ϕ
.Equation (6) is a diffusion equation with rates (7)
r
1
ϕ
=
1
-
1
∇
ϕ
.However, r
1
(
φ
)
→
-
∞ when ∇
φ
→
0, which may cause oscillation in the final LSF φ. This problem is solved by applying a new diffusion rate(8)
r
2
ϕ
=
sin
⁡
2
∏
∇
ϕ
2
∏
∇
ϕ
,
if
∇
ϕ
≤
1
1
-
1
∇
ϕ
,
if
∇
ϕ
≥
1and a constrained level set diffusion rate as (9)
r
3
ϕ
=
H
p
∇
ϕ
-
1
,where H
p
(
z
)
=
(
1
/
2
)
[
1
+
(
2
/
∏
)
a
r
c
t
a
n
(
z
/
ρ
)
].Andρ is a fixed parameter. The DRLSE methods using r
1
(
φ
)
,
r
2
(
φ
), and r
3
(
φ
) are called generalized DRLSE such as GDRLSE1, GDRLSE2, and GDRLSE3, respectively.


### 4.4. Reaction Diffusion LSS

The RD equation is constructed by adding a diffusion term into the conventional LSE equation. Such an introduction of diffusion to LSE makes LSE stable without reinitialization. The diffusion term “ε
Δ
ϕ” was added to the LSE equation (3); we get the following equation for RD:(10)
ϕ
t
=
ε
Δ
ϕ
-
1
ε
L
ϕ
,
x
ε
Ω
⊂
R
n
Subject
⟶
ϕ
x
,
t
=
0
,
ε
=
ϕ
0
x
,where ε is a small positive constant, L
(
φ
) for PDE-based LSM or L
(
φ
)
=
-
F
δ
(
φ
) for variational LSM, Δ is the Laplacian operator defined by Δ
(
·
)
=
∑
i
=
1
n
(
∂
2
·
/
∂
x
i
2
), and ϕ
0
(
x
) is the initial LSF. Equation (10) has two dynamic processes, the diffusion term ε
Δ
ϕ gradually regularizes the LSF to be piecewise constant in each segment domain Ω
i and the reaction term “-
ε
-
1
L
(
φ
)” forces the final stable solution of (10) to L
(
φ
)
=
0, which determines Ω
i.RD Level Set Segmentation Algorithm
(1)
Initialization is as follows:ϕ
n
=
ϕ
0, n
=
0.
(2)
Compute∅
n
+
1
/
2 as ϕ
n
+
1
/
2
=
ϕ
n
=
Δ
t
2
L
(
ϕ
n
).
(3)
Compute∅
n
+
1as ϕ
n
+
1
/
2
=
ϕ
n
+
Δ
t
2
L
(
ϕ
n
).
(4)
If∅
n
+
1 satisfies stationary condition, stop; otherwise, n
=
n
+
1 and return to step (2): ϕ
n
=
ϕ
n
+
1
/
2.From the analysis in DRLSE and RD, the equilibrium solution of (10) is seen t as ε
→
0
+, which is the characteristic of phase transition. On the other hand, the said equation intrinsic problem of phase transition, that is, the stiff parameter ε
-
1, makes (10) difficult to implement. In TSSM section, we propose a splitting method to implement (10) to reduce the side effect of stiff parameter ε
-
1.


#### 4.4.1. Two-Step Splitting Method (TSSM) for RD

A TSSM algorithm to implement RD has been proposed in [21] to generate the curvature-dependent motion. In [21] the reaction function is first forced to generate a binary function with values 0 and 1, and then the diffusion function is applied to the binary function to generate curvature-dependent motion, different from [21], where the diffusion function is used to generate curvature-dependent motion, in our proposed RD based LSM; the LSE is driven by the reaction function, that is, the LSE equation. Therefore, use of the diffusion function to regularize the LSF generated by the reaction function is proposed, following TSSM to solve the RD. The steps are as follows:(1)
Solve the reaction termφ
t
=
-
ε
-
1
L
(
φ
) with φ
(
x
,
t
=
0
)
=
φ
n till some time T
r to obtain the intermediate solution, denoted by φ
n
+
1
/
2
=
φ
(
x
,
T
r
).
(2)
Solve the diffusion termφ
t
=
ε
Δ
φ, φ
(
x
,
t
=
0
)
=
φ
n
+
1
/
2 till some time T
d, and then the final level set is φ
n
+
1
=
φ
(
x
,
T
d
).Although the second step faces the rush of moving the zero level set away from its original position. But this can be eliminated through the choice of a td small enough compared to the spatial resolutions.In the above two terms, by choosing smallT
r and T
d, we can discretely approximate φ
n
+
1
/
2 and φ
n
+
1 as φ
n
+
1
/
2
=
φ
n
+
Δ
t
1
(
-
ε
-
1
L
(
φ
n
)
) and φ
n
+
1
=
φ
n
+
1
/
2
+
Δ
t
2
(
ε
Δ
φ
n
+
1
/
2
), respectively, where the time steps Δ
t
1 and Δ
t
2 represent the times T
r and T
d, respectively. Obviously, we can integrate the parameter ε into the time steps Δ
t
1 and Δ
t
2 as Δ
t
1
←
Δ
t
1
(
-
ε
-
1
) and Δ
t
2
←
Δ
t
2
ε, and hence similar to the diffusion-generated or convolution-generated curvature motion, we only need to consider the two time steps Δ
t
1 and Δ
t
2 to keep numerical stability.Speed of Segmentation. For the reinitialization methods, (4) should be iterated several times to make the LSF be an SDF while keeping the zero level set stationary. This is highly time-consuming for the reinitialization methods [21]. The GDRLSE methods are computationally much more efficient than reinitialization method. Equation (8) in each iteration of the computation of GDRLSE includes two components: the regularization term and LSE term driven by force F, in each iteration of RD method; the computation also includes two similar components. The only difference is that we split the computation into two steps: first compute the LSE term, and then compute the diffusion term. Therefore, the computation complexity of RD is similar to that of GDRLSE methods and, using Xilinx system generator implemented on hardware, its XSG are shown in the Figure 6. The proposed design is implemented using Verilog HDL, configured on Vertex-2Pro FPGA and output of FPGA is applied to monitor through VGA for displaying input image and processed image as shown Figure 7.Figure 6
Curvature central, regularization, delta, signed distance, signed distance function to circle and reinitialization of signed distance of reaction, and diffusion level set function segmentation.Figure 7
The FPGA output in the monitor through VGA. The first image shows the kidney portion segmented from the US image and second image witnesses the detected stone in that image indicated with red color.
(a)
(b)


## 4.1. Conventional LSS

In the conventional level set method, consider a closed parameterized planar curve or surface, denoted byC
(
y
,
t
)
:
[
0,1
]
x
R
+
→
R
n, where n
=
2 is for planar curve and n
=
3 is for surface, and t is the artificial time generated by the movement of the initial curve or surface C
0
(
y
) in its inward normal direction N
~. The curve or surface evolution equation is as follows [21]:(2)
C
y
,
t
=
0
=
C
0
y
,
C
t
=
F
N
¯
,where F is the force function.In the above equation intrinsic drawback of interactively solving (2) lies in its difficulty to handle topological changes of the moving front, such as splitting and merging [21]. This problem can be eliminated by using the level set method (LSM) by modifying the above equation (2) by taking the derivative with respect to time t on the both sides, yielding the following equation:(3)
ϕ
t
+
∇
ϕ
·
C
t
=
ϕ
t
+
Δ
ϕ
·
F
N
¯
=
0
,
ϕ
x
,
t
=
0
=
ϕ
0
x
,where gradient operator ∇
·
=
∂
·
/
∂
x
1
,
∂
·
/
∂
x
2
,
…
,
∂
(
·
)
/
∂
x
n and ϕ
0
(
x
) is the initial LSF C
0
(
y
)
=
x
∣
ϕ
0
(
x
)
=
0.But (3) fails for too flat or too steep near the zero level set to address this issue re-initialization is introduced.


## 4.2. Reinitialization LSS

But, during evolution, the level set function (LSF) fails for too flat or too steep near the zero level set, causing serious numerical errors. Therefore, a procedure called reinitialization is periodically employed to reshape it to be a signed distance function (SDF). In reinitialization the distance signed function isϕ
(
x
)
=
1
±
d
i
s
t
2
(
x
) where d
i
s
t
(
·
) is a distance function and ± denotes the sign inside and outside the contour [21]. But it has many problems, such as expensive computational cost, blocking the emerging of new contours [21], failures when the LSF deviates much from an SDF, and inconsistency between theory and implementation using XSG shown in Figure 5. Therefore, some formulation has been proposed to regularize the variational LSF to eliminate the reinitialization and computational cost. The following reinitialization equation is given by(4)
ϕ
t
+
S
ϕ
0
∇
ϕ
=
0
,where ϕ
0
=
ϕ
0
/
ϕ
0
2
+
Δ
x
2, ϕ
0 is the initial LSF, and Δ
x is the spatial step.Figure 5
Hardware implementation of signed distance function using XSG.


## 4.3. Distance Regularized Level Set Evolution (DRLSE)

Li et al., proposed a signed distance penalizing energy functional is given by(5)
P
ϕ
1
2
∫
Ω
∇
ϕ
-
1
2
d
x
.Equation (5) measures the closeness between an LSF (ϕ) and an SDF in the domain Ω
⊂
R
n, n
=
2 or 3. By calculus of variation [21], the gradient flow of P
(
ϕ
) is obtained as (6)
ϕ
t
=
-
P
ϕ
ϕ
=
div
⁡
r
1
ϕ
∇
ϕ
.Equation (6) is a diffusion equation with rates (7)
r
1
ϕ
=
1
-
1
∇
ϕ
.However, r
1
(
φ
)
→
-
∞ when ∇
φ
→
0, which may cause oscillation in the final LSF φ. This problem is solved by applying a new diffusion rate(8)
r
2
ϕ
=
sin
⁡
2
∏
∇
ϕ
2
∏
∇
ϕ
,
if
∇
ϕ
≤
1
1
-
1
∇
ϕ
,
if
∇
ϕ
≥
1and a constrained level set diffusion rate as (9)
r
3
ϕ
=
H
p
∇
ϕ
-
1
,where H
p
(
z
)
=
(
1
/
2
)
[
1
+
(
2
/
∏
)
a
r
c
t
a
n
(
z
/
ρ
)
].Andρ is a fixed parameter. The DRLSE methods using r
1
(
φ
)
,
r
2
(
φ
), and r
3
(
φ
) are called generalized DRLSE such as GDRLSE1, GDRLSE2, and GDRLSE3, respectively.


## 4.4. Reaction Diffusion LSS

The RD equation is constructed by adding a diffusion term into the conventional LSE equation. Such an introduction of diffusion to LSE makes LSE stable without reinitialization. The diffusion term “ε
Δ
ϕ” was added to the LSE equation (3); we get the following equation for RD:(10)
ϕ
t
=
ε
Δ
ϕ
-
1
ε
L
ϕ
,
x
ε
Ω
⊂
R
n
Subject
⟶
ϕ
x
,
t
=
0
,
ε
=
ϕ
0
x
,where ε is a small positive constant, L
(
φ
) for PDE-based LSM or L
(
φ
)
=
-
F
δ
(
φ
) for variational LSM, Δ is the Laplacian operator defined by Δ
(
·
)
=
∑
i
=
1
n
(
∂
2
·
/
∂
x
i
2
), and ϕ
0
(
x
) is the initial LSF. Equation (10) has two dynamic processes, the diffusion term ε
Δ
ϕ gradually regularizes the LSF to be piecewise constant in each segment domain Ω
i and the reaction term “-
ε
-
1
L
(
φ
)” forces the final stable solution of (10) to L
(
φ
)
=
0, which determines Ω
i.RD Level Set Segmentation Algorithm
(1)
Initialization is as follows:ϕ
n
=
ϕ
0, n
=
0.
(2)
Compute∅
n
+
1
/
2 as ϕ
n
+
1
/
2
=
ϕ
n
=
Δ
t
2
L
(
ϕ
n
).
(3)
Compute∅
n
+
1as ϕ
n
+
1
/
2
=
ϕ
n
+
Δ
t
2
L
(
ϕ
n
).
(4)
If∅
n
+
1 satisfies stationary condition, stop; otherwise, n
=
n
+
1 and return to step (2): ϕ
n
=
ϕ
n
+
1
/
2.From the analysis in DRLSE and RD, the equilibrium solution of (10) is seen t as ε
→
0
+, which is the characteristic of phase transition. On the other hand, the said equation intrinsic problem of phase transition, that is, the stiff parameter ε
-
1, makes (10) difficult to implement. In TSSM section, we propose a splitting method to implement (10) to reduce the side effect of stiff parameter ε
-
1.


### 4.4.1. Two-Step Splitting Method (TSSM) for RD

A TSSM algorithm to implement RD has been proposed in [21] to generate the curvature-dependent motion. In [21] the reaction function is first forced to generate a binary function with values 0 and 1, and then the diffusion function is applied to the binary function to generate curvature-dependent motion, different from [21], where the diffusion function is used to generate curvature-dependent motion, in our proposed RD based LSM; the LSE is driven by the reaction function, that is, the LSE equation. Therefore, use of the diffusion function to regularize the LSF generated by the reaction function is proposed, following TSSM to solve the RD. The steps are as follows:(1)
Solve the reaction termφ
t
=
-
ε
-
1
L
(
φ
) with φ
(
x
,
t
=
0
)
=
φ
n till some time T
r to obtain the intermediate solution, denoted by φ
n
+
1
/
2
=
φ
(
x
,
T
r
).
(2)
Solve the diffusion termφ
t
=
ε
Δ
φ, φ
(
x
,
t
=
0
)
=
φ
n
+
1
/
2 till some time T
d, and then the final level set is φ
n
+
1
=
φ
(
x
,
T
d
).Although the second step faces the rush of moving the zero level set away from its original position. But this can be eliminated through the choice of a td small enough compared to the spatial resolutions.In the above two terms, by choosing smallT
r and T
d, we can discretely approximate φ
n
+
1
/
2 and φ
n
+
1 as φ
n
+
1
/
2
=
φ
n
+
Δ
t
1
(
-
ε
-
1
L
(
φ
n
)
) and φ
n
+
1
=
φ
n
+
1
/
2
+
Δ
t
2
(
ε
Δ
φ
n
+
1
/
2
), respectively, where the time steps Δ
t
1 and Δ
t
2 represent the times T
r and T
d, respectively. Obviously, we can integrate the parameter ε into the time steps Δ
t
1 and Δ
t
2 as Δ
t
1
←
Δ
t
1
(
-
ε
-
1
) and Δ
t
2
←
Δ
t
2
ε, and hence similar to the diffusion-generated or convolution-generated curvature motion, we only need to consider the two time steps Δ
t
1 and Δ
t
2 to keep numerical stability.Speed of Segmentation. For the reinitialization methods, (4) should be iterated several times to make the LSF be an SDF while keeping the zero level set stationary. This is highly time-consuming for the reinitialization methods [21]. The GDRLSE methods are computationally much more efficient than reinitialization method. Equation (8) in each iteration of the computation of GDRLSE includes two components: the regularization term and LSE term driven by force F, in each iteration of RD method; the computation also includes two similar components. The only difference is that we split the computation into two steps: first compute the LSE term, and then compute the diffusion term. Therefore, the computation complexity of RD is similar to that of GDRLSE methods and, using Xilinx system generator implemented on hardware, its XSG are shown in the Figure 6. The proposed design is implemented using Verilog HDL, configured on Vertex-2Pro FPGA and output of FPGA is applied to monitor through VGA for displaying input image and processed image as shown Figure 7.Figure 6
Curvature central, regularization, delta, signed distance, signed distance function to circle and reinitialization of signed distance of reaction, and diffusion level set function segmentation.Figure 7
The FPGA output in the monitor through VGA. The first image shows the kidney portion segmented from the US image and second image witnesses the detected stone in that image indicated with red color.
(a)
(b)


## 4.4.1. Two-Step Splitting Method (TSSM) for RD

A TSSM algorithm to implement RD has been proposed in [21] to generate the curvature-dependent motion. In [21] the reaction function is first forced to generate a binary function with values 0 and 1, and then the diffusion function is applied to the binary function to generate curvature-dependent motion, different from [21], where the diffusion function is used to generate curvature-dependent motion, in our proposed RD based LSM; the LSE is driven by the reaction function, that is, the LSE equation. Therefore, use of the diffusion function to regularize the LSF generated by the reaction function is proposed, following TSSM to solve the RD. The steps are as follows:(1)
Solve the reaction termφ
t
=
-
ε
-
1
L
(
φ
) with φ
(
x
,
t
=
0
)
=
φ
n till some time T
r to obtain the intermediate solution, denoted by φ
n
+
1
/
2
=
φ
(
x
,
T
r
).
(2)
Solve the diffusion termφ
t
=
ε
Δ
φ, φ
(
x
,
t
=
0
)
=
φ
n
+
1
/
2 till some time T
d, and then the final level set is φ
n
+
1
=
φ
(
x
,
T
d
).Although the second step faces the rush of moving the zero level set away from its original position. But this can be eliminated through the choice of a td small enough compared to the spatial resolutions.In the above two terms, by choosing smallT
r and T
d, we can discretely approximate φ
n
+
1
/
2 and φ
n
+
1 as φ
n
+
1
/
2
=
φ
n
+
Δ
t
1
(
-
ε
-
1
L
(
φ
n
)
) and φ
n
+
1
=
φ
n
+
1
/
2
+
Δ
t
2
(
ε
Δ
φ
n
+
1
/
2
), respectively, where the time steps Δ
t
1 and Δ
t
2 represent the times T
r and T
d, respectively. Obviously, we can integrate the parameter ε into the time steps Δ
t
1 and Δ
t
2 as Δ
t
1
←
Δ
t
1
(
-
ε
-
1
) and Δ
t
2
←
Δ
t
2
ε, and hence similar to the diffusion-generated or convolution-generated curvature motion, we only need to consider the two time steps Δ
t
1 and Δ
t
2 to keep numerical stability.Speed of Segmentation. For the reinitialization methods, (4) should be iterated several times to make the LSF be an SDF while keeping the zero level set stationary. This is highly time-consuming for the reinitialization methods [21]. The GDRLSE methods are computationally much more efficient than reinitialization method. Equation (8) in each iteration of the computation of GDRLSE includes two components: the regularization term and LSE term driven by force F, in each iteration of RD method; the computation also includes two similar components. The only difference is that we split the computation into two steps: first compute the LSE term, and then compute the diffusion term. Therefore, the computation complexity of RD is similar to that of GDRLSE methods and, using Xilinx system generator implemented on hardware, its XSG are shown in the Figure 6. The proposed design is implemented using Verilog HDL, configured on Vertex-2Pro FPGA and output of FPGA is applied to monitor through VGA for displaying input image and processed image as shown Figure 7.Figure 6
Curvature central, regularization, delta, signed distance, signed distance function to circle and reinitialization of signed distance of reaction, and diffusion level set function segmentation.Figure 7
The FPGA output in the monitor through VGA. The first image shows the kidney portion segmented from the US image and second image witnesses the detected stone in that image indicated with red color.
(a)
(b)


## 5. Lifting Scheme Wavelets Processing

The segmented image (only stone) obtained from the previous block is applied to the lifting scheme wavelet processing block. This block consists of Daubechies filter (Db12), Symlets filter (sym12), and Biorthogonal filter (bio3.7, bio3.9, and bio4.3). InDaubechies filter (Db12) the number 12 denotes the number of vanishing moments. The higher the number of vanishing moments, the smoother the wavelet (and longer the wavelet filter) and the length of the wavelet (and scaling) filter should be twice that of the number [5].Symlets filter (sym12) extracts the kidney image features and analyses the discontinuities and abrupt changes contained in the signals. One of the 12th-order Symlets wavelets is used for feature extraction.Biorthogonal filter (bio3.7, bio3.9,andbio4.4) filter’s wavelet energy signatures were considered and averages of horizontal and vertical coefficients details were calculated. Figure 8 shows different filter used in the lifting scheme that gives different energy levels or energy features. These energy features will have significant difference if there is any stone present in the particular region or location. The identification of type of stone is explained in next section.Figure 8
Wavelets filters to extract energy features.In 2D lifting scheme wavelets transformation consists of update and predictor to get Db12, sym12, bio3.7, bio3.9, and bio4.3 as shown in Figure9.Figure 9
2D lifting scheme DWT.The equations of predict and update are given by(11)
d
i
1
=
α
x
2
i
+
x
2
i
+
2
+
x
2
i
+
1
,
a
i
1
=
β
d
i
1
+
d
i
-
1
1
+
x
2
i
,
d
i
2
=
γ
a
i
1
+
a
i
+
1
1
+
d
i
1
,
a
i
2
=
δ
d
i
2
+
d
i
-
1
2
+
a
i
1
,
d
i
=
d
i
2
ε
,where x
2
i and x
2
i
+
2 are even pixels, x
2
i
+
1 is odd pixels of stone image, and α
,
β
,
γ
,
δ
,
ε are the constants.


## 6. ANN Classification

Two architectures are used in the ANN classification, namely, multilayer perceptron and back propagation, which are described in detail in the following sections.


### 6.1. Multilayer Perceptron (MLP)

A multilayer perceptron is a feedforward artificial neural network algorithm that helps in the mapping of different sets of energy and average values obtained from the wavelets subbands energy extraction shown in Table1. These energy values are given to the input layer and multiplied with initial weights. The back propagation is the modified version of linear perceptron which uses three or more hidden layers with the nonlinear activation function. The back propagation is the most extensively used learning algorithm for multilayer perceptron in neural networks and it employs gradient descent to minimize the mean squared error between the network output value and the desired output value. These error signals are taken for completion of the weight updates which represent the power of knowledge learnt by the back propagation [22]. Multilayer perceptron with back propagation (MLP-BP) is the core algorithm. Based on the literature survey, MLP-BP algorithm was found to be better than the other algorithms in terms of accuracy, speed, and performance.Table 1
Min and max features of kidney images database which are enlarged table of the table in the GUI of lifting scheme wavelets.
Min-max
Db12Dh1 average
0.0026–0.0169
Db12cVenergy
0.0011–9.0990
e
-
04
Sym12Dh1average
0.0026–0.0169
Sym12cVenergy
0.0011–9.
0990
e
-
04
rbio3.7Dh1average
0.0052–0.0255
rbio3.7cDenergy
0.0010–8.8080
e
-
04
rbio3.7cVenergy
0.0010–9.5594
e
-
04
rbio3.7Dh1average
0.0066–0.0272
rbio3.9cDenergy
0.0010–8.0706
e
-
04
rbio3.9cVenergy
0.0014–8.9330
e
-
04
rbio3.9cVaverage
0.0039–0.0061
rbio4.4Dh1average
0.0011–8.1203
e
-
04
rbio4.4cHenergy
1.4336
e
-
04–8.7336e
-
04
rbio4.4cVenergy
0.0010–9.7450
e
-
04


## 6.1. Multilayer Perceptron (MLP)

A multilayer perceptron is a feedforward artificial neural network algorithm that helps in the mapping of different sets of energy and average values obtained from the wavelets subbands energy extraction shown in Table1. These energy values are given to the input layer and multiplied with initial weights. The back propagation is the modified version of linear perceptron which uses three or more hidden layers with the nonlinear activation function. The back propagation is the most extensively used learning algorithm for multilayer perceptron in neural networks and it employs gradient descent to minimize the mean squared error between the network output value and the desired output value. These error signals are taken for completion of the weight updates which represent the power of knowledge learnt by the back propagation [22]. Multilayer perceptron with back propagation (MLP-BP) is the core algorithm. Based on the literature survey, MLP-BP algorithm was found to be better than the other algorithms in terms of accuracy, speed, and performance.Table 1
Min and max features of kidney images database which are enlarged table of the table in the GUI of lifting scheme wavelets.
Min-max
Db12Dh1 average
0.0026–0.0169
Db12cVenergy
0.0011–9.0990
e
-
04
Sym12Dh1average
0.0026–0.0169
Sym12cVenergy
0.0011–9.
0990
e
-
04
rbio3.7Dh1average
0.0052–0.0255
rbio3.7cDenergy
0.0010–8.8080
e
-
04
rbio3.7cVenergy
0.0010–9.5594
e
-
04
rbio3.7Dh1average
0.0066–0.0272
rbio3.9cDenergy
0.0010–8.0706
e
-
04
rbio3.9cVenergy
0.0014–8.9330
e
-
04
rbio3.9cVaverage
0.0039–0.0061
rbio4.4Dh1average
0.0011–8.1203
e
-
04
rbio4.4cHenergy
1.4336
e
-
04–8.7336e
-
04
rbio4.4cVenergy
0.0010–9.7450
e
-
04


## 7. Implementation and Results

The implementation of the proposed work is completed using Verilog on Vertex-2Pro FPGA and Matlab 2012a. The Graphical User Interface (GUI) is created for the developed system as shown in Figure10. From the database of the US kidney images, one kidney image is uploaded through the GUI. The uploaded image is preprocessed and is shown in the GUI. The image segmentation option present in the GUI is executed as the next step to segment the kidney stone portion from the image. The segmented image is performed with wavelet processing by picking one of the lifting wavelet filters shown in the GUI. After selecting the filter for energy level extraction, a specific wavelet code will be invoked to get the subsequent image. Then the feature extraction option is selected to get list of energy levels extracted from the segmented image.Figure 10
Lifting scheme wavelets subband decomposition, energy feature extraction, and classifications.In the GUI shown in Figure10 there is another table which lists energy levels of all the kidney images present in the database. This is performed for testing the accuracy of MLP-BP ANN system in identifying the kidney images as normal or abnormal and to mention the type of the detected stone. Essentially the database should have both normal and abnormal images in which we should have knowledge of the number of normal and abnormal image count. During the test, it is found that our system can classify the kidney images as normal and abnormal almost with accuracy of 98.8%. The energy levels of all kidney images from database are extracted and plotted and are shown in Figure 12. The maximum peak plot indicates presence of stone. All extracted features are from range of 0 to 1 as shown in Figure 11.Figure 11
Energy and average feature values are within 0 to 1 range.Figure 12
Lifting scheme filters energy and average features of database kidney images.Table1 shows the lists of energy levels extracted from the segmented image. This table is the enlarged version of the table shown in the GUI. The rows of the table give the individual energy levels of each kidney image in the database. The columns of the table give the energy level extracted from the images in the database with respect to each wavelet filter. The first two columns correspond to Daubechies filter; the third and the fourth columns correspond to symlets12. The fifth, sixth, and seventh columns correspond to Biorthogonal filter (Bio3.7). The eighth, ninth, and tenth columns correspond to Biorthogonal filter (Bio3.9). The eleventh, twelfth, and thirteenth columns correspond to Biorthogonal filter (Bio4.4).The proposed work is implemented on Vertex-2Pro FPGA and the device speed is increased by 35% with a gate delay of 3.765 ns. The number of slices LUTs are deceased by 23% in comparison with the existing methods. The design summary table of the proposed hardware model is shown in Table2.Table 2
Design summary of proposed hardware implementation.
Parameters
Used
Available
Utilization
Number of slices Register
377
28,800
1%
Number of flip flops used
376
Number of latches used
507
28,200
1%
Mo of slice LUTs
424
28,200
1%
Dynamic power
0.488 W
Gate delay
16.5 ns


## 8. Conclusion and Scope for Future Work

The proposed work is implemented on Vertex-2Pro FPGA employing level set segmentation with momentum and resilient propagation parameters. It is found to be capable of any satisfactory achievement in identifying the stones in the US kidney image. The device performance is also pleasing with very low utilization of resources. The energy levels extracted from the lifting scheme wavelet subbands, that is, Daubechies (Db12), Symlets (sym12), and Biorthogonal filterers (bio3.7, bio3.9, and bio4.4), give the perfect indication of the difference in the energy levels of the stone portion compared to that of normal kidney region. The ANN is trained with normal kidney image and classified image input for normal or abnormal conditions by considering extracted energy levels from wavelets filters. The developed system is examined for different kidney images from the database and the results are effective in classifying the types of stone successfully with the accuracy of 98.8% [23]. Thus this system can be readily utilized in the hospitals for patients with abnormality in kidney. This work proves that the combination of level set segmentation, lifting scheme wavelet filters, and multilayer perceptron with back propagation means a better approach for the detection of stones in the kidney. In the future work, the system will be designed for real time implementation by placing biomedical sensors in the abdomen region to capture kidney portion. The captured kidney image is subjected to the proposed algorithm to process and detect stone on FPGA using hardware description language (HDL). The identified kidney stone in the image is displayed with colour for easy identification and visibility of stone in monitor.



---

*Source: 581961-2015-06-02.xml*