Verifying Neutron Tomography Performance using Test Objects

Available online at www.sciencedirect.com Physics Procedia 43 (2013) 128 – 137 Verifying neutron tomography performance using test objects A.P. Kaestnera , E.H. Lehmanna , J. Hovinda , M.J. Radebeb , F.C. de Beerb , C.M. Simc a NeutronImaging and Activation Group, Paul Scherrer Institut, CH-5232 Villigen, Switzerland b RadiationScience, NECSA, P1500, P.O Box 582, Pretoria, 0001, Republic of South Africa c Korea Atomic Energy Research Institute, 1045 Daedeok-daero, Yuseong-gu, Republic of Korea Abstract In an effort to provide a standardized method to quantify the imaging capabilities of neutron imaging beam-lines worldwide, we propose a set of test objects for neutron tomography. The test objects are designed to quantify spatial resolution and material contrast in tomograms. The resolution samples aim at detecting a thin film embedded in a different material. Two samples with complementary material compositions are proposed for this purpose. The contrast sample has several insets of different materials. The measurements are proposed to be done using both radiography and tomography. The image processing methods needed to evaluate the performance of the reconstructed data are presented. The methods are automated to avoid subjective decisions by persons who evaluate the data. Experimental data to demonstrate the test objects and their analysis methods were acquired with the cold neutron imaging beam-line, ICON, in comparison with data from the thermal neutron facility, NEUTRA at Paul Scherrer Institut, Switzerland. This is a first initiative and is open for discussion among the participants to further improve the evaluation procedure. Keywords: Neutron imaging, Computed tomography, Quality assessment, Image evaluation, Test bodies 1. Introduction Neutron imaging is often used for the non-destructive investigations of samples that are unsuitable for X-ray imag- ing. Neutron imaging is based on the transmission of a radiation through an object like X-ray imaging. This means that the shadow of the semitransparent object reveals the inner structure of the sample. The difference between the two imaging type is that the attenuation coefficients are very different from each other. Computed tomography (CT) is also a possibility with neutron imaging. Similarly to X-ray CT, a tomogram is acquired by rotating the object while simultaneously acquiring radiograms from different views of the sample. This makes it possible to reconstruct the three-dimensional structure of the sample. For X-ray imaging the development of testing procedures has been driven by both the communities of medical imaging and non-destructive testing. Since neutron imaging in principle has the same requirements in repeatability and image quality similar testing procedures should be developed. There are few methods to measure the quality of the imaging capabilities of neutron imaging facilities. In general, they focus on quantifying the resolution of the radiograms. However, since about 50% of the neutron imaging beam time is used for tomography it is motivated to define a method that aims at quantifying the quality of the reconstructed tomography data. Here, two objectives are in focus: the ability the detect a thin film embedded in a sample and to investigate the contrast distribution for different well defined materials. In both cases the signal to noise ratio (SNR) is a relevant quantity since this has a direct impact on the performance of the measurements. In this paper we present two samples to measure the detectability of thin films embedded in a sample and one sample to quantify the attenuation coefficients of six different materials. In addition to the sample definition we also describe the automated image processing steps required to analyse the image data objectively. As the data is to be analysed by different operators world wide it is important to assure that the analysis is done in a repeatable manner. Email address: [email protected] (A.P. Kaestner) 1875-3892 © 2013 The Authors. Published by Elsevier B.V. Open access under CC BY-NC-ND license. Selection and/or peer-review under responsibilty of ITMNR-7 doi:10.1016/j.phpro.2013.03.016  A.P. Kaestner et al. / Physics Procedia 43 (2013) 128 – 137 129 This is an international initiative initiated by a working group of the International Atomic Energy Agency (IAEA) and involves the collaboration between Paul Scherrer Institut, Switzerland, NECSA, South Africa, and the Korea Atomic Energy Research Institute, South Korea. In section 2 we describe the proposed samples, acquisition procedure, and the required image processing steps to quantify the quality parameters from the image data. The procedures and evaluation methods are then experimentally verified by scanning the samples at the ICON beam-line, section 3. In the results section the results from the described image processing steps are presented. 2. Method 2.1. Thin film samples The intention of the thin film sample is firstly to determine the thinnest film that is possible to detect by the used imaging system. There are several factors that affect the detectability; pixel size, point spread functions of the detector system and reconstruction software, beam divergence, and scattering. When the pixel size is greater than the film thickness it may still be visible due the partial volume effect. The consequence is that the attenuation coefficient is not correctly represented for the film since each voxel contains a mix between film and sample body. The imaging system introduces a blurring of the acquired neutron shadow due to light spread in the scintillator and unsharpness in the light optical system. These effects contribute to the point spread function. The beam divergence plays a role in the blurring of the film detectability, especially when the film is thin relative to the divergence angle. Finally, sample scattering can also affect the attenuation coefficient of the film. In CT this effect appears with a cupped intensity profile across the sample. These factors makes it hard to measure the film thickness and quantify the attenuation coefficient for thin films. Still, the film may be clearly visible which may be a sufficient indicator for some experiments. In this paper we will only quantify the detectability of the film and measure the line spread of the film. The resolution is better determined using a contrast step from which the edge spread function and modulation transfer function can be derived. The thin film sample has two complementary configurations. One to determine the detectability of a thin low- contrast feature embedded in a high-contrast body. In the other configuration the contrast for film and body is reversed, i.e. high-contrast film and low-contrast body. The test object consists of two blocks between which one or more thin films are pressed, as illustrated in figure 1. The blocks are assembled with screws, this allows the most tight press of the inserted films. We chose to use the two combinations one with Fe blocks with inserted Al films and the other with Al blocks with Cu films. The film thickness in both cases is 20 μm. Using films, the gap can only be changed in discrete steps of the film thickness. Figure 1. The resolution sample. The colored lines mark the slices that will be evaluated. 2.1.1. Image acquisition The measurement procedure is divided into two steps: 1, identify the thinnest visible film using radiograms and 2, perform a computed tomography. The plots in figure 2 show that the system response to the inserted film strongly depends on the observation angle. A consequence is that if a CT scan is started at an arbitrary angle the film may or 130 A.P. Kaestner et al. / Physics Procedia 43 (2013) 128 – 137 may not produce a significant signal in the acquired projection data set. Therefore, we suggest to identify the starting angle of the CT scan using a rocking scan. This initial scan is made with a higher angular resolution than intended or needed for the following CT scan. The image that provides the greatest contrast for the film is the one where the film is most parallel aligned to the neutron beam. By starting at the position with the greatest contrast the greatest probability of film detection will also be provided. A side effect is that the reconstructed slices do not need to be rotated for the analysis, the film is already parallel to an image axis. If the thickness of the gap is less than the resolution of the imaging system a partial volume effect will be observed, i.e. the film may be detected but the thickness cannot be determined. This gives two different metrics to observe. A weak one, that essentially measures the detectability of a thin film, and a stronger one that is able to quantify the resolution of the system. In the following analysis we will only focus on the detectability. If the gap is too narrow it cannot be detected and an additional film must be inserted between the sample blocks. 1 0.9 0.8 0.7 Transmission [1] 0.6 Al−Cu 0.5 Fe−Al 0.4 0.3 0.2 0.1 0 −5 −4 −3 −2 −1 0 1 2 3 4 5 Observation angle [deg] Figure 2. Transmission through film region the two resolution test bodies at different observation angles. 2.1.2. Reconstruction of the projection data The reconstruction can be done using any available CT reconstruction software. It is however important that the output is given in floating point format or that conversion parameters from integer gray levels to linear attenuation coefficients are provided with the integer data. The best reconstruction is achieved when the center of rotation is defined to sub pixel accuracy and when possible tilt of the acquisition axis is corrected. Additionally, it is recommended to apply filters to remove typical artefacts from the projection data (ring and line artefacts). If scattering and beam-hardening correction methods are used, this will improve the performance of the evaluation. 2.1.3. Image analysis The purpose of this analysis procedure is to determine the detectability of the film and to compare the width of the film in the reconstructed slice with the known film thickness. Select a slice f from the reconstructed data. It is recommended to select slices from the three regions marked in figure 1 to verify the spatial dependency on the unsharpness. The following description assumes that the film is oriented parallel to the y-axis in the slice image. The analysis to be made on a one dimensional profile perpendicularly across the film. To improve the statistics of the profile it is computed from N lines. The first step of the profile extraction is to locate the sample center in the image. This is done directly on the gray-level image by computing the centroid in the x and y directions providing {c x , cy } as Ny Nx y=1 x=1 x f (x, y) c x = N N (1) y x y=1 x=1 f (x, y)  A.P. Kaestner et al. / Physics Procedia 43 (2013) 128 – 137 131 where N x and Ny are the number of pixels in each direction of the image. Computing cy is done analogously. With the profile  N/2 t(p) = f (p + c x , q + cy ), −u ≤ p < u (2) q=−N/2 where 2 u is the width of the extracted profile and N the number of lines to average. The width u should be selected several times greater than the width of the film. The number of lines is strongly connected to the SNR and the contrast between film and sample body. Larger N is required for low SNR images and is one indication of the detectability of the film. For this evaluation the mean (μ = E[t]) and standard deviation (σ = s[t]) are needed. If the segmented signal d given by ⎧ ⎪ ⎪ ⎨1 if 2 σ < |t(p) − μ| d(p) = ⎪ ⎪ (3) ⎩0 otherwise only contains a single cluster of adjacent pixels the film is considered to be detectable for the selected number of lines N. Count the number of detected pixels (N f ilm ) when the detectability criteria is fulfilled and compare the detected film width with the inserted film thickness (t f ilm ) using the line spread ratio N f ilm · [pixel size] LS R = (4) t f ilm The value of LSR is smaller for a high resolution system and approaches unity for an ideal system. In addition to the line spread ratio also determine the contrast between the sample body and the film as C FeAl = min(t) − μ (5) or C AlCu = max(t) − μ (6) depending on which sample configuration is used. 2.2. Contrast sample The contrast sample is a cylinder of aluminum with a concentric ring of six insets of different materials as shown in figure 3. This sample is intended for evaluatation with computed tomography. The purpose of this sample is to measure the contrast between different materials and the embedding material. The sample can also be used to quantify the attenuation coefficients of the inset materials. The sample has a diameter of 30 mm and the insets have a 6 mm diameter. The insets are uniformly distributed over a concentric ring with a diameter of 18 mm. These dimensions are chosen to fit in the field of view (FOV) at most neutron imaging facilities. Two contrast sample configurations have been made. In the first version the sample contained an inset of polyethy- lene (PE). This inset proved in the initial experiments to be too strongly attenuating which resulted in very strong starvation artefacts, see [1] for an overview of common artefacts occurring in CT slices. The revised second version has the PE inset replaced by an inset of Ti. The organization of the insets and their corresponding linear attenuation coefficients are listed in table 1. Figure 3. The layout of the contrast sample. 132 A.P. Kaestner et al. / Physics Procedia 43 (2013) 128 – 137 1 2 3 4 5 6 Configuration I Al (0.11) Cu (1.0) Pb (0.37) PE (6.4) Fe (1.2) Ni (2.1) Configuration II Al (0.11) Cu (1.0) Pb (0.37) Ti (0.59) Fe (1.2) Ni (2.1) Table 1. Inset configuration of the proposed contrast samples. The insets are listed clockwise. The numbers are theoretical linear attenuation coefficients for thermal neutrons [cm-1 ]. The data acquisition of the contrast sample is done using a typical CT setup. The number of projections required is related to the sampling theorem, i.e. the more pixels the sample covers the more projections are required. A pragmatic rule of thumb is to use as many projections as the width of the sample in pixels. With more projections, a better signal to noise ratio (SNR) will be obtained in the reconstructed data. The SNR depends on the total acquired neutron dose, hence the exposure time of each projection also counts into the relation for contrast and noise in the reconstructed data. A guideline for the number of projections N pro jections needed to achieve a slice contrast C slice with a given projection contrast C pro jection for a sample that is W sample pixels wide is given by C slice W sample = C pro jection N pro jections (7) The slice contrast relates to the number of gray levels that cover the entire image dynamics. A consequence of equation 7 is that projections acquired with low gray level dynamic can produce higher contrast when the number of projections is increased. The equation is not an absolute truth since the effects of noise and scattering are not considered. Still, it gives an indication how the quantities are related. The number of projections must also be defined by the sampling theorem otherwise aliasing artefacts will contribute to the noise level. The projection data from the contrast sample scan shall be reconstructed as described in 2.1.2. 2.2.1. Image analysis of the contrast sample The evaluation of the contrast sample is done slice-wise at three different locations of the sample – one slice at the upper 10% of the sample, the central slice, and the last at the lower 10% of the sample. These the locations are evaluated to identify spatial variations in the estimated values of the linear attenuation coefficients. The attenuation coefficient for each inset is determined as the average of the pixels in a central circular region of interest (ROI) with a diameter corresponding to 80% of the inset diameter in pixels. At the same time the standard deviation is to be determined using the same ROI. The next step is to measure the local background level. This is done using a annulus-shaped ROI centered on each material inset. The area of the annulus shall be the same as the area of the inset ROI. We propose to perform an automated evaluation of the slices to avoid biased results due to subjective choices. This is important since the data from different facilities worldwide will be evaluated by different persons. Some a priori information is required for this evaluation, most important is the pixel size of the slice. The remaining information like inset diameter and the diameter of the inset ring is given by the sample geometry. The first step of the inset localization is to filter the slice image f using a median filter (7×7) to suppress the effect of disturbing artefacts and noise. fmed = median7×7 ( f ) The blurring caused by this step will not affect the measurements in the original image at later stages. It only helps to improve the inset center estimation. To locate the insets in the slice image a detection image is computed by K convolution using a ring shaped kernel defined as Kring = x,yring K where ring ⎧  ⎪ ⎪ ⎨1, if | x2 + y2 − R| ≤ 1  Kring (x, y) = ⎪ ⎪ 2 ⎩0, otherwise The radius of the ring R is given in pixels and is determined using the known pixel size and the metric radius of the inset. The convolution (denoted by ∗) with Kring provides the detection image d = |Kring ∗ fmed |  A.P. Kaestner et al. / Physics Procedia 43 (2013) 128 – 137 133 In the detection image the centers of the insets will be indicated by cone shaped local maxima, these have to be detected. This is done using the morphological HMAX operator [6] which extracts all local maxima with a height of at least h relative to the local background. The inset centers, c, are identified by thresholding the local max image h c= < (d − HMAX(d, h)) 2 The height parameter h is determined as an intensity difference near the max intensity of d. Here we used the distance 15%R from the most intense inset center as reference point. h = d(p) − d(p + (0, 0.15 R)) where p is the position of the intensity maximum. The individual centers are identified using connected component labelling of c resulting in the image clbl . This procedure may not identify all insets due to low contrast, especially the Al inset is not distinguishable from the background, but also the Pb inset may be hard to locate. To overcome this limitation we propose to estimate the locations of these insets by parametrize the circle of insets, see appendix A, and use the known order of the identified insets to fill in the missing inset centers. The final set of centers are denoted by the coordinates ci = (xi , yi ) where i indicates the location in the sequence of insets. Once the centers of each inset are identified the measurements can be done using the following definition for the region of interest ROIr = {p = (x, y) | x2 + y2 ≤ r; −r ≤ x, y ≤ r; (x, y) ∈ Z} Then, M sample = ROI80%R represents the region inside the inset and the ring outside the inset is defined by the set difference between two discs as Mbackground = ROI136%R \ ROI110%R as shown in figure 4. With these ROIs defined we Figure 4. The ROIs used to compute the inset statistics. can compute the average intensity and standard deviation as 1  Iisample = E[ f (p)|p ∈ M sample + ci ] = f (p) (8) #M sample p∈(M sample +ci ) 1  Iibackground = E[ f (p)|p ∈ Mbackground + ci )] = f (p) (9) #Mbackground p∈(Mbackground +ci ) 1  sisample = s[ f (p)|p ∈ M sample + ci )] = ( f (p) − Ii )2 (10) #M sample − 1 p∈(M sample +ci ) The subindex i indicates the material as listed in table 1. It is now possible to compare the reconstructed attenuation co- sample sample background Ii Ii −Ii efficient with the expected (using Iisample ), compute the S NR = and contrast to noise ratio CNR = . sisample sisample The results for each material will be used in the final facility performance report. 134 A.P. Kaestner et al. / Physics Procedia 43 (2013) 128 – 137 3. Experiments The three samples were scanned at the cold neutron imaging beam-line ICON [3] at Paul Scherrer Institut. The FOV was set to 40×40 mm and an Andor DV-434 CCD camera (1024×1024) was used. A Zeiss Macro Planar 100mm f/2.0 lens was focussed on the 100 μm thick 6 LiF scintillator screen. The collimation ratio was set to L/D=342. The samples were scanned using 375 projections uniformly distributed over 360◦ . The exposure time was set to 15 s, which provides a gray level dynamic of about 32000 levels. Longer exposure times would saturate the camera chip. The resulting projection data was reconstructed using the MuhRec reconstruction software [2]. The angle at which the film appeared with the greatest contrast was identified before each scan of the resolution samples described in section 2.1. The scans were started without removing the sample from the identified position. 4. Results 4.1. Thin film samples For this evaluation we only present the results from the Fe sample with Al films. The reconstructed images were analyzed using a rectangular ROI which was located by the procedure described in section 2.1.3. With the used pixel size (37.6μm) a single film was faintly visible as shown in figure 5. This film line is hard to separate from the background by thresholding since the noise produces outliers that are detected as film. With the non-ambiguity criteria from section 2.1 this sample configuration is not resolved by the current setup. One reason is the cupping caused by sample scattering that has a negative impact on the profile statistics since the curved profile indicates a greater noise variance than is actually true. This over-estimate of the variance may make the detectability worse than it is. The profile in figure 5 shows an example when the cupping disturbs the detectability test in a slice with a single 20 μm foil. 0.93 0.92 100 0.91 200 0.9 Intensity 0.89 300 0.88 400 0.87 0.86 500 0.85 600 0.84 100 200 300 400 500 600 1 2 3 4 5 6 7 Position [mm] Figure 5. A profile from the Fe sample with a single 20 μm films inserted. Cupping makes unambiguous detection impossible. By adding a second film to the sample the film was more pronounced and the detectability increased, figure 6. The profile from N=100 lines shows that the CNR improved a lot as an effect of less partial volume data in the line. The width of the detected line using 2σ as detection criteria was compared with the thickness of the film stack. In this evaluation N=100 lines on a single slice was used. A more localized alternative would be to compute the profile from several slices but with less lines on each slice. That would provide the same statistics, but much more localized which would be more realistic since thin features are rarely flat.  A.P. Kaestner et al. / Physics Procedia 43 (2013) 128 – 137 135 0.92 100 0.9 200 0.88 Intensity 300 0.86 400 0.84 500 0.82 600 100 200 300 400 500 600 0 2 4 6 8 Position [mm] Figure 6. A profile from the Fe sample with an two 20 μm films inserted. The vertical red lines indicate the width of the detected film. The rectangle in the slice image marks the region used to compute the profile. 4.2. Contrast Three reconstructed slices were used in the evaluation of the contrast sample. The slices are shown in figure 7. The first two slices were measured at NEUTRA and should show essentially the same contrasts for the different insets, with exception for one inset that is PE in slice (a) and Ti in slice (b). It should be noted that the attenuation by the PE inset in combination with one of the insets Fe, Cu, or Ni is too strong. This causes strong artefacts in the slice. These artefacts also affect the results of the evaluation since the reconstructed attenuation coefficient is biased by the effect of the detector starvation. The last slice that was measured at ICON shows that the beam spectrum also plays a role for the evaluation. The total cross section for Fe, Cu, and Ti approach each other for a colder spectrum making them hard to identify when they are combined. (a) (b) (c) Figure 7. The reconstructed slice that were used in the evaluation. Slices a and b are from NEUTRA and slice c is from ICON. Using the automated evaluation scheme described in section 2.2 provided the results presented in table 2. Each image produced six measurements of image intensities and contrast to noise ratio. Tests to evaluate the robustness of the evaluation scheme showed that the method is insensitive to rotations, translations, and variations in image dynamics. The only parameter that is crucial for a successful evaluation are pixel size since it is used to determine the inset diameter. Mostly the inset detection finds four of the six insets, while the two weakest (Al and Pb) have too low contrast to be detected. In general, three insets are sufficient to find the center of the inset ring. 136 A.P. Kaestner et al. / Physics Procedia 43 (2013) 128 – 137 Al Cu Pb PE Fe Ni Ti NEUTRA 0.10 (0.16) 0.88 (13.4) 0.30 (4.43) 2.61 (11.7) 0.79 (15.4) 1.43 (25.5) – NEUTRA 0.11 (0.54) 0.78 (18.0) 0.24 (4.03) – 0.93 (19.4) 1.46 (26.2) 0.52 (11.6) ICON 0.09 (0.34) 1.00 (13.1) 0.16 (1.30) – 0.98 (12.1) 1.71 (18.7) 0.86 (10.9) Table 2. Measured attenuation coefficients in cm-1 from contrast sample I and II. The value in parentheses is the contrast to noise ratio. The data was acquired at NEUTRA and ICON. 5. Discussion In this paper we presented two sample types for the characterization of the configuration of a neutron tomogra- phy setup. In addition to the samples we also describe the automated image analysis procedures that improves the repeatability of the image evaluation. First experiments with measured data from the two neutron imaging beam lines NEUTRA and ICON at Paul Scherrer Institut have shown that the first concept still requires some refinements. Further experimental results are presented in [5], in this work an outline for a standard measurement procedure is presented. The geometry of the test objects in the current set exclude experiments with some instrument setups at the beam- lines. In particular those with high resolution detectors since there the field of view is mostly smaller than the samples and can not include sample in a single view. A proposal for a new version would be to reduce the greatest sample width to 25 mm for all samples. Ideally, two sets should be defined to be able to handle two levels of resolution. In addition to the geometry of the samples, it has also come clear that the neutron spectrum plays an important role in the evaluation and comparison of the results. This is most important for the contrast sample. The original design was made with thermal neutrons in mind. With thermal energies the attenuation coefficients of the insets are well distributed, while the attenuation coefficients are more clustered with a cold spectrum. Determining the optimal number of projections for the different experiments may require many scans in the testing phase of the routine. To save beam time for this investigation we propose to use a scanning protocol based on angular increments using the Golden ratio[4]. Using this scanning protocol a single scan can provide all data needed for a reconstruction for any number of projections up the acquired total number of projections. The samples can also be used to verify the performance of reconstruction software and artefact correction methods. One example is spot cleaning algorithms detect outliers in the projection data. This is a problem with the thin film of the resolution sample since the film could be identified as a spot with devastating results on the reconstructed images. The proposed set of test objects and automated analysis procedure are a part of the development toward a standard on how to determine the performance of the used experimental conditions. A round robin procedure to send sets of samples to different neutron imaging facilities world wide has been started. The results from these experiments will be compiled and reported to the International Atomic Energy Agency (IAEA) when all involved institutes have returned their results. First feedback and results have already returned from contributing facilities. From this information it has become clear that a very strict testing procedure is difficult to define since there are so great variations in the facility configurations. These variations include neutron spectrum and different boundary conditions for the detector/image acquisition systems. References [1] J. Barret and N. Keat. Artifacts in ct: Recognition and avoidance. RadioGraphics, 24(6):1679–1691, 2004. [2] A.P. Kaestner. MuhRec – a new tomography reconstructor. Nuclear Instruments and Methods A, 651(1):156–160, September 2011. [3] A.P. Kaestner, S. Hartmann, G. Kuehne, G. Frei, C. Gruenzweig, L. Josic, F. Schmid, and E.H. Lehmann. The ICON beamline - A facility for cold neutron imaging at SINQ. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 659(1):387 – 393, 2011. [4] A.P. Kaestner, B. Münch, P. Trtik, and L.G. Butler. Spatio-temporal computed tomography of dynamic processes. Optical Engineering, 50(12):123201, 2011. [5] M.J. Radebe, F.C. de Beer, A. Kaestner, E. Lehmann, C.M. Sim, and E. Sideras-Haddad. Evaluation procedures for spatial resolution and contrast standards for neutron tomography. Physics Procedia, This issue, 2012. [6] P. Soille. Morphological image analysis. Springer Verlag, 2nd edition, 2002.  A.P. Kaestner et al. / Physics Procedia 43 (2013) 128 – 137 137 Appendix A. N point circle parameter estimate (x1 − c x )2 + (y1 − cy )2 = R2 (x2 − c x )2 + (y2 − cy )2 = R2 .. (A.1) . (xN − c x )2 + (yN − cy )2 = R2 Subtract the last equation from the other N − 1 eliminates R (x1 − c x )2 + (y1 − cy )2 − (xN − c x )2 − (yN − cy )2 = 0 (x2 − c x )2 + (y2 − cy )2 − (xN − c x )2 − (yN − cy )2 = 0 .. (A.2) . (xN−1 − c x )2 + (yN−1 − cy )2 − (xN − c x )2 − (yN − cy )2 = 0 Expand x12 − 2x1 c x + c2x + y21 − 2y1 cy + c2y − (x2N − 2xN c x + c2x ) − (y2N − 2yN cy + c2y ) = 0 x22 − 2x2 c x + c2x + y22 − 2y2 cy + c2y − (x2N − 2xN c x + c2x ) − (y2N − 2yN cy + c2y ) = 0 .. (A.3) . x2N−1 − 2xN−1 c x + c2x + y2N−1 − 2yN−1 cy + c2y − (x2N − 2xN c x + c2x ) − (y2N − 2yN cy + c2y ) = 0 Clean up x12 − 2x1 c x + y21 − 2y1 cy − (x2N − 2xN c x ) − (y2N − 2yN cy ) = 0 x22 − 2x2 c x + y22 − 2y2 cy − (x2N − 2xN c x ) − (y2N − 2yN cy ) = 0 .. (A.4) . x2N−1 − 2xN−1 c x + y2N−1 − 2yN−1 cy − (x2N − 2xN c x ) − (y23 − 2y3 cy ) = 0 Simplify x12 − x32 + y21 − y23 + 2(x3 − x1 )c x + 2(y3 − y1 )cy = 0 x22 − x32 + y22 − y23 + 2(x3 − x2 )c x + 2(y3 − y2 )cy = 0 .. (A.5) . x2N−1 − x2N + y2N−1 − y2N + 2(xN − xN−1 )c x + 2(yN − yN−1 )cy = 0 Put in matrix form ⎡ ⎤ ⎡ ⎤ ⎢⎢⎢ 2(xN − x1 ) 2(yN − y1 ) ⎥⎥ ⎥⎥⎥  ⎢⎢⎢ x2N − x12 + y2N − y21 ⎥⎥⎥ ⎢⎢⎢ 2(x − x ) 2(yN − y2 ) ⎥⎥ c x  ⎢ ⎢⎢⎢ x2N − x22 + y2N − y22 ⎥⎥⎥ ⎢⎢⎢ N 2 ⎥ ⎥ ⎢ ⎢ ⎥⎥⎥ ⎢⎢⎢ .. .. ⎥⎥⎥ = ⎢⎢⎢ .. ⎥⎥⎥ (A.6) ⎢⎢⎢ . . ⎥ c ⎥⎥⎦  ⎢⎢⎣ y ⎢ . ⎥⎥⎥ ⎣ ⎦ 2(xN − xN−1 ) 2(yN − yN−1 ) θ̂ xN − xN−1 + yN − yN−1 2 2 2 2   H a Solve Hθ̂ = a for θ̂  −1 θ̂ = HT H HT a (A.7) Provides the least square estimate of the center of the circle. The radius can be computed by inserting θ in any of the equations.