Automatic estimation of the cross‐sectional area of the waist of the nerve fibre layer at the optic nerve head

Glaucoma leads to pathological loss of axons in the retinal nerve fibre layer at the optic nerve head (ONH). This study aimed to develop a strategy for the estimation of the cross‐sectional area of the axons in the ONH. Furthermore, improving the estimation of the thickness of the nerve fibre layer, as compared to a method previously published by us.


| I N T RODUC T ION
In glaucoma, the ganglion cell axons forming the nerve fibre layer of the retina gradually degenerate, causing loss of vision.Ophthalmologists evaluate the capacity of the visual system to transfer signals to the brain in glaucoma qualitatively by observing the structure of the optic nerve head and quantitatively by observing the loss of the visual field.Segmentation of the optic nerve head (ONH) in OCT captures allows objective quantification of the retinal nerve fibre layer.There are suggestions to measure the thickness of the nerve fibre layer with OCT around a circle close to the ONH (Greaney et al., 2002;Jaffe & Caprioli, 2004) integrated over a surface close to the ONH (Hwang et al., 2016), or alternatively as the waist of the nerve fibre layer at the ONH measured as average thickness or cross-sectional area integrated over an angle in the facial plane (Považay et al., 2007).Považay et al. (2007) suggested defining the cross-section of the waist of the nerve fibre layer at the ONH as the minimal distance between the ONH Pigment epithelium Central Limit (OPCL) and the Inner limit of the Retina Closest Point (IRCP; Figure 1).
The minimal distance was suggested by Považay et al. to avoid tilt angle dependence of the measured quantity.
Typically, OPCL and IRCP are annotated discretely in the captured OCT volume.Algorithms for automation of the annotation of OPCL and IRCP have been published (Belghith et al., 2014;Miri et al., 2017).We developed a semi-automatic algorithm for 3D scans of the ONH requiring manual annotation of OPCL and automatic detection of IRCP (Sandberg Melin et al., 2019a).Later, we automated the annotation of OPCL in 500 radial representations of the volume at equidistant angles in the frontal plane, as well as the 3D representation of the inner limit of the retina, using deep learning algorithms (Brusini et al., 2021;Carrizo et al., 2020).For each OPCL, the corresponding closest point on the 3D representation of the inner limit of the retina, IRCP, is automatically detected.The distance between OPCL and IRCP, the Pigment epithelium central limit-Inner limit of the retina Minimal Distance (PIMD; Figure 1) is thus resolved in 500 angles in the frontal plane.PIMD at each angle α was defined as PIMDα.Our previously published strategy for PIMD estimation (Sandberg Melin et al., 2019a) will be referred to as PIMDv1 in the current paper.
The commercially available algorithm (Spectralis, Heidelberg Engineering, Germany) for estimating the waist of the nerve fibre layer as thickness (Bruch's Membrane Opening Minimal Rim Width, BMO-MRW [Chauhan et al., 2013]) or alternatively as cross-sectional area (Bruch's Membrane Opening Minimal Rim Area, BMO-MRA [Gardiner et al., 2014]) integrated over an angular segment from radial scans of the ONH is described elsewhere (An et al., 2021;Enders et al., 2016;Gardiner et al., 2014).The cross-sectional area was estimated as the sum of surface elements on a conical frustum.
On the assumption that healthy human eyes have roughly a similar amount of ganglion cell axons and that in the vicinity of the ONH, the ganglion cell axons are of even diameter, it is expected that the average thickness of the waist is smaller in large ONHs (Považay et al., 2007).This is because the ganglion cell axons will be distributed over a larger surface in a large ONH.On the contrary, the total area is expected to be independent of ONH size.Hence, estimates of cross-sectional thickness in subjects are expected to vary more than estimates of cross-sectional area.Consequently, the limits for a non-pathological ONH will be wider for cross-sectional thickness than for cross-sectional area.This would render measurements of cross-sectional thickness less sensitive for the detection of glaucomatous change.Indeed, the estimated cross-sectional thickness was found larger in micro ONHs than in macro ONHs, while the crosssectional area was independent of ONH size (Enders et al., 2016).
The main purpose of the current paper was to develop an algorithm for the estimation of the cross-sectional area of the waist of the nerve fibre layer at the ONH from a voxel representation of the ONH, considering local deviations of radius in the frontal plane.In addition, a new strategy for estimating the thickness of the waist of the nerve fibre layer at the ONH was planned to be compared to our previously published strategy (Sandberg Melin et al., 2019a).

| Subjects
Subjects recorded in the present study are young adults without eye disease, aged [20;30] years with equal gender distribution and a spherical equivalent refractive error within [−5;+5] D. Ethical approval with diary no.2019-05407 was obtained from Etikprövningsmyndigheten, Sweden.

| Anatomic space definition
Space is defined along the three anatomical dimensions, the sagittal axis perpendicular to the facial plane, the frontal axis perpendicular to the sagittal plane, and the longitudinal axis perpendicular to the transverse plane (Figure 2).

| OCT procedure
Before OCT capture, the eye was dilated with tropicamide 5 mg × mL −1 .All OCT volumes were captured with OCT-Triton (Topcon, Japan) in a low-level light environment.The volume captured is, according to the manufacturer, 2.6 mm resolved in 992 px along the sagittal axis, 6 mm resolved in 512 px along the frontal axis, and 6 mm resolved in 256 px along the longitudinal axes.This provides a resolution of 3 μm along the sagittal axis, 12 μm along the frontal axis, and 23 μm along the longitudinal axis.

| Image analysis
Two separate deep learning algorithms were used for automatically annotating OPCL and delineation of the inner limit of the retina, respectively.Both algorithms were fully convolutional networks (fCNN) utilizing 2D U-Nets (Ronneberger et al., 2015).Images used for teaching the algorithm were not used in data analysis.The details of the location of OPCLs and the inner limit of the retina are provided elsewhere (Brusini et al., 2021;Carrizo et al., 2020).In short, the backscattered intensities within the OCT representation of the ONH in Cartesian coordinates were resampled to polar coordinates, with the origin (first reference origin) in the centre of the ONH at an arbitral height on the sagittal axis.Then, the OCT volume of the ONH was represented on 500 equidistant radii, approximately perpendicular to the frontal plane.The location of the inner limit of the retina was defined in space with an AI model.The location of the OPCL (Figure 1) in each radius was automatically detected with a second AI model.Finally, the location of an IRCP corresponding to each OPCL was automatically detected.
Following the automatic annotation, OPCL locations are least-square fitted to a plane in 3D space (Figure 3).
Subsequently, the best-fitting plane of the OPCL coordinates together with all OPCL and IRCP locations are rotated in space with an unchanged spatial relationship between them (Figure 3a, thick lines).Thus, the normal vector to the best-fitting plane of the OPCL coordinates is parallel with the sagittal axis.The new translocated OPCL and IRCP coordinates are named the OPCL′ and IRCP′ coordinates.
Then, a circle with a plane parallel to the best-fitting plane is fitted to the OPCL′ coordinates to establish the origin of the circle, secondary reference origin (Figure 3b).
In the next step, all OPCL′ and IRCP′ coordinates are translated so that the new OPCL′ coordinates (OPCL″ coordinates) lie in the best-fitting plane with an unchanged spatial relationship between each OPCL″ coordinate and the corresponding new IRCP′ coordinate (IRCP″ coordinate).
For each OPCL″ coordinate, a radial OPCL plane parallel to the sagittal axis is created (Figure 4a).
Thereafter, the IRCP″ locations around the radial plane corresponding to a specific OPCL″, OPCL″ location of interest, are located.Finally, the cross-over between the radial plane and the straight line running through adjacent IRCP″-locations is established mathematically (Figure 4b).The cross-over point is considered the IRCP location of interest.The distance between IRCP″ and OPCL″ location of interest and the corresponding IRCP location of interest is defined as PIMDv2 (Figure 4b).

| Resolution
When plotting the transformed OPCL″ and IRCP″ locations, it was apparent that their locations in 3D space sometimes are more discontinuous than expected from the structure of the ONH.The discontinuities are caused by errors in the spatial localization of each OPCL and IRCP due to limited resolution in the OCT capture.The discontinuities are visualized when plotting the distance in space between the secondary reference origin and the OPCL″ location as a function of angle (Figure 5a, red).
The errors in spatial localization of OPCL and IRCP were filtered out by applying a low-pass filter.The lowpass filter consisted of averaging pixels in each of the three dimensions in space.The filtering window size was set to 11 pixels in all three dimensions.
A Fourier transform of the length of the OPCL″ vector revealed that most of the low-frequency content describing the structural variation of the length of the OPCL″ vector is below 40 cycles/2π (Figure 5b).
This corresponds to a period of 1/40 radians.For this reason, the angular resolution is limited to approximately slightly more than a quarter of a clock hour.
When low-pass filtering the three-dimensional OPCL″ vector, most of the high-frequency noise disappeared (Figure 5a, black).

| Waist area
We suggest that the cross-sectional area of the waist of the nerve fibre layer in the ONH is defined as Pigment epithelium Inner Limit of the retina Minimal Area (PIMA).
The surface element of PIMA at a specific angle in the frontal plane, PIMAα (Figure 4b), is defined as the surface enclosed by two adjacent OPCL″s and the corresponding adjacent IRCP locations of interest.The angle in the frontal plane (α) is defined with 0 radians at the nasal direction incrementing towards the superior direction.Then, the surface area of the PIMA element is a trapezoid estimated as the sum of two triangles (Figure 4b).
Adding all surface elements within an angular segment incrementing from α 1 to α 2 provides the total surface for the angular segment, PIMA-[α 1 ; α 2 ], and the total sum over 2π radians, PIMA-2π, provides an estimate of the total minimal cross-sectional area of the waist of the nerve fibre layer.
Variation in OPCL and IRCP locations introduces the Coastline paradox.If the radius of a boundary opposite to the origin of the radius varies as a function of angle, the boundary length varies proportionally.When strictly increasing radial change per angle change, the boundary length increases.Both the OPCL locations and the IRCP locations are estimated with a random error.Consequently, PIMA-2π will strictly increase with increasing resolution of measurement error.Therefore, it is necessary to reject nonstructural information in the boundary by low pass filtering of the three-dimensional locations of OPCL and IRCP.
To verify the consequences of filter window size, PIMA integrated over the circumference of the ONH, PIMA-2π, as a function of filtering window size, was examined (Figure 6).
Below a filtering window of 11 pixels, PIMA-2π drops quickly, probably due to filtering of random noise caused by limitations in the 3D resolution of the localization of OPCL and IRCP.However, above a filtering window of 10, PIMA-2π decreases strictly at a low rate, probably due to increasing loss of relevant structural information.

| Experimental design
One eye from each of the 16 subjects, 8 males and 8 females, was analysed for PIMD-angle and PIMA-angle.

| Statistical parameters
The confidence coefficient and the significance limit were set to 0.95 and 0.05, respectively.
When resolving PIMA into angular segments, the typical superior and inferior peak thickness of the nerve fibre layer is reflected (Figure 7).
When resolving PIMDv2 into angular segments, the typical superior and inferior peak thickness of the nerve fibre layer is reflected (Figure 8).

| DI SC US SION
The current paper introduces a strategy to estimate the cross-sectional area of the waist of the nerve fibre layer at the ONH in a 3D grid representation of the ONH, captured with OCT.Further, a computational method for estimating the cross-sectional thickness is presented.
The current paper was designed to develop an algorithm that allows the estimation of minimal crosssectional surface area of the nerve fibre layer.The developed algorithm was applied on all available data from an ongoing study for the determination of the age dependence of PIMD-2π.Differences in regard to ONH disc size remain to be established.To verify clinical validity, the presented algorithm should be applied on captures of glaucomatous eyes.The sample size is limited but sufficient to provide approximate estimations of levels of quantities measured.
The OCT device was selected to allow analysis in OCT volumes generated by a commercially available standard device in ophthalmological clinical practice.The fact that different axes are collected with different spatial resolutions makes extraction of distance complex and poses a risk for image artefacts.In the present system, the angular resolution appears to be limited to approximately 10/500 × 2π radians, or slightly less than a quarter of a clock hour.
Originally, the inner limit of the tissue surrounding the nerve fibres in the ONH was defined as BMO (Reis et al., 2012).The thickness of Bruch's membrane is below the resolution limit while the pigment epithelium is clearly resolved in current commercial OCT devices (Sandberg Melin et al., 2019a).Therefore, Sandberg et al. proposed to use the central limit of the pigment epithelium as the delimiter to the tissue surrounding the nerve fibres in the ONH.On the condition that pigment epithelial cells always are associated with a basement membrane, point estimations of BMO and OPCL are anatomically identical.
In our original publications (Sandberg Melin et al., 2019a, 2019b), PIMD was estimated as the spatial distance between OPCL and IRCP.Due to the limitation in the spatial resolution, often, the same IRCP corresponded to consecutive OPCL coordinates.This caused a slightly slanted PIMD in relation to a best-fit plane of the OPCL coordinates.Therefore, the spatial distance was overestimated.In the current estimation, the best-fit plane of the OPCL coordinates was used as a reference plane (Figure 3a).The IRCP corresponding to a specific OPCL coordinate is located between adjacent IRCPs perpendicular to the best-fit plane (Figure 4b) by assuming a straight line in space between adjacent IRCPs (Figure 4b).Thus, PIMDv2 is always perpendicular to the best-fit plane of the OPCLs and therefore expected to be closer to the truth and, on average, shorter than PIMDv1.The currently estimated confidence interval for the mean difference between PIMDv1-2π and PIMDv2-2π did not indicate any difference.The lack of difference is probably due to the angle between PIMDv1 and the sagittal axis, on average, being very small.
Assume that nerve fibres have a constant diameter in the proximity of the ONH.Then, a widening of the central limit of the pigment epithelium over time in a patient without loss of nerve fibres is expected to decrease the thickness of the waist of the nerve fibres.This is due to the distribution of the nerve fibres over a larger surface.In contrast, the total cross-sectional area is expected to be constant.The diameter of the optic nerve head is known to be very variable among individuals.Thus, a potential Gaussian distribution of estimates of the waist of the nerve fibre layer in the ONH in individuals without pathology, normalized to the mean, is expected to be wider for thickness estimates than for estimates of the cross-sectional area.Therefore, if the frequency distribution for non-pathological individuals is used to set the threshold for pathological loss, it is expected that estimates of the waist of the nerve fibre layer in the ONH expressed in thickness are less efficient for the detection of pathological loss.Alternatively, the threshold for pathological loss can be set based on the frequency distribution for that individual.Then, for estimates of the waist of the nerve fibre layer expressed in thickness, a change of the diameter in that individual over time is expected to cause a widening of the frequency distribution over time.If estimates of the waist of the nerve fibre layer are instead expressed as surface area, the frequency distribution is expected to remain the same.On the contrary, it is known that the fibre diameter of retinal ganglion cell axons is variable.This may contribute to increased variability in area measurements (Sanchez et al., 1986).
Already Považay et al. presented both average thickness and average cross-sectional area of the waist of the nerve fibre layer (Považay et al., 2007).Several other authors published both quantities.The sensitivity and specificity of parameters measuring waist thickness and waist area when comparing non-glaucomatous and glaucomatous individuals have been evaluated in two studies with contradicting results (Enders et al., 2019;Li et al., 2021).For a specificity of 90%, Enders found that the sensitivity of BMO-MRA was slightly higher than for BMO-MRW (72.7% and 71.9%, respectively) when adjusting for different ONH sizes between the glaucomatous and the control group (Enders et al., 2019) while Li concluded the opposite when adjusting for ONH size (75.6% and 82.9% respectively) (Li et al., 2021).For both studies, Area Under the Curve (AUC) differences between BMO-MRA and BMO-MRW were not significant.When comparing the glaucomatous and control group without adjusting for different ONH sizes, Enders found a larger difference in sensitivity at 90% specificity between the parameters BMO-MRA and BMO-MRW (70.1% and 66.0%, respectively).This supports the idea that cross-sectional waist area might be a more efficient parameter.Recently, it was found that there is no difference in efficiency between PIMD-2π and PIMA-2π at 0.95 specificity considering averages among subjects (Söderberg et al., 2022).The lack of difference is probably due to substantial variability among individuals (Sandberg Melin et al., 2019b).Within-subject comparison may result in different efficiencies.
In the Spectralis SD-OCT system (Heidelberg Engineering), the ONH is represented by radial OCT scans.In each radial scan, BMO (corresponding to OPCL) and the inner limit of the retina are delineated.Then, the BMO locations are best fitted to a BMO plane.BMO-MRW is then defined as the shortest distance perpendicularly to the BMO-plane between the BMO plane and the inner limit of the retina in the same radial scan.BMO-MRA is defined as the trapezoidal surface enclosed by consecutive BMO-MRWs (Enders et al., 2016;Gardiner et al., 2014).
The currently presented algorithm, PIMA-2π estimated from a 3D-grid representation of the ONH, corresponds anatomically to BMO-MRA.The estimated PIMA-2π was of the same magnitude but slightly higher than similar estimates of the cross-sectional area of the waist of the nerve fibre layer as BMO-MRA (An et al., 2021;Enders et al., 2016;Gardiner et al., 2014;Li et al., 2021;Považay et al., 2007).BMO-MRA is estimated as the sum of trapezoidal surfaces with the base centred on BMO and potentially underestimates the real minimum waist cross-sectional area (Enders et al., 2016).Our voxel representation of the 3D structure implies that the cross-sectional area of the waist undulates around the centre of rotation.PIMA-2π is the sum of surface elements along the undulating cross-section and, therefore, is expected to be slightly higher.The frequency analysis in Figure 6 indicates that sampling bias was limited, while the integrity of real undulation was preserved.It cannot be excluded that some of the difference between PIMA-2π and BMO-MRA measurements is due to differing calibration of the OCT devices used.To prove the clinical significance of recording undulations, comparative measurements between the present algorithm and the radial scan algorithm (Gardiner et al., 2014) are required.
The capacity of the deep learning model to correctly annotate ONH's with abnormal anatomy (e.g.peripapillary atrophy, macro-and microdiscs, coloboma) can be achieved by training the algorithm with such samples.If pathological samples are included for training the current segmenting deep learning algorithm, the reliability for abnormal anatomy in the ONH needs to be proven for each condition.
Today, retinal Nerve Fibre Layer Thickness (RNFLT) is clinically used to monitor glaucoma.A recent comparison between BMO-MRW and RNFLT showed a slightly lower sensitivity for BMO-MRA (Li et al., 2021).A reason for this may be an underestimation (Enders et al., 2016).
The angular variation of PIMA (Figure 7) and PIMD (Figure 8) along the frontal plane around the ONH agreed with previous findings when estimating the average nerve fibre layer thickness close to the ONH; by circular scanning around the ONH (Jaffe & Caprioli, 2004), when estimating the average thickness of the nerve fibre layer close to the optic nerve head with polarized light (Cense et al., 2004), and when estimating the waist of the nerve fibre layer by AIsupported segmentation of OCT volumes of the ONH (Carrizo et al., 2020).
The confidence interval for the difference between PIMDv1 and PIMDv2 estimated in the current study did not elucidate a potentially higher PIMDv1 than PIMDv2.Considering that PIMDv1 corresponding to a defined OPCL often is slanted, decreased thickness of the nerve fibre layer will result in an overestimated PIMDv1.Therefore, PIMDv2 may be more appropriate for detecting local deviations of nerve fibre layer thickness in early glaucoma.In the current sample of healthy young eyes, the nerve fibre layer was probably thick enough to abolish the effect of the slant angle.The difference between PIMDv1 and PIMDv2 in thin nerve fibre layer sectors in glaucoma patients remains to be settled.

F
I G U R E 2 Schematic image with anatomical space definitions.F I G U R E 3 Transformation of ONH pigment epithelium central limit (OPCL) and inner limit of the retina closest point (IRCP) coordinates.(a) Red line: OPCL original coordinates before (upper thin dots) and after translation (lower thick dots).Red surfaces are best-fit planes to the corresponding OPCL coordinates.Blue line: IRCP original coordinates before (upper thin dots) and after translation (lower, thick dots).(b) After fitting translated OPCL coordinates (red dots) to a circle (black interrupted line).

F
I G U R E 4 (a) Visual 3D representation of the cross-section of the waist of the nerve fibre layer in the optic nerve head.The green surface is a radial plane perpendicular to the best-fitting plane of the ONH pigment epithelium central limit (OPCL) locations in space (red) transecting an OPCL location of interest.(b) Schematic of one trapezoidal surface element of PIMA constituted by two triangles (striped areas).The green surface perpendicular to the PIMA surface element corresponds to the green surface in (a).a) Red line: Length of ONH pigment epithelium central limit (OPCL)″ vector originating from (0.0) as a function of angle.Black line: OPCL″ vector after application of a low-pass filter excluding frequencies above 40 cycles/2π.(b) Fourier transform of the length of OPCL″ vector as a function of number of frequencies.Red: Frequency content of OPCL″ vector as a function of angle before low-pass filtering.Black: Frequency content after low pass filter excluding frequencies above 40 cycles/2π.

F
I G U R E 6 PIMA-2π as a function of filter size.Sample from one subject.F I G U R E 7 PIMA-2π resolved into angular segments.Each bar represents the total surface area of 10 consecutive angles (0.13 radians) in one subject.F I G U R E 8 PIMDv2-2π resolved into angular segments.