Calibration Alignment Sensitivity in Corneal Terahertz Imaging

Improving the longitudinal modes coupling in layered spherical structure contributes significantly to corneal terahertz sensing, which plays a crucial role in the early diagnosis of cornea dystrophies. Using a steel sphere to calibrate reflection from the cornea sample assists in enhancing the resolution of longitudinal modes. The requirement and challenges toward applying the calibration sphere are introduced and addressed. Six corneas with different properties are spotted to study the effect of perturbations in the calibration sphere in a frequency range from 100 GHz to 600 GHz. A particle-swarm optimization algorithm is employed to quantify corneal characteristics considering cases of accurately calibrated and perturbed calibrated scenarios. For the first case, the study is carried out with signal-to-noise values of 40 dB, 50 dB and 60 dB at waveguide bands WR-5.1, WR-3.4, and WR-2.2. As expected, better estimation is achieved in high-SNR cases. Furthermore, the lower waveguide band is revealed as the most proper band for the assessment of corneal features. For perturbed cases, the analysis is continued for the noise level of 60 dB in the three waveguide bands. Consequently, the error in the estimation of corneal properties rises significantly (around 30%).


Introduction
There has been a surge in THz application in medical imaging and sensing. One of these promising areas is using this spectrum for quantification of the human cornea water content and thickness, which contributes to the early diagnosis of cornea dystrophies [1][2][3][4][5][6]. The nominal cornea dimension (0.5 mm corneal phantom sitting on a 7.5 mm aqueous core) manifests as a lossy thin film at THz frequencies, hence the extraction of corneal central thickness (CCT), corneal anterior water content (CAWC), and corneal posterior water content (CPWC) is achievable through the coupling of longitudinal modes via broadband reflectometry [7]. Given the cornea anatomy acting as a half-space structure in THz band, back-scattered field from the cornea is achievable by stratified medium theory [8][9][10], although some challenges such as phase matching in peripheral areas needed to be addressed. A study based on Fourier optics and vector spherical harmonics enables accurate calculation of back-scattered field from the cornea, by modelling cornea as a layered sphere under Gaussian illumination [11,12]. In practice, the two-way propagation through the reflectometer quasioptics is frequency-dependent, and calibrator targets with known reflectivity are necessary. The studies suggest using a PEC sphere identical to the cornea radius of curvature (ROC) and center of curvature (COC) for calibration [9,12], but reaching this identical situation comes with challenges.
Extracting cornea properties from the back-reflected beam (computed by stratified medium theory) is proposed using a particle-swarm optimization (PSO) algorithm and suggests the WR-5.1 frequency band is more compatible for corneal sensing in [13,14].
In this paper, the theory to compute the coupling efficiency as criteria of resolution of longitudinal modes is briefly described. Six corneal properties which are targeted for study are introduced and compared with the plane-wave model after correct calibration. Different scenarios that can perturb the calibration with the PEC sphere in a realistic situation are introduced and the consequences of these perturbations are addressed by comparing with the plane-wave model. Next, the PSO algorithm is applied to correctly calibrated simulation with added signal to noise ratios (SNRs), to quantify target properties in different waveguide bands. Finally, the root-mean-square deviation (RMSD) of six targets parameters after perturbed calibration in different waveguide bands are reported in the tables.

Theory Based on Fourier Optic and Vector Spherical Harmonics
To solve the scattering problem of a layered sphere (representing a cornea) under Gaussian illumination, the methodology introduced in [15][16][17] appears to be the best choice since the corneal dimensions are of the order of the wavelength of the incident beam. The incident field E i and scattered field E s are obtained as: where a e , a o , b e , and b o are the incident field coefficients. The f e , f o , g e , and g o are scattered field coefficients. The D is a normalization factor. The n and m represent radial and azimuthal mode numbers. The M 1 , N 1 , M 3 , and N 3 are vector spherical harmonics of first and third kind [18]. A through details of the theory can be found in [11,12,[15][16][17]19]. Coupling between the incident and back-scattered field (along optical axis z) are calculated by: in a plane perpendicular to the optical axis, where * denotes the complex conjugate and the incident field E i and scattered field E s described by Equations (1) and (2).

Correctly Calibrated Cornea
Cornea modelled as a layered corneal phantom sitting on an aqueous core with different characteristics listed in Table 1. The cornea ROC is set as 7.5 mm at the anterior and the cornea COC is fixed on x 0 = 0, y 0 = 0, and z 0 = 0. Corneal phantom thickness is considered either 580 µm or 680 µm corresponding to posterior radius of 6.92 µm and 6.82 µm, as shown in Figure 1. Aqueous-humour permittivity is derived from the double-Debye model [20].  For the corneal phantom, the gelatin hydrogel characterized in [21] is used, where the complex permittivity is a function of water content. The permittivity model is based on the effective medium theory and constitutes two main components: free water and solid content; the latter is a combination of dry gelatin and bound water. Furthermore, it is assumed that the total water content varies linearly from the drier anterior surface to the more hydrated posterior surface. The cornea is discretized at most 20 µm-thick layers, whose permittivity depends on the linear water gradient. Hence, for the phantom thickness of 580 µm, cornea modeled as a core with 29 number of equally distanced layers and for the thickness of 680 µm, cornea modelled as a core with 34 number of equally distanced layer.
Each cornea coupling coefficient is computed with Equation (3) and calibrated against the coupling coefficient of the same size steel sphere, which is accurately located in the cornea location (called Cal 0 ), for the frequency range from 100 GHz to 600 GHz. Targets are illuminated by a Gaussian beam in a way that the beam focus on the sub-confocal point introduced in [10,12], shown in Figure 1. Figure 2 compares calibrated efficiency with the stratified medium theory which models the cornea as a planar half-space illuminated by plane wave [8] at the frequency range of 100-600 GHz. Due to the proper matching between plane-wave phase front and planar boundary, the plane-wave calculation is considered as a reference offering the maximum coupling between the incident and back-scattered field.  Figure 2 shows two important points. First, correct calibration (calibration sphere has identical ROC and COC to the target) delivers a very close coupling coefficient to reference plane-wave model coupling, leading to a high resolution of longitudinal modes. Second, the difference between each target is more distinguishable in lower frequencies.

Perturbed Calibrated Cornea
As mentioned in the previous subsection, an accurate calibration leads to reaching a resolution close to maximum coupling, although there are challenges to attaining a precise calibration. In a realistic situation, the PEC sphere calibration can be perturbed in various ways. Eight different scenarios including calibration COC misalignment and error in calibration ROC are described for further investigation. Perturbation scenarios are divided into 4 categories including 1-calibration ROC variation, fixed COC, 2-calibration ROC variation, fixed apex, 3-calibration transverse misalignment, and 4-calibration axial misalignment.
In category 1, ∓2 mm of error in calibration ROC is considered in the way that calibration sphere COC and cornea COC co-located, called Cal 1 and Cal 2 , respectively. In category 2, an error of ∓0.5 mm in calibration ROC is considered, and calibration COC moved in the way that calibration sphere apex and cornea apex co-located, called Cal 3 and Cal 4 , respectively. In category 3, calibration ROC is kept identical with the cornea and 0.5 mm misalignment of calibration COC in transverse location is considered, once the calibration COC is located in x 0 = 0.5 mm and once the calibration COC is located in y 0 = 0.5 mm, called Cal 5 and Cal 6 , respectively. Finally, in category 4, the calibration COC misalignment is considered in the optical axis (z 0 = ∓0.5), called Cal 7 and Cal 8 , respectively. Figure 3a-h demonstrates these scenarios and compares them with the accurate ROC and COC for calibration. In Figure 4, cornea 1 is used as a target and calibrated with various calibration scenarios at a frequency range of 100-600 GHz. A sawtooth behaviour in the phase of perturbed case 1, 2, 7, and 8 is observed which is an indicator of COC misalignment toward the optical axis or error in calibration ROC. The Cal 1 and Cal 2 phase is corrected by considering ∓2 mm phase shift corresponds to an error in calibration ROC. Furthermore, Cal 7 and Cal 8 phase is corrected by considering the ∓0.5 mm phase shift corresponds to an error in the axial location of calibration COC and is reported in Figure 4b. Phase correction revises the phase consistent with the accurate calibration phase. Interestingly, Cal 3 and Cal 4 did not need any phase correction because the ROC and COC change together, in a way that they cancel out the sawtooth behavior of each other and reveal a phase close to accurate calibration. The Cal 5 and Cal 6 , transverse COC misalignment, affects the amplitude and phase stronger than the rest of the scenarios, showing an error of around 10% to more than 70% over the band with respect to the Cal 0 . Among different calibration scenarios, Cal 1 and Cal 2 are less sensitive to the perturbation and act closer to the Cal 0 after phase correction.

Extraction of Corneal Features
To obtain corneal properties from the reflection spectrum, including CAWC, CPWC and CCT, the PSO algorithm is applied [14]. The Additive white Gaussian noise (AWGN) was added to the cornea and PEC sphere coupling efficiencies to build a real-life situation. The noise power was chosen to mimic scenarios in which the SNR was either 40, 50, or 60 dB. For calibration, the noisy cornea coupling efficiency was divided by the noisy PEC coupling efficiency to make it comparable to an analogous plane-wave model. This operation was repeated for six cornea targets. Then, the PSO algorithm was used to extract the target parameters by searching for the corresponding plane-wave model whose reflectivity minimizes the difference with the noisy simulated reflectivity. The adopted merit function is [22,23]: where e ikz is added to compensate for the deviation of the calibration COC along optical axis z or calibration ROC deviation. Γ Cal is the noisy calibrated reflection coefficient, Γ pw is the plane-wave model reflection coefficient, and N is the number of frequency points. The merit function is a sum of the mean squared error normalized with the average power of the merit function metric.

PSO Analysis for Correctly Calibrated Cornea
For the six cornea properties defined in Table 2, an SNR of 40 dB, 50 dB, and 60 dB were added to target and calibration coupling efficiency and was used as an input to the PSO algorithm to extract the corneal properties. The standard frequency range of waveguide WR 5.1 (140-220 GHz) was considered. The extracted values for cornea CAWC, CPWC, and CCT are shown in Figure 5. The different corneas are plotted with different colours. Nominal values are shown with yellow spots. It is obvious that higher SNR contributes to a more accurate estimation of variables.
In the following, the standard waveguide bands of WR-3.4 (220-330 GHz) and WR-2.2 (330-500 GHz) were compared for SNR of 60 dB, repeating a PSO analysis for six noisy correctly calibrated cornea, and the comparison is depicted in Figure 6. The WR-3.4 provides a better estimation of variables compared to WR-2.2, yet worse than WR-5.1.

PSO Analysis in Case of Perturbation
For more exploration, PSO analysis for perturbed calibrated cornea was performed. The 60 dB SNR was added to six corneas described in Table 1, then calibrated by various noisy perturbation scenarios introduced in Figure 3. Tables 2-4 reports the root-mean-square deviation (RMSD) of estimated CCT, CAWC, and CPWC from their nominal values at the frequency band of WR 5.1, WR 3.4 and WR 2.2, respectively. The consequence of the error in calibration sphere ROC and COC is elucidated in Tables 2-4. The error in parameter estimation rises significantly (around 30%), especially for transverse misalignment calibration, Cal 5 and Cal 6 , consistent with the result in Figure 4. Table 2. RMSD of estimated CCT, CAWC, and CPWC for six noisy cornea calibrated with various noisy perturbed PEC sphere at frequency band of WR 5.1. The SNR 60 dB is considered.

Discussion
This study suggests extending the corneal spectroscopy from the previously studied WR-3.4 band to the WR-5.1 using an accurate computational method based on Fourieroptics and vector spherical harmonics and introducing an acceptable noise level of 60 dB.
Furthermore, in terahertz cornea sensing and imaging, it is suggested to calibrate the cornea reflection with a steel sphere (identical to the radius and location of the target) reflection that contributes to an enhancement in longitudinal modes resolution. The challenge to reach an optimal calibration was discussed and addressed. Perturbed scenarios affect phase significantly. Adding a phase correction coefficient (coincidence to the misalignment in optical axis or ROC calibration) contributed to the proper compensation of the phase, yet phase correction was not adequate for the transverse misalignment.
This study assists in analyzing measured reflection in the laboratory and offers some insights such as avoiding the transverse dislocation of the calibration sphere due to a higher level of the error, also mitigating the misalignment in the optical axis or the deviation of the calibration ROC through distance estimation in the PSO algorithm.

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