Ultrasound Prostate Segmentation Using Adaptive Selection Principal Curve and Smooth Mathematical Model

Abstract

Accurate prostate segmentation in ultrasound images is crucial for the clinical diagnosis of prostate cancer and for performing image-guided prostate surgery. However, it is challenging to accurately segment the prostate in ultrasound images due to their low signal-to-noise ratio, the low contrast between the prostate and neighboring tissues, and the diffuse or invisible boundaries of the prostate. In this paper, we develop a novel hybrid method for segmentation of the prostate in ultrasound images that generates accurate contours of the prostate from a range of datasets. Our method involves three key steps: (1) application of a principal curve-based method to obtain a data sequence comprising data coordinates and their corresponding projection index; (2) use of the projection index as training input for a fractional-order-based neural network that increases the accuracy of results; and (3) generation of a smooth mathematical map (expressed via the parameters of the neural network) that affords a smooth prostate boundary, which represents the output of the neural network (i.e., optimized vertices) and matches the ground truth contour. Experimental evaluation of our method and several other state-of-the-art segmentation methods on datasets of prostate ultrasound images generated at multiple institutions demonstrated that our method exhibited the best capability. Furthermore, our method is robust as it can be applied to segment prostate ultrasound images obtained at multiple institutions based on various evaluation metrics.

Appendix

• Used symbols in this work

Appendix table used symbols in this work

	Description	Symbols
Global variables	D-dimensional real number set space	RD
	Number of points in dataset	n
	X-axis coordinate of data point	x
	Y-axis coordinate of data point	y
ASPC	Principal curve	f
	Data point set/data point	Pn/pn
	Vertex/segment subset of principal curve	Vi = {v1, v2, …, viv}/Si = {s1, s2, …, sis}
	Vertex/segment of principal curve	v/s
	Number of vertices/segments of principal curve	iv/is
	Projection index	t
	Mean shift vector	m(•)
	Normalization factor	z
	Kernel bandwidth	h
	Radially symmetric kernel	K(·)
	Kernel density estimator	S(·)
	Derivative of the kernel profile	L(·)
	Bright/indeterminate/non-bright pixel point sets	TC/IC /FC
	Intensity value at data point p in the image	gd
	Gradient value at data point p in the image	Gd
	Kernel function of the indeterminacy filter	GIc
	Standard deviation of kernel function	$\sigma_{I}$
	Adjustment parameters in the linear function	ap1/ap2
	Average indeterminacy value of the current cluster point	Icavg
FBNNL	Number of neurons of input layer	I
	Number of neurons of hidden layer	H
	Number of neurons of output layer	O
	Weight from input layer to the hidden layer	w1
	Weight from hidden layer to the output layer	w2
	Threshold of the hidden neuron	a
	Threshold of the output neuron	b
	Learning rate from input layer to the hidden layer	μ1
	Learning rate from hidden layer to the output layer	μ2
	Adjustment parameter	ap
	Fractional order	α
	Iteration number	g
	Output of output units	c(•)
	Total error	E
Evaluation metrics	Dice similarity coefficient	DSC
	Jaccard similarity coefficient	OMG
	Accuracy	ACC

• Values of important neural parameters

After training, we can determine the optimal parameters of our model to determine the experimental result contour according to Eqs. (14) and (15). Table 5 represents the values of important parameters of our model for expressing the contours in three randomly selected cases, and the qualitative results of these cases are shown in Fig. 8.

Table 5 Values of neural parameters after model’s training. In this work, due to the projection index and vertices’ coordinates being the input and output of FBNNL, we set the number of input (I) and output (O) neurons to be one and two, respectively. In addition, we set the number of hidden neurons (H) be ten

• MSC method

To search the cluster of data, Cheng et al. [27] proposed the MSC method, whose details are summarized as follows:

Step 1: Compute the mean shift vector m_h(p_i), shown as follows.

To perform kernel density estimation, the kernel density estimator S(•) for p_i is denoted as

$$S(p)=\frac{z}{n\times {h}^{D}}\sum\limits_{i=1}^{n}K({\Vert \frac{p-{p}_{i}}{h}\Vert }^{2})$$

(16)

where h is kernel bandwidth with a constant, K(·) is the radially symmetric kernel, and z is a normalization factor.

The density gradient estimator is obtained as the gradient of the density estimator based on Eq. (16), shown as below:

$$\nabla S\left(p\right)=\frac{2z}{n\times h^{D+2}}\left[\sum_{i=1}^nL\left({\Arrowvert\frac{p-p_i}h\Arrowvert}^2\right)\right]\underbrace{\left[\frac{\sum_{i=1}^np_i\times L\left({\Arrowvert\frac{p-p_i}h\Arrowvert}^2\right)}{\sum_{i=1}^nL\left({\Arrowvert\frac{p-p_i}h\Arrowvert}^2\right)}-p\right]}_{mean\;shift\;vector}$$

(17)

where L(p) = K' (p) shows the derivative of the selected kernel profile. The first term of Eq. (17) is a constant, and the second term of Eq. (17) is called the mean shift vector m(p), which points to the direction of the greatest increase in density, shown below:

$${m}_{h}\left(p\right)=\frac{\sum\limits_{i=1}^{n}{p}_{i}\times L\left({\Vert \frac{{y}_{t}-{p}_{i}}{h}\Vert }^{2}\right)}{\sum\limits_{i=1}^{n}L\left({\Vert \frac{{y}_{t}-{p}_{i}}{h}\Vert }^{2}\right)}$$

(18)

Step 2: Update each initial point p_i on the m(p) direction according to p_i = + m_h(p_i).

Step 3: Iterate Step 1 and Step 2 until convergence.

Step 4: Determine the cluster point.

• NAMSC method

According to the traditional MSC method, Guo et al. [21] proposed the NAMSC method, which combined an NS-based filter into the MSC method. The details of the NAMSC method are shown as follows:

Step 1: To unify the dataset, the initial data points P_n are normalized into the range of {(− 1, − 1) ~ (1,1)}.

Step 2: Map the P_n to each channel of the neutrosophic domain, where Tc(P_n), Ic(P_n), and Fc(P_n) represent the bright, indeterminate, and non-bright pixel point sets, respectively.

$${T}_{C}\left(x,y\right)=\frac{gd\left(x,y\right)+g{d}_{\mathrm{min}}}{g{d}_{\mathrm{max}}-g{d}_{\mathrm{min}}}$$

(19)

$${I}_{C}\left(x,y\right)=\frac{Gd\left(x,y\right)+G{d}_{\mathrm{min}}}{G{d}_{\mathrm{max}}-G{d}_{\mathrm{min}}}$$

(20)

$${F}_{C}\left(x,y\right)=\frac{g{d}_{\mathrm{max}}-gd\left(x,y\right)}{g{d}_{\mathrm{max}}-g{d}_{\mathrm{min}}}$$

(21)

where gd(x, y) and Gd(x, y) are the intensity value and gradient value at the position of (x, y) in the image, respectively.

Step 3: Convolute each channel using the indeterminacy filter.

$${\sigma }_{I}\left(x,y\right)=a{p}_{1}\times {I}_{C}\left(x,y\right)+a{p}_{2}$$

(22)

$${G}_{I{\text{c}}}(u,v)=\frac{1}{2\times \pi \times {\sigma }_{I}^{2}}\mathrm{exp}\left(-\frac{{u}^{2}+{v}^{2}}{2\times {\sigma }_{I}^{2}}\right)$$

(23)

where $\sigma_{I}$ is the standard deviation, which is used to determine the shape of the kernel function. G_Ic is a kernel function of the indeterminacy filter. ap₁ and ap₂ are the adjustment parameters in the linear function, which is used to transform the indeterminate value to filter’s parameter value.

Step 4: Compute the indeterminate values of the channels in the neutrosophic domain.

$${Tc}^\prime\left(x,y\right)=\sum\limits_{v=y-\frac{fc}2}^{y+\frac{fc}2}\sum\limits_{u=x-\frac{fc}2}^{x+\frac{fc}2}Tc\left(x-u,y-v\right)\times G_{Ic}\left(u,v\right)$$

(24)

where Tc' is the result after indeterminate filter on Tc and fc is the filter size.

Step 5: Select a point p_i (x, y) and compute the bandwidth h according to the corresponding indeterminate value.

$$h(x,y)=Ic_\text{avg}(x,y)\times\left(Tc^{\prime}_{\max}-Tc^{\prime}_{\min}\right)$$

(25)

where Ic_avg is the average indeterminacy value of the current cluster point. $Tc^{\prime}_{\max}$ and $Tc^{\prime}_{\min}$ are the maximum and minimum of Tc at the current cluster among all the channels of the neutrosophic domain, respectively.

Step 6: Compute the mean shift vector m(p) (same with Step 2 of the MSC method).

Step 7: Transfer each initial point p_i on the m(p) direction (same with Step 3 of the MSC method).

Step 8: Go to Step 6 until it meets the condition $\nabla S\left(p\right)=0$ (same as Step 4 of the MSC method).

Step 9: Update the bandwidth h using the average value of the indeterminate values in the current cluster.

Step 10: Go to Step 6 until the mean of the current cluster point become unchanged.

Step 11: Go to Step 5 until all the data points are clustered.

Step 12: Exit the loop and output the cluster points.