*Proceeding Paper* **Comparing two Fitting Algorithms for Determining the Cole–Cole Parameters in Blood Glucose Problems †**

**Roberto Dima \*, Giovanni Buonanno and Raffaele Solimene**

† Presented at the 2nd International Electronic Conference on Applied Sciences, 15–31 October 2021 ; Available online: https://asec2021.sciforum.net/.

**Abstract:** This paper addresses the nonlinear inverse problem of estimating the parameters of the Cole–Cole model used to describe the behavior of the complex permittivity of blood samples. Such a model provides an efficient and accurate representation of biological tissues in the entire frequency band considered and reduces the complexity of the experimental data to a few parameters. In this way, it is possible to extract a "synthetic view" of the dielectric properties of tissues in such a way that more information on the glucose concentration can be derived, in addition to the resonance peak or phase shift. In order to perform the fitting of the Cole–Cole model, two different algorithms were used and compared: the Levenberg–Marquardt and the variable projection algorithms. Synthetic data present in the literature were used to evaluate the performances obtainable with these methods. In particular, Monte Carlo analysis was used in order to evaluate the accuracy and the precision that these two methods provide in the process of estimating the parameters involved, with respect to the starting points of the parameters. The results obtained showed that the variable projection algorithm always outperformed the Levenberg–Marquardt one, although the former has a greater computational burden than the latter.

**Keywords:** glucose measurement; Cole–Cole model; Levenberg–Marquardt algorithm; variable projection algorithm; blood dielectric properties; nonlinear fitting problem

#### **1. Introduction**

Diabetes is a metabolic disorder that afflicts millions of people in the world [1]. It degrades the cell's ability to absorb glucose from the bloodstream because of the improper regulation of the insulin hormone. For this reason, great efforts have been dedicated to the development of non-invasive glucose monitoring devices, which may considerably improve the quality of life for diabetics [2].

This work is particularly concerned with microwave sensor technology that relies on the change in the dielectric and conductivity properties of blood plasma as a function of the glucose concentration in order to track such a change.

In this framework, developing accurate and precise fitting methods for blood models, at different glucose concentrations, is essential for the development of robust electromagnetic (EM)-based techniques that could be employed for non-invasive, continuous glucose monitoring. Indeed, accurate electromagnetic tissue modeling is of paramount importance since it affects the simulation stage required for sensor design [3]. Moreover, extracting a "synthetic view" (in terms of a few parameters) of the sensor response data is essential for analyzing patterns and possibly extracting more information, besides the resonance peak or phase shift, about glucose concentration.

In this paper, the aim was just to address the fitting problem. More in detail, starting from the dielectric spectrum, which is assumed known over a certain number of frequencies, we aimed at estimating the parameters of a single-pole Cole–Cole model. As is well

**Citation:** Dima, R.; Buonanno, G.; Solimene, R. Comparing two Fitting Algorithms for Determining the Cole–Cole Parameters in Blood Glucose Problems. *Eng. Proc.* **2021**, *11*, 45. https://doi.org/10.3390/ ASEC2021-11188

Academic Editor: Roger Narayan

Published: 15 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Department of Engineering, University of Campania, 81031 Aversa, Italy; giovanni.buonanno@unicampania.it (G.B.); raffaele.solimene@unicampania.it (R.S.) **\*** Correspondence: roberto.dima@unicampania.it

known, this entails solving a nonlinear inverse problem, which here was addressed by two different methods: the classical Levenberg–Marquardt method [4,5] and the variable projection algorithm [6]. We evaluated how sensitive the two methods are with respect to the starting points of the parameters and with what accuracy and precision these parameters can be estimated. In order to check the two methods, we first generated synthetic relative permittivities by employing a single-pole Cole–Cole model, using data from the literature [7] as true values for its parameters, and then solved an inverse problem in order to trace these values by resorting to the two aforesaid methods.

Although higher-order models could perform better, we considered a first-order model to perform the comparison in the simplest possible case.

#### **2. Methods**

#### *2.1. Cole–Cole Model*

The Cole–Cole model [8] is widely used to describe the complex relative permittivity of biological tissues, *ε*r(*ω*) = *ε*(*ω*)/*ε*0, and its equation is:

$$\varepsilon\_{\rm r}(\omega) = \varepsilon\_{\infty} + \sum\_{n=1}^{N} \frac{\varepsilon\_{\rm s\rm t} - \varepsilon\_{\infty}}{1 + (j\omega \tau\_n)^{1 - a\_n}} - \frac{\sigma}{j\omega \varepsilon\_0} \tag{1}$$

in which *N* is the number of poles and thus the order of the model, *ε*<sup>∞</sup> = lim*ω*→<sup>∞</sup> *ε*r(*ω*) is the permittivity at high frequencies, *σ* is the static ionic conductivity and *ε*s*<sup>n</sup>* = lim*ω*→<sup>0</sup> *ε*r(*ω*), and *τn*, and *α<sup>n</sup>* are the static permittivity, the relaxation time constant, and the so-called distribution parameter of the *n*-th addend of the summation, respectively. Such a model incorporates the Debye model [9]. Indeed, the main difference between the Debye and the Cole–Cole models is that the latter includes the exponent 1 − *α*, with 0 ≤ *α* ≤ 1. When the exponent becomes smaller, the relaxation time distribution becomes broader, i.e., the transition between low- and high-frequency values becomes wider and the peak on the imaginary part of the spectrum also becomes wider.

The complexity of both the structure and composition of biological material is such that the dispersion region of each pole may be broadened by multiple contributions to it. The broadening of the dispersion could be empirically accounted for by using the Cole–Cole model [10], which is expected to give more accurate dielectric spectrum curve fitting.

#### *2.2. Curve Fitting Algorithms*

Let be *x* the vector of model parameters and *P* its length and *M* the number of frequency points for which the measures are taken. We define the data vector ( stands for transposition) *y* = *y*(*ω*1)... *y*(*ωm*)... *y*(*ωM*) in which the *m*-th component of the vector *y* is the observed value *y*(*ωm*). Let also *ε*<sup>r</sup> = *ε*r(*ω*1; *x*)... *ε*r(*ωm*; *x*)... *ε*r(*ωM*; *x*) be the model vector, here given by Equation (1)), with *ε*r(*ωm*; *x*) being the estimation at *ωm*.

Solving the least-squares problem means finding *x***ˆ** such that:

$$\mathfrak{X} = \arg\min\_{\mathbf{x} \in \mathbb{R}^P} \left\{ \frac{1}{2} \left\| \boldsymbol{\varepsilon}\_{\mathbf{r}}(\mathbf{x}) - \boldsymbol{y} \right\|\_2^2 \right\} \tag{2}$$

in which the function to minimize, Ψ = <sup>1</sup> <sup>2</sup> *ε*r(*x*) <sup>−</sup> *<sup>y</sup>*<sup>2</sup> <sup>2</sup>, is the <sup>2</sup> quadratic norm of the misfit *<sup>r</sup>* <sup>=</sup> *<sup>ε</sup>*r(*x*) <sup>−</sup> *<sup>y</sup>*, which is a nonlinear function such that *<sup>r</sup>* : <sup>R</sup>*<sup>P</sup>* → <sup>R</sup>*<sup>M</sup>* with *<sup>P</sup> <sup>M</sup>*.

We addressed the nonlinear fitting problem with two methods: the Levenberg– Marquardt algorithm (LMA) and the variable projection algorithm (VPA).

The Levenberg–Marquardt algorithm [4,5] acts more as a gradient-descent method when the parameters are far from their optimal value and acts more as the Gauss–Newton method when the parameters are close to their optimal value [11]. The equation for the step *h* at the *k*th iteration is:

$$\left(J(\mathbf{x}\_k)^\top J(\mathbf{x}\_k) + \lambda\_k I\right)\mathbf{h} = -J(\mathbf{x}\_k)^\top f(\mathbf{x}\_k) \tag{3}$$

where *J* is the Jacobian of *f* and *λ<sup>k</sup>* is the damping parameter. It controls both the magnitude and direction of *h*, and it is chosen at each iteration. It can be shown [5] that, at each iteration, Equation (3) solves the minimization problem over a reduced set of admissible solutions, i.e., those that satisfy *h* <sup>≤</sup> *<sup>R</sup>*(*λ*), limiting the correction step to within a region near *xk*. The radius of the trust region *R* = *R*(*λ*) is a strictly decreasing function with lim*λ*→<sup>∞</sup> *<sup>R</sup>*(*λ*) = 0. When *<sup>λ</sup><sup>k</sup>* <sup>=</sup> 0, the step *<sup>h</sup>* is identical to that of the Gauss–Newton method, i.e., the same direction and maximum magnitude. As *<sup>λ</sup>* <sup>→</sup> <sup>∞</sup>, *<sup>h</sup>* tends towards the steepest descent direction, with the magnitude tending towards 0.

Based on the above, we inferred the qualitative update rule for *λk*<sup>+</sup>1: If Ψ(*x<sup>k</sup>* + *h*) < Ψ(*xk*), then the quadratic approximation works well, and we can extend the trust region, i.e., it will be *λk*+<sup>1</sup> < *λk*. Otherwise, the step is unsuccessful and we reduce the trust region, i.e., it will be *λk*+<sup>1</sup> > *λk*; in this way, the next step tends towards the negative gradient method, and a lower value of Ψ is more likely to be found.

The MATLAB implementation was used, in particular the lscurvefit function with the Levenberg–Marquardt option.

The variable projection algorithm [6] is a method used to solve separable nonlinear least-squares problems. The least-squares problem is said to be separable when the model parameters can be separated into two sets of parameters, one that enters linearly into the model, *c* = [*c*1, ... , *ck*], and another set of parameters that enter the model nonlinearly, *a* = [*a*1, ... , *al*], so that *x* = [*c*, *a*]. For each observation *ym* of a separable nonlinear leastsquares problem, the model is a linear combination of nonlinear functions that depend on nonlinear parameters, and the model function can be written as *ε*r(*ω*) = ∑*<sup>k</sup> <sup>j</sup>*=<sup>1</sup> *cj <sup>φ</sup>j*(*ω*; *<sup>a</sup>*).

The functional Ψ is written in terms of residual vector *r* as:

$$\Psi(a,\mathbf{c}) = \frac{1}{2} \left\| \mathbf{y} - \Phi(a)\mathbf{c} \right\|^2 \tag{4}$$

in which the columns of the matrix **Φ** are the nonlinear functions *φj*(*ω*; *a*). The linear parameters *c* could be obtained from the knowledge of *a*, by solving the linear least-squares problem:

$$\mathfrak{c} = \Phi(a)^\dagger y \tag{5}$$

which stands for the minimum-norm solution of the linear least-squares problem for fixed *a*, where **Φ**(*a*)† is the Moore–Penrose generalized inverse of **Φ**(*a*). By replacing this in Equation (4), we obtain the variable projection functional:

$$\Psi\_{\rm VP}(a) = \frac{1}{2} \left\| y - \Phi(a)\Phi(a)^\dagger y \right\|^2 \tag{6}$$

The variable projection algorithm consists of two steps: first minimizing Equation (6) with an iterative nonlinear method and then using the optimal value found for *a* to solve for *c* in Equation (5) [12]. The principal advantage is that the iterative nonlinear algorithm used to solve the first minimization problem works in a reduced space and less initial guesses are necessary. A robust implementation in MATLAB, called VARPRO [13], was adapted and used to deal with complex-valued problems, choosing the Levenberg–Marquardt option for the solution of Equation (6).

#### *2.3. Numerical Simulations*

The generation of the synthetic complex relative permittivity of blood plasma relied on quadratic fits to glucose-dependent Cole–Cole parameters reported in [7]; in particular, we considered two different concentrations, 100 mg/dL and 250 mg/dL, one normal and the other of severe diabetes, respectively, in accordance with the diagnostic criteria in [1]. The data vector consisted of *M* = 1000 points in the frequency range 500 MHz–20 GHz.

In gradient-like algorithms (such as those used in this paper), the choice of the initial point is a crucial factor for the convergence of the procedure. For the single-pole model case, it is fairly easy to exploit the physical meaning of the parameters to infer an initial estimate. However, since the noise can invalidate the initial estimate, we propose to study the robustness of the two algorithms with respect to the initial point. To this end, a Monte Carlo analysis was performed, iteratively evaluating the deterministic algorithms using a set of N = 1000 uniformly distributed random initial points arranged in a 5D hypercube of the parameter space in order to statistically characterize the results. Each side of the hypercube represents an interval containing the range of variation of each parameter for the glucose concentrations considered.

We ran simulations on a machine with Intel i9-10850K (10 physical cores), 32 GB RAM, and Ubuntu 21.04, and we took advantage of the Parallel Computing Toolbox, using the parfor loop for running the 1000 simulations.

The intervals for generating the random initial value for each parameter (of the Cole– Cole model) were chosen from the data tabulated in [7]. In particular, the widths of these intervals were the same for each glucose concentration and were: [1, 5] for *ε*∞, [1, 150] for *<sup>ε</sup>s*, [<sup>1</sup> <sup>×</sup> <sup>10</sup>−14, <sup>1</sup> <sup>×</sup> <sup>10</sup>−11] for *<sup>τ</sup>*, [0.1 <sup>−</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−9, 0.1 <sup>+</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−9] for *<sup>α</sup>*, and [0, 5] for *<sup>σ</sup>*. These intervals are relatively large compared to the values taken from [7] in order to test the two algorithms in sufficiently stressful situations. Only the range of variation of *α* is extremely small because the model used in [7] practically fixes it a priori to 0.1. Obviously, it must be taken into account that the VPA requires only generating the values for *τ* and *α*.

For an initial bland qualitative assessment, we established evaluation intervals (the same for generations) for the estimated parameters so that we could assert that a reconstruction is "good" if it falls within these ranges, and "wrong" otherwise.

Now, let *<sup>x</sup>*<sup>ˆ</sup> (*i*) = [*<sup>ε</sup>* (*i*) <sup>∞</sup> ,*<sup>ε</sup>* (*i*) *<sup>s</sup>* , *<sup>τ</sup>*(*i*), *<sup>α</sup>* (*i*), *<sup>σ</sup>* (*i*)] be the vector of the parameter estimates returned by the two algorithms at the *i*-th simulation, and let *x*ˆ(*i*) denote one of its five elements. Moreover, let *x* = (1/*N*sim) <sup>∑</sup>*N*sim *<sup>i</sup>*=<sup>1</sup> *x*ˆ*<sup>i</sup>* and *σ<sup>x</sup>* = [1/(*N*sim <sup>−</sup> <sup>1</sup>)] <sup>∑</sup>*N*sim *<sup>i</sup>*=<sup>1</sup> <sup>|</sup>*x*ˆ*<sup>i</sup>* <sup>−</sup> *x*|<sup>2</sup> be the sample mean and standard deviation, respectively, calculated for each parameter.

For a quantitative evaluation of the performance of the two algorithms, we then define multiple figures of merit for statistically characterizing the results of Monte Carlo analysis. For each parameter, Equations (7) and (8) define measures of accuracy and precision, respectively, defined over the entire set of reconstructions. However, such measures can be greatly affected by estimates that are very far from the true value, *xtrue*, where the latter represents one of the five elements of the vector of reference values *<sup>x</sup>*ˆ*true* = [*<sup>ε</sup>*∞*true*,*<sup>ε</sup>strue*, *<sup>τ</sup>true*, *<sup>α</sup>true*, *<sup>σ</sup>true*]. For this reason, we also introduced <sup>A</sup>*cut* and <sup>P</sup>*cut* in Equations (9) and (10) in order to define the accuracy and precision measures, respectively, which instead dampen the effect of the above isolated events. The subscript *cut* means that they were calculated on a subset obtained by eliminating *ζ*% of reconstructions with lower values and *ζ*% of reconstructions with higher values, in which 0 < *ζ* < 50.

$$\mathcal{A} = \left| \frac{\mathbf{x}\_{true} - \langle \mathbf{x} \rangle}{\mathbf{x}\_{true}} \right| \times 100\% \tag{7}$$

$$\mathcal{P} = \frac{\sigma\_{\text{x}}}{\langle \mathbf{x} \rangle} \times 100\% \tag{8}$$

$$\mathcal{A}\_{\rm cut} = \left| \frac{\chi\_{\rm true} - \langle \mathbf{x} \rangle\_{\rm cut}}{\chi\_{\rm true}} \right| \times 100\% \tag{9}$$

$$\mathcal{P}\_{\rm cut} = \frac{\sigma\_{\rm x\_{\rm cut}}}{\langle \mathbf{x} \rangle\_{\rm cut}} \times 100\% \tag{10}$$

#### **3. Results**

We conducted many numerical simulations by widening more and more the generation intervals. In this paper, we report the case where the generation intervals are very large except for the *α* interval, due to the above explanation. In all these experiments, following the qualitative criterion mentioned above, the VPA always provided good estimations, while the same was not true for the LMA. In particular, in the case considered, the LMA provided the wrong estimates in about 260 simulations out of 1000, for each of the two

glucose concentrations, while no wrong estimate was returned by the VPA. Here, by the "wrong" estimate, we mean that at least one component of the parameter vector *x* had a value that was outside its generation range. Graphical representations of those qualitative results are provided in Figure 1.

**Figure 1.** Graphical representations of the convergence of the LMA and VPA for 1000 simulations and 2 different glucose concentrations: (**a**) the LMA for 100 mg/dL, (**b**) the LMA for 250 mg/dL, (**c**) the VPA for 100 mg/dL, and (**d**) the VPA for 250 mg/dL.

Consistent with the qualitative results, considering the whole set of 1000 estimates, the VPA exhibited excellent accuracy and precision, while this was not the case for LMA, as can be observed in Tables 1 and 2. For this reason, for the VPA only, the results calculated by means of Equations (7) and (8) are reported, whilst for the LMA, the results derived from Equations (9) and (10), obtained by cutting 25% of the lowest values and 25% of the highest values, were also considered.

Regarding execution times, the LMA took 7 s to finish 1000 simulations, while the VPA took around 170 s.


**Table 1.** Figures of merit (in %) of the two algorithm in which the glucose concentration is 100 mg/dL.

**Table 2.** Figures of merit (in %) of the two algorithm in which the glucose concentration is 250 mg/dL.


#### **4. Discussion**

In this paper, we faced the problem of fitting the dielectric spectrum of a blood sample in order to estimate the parameters of the single-pole Cole–Cole model. In particular, we compared the performance of two different algorithms, the LMA and VPA, in terms of the accuracy and precision with respect to the starting points of the parameters.

For the parameter range considered, the VPA outperformed the LMA in robustness with respect to the initial point of the algorithm. However, analyzing the figures of merit related to the LMA, it became clear that there were erroneous reconstructions so far from the true value such that they heavily deteriorated the (standard) accuracy and precision, while the cut versions did not suffer from this problem. In fact, once the erroneous ones were removed, in all other simulations, the algorithm converged to the true values.

On the other hand, the VPA gave these good results because less initial guesses were necessary and because the iterative nonlinear algorithm used to solve the first minimization problem works in a reduced space. The big disadvantage of the VPA, or at least of the implementation used in this paper, was the execution time. This was certainly due to the numerous singular-value decompositions (SVDs) that the algorithm calculates in its runtime.

The results are promising, and the research will continue by evaluating the algorithms in increasingly realistic scenarios, including adding noise to the synthetic data.

**Author Contributions:** Conceptualization, R.S.; methodology, R.D. and G.B.; software, R.D.; validation, R.D. and G.B.; formal analysis, R.S., R.D. and G.B.; investigation, R.D. and G.B.; resources, R.S.; data curation, R.D.; writing—original draft preparation, R.D.; writing—review and editing, G.B. and R.D.; visualization, R.D. and G.B.; supervision, R.S. and G.B.; project administration, R.S.; funding acquisition, R.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by MIUR, Project Grant PRIN 2017 "Microwave Biosensors: Enhanced Non-Invasive Methodology for Blood Glucose Monitoring". The APC was funded by the above Project Grant PRIN 2017.

**Data Availability Statement:** The data presented in this study are available upon request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Proceeding Paper* **Drone Polariscopy—Towards Remote Sensing Applications †**

**Soon Hock Ng 1,\*, Blake Allan 2, Daniel Ierodiaconou 2, Vijayakumar Anand 1, Alexander Babanin <sup>3</sup> and Saulius Juodkazis 1,4,\***


**Abstract:** Remote sensing is critical for a wide range of applications, including ocean and wave monitoring, planetary exploration, agriculture, and astronomy. We demonstrate a polariscopy concept that is able to determine orientation of patterns below the optical resolution limit of a system. This technique relies on measuring at least four different polarisation angles and calculating the orientation from this set of intensity information. It was initially demonstrated on the Infrared Microspectroscopy Beamline at the Australian Synchrotron using IR light in transmission. Using a monochrome polarising camera mounted onto a drone as a remote sensing platform analogue, orientation information was extracted from 3D-printed targets in reflection. The images were taken at an altitude where conventional imaging could not resolve the test patterns. The system had a 3.33 mm ground resolution. Patterns consisting of 0.5 mm lines spaced 0.5 mm apart were detected using the method, demonstrating the capability of detecting features over six times smaller than the resolution limit. In the interest of moving towards high-speed data acquisition and processing, two methods for processing the image are compared—an analytical and a curve fitting method.

**Keywords:** polariscopy; remote sensing; drone; image processing

Ierodiaconou, D.; Anand, V.; Babanin, A.; Juodkazis, S. Drone Polariscopy— Towards Remote Sensing Applications. *Eng. Proc.* **2021**, *11*, 46. https:// doi.org/10.3390/ASEC2021-11161

**Citation:** Ng, S.H.; Allan, B.;

Academic Editors: Nunzio Cennamo and Nicholas Vassiliou Sarlis

Published: 15 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

The ability to detect periodic features below the resolution limit of a system is valuable to all forms of imaging. This is especially important in areas such as remote sensing, where the cost of launching a satellite is dependent on its mass. Polarisation information adds an additional dimension to image data and has shown value in machine vision, where it is utilised for analysis of reflections [1]. In remote sensing, the orientation detection of waves is studied using synthetic aperture radar (SAR). The polarisation orientation has been used on an existing resolved image [2,3]. These satellites need large, heavy antennas to achieve their resolutions. However, in such a contexts, polarisation has not been used for diffraction-limited imaging of anisotropy. A method using Fourier transform infrared spectroscopy was developed and transferred to the Infrared Microspectroscopy Beamline at the Australian Synchrotron, which demonstrated anisotropy recognition and mapping below the spatial resolution limit for the first time [4,5]. The method described—the 4-pol method—involves measuring the sample at four polarisation angles (0◦, 45◦, 90◦, and 135◦ (−45◦)). In addition to molecular anisotropy, the orientation of a circular grating with a 200 nm pitch and 100 nm line width was detectable with a system resolution of 5 μm, a 25 × difference [6]. To demonstrate the ability of the 4-pol method to move beyond the IR wavelengths and microscopy scale, a 4-pol visible light camera was used to show orientation detection in the visible wavelength range [7]. In this preliminary study, we extended the method even further and mounted the camera on a drone, which functioned as a remote sensing platform analogue and operated in reflection mode.

#### **2. Materials and Methods**

#### *2.1. Drone Flight with Polarisation Camera*

The CS505MUP monochrome polarisation camera (Thorlabs Inc.) has a 2448 × 2048 pixel sensor with micro-polarisers directly over each 3.5 μm pixel. Each 2 × 2 group of pixels provides 0◦, 45◦, 90◦, and 135◦ intensity information for a total of 1224 × 1024 for each polarisation angle. A Navitar MVL8M23 f/1.4, 8 mm lens was attached to the camera, which was then mounted on a custom mounting rig which provided vibration damping, a battery, a LattePanda computer with a receiver, and a GoPro video camera (Figure 1a). This rig was then bolted under a Kraken 130V2 octocopter drone platform (Figure 1b). After takeoff, the acquisition was remotely triggered to conserve power and storage space. To account for different lighting conditions, the camera was programmed to record a bursts of 10 images at varying exposure times, every 30 s. The drone was flown up in steps of 20 m and hovered for 1 min to ensure a full set of images were acquired at each altitude. 3D-printed targets (Figure 2) consisting of 0.5 mm lines with a 1 mm pitch (one circular with an Archimedean spiral and one rectangular with different orientations in each quadrant) were set up on the ground as points of known orientation. The diameter of the circle was 20 cm and the dimensions of the rectangle were 24 × 20 cm.

**Figure 1.** (**a**) Photograph of the polarisation camera in its mounting rig along with a remote control. (**b**) Preflight photo of the octocopter with camera module attached underneath. Inset shows the drone and camera just after takeoff.

**Figure 2.** (**a**) Intensity image of 3D-printed targets on a grass field at an altitude of ∼4 m. Inset shows a photo of the centre of the rectangular target showing the orientation in each quadrant. (**b**) Corresponding azimuth map with the colours indicating azimuth angle. (**c**) Enlarged intensity image of the targets. (**d**) Enlarged azimuth image of the targets with black arrows indicating approximate calculated azimuth orientation.

#### *2.2. Data Processing*

The data were subsequently converted to intensity and orientation information using two methods—analytical and curve fitting [4,8]. For the analytical solution, the intensity is given by

$$I = \frac{I\_{\rm O^{\circ}} + I\_{45^{\circ}} + I\_{90^{\circ}} + I\_{135^{\circ}}}{2},\tag{1}$$

where *I* is the image intensity and *Iθ*◦ is the intensity at a specific polarisation angle. The orientation (azimuth) is given by

$$\psi = \frac{1}{2} \arctan2\left(\frac{I\_{45^\circ} - I\_{135^\circ}}{I\_{0^\circ} - I\_{90^\circ}}\right),\tag{2}$$

where arctan2 is the four-quadrant inverse tangent. The equation used for fitting was

$$I(\theta) = A\cos(2\theta - 2\psi) + c,\tag{3}$$

where *I*(*θ*) is the intensity at a specific polariser angle, *A* is a factor proportional to the degree of polarisation, *θ* is the polariser angle, *ψ* is the azimuth, and *c* is a factor proportional to *<sup>I</sup>*. The bounds of the fit were 0 <sup>&</sup>lt; *<sup>A</sup>* <sup>&</sup>lt; 5, <sup>−</sup>*<sup>π</sup>* <sup>2</sup> <sup>&</sup>lt; *<sup>ψ</sup>* <sup>&</sup>lt; *<sup>π</sup>* <sup>2</sup> , and −5 < *c* < 5. The initial guesses were *A* = 1, *ψ* = 0, and *c* = 1. In all cases, these were the values used for all pixels.

#### **3. Results and Discussion**

#### *3.1. Orientation Determination*

Figure 2a shows an intensity image, calculated with Equation (1), at an altitude of approximately 4 m, and the inset shows a photo of 3D-printed lines on the rectangular target. This height is equivalent to 3 pixels/cm in the image or 6 pixels/cm over the full sensor (calculated using the 20 cm diameter circle as a reference). The highest theoretical resolution of a sensor is when two points are imaged onto separate neighbouring pixels. Naively for the above sensor with a pixel size of 3.5 μm, the sensor resolution *R*sensor can be considered 3.5 μm, although physical limitations preclude this. The polarisation image (equivalent pixel size of 6.9 μm) has half the resolution of the full sensor. The magnification factor

$$m = \frac{\text{Sensor width}}{\text{Image width}}$$

between the object space and sensor can be used to calculate the maximum object resolution *R*object. For the image in Figure 2a,

$$m = \frac{8.445\text{ mm}}{4080\text{ mm}} = 0.00207.5$$

Hence,

$$R\_{\text{object}} = \frac{3.45 \text{ } \text{\textmu m}}{0.00207} = 1.67 \text{ mm}$$

for the full sensor and 3.33 mm for the polarisation intensity image.

Since the spaces between the lines were 0.5 mm, the lines could not be resolved in the image. Applying Equation (2) resulted in Figure 2b, which shows colours representing the azimuth angle with 0◦ defined as shown in the bottom left. It is evident that each quadrant of the rectangular target shows a different orientation, and the circular target shows radial changes. Figure 2c,d shows enlarged versions of the region of interest. Arrows in Figure 2d show approximate azimuth orientation. The azimuth is orthogonal to the orientation of the printed lines for both rectangular and circular targets. This is in contrast to transmission mode, where the calculated azimuth is parallel to the alignment [6]. However, the fact that the azimuth is consistent still indicates that the 4-pol method detects orientation that cannot be resolved by the optics in reflection. This system, with a theoretical 3.33 mm resolution, was able to determine the orientation of 0.5 mm lines spaced 0.5 mm apart, demonstrating

detection of patterns over six times smaller than the resolution limit. It also shows that it works in cases where either the detected features are much smaller or much larger than the wavelength of the probing light. The difference in azimuth alignment might have been due to the illumination from the sky, since polarisation varies with angle from the sun [9].

#### *3.2. Data Processing Method*

There are cases wherein there is a low signal-to-noise ratio, the angle of a polariser cannot be set precisely, or the polarisation combinations are not 45◦ apart. The analytical method cannot accommodate these, and so curve fitting can be used to determine the azimuth and intensity. This comes at a time cost, however. While Equations (1) and (2) can be vectorised so the calculation can be performed on an entire image at once, curve fitting needs to be done for each pixel. As a single-thread process, the analytical method takes <1 s to process 1 frame. The same frame using curve fitting takes over 3 h.

While the fitting method is slow, it ultimately returns the very similar results to the analytical method, when the bounds are set correctly. Figure 3a shows an enlarged cropped section from Figure 2b, and Figure 3b shows the same area but calculated by fitting each pixel. At a glance, they are indistinguishable; however, Figure 3c shows that there are several pixels which are different. These differences are not very significant, as the number of pixels that do not agree are few. The time disparity demonstrates that the fitting method should only be used when a great number of polarisation measurements are needed—when the signal-to-noise ratio is low.

**Figure 3.** (**a**) Azimuth image calculated using Equation (2). (**b**) Azimuth image calculated by fitting Equation (3). (**c**) Images (**a**,**b**) overlaid with the difference image blend mode showing pixels that are different in colour and pixels which are the same as black. Arrows are eye guides to one pixel which is different.

#### **4. Conclusions**

The 4-pol method was successfully transferred from the Infrared Microspectroscopy Beamline to a drone platform. It was shown that orientation could be determined even when the aligned features (0.5 mm lines 0.5 mm apart) could not be spatially resolved (3.33 mm resolution). This demonstrated that patterns over six times smaller than the system's resolution limit could be extracted from the image. The calculated azimuth was orthogonal to the alignment direction, and the reason for this requires further study. This is a stepping stone towards satellite-based remote sensing using the same technique, and opens up possibilities for ocean monitoring, planetary observation, and astronomy, where orientation information is important.

**Author Contributions:** Conceptualisation, S.J. and S.H.N.; methodology, S.H.N.; software, S.H.N., V.A.; investigation, S.H.N., B.A., D.I., V.A., A.B. and S.J.; resources, B.A., D.I. and A.B.; data curation, S.H.N.; writing—original draft preparation, S.H.N. and S.J.; writing—review and editing, S.H.N., S.J., B.A., D.I. and V.A.; visualization, S.H.N.; supervision, D.I., A.B. and S.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Australian Research Council, grant number LP190100505.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

