## Abstract

Measurement of optic nerve head (ONH) deformations could be useful in the clinical management of glaucoma. Here, we propose a novel three-dimensional tissue-tracking algorithm designed to be used *in vivo*. We carry out preliminary verification of the algorithm by testing its accuracy and its robustness. An algorithm based on digital volume correlation was developed to extract ONH tissue displacements from two optical coherence tomography (OCT) volumes of the ONH (undeformed and deformed). The algorithm was tested by applying artificial deformations to a baseline OCT scan while manipulating speckle noise, illumination and contrast enhancement. Tissue deformations determined by our algorithm were compared with the known (imposed) values. Errors in displacement magnitude, orientation and strain decreased with signal averaging and were 0.15 µm, 0.15° and 0.0019, respectively (for optimized algorithm parameters). Previous computational work suggests that these errors are acceptable to provide *in vivo* characterization of ONH biomechanics. Our algorithm is robust to OCT speckle noise as well as to changes in illumination conditions, and increasing signal averaging can produce better results. This algorithm has potential be used to quantify ONH three-dimensional strains *in vivo*, of benefit in the diagnosis and identification of risk factors in glaucoma.

## 1. Introduction

Glaucoma is the most common cause of irreversible blindness, affecting about 60 million people worldwide [1]. The central event in glaucoma is irreversible damage of retinal ganglion cell axons, which carry visual information from the retina to the brain. The presumed primary location of retinal ganglion cell damage is within the optic nerve head (ONH), the site where retinal ganglion cell axons converge to leave the back of the eye in transit to the brain (figure 1). Glaucoma is diagnosed clinically in patients who have characteristic changes (increased ‘cupping’ or excavation) visible at the ONH (glaucomatous optic neuropathy) associated with typical glaucomatous patterns of visual field loss. It is well established that increased intraocular pressure (IOP) is associated with increased prevalence [3] and incidence [4] of glaucoma. Furthermore, reduction in IOP, either pharmacologically or surgically, is the only treatment proved to slow the progression of the disease [5,6]. In addition, glaucoma occurs nearly as often in patients with normal IOP as in those with elevated IOP, and does so without distinct aetiology [7], although in a significant proportion IOP-lowering still retards disease progression. Finally, some subjects have a statistically elevated IOP (more than 21 mmHg) but will never develop glaucomatous optic neuropathy or visual field loss; these subjects are referred to as ‘ocular hypertensive’.

The above-mentioned observations demonstrate that our current understanding of glaucoma is insufficient: IOP is important, but risk factors other than IOP must play important roles in the development and progression of this pathology. Unfortunately, we currently lack the tools to measure and assess other risk factors.

To explain the above-mentioned clinical observations, it has been suggested that glaucoma patients with normal IOP have an inherently ‘weak’ ONH, whereas ocular hypertensive patients without glaucoma do not develop the disease owing to an inherently ‘strong’ ONH. In other words, the biomechanics of an individual's ONH may dictate the IOP level at which its homeostatic mechanisms are effective [8]. Therefore, above an individual-specific threshold level of IOP, a series of cellular events may be initiated which eventually lead to glaucomatous damage. There is strong evidence to suggest that this is a valid hypothesis [9–26]; however, clear proof has yet to be obtained, particularly as no direct investigation of ONH biomechanics has yet been undertaken *in vivo* in human subjects.

Our long-term research aim is to evaluate the *in vivo* biomechanical characteristics of human subjects' ONHs by imaging deformation of the ONH using optical coherence tomography (OCT). Quantifying biomechanics of the ONH *in vivo* using this novel approach could provide a crucial understanding of the underlying mechanisms involved in glaucoma, and may eventually lead to improved diagnosis of individuals at risk for glaucomatous vision loss. In this study, as the first step towards our long-term research aim, we developed and tested a three-dimensional tracking algorithm that can extract deformations from artificially deformed OCT image volumes of patients' ONHs captured at different IOPs.

The proposed algorithm could have wide applicability and could be used to study other biomechanically related disorders accessible by OCT such as papilloedema [27], atherosclerosis [28,29], keratoconus, macular degeneration and myopia.

## 2. Methods

We wished to develop an algorithm that can track displacements and strains of the ONH following a change in IOP. This algorithm requires two OCT volumes of the ONH: one is captured before a change in IOP, and is referred to as the ‘undeformed’ volume; the other is captured after a change in IOP and is referred to as the ‘deformed’ volume. Each OCT volume is a stack of images but is considered here as a group of voxels (i.e. intensity values are distributed on a three-dimensional rectangular grid).

In the following sections, we first present our three-dimensional tracking algorithm. Its accuracy is tested, using artificially deformed volumes. Its robustness to speckle noise (the noise inherent to OCT) is also tested.

### 2.1. Development of an optical coherence tomography-based three-dimensional tracking algorithm

#### 2.1.1. Principle

Our three-dimensional tracking algorithm was developed using the principle of digital volume correlation [30–34], which has not, to the best of our knowledge, been applied to OCT data previously. In order to track displacements, a region of interest (ROI) is created in the undeformed volume by centring a cuboidal subvolume at each material point where a displacement vector is sought (figure 2). Each ROI (or subgroup of *n _{i}* ×

*n*×

_{j}*n*voxels) is expected to deform with IOP, and its deformation can be mathematically expressed using continuum mechanics concepts as 2.1where

_{k}*x*and

*X*are the position vectors of a material point that belongs to the ROI in the undeformed (reference) and deformed (current) configurations, respectively.

**F**is the deformation gradient tensor (nine unique components), and

*t*is the translation vector (three unique components), both of which are unknown quantities that need to be determined for each ROI. With equation (2.1), we can describe the ROI deformations as a combination of four possible affine transformations (rigid translation with vector

*t*; rigid rotation, stretch and shear with tensor

**F**).

Using this approach, the displacement vector of the ROI centroid, *u*, which is the actual displacement vector we are seeking at each material point, can be expressed as
2.2where *x** and *X** are the position vectors of the ROI centroid in the undeformed and deformed configurations, respectively. By combining equations (2.1) and (2.2), we derive
2.3where **F** and *u* are the quantities to be determined (total of 12 parameters) at each material point.

#### 2.1.2. Global optimization

Extracting **F** and *u* for each material point requires solving a global optimization problem. Put simply, for each ROI, we want to find the ROI deformation (characterized by **F** and *u*) for each couple (*X*, *x*) so that
2.4where *I*_{2} is the voxel intensity function of the deformed OCT volume, and *I*_{1} is the corresponding function of the undeformed OCT volume. Numerically, we chose to minimize the zero-normalized cross-correlation coefficient (ZNCC) [35] between the intensities, which is expressed as:
2.5where
2.6and
2.7where *X ^{ijk}* is the reference position vector of a voxel centroid.

*I*

_{1}(

*X*) is the voxel intensity assumed to be defined at its centroid. In general, the vector

^{ijk}*x*(derived with equation (2.3)) is not required to point to any voxel centroid in the deformed OCT volume.

^{ijk}*I*

_{2}(

*x*) was therefore obtained using a linear interpolation from the intensity values of the eight closest voxel centroids (in the deformed volume) of where the vector

^{ijk}*x*was pointing. Note that the use of such linear interpolation is what allows subvoxel resolution of the proposed three-dimensional tracking algorithm. Therefore, our three-dimensional tracking can deal with conditions where the axial voxel resolution is different from the lateral voxel resolutions without further compensation techniques.

^{ijk}Several other definitions can be used instead of the ZNCC criterion, such as a least-squares criterion. However, the ZNCC criterion has been shown to be very robust to noise and insensitive to any offset or scaling of the illumination conditions [35]—an important aspect for OCT imaging where illumination variations can occur from one acquisition to the next.

In order to minimize equation (2.5), we used a genetic optimization algorithm known as differential evolution [36]. We have used such an algorithm in the past for other global optimization problems such as simple curve fitting [37] and inverse finite element problems [10,11,38]. As control parameters, we used a weighting factor that varied randomly from 0.5 to 1, a crossover constant of 0.9 and a number of trial vectors of 40. Such parameters were chosen using the knowledge we developed by solving previous global optimization problems [10,11,38]. The reader is referred to the work of Price *et al.* [36] for a better understanding of the differential evolution algorithm. To the best of our knowledge, the differential evolution algorithm has not been previously combined with a digital volume correlation algorithm; it is advantageous because it has been proved to be more robust than other genetic and non-genetic optimization algorithms in other contexts [36]. Owing to its demonstrated robustness, we have opted for not performing a pre-search (using translation as the only mode of deformation for each ROI), as other digital volume correlation implementations commonly require [30].

#### 2.1.3. Computational speed improvements

Our three-dimensional tracking algorithm was written in C++ to take advantage of software/hardware performance enhancements. Because tracking the displacement of each point *X* is carried out independently from all other points, our algorithm can be parallelized safely, which was performed using openMP—an open interface for multiprocessing programming. Our simulations were run on a 12-core machine and therefore the tracking of 12 points could be performed simultaneously. We typically required 15 h of computational time to track the displacements of 10 000 points. Note that speed can be affected by factors such as optimization parameters, data noise and ROI size.

#### 2.1.4. Displacement filtering

Although the differential evolution algorithm is robust, it can sometimes fail, resulting in erroneous displacement vectors. To remove such ‘bad’ vectors, we assumed that the ONH was a continuum and we used a displacement filter. Briefly, each displacement vector was compared with its neighbours within a 100 µm radius. Specifically, if the vector's magnitude deviated by more than 15% from the median magnitude in its neighbourhood, or if its direction deviated by more than 15° from the median direction in its neighbourhood, then this displacement vector was considered erroneous and removed from the vector list. These values (15% and 15°) were chosen empirically but represent local deformations that the ONH is unlikely to exhibit *in vivo* [39]. Such an empirical approach was found to be robust for all cases.

Following filtering, displacement vectors were smoothed in three dimensions using a cubic smoothing spline to further remove noise (function csaps in Matlab; smoothing parameter = 0). Vectors removed in the filtering step were replaced by results from such interpolation.

#### 2.1.5. Strain calculations

Following the displacement filtering and smoothing steps, we computed the Green–Lagrange strain tensor at each material point. The effective strain [40], a single index that conveniently summarizes the three-dimensional state of strain, can be extracted as
2.8Here, *E*_{1}, *E*_{2} and *E*_{3} are the principal components of the Green–Lagrange strain tensor [41].

### 2.2. Artificial deformations to verify our algorithm

Our three-dimensional tracking algorithm was applied to three sets of ONH volumes (one undeformed and one deformed per set) in order to test its accuracy. The undeformed OCT volume (same for all sets) was obtained from a baseline OCT scan of a myopic (−6 diopter refraction), non-glaucomatous patient's ONH (size: 768 × 496 × 97 voxels; resolution: 6.25 × 3.87 × 33.3 µm; Spectralis OCT, Heidelberg Engineering, Germany). The corresponding deformed volume was obtained by application of a predetermined, artificial deformation to the baseline scan, using a morphing technique (i.e. inverse mapping). In other words, the baseline volumes are real OCT volumes, but the deformed volumes are simulated. The imposition of known deformations in this way is a standard verification technique, permitting us to assess the accuracy our three-dimensional tracking algorithm and is required before future studies using clinical deformations following a change in IOP. For this section, we considered the following four types of artificial deformations as illustrated in figure 3 as they are consistent with the deformations predicted by biomechanical simulations [19,42–44].

(1) A rigid body motion (translations: 11.1, −22.2 and 33.3 µm along the

*x*,*y*and*z*axes, respectively; rotation: −2.86°, 2.86° and 2.86° about the*x*,*y*and*z*axes, respectively. This deformation is characterized by a mean displacement magnitude of 117 µm and a mean effective strain of 0.(2) Three-dimensional radial expansion of the volume with a 1000 µm radius to mimic eye globe expansion induced by increased IOP. The centre of this latter deformation was ‘virtually’ located 12.5 mm away from the ONH to represent the centre of the eye globe. A rigid translation of 1000 µm along the axial direction (directed towards the centre of the eye) was also applied, so that the volume remains within the field of view. This deformation is characterized by a mean displacement magnitude of 107 µm and a mean effective strain of 0.08.

(3) The same three-dimensional radial expansion of the volume as in the item (2) but with a varying radius from 1000 to 500 µm to mimic increased IOP-induced eye globe expansion accompanied with tissue thinning. A rigid translation of 750 µm along the axial direction (directed towards the centre of the eye) was also applied, so that the volume remains within the field of view. This deformation is characterized by a mean displacement magnitude of 103 µm and a mean effective strain of 0.28.

(4) The same three-dimensional radial expansion of the volume as in the item (3) but with a radius varying nonlinearly (exponential behaviour) from 1000 to 750 µm to provide more realistic tissue thinning (higher in the retina, pre-laminar tissue and choroid than in the sclera and lamina cribrosa [45]). In addition, posterior bowing of the lamina cribrosa (as observed in [23]) was modelled using a three-dimensional Gaussian-type function applied to the centre of the three-dimensional volume (with a variance

*σ*^{2}of 500 µm^{2}and a bowing magnitude*R*_{magnitude}of 150 µm). In other words, a radius of was added to the entire three-dimensional expansion deformation to simulate bowing, where*x*and_{c}*y*are the lateral coordinates of the optic disc centre. Finally, a rigid translation of 875 µm along the axial direction (directed towards the centre of the eye) was also applied, so that the volume remained within the field of view. This deformation is characterized by a mean displacement magnitude of 137 µm and a mean effective strain of 0.19. The performance of our tracking algorithm can then be evaluated by comparing the extracted deformation with the imposed (known) deformation for all four cases._{c}

### 2.3. Modelling optical coherence tomography speckle noise to test robustness

It has been shown that noise in OCT images arises from two main sources: the shot (electronic) noise, which is additive and can be described as a Gaussian white noise, and the speckle noise, which is multiplicative [46] and by far the most important in OCT. In this study, we aimed to evaluate the robustness of our three-dimensional tracking algorithm to speckle noise, which has not been previously undertaken in a digital volume correlation algorithm. Specifically, speckle noise was added to each deformed volume such that
2.9in which *I*_{3} is the intensity of the newly obtained volume (deformed with speckle noise). The *n* values, which are different for each voxel, are randomly generated numbers that follow a negative exponential distribution, i.e.
2.10This model has been shown to represent speckle noise in OCT images relatively well [46] as illustrated in figure 4. In addition, note that this model has no control parameters.

### 2.4. Speckle noise reduction through signal averaging

OCT can acquire data in a remarkably fast manner (from 40 000 A-scans per second for spectral-domain OCT to 100 000 A-scans per second for swept source OCT). Most commercial devices are now able to perform signal averaging (i.e. averaging repeated acquisitions of the same object) in order to reduce speckle noise. In this study, we mimicked signal averaging by
2.11where *I*_{4} is the intensity of a new OCT volume following signal averaging derived from *I*_{2}, *A* is the total number of scans to create an average and *n*_{a}s are random numbers that follow a negative exponential distribution as in equation (2.10). We further tested our three-dimensional tracking algorithm with different levels of OCT speckle noise that can be achieved through signal averaging. We performed 1× (no averaging), 10×, 20× and 30× scan averages (figure 4). These signal-averaging values were chosen, because our group commonly uses them when performing raster scans in a clinical setting.

### 2.5. Modelling optical coherence tomography illumination variations to test robustness

Repeated OCT scans of the same eye may not exhibit identical illumination profiles owing to slight changes in eye orientation with respect to the OCT light source (e.g. subject's eye and head movements). Thus, here we aimed to verify the robustness of our three-dimensional tracking algorithm owing to changes in illumination. To this end, a clinical technician consecutively acquired five ONH volumes of the left eye of the aforementioned subject, over a 30 min time period, with identical acquisition parameters for all scans. All five volumes were then manually realigned. The first ONH volume was considered as the reference volume, from which the voxel intensity of all other volumes was assumed to vary owing to changes in illumination as: 2.12

where *I*_{1} is the voxel intensity of the first ONH volume, and *I _{i}* that of volume

*i. k*is the illumination mask for volume

_{i}*i*, which contains information about illumination variations. Our goal was then to identify the four illumination masks from equation (2.12).

For each ONH volume, we computed the averaged voxel intensity (within 40 × 40 × 40 voxel ROIs) at the same eight corner locations of each ONH volume. Note that averaging was necessary to reduce influence from speckle noise. These experimental data were then used to derive the four illumination mask functions (*k _{i}*(

*X*) from equation (2.12)). Through linear regression, we found that the four illumination masks could be expressed as linear functions of the Cartesian coordinates (

*x*,

*y*,

*z*) (

*R*

^{2}-values of 0.94, 0.72, 0.78 and 0.91, respectively). These derived illumination masks were then applied to each deformed volume to test the robustness of our three-dimensional tracking algorithm to illumination variations (case 4 only).

### 2.6. Shadow removal, increased tissue visibility and contrast enhancement through adaptive compensation

OCT images of the ONH are usually hampered by the presence of shadow artefacts (cast by blood vessels) and by poor tissue visibility in deep tissue layers (i.e. sclera and lamina cribrosa). Both phenomena can impact the quality of the three-dimensional tracking. We have recently developed algorithms (standard [28,29,47,48] and adaptive compensation [49]) that were demonstrated to (i) significantly reduce shadow artefacts; (ii) improve tissue contrast; and (iii) increase the visibility of the deepest tissue layers (see example in figure 5). In this section, we applied adaptive compensation ([49]; contrast exponent = 2; energy threshold = 0.001) to baseline and deformed volumes (case 4 only) to verify the impact of shadow artefacts and low visibility on three-dimensional tracking.

### 2.7. Application, additional sensitivity studies and error definition

For each test case, displacements were tracked at 1210 points in two areas of the ONH (605 points per area that formed a three-dimensional rectangular grid). One area was defined within the optic disc margin and incorporated the lamina cribrosa and the prelaminar tissue. The other area was defined in the peripapillary (outside the disc margin) region surrounding the ONH, including the retina, the choroid and the sclera.

We also tested the performance of our three-dimensional tracking algorithm by varying the ROI size (11 × 11 × 11, 21 × 21 × 21 or 31 × 31 × 21 voxels), and the number of optimization parameters used (with either three to account for rigid translation, six to account for rigid translation and rotation, or 12 to account for rigid translation, rigid rotation, shear and stretch).

For each case, we computed the total error (imposed versus extracted) in displacement magnitude as
2.13where DM_{imposed} and DM_{extracted} are the imposed (three cases above) and the extracted (using our tracking algorithm) displacement magnitudes, respectively. NP is the total number of material points. Similarly, the total error in effective strain was defined as
2.14where *E*_{eff_imposed} and *E*_{eff_extracted} are the imposed and the extracted effective strains. Finally, the total error in displacement direction, which reflects how much displacement vectors deviate in direction, was defined as
2.15where *u*_{imposed} and *u*_{extracted} are the imposed and the extracted three-dimensional displacement vectors. The function atan2 is the four-quadrant arctangent function.

## 3. Results

Total errors in displacement magnitude, displacement direction and in effective strain are depicted as functions of speckle noise amount (decreased through signal averaging) in figure 6 (rigid body motion), figure 7 (three-dimensional radial expansion), figure 8 (three-dimensional radial expansion with tissue thinning) and figure 9 (three-dimensional radial expansion with nonlinear tissue thinning and laminar bowing). For figures 6⇓–8, the effect of varying ROI size (left column) and the total number of optimization parameters (right column) was tested. For figure 9, the effect of varying illumination condition (left column) and applying adaptive compensation (right column) was tested.

### 3.1. Effect of speckle noise

When no speckle noise was added to the deformed volume, we found that our three-dimensional tracking algorithm performed extremely well with a best-case total error in displacement magnitude of 0.0004 µm, in displacement direction of 0.0002° and in effective strain of 0.000002; with added speckle, noise errors were markedly increased, to 4.2 µm, 4.5° and 0.06, respectively.

### 3.2. Effect of signal averaging

Overall, we found that signal averaging greatly improved the accuracy of the tracking algorithm as total errors decreased consistently with increased signal averaging in most cases (figures 6⇑–8). We also observed that 10× signal averaging usually produced considerable improvements with respect to no signal averaging. For 10× signal averaging, best-case total errors were 0.45 µm, 0.36° and 0.0042. Further signal averaging (20× and 30×) produced quasi-linear and smaller improvements consistently. For 30× signal averaging, best-case total errors were reduced to 0.15 µm, 0.15° and 0.0019 (figure 10).

### 3.3. Effect of illumination variations

From five successive acquisitions of the same subject's ONH, we found that tissue reflectivity (or OCT pixel intensity) varied on average by only 5.7% through comparison of the illumination masks. Such changes were almost indiscernible by visual examination of the OCT images. Overall, we found that the three-dimensional tracking algorithm was insensitive to varying illumination masks as shown in the left column of figure 8; furthermore, this conclusion held for any level of speckle noise that we tested.

### 3.4 Effect of shadow removal, increased tissue visibility and contrast enhancement

We found that application of adaptive compensation to the baseline and deformed ONH volumes was able to improve three-dimensional tracking (smaller errors in displacement magnitude/direction and effective strain; figure 9) when no signal averaging was performed. When signal averaging was used, we found that adaptive compensation gave slight improvements in displacement magnitude error, but slight impairments in displacement direction and effective strain errors (figure 9).

### 3.5. Effect of region of interest size

In the presence of speckle noise without averaging, increasing the ROI size could reduce total errors with deformation cases 2 and 3 (figures 7 and 8), whereas the opposite was observed with deformation case 1 (figure 6). In the presence of both speckle noise and signal averaging, increasing the ROI size yielded either no or slight improvements for all deformation cases (figures 6⇑–8).

### 3.6. Effect of optimization parameters

When using only three optimization parameters (i.e. allowing each ROI to exhibit a rigid translation deformation only), we found that the accuracy of our algorithm was the poorest (figures 6⇑–8). Best-case total errors were 2.4 µm in displacement magnitude, 1.8° in displacement direction and 0.02 in effective strain.

Increasing the optimization parameters from 3 to 6 (i.e. allowing each ROI to exhibit rigid translation and rotation, but no deformation), resulted in substantial improvements in deformation case 1 (figure 5) but not in deformation cases 2 and 3 (figures 7 and 8). Such behaviour is readily understandable as case 1 is itself a rigid body transformation (translation + rotation), whereas cases 2 and 3 do not contain rigid rotations.

Likewise, increasing the optimization parameters from 6 to 12 (i.e. using rigid translation, rigid rotation, stretch and shear as the permissible deformations for each ROI) resulted in major improvements in deformation cases 2 and 3 (figures 7 and 8) but not with case 1 (figure 6). This is to be expected as case 1 contains neither shear nor stretch, unlike cases 2 and 3.

## 4. Discussion and conclusion

In this study, we have developed an OCT-based three-dimensional tracking algorithm that could be used to map *in vivo* deformations of ONH tissues following a change in IOP. We have tested its accuracy using OCT images with predetermined and artificial deformations and assessed its robustness to speckle noise and illumination. Our algorithm was found to be accurate and can extract three-dimensional displacements with errors on the order of 1 µm, well below the OCT lateral resolutions used here (6.25 and 33.3 µm). Furthermore, in the presence of speckle noise, signal averaging was shown to produce excellent results. Changes in illumination conditions had negligible impact on tracking performance. We are now applying this algorithm to the tracking of *in vivo* deformations of the ONH following surgical IOP-lowering, and IOP increase *in vivo* and *ex vivo*.

The ability of our algorithm to track known displacements was found to be excellent with best-case total errors in displacement magnitude of 0.0004 µm, in displacement direction of 0.0002°, and in effective strain of 0.000002. When speckle noise was added to the images, our algorithm was still able to track displacements (figure 10); however, its performance was at its worst. Signal averaging, a common feature of most OCT scanners, was modelled and was shown to drastically improve the robustness of our algorithm against speckle noise. With 10× signal averaging, best-case total errors were 0.45 µm (displacement magnitude), 0.36° (displacement direction) and 0.0042 (effective strain). We currently observe strain levels exceeding 0.1 following IOP-lowering by trabeculectomy in glaucoma subjects [50], which suggests that the strain error in our algorithm will be acceptable with at least 10× scan averaging.

Total errors were further decreased with signal averaging. For 30 × signal averaging, best-case total errors were 0.15 µm, 0.15° and 0.0019, respectively. Overall, signal averaging produced nonlinear improvements, with large improvements observed with 10× signal averaging only, and smaller improvements observed with more than 10× signal averaging. These results suggest that high amounts of signal averaging (more than 10×) may not be required to obtain the necessary accuracy for quantification of ONH biomechanics. This is an advantage as an increase in signal averaging results in a proportional increase in OCT acquisition time.

We found that our three-dimensional tracking algorithm was extremely robust to changes in illumination conditions (figure 9). These changes were estimated and modelled using repeated scan acquisitions of the same subject's ONH. These results are consistent with the use of the ZNCC criterion, which is insensitive to any offset or scaling of the illumination conditions. It should also be emphasized that most commercial OCT systems will only acquire data if relatively constant illumination is achieved within the entire image. Thus, illumination variations should not pose a problem to three-dimensional tracking as long as longitudinal acquisitions are performed by suitably trained clinical technicians.

Our algorithm was also found to be more robust when adaptive compensation was applied to both the undeformed and deformed volumes (no speckle noise averaging; figure 9). This is likely due to the fact that contrast is enhanced and tissues are more visible through the compensation step, thus helping the global optimization procedure. However, improvements were not consistent when signal averaging was first applied (i.e. better error in displacement magnitude, but worse errors in displacement direction and effective strain). Because adaptive compensation has a tendency to increase noise (owing to signal amplification), it is likely that this compensated generated noise could affect tracking performance and a trade-off between compensation and signal averaging is observed.

It should be emphasized that our algorithm is able to track displacements with subvoxel resolution. For instance, when using signal averaging, most total errors were less than 1 μm even though OCT axial and lateral resolutions were 3.87, 6.25 and 33.3 µm. This is a considerable advantage comparing with manual displacement tracking methods used by other investigators [45,51] that cannot provide subvoxel accuracy. Because IOP-induced displacements of the ONH tissues are expected to be small [52,53], manual tracking techniques are unlikely to provide the necessary resolution required for proper biomechanical quantification of the ONH.

In this study, we also found that using 12 optimization parameters was necessary in order to obtain optimal tracking performance. For example, total error in effective strain was reduced from a high of 0.06 (three optimization parameters) to 0.0019 (12 optimization parameters) for case 3 (with 30× signal averaging). Because small IOP-induced deformations are expected in the ONH, in order to minimize the error, we strongly recommend the use of 12 optimization parameters rather than only 3, 6 or 9 as done in several studies that implemented digital volume correlation [32,34,54].

We also found that increasing the ROI size could further reduce total errors when signal averaging was used. However, the improvements were relatively small, suggesting that using a large ROI may not be necessary for three-dimensional tracking. This is an advantage as the use of smaller ROIs can result in a considerable increase in computational speed.

Several limitations exist in our study. First, we used only affine transformations to represent local deformations of the ONH tissues (i.e. local ROI deformations) rather than more complex deformations (e.g. quadratic) [33]. However, because each ROI is small with respect to the OCT volume, affine transformations at each ROI can be sufficient to fully characterize the ONH tissue deformations.

Second, we have evaluated only the algorithm's robustness to speckle noise, the dominant noise in OCT, and to changes in illumination conditions. In previous preliminary studies, we found that the shot noise (or Gaussian white noise) had a very small impact on three-dimensional tracking (data not shown). This is consistent with the fact that shot noise is usually influential when signal levels are very low. This is not the case for the ONH, as most of the ONH layers of biomechanical importance in glaucoma can be seen as bright or very bright. Other noise sources may include noise due to digitization and thermal detector noise, but like shot noise, their influence is usually small. Noise can also come from patients’ head movements and OCT registration error. However, the current OCT system we are using (Spectralis OCT from Heidelberg Engineering) is equipped with an additional imaging modality (confocal scanning laser ophthalmoscopy) that is used to reference the OCT images and allows ‘true’ anatomy measurements despite patient head movements. Potential changes in tissue optical properties from one OCT acquisition to the subsequent one could also be seen as a source of error; however, no study has yet reported reflectivity changes in glaucomatous connective tissues. Finally, birefringence artefacts [55] could potentially cause a problem, but they also have yet to be characterized properly.

Third, it is possible that errors under ‘real’ clinical deformation conditions could prove higher than those indicated here. As discussed, strain and displacement errors could arise from speckle noise, differing illumination profiles, light attenuation, changes in optical properties, birefringence artefacts, shot noise, image registration, patient movements or even from complex interactions among all these factors. While we have tried our best to verify the accuracy of our algorithm under multiple conditions, a validation step is still required for clinical acceptance of the methodology. Unfortunately, to date, there are no alternative methods to experimentally map displacement and strain within the ONH tissues, thus hampering direct validation of the proposed method. Further validation of the method will include manual tracking of *in vivo* deformations, or the use of alternative technologies when those become available. It is worth adding that improvement in *in vivo* ONH imaging resolution from adaptive optics [48] or μOCT [56] technology has strong potential to provide further credibility to three-dimensional tracking—a direction that we are currently investigating.

Our three-dimensional tracking algorithm is now being tested with OCT data captured in patients at different levels of IOP [50]. Our goal is to verify whether the biomechanical properties of each individual tissue layer of a patient's ONH can indicate susceptibility to glaucoma, and whether a correlation exists between IOP-induced strain and loss of visual function. Our long-term research goal is to develop tools to quantify the biomechanics of the ONH *in vivo*; this may have a strong potential to be used for the diagnosis of ONH disorders such as glaucoma.

The views expressed are those of the authors and not necessarily those of the NHS, the NIHR or the Department of Health.

## Funding statement

This work was supported by the Ministry of Education, Academic Research Fund, Tier 1, Singapore (Girard), by an Imperial College Junior Research Fellowship, London, UK (Girard), by a Royal Society Wolfson Research Merit Award (Ethier) and by the National Institute for Health Research (NIHR) Biomedical Research Centre based at Moorfields Eye Hospital NHS Foundation Trust and UCL Institute of Ophthalmology (Strouthidis).

- Received May 20, 2013.
- Accepted July 2, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.