Circle Fitting Using a Virtual Source Localization Algorithm in Wireless Sensor Networks

Abstract

A novel circle fitting algorithm is proposed in this paper. The key points of this paper are given as follows: (i) it formulates the circle fitting problem into the special source localization one in wireless sensor networks (WSN); (ii) the multidimensional scaling (MDS) analysis is applied to the data points, and thus the propagator-like method is proposed to represent the circle center parameters as the functions of the circle radius; (iii) the virtual source localization model can be rerepresented as special nonlinear equations of a unique variable (the circle radius) rather than the original three ones (the circle center and radius), and thus the classical fixed-point iteration algorithm is applied to determine the radius and the circle center parameters. The effectiveness of the proposed circle fitting approach is demonstrated using the simulation and experimental results.

1. Introduction

Circle fitting receives considerable attention because it plays an important role in computer vision, observational astronomy, structural geology, industry inspection, medical diagnosis, Iris recognition, military, security, and so forth [1–8]. For instance, to meet the increasing demand for manufacturing automation, the circle fitting technique is often applied to measure the diameter of the processing product in the manufacturing systems.

The fitting problem can be viewed as follows: estimate the parameters of a circle from a set of coplanar points. Several classical approaches [1–8], including the Hough transform (HT) methods [4, 5] and the least square (LS) approaches [6–8], have been developed to solve this problem. The former are actually to carry out a voting procedure in a three-dimensional (3D) Hough accumulator space, where every point represents a circle of a certain size. The corresponding coordinate of the local maxima is obtained as the estimated parameters of the circle. In comparison, the latter attempt to find the parameters of a circle by minimizing an error metric between the primitive and the data points.

In this paper, we develop a novel circle fitting approach by borrowing the idea from source localization in wireless sensor networks (WSN) [9, 10]. It is worthwhile to highlight the main contributions of this paper here.

• (i)

It formulates the circle fitting problem into special source localization one in WSN, where each data point should be understood as an abstract sensor node in sensor networks, and the circle center represents the localized target. However, the propagation delays are unknown, and thus the existing source localization algorithms in WSN cannot be applied to solve the special source localization problem.

(ii)

The multidimensional scaling (MDS) analysis [11] is applied to the data points, and a special covariance-like matrix is constructed. Thus, we propose the propagator-like method to represent the circle center as the functions of the circle radius.

(iii)

The virtual source localization model can be rerepresented as special nonlinear equations, where the radius is the unique variable rather than the original three ones (the circle center and radius), and thus the classical fixed-point iteration algorithm [12] is applied to determine the radius and circle center.

The rest of this paper is organized as follows. The circle fitting problem is described in Section 2. A novel circle fitting approach is developed in Section 3. Simulated and experimental results are presented in Section 4. The paper is concluded in Section 5.

2. Problem Formulation

The equation for a circle centered at $(x_{\mathrm{0}},y_{\mathrm{0}})$ with a radius r in $(x,y)$ coordinates has the following form:

$$\begin{array}{cc}& \frac{{\left(x-x_{\mathrm{0}}\right)}^{\mathrm{2}}}{r^{\mathrm{2}}}+\frac{{\left(y-y_{\mathrm{0}}\right)}^{\mathrm{2}}}{r^{\mathrm{2}}}=\mathrm{1}.\end{array}$$

(1)

The circle fitting problem (CFP) [1–8] can be described in Figure 1, that is, given data points $(x_i,y_i)$ , $i=\mathrm{1},\dots ,I$ , the objective of circle fitting is to estimate circle parameters $(r,x_0,y_0)$ that best fit to these data points. However, since in the actual application noise is introduced by some operations (e.g., the segmentation and edge detection operations in the image processing application), these obtained points are not completely precise, that is,

$$\begin{array}{cc}& \frac{{\left(x_i-x_{\mathrm{0}}\right)}^{\mathrm{2}}}{r^{\mathrm{2}}}+\frac{{\left(y_i-y_{\mathrm{0}}\right)}^{\mathrm{2}}}{r^{\mathrm{2}}}=\mathrm{1}+v_i,\end{array}$$

(2)

where $v_i$ is the introduced noise. The objective of this paper is to estimate the circle parameters $(r,x_{\mathrm{0}},y_{\mathrm{0}})$ from the given data points $(x_i,y_i)$ ,   $i=\mathrm{1},\dots ,I$ .

Figure 1

Circle fitting problem description.

3. Proposed Algorithm

In this section, we first reformulate CFP into a virtual source localization problem in wireless sensor networks (WSN) [9, 10] and then develop a novel circle fitting algorithm in this framework.

Let us review the source localization model in WSN [9, 10]:

$$\begin{array}{cc}& \frac{\mathrm{1}}{c}\parallel \mathbf{s}-{\mathit{x}}_i\parallel ={\tau }_i+w_i={\tilde{\tau }}_i,\end{array}$$

(3)

where c is the propagation speed of light; s is the unknown location of the emitter; ${\mathit{x}}_i$ is the location of the $i\text{th}$ sensor node in wireless sensor networks; ${\tau }_i$ is the (ideal) signal propagation delay from the target to the $i\text{th}$ sensor node, but ${\tilde{\tau }}_i$ is the available delay measurement, containing noise $w_i$ . In addition, $\parallel •\parallel$ denotes the Euclidean norm of vector •. The objective of source localization in WSN is to estimate s from the given measurement delays ${\tilde{\tau }}_i$ ,   $i=\mathrm{1},\dots ,I$ .

To reformulate CFP into the source localization problem in WSN, we rewrite (2) in another form as

$$\begin{array}{ll}& {\left(x_i-x_{\mathrm{0}}\right)}^{\mathrm{2}}+{\left(y_i-y_{\mathrm{0}}\right)}^{\mathrm{2}}=r^{\mathrm{2}}+r^{\mathrm{2}}{v}_i\\ \hspace{1em}\Rightarrow \parallel \mathbf{s}-{\mathbf{x}}_i\parallel =r+n_i,\hspace{1em}i=\mathrm{1},\dots ,I,\end{array}$$

(4)

where $n_i=r(\sqrt{\mathrm{1}+v_i}-\mathrm{1}),\hspace{0.17em}\mathbf{s}=[x_{\mathrm{0}}\hspace{0.17em}\hspace{0.17em}{y}_{\mathrm{0}}],\hspace{0.17em}{\mathbf{x}}_i=[x_i\hspace{0.17em}\hspace{0.17em}{y}_i]$ , and $i=\mathrm{1},\dots ,I$ .

The source localization model in (3) is quite similar to the circle fitting model in (4), especially when $c=\mathrm{1}$ : (i) each data point should be understood as a virtual “sensor node” in sensor networks; (ii) the circle center represents the virtual “emitter” or localized “target”; and (iii) the circle radius r is the virtual “propagation delay”, which are clearly described in Figure 2.

Figure 2

Virtual source localization in WSN.

By comparing (3) with (4), we can easily observe their differences; that is, all “propagation delays” from the “target” to “sensor nodes” ${\tilde{\tau }}_i$ or ${\tau }_i$ are unknown. Therefore, the existing source localization algorithms cannot solve $(x_{\mathrm{0}},y_{\mathrm{0}})$ since they require the knowledge of “propagation delays” r, ${\tilde{\tau }}_i$ , or ${\tau }_i$ .

In the rest of this section, we will develop a novel algorithm for estimating the “target position” $(x_{\mathrm{0}},y_{\mathrm{0}})$ and the “propagation delays” r of all sensor nodes.

Let

$$\begin{array}{cc}& d_{i,j}=\parallel {\mathbf{x}}_i-{\mathbf{x}}_j\parallel ,\hspace{1em}i,j=\mathrm{1},\dots ,I.\end{array}$$

(5)

And then define an $I\times I$ similarity matrix B [11]:

$$\begin{array}{cc}& \mathbf{B}=\left[\begin{array}{c}\mathbf{s}-{\mathbf{x}}_{\mathrm{1}}\\ ⋮\\ \mathbf{s}-{\mathbf{x}}_I\end{array}\right]{\left[\begin{array}{c}\mathbf{s}-{\mathbf{x}}_{\mathrm{1}}\\ ⋮\\ \mathbf{s}-{\mathbf{x}}_I\end{array}\right]}^T,\end{array}$$

(6)

the $(i,j)\text{th}$ element of which can be represented as

$$\begin{array}{ll}& \mathit{B}\left(i,j\right)=\mathrm{0.5}\left({\parallel \mathbf{s}-{\mathbf{x}}_i\parallel }^{\mathrm{2}}+{\parallel \mathbf{s}-{\mathbf{x}}_j\parallel }^{\mathrm{2}}-{\parallel {\mathbf{x}}_i-{\mathbf{x}}_j\parallel }^{\mathrm{2}}\right)\\ =\mathrm{\hspace{0.17em}\hspace{0.17em}}{r}^{\mathrm{2}}-\mathrm{0.5}{d}_{i,j}^{\mathrm{2}}+r\left(n_i+n_j\right)+\mathrm{0.5}\left(n_i^{\mathrm{2}}+n_j^{\mathrm{2}}\right).\end{array}$$

(7)

Under the ideal (without noise) case,   $\mathit{B}(i,j)=r^{\mathrm{2}}-\mathrm{0.5}{d}_{i,j}^{\mathrm{2}}$ . Note that r is unknown, and thus B is actually the function of r, that is, $\mathit{B}(r)$ .

Since the rank of $\left[\begin{array}{c}\mathbf{s}-{\mathbf{x}}_{\mathrm{1}}\\ ⋮\\ \mathbf{s}-{\mathbf{x}}_I\end{array}\right]$ equals 2, the rank of   $\mathit{B}(r)$ is also 2. From $\mathit{B}(r)$ , we introduce the following partition:

$$\begin{array}{cc}& \mathit{B}\left(r\right)=\left[{\mathit{B}}_{\mathrm{1}}\left(r\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\mathit{B}}_{\mathrm{2}}\left(r\right)\right],\end{array}$$

(8)

where ${\mathit{B}}_{\mathrm{1}}(r)$ and ${\mathit{B}}_{\mathrm{2}}(r)$ are the first two and last $(I-\mathrm{2})$ columns of $\mathit{B}(r)$ , respectively.

Similar to the conventional propagator method [13], we define the propagator

$$\begin{array}{cc}& \mathit{P}\left(r\right)={\left({\mathit{B}}_{\mathrm{1}}{\left(r\right)}^T{\mathit{B}}_{\mathrm{1}}\left(r\right)\right)}^{-\mathrm{1}}{\mathit{B}}_{\mathrm{1}}{\left(r\right)}^T{\mathit{B}}_{\mathrm{2}}\left(r\right),\end{array}$$

(9)

which satisfies

$$\begin{array}{cc}& \mathbf{P}\left(r\right)=\underset{\mathit{P}}{\min }\parallel {\mathbf{B}}_{\mathrm{2}}\left(r\right)-{\mathbf{B}}_{\mathrm{1}}\left(r\right)\mathit{P}\parallel .\end{array}$$

(10)

Let $\mathbf{x}(:,\mathrm{1})$ and $\mathbf{x}(:,\mathrm{2})$ stand for the first and second column of ${[{\mathbf{x}}_{\mathrm{1}}^T\hspace{0.17em}\hspace{0.17em}{\mathbf{x}}_{\mathrm{2}}^T\hspace{0.17em}\hspace{0.17em}\cdots \hspace{0.17em}\hspace{0.17em}{\mathbf{x}}_I^T]}^T$ . Based on the propagator method [13], we have

$$\begin{array}{cc}& {\left[x_{\mathrm{0}}\times \mathrm{1}-\mathbf{x}(:,\mathrm{1})\right]}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]=\mathrm{0},\\ {\left[y_{\mathrm{0}}\times \mathrm{1}-\mathbf{x}(:,\mathrm{2})\right]}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]=\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{0},\end{array}$$

(11)

where 1  denotes the $I\times \mathrm{1}$ vector with all elements 1.

From (11), we can solve $x_{\mathrm{0}}$ and $y_{\mathrm{0}}$ as follows:

$$\begin{array}{cc}& {\widehat{x}}_{\mathrm{0}}\left(r\right)=\frac{\left({\mathrm{1}}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]\right){\left(\mathit{x}{\left(:,\mathrm{1}\right)}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]\right)}^T}{\left({\mathrm{1}}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]\right){\left({\mathrm{1}}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]\right)}^T},\\ {\widehat{y}}_{\mathrm{0}}\left(r\right)=\frac{\left({\mathrm{1}}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]\right){\left(\mathit{x}{\left(:,\mathrm{2}\right)}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]\right)}^T}{\left({\mathrm{1}}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]\right){\left({\mathrm{1}}^T\left[\begin{array}{c}\mathbf{P}\left(r\right)\\ -\mathbf{I}\\ \end{array}\right]\right)}^T}.\end{array}$$

(12)

Note that $\mathbf{P}(r)$ depends on r and thus $x_{\mathrm{0}}(r)$ and $y_{\mathrm{0}}(r)$ are not determined from (12) directly and can only be represented as the functions of the unknown “propagation delay” r.

Plugging (12) (i.e., $x_{\mathrm{0}}(r)$ and $y_{\mathrm{0}}(r)$ ) into (4), we can obtain

$$\begin{array}{cc}& \parallel \mathbf{s}\left(r\right)-{\mathbf{x}}_i\parallel =r+n_i,\hspace{1em}i=\mathrm{1},\dots ,I,\end{array}$$

(13)

which implies that r is the root of the equations previously mentioned in absence of noise.

According to the fixed-point iteration theory [12], r is the fixed point of the function $\parallel \mathbf{s}(r)-{\mathbf{x}}_i\parallel$ , that is, the value that the function $\parallel \mathbf{s}(r)-{\mathbf{x}}_i\parallel$ “locks onto” in the iterative process. Thus, we can compute “propagation delay” r using the following iteration procedure:

$$\begin{array}{cc}& r\left(k+\mathrm{1}\right)=\parallel \mathbf{s}\left(r\left(k\right)\right)-{\mathbf{x}}_i\parallel ,\hspace{1em}k\ge \mathrm{0},\hspace{0.17em}i=\mathrm{1},\dots ,I,\end{array}$$

(14)

which are combined to yield the following iteration process:

$$\begin{array}{cc}& r\left(k+\mathrm{1}\right)=\frac{\mathrm{1}}{I}\sum _{i=\mathrm{1}}^I‍\parallel \mathit{s}\left(r\left(k\right)\right)-{\mathbf{x}}_i\parallel ,\hspace{1em}k\ge \mathrm{0}.\end{array}$$

(15)

Once circle radius r is obtained from the previously mentioned iterative procedure, the circle center $(x_{\mathrm{0}},y_{\mathrm{0}})$ can be solved from (12).

4. Simulation Results

In this section, some experiments are conducted to evaluate the performance of the proposed method. For comparison, we simultaneously implement the HT method [4, 5] and the LS approach [6–8].

4.1. Experiment 1

The first experiment is implemented on $I=\mathrm{50}$ data points of a circle shown in Figure 3, where the noise is stochastic additive white Gaussian noise with zero mean and variance 0.01, and the true circle parameters are $r=\mathrm{10},\hspace{0.17em}\hspace{0.17em}{x}_{\mathrm{0}}=\mathrm{8}$ , and $y_{\mathrm{0}}=\mathrm{15}$ . The initial value of $r(\mathrm{0})$ is 0.1. Figure 4 shows the realized $r,x_{\mathrm{0}}$ , and $y_{\mathrm{0}}$ using the proposed algorithm with 50 iterations. We can see from Figure 4 that the realization of $r,\hspace{0.17em}\hspace{0.17em}{x}_{\mathrm{0}}$ , and $y_{\mathrm{0}}$ generally stays around 10, 8, and 15, respectively, after 3 iterations in this problem, which shows that the proposed algorithm converges rapidly. Table 1 gives the estimated results using the proposed algorithm, the LS method, and the Hough transform approach. We can see from Table 1 that the proposed algorithm has the approximate estimation accuracy as those of HT and LS and is slightly more accurate in estimating $(r,x_{\mathrm{0}})$ than the LS method. Figure 5 shows the fitting results using different algorithms, which further shows that the fitting result obtained by the proposed algorithm approaches the true circle.

Table 1

Estimated results using different algorithms (Experiment  1).

Parameters	Hough	LS	Proposed
$r(10)$	10	9.9177	9.9204
$x_0(8)$	8	8.0354	8.0332
$x_0(15)$	15	14.9800	14.9703

Figure 3

Data points used in the first experiment.

Figure 4

Estimated circle parameters versus iteration number (Experiment 1).

Figure 5

Fitting results using different algorithms.

4.2. Experiment 2

In this experiment, the proposed algorithm is applied to the real data. Figure 6 (resolution $\mathrm{300}\times \mathrm{244}$ ) shows an example which computes the diameter of the bright hole after orientating the clouds of points to get the hole plane parallel to the projection plane (Available from http://www.aqsense.com/docs/docu/Compatibility.html). We can obtain the inner edge points of the hole via threshold segmentation, edge detection, and spectral clustering and then fit these points. Figure 7 shows the realized $r,\hspace{0.17em}\hspace{0.17em}{x}_{\mathrm{0}}$ , and $y_{\mathrm{0}}$ using the proposed algorithm with 50 iterations. Table 2 lists the estimation results using the proposed algorithm, the HT method, and the LS approach.

Table 2

Estimated results using different algorithms (Experiment  2).

Parameters	Reference	Hough	LS	Proposed
$\mathrm{2}\times r$ (mm)	45.07	45.0700	44.9838	44.9874
$x_{\mathrm{0}}$ (pixel)	—	161	160.5635	160.5635
$y_{\mathrm{0}}$ (pixel)	—	123	122.8443	122.8443

Figure 6

Experimental data used in Experiment 2.

Figure 7

Estimated circle parameters versus iteration number (Experiment 2).

4.3. Experiment 3

Iris recognition is a biometric identification technique based on images of the irides of an individual's eyes. Since the Iris area lies between the pupil region (a dark ellipse with the lowest intensity) and limbus region, determining the pupil region is an important preprocessing step of Iris localization. In the third experiment, we implement the proposed algorithm on the Iris image, as shown in Figure 8. Via the thresholding segmentation and Sobel edge detection, edge points are given in Figure 9. Table 3 gives the estimation results from three different methods. Figure 10 displays the fitting results by the proposed algorithm that are marked by red points, and Figure 11 shows the realized $r,\hspace{0.17em}\hspace{0.17em}{x}_{\mathrm{0}}$ , and $y_{\mathrm{0}}$ using the proposed algorithm with 50 iterations, which shows that the proposed algorithm can fit the Pupil's boundary well.

Table 3

Estimated results using different algorithms (Experiment  3).

Parameters	Hough	LS	Proposed
r (pixel)	23	23.2277	23.2237
$x_{\mathrm{0}}$ (pixel)	120	119.7545	119.7542
$y_{\mathrm{0}}$ (pixel)	156	156.4885	156.4885

Figure 8

Iris Image used in Experiment 3.

Figure 9

Edge points used in Experiment 3.

Figure 10

Fitting result using the proposed algorithm.

Figure 11

Estimated circle parameters versus iteration number (Experiment 3).

Although the HT method is of the highest estimation accuracy, it needs to be pointed out that the HT method requires (i) quantizing the three-dimensional space finely enough; otherwise the peaks in the transform plane will be broadened and (ii) the overwhelming burden of the three-dimensional search in the $r,\hspace{0.17em}\hspace{0.17em}{x}_{\mathrm{0}}$ , and $y_{\mathrm{0}}$ plane.

5. Conclusion

In this paper, we propose a novel circle fitting algorithm by borrowing the idea from source localization in wireless sensor networks. Since the virtual propagation delays of all sensor nodes are unknown, the existing source localization algorithms cannot be applied. This paper formulates the virtual source localization model of three unknown parameters $(r,x_{\mathrm{0}},y_{\mathrm{0}})$ into special nonlinear equations of a unique parameter, that is, the circle radius r, using the MDS analysis and propagator method, and then it employs the classical fixed-point iteration theory to determine the circle radius and circle center.