Maximum likelihood estimation for disk image parameters

TL;DR

This paper introduces a maximum likelihood method for estimating disk parameters from 2D images, leveraging edge pixels and intensity gradients, offering computational efficiency and clear edge distinction.

Contribution

The paper presents a new likelihood-based approach with closed-form solutions for disk parameter estimation, improving efficiency and edge differentiation over existing methods.

Findings

Effective on synthetic and real data

Requires less computational resources

Accurately distinguishes inner and outer edges

Abstract

We present a novel technique for estimating disk parameters (the centre and the radius) from its 2D image. It is based on the maximal likelihood approach utilising both edge pixels coordinates and the image intensity gradients. We emphasise the following advantages of our likelihood model. It has closed-form formulae for parameter estimating, requiring less computational resources than iterative algorithms therefore. The likelihood model naturally distinguishes the outer and inner annulus edges. The proposed technique was evaluated on both synthetic and real data.

Figures

Tables

Table 1. Table 1: The algorithm evaluation results for the disk centre point x 0 , y 0 subscript 𝑥 0 subscript 𝑦 0 x_{0},y_{0} . 25 % percent 25 25\% , 50 % percent 50 50\% , and 75 % percent 75 75\% percentiles are given for the absolute difference between the true and estimated parameters in units of pixels.

	Points	30			60			120			240			320
	Percentile	25%	50%	75%	25%	50%	75%	25%	50%	75%	25%	50%	75%	25%	50%	75%
Technique	Noise level
This paper	1	0.49	1.21	2.69	0.35	0.87	1.97	0.24	0.64	1.48	0.17	0.46	1.13	0.15	0.40	1.02
	256	0.55	1.32	2.99	0.38	0.96	2.18	0.27	0.71	1.67	0.20	0.52	1.33	0.17	0.45	1.22
	1024	0.75	1.87	4.25	0.51	1.35	3.23	0.37	1.00	2.61	0.28	0.76	2.29	0.24	0.68	2.21
Li, et al.	1	0.11	0.25	0.46	0.08	0.18	0.35	0.06	0.13	0.27	0.05	0.10	0.22	0.04	0.09	0.21
	256	0.11	0.24	0.46	0.08	0.18	0.34	0.06	0.13	0.26	0.04	0.10	0.21	0.04	0.09	0.21
	1024	0.11	0.25	0.48	0.08	0.18	0.36	0.06	0.14	0.29	0.05	0.11	0.25	0.04	0.10	0.24

Table 2. Table 2: The algorithm evaluation results for the disk radius R 𝑅 R . The table layout and units are the same as for Table  1 .

	Points	30			60			120			240			320
	Percentile	25%	50%	75%	25%	50%	75%	25%	50%	75%	25%	50%	75%	25%	50%	75%
Technique	Noise level
This paper	1	0.33	0.80	2.19	0.32	0.71	1.86	0.33	0.70	1.71	0.36	0.72	1.63	0.37	0.73	1.62
	256	0.38	0.93	2.54	0.38	0.85	2.16	0.42	0.88	2.11	0.45	0.90	2.04	0.47	0.94	2.07
	1024	0.62	1.50	4.00	0.64	1.45	3.76	0.72	1.50	3.54	0.80	1.59	3.70	0.81	1.63	3.65
Li, et al.	1	0.09	0.21	0.42	0.07	0.15	0.31	0.05	0.11	0.25	0.04	0.08	0.21	0.03	0.07	0.20
	256	0.09	0.20	0.41	0.07	0.14	0.30	0.05	0.11	0.23	0.04	0.08	0.20	0.03	0.07	0.21
	1024	0.09	0.21	0.43	0.07	0.15	0.34	0.05	0.11	0.29	0.04	0.09	0.25	0.03	0.08	0.25

Equations

$$\ln p(x_i,y_i|\theta )=-\frac{1}{2}\ln 2\pi \sigma -\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left(\sqrt{(x_i-x_0{)}^2+(y_i-y_0{)}^2}-R\right)}^2,$$

$$\ln p(x_i,y_i|\theta )=-\frac{1}{2}\ln 2\pi \sigma -\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left((x_i-x_0{)}^2+(y_i-y_0{)}^2-R^2\right)}^2,$$

$$\ln p(x_i,y_i,n_{x,i},n_{y,i}|\theta )$$

$$+{\left(n_{y,i}\left(x_0-x_i\right)-n_{x,i}\left(y_0-y_i\right)\right)}^2],$$

$$n_{x,i}$$

$$n_{y,i}$$

$$d_l$$

$$d_{\tau }$$

$$\mathrm{E}\left[d_l\right]$$

$$\mathrm{E}\left[d_{\tau }\right]$$

$$\ln p(x_i,y_i,n_{x,i},n_{y,i}|\theta )$$

$$\displaystyle\Bigl{.}-2Rn_{x,i}\left(x_{0}-x_{i}\right)-2Rn_{y,i}\left(y_{0}-y_{i}\right)\Bigr{]},$$

$$R=\frac{\frac{{\sum }_{i=1}^N{n}_{x,i}}{N}\frac{{\sum }_{i=1}^N{x}_i}{N}+\frac{{\sum }_{i=1}^N{n}_{y,i}}{N}\frac{{\sum }_{i=1}^N{y}_i}{N}-\frac{{\sum }_{i=1}^N{x}_i{n}_{x,i}}{N}-\frac{{\sum }_{i=1}^N{y}_i{n}_{y,i}}{N}}{1-{\left(\frac{{\sum }_{i=1}^N{n}_{x,i}}{N}\right)}^2-{\left(\frac{{\sum }_{i=1}^N{n}_{y,i}}{N}\right)}^2},$$

$$x_0=\frac{{\sum }_{i=1}^N{x}_i}{N}+R\frac{{\sum }_{i=1}^N{n}_{x,i}}{N},$$

$$y_0=\frac{{\sum }_{i=1}^N{y}_i}{N}+R\frac{{\sum }_{i=1}^N{n}_{y,i}}{N},$$

$${\sigma }^2$$

$$\displaystyle\Bigl{.}+\left(n_{y,i}\left(x_{0}-x_{i}\right)-n_{x,i}\left(y_{0}-y_{i}\right)\right)^{2}\Bigr{]}.$$

$$K_{x,pq}$$

$$Q\left(\theta |{\theta }^{\prime }\right)=\sum _{j=1}^M\sum _{i=1}^N{T}_{i,j}^{\prime }\ln p_j({\mathbf{x}}_i|\theta ),$$

$$T_{i,j}^{\prime }$$

$$\theta =\arg \underset{\theta }{\min }Q\left(\theta |{\theta }^{\prime }\right).$$

$$p_1(x_i,y_i,n_{x,i},n_{y,i}|\theta )$$

$$\displaystyle+\biggl{.}\left.\left(n_{y,i}\left(x_{0,1}-x_{i}\right)-n_{x,i}\left(y_{0,1}-y_{i}\right)\right)^{2}\right)\biggr{)},$$

$$p_2(x_i,y_i|\theta )=\frac{1-\tau }{2\pi {\sigma }_2^2}\exp \left(-\frac{1}{2{\sigma }_2^2}\left({\left(x_{0,2}-x_i\right)}^2+{\left(y_{0,2}-y_i\right)}^2\right)\right).$$

$$x_0\sum _{i=1}^N{n}_{y,i}^2-y_0\sum _{i=1}^N{n}_{x,i}{n}_{y,i}$$

$$x_0\sum _{i=1}^N{n}_{x,i}{n}_{y,i}-y_0\sum _{i=1}^N{n}_{x,i}^2$$

$$x_{0,1}=x_{0,2}$$

$$y_{0,1}=y_{0,2}$$

$${\sigma }_2=R$$

$${\sigma }_1$$

Full text

The proposed technique was evaluated on both synthetic and real data.

1 Introduction

Circles and disks are among the most basic geometric primitives. Therefore, the detection and fitting problems are widespread over different fields: microwave engineering [1], particle physics [2, 3], pattern recognition [4], quality control [5], robotic systems [6, 7, 8], and others. Since a disk edge is a circle, detecting a circle and detecting a disk are closely related problems. In this paper, we apply disk fitting in the experimental astronomy field.

The techniques may be coarsely divided into the two following classes. First, convolutional image-based techniques such as the circle Hough transform [9, 10], or the phase-coded annulus [11]. They work with a 2D image represented as an array of pixels. Second, point-based statistical techniques, for instance, maximal likelihood estimation. Here, input data is a list of points presumable residing on a circle. A review of different approaches is given, for instance, in [12].

In this work, we limit to fit only a single one disk or annulus per an image frame, so that the 2D image array can be straightforwardly transformed to the list of points by a gradient filtering followed by thresholding such as Otsu thresholding [13]. As soon as the maximal likelihood fitting follows the gradient filtering, we also use normalised gradient vectors during the likelihood fitting.

The paper is organised as the following. A gradient-based maximal likelihood estimation with an analytic solution is derived in Section 2. In Section 3 we describe full detection pipeline and evaluate it on synthetic data. In Section 4 we show how the expectation minimisation technique applied on top of maximal likelihood can recognise an astronomical telescope entrance pupil. The major results are discussed in the Conclusion.

2 Gradient-based likelihood model

When a table of noisy $(x,y)$-pairs is known, under a specific natural assumption we may derive the following log-likelihood loss function to estimate the centre and the radius for a circle:

[TABLE]

where $\theta$ is the set of distribution parameters: the centre $x_{0},y_{0}$, the radius $R$. The measured points are denoted by $x_{i},y_{i}$, $N$ is the total points number. For instance, this approach is used in [14, 15, 2].

However, since there is no analytic solution for minimisation of (1), one needs either to use iterative numerical methods or to change the loss function:

[TABLE]

as it is done in [1, 16].

We propose the following likelihood model:

[TABLE]

where $n_{x,i},n_{y,i}$ is the measured normalised gradient at the point $x_{i},y_{i}$. The sign before $R$ is negative for outer disk edges and positive for inner annulus edges. In essence, the gradient $n_{x,i},n_{y,i}$ is the normal vector at the circle point $x_{i},y_{i}$.

Let us show how equation (3) is derived. We have the following four quantities for each point: $x_{i},y_{i},n_{x,i},n_{y,i}$. Since the points are supposed to be located at the circle with some precision, the quantities are not independent random variables. It is known that the vectors normal to the circle are crossed at its centre, so we assume that the following relations are held:

[TABLE]

where $\epsilon_{x,i}$ and $\epsilon_{y,i}$ are random variables with zero mean. Consider the following two new random variables which geometrical meaning is explained in Fig. 1:

[TABLE]

Using equations (4) and (5), one may find mean values for $d_{l}$ and $d_{\tau}$:

[TABLE]

where ${\mathrm{E}}\left[\cdot\right]$ denotes an expectation operator. Here we assume that $\epsilon_{x,i}$ and $\epsilon_{y,i}$ are independent from $x_{i}$ and $y_{i}$.

There can also be an alternative geometric interpretation for (3). It is easy to show that every circle is located at some two-dimensional plane in the four-dimensional space of $x,y,n_{x},n_{y}$. Equation (3) suggests minimal distances between the plane and the measurement points in this four-dimensional space.

The model has two important properties. First, it can be reduced to the exact model (1). Indeed, expanding and collecting the terms in (3) we have the following:

[TABLE]

which is the same as (1), when equations (4) and (5) are substituted in (10), and $\epsilon_{x,i}=\epsilon_{y,i}=0$. That corresponds to the case of unobservable normal vectors.

Second, it has the following simple analytic solution for the maximal likelihood estimators:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In [17], a similar approach is used, authors implicitly constructed normal vectors from the model parameters and the measured point positions. However, since they didn’t measure the normal vectors, their disk parameters estimates are obtained through an iterative procedure.

3 Maximal likelihood disk fitting

Let us recall that we are interested in locating a single disk on a 2D image. We suppose that the image is represented as a 2D grey scale array. The following two steps are performed on the data to convert the array to the list of circle points coordinates suitable for applying equation (3).

First, the gradient filters are applied to produce the horizontal gradient $g_{x,nm}$ and vertical gradient $g_{y,nm}$ component maps. The following convolution kernels are used to evaluate the gradients:

[TABLE]

where $s^{2}$ is selected to be $2$ to fit the kernel into five pixels. We also ignore the Gaussian norm in the kernel, since we are interested in normalised gradient vectors rather than absolute gradient values. Kernels (15) and (16) are separable111i.e. a product of two expressions each depending on its own index., so the filtering can be efficiently applied as a sequence of two 1D convolutions. When the both gradient projections are known, the gradient norm $g_{nm}=\sqrt{g_{x,nm}^{2}+g_{y,nm}^{2}}$ map can be evaluated.

Second, Otsu threshold selection method [13] is applied to the gradient norm map to detect pixels belonging to the circle (disk edge). The thresholding technique relies on the data distribution itself, so the Gaussian norm in kernels (15)–(16) may be safely ignored here too. All the image pixels are divided into the circle and the background using the threshold. There may be too many first class pixels, so the smaller fixed-size subset is chosen at random. Adjusting subset size helps to find a balance between computational expenses and accuracy. Finally, for every pixel from the subset we know its integer coordinates $x_{i}=n,y_{i}=m$ and normalised gradient vector $n_{x,i}=g_{x,nm}\cdot g_{nm}^{-1},n_{y,i}=g_{y,nm}\cdot g_{nm}^{-1}$. Then, the gradient-based maximum likelihood algorithm is applied to the list to estimate circle parameters.

To evaluate the algorithm performance, we generate $10000$ sample disk images $640\times 480$ pixels each. Each disk has a random size and a random position. The radius is distributed uniformly between $30$ and $270$ pixels. The position is distributed uniformly between zero and the maximal dimension to allow existing of partial disk segments. The disk pixels have a value of $255$, the background is [math], and then we apply Poisson noise of different scales $\lambda$ to the pixels of this input data set. A sample is shown at Fig. 2. For each combination of the noise level $\lambda$ and the points number each input image is processed as described above, the circle parameters are estimated using equations (11)–(14), and an a-posteriori estimation error is evaluated. The maximal absolute parameter difference is used as an accuracy measure. The disk parameters are expressed in units of pixels.

The proposed algorithm is also compared to the algorithm of Li et al. (referred as Li’s algorithm below) presented in [17]. Since that algorithm is iterative, the termination criterion is chosen to be so that relative error is less than $10^{-4}$ to ensure that the accuracy is limited by the data properties rather than an insufficient number of iterations.

The evaluation results are given in Table 1 and Table 2. We may see that the proposed algorithm has an important drawback. It has a considerably lower accuracy than the algorithm from [17] has. The origin of this additional error is clear. When we add normalised vectors $(n_{x,i},n_{y,i})$ into consideration, we also introduce additional measurements error. On the other hand, it takes only about $60\times 10^{-6}s$ for the proposed algorithm to process one circle, while it takes $2\times 10^{-3}s$ for Li’s algorithm to reach the obtained accuracy. The execution times were measured at the Intel i5-7200U CPU for Python based implementation. They are consistent within this paper but should be thought of as relative values.

The execution times can be made comparable by limiting the maximal iteration number for Li’s algorithm. The median error for Li’s algorithm is about $20$ pixels in that case. It is natural to try using the results of the proposed algorithm as initial estimates for iterative Li’s algorithm. In that case, Li’s algorithm execution time is reduced to $0.9\times 10^{-3}s$ keeping the same good accuracy as in Tables 1 and 2.

The evaluation results make us optimistic about the algorithm applications. Since the proposed algorithm is not an iterative one, its major advance is the low and predictable execution time.

4 Astronomical telescope pupil detection

Simplicity of equations (11)–(14) providing optimal parameters for likelihood function (3) allows us to construct more complex models on top of it. The expectation maximisation (EM) algorithm [18] is usually considered as an elegant and effective way to maximise complex likelihood functions for so-called mixture models, when the probability density function for every observation is expressed as a linear combination of the single-class probability density functions each having its own parameters to be estimated.

When $N$ independent observations are given, the EM algorithm iteratively optimises the following goal function:

[TABLE]

where $\theta$ is a distribution parameters vector to be estimated, ${\mathbf{x}_{i}}$ is a vector of $i$-th observation, index $j$ enumerates the class, and weights $T^{\prime}_{i,j}$ are evaluated based on the previously estimated distribution parameters $\theta^{\prime}$ as:

[TABLE]

The next step parameters estimate $\theta$ is evaluated as a solution for the following optimisation problem:

[TABLE]

It is important to note that the equations for $\theta$ can be easily derived if a solution for every $j$-th single-class likelihood problem is known. This property follows from the linearity of (17) and forms the key part of EM technique success. If $p_{j}$ is described by (3) for some $j$, then the corresponding parameters $x_{0},y_{0},R,\sigma$ are evaluated using modified equations (11)–(14) where $\sum^{N}_{i=1}T^{\prime}_{i,j}\left(\cdot\right)_{i}$ is substituted for each $\sum^{N}_{i=1}\left(\cdot\right)_{i}$, and $\sum^{N}_{i=1}T^{\prime}_{i,j}$ is substituted for $N$.

Let us consider the following real-data case as an example for applying the EM algorithm on top of the proposed disk parameters estimation algorithm.

An astronomical telescope focuses parallel beams of the light coming from stars at a focal plane. Usually, two-mirror optical schemes are common nowadays. A typical optical scheme for a modern optical astronomical telescope is given at Fig. 3. Note, that the part of the parallel light beam is dissipated by the secondary mirror back surface. The outer size of the beam is limited by the primary mirror size. An optical scheme part limiting the amount of light passing through the whole system is called an actual aperture stop. We may imagine that the aperture stop splits the whole optical system into the two parts. An image of the actual aperture stop produced by the back of the optical system is called an exit pupil. An entrance pupil is the image of the actual aperture stop produced by the front optical system part in the reverse direction. The entrance pupil is identical to the actual aperture stop in case of typical astronomical telescope optical scheme given at Fig. 3. The pupil image can be obtained by putting an image detector at the exit pupil. The sample for the real telescope is given at Fig. 4.

Installing additional equipment in the exit pupil is used, for instance, while examining the telescope optics quality [19, 20] or optical turbulence profiling [21]. Since the telescope optical system geometry is always a subject of thermal expansions and gravitational bending, the exact position for the exit pupil is known with limited accuracy. In real conditions, an additional adjustment would be required. The sharpness of the outer pupil edge may be controlled to achieve the goal. Note, that the inner pupil edge can be unsharp at the same moment since it is produced by the secondary mirror edge which is located at some distance from the primary mirror. This is why the outer edge locating, followed by spatial filtering, is required before sharpness detection.

The proposed algorithm naturally distinguishes between an annulus outer edge and inner one. The outer edge has gradients directed towards the circle centre where they are crossed. At the same time, the inner edge has gradients directed in the reverse direction. We could craft an algorithm detecting inner edges by changing the sign for $R$ in equation (3), assuming that $R$ is always positive.

Two classes of edge points are obtained by the edge detection algorithm in the case of images similar to Fig. 4. The first class is consisted from the outer edge points. The second class consists of the inner edge points and the spiders edge points. We assume that the gradient vectors direction is distributed uniformly for the latter class. Since the class is a hidden variable, it seems to be natural to apply mentioned expectation maximisation algorithm [18] on the following two-component model. The probability density functions are defined as follows:

[TABLE]

[TABLE]

Here, $p_{1}$ follows from equation (3) and represents the outer edge point distribution. The vector ${\theta}$ stands for the probability density function parameters to be determined: $\tau,x_{0,1},y_{0,1},R,\sigma_{1},x_{0,2},y_{0,2},\sigma_{2}$. The parameters are evaluated using modified equations (11)–(14) and ordinarily equations for the Gaussian mixture EM model.

The initial parameter estimation is required to start the iteration process. The following has been found to be working satisfactory. The pupil centre is estimated as the solution for the following linear equation set:

[TABLE]

while other initial parameters are evaluated as follows:

[TABLE]

We examined two data sets of real data: one for the focused pupil image and another for the defocused pupil. Each data set consists of $1000$ individual images obtained in a series under the same conditions. The pupil parameters estimator accuracy can be evaluated by running the algorithm for every individual image frame followed by evaluating the parameters standard deviation.

In the focused pupil case, we have found that the standard deviation for the pupil centre $\sigma_{c}=0.85$ pix. and the standard deviation for the outer pupil radius $\sigma_{R}=0.42$ pix. Given the average outer radius $R=113.5$ pix, we have a relative accuracy of $0.4\%$ here. For the defocused pupil case, $\sigma_{c}=1.2$ pix. and $\sigma_{R}=0.82$ pix, while the average outer radius remains almost the same $R=113.4$ pix.

It takes approximately $120\times 10^{-3}s$ to process a single $659\times 493$ frame using Python programming language in both cases. Most of the time is spent at the gradient filtering stage. However, our production C++ based algorithm implementation requires only $6\times 10^{-3}s$ that allows us to process approximately $160$ frames per second.

Hence, we have found that the proposed algorithm is able to cope with real data.

5 Conclusion

In this paper, the new likelihood disk fitting model was presented. We showed how using additional information about image gradients can change the model properties. We obtained likelihood model (3) that has simple equations (11)–(14) for its solution and can be reduced to the exact likelihood model (1) at the same time. The likelihood model parameters evaluating doesn’t require iterative techniques because of the equations’ simplicity. So, compared to iterative techniques, less computational resources are generally required, and the computational time is more predictable, which may be important in real-time image processing.

At the same time, the proposed model also has the following considerable drawbacks. First, additional input data are required. However, when gradient-based edge detection algorithms are applied during the image processing pipeline, the required data is usually already available. Second, it can only apply to disks or annuluses. Fitting real circles may require further algorithm complication, such as considering a circle as a thin annulus and using expectation maximisation technique to fit the annulus parameters using both its outer and inner edges. Third, it has a lower accuracy than can be reached by other algorithms under the same conditions.

In Sections 3 and 4, we demonstrated that the proposed likelihood model may be efficiently used both solely and for estimating initial parameters for iterative circle fitting algorithms.