© The Author(s), published by EDP Sciences, 2025

1 Introduction

Refractive surgery is a broad term encompassing various surgical procedures to correct vision. The primary goal is to reduce or eliminate the need for eyeglasses or contact lenses and to enhance the visual acuity [1–3]. Femtosecond lasers have advanced the field of refractive surgery by enabling precise and reproducible tissue cutting, as supported by several studies (e.g., [4, 5]).

In refractive surgery, distinct substructures are created depending on the procedure. For instance, the optical breakdown procedure generates an intrastromal lenticule that is extracted through a small incision, while LASIK involves the creation of a corneal flap, followed by stromal ablation with an excimer laser. While the mechanical characteristics of the substructures differ, understanding their properties is critical for surgical outcomes.

One of the earliest applications of femtosecond lasers in refractive surgery is flap creation for LASIK. Compared to mechanical microkeratomes, femtosecond lasers offer several advantages, including a reduced risk of severe flap breaks during surgery and the ability to customize flap shape [6, 7].

Lenticule extraction is a well-established alternative for laser vision correction. This procedure uses an ultrashort pulse laser system to define the shape of a tissue volume that will be removed to achieve refractive correction [8, 9].

A precise lenticule design ensures it can be removed without damaging the surrounding tissue, thus lowering risks such as corneal dryness, infection, or refractive errors [10–12]. Flap dimensions such as size and thickness can influence visual recovery [13]. Flap morphology characterization can reduce potential risks of postoperative flap striae and displacements [14].

The geometry characterization of flap and lenticule cuts is an associated procedure playing a crucial role in refractive surgeries. The geometric parameters such as flap diameter and thickness are of importance for achieving the desired refractive outcomes. For instance, the semi-automated measurement of central values was studied by Lwowski et al. [15].

The geometry of lenticule cuts determines the amount of corneal tissue removed or reshaped. Accurate characterization ensures the desired refractive correction, minimizing postoperative refractive errors and enhancing predictability. Characterizing the geometry reduces the risks associated with surgical interventions as well as explaining probable post-operative issues originating from imperfect residuals. Previous systematic studies predominantly relied on manual expert assessment for geometry characterization [16–20]. These approaches often involved subjective interpretation and significant variability due to human factors. Consequently, they could be less consistent and reproducible compared to automated methods. Our proposed approach aims to address these limitations through advanced image processing techniques.

Optical Coherence Tomography (OCT) is an essential imaging technology that has significantly assisted the field of medicine [21–24]. By providing non-invasive, high-resolution cross-sectional images, OCT images enable medical professionals to accurately analyze internal structures.

OCT imaging enables high-resolution visualization of the eye’s microstructures, establishing it as an indispensable tool for both diagnostic and treatment planning purposes [25–27]. Its ability to deliver real-time imaging with micrometer-level precision has elevated the standard of care in ophthalmology [28–32].

In refractive surgery, OCT imaging plays a pivotal role throughout the surgical process, from preoperative planning to postoperative assessment [2]. It provides critical data on corneal thickness, curvature, and epithelial profile. During surgery, OCT assists in verifying lenticule cuts or flap quality, while postoperatively, it aids in detecting and managing complications [33].

By providing high-resolution imaging of ocular structures, the use of OCT enables clinicians to make more informed decisions, which can contribute to improved surgical outcomes [34, 35]. However, challenges such as the presence of noise have been addressed through various advancements in imaging techniques, leading to enhanced image quality and more accurate assessments [36–39].

Speckle noise along with the scattering and small contrast through intrastromal OCT images complicates accurate segmentation (characterization) of peri-operative cuts [40–42].

Currently, manual measurements serve as the reference, but they can potentially introduce systematic inaccuracies. These arise from user precision and arbitrary marker positioning relative to intrastromal surfaces, leading to significant variations in similar measurements.

We propose a combined approach that involves denoising images and optimization to accurately determine the desired segments. The Sobel gradient is employed to enhance features in the captured images, which aids in localizing the intrastromal cuts. Bayesian optimization is then used to fine-tune the hyper-parameters of the segmentation algorithm.

A pioneering work of Li et al. provides a framework for corneal interference segmentation to detect flap interfaces and the other anatomical corneal layers [43]. The proposed algorithm utilizes a transverse average filter followed by an intensity scan to identify local maxima. While average filtering demonstrates reliable performance, it can blur edges and fine details. To address this, we applied Non-Local means (NL Means) denoising, which effectively preserves details while reducing noise. Advances in computational image processing have made NL Means approaches the preferred method for denoising and maintaining fine structures. However, the choice of denoising strategy depends heavily on image registration and the technology applied.

Additionally, an active contour model (commonly referred to as a “snake”) was used to optimize boundary detection [44–46]. This model incorporates image gradients and a smoothness constraint to guide the contour towards sharp boundaries.

In this study, we initially applied a point-density scanning routine, which identifies the boundaries of flap or anterior lenticule cuts. Furthermore, a classifier and Bayesian optimization were employed to distinguish the boundaries of lenticule cuts. Both the “snake” model and Bayesian optimization minimize a cost function, with the strategy varying based on the specific problem. However, unlike the “snake” model, the proposed optimization does not require progressive filtering. On the other hand, our utilized optimization routine constantly assumes a smooth and non-singular underlying function.

As part of our deterministic study, we analyzed laser-treated, peri-operative cuts on porcine eyes to demonstrate the effectiveness of our method. Our work not only targets peri-operative characterizations, but also can be of use for postoperative residual monitoring.

This study is conducted in a laboratory setting, with a primary focus on validating the geometry of laser-generated corneal cuts prior to any surgical manipulation. By isolating laser-induced geometries from surgical effects, the performance of laser systems in a more controlled and precise manner can be evaluated.

We compared the results of our proposed method against manual measurements of identical substructures to validate its accuracy and reliability. In addition to matching the central measurements obtained by traditional approaches, our method demonstrated a significant advantage by enabling the characterization of key geometric properties, including lenticule power, transition zones, and incision angles. These characteristics could not be determined using previously developed methods that focused solely on central substructure measurements.

Furthermore, despite the ex vivo nature of the experiment, we observed a strong correlation between the intended design parameters and the automated measurements produced by our method. Such capabilities not only extend the scope of preoperative and peri-operative analysis but also pave the way for more comprehensive postoperative evaluations, ultimately contributing to improved surgical predictability and patient outcomes.

Our approach can potentially fill a current gap in precise reassessment of intrastromal laser cuts where manual expert measurements were seen as the conventional (user-based) standard. Our method can realize the integration of advanced fit routines or machine learning approaches for comprehensive live characterization.

We note that the primary objective of this work is to demonstrate and validate the generalizability and robustness of our proposed approach, rather than to derive specific clinical conclusions.

2 Materials and methods

2.1 Experiment

OCT images were obtained with a Thorlabs GAN111 OCT consisting of a Ganymede GAN111 as the base unit and an OCTG9 (with OCT-LK4-BB scan lens) as the scan head unit [47]. With a wet lab handling, the epithelium of each eye was first removed, and the eyes were placed under the laser for further processing [48]. The treatments were performed with the SCHWIND ATOS femtosecond laser system (SCHWIND eye-tech-solutions, Germany). The OCT images were collected from 83 ex vivo porcine eyes underwent laser treatments leading to either flap or lenticule creation. Due to the incomplete cuts or faint traces, 23 eyes were excluded from the study. The remaining cohort of ex vivo treated eyes was divided into two groups of 24 (lenticule) and 36 (flap) eyes. The intended substructure geometries for all eyes are summarized in Table 1.

For both flap and lenticule creations, a combination of key parameters was employed to achieve versatility through the generation of customized intrastromal substructures. The values were chosen to represent the full spectrum of possibilities, including the most extreme cases, for both flap and lenticule creation. The intended cap thickness for all lenticule substructures was set to 150 μm.

All animal specimens used in this study were procured as a by-product from a local approved slaughterhouse, and the handling of these specimens was not part of the study, since eyes were enucleated post-mortem. All the treatments were performed on ex vivo porcine eyes within five hours of procuring the eyes. The eyes were stored in Normal Saline solution at room temperature (≈23 °C), during transportation and before treatment ensuring cornea’s health and transparency. For all treatments, the eyes were mounted on a holder and de-epithelialized manually prior to treatments (as a standard routine) using an Amoils brush, to have consistent smooth corneal surfaces, which is important for creating exact substructures. Intraocular pressure was controlled using a specialized ex vivo eye holder to ensure consistent experimental conditions. These conditions could preserve optimal condition for corneal integrity and ensure proper standardization for laser cutting. The study did not involve live animal experimentation and used the eyes ex vivo as a suitable animal model.

Each eye then was docked with the disposable patient interface and held with a vacuum level of 250 mmHg. The treatment specifications were set before each treatment to minimize the handling time for the treatment and image acquisition. The porcine eyes were treated and, immediately (peri-operatively) after treatment, placed under the OCT instrument to capture images. Neither flap cuts were lifted nor lenticule volumes were extracted after imaging.

Subsequently, two independent users manually performed measurements using the Thorlabs OCT software, where they loaded images and employed the integrated measurement tool. In these measurements, the refractive index was set to 1.37 to ensure consistency. The users’ expertise ensured that the measurements were conducted with a high degree of precision.

Emphasis was placed on establishing a systematic process to ensure tissue homogeneity throughout the experiment, from preparation to OCT scans. This approach aimed to minimize variability and ensure consistent and reliable results across all stages of the procedure.

Representative images from each category of treatments are provided in Figure 1. All scans (raw data) were imported using ThorImageOCT – Version: 5.4.8 – software into our Python routine for characterization.

Figure 1 (a, b) Illustrate the aerial and OCT scan of the ex vivo intrastromal flap cut and (c, d) the lenticule cuts with SCHWIND ATOS on porcine eyes. Note the edge cut of the flap delineates as a circle in (a) and the reflection of plasma bubble layers (anterior and posterior) can be seen in (c).

The OCT GAN111 is a configurable high-resolution spectral-domain OCT (SD-OCT) system equipped with an 880 nm center wavelength [47]. It supports A-scan rates of up to 20 kHz, enabling fast and precise imaging. The system offers imaging depths in air(eye) of up to 3.4 mm (2.5 mm), with axial resolutions of 6.0 μm in air (4.4 μm in eye). Its sensitivities range from 96 dB (at 20 kHz) to 106 dB (at 1.5 kHz), allowing for the cornea and anterior imaging (Table 2).

We used a corneal refractive index of 1.37 and the correction model proposed by [49] to justify thickness measurements.

The effects of optical refraction on the OCT images might affect lateral measurements primarily by overestimating lateral extents [50]. However, due to the shallow cuts and the naturally flat shape of porcine corneas, optical refraction-induced distortion was minimized. Based on Snell’s law, the estimated lateral distortion due to refraction at a depth of 150 μm is predicted to be less than 1%, as discussed in [51]. This small percentage reflects a very minimal distortion in the OCT image, making this effect negligible in our practical application.

2.2 Image processing

Given a discrete gray scale and noisy image, we utilized the NL Means algorithm for denoising [52–54]. The intensity of a pixel at position q in the noisy image thus is given,$$I\left(p\right)=\sum _p\mathrm{}W(p,q) I\left(q\right),$$(1)where W(p, q) is a non-local weighting function as follows,$$\begin{array}{l}W(p,q)=\frac{1}{\sum \exp \left(L_{\mathrm{HUB}}\left[-\frac{E\left(p,q\right)}{h^2}\right]\right)} \exp \left(L_{\mathrm{HUB}}\left[-\frac{E(p,q)}{h^2}\right]\right),\\ E(p,q)=\left|\right|P_p,P_q\left|\right|.\end{array}$$(2)

Here E(p, q) is the Euclidean distance of pixel patches centered at p and q · h² is a smoothing parameter, which controls the sensitivity of the weighting function to differences in E(p, q).

To measure how different two patches are, W(p, q) calculates their distance. Further, W(p, q) incorporates the Huber loss as the robust penalty operator like the conventional L₁ and L₂ norms. The Huber loss is given,$$L_{\mathrm{HUB}}\left(f\right)=\{\begin{array}{ll}\frac{1}{2}·f^2& \in \left|f\right|\le \delta \\ \delta ·\left(\right|f|-\frac{1}{2}·f)& \mathrm{otherwise}\\ & \end{array}$$(3)where f is the residual and δ is a threshold parameter that determines the switch between the L₂ and L₁ norms [55, 56]. The L₂ and L₁ terms are two commonly used metrics to quantify differences between two entities, such as image patches in this context [57]. The Huber loss function transitions smoothly between the quadratic behavior of L₂ norms and the linear behavior of L₁ norms, depending on the value of the input f:

• For |f| ≤ δ, the quadratic term behaves like the L₂ norm and penalizes small differences strongly.

• Otherwise, the linear term behaves like L₁ norm and provide robustness to large outliers.

By combining these two approaches, the Huber loss provides a hybrid solution. It behaves like L₂ for small differences, preserving fine details, while transitioning to L₁ for larger differences, ensuring robustness against noise and outliers. This balance is particularly advantageous in image processing tasks, where maintaining important features while reducing noise is critical.

The estimation of δ mostly relies on the noise standard deviation of images (σ) commonly set to a range of 1.5 σ–2.5 σ, a heuristic that has been followed in our work. Thus, given average σ = 5 pixels throughout the collected dataset, δ = 10 pixels was set. The optimal value of h was empirically determined to be 10 which ensures balancing noise reduction and structural integrity.

The Sobel operator is used to find the approximate gradient of the image intensity function [58, 59]. It consists of two convolution kernels (or masks) that estimate the gradients in the horizontal (x) and vertical (y) directions. Given an input image I, the gradients in the x and y directions (denoted as I _x and I _y ) are calculated using the convolution of the image with the Sobel kernels,$$I_i=G_i\otimes I : i=\left(x,y\right).$$(4)

The maximum amplification of gradients through OCT images can be obtained with the Sobel’s kernel (G) in the horizontal and vertical directions as follows,$$G_x,G_y=\left|\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right|,\left|\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right|$$(5)where the combination resembles Robert Cross operator [60].

The above chain of processing leads to collect a 2D cloud of points encompassing the cornea, and intrastromal cuts. A vertical peak finder algorithm was used to specify points belonging to the topmost intensified feature that forms the cornea [61].

Additionally, fitting a parabolic polynomial classifier (PPC) to the entire cloud of points distinctly splits the cloud of points such that,$$\{\begin{array}{ll}{x}_i\in \mathrm{Corneal} \mathrm{Segment}& y_i-\mathrm{PPC}\left(x_i\right)>ϵ\\ x_i\in \mathrm{Intrastromal} \mathrm{Cuts}& y_i-\mathrm{PPC}\left(x_i\right)<ϵ\end{array}$$(6) ϵ is a threshold to assure that points are in proximity to the desired regions. The PPC aims to minimize the fitting error, ultimately leading to a fit positioned between the cornea and the intrastromal structure. The epsilon parameter, therefore, helps determine the layer to which each peak belongs. ϵ was considered to be zero for splitting two layers, whether the cornea-intrastromal or anterior-posterior, as vertical upper and lower layers.

The R² score was used to identify whether the intrastromal points would determine a flap or lenticule cut. The R² of the PPC in flap cuts reaches a higher value compared to lenticule cuts. A lenticule substructure typically exhibits more vertical peaks at lateral positions compared to the flap. Therefore, when fitting a PPC through additional corneal peaks, the resulting model exhibits higher variance. Consequently, the R² score for lenticules is typically at least 10% lower than that for flap substructures.

A similar methodology can be applied to distinguish the anterior and posterior segments of a lenticule cut. However, the PPC approach is very sensitive to sparse noises, which refers to occasional outliers that occur at isolated points within the peak data. It may lead to extended flap or posterior segments since the classifier may cross the sidebands of anterior. Thus, applying the equation 6 would not solely suffice for a precise segmentation. Hereafter, the sideband noise (or sideband falsely detected peaks) refers to peaks that appear outside the boundaries of the lenticule or flap extension. These falsely recognized peaks can lead to overestimation of the diameter.

For flap cases, a solution to filter out falsely recognized points involves scanning the point density starting from the center and extending toward the sidebands. As the scan progresses, the point density reaches a threshold where it no longer exhibits significant variation compared to the earlier steps. This method effectively identifies and excludes false peaks detected in the sidebands, improving the accuracy of the layer recognition. This approach can be seen as a density-point scan routine for refining flap boundary identification.

For posterior segments of lenticules, a Bayesian Optimization (BO) approach was implemented [62–64]. Empirically, it was observed that the point-density scanning routine, while effective for flap cases, is not yet sufficiently accurate to handle lenticule structures.

A mean-squared-error (MSE) was set as the cost function to be minimized by shortening (if needed) the lateral extent of posterior segments. We have utilized GP-Minimized package from Scikit-optimize (open-source) library which performs sequential model-based optimization. GP-Minimized algorithm has certain advantages. However, it should be noted that GP-Minimized assumes a smooth and non-singular underlying function [65].

The pseudo-algorithm to filter out noises from posterior cuts can be summarized as Algorithm 4 where (x_min, x_max) specifies the range obtained by the initial PPC, α is the set of weights for the polynomial fit on the posterior segment, and T is the maximum number of iterations, thereby the optimization reaches a monotonic behavior. The algorithm was initialized with α = 0 and N_max = 12. The initialization of α was made to provide a neutral starting condition, free from any prior assumptions or bias that could influence the optimization process. Additionally, it was anticipated that all layers could be accurately described by a 12th-order polynomial fit. This assumption was based on the need for a model that could capture the complex curvature and variations within the corneal and substructural layers.

Algorithm 1 iteratively filters out points outside a boundary while fitting polynomials. It works by proposing new candidate parameters, such as the bounds of the data and the polynomial fit order, fitting a model to the selected data points, and checking how well the model performs. The performance is evaluated using an objective function, typically the MSE, which is minimized during the process. Based on the performance, the algorithm updates its search strategy to focus on areas that are likely to lead to a better fit, refining the boundary to exclude irrelevant points. This process is repeated several times, continuously improving the fit by selecting the most relevant data points. Thus, the sideband noises are effectively filtered out, allowing the algorithm to accurately determine the true(x_min, x_max) as the boundaries of the posterior layers. After determining the boundaries of a cut, a cord line connects both ends. This line can be seen as the possible tilt reference.

Following the application of the PPC and Bayesian Optimization, the corneal and sub-structural layers are identified – a morphology as depicted in Figure 2 – which enables the definition of key geometries for further characterization. This foundational step allows for the automatic measurement of local thickness, which is essential for understanding the structural properties of the cornea.

Figure 2 (Left) Flap and (Right) Lenticule geometry to characterize as equation 7. FT and FD stand for the flap thickness and dimeter while CT, LT, CD and LD refer to the cap (lenticule) thickness and dimeter. The bed and hing cuts (for flaps and lenticules) are shown as the gray lines.

The Euclidean distance between two corresponding points at two segments estimates the local thickness. Those pairs of points are found as the cross sections between a line, normal to the cord, that crosses two segments. Thus, the following thicknesses can be determined (see Figure 2),$$\{\begin{array}{ll}\mathrm{FT}=\left|\right|C(x_c,y_c)-A(x_a,y_a)\left|\right|& \in \mathrm{flap}\\ \mathrm{LT}=\left|\right|A(x_a,y_a)-P(x_p,y_p)\left|\right|, \mathrm{CT}=\left|\right|A(x_a,y_a)-C(x_c,y_c)\left|\right|& \in \mathrm{lenticule}\end{array}$$(7)where FT, CT, and LT refer to the flap, cap, and lenticule thickness.

Due to the variation around center, all LT-CT-FT values are calculated as an average within close proximity with the center (≈10 pixels) and represented as Average ± StdDev. A positioning standard deviation is also to be considered for reporting thickness measurements.

The lenticule and flap diameters (LD, FD) are measured as the extent of the left and right boundaries (i.e. the cord length).

Figure 3 illustrates the impact of our post-processing approach. It represents the post-processed version of Figure 1, obtained by applying the steps outlined in the Materials section, which are described in detail below.

Figure 3 (Left) – Raw images of Figure 1. (Right) – The proposed image processing pipeline segments intrastromal lenticule and flap cuts by first enhancing the structures (a, b), then applying a segmentation algorithm to identify corneas (c, d) and consequently the domains of cuts (e, f). Panels in (e) illustrate the presence of sideband noise and its removal to achieve the correct flap geometry. Panel (f) demonstrates the optimization process for detecting the correct posterior (lenticule) cut (orange dots) using Algorithm 1.

To ensure optimal performance, the initial parameters of our algorithms were determined via a preliminary tuning. This trial involves a 5% subset of the dataset from both lenticule and flap intrastromal cuts. These parameters were set to maximize performance across the entire dataset. Subsequent to fine-tuning on this subset, the algorithm operated unsupervised on the complete dataset. Thereby, the output variability is minimized. All parameters were derived from this initial optimization phase. As is customary in pre-tuning for unsupervised tasks, this subset was selected to optimize algorithm performance. Notably, the algorithm processes the remaining data in its entirety, unless incomplete cuts or faint traces are detected.

Furthermore, the proposed approach enables the determination of the optical zone. Characterizing transition zones would determine the quality of the created corneal interface. The design of the transition zone can influence residual aberrations, such as coma and spherical aberration [66]. Numerically, the transition zone is determined as the region where the posterior curvature shifts from convex to concave. This definition corresponds to a sign change in the second-order gradient of the posterior fits (P(x _i , y _i ) as shown in Fig. 2).

The same rationale can be applied to flap cuts for determining edge cut angles. However, this method necessitates continuous and complete flap cuts to ensure a high-quality fit. To estimate the edge cut angle, peaks within a certain window surrounding a boundary can be selected. Then, two vectors pointing from the boundary point towards the most left ($\vec{A}$) and right ($\vec{B}$) peaks are formed. Thus, an inner product determines the inner angle as follows,$$\theta =\mathrm{arcosin}\left(\frac{\vec{A}·\vec{B}}{\left|A\right|\left|B\right|}\right),$$(8)where || is the norm of vectors. The deviation may root back to monomial fits that always retrieve a smoother edge compared to practical (planned) sharp edges. Thus, a widening is to account for.

Moreover, the proposed characterization approach enables scanning the lenticule thickness profile. The thickness profile can be described with a parabolic equation which determines the lenticule extraction curvature (R).

Given the refractive indices of air and porcine eyes, the lenticule power (LP) formulates as$$\mathrm{LP}=\left(\delta n\right)·\frac{1}{R}.$$(9)

The outlined steps form the proposed numerical pipeline for characterizing intrastromal substructures with a high level of morphological detail.

3 Results

Figures 3a and 3b illustrate the denoising of the OCT scans (corresponding to Fig. 1). By initial use of the peakfinder algorithm, the corneal segment was recognized as shown in Figures 3c and 3d. For flap cases, the intrastromal peaks identified with the PPC may often lead to wider boundary detections. Applying the point-density routine would exclude the sideband peaks as illustrated in Figure 3e.

In Figure 3f, through approximately 8 iterations of Algorithm 1, the optimized posterior (orange line) was obtained. The MSE of Algorithm 1 typically decreases by an order of magnitude after the first four iterations and then continues to decrease monotonically. Additionally, the same point-density scan routine, as Figure 3e, was applied to the lenticule cuts for filtering out falsely detected intrastromal peaks. The point-density scan determines the extent of anterior cuts but may also extend the posterior boundaries. Consequently, the BO (Bayesian optimization) is required to accurately adjust the width of the posterior cuts. Generally, controlled optimization algorithms such as Gradient Descent or RMSprop [67] which rely on the gradient of the objective function can be utilized. However, the performance analysis can form the scope of future work.

The optimized range for peak prominence was determined by implementing a feedback loop that monitors the R² score for each segmentation. Typically, the fine-tuning converges more rapidly for corneal segments compared to intrastromal cuts. An example of the prominence range setting is illustrated in Figure 4a, which is applied to the data shown in Figure 3e. Figure 4b demonstrates the MSE drop over iterations.

Figure 4 (a) Peak prominence tuning and MSE drop illustration (b).

In Figures 5f and 5e, the segmentation of Figures 3a and 3b have been illustrated. In both cases, the intended flap, cap, and lenticule thicknesses were set to values higher than 100 μm. Manual characterization (segmentation) of thicknesses below than 100 μm can be particularly challenging. The challenges originate from the uncertainties such as the lack of corresponding points at the anterior or posterior segments.

Figure 5 (a–d) Raw OCT images and (e–h) their corresponding segmentations. (a, b) and (e, f) display thick cuts, including the flap and lenticule, whereas (c, d) and (g, h) represent thin cuts.

The proposed approach can systematically, faster and reproducibly, resolve these problems as shown in Figures 5g and 5h. The intended FTs in Figures 5f and 5h were 115 μm and 80 μm which were estimated as 120.4 μm and 81.9 μm. The treatment parameters in Figures 5e and 5g are (Sphere = −10D, Cylinder = −6D) and (Sphere = −2D, Cylinder = −2D) with the intended LTs of 157 μm and 54 μm. The LTs were estimated 155.6 μm and 59.9 μm. The lenticule, cap, and flap thicknesses were calculated as the average within a central region of 10 pixels. This approach ensures a reliable estimation by smoothing out local variations and potential noise in the calculations.

For thicknesses below 100 μm, larger deviations from the intended values were observed in the dataset. This discrepancy arises because both anterior and posterior intrastromal cuts have a certain width, which affects the performance of the peakfinder. Consequently, the estimated thickness may fluctuate by up to the average cut width, leading to variations in measurements for thinner thicknesses (Fig. 6). The two white lines per layer represent the boundaries of the laser tissue disruption, i.e. the vertical extent of tissue damaged by the irradiation. This disruption occurs where laser energy exceeds the breakdown threshold of the tissue. Bubbles are caused by localized plasma formation due to high-intensity laser pulses, leading to tissue disruption. The breakdown process scatters light significantly, making these regions appear as bright spots on the OCT images.

Figure 6 Illustrating intrastromal cuts width within thickness measurements may vary.

As a result, whether measurements are obtained manually or automatically, a positioning error of approximately 10 μm should be considered when comparing the measurements to the intended (thickness) values.

This positioning error adds up to StdDev calculated to present thickness values in proximity with the center. This error margin may vary with pulse energies exceeding 100 nJ.

Considering the proposed pipeline, the automated results align with the gold standard manual measurements, given correlation coefficients greater than 92% as shown in Figure 7. This validation, crucial for establishing the reliability of our approach, highlights the automated system’s ability to replicate manual measurements while offering significant advantages. Notably, automated measurements exhibit significantly less variation compared to manual measurements, suggesting improved precision and reduced inter-observer variability.

Figure 7 Direct comparison of the automated results to the manual results. (a, b) illustrate the comparison for the flap and (c) is an instance of the lenticule morphology.

Furthermore, our approach offers additional insights beyond basic measurements. For instance, it enables the analysis of lenticule power (LP), as well as the optical transition zone and the incision angles.

The LPs were estimated and compared with the intended LPs in Figure 8a. With a R² = 0.96, the values signify a correlation. However, the mild corrections (−2D and −4D) deviate more compared to the higher corrections. This observation may route back to the overall systematic errors associated with the tiny lenticule estimations. In Figure 8b, the powers associate with the larger planned diameters, and large powers usually correspond to the treatment planned with smaller diameters.

Figure 8 (a) illustrate the correlation between the intended and estimated LPs. (b) displays the different geometries (23 data points) versus the intended LPs annotated with the estimated LPs. Myopic corrections were intended. Minor corrections associate with the larger planned diameters, and vice versa.

All intended treatments were planned with corrections for myopic or compound myopic astigmatism. Without a fully automated characterization approach, inferential analyses may not be achievable as addressed in [15].

The estimated LPs in Figure 8a determine approximately 0.5D less correction than intended. The observed discrepancy might be attributed to the diameter shortening by the BO routine, the ex vivo optical breakdown and measurements on porcine eyes [49].

Overall, despite the influence of the discussed factors, the estimated values still remain reasonably close to the intended. This demonstrates that the observed discrepancies are within an acceptable range.

The optimization algorithm, combined with higher-order Taylor’s monomials, allows for a comprehensive capture of the posterior segments of lenticule cuts. This capability is particularly valuable as it enables the differentiation between the optical zone and the transition zone.

An example of transition zone determination has been shown in Figure 8b. On average, the transition zone could be determined in 8 lenticule cases as 382 μm ± 15 μm versus the expected 400 μm. This small subset corresponds to the cases with strong astigmatism corrections which require the greatest bend in the transition zone.

To accurately characterize possible transition zones, Taylor’s monomials of at least the 10th order should be fitted. Otherwise, the true shape of posterior cuts may not be properly reconstructed.

Incision angles for flap cuts can be determined using the same method. However, it is crucial to make continuous and complete flap cuts to ensure a precise fit.

An example of edge detection is shown in Figure 9c. Due to the small width of edge cuts, angle estimations may have an error margin of up to 10°. For a subgroup with an intended edge cut of 135°, the algorithm estimated an angle of 150° across 10 flap cuts, resulting in a total deviation of 15°.

Figure 9 (b) identifies the transition zones of (a), as the regions within the posterior fit curvature alters from convex to concave corresponding to a sign change of the 2nd order gradient. (c) illustrate the hinge angle detection (on flap cuts) using a similar rationale. The black dashed lines illustrate $\vec{A}$ and $\vec{B}$ to determine θ.

Moreover, full 2D segmentations would facilitate the application of 3D analysis such as determination of Keratometry reading or analysis of anterior aberrations [68, 69].

Despite the potential influence of factors such as hydration on the peri-operative eyes, the analysis demonstrated a high correlation between intended and automatically measured values. It might suggest that ambient conditions did not significantly alter the morphologies.

Figure 10 provides a summary of the linear regression analysis for the estimations of FD, FT, LD, and LT. The proposed approach estimates flap and lenticule diameters, as well as thicknesses, without additional computational overhead. The regression analysis indicates that the estimations on flap cuts show a stronger correlation with their intended values compared to lenticule cuts. Specifically, the estimated flap thicknesses (FTs) exhibit the smallest discrepancies from the intended values. The lenticule characteristics generally demonstrate greater deviations. Notably, the highest deviation is observed in the subgroup with an intended lenticule thickness (LT) of 150 μm. This deviation from regression might originate from the statistically insufficient sample size. A larger sample size may better determine if measured lenticule thicknesses will be consistently larger than intended values. It was found that the measured LT was significantly greater than the intended LT. This discrepancy was observed in the overall scheme rather than a specific thickness range [15].

Figure 10 Regression representation comparing intended lenticule-flap diameter and thickness with the intended values. Thicknesses are in μm and diameters in mm.

The estimated LD values are approximately 0.8 mm smaller than the intended values (Fig. 10c). These discrepancies can be partially attributed to the BO that numerically compresses the posterior peaks. In the absence of an optimization routine, these values could significantly deviate towards much larger values. Thus, a scope of future work would integrate the performance of versatile optimization approaches to figure out the best combination for posterior characterization. Additionally, associated difficulties to determine transition zone within LD can potentially contribute to the observed discrepancies. Moreover, the estimations (or even manual measurements) were performed ex vivo on porcine eyes unlike human living subjects which also can trigger shorter estimations than planned [49]. Complementarily, Figure 10a illustrates that the estimated values of FD lag behind the intended values, approximately 0.5 mm. The discrepancies can root to a probable combination of edge compression by the BO and ex vivo measurements.

For all panels in Figure 10, the regression analysis shows R² > 0.96 and p < 0.001 which indicates a strong statistical relationship between the estimated and intended values.

The intrastromal cuts were scanned across four meridians with a 45° angular step size (referred to as bunches). All scans from the bunches were processed by the algorithm for characterizations. No significant fluctuations in thickness or diameter were estimated within these bunches. For each quantity (such as FT, FD, LT, and others) within each bunch, a StdDev was calculated. In Table 3, we report the maximum StdDev observed for each quantity across all bunches. The highest variations in estimates, although still minor, were observed in bunches with larger diameters or thinner thicknesses. It can characterize the highest level of variability observed for each particular quantity among all the bunches. The repeatability of the overall characterization pipeline (estimated for the maximum within bunch StdDev values) results in a 95% confidence interval of ±200 μm for determining diameters.

Since the algorithm is deterministic, multiple runs of the algorithm on the same scan always produce an identical segmentation (Dice similarity of 1 [70]). In contrast, the manual measurements performed by independent (skilled) users often struggled to converge to a unique solution. This discrepancy is particularly evident in thickness measurements.

As shown in Figure 11, two skilled users performed manual measurements on a subset of data. The manual measurements revealed a higher deviation in LT compared to FT, due to the smaller thicknesses involved. Overall, comparing the average of linear regressions between Figures 10b and 10d and Figures 11a and 11b reflect the same underlying trend. Despite these challenges, the manual measurements and algorithmic characterizations are correlating, with a positioning error of approximately 10 μm.

Figure 11 User based measurements of flap and lenticule thicknesses. Significant variation on lenticule based measurements were observed.

4 Discussion

While Canny and Hildreth methods could potentially outperform the Sobel operator, they require significant tuning [71]. In particular, the Canny method, with additional adjustments, can segment both ends of a single cut. Nevertheless, for our initial approach, this level of precision was not critical. The Canny edge detection method, with its Non-Maxima suppression, produces smoother edges. However, the rougher edges detected by the Sobel operator were sufficient for the requirements of this particular application. Future developments may necessitate more sophisticated modeling to capture finer substructures. Further advancements may, nonetheless, necessitate the incorporation of more sophisticated modeling to resolve finer substructures.

One limitation of our approach is the reliance on conventional heuristics for setting parameters such as δ and h. This legitimacy directly impacts the noise removal algorithm and, consequently, the performance of the Sobel gradient. The current method involves fine-tuning of these parameters using a small set of data. While this approach has not shown significant variation across the entire dataset, it limits the effectiveness of patch-specific fine-tuning. To streamline the image processing pipeline, a standardized capture routine could be beneficial. Computationally, more sophisticated techniques such as Grid and Elastic search could be explored for parameter tuning [72, 73]. These methods could potentially reduce the manual effort involved in parameter optimization and improve the overall robustness of the algorithm.

For all collected datasets of porcine eye the peak prominence range was identical. No case showed a fluctuating R² versus prominence or a monotonic behavior since the start. The BO also was of necessity for all cases, though the drop down of MSE might occur earlier than four iterations.

Three major parameters conventionally specify the quality of the flap and lenticule intrastromal cuts as the spot distance, track distance and the pulse energy (nJ). Certain combinations of these parameters would lead to improve intrastromal cuts. However, incomplete cuts are probable to observe (Fig. 12). Those intrastromal patterns (as Fig. 12a) lead to low values of R² score (of fit) that rejects proceed of segmentation. In contrary, although the flap in Figure 12b was not performed perfectly, the R² score still gains a value higher than 95%. A shorter flap diameter compared to the planned geometry could qualify the performance of the treatment system on creating intrastromal cuts.

Figure 12 (a, b) Examples of instrastromal flap (incomplete) cut at 75 nJ laser pulse energies.

In our model, a refractive index of 1.37 was employed for vertical measurements, which is a commonly accepted approximation. Nevertheless, it is important to acknowledge that the accuracy of these measurements is influenced by the refractive index, which can vary substantially with the state of the corneal hydration. Since these are ex vivo measurements where tissue hydration is not strictly controlled, a more precise approach could involve performing the lenticule extraction (or flap cut) followed by direct measurement with a micro-meter. This consideration will be addressed in future studies to further enhance accuracy.

Through analysis, the error propagation was generally low across all the parameters which may indicate that the uncertainties introduced by individual measurements had a minimal impact on the final results. Specifically, for parameters such as the flap thickness, the low propagation of error suggests that the precision of the measurement system is well-maintained throughout the process.

For instance, when calculating the axial thickness, we found that the uncertainty from the resolution of boundary measurements propagated with only a small increase in overall uncertainty, as follows.

Considering Table 3, the FT measurement determines a σ value of 1.93 which results in approximately ±3.6 μm uncertainty for axial thicknesses.

The error propagation of uncertainties, when adding independent measurements, is given as follows,$${\sigma }_{\mathrm{results}}=\sum \mathrm{}{\sigma }_i^2.$$(10)

Here, for two boundaries with equal resolution uncertainty (σ _r ), the axial thickness uncertainty (σ _t ) becomes,$${\sigma }_t=\sqrt{{\sigma }_r^2+{\sigma }_r^2}=\sqrt{2}{\sigma }_r.$$(11)

Substituting σ _r  = 2.46 μm, equation 11 yields approximately 3.47 μm which matches the expected uncertainty of the nominal resolution.

The primary objective of this work was to demonstrate the applicability of our proposed method. However, its feasibility and generalizability can potentially provide deeper insight into new measurements via another diagnostic device, incorporating another femtosecond laser and a different eye model. In light of that, we envisage further work that incorporates the SCHWIND MS-39 OCT system [74], the goat eye model and two different femtosecond lasers (SCHWIND ATOS and VisuMax 500). This investigation further supported our aim of developing a generalized approach with consistent performance across various combinations [68].

For instance, with the MS-39 instrument, each eye can be scanned across 12 full meridians. The application of multi-meridian scans (with a minimum of 3) would enhance the capability to fit Zernike polynomials or the Gatinel-Malet polynomials to the detected corneal heights [69, 75–77]. This advancement could significantly improve the determination of corneal properties, including wavefront representation and aberration analysis. While Taylor’s monomials are effective for describing 2D slices, they do not support comprehensive inferences related to aberration or refraction, which are critical for advanced diagnostic and corrective applications. Integrating these polynomial fits could thus offer a more detailed and accurate assessment of corneal topography and optical performance.

Notably, our proposed method shows promise not only as a standalone solution but also as an annotation tool for training machine learning segmentation architectures, such as Image-to-Image translation models [78]. As part of future research, we envision integrating knowledge of corneal layer sequencing into our model, as outlined in [43].

These findings have direct applications in the design and calibration of laser systems. Also, it facilitates assessing the repeatability and reproducibility (accuracy and precision) of laser-generated corneal substructures.

5 Conclusion

Intrastromal peri-operative (intra-OP) segmentations allow for precise measurement and visualization of the cuts made by the laser. This process ensures that the intended refractive corrections are accurately implemented.

We proposed a generalized routine to accomplish intrastromal flap and lenticule segmentation with an accuracy comparable with the manual expert level segmentations. Our method also characterizes the geometry of peri-operative intrastromal flap and lenticule cuts. This is particularly valuable for ultra-thin flaps and lenticules, where manual measurements can be highly variable and challenging. Additionally, our approach exceeds the accuracy of manual measurements through complex characterizations where user-based marking can be significantly doubted. Specifically, for the purpose of evaluational (comparison) measurements, through different refractive correction systems, an automated approach would essentially substitute manual measurements.

The proposed computational algorithm significantly accelerates characterization tasks. It reduces the processing time from several minutes to just a few seconds and maintains the accuracy and solution uniqueness. After fine-tuning the parameters over a small batch, the algorithm runs unsupervised across the dataset.

Our approach can facilitate feedback assessment routines to evaluate the peri-operative situations. It also can be used during (or ex-situ of) a refractive surgery treatment to reflect the influence of versatile parameters on the cuts’ quality.

Our effort was concentrated on demonstrating the potential of the proposed work, particularly its ability to achieve comprehensive intrastromal characterization. This capability, to the best of our knowledge, has been showcased for the first time. Consequently, statistical analysis of subgroups was not prioritized. Our primary goal was not to derive specific clinical findings or conclusions, but rather to establish the generalizability and robustness of our approach.

Acknowledgments

Dr. Shwetabh Verma from SCHWIND eye-tech-solutions and Francesco Versaci from CSO, Italy are acknowledged for the handling and preparing data.

Funding

This research received no external funding.

Conflicts of internet

The authors have nothing to disclose.

Data availability statement

The source code supporting the findings of this study are available from the corresponding author upon reasonable request. Data supporting the findings of this article are not publicly available due to clinical restrictions.

Author contribution statement

Conceptualization, S.A.M. and K.D.; Methodology, M.M.; Software, M.M., F.B; Experiment, A.P. L.K; Validation, M.M., S.A.M; Investigation, M.M., A.P.; Data Curation, M.M.; Writing – Original Draft Preparation, M.M.; Writing – Review Editing, M.M., S.A.M.; Visualization, M.M.; Supervision, K.D. and S.A.M; Project Administration, A.P..

References

All Tables

All Figures

	Figure 1 (a, b) Illustrate the aerial and OCT scan of the ex vivo intrastromal flap cut and (c, d) the lenticule cuts with SCHWIND ATOS on porcine eyes. Note the edge cut of the flap delineates as a circle in (a) and the reflection of plasma bubble layers (anterior and posterior) can be seen in (c).
In the text

	Figure 2 (Left) Flap and (Right) Lenticule geometry to characterize as equation 7. FT and FD stand for the flap thickness and dimeter while CT, LT, CD and LD refer to the cap (lenticule) thickness and dimeter. The bed and hing cuts (for flaps and lenticules) are shown as the gray lines.
In the text

Figure 3 (Left) – Raw images of Figure 1. (Right) – The proposed image processing pipeline segments intrastromal lenticule and flap cuts by first enhancing the structures (a, b), then applying a segmentation algorithm to identify corneas (c, d) and consequently the domains of cuts (e, f). Panels in (e) illustrate the presence of sideband noise and its removal to achieve the correct flap geometry. Panel (f) demonstrates the optimization process for detecting the correct posterior (lenticule) cut (orange dots) using Algorithm 1.

In the text

	Figure 4 (a) Peak prominence tuning and MSE drop illustration (b).
In the text

	Figure 5 (a–d) Raw OCT images and (e–h) their corresponding segmentations. (a, b) and (e, f) display thick cuts, including the flap and lenticule, whereas (c, d) and (g, h) represent thin cuts.
In the text

	Figure 6 Illustrating intrastromal cuts width within thickness measurements may vary.
In the text

	Figure 7 Direct comparison of the automated results to the manual results. (a, b) illustrate the comparison for the flap and (c) is an instance of the lenticule morphology.
In the text

	Figure 8 (a) illustrate the correlation between the intended and estimated LPs. (b) displays the different geometries (23 data points) versus the intended LPs annotated with the estimated LPs. Myopic corrections were intended. Minor corrections associate with the larger planned diameters, and vice versa.
In the text

	Figure 9 (b) identifies the transition zones of (a), as the regions within the posterior fit curvature alters from convex to concave corresponding to a sign change of the 2nd order gradient. (c) illustrate the hinge angle detection (on flap cuts) using a similar rationale. The black dashed lines illustrate $\vec{A}$ and $\vec{B}$ to determine θ.
In the text

	Figure 10 Regression representation comparing intended lenticule-flap diameter and thickness with the intended values. Thicknesses are in μm and diameters in mm.
In the text

	Figure 11 User based measurements of flap and lenticule thicknesses. Significant variation on lenticule based measurements were observed.
In the text

	Figure 12 (a, b) Examples of instrastromal flap (incomplete) cut at 75 nJ laser pulse energies.
In the text