High-Index Glass Ball Retroreflectors for Measuring Lateral Positions

This paper is concerned with backscattered luminous signals, coming from a particular class of dielectric spheres illuminated by a coherent source. The purpose is to measure the lateral position of the sphere serving as an optical target, to achieve an overall contactless sensor of lateral position in space. Traditional approaches and theories such as ray-tracing and Mie scattering—as implemented in dedicated software—are applied to investigate their fitness for purpose in this application. No previous literature was found dealing with this specific case. Unfortunately, our observations did not match the theories’ predictions to an acceptable degree, and these approaches proved to be unsatisfactory. The rest of the paper focusses then on the development and comparison of suitable algorithms to compute the image coordinates of a representative point, which was in fact the true motivation of this work. Two original algorithms are proposed and discussed. Their robustness and repeatability are benchmarked under noisy conditions and at different distances from the target, with simulated as well as real images. Both resulted capable of sub-pixel accuracy.


Introduction
High-index ball lenses are widely used in optics-oriented applications, e.g., in coupling fibers and-more generally-in focusing or collimating light beams. This paper is centred on their peculiar property of retroreflecting the light. Several alternative retroreflector types ara available based on different principles, such as glass or hollow corner cubes [1,2], cat's eyes, and Spherically Mounted Retroreflectors (SMRs). One of the main disadvantages of these retroreflectors is a rather limited acceptance angle, restricing the illuminating beam to impinge from a narrow solid angle. To overcome this, glass spheres were introduced in the 1960s [3]. More recently, glass spheres with refracrive index n = 2 and thin supporting stems were introduced [4], to enlarge the acceptance angle to a nearly full solid angle (apart from the stem) and to exploit their (nominal) property of ideal retroreflection of light (at least under certain circumstances, see the next section for more details). However, they suffer poor efficiency as most of the impinging light is refracted through and is wasted on the rear: as little as approximately 15% of the light power is retroreflected in fact. Such spheres are available on the market up to a limited size (typically 8 to 10 mm diameter) as they are usually machined out of a single glass sheet. The introduction of these spheres triggered many studies available in literature, among which a detailed ray-tracing interpretation of their retroreflection [5].
Retroreflectors can be used for two different purposes: to measure either the longitudinal distance or the later position of a target to a reference point. Interferometric measurements are a typical example of the former use, while laser trackers [6,7] combine the two, effectively measuring the 3D spherical coordinates of a target point. The mentioned studies focussed on the longitudinal distance [4], while none investigated in full the ultimate limit of the lateral position detection. In addition, the actual distribution of the retroreflected light-detected as an image by a camera-has not been investigated yet.
The work described in this paper was intended to improve the performance of a novel coordinate measurement concept referred to as InPlanT [8]. In it, a n = 2 glass sphere is used as retroreflecting target, whose image is acquired by a camera. The sphere lateral position is detected and measured by the 2D (local) coordinates of the image in the camera system of coordinates. The algorithmical reduction of the image to a representative point is key and any improvement in the algorithm would impact straight to the overall InPlanT measurement accuracy. An apparently promising approach was by fitting the actual image against a template image as predicted by a model equation, which was not available though. This motivated a theoretical/experimental study to benchmark available theories in literature, from classical ray-tracing to the exotic theory of the glory. Unfortunately, none proved adequate for this specific case. As the original intent was to improve a localisation algorithm, we eventually forewent the image template approach and focussed on the only available first-principles characteristic, namely the intrinsic image central symmetry naturally stemming from a sphere.
In the end, this paper describes a method for measuring the lateral position of a n = 2 sphere target. Such measurement originated from a specific project (where the x-y coordinates of the center of the sphere needed to be carefully estimated to close the control loop of a tracking system), but it is of wider applicability, i.e., in non contact measurements in a view plane. This is a common problem to all camera based sensors. The knowledge of the actual beam pattern retroreflected by the target may be beneficial also when the light detection is done by means other than cameras, as it helps designing the optical set up.
The next section decribes the optical set up used in our investigation. Section 3 reviews exisiting theories in the attempt of predicting the actual image. Section 4 describes the developed algorithms. Some conclusions follow in Section 5.

Experimental Set Up
We used 16 mm diameter spheres made of Ohara S-LAH79 n = 2 glass (Ohara Corporation, Rancho Santa Margarita, CA, USA), illuminated by collimated coherent light; the spheres were kindly provided by the National Physical Laboratory (NPL, Teddington, Middlesex, UK); the complete optical characteristics of their glass can be found in the Ohara reference page, while the refractive index was verified on the site refractiveindex.info). The setup used to study the physics of the n = 2 spheres is shown in Figure 1. A 20 mW fiber coupled red laser crosses a polarizer and is collimated by means of a converging lens, then directed towards the sphere; the backscattered beam is then deviated at 90 • through a non-polarizing cubic beam splitter, and impinges on the active area of a CMOS camera; the acquired images, afterward, are post-processed to extract the information needed. In Figure 2 a typical pattern of the backscattered signal is shown: it is generated by a 8 mm radius sphere located at a distance of 254 mm far from the camera used to collect the images; our main purpose was to define a viable and repeatable method to gather information on the lateral positions of the sphere itself.
Also keeping the kind of glass fixed, the ring-shaped pattern revealed itself to be strongly dependent on some geometrical and physical factors, accounting for the sphere diameter, the reciprocal distance between sphere and camera, the accuracy of the collimation of the impinging beam, the portion of the sphere that is illuminated, the shape of the wavefront etc. More specifically, the thickness and the angular spacing of the rings is subject to dramatic variation when the sphere is moved along the direction joining the source of the signal to the center of the sphere; correlating such morphological parameters to the position of the sphere in space is the key factor for developing a metrological approach when using this kind of retroreflector in applications where the lateral displacement of a target is the measurand.

Review and Investigation of Existing Models
Some different strategies have been adopted to interpret the signals and infer the needed information: a classical ray-tracing approach, the more general Mie scattering theory and the models developed upon the Van de Hulst dissertation on light scattering by particles; anyway, all of them turned out to be unsatisfactory in explaining what happens in our particular situation, so a phenomenological and quantitative study of the peculiar pattern of the signals through two newly developed algorithms has been carried out, with the aim of extracting the localization of the target point.

Ray-Tracing
Most of the calculations are based on a few but elegant geometric optics concepts easily found in literature [5]. In the following Figures 3-6, the main optical properties of the spheres can be inferred, basing the calculations on the same principles. In Table 1, a list of the fundamental equations accounting for a geometrical interpretation of the backscattered signal is provided (same nomenclature as shown in Figure 3):

Review and Investigation of Existing Models
Some different strategies have been adopted to interpret the signals and infer the needed information: a classical ray-tracing approach, the more general Mie scattering theory and the models developed upon the Van de Hulst dissertation on light scattering by particles; anyway, all of them turned out to be unsatisfactory in explaining what happens in our particular situation, so a phenomenological and quantitative study of the peculiar pattern of the signals through two newly developed algorithms has been carried out, with the aim of extracting the localization of the target point.

Ray-Tracing
Most of the calculations are based on a few but elegant geometric optics concepts easily found in literature [5]. In the following Figures 3-6, the main optical properties of the spheres can be inferred, basing the calculations on the same principles. In Table 1, a list of the fundamental equations accounting for a geometrical interpretation of the backscattered signal is provided (same nomenclature as shown in Figure 3):

Review and Investigation of Existing Models
Some different strategies have been adopted to interpret the signals and infer the needed information: a classical ray-tracing approach, the more general Mie scattering theory and the models developed upon the Van de Hulst dissertation on light scattering by particles; anyway, all of them turned out to be unsatisfactory in explaining what happens in our particular situation, so a phenomenological and quantitative study of the peculiar pattern of the signals through two newly developed algorithms has been carried out, with the aim of extracting the localization of the target point.

Ray-Tracing
Most of the calculations are based on a few but elegant geometric optics concepts easily found in literature [5]. In the following Figures 3-6, the main optical properties of the spheres can be inferred, basing the calculations on the same principles. In Table 1, a list of the fundamental equations accounting for a geometrical interpretation of the backscattered signal is provided (same nomenclature as shown in Figure 3):    [5]).
Relation between ray (coming from the right side of the sphere) impact parameter b and impact angle φ (with R = radius of the sphere):

= sin
Deflection angle of the ray exiting the sphere ( Figure 4): Position r of the point (on a hypothetical screen at distance z from the sphere bigger than its radius) mapping the input ray: χ, ≅ sin + χ + sinχ Maximum deflection angle corresponding to the impact angle φm (when z ≫R): Condition of perfect (when χ = 0) back-reflection: Beam intensity distribution on a screen at distance z from the sphere:  . Transfer function for incident and reflected rays of different refractive indices n (the red curve is the most representative for our situation); as soon as a ray moves away from Z axis, increasing its impact parameter b, the deflection angle increases in this way: from a null value, it reaches a maximum in correspondence of , then it starts decreasing until a null value (parallel to itself but mapped in the other half-space) when = . Image inspired by [5] (Figure 3).

Ray-Optics on the Sphere at our Disposal
In our specific situation (glass used: Ohara S-LAH79, with n = 1.9957745 @ 635 nm, T = 20 °C and P = 1 Atm), numerical analyses with realistic refraction data have been carried on. Adopting the aforementioned nomenclature, the following relations can be derived:
From Figure 5b it can be argued that, according to the ray-optics, three distinct regions in the codomain of the mapping function r(b, z) exist: the first ("3-rays region"), mapped by the impact parameter range [-bm, bm], in correspondence of which, at a certain ray mapped on the screen, three different rays (with three distinct impact parameters) can impinge on the sphere giving the same r(χ, z); the other two regions, symmetrical to the central field, where there is a 1:1 correspondence between a ray impacting on the sphere and a scattered ray mapped on the screen. In Figure 6 a 3D . Transfer function for incident and reflected rays of different refractive indices n (the red curve is the most representative for our situation); as soon as a ray moves away from Z axis, increasing its impact parameter b, the deflection angle increases in this way: from a null value, it reaches a maximum in correspondence of φ m , then it starts decreasing until a null value (parallel to itself but mapped in the other half-space) when φ = φ 0 . Image inspired by [5] (Figure 3). codomain of the mapping function r(b, z) exist: the first ("3-rays region"), mapped by the impact parameter range [-bm, bm], in correspondence of which, at a certain ray mapped on the screen, three different rays (with three distinct impact parameters) can impinge on the sphere giving the same r(χ, z); the other two regions, symmetrical to the central field, where there is a 1:1 correspondence between a ray impacting on the sphere and a scattered ray mapped on the screen. In Figure 6 a 3D representation of the same concept (changing χ with b) can be seen.   Table 2 shows the geometrical properties of the ray-tracing simulation, in input and in output. Figure 7b justifies the choice of S-LAH79 glass for our purposes: the plot clearly shows that the deflection of a beam impacting on a sphere made with this glass changes quite slowly with respect to the impact parameter, inside a reasonable range (low ray deflection zone: | | < 10 mrad for b < 0.35•R), without requiring an unpractically low b/R ratio, whereas the other glasses suffer from a greatest angular spread of the exiting beam, at the same radius of the sphere.   Table 1. Fundamental equations for tracing rays of light impacting on a sphere (taken from [5]).
Relation between ray (coming from the right side of the sphere) impact parameter b and impact angle φ (with R = radius of the sphere): Deflection angle of the ray exiting the sphere ( Figure 4): Position r of the point (on a hypothetical screen at distance z from the sphere bigger than its radius) mapping the input ray: Maximum deflection angle corresponding to the impact angle ϕ m (when z R): Condition of perfect (when χ = 0) back-reflection: Beam intensity distribution on a screen at distance z from the sphere: Needed aperture radius of the impinging beam to have the maximum acceptance of the sphere: Ray-Optics on the Sphere at our Disposal In our specific situation (glass used: Ohara S-LAH79, with n = 1.9957745 @ 635 nm, T = 20 • C and P = 1 Atm), numerical analyses with realistic refraction data have been carried on. Adopting the aforementioned nomenclature, the following relations can be derived:
From Figure 5b it can be argued that, according to the ray-optics, three distinct regions in the codomain of the mapping function r(b, z) exist: the first ("3-rays region"), mapped by the impact parameter range [-b m , b m ], in correspondence of which, at a certain ray mapped on the screen, three different rays (with three distinct impact parameters) can impinge on the sphere giving the same r(χ, z); the other two regions, symmetrical to the central field, where there is a 1:1 correspondence between a ray impacting on the sphere and a scattered ray mapped on the screen. In Figure 6 a 3D representation of the same concept (changing χ with b) can be seen. Table 2 shows the geometrical properties of the ray-tracing simulation, in input and in output. Figure 7b justifies the choice of S-LAH79 glass for our purposes: the plot clearly shows that the deflection of a beam impacting on a sphere made with this glass changes quite slowly with respect to the impact parameter, inside a reasonable range (low ray deflection zone: |χ| < 10 mrad for b < 0.35·R), without requiring an unpractically low b/R ratio, whereas the other glasses suffer from a greatest angular spread of the exiting beam, at the same radius of the sphere.
Scattering angles at which the intensity is infinite, according to geometrical optics [9]; impact angles φ determininig rainbows, after undergoing n r reflections inside the sphere:

Mie Scattering: A Simple Software Approach
On the web, many software programs aimed to solving the problem of scattering of electromagnetic signals by spheres/spheroids can be found [10,11] and tested; among them, we decided to give preference to the package PyMieScatt [12], implemented in Python, and to the software MiePlot [13]. The angular functions of interest were, respectively, the parallel, perpendicular, and unpolarized scattered intensities SR(θ), SL(θ) and SU(θ), with a focus on the backscattering region (θ~180°), that is the region where the phenomenon of the glory makes its spectacular appearance (Figures 8 and 9). The literature about glories is quite rich [14][15][16], but the most recent outcomes on the topic derive from postulating the existence of surface waves travelling around the outside of the sphere [17,18]. In our case, we deal with a sphere with a refractive index very close to 2, and this means that backscattering is something correlated by geometric rays passing through a huge (with respect to the impinging wavelength) sphere and then suffering at least one internal reflection before travelling back towards the light source; so, simply speaking, the empirical approach does not require us to threat the phenomenon as a glory in a strict sense. In this paper, indeed, we are not interested in elaborating a mathematical model capable of explaining the observable scattering patterns produced by retro-reflecting spheres, but rather to develop a viable artifice to interpret our peculiar class of patterns-and consequently infer metrological information about the coordinates of a target-that, to our knowledge, have not been described before in literature.

Mie Scattering: A Simple Software Approach
On the web, many software programs aimed to solving the problem of scattering of electromagnetic signals by spheres/spheroids can be found [10,11] and tested; among them, we decided to give preference to the package PyMieScatt [12], implemented in Python, and to the software MiePlot [13]. The angular functions of interest were, respectively, the parallel, perpendicular, and unpolarized scattered intensities SR(θ), SL(θ) and SU(θ), with a focus on the backscattering region (θ~180 • ), that is the region where the phenomenon of the glory makes its spectacular appearance (Figures 8 and 9). The literature about glories is quite rich [14][15][16], but the most recent outcomes on the topic derive from postulating the existence of surface waves travelling around the outside of the sphere [17,18]. In our case, we deal with a sphere with a refractive index very close to 2, and this means that backscattering is something correlated by geometric rays passing through a huge (with respect to the impinging wavelength) sphere and then suffering at least one internal reflection before travelling back towards the light source; so, simply speaking, the empirical approach does not require us to threat the phenomenon as a glory in a strict sense. In this paper, indeed, we are not interested in elaborating a mathematical model capable of explaining the observable scattering patterns produced by retro-reflecting spheres, but rather to develop a viable artifice to interpret our peculiar class of patterns-and consequently infer metrological information about the coordinates of a target-that, to our knowledge, have not been described before in literature.

Mie Scattering: A Simple Software Approach
On the web, many software programs aimed to solving the problem of scattering of electromagnetic signals by spheres/spheroids can be found [10,11] and tested; among them, we decided to give preference to the package PyMieScatt [12], implemented in Python, and to the software MiePlot [13]. The angular functions of interest were, respectively, the parallel, perpendicular, and unpolarized scattered intensities SR(θ), SL(θ) and SU(θ), with a focus on the backscattering region (θ~180°), that is the region where the phenomenon of the glory makes its spectacular appearance (Figures 8 and 9). The literature about glories is quite rich [14][15][16], but the most recent outcomes on the topic derive from postulating the existence of surface waves travelling around the outside of the sphere [17,18]. In our case, we deal with a sphere with a refractive index very close to 2, and this means that backscattering is something correlated by geometric rays passing through a huge (with respect to the impinging wavelength) sphere and then suffering at least one internal reflection before travelling back towards the light source; so, simply speaking, the empirical approach does not require us to threat the phenomenon as a glory in a strict sense. In this paper, indeed, we are not interested in elaborating a mathematical model capable of explaining the observable scattering patterns produced by retro-reflecting spheres, but rather to develop a viable artifice to interpret our peculiar class of patterns-and consequently infer metrological information about the coordinates of a target-that, to our knowledge, have not been described before in literature.

Novel Algorithms for Center Detection
By using both of the computational approaches described in the previous paragraph, we were not able to find a useful correspondence with our intensity patterns: on the one hand, ray-tracing concepts help in predicting the existence of an "accumulation" zone-only one ring where we have many of them-for the rays scattered by the sphere very close to the 180° region (the primary rainbow, as stated in Table 2, corresponding to nr = 1); on the other hand, the Mie approach has shown to be affected by two issues: 1) The simulated backscattering pattern (Figure 8b) clearly shows a monotonically decreasing envelope (with a hyperbolic trend, according to [17] (ch.13)) of the ripples, starting from the 180° zone, whereas we can infer the presence of a sort of modulating function that breaks the regularity of the envelope (see plot in Figure

Novel Algorithms for Center Detection
By using both of the computational approaches described in the previous paragraph, we were not able to find a useful correspondence with our intensity patterns: on the one hand, ray-tracing concepts help in predicting the existence of an "accumulation" zone-only one ring where we have many of them-for the rays scattered by the sphere very close to the 180 • region (the primary rainbow, as stated in Table 2, corresponding to n r = 1); on the other hand, the Mie approach has shown to be affected by two issues: (1) The simulated backscattering pattern (Figure 8b) clearly shows a monotonically decreasing envelope (with a hyperbolic trend, according to [17] (ch.13)) of the ripples, starting from the 180 • zone, whereas we can infer the presence of a sort of modulating function that breaks the regularity of the envelope (see plot in Figure   For these reasons, we decided to adopt a strictly computational strategy to infer the "centers" of the acquired images based on the only first principle property of the image we could rely on, i.e., the central symmetry of a spherical image. To do so, the cartesian space is convertedto a polar space. The transformed image is then evaluated according to two different criteria implemented in Python, as described in the following sections. Both proved to achieve sub-pixel accuracy.

Variance-Based Algorithm
This version of software is built on the concept of the minimum value of the sum of the variances per column of the pixel intensities in the polar-converted input images of the backscattered signal (block diagram shown in Figure 10). transformed image is then evaluated according to two different criteria implemented in Python, as described in the following sections. Both proved to achieve sub-pixel accuracy.

Variance-Based Algorithm
This version of software is built on the concept of the minimum value of the sum of the variances per column of the pixel intensities in the polar-converted input images of the backscattered signal (block diagram shown in Figure 10).
The algorithm branches off into two different methods, labeled as "Automatic" or "Forced" procedures. In the former modality, through a subroutine that helps in recognizing the central disc of the acquired image (1A), it starts finding the coordinates (y_start, x_start) of the point of maximum brightness (2A), that locates the region where the maximum of the scattered energy is stored on the detector; some explicative and simplified lines of code will follow: Figure 10. Block diagram of the variance-based approach to the problem. The next step (3A) is determining a squared region on interest ("ROI", whose odd side length is equal to "span"), centered on the brightest pixel of the image, that will be the "envelope" of the points defining the centers (subsequently cycled in a double loop) of what we have defined "sensitive areas" (4), exploring the whole image: The algorithm branches off into two different methods, labeled as "Automatic" or "Forced" procedures. In the former modality, through a subroutine that helps in recognizing the central disc of the acquired image (1A), it starts finding the coordinates (y _ start, x _ start) of the point of maximum brightness (2A), that locates the region where the maximum of the scattered energy is stored on the detector; some explicative and simplified lines of code will follow: img64 _ float = source.astype(np.float64) a = np.array(img64 _ float) data = np.unravel _ index(a.argmax(), a.shape, order='C') The next step (3A) is determining a squared region on interest ("ROI", whose odd side length is equal to "span"), centered on the brightest pixel of the image, that will be the "envelope" of the points defining the centers (subsequently cycled in a double loop) of what we have defined "sensitive areas" (4), exploring the whole image: span = int(np.ceil(kp.size/(4 * np.sqrt(2))) * 2 + 1) In the automatic procedure, kp.size is the diameter of the blob found by the OpenCV [19] SimpleBlobDetector class: in this formulation, span is equal to half of the side of the maximum square inscribed in the central bright disc of the image; in the forced procedure (1B), this parameter is simply replaced by an arbitrary quantity, evalutaed in relation with the central disc size (this sometimes reduced the computation time during the testing process). An appropriate cropping of the original image, then, avoids inconsistencies when dealing with situations when the pattern of the signal is very asymmetrical inside the frame: The iteration order n (with 1 < n < span 2 ) where the minimum of the sum of the variances per columns is found, finally, identifies the pole of the polar transformation corresponding to the best candidate as "center" of the original image (7). The sub-pixel accuracy is obtained by least squares fitting the 3D matrix of the variances inside the ROI area with a rotated paraboloid (Figure 11a): matrix _ variances = np.transpose(np.resize(results, int(span) ** 2).reshape(int(span), int(span))) def parab(X, amp, shift _ x, shift _ y, cost): x, y =~X return amp * (x ** 2 + y ** 2) + shift _ x * x + shift _ y * y + cost shift _ x = popt [1] shift _ y = popt [2] cost = popt [3] Sensors 2019, 19

FFT-Based Algorithm
This second version of the software is built on the concept of the minimum value of the sum per column of the pixel intensities of the FFT-transformed pictures of the polar-converted original images (block diagram shown in Figure 12).
In this more sophisticated arrangement, the strategy followed by the FFT-based algorithm is the same as before in its first steps (1-5), but the basic concept, here, is computing the 2D discrete Fourier Transform of a bidimensional array (our polar-transformed image) by means of the Fast Fourier Transform (Numpy function fft.fft2) (

FFT-Based Algorithm
This second version of the software is built on the concept of the minimum value of the sum per column of the pixel intensities of the FFT-transformed pictures of the polar-converted original images (block diagram shown in Figure 12).
In this more sophisticated arrangement, the strategy followed by the FFT-based algorithm is the same as before in its first steps (1-5), but the basic concept, here, is computing the 2D discrete Fourier Transform of a bidimensional array (our polar-transformed image) by means of the Fast Fourier Transform (Numpy function fft.fft2) (6): for x, y in itertools.product(x, y): The hearth of this portion of the algorithm is the removal of the central row of the raw power spectrum, that is the DC component of the FFT: choosing a "bad" pole, in the double loop, for the cartesian to polar transformation, in fact, induces a periodic perturbation of the vertical bands of the polar image (adding a spurious component to the spectrum of the frequencies of the image), as much important as the further we go away from the point of maximum symmetry. The choice of the point that minimizes the sums of the squared intensities per column-the probe of this version of the algorithm-of the raw power spectrum (deprived of the central row) of the polar-transformed image defines, therefore, the final coordinates of our "center" (Figure 13). The sub-pixel accuracy is then gained like in the previous version of the algorithm.

Analysis of the Pattern of the scattered Signal
After determining the centers of the acquired images (Figure 14), we studied the radial profile of the backscattered signal. It showed luminous and dark rings, with a central bright disc concentrating about 60% of the luminous power ( Figure 15).

Analysis of the Pattern of the scattered Signal
After determining the centers of the acquired images (Figure 14), we studied the radial profile of the backscattered signal. It showed luminous and dark rings, with a central bright disc concentrating about 60% of the luminous power ( Figure 15).

Analysis of the Pattern of the scattered Signal
After determining the centers of the acquired images (Figure 14), we studied the radial profile of the backscattered signal. It showed luminous and dark rings, with a central bright disc concentrating about 60% of the luminous power ( Figure 15).  We then applied the developed two algorithms to several images, computed the centers of the patterns, and inferred the radial distances and angular separation (as seen from the sphere) of the dark rings ( Figure 16). We observed that:

•
At each distance, a slow envelope function modulates a higher frequency oscillation ( Figure 15). This latter turns out not fully periodical, as the ring mutual distances decrease outwards.

•
The pattern increases in size with increasing distances, revealing a divergence of the backscattered beam, at least in the far field.
• This divergence does not exhibit the same angle (for a specific ring) at all distances. Figure 16b shows an increase of the dark ring angular separation with the distance. If each dark ring is regarded as the base of a cone with apex in the target sphere, then the apex angle increases with the sphere distance, as in a sort of angular amplifier. It would remain constant instead if the backscattered beam were shaped as a series of coaxial cones. This is in contrast with one of the most accepted models for the glory found in literature [17]. This theory asserts that the backscattered wavefront has a toroidal shape whose intensity is a function (linear combination of squared Bessel functions of first kind) of a parameter = π , where a is the radius of the sphere, n the sphere refractive index and is supplementary to the backscattering angle ( Figure  17). There is no dependance on the distance, so the divergency angle u is expected to be the same at all distances with no near/far-field separation.  We then applied the developed two algorithms to several images, computed the centers of the patterns, and inferred the radial distances and angular separation (as seen from the sphere) of the dark rings ( Figure 16). We observed that:

•
At each distance, a slow envelope function modulates a higher frequency oscillation ( Figure 15). This latter turns out not fully periodical, as the ring mutual distances decrease outwards.

•
The pattern increases in size with increasing distances, revealing a divergence of the backscattered beam, at least in the far field.
• This divergence does not exhibit the same angle (for a specific ring) at all distances. Figure 16b shows an increase of the dark ring angular separation with the distance. If each dark ring is regarded as the base of a cone with apex in the target sphere, then the apex angle increases with the sphere distance, as in a sort of angular amplifier. It would remain constant instead if the backscattered beam were shaped as a series of coaxial cones. This is in contrast with one of the most accepted models for the glory found in literature [17]. This theory asserts that the backscattered wavefront has a toroidal shape whose intensity is a function (linear combination of squared Bessel functions of first kind) of a parameter = π , where a is the radius of the sphere, n the sphere refractive index and is supplementary to the backscattering angle ( Figure  17). There is no dependance on the distance, so the divergency angle u is expected to be the same at all distances with no near/far-field separation. We then applied the developed two algorithms to several images, computed the centers of the patterns, and inferred the radial distances and angular separation (as seen from the sphere) of the dark rings ( Figure 16). We observed that:

•
At each distance, a slow envelope function modulates a higher frequency oscillation ( Figure 15). This latter turns out not fully periodical, as the ring mutual distances decrease outwards.

•
The pattern increases in size with increasing distances, revealing a divergence of the backscattered beam, at least in the far field. • This divergence does not exhibit the same angle (for a specific ring) at all distances. Figure 16b shows an increase of the dark ring angular separation with the distance. If each dark ring is regarded as the base of a cone with apex in the target sphere, then the apex angle increases with the sphere distance, as in a sort of angular amplifier. It would remain constant instead if the backscattered beam were shaped as a series of coaxial cones. This is in contrast with one of the most accepted models for the glory found in literature [17]. This theory asserts that the backscattered wavefront has a toroidal shape whose intensity is a function (linear combination of squared Bessel functions of first kind) of a parameter u = 2πan λ γ, where a is the radius of the sphere, n the sphere refractive index and γ. is supplementary to the backscattering angle ( Figure 17). There is no dependance on the distance, so the divergency angle u is expected to be the same at all distances with no near/far-field separation.

Comparison with "Classical" Software and Performance Analysis
We compared the developed algorithms with estabilished software (Classes cv2.moments and scipy.ndimage.center_of_mass and ImageJ with its measurement tool.). The derived image center after segmentation was the base for comparison. We carried out some tests: 1) Center repeatability vs. distance of the camera to the sphere. This test is based on simulation only. We used Zemax [20] as simulation software. Figure 18 shows the used layouts, Figure 19 some examples of the simulated backscattered images, Figure 20 the size of a surrounding ring as a function of the distance to the sphere. Figure 20 shows a minimum, which separates the converging (near-field) from the diverging (far-field) portions of the beam, thus being a "waist". The derived centers were compared with the theoretical one imposed by simulation. The discrepancy ∆ = δ + δ is a measure of the error incurred by each algorithm at each distance ( Figure 21).

Comparison with "Classical" Software and Performance Analysis
We compared the developed algorithms with estabilished software (Classes cv2.moments and scipy.ndimage.center_of_mass and ImageJ with its measurement tool.). The derived image center after segmentation was the base for comparison. We carried out some tests: 1) Center repeatability vs. distance of the camera to the sphere. This test is based on simulation only. We used Zemax [20] as simulation software. Figure 18 shows the used layouts, Figure 19 some examples of the simulated backscattered images, Figure 20 the size of a surrounding ring as a function of the distance to the sphere. Figure 20 shows a minimum, which separates the converging (near-field) from the diverging (far-field) portions of the beam, thus being a "waist". The derived centers were compared with the theoretical one imposed by simulation. The discrepancy ∆ = δ + δ is a measure of the error incurred by each algorithm at each distance ( Figure 21). Figure 17. Intensity of the glory as predicted in [17] (with C1 = 0, C2 = 1).

Comparison with "Classical" Software and Performance Analysis
We compared the developed algorithms with estabilished software (Classes cv2.moments and scipy.ndimage.center_of_mass and ImageJ with its measurement tool.). The derived image center after segmentation was the base for comparison. We carried out some tests: (1) Center repeatability vs. distance of the camera to the sphere. This test is based on simulation only.
We used Zemax [20] as simulation software. Figure 18 shows the used layouts, Figure 19 some examples of the simulated backscattered images, Figure 20 the size of a surrounding ring as a function of the distance to the sphere. Figure 20 shows a minimum, which separates the converging (near-field) from the diverging (far-field) portions of the beam, thus being a "waist". The derived centers were compared with the theoretical one imposed by simulation. The discrepancy ∆ = δx 2 + δy 2 is a measure of the error incurred by each algorithm at each distance ( Figure 21). The simulation allowed to move the target over an ideal straight line parallel to the camera optical axis, theoretically resulting in the same image coordinates of the centers. This would have been impossible in practice-or very difficult and expensive-over distances of several meters. The pooled standard deviations of the derived center x and y coordinates over the full distance range (100-5000) mm were: 0.52 and 0.51 pixel for variance-based and FFT-based algorithms respectively; 0.69 pixel for the scipy.ndimage.center_of_mass algorithm; and an impressive 0.03 pixel for the ImageJ software. In ideal conditions (full pattern captured with no border cropping and no noise), ImageJ outstands.      The simulation allowed to move the target over an ideal straight line parallel to the camera optical axis, theoretically resulting in the same image coordinates of the centers. This would have been impossible in practice-or very difficult and expensive-over distances of several meters. The pooled standard deviations of the derived center x and y coordinates over the full distance range (100-5000) mm were: 0.52 and 0.51 pixel for variance-based and FFT-based algorithms respectively; 0.69 pixel for the scipy.ndimage.center_of_mass algorithm; and an impressive 0.03 pixel for the ImageJ software. In ideal conditions (full pattern captured with no border cropping and no noise), ImageJ outstands. 2) Center accuracy vs. lateral displacement of the sphere. We then investigated the effect of eccentriciy of the image to the camera and consequent border cropping of the signal; in other words, when the sensor is not used merely for zero detection but for a full measurement of the lateral position. This might be the case e.g., when alignment or tracking is being sought. A same acquired (not simulated) image was shifted in steps to obtain other images corresponding to as many shifts of the sphere. In this case, the developed algorithms are the best choice for achieving the information needed: they are reliable and robust, since they can correctly track the center also when the symmetry of the signal is seriously impaired or totally broken (e.g.,: when an important part of the peripheral rings is cut off the frame); in these situations, the coordinates of the center determined by ImageJ are manifestly wrong. To show this concept, a set of 11 images were cropped along a specific direction (x-axis), with each new image separated from the previous one by an arbitrary number of pixel; the error in the determination of the x-coordinate of the center has then been calculated and plotted for comparison. In the total absence of a thresholding strategy, when applying ImageJ measurement tool to the cropped (rough) image sequence, the accuracy of VAR and FFT-based algorithms outperforms ImageJ capabilities (Figure 22a). When the correct ImageJ thresholding procedure selects the portion of the signal  The simulation allowed to move the target over an ideal straight line parallel to the camera optical axis, theoretically resulting in the same image coordinates of the centers. This would have been impossible in practice-or very difficult and expensive-over distances of several meters. The pooled standard deviations of the derived center x and y coordinates over the full distance range (100-5000) mm were: 0.52 and 0.51 pixel for variance-based and FFT-based algorithms respectively; 0.69 pixel for the scipy.ndimage.center_of_mass algorithm; and an impressive 0.03 pixel for the ImageJ software. In ideal conditions (full pattern captured with no border cropping and no noise), ImageJ outstands. 2) Center accuracy vs. lateral displacement of the sphere. We then investigated the effect of eccentriciy of the image to the camera and consequent border cropping of the signal; in other words, when the sensor is not used merely for zero detection but for a full measurement of the lateral position. This might be the case e.g., when alignment or tracking is being sought. A same acquired (not simulated) image was shifted in steps to obtain other images corresponding to as many shifts of the sphere. In this case, the developed algorithms are the best choice for achieving the information needed: they are reliable and robust, since they can correctly track the center also when the symmetry of the signal is seriously impaired or totally broken (e.g.,: when an important part of the peripheral rings is cut off the frame); in these situations, the coordinates of the center determined by ImageJ are manifestly wrong. To show this concept, a set of 11 images were cropped along a specific direction (x-axis), with each new image separated from the previous one by an arbitrary number of pixel; the error in the determination of the x-coordinate of the center has then been calculated and plotted for comparison. In the total absence of a thresholding strategy, when applying ImageJ measurement tool to the cropped (rough) image sequence, the accuracy of VAR and FFT-based algorithms outperforms ImageJ capabilities (Figure 22a). When the correct ImageJ thresholding procedure selects the portion of the signal (2) Center accuracy vs. lateral displacement of the sphere. We then investigated the effect of eccentriciy of the image to the camera and consequent border cropping of the signal; in other words, when the sensor is not used merely for zero detection but for a full measurement of the lateral position. This might be the case e.g., when alignment or tracking is being sought. A same acquired (not simulated) image was shifted in steps to obtain other images corresponding to as many shifts of the sphere. In this case, the developed algorithms are the best choice for achieving the information needed: they are reliable and robust, since they can correctly track the center also when the symmetry of the signal is seriously impaired or totally broken (e.g.,: when an important part of the peripheral rings is cut off the frame); in these situations, the coordinates of the center determined by ImageJ are manifestly wrong. To show this concept, a set of 11 images were cropped along a specific direction (x-axis), with each new image separated from the previous one by an arbitrary number of pixel; the error in the determination of the x-coordinate of the center has then been calculated and plotted for comparison. In the total absence of a thresholding strategy, when applying ImageJ measurement tool to the cropped (rough) image sequence, the accuracy of VAR and FFT-based algorithms outperforms ImageJ capabilities (Figure 22a). When the correct ImageJ thresholding procedure selects the portion of the signal that keeps its best symmetry and integrity during the lateral displacement of the sphere, the performances of the three methods become comparable ( Figure 22b); here, our algorithms still work on the original, not pre-processed images); anyway, also in this case, the variance of the plotted parameter obtained by means of ImageJ is about one order of magnitude greater than the variance of the same parameter obtained using our algorithms: this means that the latter generate coordinates of the centers that are stable and reliable also in presence of images showing patterns with a very low level of simmetry.
The artifice that can be used to improve the ImageJ results is, therefore, thresholding the signal with the most appropriate method (that can be different from time to time, according to the pattern, the noise, the number, the thickness and the spacing of the rings . . . ) to isolate its most symmetric portion inside the frame; however, the price, in terms of the amount of the power that is lost in this way, can be high, and the information acquired about the center inaccurate. Without such pre-processing of the eccentric images, on the other hand, ImageJ is very likely to generate macroscopic inconsistencies in the determination of the needed information ( Figure 23).
that keeps its best symmetry and integrity during the lateral displacement of the sphere, the performances of the three methods become comparable ( Figure 22b); here, our algorithms still work on the original, not pre-processed images); anyway, also in this case, the variance of the plotted parameter obtained by means of ImageJ is about one order of magnitude greater than the variance of the same parameter obtained using our algorithms: this means that the latter generate coordinates of the centers that are stable and reliable also in presence of images showing patterns with a very low level of simmetry. The artifice that can be used to improve the ImageJ results is, therefore, thresholding the signal with the most appropriate method (that can be different from time to time, according to the pattern, the noise, the number, the thickness and the spacing of the rings…) to isolate its most symmetric portion inside the frame; however, the price, in terms of the amount of the power that is lost in this way, can be high, and the information acquired about the center inaccurate. Without such pre-processing of the eccentric images, on the other hand, ImageJ is very likely to generate macroscopic inconsistencies in the determination of the needed information ( Figure 23).   that keeps its best symmetry and integrity during the lateral displacement of the sphere, the performances of the three methods become comparable ( Figure 22b); here, our algorithms still work on the original, not pre-processed images); anyway, also in this case, the variance of the plotted parameter obtained by means of ImageJ is about one order of magnitude greater than the variance of the same parameter obtained using our algorithms: this means that the latter generate coordinates of the centers that are stable and reliable also in presence of images showing patterns with a very low level of simmetry. The artifice that can be used to improve the ImageJ results is, therefore, thresholding the signal with the most appropriate method (that can be different from time to time, according to the pattern, the noise, the number, the thickness and the spacing of the rings…) to isolate its most symmetric portion inside the frame; however, the price, in terms of the amount of the power that is lost in this way, can be high, and the information acquired about the center inaccurate. Without such pre-processing of the eccentric images, on the other hand, ImageJ is very likely to generate macroscopic inconsistencies in the determination of the needed information ( Figure 23).   (3) Repeatability of center determination vs. noise of the image. In this case, the experimentally acquired images were artificially altered by means of Gaussian noise with a gradually increasing variance: a set of 10 images of the sphere at a fixed distance from the sphere, with decreasing peak signal to noise ratio was used in a test-like environment, starting from the image with the lowest level of noise (very close to the one effectively acquired by the CMOS camera) to the most noisy one. Figure 24 shows the performances of the developed algorithms in a short sphere-camera distance situation. The developed algorithms do not perform remarkably better than the existing ones for all the investigated distances: orange and grey bars-representative of the σ-pooled of the centers found by using VAR and FFT-based algorithms-are, apart from a couple of occurrences, higher than yellow and blue ones-representative of the competitor algorithms-σ = 40 is the maximum value of noise tolerated by them; values of σ > 40 deteriorate too much the images and generate inconsistencies in the determination of the coordinates of the center. Anyway, it is useful to point out that, in some situations, the centroid algorithm completely failed the test (this is the reason why it is not represented in the plot), while both VAR and FFT algorithms always succeeded in the task of determining the center, even in the cases of heavily degraded images.
a set of 10 images of the sphere at a fixed distance from the sphere, with decreasing peak signal to noise ratio was used in a test-like environment, starting from the image with the lowest level of noise (very close to the one effectively acquired by the CMOS camera) to the most noisy one. Figure 24 shows the performances of the developed algorithms in a short sphere-camera distance situation. The developed algorithms do not perform remarkably better than the existing ones for all the investigated distances: orange and grey bars-representative of the σ-pooled of the centers found by using VAR and FFT-based algorithms-are, apart from a couple of occurrences, higher than yellow and blue ones-representative of the competitor algorithms-σ = 40 is the maximum value of noise tolerated by them; values of σ > 40 deteriorate too much the images and generate inconsistencies in the determination of the coordinates of the center. Anyway, it is useful to point out that, in some situations, the centroid algorithm completely failed the test (this is the reason why it is not represented in the plot), while both VAR and FFT algorithms always succeeded in the task of determining the center, even in the cases of heavily degraded images. Another way of estimating the performance of the algorithms is to test the repeatability of center determination with respect to the choice of the ROI size: the algorithms turn out to be as much robust as they are more capable of generating "similar" center coordinates Cxi and Cyi, independently on the ROI size. To verify this aspect (directly correlated to the achievable uncertainty), the quantity  Another way of estimating the performance of the algorithms is to test the repeatability of center determination with respect to the choice of the ROI size: the algorithms turn out to be as much robust as they are more capable of generating "similar" center coordinates Cx i and Cy i , independently on the ROI size. To verify this aspect (directly correlated to the achievable uncertainty), the quantity has been plotted for a randomly chosen test image (Figure 25). The distribution of the results obtained seem to limit the best achievable uncertainty to a range of ~±0.5 pixel for the variance-based algorithm and ±0.3 pixel for the FFT-based algorithm. It is remarkable the fact that a small ROI size (compared to the diameter of the central disc, equal to 85 pixel in the image studied) seems adequate to ensure a small value of uncertainty: that allows a fast computation time to determine the center coordinates, even in the case of high-resolution images. The distribution of the results obtained seem to limit the best achievable uncertainty to a range of~±0.5 pixel for the variance-based algorithm and ±0.3 pixel for the FFT-based algorithm. It is remarkable the fact that a small ROI size (compared to the diameter of the central disc, equal to 85 pixel in the image studied) seems adequate to ensure a small value of uncertainty: that allows a fast computation time to determine the center coordinates, even in the case of high-resolution images.

Some Words on Longitudinal Metrology
As stated before, the pattern of the backscattered signal captured by the camera shows a strong dependence on the distance of the camera from the sphere: the ring spacing thickens when the camera is close to the sphere, while it thins out when the camera is far from the sphere. In first approximation, one could be tempted by the possibility of exploiting this behavior to infer some kind of longitudinal metrology on the effective distance of the sphere: let's see with which level of accuracy. Adding to the algorithms the capability of averaging-along columns-the pixel intensity of the polar-transformed image of the sphere, and then performing the FFT on the row vector of these values, the coordinates of the peaks reveal the average spacing of the bright rings of the original image. The subsequent step is fitting the positions of these peaks with respect to the distance of the sphere, as shown in Figure 26. The distribution of the results obtained seem to limit the best achievable uncertainty to a range of ~±0.5 pixel for the variance-based algorithm and ±0.3 pixel for the FFT-based algorithm. It is remarkable the fact that a small ROI size (compared to the diameter of the central disc, equal to 85 pixel in the image studied) seems adequate to ensure a small value of uncertainty: that allows a fast computation time to determine the center coordinates, even in the case of high-resolution images.

Some Words on Longitudinal Metrology
As stated before, the pattern of the backscattered signal captured by the camera shows a strong dependence on the distance of the camera from the sphere: the ring spacing thickens when the camera is close to the sphere, while it thins out when the camera is far from the sphere. In first approximation, one could be tempted by the possibility of exploiting this behavior to infer some kind of longitudinal metrology on the effective distance of the sphere: let's see with which level of accuracy. Adding to the algorithms the capability of averaging-along columns-the pixel intensity of the polartransformed image of the sphere, and then performing the FFT on the row vector of these values, the coordinates of the peaks reveal the average spacing of the bright rings of the original image. The subsequent step is fitting the positions of these peaks with respect to the distance of the sphere, as shown in Figure 26.  The most accurate available points (in a limited range of distances) have been used to fit the data, but, also in this way, the slightly increasing slope of the fitting line prevents from getting a very useful information: solving the fitting equation with respect to the variable d sphere-camera , the equation: is obtained, where A = intercept and B = slope of the line, the distance cannot be resolved at a high degree (e.g.,: submillimetric) of accuracy; indeed, the sensibility 1/B of the whole system is approximately equal to 6.25 mm/px, too low to confer a metrological valence to this approach of image analysis. Better results could be achieved in presence of a camera with a higher resolution and a bigger active area of the sensor, in order to collect a useful part of the diverging signal at greater distances.

Conclusions
With this work, the possibility of extracting accurate metrological information about the position of glass spheres, used as targets, with refractive index almost equal to 2 was investigated. More specifically, two completely new algorithms were written to generate the coordinates of the centers of the same spheres through the analysis of their backscattered signals: the complexity of these signals, combined to the lack of a satisfying theoretical approach to interpret those such peculiar patterns observed, forced us to limit our research in extracting quantitative data from empirical observations rather than applying ready-made models. The main strength that these algorithms demonstrated in the validation tests is their ability of working on images with a poor level of symmetry of the pattern caught on the active area of the sensor, a condition that experimentally is quite common, and inferring the desired information without a specific pre-processing (tipically starting from an arbitrary and dangerous thresholding) of the images themselves. Moreover, they proved to be robust and repeatable, also in conditions of heavy noise and degradated images: their accuracy is comparable to the classical software used during the validation tests.
Future developments of this work will concern partly the improvement of the codes efficiencyhopefully taking advantage of the parallel computation through optimized toolkits (like CUDA [21]), very useful when dealing with high resolution images-and on their integration into a real-time acquisition system, partly the efforts to build a suitable mathematical model in order to explain, in a rigorous way, the behavior of the specific retroreflectors we used for target acquisition.  [23]. This project has received funding from the EMPIR programme co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme. image analysis. Better results could be achieved in presence of a camera with a higher a bigger active area of the sensor, in order to collect a useful part of the diverging si distances.

Conclusions
With this work, the possibility of extracting accurate metrological information abo of glass spheres, used as targets, with refractive index almost equal to 2 was inve specifically, two completely new algorithms were written to generate the coordinate of the same spheres through the analysis of their backscattered signals: the comp signals, combined to the lack of a satisfying theoretical approach to interpret those patterns observed, forced us to limit our research in extracting quantitative data observations rather than applying ready-made models. The main strength that th demonstrated in the validation tests is their ability of working on images with a symmetry of the pattern caught on the active area of the sensor, a condition that ex quite common, and inferring the desired information without a specific pre-proce starting from an arbitrary and dangerous thresholding) of the images themselves. M proved to be robust and repeatable, also in conditions of heavy noise and degradate accuracy is comparable to the classical software used during the validation tests.
Future developments of this work will concern partly the improvement of the cod hopefully taking advantage of the parallel computation through optimized toolkits (li very useful when dealing with high resolution images-and on their integration i acquisition system, partly the efforts to build a suitable mathematical model in order rigorous way, the behavior of the specific retroreflectors we used for target acquisitio