Image Processing and Concentric Ellipse Fitting to Estimate the Ellipticity of Steel Coils

1. Introduction

In this chapter a particular image processing problem existing in quality assurance for steel industry is addressed. Since the quality standards in steel industry are ever increasing, even steel industry starts to use machine vision systems to monitor their product quality. In case of steel coils a particular figure of merit in quality is the roundness (ellipticity) or lack thereof of the coil. The mathematical correct measure of the roundness is the eccentricity$e=\frac{\sqrt{a^2-b^2}}{a}$, with the major semi axis a and the minor semi axis b. But the fastest way for the quality inspector is just to measure a and b of the coil and calculate their ratio$e=\frac{a}{b}$. If ( exceeds a set limit (≈1.03 measured at the innermost coil layer), the coil is considered as cull. As discussed by Schleicher & Zagar (2008; 2009a) the extraction of elliptic arc segments for the ellipticity estimation and the basic concentric ellipse fitting algorithm are the main parts of the algorithm estimating the coils ellipticity. Comparing to ordinary ellipse fitting tasks that require input data distributed over the whole circumference, the estimation of the ellipse parameters using only data forming low length arcs is a much more challenging problem. To achieve a precise estimation of the ellipticity the fact of the common center of the elliptic arcs is used by the concentric ellipse fitting algorithm to estimate the center of the coil precisely. The functionality of the algorithm is show based upon of real coil images.

Indeed some error sources reduce the accuracy, which are discussed at the end of the chapter. The chapter is organized as follows: First the image acquisition system is described. Then the algorithm is explained consecutively in Sec. 3. Further sources of different errors are investigated in Sec. 4. Finally major results are presented (Sec. 5) and a conclusion (Sec. 6) is given.

2. Image acquisition system

The image acquisition setup shown in Fig. 1a consists of a camera placed on a tripod, a laser trigger unit to initiate the acquisition and a laptop to store, process and view the images. A single image frame $$I_{raw}\left(x,y\right)$$of the coil is acquired, when the coil transported by a lifting ramp reaches its end position triggering a light barrier. The camera is a 1/2” monochrome firewire camera (AVT Guppy 146–B) with 1392×1040 pixels. Its optical zoom is adjusted so as to image the biggest coil full frame, which has an outer diameter of 2 m. The sheet metal layers are 0.5 to 2.5 mm thick and the inner diameter of each coil is known to be 0.6 m. To ensure a circular image of a truly circular coil, the camera chip plane must be oriented parallel to the coils front side (Scheimpflug condition (Luhmann, 2003, p.147)). Figure 2 shows the camera coordinate system $\left\{{\overrightarrow{e}}_x,{\overrightarrow{e}}_y,{\overrightarrow{e}}_c\right\}$ in relation to the world coordinate system$\left\{{\overrightarrow{e}}_{\xi },{\overrightarrow{e}}_{\eta },{\overrightarrow{e}}_{\zeta }\right\}$, where the$$\xi ,\eta$$–plane (coil front side) is kept parallel to the${\overrightarrow{e}}_x,{\overrightarrow{e}}_y$–plane (camera chip). To assure the correct alignment of the camera, a spirit level was used to do the horizontal alignment (${\overrightarrow{e}}_x$,${\overrightarrow{e}}_z$–plane). Furthermore the horizontal edges of the coils carrier were used to align the rotation about the${\overrightarrow{e}}_y$–axis that the coils front side is parallel to the camera chip plane. With this setup the image $$I_{raw}\left(x,y\right)$$is acquired and processed by the algorithm in the next section.

Figure 1.

Image acquisition setup (a) and typical captured image (b).

Figure 2.

The cameras${\overrightarrow{e}}_x$,${\overrightarrow{e}}_y$–plane is oriented parallel to the coils front side ${\overrightarrow{e}}_{\xi \text{'}}{\overrightarrow{e}}_{\eta }$ – plane.

3. Algorithm

The block diagram of the algorithm, shown in Fig. 3, reflects the structure of this section. The preprocessing, circular hough transform and extraction of all curve segments $${\text{C}}_a$$are applied on the pixel level and then a reduced dataset of selected curve segments $${\text{C}}_s$$is used for the ellipticity estimation.

Figure 3.

Block diagram of the algorithm to estimate ellipticities. The preprocessing is done at pixel level prior to the ellipticity estimation done at curves level.

3.1. Preprocessing

The first preprocessing step is to crop $$I_{raw}(x,y)$$to the region of interest, which fits the biggest coil to capture. Further the unsharp masking sharpening algorithm Eqn 1 (Jain, 1989, pp.249– 251) is used to boost the edges of the coil.

$$I\left(x,y\right)=Iraw\left(x,y\right)+\lambda G\left(x,y\right)$$E1

where G(x,y) is the discrete Laplacian

$$G\left(x,y\right)=Iraw\left(x,y\right)-\frac{1}{4}\left[Iraw\left(x-1,y\right)+Iraw\left(x,y-1\right)+Iraw\left(x+1,y\right)+Iraw\left(x,y+1\right)\right]$$E2

of the raw image. The parameter λ determines the sharpening power of the algorithm. Setting λ too large leads to a very noisy image result, since the Laplacian is very sensitive to noise. Since the sharpening process is similar to a high pass filter, any influence of the inhomogeneous natural lighting is reduced.

The next important preprocessing step is the segmentation of the image into background and foreground pixels (front side highlighted in Fig. 1a). Usually this is done by comparing the actual image to a background image only (Zivkovic, 2004). Since only a single image of the coil is available this procedure can not be applied and the results of the circular hough transform are used to first find the coils center and then to segment the coils front side, which is described in the next section.

3.2. Circular Hough Transform (CHT)

The circular Hough transform (Illingworth & Kittler, 1987; Peng, 2007; Peng et al., 2007) is a transform mapping points in the x,y–plane belonging to circular slopes into parameter space. Thus filling a so called accumulator. A circle in x,y–space can be described by

$${\left(x-x_c\right)}^2+{\left(y-y_c\right)}^2=r^2$$E3

In parameter space it is then$$\left(x_c,y_c,r\right)$$. This 3D parameter space can be reduced to a 2D with$$\left(x_c,y_c\right)$$, we incorporate the constraint that the vectors, which are normal to the circle boundary, must all intersect at the circle center$$\left(x_c,y_c\right)$$. Figure 4 illustrates this relation between the circle and its center. These normal directions can be obtained by local gray–level detection operators. The normal direction can also be expressed by

$$\tan \theta =\frac{y-y_c}{x-x_c}$$E4

Figure 4.

Circle parameters for the circular Hough transform. The points on the circle (x,y) in image space is mapped onto a line in the accumulator.

Rearranging Eqn 4 yields

$$y_c=x_c\tan \theta +\left(y-x\tan \theta \right)$$E5

and transfers the problem in a form similar to the well known (line) Hough transform (HT). The accumulator of the HT, where the circles centers $\left({\widehat{x}}_c,{\widehat{y}}_c\right)$ is determined by the peak, is shown in Fig. 5a. Since the determined center is only a rough estimate (denoted by ˆ), the center$$\left(x_c,y_c\right)$$. After incorporating the circle center $\left({\widehat{x}}_c,{\widehat{y}}_c\right)$ in Eqn 3, the radius can be found by histogramming,$r^2={\left(x-{\widehat{x}}_c\right)}^2+{\left(y-{\widehat{y}}_c\right)}^2$. This histogram, also called signature curve, shown in Fig. 5b, is used to obtain an estimate for the inner radius r_i aswell as the outer:

Figure 5.

Accumulator and signature curve of the CHT of sample image of Fig. 1b.

radius r_o (marked by a red asterisk). Local maxima in the signature curve are marked with a blue asterisk and signals the presence of a circular shape detected in the front side image.

To reduce the execution time of the CHT algorithm the preprocessed image $$I\left(x,y\right)$$is resized (down sampled) to 150×170 pixels.

Due to the reduced image size the center estimation is only very rough, though fast. Furthermore the two stage concept causes the error of the center estimate to propagate further to the radius estimation in the signature curve. That’s why the signature curve is not very reliable to return the proper real inner and outer radius, respectively. When no proper radii are found, the maximum coil size is chosen to segment the coil front side, since subsequent selection criteria will be applied to the curves before they are used for the fitting process.

3.3. Curve extraction

To extract arc curve segments at least of some of the sheet metal layers the Canny edge detector (Canny, 1986) is applied to$$I\left(x,y\right)$$. As a result this edge detector creates a binarized image $$I_t\left(x,y\right)$$with thinned out edges. Due to the hysteresis threshold technique (Nixon & Aguado, 2002) implemented within the Canny edge detector, even weak edges are reliably extracted. As a drawback however, some edges, especially those resulting from rather weak gradients, have ramifications and squiggles that need to be taken care of.

Figure 6.

Different 3 ×3 erosion masks M₁... (a), M_2... (b), M₃... (c), are used to erode all possible types of ramifications of a thinned binary image.

First the ramifications are removed by erosion (Gonzales & Woods, 2002, pp. 525–527) with the templates shown in Fig. 6 and their permutations by 0 , 90 , 180 and 270 .

$$I_c\left(x,y\right)=I_t\left(x,y\right)-{\displaystyle \sum _{M=\left\{M_1,M_2,M_3\right\}}{\displaystyle \sum _{\varphi =\left\{0^0,{90}^0,{180}^0,{270}^0\right\}}{I}_t}}\left(x,y\right)\Theta rot\left(M,\varphi \right)$$E6

Figure 7.

a) is a cutout of the thinned binary edge image$$I_t(x,y)$$, which is eroded by the templates shown in Fig. 6. (b) shows only points which are removed. (c) is the cutout of the resultant image $$I_c(x,y)$$without ramifications (a)–(b)

After the erosion, the resultant image $$I_c\left(x,y\right)$$contains only disrupted curve segments, which are remapped from the image matrix onto a curve set$${\text{C}}_a=c_{\text{1}},...,c_j,...,c_J$$, which contains all extracted curve segments denoted by $$c_j$$(See Fig. 7). The curve segment c_j is a $$M\times \text{2}$$vector containing the curve pixels coordinates.

3.4. Curve selection

Since the extracted curve set $${\text{C}}_a$$also contains curves derived from the edges of the coils tapes or packaging labels, which must not be processed in the ellipticity estimation algorithm, these curves must also be remove to increase the quality of the estimators input data. Therefore figures of merits are calculated to select only those that are kept for further processing $$({\text{C}}_s)$$and reject straight lines, rather short and very squiggled curves. Useful figures of merit are the curve length and the extend, which is the ratio between curve bounding box and curve length. Furthermore closed curve segments, with an eccentricity $$e$$of the ellipse that has the same second-moments as the region ($$e=0$$is actually a circle;$$e=\text{1}$$ is an ellipse degenerated to a line segment) over a specific threshold (in this application 0.1 turned out to be useful), are also rejected (Rafael C. Gonzalez, 2003).

Figure 8.

The tangential vector of a circular arc is normal to the radius vector. (1) acts as a seed to an elliptic curve, (2) will be removed, (3) will be part of an elliptic curve, if it is connected to a seed point.

To select only elliptic arc segments which are part of the coil the scalar product of$$n_r=r/\left|r\right|$$, the normalized radius vector, and t, the tangential vector, is calculated for each point on each arc curve segment c_j. Figure 8 shows two curves with typical cases indicated. Concerning an ideal round coil, the scalar product $$n_r\cdot t=0$$for each arc segment. Since n_r and t are normalized, the range of the result of the scalar product is [0,...,1]. The subsequently applied hysteresis thresholding technique (Nixon & Aguado, 2002) returns only reasonable curve segments for the final ellipticity estimation algorithm.

3.5. Ellipticity estimation

To estimate the ellipticity an ellipse fit is applied to the dataset. Curve fitting is an extensive topic and even ellipse fitting is discussed in literature very detailed (Gander et al., 1994; Matei & Meer, 2000b; Leedan & Meer, 2000; Halir, 1998; Ahn et al., 2001; Nair & Saunders, 1996; Liu et al., 2006; Lei &Wong, 1999; Fitzgibbon et al., 1999; Halir, 1998).

The paper written by Gander et al. (1994), was one of the first that covered circle and ellipse fitting using least squares and iterative methods.

Figure 10 shows typical fitting results when applying the general ellipse fit algorithm to the previously extracted arc segment dataset$${\text{C}}_s$$. Obviously, the common center of the coil can’t be reliably determined, since the general ellipse fit algorithm is executed only for one arc segment at a time. This renders it necessary to consider additional constraints for the fitting algorithm and adopt it to use the whole dataset with all arc segments simultaneously.

In the following sections first the chosen curve fitting algorithm is described on a general way. Then the general, the simple and the new concentric ellipse fit, which incorporates the fact of the common coil center, is described in more detail.

Figure 9.

Plot of the selected elliptic curve segments $${\text{C}}_s$$used for the subsequent ellipticity estimation.

Figure 10.

Ellipse fitting of curve segments without any constraints leads to unusable results for the ellipticity estimation problem of a coil.

3.5.1. Curve fitting

Before doing the ellipse fit, a short introduction into curve fitting will be given. Therefore the same notation as proposed by Chernov (2007) will be used.

The principal equation for algebraic curve fitting is

$$u{\left(x\right)}^T\Theta =0$$E7

with$$u\left(x\right)={\left[u_{\text{1}}\left(x\right),...,u_L\left(x\right)\right]}^T$$, which is the data transformed in a problem-dependent manner and$$Q={[q_{\text{1}},...,q_L]}^T$$, which contains the L unknown variables. The measured data vector $$x={\left[x_{\text{1}},...,x_I\right]}^T{}_{}$$is an I element long vector, assumed to represent the true data ${\overline{x}}_i$ corrupted by some random noise $$\delta x_i$$(here assumed to be normally distributed with$$N(0,{\sigma }^{\text{2}}I))$$.

$$x_i={\overline{x}}_i+\delta x_i,i=1,\mathrm{...},I$$E8

To estimate Θ the approximate maximum likelihood (AML) method (Chernov, 2007) is used.

$$F_{AML}={\displaystyle \sum _{i=1}^I\frac{{\Theta }^{T}u\left(x_i\right)u{\left(x_i\right)}^T\Theta }{{\Theta }^T\left[{\nabla }_{x}u\left(x_i\right)\right]{\left[{\nabla }_{x}u\left(x_i\right)\right]}^T\Theta }}=\frac{{\Theta }^T{A}_i\Theta }{{\Theta }^T{B}_i\Theta }\to \min$$E9

The square matrices $$A_i$$and $$B_i$$depend only on x_i and have the size$$L\times L$$.

$$A_i=u\left(x_i\right)u{\left(x_i\right)}^T$$E10

$$B_i=\left[{\nabla }_{x}u\left(x_i\right)\right]{\left[{\nabla }_{x}u\left(x_i\right)\right]}^T$$E11

Using these matrices for the representation of the input data many popular schemes can be used to solve the minimization problem:

1. (TAU) Taubin’s fit (Taubin, 1991)

2. (HEIV) Heteroscedastic errors–in–variables (Leedan & Meer, 2000; Matei & Meer, 2000a)

3. (FNS) Fundamental numerical scheme (Chojnacki et al., 2001; 2005; 2004)

4. (RED) Reduced scheme (Kanatani, 2005)

5. (REN) Renormalization procedure (Kanatani, 2005; 1997; 2008)

All of the above algorithms were tested for the concentric ellipse fitting problem and only REN showed unstable behavior. The other algorithms only vary in execution time and achieved accuracy, which is also dependent on the maximum number of iterations and the set tolerances. For additional information about these schemes we refer to Chernov (2007). In the following sections the way to obtain the $$A_i$$and $$B_i$$matrices, $${Ñ}_{x}u\left(x_{\text{i}}\right)$$and$$u\left(x_{\text{i}}\right)$$, respectively, will be described.

3.5.2. General ellipse fitting

A general ellipse can be described by the equation for the algebraic distance

$$F_{general}\left(x,y\right)=A{x}^2+Byx+C{y}^2+Dx+Ey+F=0$$E12

To obtain the unknown parameters A,..., F of Eqn 13 the AML method is used. Therefore the equation is written in terms u(x) and Θ.

$$x={\left[x,y\right]}^T$$E13

$$u\left(x\right)={\left[x^2,xy,y^2,x,y,1\right]}^T$$E14

$$\Theta ={\left[A,B,C,D,E,F\right]}^T$$E15

Further we obtain

$$A_i=u\left(x_i\right)u{\left(x_i\right)}^T=\left(\begin{array}{l}{x}_i^4,x_i^3{y}_i,x_i^2{y}_i^2,x_i^3,x_i^2{y}_i,x_i^2\\ x_i^3{y}_i,x_i^2{y}_i^2,x_i{y}_i^3,x_i^2{y}_i,x_i{y}_i\\ x_i^2{y}_i^2,x_i{y}_i^3,y^4,x_i{y}_i^2,y_i^3,y_i^2\\ x_i^3,x_i^2{y}_i,x_i{y}_i^2,x_i^2,x_i{y}_i,x_i^{}\\ x_i^2{y}_i,x_i{y}_i^2,y_i^3,x_i{y}_i,y_i^2,y_i\\ x_i^2,x_i{y}_i,y_i^2,x_i^{},y_i,1\end{array}\right)$$E16

And

$$B_i=\left[{\nabla }_{x}u\left(x_i\right)\right]{\left[{\nabla }_{x}u\left(x_i\right)\right]}^T=\left(\begin{array}{l}4x_i^2,2x_i{y}_i,0,2x_i,0,0\\ 2x_i{y}_i,x_i^2+y_i^2,2x_i{y}_i,1,0,0\\ 0,2x_i{y}_i,4y_i^2,0,2y_i,0\\ 2x_i,1,0,1,0,0\\ 0,0,2y_i,0,1,0\\ 0,0,0,0,0,0\end{array}\right)$$E17

With

$${\nabla }_{x}u\left(x\right)={\left[\begin{array}{l}2x,y,0,1,0,0\\ 0,x,2y,0,1,0\end{array}\right]}^T$$E18

3.5.3. Simple ellipse fitting

The denotation simple ellipse fitting (SEF) is used when the simplest form of an ellipse, with the center at the origin (0,0) and no rotation, is to be fitted. Equation 12 is then simplified to

$$F_{simple}\left(u,v\right)=A{u}^2+C{v}^2+F=0=\frac{v^2}{a^2}+\frac{u^2}{b^2}=0$$E19

Here a and b are the major and minor semi axes and u and v are the coordinates related to the semi axis of the ellipse. To obtain u and v the ellipse must be centered in the origin thus:

$$v=x-x_c$$E20

$$u=y-y_c$$E21

For $$A_i$$and $$B_i$$calculated by Eqn 12 we further need

$$x={\left[u,v\right]}^T$$E22

$$u\left(x\right)={\left[u^2,v^2,1\right]}^T$$E23

$$\Theta ={\left[A,C,F\right]}^T$$E24

and

$${\nabla }_{x}u\left(x\right)=\left[\begin{array}{l}2u,0\\ 0,2v\end{array}\right]$$E25

to calculate $$A_i$$and$$B_i$$.

3.5.4. Concentric ellipses fitting

The concentric ellipse fit takes the whole dataset $${\text{C}}_s$$of the coils elliptic arc segments into account, in contrast to the previously discussed algorithms, which use only one curve segment at a time. Thereby a highly precise center estimation is possible. To achieve this, some simplifications needs to be made:

1. The rotation of the ellipses is zero.

2. The ellipticity of all elliptic curve segments is equal.

For a typical coil the first assumption is justifiable since mainly the acting force of gravity determines the orientation of the semi major axis to be horizontal. Further the requirement of constant ellipticity for all elliptic curve segments doesn’t hold here exactly and will cause a small error, but the center estimation is not affected by this fact (Schleicher & Zagar, 2009b).

First Eqn 12 is rewritten using the center coordinate $$\left(x_c,y_c\right)$$and the semi axes a and b.

$$\begin{array}{l}{F}_1\left(x,y\right)=\frac{1}{a^2}{\left(x-x_c\right)}^2+\frac{1}{b^2}{\left(y-y_c\right)}^2-1=0\\ =\underset{A}{\underbrace{\frac{1}{a^2}}}{x}^2+\underset{C}{\underbrace{\frac{1}{b^2}}}{y}^2\underset{D}{\underbrace{-\frac{1}{a^2}2x_c}}x\\ \underset{E}{\underbrace{-\frac{1}{b^2}2y_c}}y+\underset{F}{\underbrace{\frac{1}{a^2}{x}_c^2+\frac{1}{b^2}{y}_c^2-1}}=0\end{array}$$E26

Introducing $e=\frac{a}{b}$ and factorizing $\frac{1}{b^2}$ leads to

$$\begin{array}{l}{F}_1\left(x,y\right)=\frac{1}{b^2}(\underset{A^{\prime }}{\underbrace{\frac{1}{e^2}}}{x}^2+\underset{C^{\prime }}{\underbrace{1}}{y}^2\underset{D^{\prime }}{\underbrace{-\frac{1}{e^2}2x_c}}x\\ \underset{E^{\prime }}{\underbrace{-2y_c}}y+\underset{F^{\prime }}{\underbrace{\frac{1}{e^2}{x}_c^2+y_c^2-b^2}})\end{array}$$E27

Concerning Eqn 27 the newly introduced parameters $$A^{\prime },...,E^{\prime }$$are independent on b, while $$F^{\prime }$$contains the curve dependent parameter b. Therefore $$\Theta$$ must now be written as

$$\Theta ={\left[A^{\prime },C^{\prime },D^{\prime },E^{\prime },{F^{\prime }}_1,\mathrm{...},{F^{\prime }}_j,\mathrm{...},{F^{\prime }}_J\right]}^T$$E28

using all curve segments in the dataset _s, which is indexed as follows: The dataset contains J curves indexed by j, where curve j has M_j points. A point on curve j is denoted by$$\left(x_j{}_{,m},y_j{}_{,m}\right)$$, where m is the index of the point (m = 1,...,M). $$X_i={\left[x_{j,m},y_{j,m}\right]}^T$$

When concatenating the curves to the input data vector x a point is addressed by index i.

$$\begin{array}{l}x={\left[\begin{array}{l}{x}_{1,1},x_{1,2},\mathrm{...},x_{1,M_1-1},x_{1,M_1},\mathrm{...},x_{j,1},\mathrm{...},x_{j,M_j},\mathrm{...},x_{J,1},\mathrm{...},x_{J,M_J}\\ y_{1,1},y_{1,2},\mathrm{...},y_{1,M_1-1},y_{1,M_1},\mathrm{...},y_{j,1},\mathrm{...},y_{j,M_j},\mathrm{...},y_{J,1},\mathrm{...},y_{J,M_J}\end{array}\right]}^T\\ =\left[\begin{array}{l}{x}_1,\mathrm{...},x_i,\mathrm{...},x_I\\ y_1,\mathrm{...},y_i,\mathrm{...},y_I\end{array}\right]\end{array}$$E29

Furthermore u(x) must also be extended to

$$u(x_i)={\left[x_{j,m}^2,y_{j,m}^2,x_{j,m},y_{j,m},\mathrm{0...},\underset{1at\overline{)}index(j+4)}{\underbrace{1}},\mathrm{...},0\right]}^T$$E30

In the general sense we can write

$$u(x_i)={\left[x_{j,m}^2,y_{j,m}^2,x_{j,m},y_{j,m},\delta (1-j),\mathrm{...},\delta (J-j)\right]}^T$$E31

For the gradient we obtain

$${\nabla }_{x}u(x_i)={\left[\begin{array}{l}2x_i,0,1,0,0,\mathrm{...},0\\ 0,2y_i,0,1,0,\mathrm{...},0\end{array}\right]}^T$$E32

Having $$u\left(x_i\right)$$given the matrices $$A_i$$and $$B_i$$are calculated and $$\Theta$$ is estimated using one of the estimators described in Sec. 3.5.1.

Out of all parameters in $$\Theta$$ jointly estimated the center coordinate (x_c,y_c) and parameter and b are calculated.

First the parameter vector $$\Theta$$ is normalized so that $$C$$is set 1. By comparison of coefficients the center coordinate $$\left(x_c,y_c\right)$$are obtained:

$$x_c=\frac{D^{\prime }}{2A^{\prime }}$$E33

$$y_c=-\frac{E^{\prime }}{2}$$E34

Since the other parameters and b

$$e=\frac{1}{\sqrt{A^{\prime }}}$$E35

$$b=\sqrt{A^{\prime }{x}_c^2+y_c^2-F^{\prime }}$$E36

$$=\sqrt{A^{\prime }\frac{{D^{\prime }}^2}{2{A^{\prime }}^2}+\frac{{E^{\prime }}^2}{4}-F^{\prime }}$$E37

are biased due to the (invalid) assumption of equal ellipticity of all segments the simple ellipse fit (Sec.3.5.3) is done for each curve separately, using the previously estimated center coordinate$$\left(x_c,y_c\right)$$.

The result of the concentric ellipse fit followed by the simple ellipse fit, where the ellipticity $$e_j$$for each curve segment is calculated (shown in Fig.13) and discussed after the error analysis in Sec. 5.

4. Error analysis

In the following sections the main error sources which influence the accuracy of the estimated ellipticity will be discussed.

4.1. Aliasing

When sampling in spatial coordinates the aliasing effect (Jähne, 2002; Burger & Burge, 2006) is similar to consider as the case in the time domain. To be able to reconstruct a sampled signal the sampling frequency $$f_s$$must be larger than twice the largest occurring frequency component$$f_c$$. In case of the coil image the thinnest sheet metal layer is 0.5 mm, which means$$f_c=0.\text{5 mm}/\text{layer}$$. Since the Nyquist criterion is also valid for any orientation the spatial frequencies $$f_c$$must be multiplied by $\sqrt{2}$

$$f_{cd}=f_c\sqrt{2}$$E38

For a coil with a diameter of 2 m a minimum camera chip size of

$$\frac{2000mm}{0.5mm/layer}\sqrt{2}\cdot 1pixel/layer\approx 5657pixel$$E39

in horizontal and vertical direction, respectively, is needed. This corresponds to a 32 Mpixel camera chip. Even if we neglect $\sqrt{2}$ –factor for the diagonal layers, we still would require a camera chip with 4000 × 4000 = 16 Mpixel.

Luckily the small axial displacement of the sheet metal layers with respect to each other causes some shadings which can be captured very good avoiding aliasing. Also the roughly–textured cut surface of most thick sheet metal layers improves the possibility to recognize a single layer. In the worst case when a coil has thin layers, these layers appear blurred in the image.

Further the goal is still to estimate the ellipticity and not to capture each layer. Assuming that the arc segment is recognized only by its shadow which extends from one layer to the next (one pixel = 2.5 mm in this setup), the relative error for the ellipticity is bounded by

$$f_e=\frac{\Delta e}{e}=\frac{\frac{a+\Delta a}{b}-\frac{a}{b}}{\frac{a}{b}}=\frac{\Delta a}{a}=0.67\%$$E40

with the major semi axis $$a=0.\text{3 m}$$and the arc displacement$$\Delta a=\text{2}.\text{5 mm}$$. Considering $$f_e$$this error is quite large for high precise ellipticity estimation.

4.2. Camera orientation

One possible error source is the miss–orientation of the camera. As described previously the camera chip plane should be oriented strictly parallel to the coils front side. Even when the setup is build up very carefully using a spirit level and using edges of the coils carrier for the adjustment, an error free orientation can not be obtained. For a circular coil, with radius r and a misalignment angle$$\phi$$, a rough error estimate can easily be obtained:

$$\Delta e=\frac{r}{r\cos \phi }-1=\frac{1}{\cos \phi }-1$$E41

The angle $$\phi$$ can have any direction, since the semi major axis of the resulting ellipse is always $$r$$and the semi minor axis is$$r\text{cos}\phi$$. Assuming $$\phi =\text{ 1}°$$leads to an error of 0.0152% and even a larger angle of 5 , which should be possible to underbid when adjusting the camera, leads only to 0.55%.

Figure 11.

Between each $$\Delta s$$thick sheet metal layer is a small oil gap $$\Delta y$$and $$\Delta x$$in x and y direction. a_n and b_n are the semi major and minor axes of the n^th layer.

4.3. Ellipticity–radius dependency

The assumption of equal ellipticity for all arc segments made in Sec. 3.5.4 leads to a erroneous . Knowing the dimensions of the coil a simulation of the ellipticity progression can be done using the following model.

Assume an ideal circular coil with the inner radius$$r_i$$, outer radius$$r_o$$, metal sheet thickness Δs and initial oil gap $$\Delta r$$between the layers (see Fig.11). Approximating the layers of the coil, which in reality are forming an so–called Euler’s–spiral, by concentric circular layers, the radius of the $$n^{\text{th}}$$layer $$r\left(n\right)$$is

$$r(n)=r_i+n(\Delta s+\Delta r)$$E42

The assumed maximum number of turns on the coil is

$$n_{\max }=(r_0-r_i)/(\Delta s+\Delta r)$$E43

When applying a virtual vertically acting force (e.g. gravity) on this coil, it is compressed somewhat flat, the layers get elliptic and the oil gap $$\Delta y\left(n\right)$$ decreases in vertical direction while $$\Delta x\left(n\right)$$increases in horizontal direction. Optionally, the oil gap can be modeled to be smaller at the middle layers, since there is a higher inner pressure caused by the winding action in production, by following equation:

$$\Delta y(i)=\Delta s\left[\alpha +(1-\alpha )\frac{4{(i-n_{\max }/2)}^2}{n_{\max }^2}\right]$$E44

Or, since the oil gap is very small in ratio to the metal sheet thickness, to be a constant setting$\alpha =1$. For the vertical minor semi axis of the coil we obtain:

$$b(n)=r_i/\beta +n\Delta s+{\displaystyle \sum _{i=0}^n\Delta y(i)}$$E45

where $\beta =\frac{r_i}{b(0)}$ is the compression ratio of the inner radius.

Since the mass of the circular coil m_c is equal to the mass of the deformed coil$$m_d$$, also the cut surface of the metal sheet layers must be the same ($$ñ$$... density of steel,$$w$$... width of the coil).

$$m_c=m_d$$E46

$$\tilde{n}w{A}_c=\tilde{n}w{A}_e$$E47

The areas are calculated by

$$A_c(n)=\pi ({(r(n)+\Delta s)}^2-r{(n)}^2)$$

$$A_e(a(n),n)=\pi ((b(n)+\Delta s)\cdot (a(n)+\Delta s)-b(n)a(n))$$E48

Resolving $$a\left(n\right)$$using Eqns 48, 49 and 50 we obtain

$$e(n)=\frac{a(n)}{b(n)}=f(r_i,r_0,\Delta s,\Delta y(n),\beta ,n)$$E49

Figure 12 shows two simulations of the ellipticity dependency over the turn number for various coil parameter. The layer dependent oil gap $$\Delta y\left(i\right)$$has the most effect on coils with thin metal sheets (see Fig. 12c vs. Fig. 12d). Generally the ellipticity decreases strongest within the first third of the radius. Setting $$\Delta y\left(i\right)=0.00\text{1 mm}$$and smaller gives similar results like shown in Fig. 12a. Since the exact conditions of the pressure distribution in the oil gap is unknown and no proper literature could be found considering this problem, this simulation can only give a qualitative evidence for the ellipticity progression.

4.4. Estimation errors

A comprehensive analysis of the estimation errors was done by Schleicher & Zagar (2009b), whose results are summarized in this section shortly. Using synthetic curve datasets, similar to $${\text{C}}_s$$the different estimators where evaluated. First the CEF algorithm was applied onto 1500, 5000 and 20000 synthetic datasets. Then the Euclidean distance between d the true

Figure 12.

Simulation of the ellipticity $$e=f\left(n\right)$$for different inner ellipticities. (Δy = 0.01 mm,r_i = 0.3 m,r_o = 1 m, β = (1.005, 1.010, 1.015, 1.020, 1.025,1.030))

center $${\overline{x}}_c=\left({\overline{x}}_c,{\overline{y}}_c\right)$$ and estimated center $$x_c=\left(x_c,y_c\right)$$was calculated, which corresponds to the error of the center estimation.

$${\mu }_d=E\left\{{\Vert \left({\overline{x}}_c-x_c|{}_{A\lg .}\right)\Vert }_2\right\}$$E50

$${\sigma }_d=\sqrt{Var\left\{{\Vert \left({\overline{x}}_c-x_c|{}_{A\lg .}\right)\Vert }_2\right\}}$$E51

With increasing number of datasets $${\mu }_d$$decreases, justifying and proofing that the CEF algorithm has no bias. Using the best performing HEIV algorithm an error of 0.00251±0.251 pixel distance to $${\overline{x}}_c$$ the was obtained. The variance was very similar for all algorithms and is probably limited by the quantisation of the dataset itself.

Furthermore the SEF algorithm was examined, concerning an erroneous center coordinate as input and arc segments with various center angles and lengths. The result was, that for a center displacement of (1,1) (which is very large in relation to$${\sigma }_d$$) the relative ellipticity error $f_e=\frac{\Delta e}{e}$ was smaller than 0.2 ± 0.3% for all algorithms (except Taubins) when curve segments with a length longer than 1/25 of their circumference where used. Furthermore arcs symmetric around the ellipses apex lead to larger errors than others.

Figure 13.

Elliptic arc segments are extracted from the image of the coils front side, and used to estimate the ellipticity.

Figure 14.

The estimated ellipticities of the arc segments are plotted over the major semi axis.

5. Results

Figure 13 shows the results of the algorithm described in Sec. 3. The estimated parameters of the extracted curve segments are used to overlay the ellipses onto the image of the coils front side. Furthermore the estimated ellipticities $$e_j$$are plotted over their major semi axis $$a_j$$in Fig. 14. Moreover some possible figures of merits, which could characterize the coils quality, are plotted. Depending on the extracted curve segments, the model of ellipticity progression derived in Sec. 4.3 fits more or less perfect. The main error source results from the low camera resolution, which might be changed in future when the cost of high resolution cameras decrease.

6. Conclusion

This chapter introduced the problem of ellipticity estimation on hand of a practical problem. To obtain a good ellipticity estimation proper datasets must be obtained. These datasets are generally not very easy to obtain, for which reason also the preprocessing to extract the datasets was discussed in detail. Since only arc segments could be used for the ellipticity estimation, normal ellipticy fits fail. Hence the estimation was split into two parts, first to estimate the center of the coil with the concentric ellipticity estimator and then to estimate the ellipticity of the arc segments. Using the general approach for the ellipse fit, enables also to use different solver implementations for fitting problems. When discussing the errors, it was accented that a lot of error sources affect the final result. The correct understanding of the partial errors helps to optimize the corresponding section of the algorithm.

Acknowledgments

At first the authors like to thank N. Chernov and Tao Peng for their support. Furthermore the authors gratefully acknowledge the partial financial support for the work presented in this paper by the Austrian Research Promotion Agency under contract grant 814278 and the Austrian COMET program supporting the Austrian Center of Competence in Mechatronics (ACCM). Last but not least the authors thank the company Industrie Logistic Linz (ILL) for its financial and technical support.