1. Introduction

The development of non-invasive methods to estimate ocular rigidity (OR) will have profound implications for research into ocular disease. Importantly, glaucoma remains a major cause of blindness due to formidable challenges in both its diagnosis and treatment, and its pathogenesis is poorly understood. Reducing intraocular pressure (IOP) is the most widely used clinical method for halting the progression of open angle glaucoma (OAG). However, the link between IOP and development of OAG is not straightforward [1–7]. Considerable recent evidence from experimental studies in primates and from mathematical modeling suggests that ocular biomechanics may play a major role in glaucoma pathogenesis [8–17]. According to finite element modeling, major determinants of optic nerve head stress and strain leading to glaucoma damage include IOP, but also scleral elasticity as well as other biomechanical factors. In fact, scleral elasticity is considered to be the most important determinant of optic nerve head stress and strain, more important than IOP [18] and it is clear that additional factors, such as ocular biomechanics, must play an important role.

Additionally, several investigations into age-related macular degeneration (AMD) have led to both mechanical and ischemic theories of pathophysiology related to OR, particularly in neovascular AMD [20, 21], but it remains unknown as to whether changing rigidity plays a role in the pathophysiology of this disease. Reduced scleral rigidity is also an important feature of pathological myopia [19].

The rigidity of the eye can be derived from the Friedenwald’s empirical function that estimates the change in IOP produced by a modification of the ocular volume V, according to:

The ocular volume fluctuates due to the pulsatile vascular filling. Since 90% of blood flow into the eye is though the choroid [25], we propose to model the fluctuations of ocular volume by estimating the total choroidal volume change over time.

Although to date investigations into the elastic properties of the eye have produced values for this rigidity constant, they have either been on post-mortem testing [4, 12, 16], they make use of invasive cannulation in order to control the volume change [19], estimate indirectly the ocular volume change [26, 27] or use measurements that intrinsically depend on ocular rigidity to estimate OR [27, 28]. For clinical applications, ocular volume change (and thus OR) needs to be measured non-invasively, but no technology is available to measure them in a clinical environment.

Recent advances in optical coherence tomography (OCT), specifically with Enhanced Depth Imaging (EDI), have improved the signal to noise ratio in deeper tissues to the point that the choroidal-scleral interface (CSI) can now be distinguished. Segmentation methods exist to delineate this interface from high quality images [29–32], but the fast acquisition required to assess the changes of choroidal volume with pulsatile blood flow limits the imaging signal to noise ratio (SNR), thereby complicating the segmentation.

In this paper we propose a novel method for automated choroid segmentation in sequential FD-OCT images that is relatively robust to noise and low image quality, and which allows us to estimate the volumetric changes of the eye due to choroidal pulsations. These measurements, in combination with intraocular pressure measurements and biometry, allow the first non-invasive, direct estimation of OR.

2. Method

In order to track the pulsatile volume changes due to choroidal filling, B-scans have to be acquired faster than the heart rate, which renders the CSI nearly indistinguishable from noise. Our method computes the area between the posterior part of the RPE layer and the CSI to extrapolate choroidal volume changes based on a simple model detailed below. We have found that previous choroidal segmentation methods are not robust enough for our application; hence a new approach is required.

We combined a robust contour-detection method with a graph search based on a novel weighting scheme to develop a segmentation algorithm that boosts the reliability of CSI delineation, as described below.

2.1 Data collection

Images of the choroid were acquired using a FD-OCT (Spectralis OCT Plus, Heidelberg Engineering, Germany) system whose software was modified to provide time series where each frame results from the average of an adjustable number of B-scans. The acquisition was set to high-speed mode (496 pixels per A-scan x 768 A-scans images), enhanced depth imaging (EDI), using a class 1 laser at 870nm, and 30° wide (8-9 mm on the retina, optimized for each subject). With these settings, B-scans can be acquired at a maximum of 40Hz, at 3.9 μm axial and 11 μm of lateral resolution, and at 400 images the memory buffer is filled. The number of B-scans per frame is determined by the required time resolution of the series but it also impacts the quality of individual frames. Averaging 5 scans per image (8 Hz sampling) was an acceptable trade-off. Measurements were centred on the fovea and the azimuthal angle was chosen to maximize the visibility and continuity of the CSI for each subject. While a full 400-frame movie was acquired the subject’s heart rate was measured with a finger oximeter.

The system is equipped with an eye-tracker to keep the scanning beam in place, but this feature also introduces pauses into the acquisition, producing fluctuations both in the number of averaged B-scans per frame and in the acquisition rate. Since the resulting time series are unequally spaced, the image’s time-stamp was used when computing the frequency spectrum. Only frames with Spectralis’ quality parameter above 20 were kept in the time series.

Immediately after imaging, the intraocular pressure is measured with a Pascal Dynamic Contour Tonometer, which is not dependant on corneal rigidity [33]. The average of three measurements having a Quality Index (Ziemer proprietary algorithm) not below 3 is computed. This provides two values, the intraocular pressure (IOP) at diastole and the ocular pulse amplitude (OPA), which represents the difference between diastolic and systolic pressure.

2.2 Preprocessing

Depending on eye size, a variable portion of the retina may be visible in each movie. Since we are only interested in estimating the volumetric changes of the eye due to choroidal filling (extrapolated from 2D images) we discard A-scans near the optic disk, where the choroid is absent (Fig. 1(A)). Every image in the time-series is then aligned to the first one, using Matlab (The MathWorks, Natick, MA) imregister function with mean squared error metric and one-plus-one evolutionary optimizer. The registration is limited to rigid transformations (i.e. no shear or dilation) to prevent biasing the measurements with artificial image distortions. Each frame is then analyzed independently, provided the Spectralis’ metric of Quality is not lower than 20, and at least two raw scans are averaged to create the frame used for further analysis.

 

Fig. 1 Automated segmentation of retinal layers of interest from OCT images. The image in A is a typical frame from the video series. A) A-scans where the choroid is absent (highlighted in green) are discarded from all frames. B) Segmentation of the outmost layers of the retina: RVI is indicated in green, anterior RPE in red, and posterior RPE in blue. The A-scans are shifted so that the blue layer appears flattened (C).

Download Full Size | PDF

To enhance the likelihood of correctly delineating layers of interest, despite variations in individual retinal shapes, the algorithm identifies several retina layers sequentially, in an anterior to posterior direction. We first normalize each A-scan to its maximum intensity, remove noise with a 5x5 pixels Wiener filter, and apply a Canny edge detector with thresholds of 0.01 and 0.3 The Gaussian filter size of the edge detector (σ = 4) is chosen large enough to ensure the top layer, the retina-vitreous interface (RVI), is continuous. This layer is segmented by joining the topmost 8-connected edge segments that are wider than 50 A-scans with cubic splines (Fig. 1(B), green line).

The next two segmented layers, the anterior and posterior interfaces of the retinal pigment epithelium (RPE), have been segmented using a previously published strategy [31], where a few modifications have been added to improve robustness. We profit from the positive intensity gradient that separates the neuroretina and the RPE, and the negative gradient between RPE/Bruchs membrane (BM) and the choroid to delimit these two interfaces more robustly by finding them together. After computing the gradient of the raw image and smoothing with a Gaussian kernel (σ = 3pix horizontally and σ = 0.5pix vertically), we search for its positive maximum in each A-scan between 39µm and 780µm below the RVI. The resulting points are connected with a local 2nd degree polynomial least squares weighted fit, to render the RPE (Fig. 1(B), red line). Analogously, the negative gradient minima found between 39µm and 117µm below the anterior RPE, are connected to give the posterior RPE (Fig. 1(B), blue line)

Since we use a graph search method to segment the CSI, each image is 'flattened' with respect to the RPE in order to eliminate erroneous paths introduced by the curvature or tilt of the image. This is done by shifting and zero-padding each A-scan until the pixels that describe the posterior RPE are vertically aligned [31] (Fig. 1(C)).

2.3 CSI segmentation

We developed a segmentation method for the CSI to match the particular requirements of our application. In EDI-OCT scans, the CSI is a remarkably heterogeneous boundary consisting of fragments of blood vessel cross-sections, which cannot be segmented with usual edge detection approaches. Graph search edge detection is especially well suited for this problem as has already been shown [31]. Briefly, this approach associates pixels that loosely describe the target interface to nodes in a graph, and minimizes the path across the nodes based on weights assigned to each connection. The reliability of the method depends strongly on the choice of nodes and the weights. Previous implementations of graph search to segment the CSI do not suit the present application where subtle changes are essential for accurate measurements of ocular rigidity.

We propose a novel approach that combines graph search with a robust contour detection method, which additionally profits from information gathered in time to boost the reliability of the segmentation.

Node locations

In earlier work, graph nodes were found using variations of image intensity along each A-scan [31]. Tian et al. looked for local maxima and minima of intensity, and used their valley pixels (the local minima) as nodes. Our implementation of this approach is unsuitable for this work for two reasons. Firstly, with this approach local minima are placed in the center of visible blood vessels rather than their bottom border, where the real CSI lies. Secondly, due to high noise, node detection is unreliable. Our approach provides a much improved weight function that favours nodes located at regions of higher local contrast and also profits from time-series information.

We find nodes using the smoothed first and second gradient of image intensity along each A-scan. The input in both cases is the A-scans of a preprocessed image, which has been smoothed with a span of 10 pixels. Pixels in which the first derivative exceeds a positive threshold of 0.7 identify the dark-to-light transition characteristic of the deepest interface of the choroid blood vessels. Additionally, those pixels whose second derivative absolute values are smaller than a near-zero threshold (10⁻¹⁶) mark the inflection point of intensity on the lower extremity of the transition. A binary image meeting both conditions undergoes a sequence of morphological operations that reduce the detected regions to isolated pixels. First, regions are cleaned to eliminate single-pixel regions. Then, regions are filled to eliminate holes, prior to being skeletonized in order to retain only the central pixels of each connected region. Next, all pixels in every other column are eliminated to reduce the number of nodes the graph search must include. Images are shrunk to ensure any remaining regions are single pixels. Then extended-minima transformation is computed for the original preprocessed intensity image with a threshold of 10 pixels, to find the rough central shadow of large blood vessels, and any nodes that intersect these shadows are eliminated. Finally, those pixels at a depth greater than 150 pixels (585 µm) from the BM are also discarded (green dashed line in Fig. 2 E), since the CSI is unlikely to go this deep. The remaining pixels pinpoint the nodes fed to the graph search (yellow pixels in Fig. 2 E).

 

Fig. 2 A) The uncompensated sub-RPE region. B) The same region, compensated and contrast enhanced, to which the oriented gradient algorithm will be applied. Overlaid is an example disk of radius 30 pixels. The red and green regions correspond to their respectively coloured histograms (D). C) The oriented gradient image, composed of the combination of the X² images of different scales and orientations, as described. The heat map shows pixels which are very likely to lie on a boundary. Even weak boundaries can be detected while excluding noisy regions with this method. E) Overlay of the oriented gradient image (heatmap), node locations (yellow x’s), and the CSI found using these two inputs to the graph search (redline), onto the flattened B-scan. The green dashed line shows the limit of 585 µm below the Bruchs beyond which nodes are discarded. F) The original B-scan overlaid with the RPE (blue), CSI (yellow) and the mean RPE-CSI distance or CT (red dotted line). This distance CT is what is tracked from frame to frame.

Download Full Size | PDF

Graph search

The graph is constructed by connecting each node to all other nodes in the neighbourhood delimited by C_max columns to the right and R_max rows above and below it. C_max must be sufficiently large to ensure the resulting graph is connected, even across dim regions with sparse nodes [31]. Connections between each pair of nodes a and b are assigned weights according to:

Indeed, the CSI is most likely smooth, but the minimal distance condition does not guarantee a smooth interface. In fact, the shortest path between any two nodes is a straight line that avoids intermediate nodes, even if they are close together. Actually, the hard thresholds C_max and R_max prevent the graph search from just tracing a straight line between start and end nodes, because they limit the maximum length of edge segments. Unfortunately, these parameters cannot be reduced arbitrarily since they have to be long enough to overcome gaps produced by missing nodes. Therefore, in order to improve the edge smoothness while preventing gaps, we added a horizontal weight term

Due to the inhomogeneous nature of the CSI, combined with the low SNR resulting from high-speed acquisition, the node locations retrieved with the method above are not reliable enough for the graph search to delineate the right interface. To reinforce the reliability of the segmentation we compute a boundary probability which we use to compute a connection weight w_Affin that favors paths through the most likely boundary.

An excellent choice for this is a multi-scale, multi-orientation approach such as the contour-detection algorithm by Arbelaez et al. [34]. This method computes the posterior likelihood X² that each pixel (i,j) in the image belongs to a boundary of scale σ at orientation θ. The computation uses the histograms g(I) and h(I) of pixels intensities in the two halves of a disk of radius σ, centered on (i,j), and divided along its diameter at an angle θ (Fig. 2(B) and 2(D)), according to the formula:

Even with EDI, the OCT signal from deep sections of the eye is often weak and inhomogeneous due to the presence of attenuating structures above them, including blood vessels. In order to improve the reliability of P_b, we used the adaptive image enhancement proposed by Mari et al, to increase the contrast in the region extending from 5 pixels below the RPE to the bottom (Fig. 2(B)) [29, 35]. The weighting of the graph search, node locations and connectivity matrix are fundamentally different from previous approaches that used such compensation [29].

Since the computation of P_b is time intensive, we restricted it to a region no deeper than 150 pixels below the RPE, as depicted by the green dashed line in (Fig. 2(E)), where the CSI is most likely found.

Finally, the term w_Affin is computed as

For the start and end nodes, two 'virtual' nodes are added before the first and after the last columns, and are connected to the nodes inside the image as per the restrictions on Cmax. The graph search then uses Dijkstra's algorithm [36] to find the minimum-weight path between these virtual nodes. The resulting path is finally interpolated and smoothed to render the CSI boundary.

Due to high noise, delineation errors may occur in some frames, yielding unrealistic CSI profiles. Assuming that the CSI should not undergo significant change in shape during the cardiac cycle, we use the contours computed on all frames to correct for these outliers. First, we compute the mean CSI curve over all frames, and we measure the correlation of each individual frame CSI to the mean curve. We recomputed the frames whose CSI correlation fall below the mean correlation value, those in which the total area (in pixels) enclosed between the posterior RPE and CSI differs more than 3000 from the median, and those whose depth of the first or last pixel of the CSI deviate by more than 15 pixels from their respective means. For this second computation an additional weight term is included, as:

2.4 Computation of ocular rigidity

Once the RPE and CSI delineations are determined, the mean RPE-CSI distance across A-scans is computed in each frame, giving a time series of sub-macular choroidal thickness (CT), as shown in Fig. 3(A) (bottom) as a black line. As expected, the frequency spectrum analysis on most time series also revealed high frequency components coincident with the first and second harmonics of the heart rate frequency F_H (see Fig. 3(B), black line), which we measured independently from the oximeter signal (Fig. A and B, top). This correlation proves that CT fluctuations in time are, at least in part, due to the pulsatile blood flow. For spectral analysis we used the Lomb-Scargle periodogram [37] instead of the popular Fourier Transform because CT time series are unequally sampled. This is due to images of Spectralis quality metric below 20 being omitted from the time series or to pauses in the acquisition due to the eye tracker.

 

Fig. 3 Spectrum analysis of CT fluctuations in time. A) Top: Oximeter signal. Bottom: Raw fluctuations of CT versus time (black). Overlaid in red is the band-pass filtered CT signal (red). B) Frequency spectrum of the oximeter signal (top), and CT signal (bottom), where the offset component has been omitted. The filtered frequency band for the CT spectrum is shown in red..The dashed blue lines indicate the two first harmonics of the measured heart rate which are observed in both signals.

Download Full Size | PDF

In order to extract the CT fluctuations associated mainly to the heart rate and discard respiration, head movement, saccades and segmentation noise, we filter out frequency components below ½ F_H, above 3F_H, and those with values below 10% of the maximum peak within this range (Fig. 3(B), red line). The inverse Fourier Transform was used to retrieve the filtered signal (Fig. 3(A), red line).

The pulsatile fluctuation of CT is obtained by using a windowed peak-to-valley algorithm, which ignores peaks and valleys that are spaced in time less than 1/6 of the heart period (T_H = 1/F_H), or that are smaller than 10% of the maximum peak. All sequential peak-to-valley gaps, which are greater than the vertical resolution of the Spectralis OCT (4 µm for the used settings), are averaged to render the final ΔCT.

The ocular volume change ΔV = V-Vo is derived from ΔCT using a first order approximation of a spherical eye model. In it, the choroid is modeled as the volume between two spheres shifted by ΔCT, as:

2.5 Subjects

In order to test our method, we enrolled 45 subjects from the Ophthalmic Clinic. Since the aim of the study was to test the ability of our approach to measure OR, the sole eligibility criteria imposed on the subjects was that eyes had to be clear enough media to allow distinguishing the choroid-sclera interface from the OCT images. The study protocol adhered to standards outlined in the Declaration of Helsinki. All participants were informed of the nature and objective of the study, the procedures that would be performed and the risks, and gave their informed consent. The local research ethics committee approved the protocol.

3. Results

This work introduces a novel method for assessing the Ocular Rigidity non-invasively from OCT time series and usual biometric measurements. One important component of the method is the algorithm capable of segmenting the CSI from low SNR images. In the next subsections we describe experiments that demonstrate the improvements provided by this novel CSI segmentation, as well as OR measurements on our cohort of subjects that validate the methodology.

3.1 Evaluation of CSI Segmentation

In contrast to some segmentation problems, there is no gold standard for the choroid contour to compare with; therefore we need to use an alternative approach to validate the new algorithm. Ophthalmologists and other eye specialists typically assess OCT images depicting the choroid qualitatively; hence they are the best-suited people to evaluate the CSI segmentation efficiency. We presented a set of 25 OCT images to 5 independent specialists using a tablet equipped with a stylus to delineate the choroidal-scleral interface on top of the OCT images. We retrieved the manual contours and compared them with the automated traces yielded by our method and Tian’s [31] (for example, see Fig. 4(A) and 4(B)). For each image, we calculated the average manually segmented CSI based on the 5 independent traces and compared the performance of each specialist and both automated methods. Histograms of deviation from the mean trace for every A-scan of all images are shown as violin plots in Fig. 4(C). Such histograms illustrate the improved capacity of our method, which is virtually indistinguishable from the specialist’s.

 

Fig. 4 Segmentation results for the proposed method (blue), Tian’s method (red) and manual tracings (green). A set of 25 OCT images was manually segmented by 5 independent specialists and used for benchmarking our method. A, B) Examples illustrating the performance. C) We calculated the average manually segmented CSI and compared the performance of each specialist and both automated methods. Histograms of deviation from the mean trace for every A-scan of all images are shown as violin plots

Download Full Size | PDF

The image shown in Fig. 4(A) illustrates one of the main problems of Tian’s method mentioned above. The nodes are mostly located in the middle of the vessels, rather than on their bottom edge, and this causes the overall offset that can be observed in Fig. 4(C). Furthermore, noise can cause the image to be crowded with intensity minima, which renders a purely Euclidean shortest path to be located blow the CSI (see Fig. 4(B)).

3.2 Measurements of ocular rigidity

In order to test our approach, we measured the ocular rigidity on a group of subjects, as described above. The mean CT was 263.8 (SD = 78.4) µm. The mean magnitude of thickness change at the macula ΔCT was 16.7 (SD = 10.9) µm. The pulse volume ΔV was 7.8 (SD = 4.9) µL, and this was used to estimate the pulsatile ocular blood flow (POBF) 595.6 (495.1) µL/min, by multiplying with heart rate (HR) measured while OCT images were acquired. Finally, the mean OR constant in the tested set was 0.028 (SD = 0.022) 1/µL.

A positive correlation was observed between OR and OPA (Fig. 5(A)). Although small, this correlation is consistent with the reasoning that for a more rigid scleral shell, larger pressure pulses are expected. More importantly, we observed a statistically significant negative correlation between OR and the axial length (AL) (see Fig. 5(A)) which agrees with a recent study that uses cannulation to artificially modify the eye volume [19].

 

Fig. 5 Ocular rigidity measurements. A) Correlations of OR with OPA and AL. B) Reproducibility results (4 subjects). Along with ICC value, its lower and upper bonds for alpha = 0.05, are indicated.

Download Full Size | PDF

Finally, four subjects were measured repeatedly to assess the reproducibility of the method. The results plotted in Fig. 5(B) show that the OR assessment is reproducible across the whole range of rigidity values observed in this study. In fact, these measurements render an intraclass correlation coefficient of 0.96, with a lower confidence band (at α = 0.05) of 0.63. Altogether, these results validate our method to measure the ocular rigidity.

4. Discussion

The Friedenwald equation is the currently accepted conceptual framework for investigating the pressure-volume relationship in the eye. It synthesizes the pressure-volume curve in a single value, the ocular rigidity constant, which simplifies the study of relationships to other variables [19, 26, 27, 38]. It is important, however, to ensure that the quantities involved in this formula represent the real physical changes occurring in the eye. Given that the primary physiological cause of pressure fluctuations in the eye is the pulsatile choroidal volume change, the method we present offers the most representative determination of ocular volume change through non-invasive direct quantification, rather than through indirect variables such as fundus pulse amplitude (FPA) [27] or laser Doppler flowmetry of pulsatile choroidal blood flow [11]. As opposed to FPA, our method directly measures the expansion of the choroid produced by blood inflow.

The method used to calculate ocular volume change from change of choroidal thickness at the macula is a first order approximation. The fluctuation ΔCT is found in a relatively small (~9mm) section at the fundus, which amounts to between 15 and 20% of the average circumferential length of the choroid of the eyes included in this study, but accounts for a greater percentage in terms of total choroidal blood flow since the macular region has the greatest perfusion. Choroidal thickness has been shown to vary not only in the temporal-nasal directions [39], but more generally in all directions [40]. However, mean CT found using our automated segmentation (263.8 µm) agrees closely with other studies centered on the fovea making use of manual CSI segmentation [39, 40]. Additionally, the pulse volume ΔV and POBF agree with estimated values obtained both using commercial devices that assume a given ocular rigidity [41], and using the slope of pressure-volume curves and OPA [19]. Together these results provide good evidence for not only the proper segmentation of the choroid, but also that the dynamic change in volume is being well computed. Importantly, as opposed to some other approaches to measure choroidal blood flow, the results obtained with our method based on OCT segmentation is neither sensitive nor biased by axial head movement.

Alternative equations to refine the volume estimation as a first order approximation in ΔCT can be used, but the final results for OR will only differ by some common factor, and comparisons among populations or the study of trends would remain legitimate. What is more noteworthy is that this first-order approximation leads not only to OR values that are the right order of magnitude, but also in the same range of values reported in earlier in vivo studies [19, 21, 24, 27]. A future refinement of the choroid volume estimation may be to use a wide angle OCT for collecting data from a broader area of the fundus, or even volume reconstructions. Nevertheless, such volumes must be obtained significantly faster than the heart rate, which is not possible using technology that is currently in widespread clinical use . Of note, the amount of time and processing required to analyze time series of volume reconstructions would be increased by approximately two orders of magnitude.

From in vivo studies, Dastiridou et al. estimated ocular rigidity by cannulating the anterior chamber of patients undergoing cataract surgery and recording pressure change for known volume infusion [19]. These results suggest the existence of a negative correlation between the axial length and rigidity, which we also observe from our results, providing further evidence to the validity of our method.

5. Conclusion

We have presented a novel approach to measure choroidal blood flow and ocular rigidity. To the best of our knowledge, this is the first non-invasive method that allows calculation of the true OR parameter as defined by Friedenwald, as it is based on directly quantifying ocular volume changes rather than estimating it based on FPA or a Doppler flowmetry signal. The method relies on measuring IOP and OPA using dynamic contour tonometry, as this technology enables the most precise estimation. We strongly believe the combination of deep penetration dynamic OCT imaging and the powerful automated image segmentation we present is seminal to further understanding of key biomechanical determinants to ocular disease, and will become clinically invaluable as these measuring devices become more accurate.