Real-time correction of geometric distortion artefacts in large-volume optical coherence tomography

Optical coherence tomography (OCT) is a non-contact, non-invasive interferometric technique which allows micrometre resolution imaging of sample microstructures up to several millimetres below the surface of the observed object [1]. The development of spectral domain (SD)-OCT has enabled larger volumes to be scanned at higher speeds [2–3] by treating the interference signal as a modulation in wavenumber k. Large-volume OCT setups typically employ scanning mirrors and suffer from nonlinear geometric distortion artefacts in which the degree of distortion is determined by the maximum angles over which the mirrors rotate. To compensate for this effect, numerous authors have reported the implementation of ray-tracing algorithms [4–7], programmatically applying Fermat's principle. However, as ray-tracing schemes require prior knowledge of the sample microstructure and typically demand high processing power, they are fundamentally unable to apply these corrections in real time. Currently, this causes operators to rely on post-processing to calculate the actual shape of the observed sample. Geometric distortion correction was for instance well described by Zawadzki et al [8] who stated that it took 7–9 min on average to segment the surface and correct the artefacts on a six-core 2.8 GHz processor. The essential advantage of our method is that it is suitable for parallelization and can therefore be implemented in real time using a graphics processing unit (GPU).

A broadband light source, based on a ytterbium-doped fibre (Multiwave Photonics, Porto, Portugal) with a central wavelength at λ = 1050 nm and bandwidth Δλ = 70 nm was used in the SD-OCT setup (figure 1). The optical signal originating from the light source is divided into reference and sample arm by a directional coupler (DC, splitting ratio 80/20). In the sample arm, the beam is collimated by MO1 and diverted towards the sample via the scanning unit GXY and the microscope objective achromatic lens L1 with a working distance of 20 mm.

Zoom In Zoom Out Reset image size

The scanning unit GXY consists of two galvo-scanning mirrors, S_x and S_y, whose respective axes of rotation are perpendicular to each other. Backscattered light from the sample is collected and guided towards the DC where it interferes with the light returning from the reference arm. The interferometric signal is conveyed towards a spectrometer which comprises two achromatic lenses (L2–3), a transmission grating (TG, 1450 lpm, Wasatch Photonics) and an InGaAs line-scan camera (SU-LDH, Goodrich-SUI), operating at 46.9 kHz. Data from the linear CCD are transferred to a PC via a Cameralink cable and a high-speed PCI Express frame grabber (NI-1427 National Instruments). At the working distance of L1, our experimental setup has a lateral field of view in the X- and Y-directions of 12.35 and 10.13 mm, respectively, covering an axial scanning depth of 2.36 mm.

An enlarged detail of the scanning configuration at the end of the sample arm is shown in figure 2 during the recording of a flat mirror (FM). The scanning mirrors S_x and S_y arch the OCT probing beam path, causing the resulting cross-sectional profiles of the FM to become curved. As S_x and S_y reflect the probing beam in succession, their respective beam reflecting points have different optical path lengths to the surface of the sample, and a difference in curvature between the mutually perpendicular transversal scans (or B-scans) B_x and B_y of the same FM can be noticed. Normally, if the distance between the pivot of the galvoscanner and the collimating lens is equal to the focal length of the lens, then the wavefront behind the lens should be flat. However, this distance cannot be adjusted precisely and therefore the wavefront behind the lens deviates from a flat surface. There is also a distance of 1 cm between the pivots of the two mirrors. This is rather large in comparison with the focal length, f = 2.5 cm, of the lens. This means that if the pivot of one galvoscanner is placed at f, then the pivot of the other galvoscanner is at a different distance and the wavefront cannot be flat for both. In this regard, when the OCT setup is configured similarly to the one depicted in figure 2, the amount of angular distortion imposed onto B_y by scanning mirror S_y will be significantly less than that imposed onto B_x by scanning mirror S_x as the origin of distortion of B_y lies further away from the object surface than that of B_x.

Zoom In Zoom Out Reset image size

To calibrate the full 3D stack of data, a series of such full-field OCT images of the FM was recorded at 100 different angles in both scanning directions and with the mirror positioned at several depths within the observed volume. First, threshold filtering was applied to the greyscale OCT images of the observed mirror to infer its topography. Next, the analytical equation of a circle was fitted to these data points using least-squares fitting. In this way, the location of the perceived origin of radial distortion could be determined for each B-scan. The virtual scanning radii R_x and R_y are defined as the shortest distances between the perceived point source of distortion and the surface of the observed FM, in the respective scanning directions. By scanning in the X-direction whilst keeping the Y-scanner fixed and vice versa, the geometric distortion effects imposed onto the OCT volume by the two galvoscanners can be considered to be decoupled. This allows gauging of the distortion characteristics of the two scanners separately and compensating for them sequentially.

As a first step in the data calibration procedure, the 2D Cartesian coordinate points (x, z) on every B_x-scan in the decoupled Y-direction are associated with their polar coordinate system equivalents (R, θ) originating from the virtual origin of distortion (x_c, z_o). After having determined—in post-processing—the full set of radii R(x_c, z) at the central column (x = x_c, z) of each B_x-scan and for every z, the angular component θ(x, z) of the 2D polar coordinates can be calculated per B_x-scan using the basic geometry:

$$\begin{equation} \forall z \left\{ {\begin{array}{@{}l@{}} {R( {x,z}):= R( {x_c ,z})} \\ {\theta ( {x,z}) = {\rm tan}^{ - 1} \left( {\frac{{x - x_c }} {{R( {x_c ,z})}}} \right)} \\ \end{array}}.\right. \end{equation} \tag{ 1 }$$

By assigning the same radius value R(x_c, z) to all coordinate points that have identical z-values, whilst keeping the respective angles θ(x, z) of all coordinate points fixed, the data block is reshaped, so that the distortion effect caused by the X-scanning direction is corrected for, as illustrated in part II of figure 3. This polar coordinate system is reconverted to its Cartesian equivalent by applying the following transformation to all (R, θ)-coordinate couples:

$$\begin{equation} \forall x,z \left\{\begin{array}{@{}l@{}} {x' = R( {x_c ,z}){\rm sin}(\theta ( {x,z}))} \\ {z' = R( {x_c ,z}){\rm cos}(\theta ( {x,z})),} \\ \end{array}\right. \end{equation} \tag{ 2 }$$

creating the Cartesian coordinate data block D(x', y, z') where the data points in the X- and Z-directions have been relocated, so that the distortion artefacts in the X-direction are corrected for. Similar to the above, the set of radii R'(y_c, z') is determined in post-processing for each central column (y = y_c, z') of each B_y-scan for every z'. The angular component φ(y, z') of the 2D coordinate couples (R'φ) can be calculated per B_y-scan as

$$\begin{equation} \forall z'\left\{\begin{array}{@{}l@{}} {R'( {y,z'} ): = R' ( {y_c ,z' } )} \\ \varphi{( {y,z'} ) = {\rm tan}^{ - 1} \left( {\frac{{y - y_c }} {{R'(y_c ,z')}}}\right)}. \end{array}\right. \end{equation} \tag{ 3 }$$

By assigning the same radius value R'(y_c, z') to all coordinate points that have identical z'-values, and by keeping the respective angles φ(y, z') of all coordinate points fixed, the data block is reshaped so that the distortion effects caused by the Y-scanning direction are now corrected for as well, as illustrated in part III of figure 3.

Finally, this polar coordinate system is reconverted to its Cartesian equivalent by applying the following transformation to all (y, z')- coordinate couples:

$$\begin{equation} \forall y,z' \left\{\begin{array}{@{}l@{}} {y' = R'(y_c ,z' ) {\rm sin} (\varphi(y,z'))} \\ {z'' = R'(y_c ,z'){\rm cos}(\varphi (y,z')),} \end{array} \right. \end{equation} \tag{ 4 }$$

creating the final coordinate location points (x', y', z'') that have been relocated to floating point coordinates, so that the distortion artefacts in both scanning directions are now compensated for. After interpolating between these floating point coordinates at the predefined Cartesian coordinate grid $\left( {i,j,k} \right) \in \mathbb{N}^3$:

$$\begin{equation} D( {i,j,k}) = I'_{x,y,z} ,\end{equation} \tag{ 5 }$$

the new intensity values $I^{\prime} _{x,y,z}$ can be correctly displayed.

Zoom In Zoom Out Reset image size

To verify the accuracy of the presented method, an OCT data stack of 12.35×10.13×2.36 mm³ of a human eardrum was obtained and distortion corrected using the above-described calibration procedure. The results of both the corrected and the uncorrected OCT volumes were compared quantitatively against a high-resolution optical moiré topogram. Moiré profilometry is a non-contact optical measurement technique which allows measuring the shape of a 3D surface using projection and triangulation of structured light patterns. The employed moiré setup has a height accuracy better than 15 µm and an in-plane accuracy of 9.35 µm px⁻¹ [9]. The RMS value of the distance between the two obtained surfaces (after alignment) is an indicator of the surface alignment in three dimensions. The proposed calibration algorithm improves the RMSu = 45.24 µm for the uncorrected model to RMSc = 13.49 µm for the corrected model.

Sections through the 3D surfaces are plotted in figure 4. To highlight the deviation from the moiré surface (black full line), a difference plot is included in the bottom graph of figure 4 for both the uncorrected (blue dashed line) and the corrected (red vertically dashed line) models. Near the edges, the uncorrected surface model is misaligned with the moiré model by almost 200 µm, whereas the corrected model misalignment stays within the 20 µm range throughout the full length of the scan.

Zoom In Zoom Out Reset image size

As we have demonstrated that a succession of coordinate transformations suffices to accurately remap the 3D data acquired with a standard OCT setup, the technique described in [10] can be used to create a look-up table of floating point coordinates at which the recorded OCT data blocks need to be interpolated. These calculations can be implemented very effectively on dedicated memory regions in standard GPUs. To assess the processing time delays caused by the proposed technique, we have benchmarked our GPU-based software using our system, containing a 3.33 GHz Intel i5 quad-core CPU and a GTX570 graphics card with 480 CUDA processing cores and 1280 MB of on-board memory. For an imaging volume of 768×768×512 pixels, memory transfer between CPU and GPU takes 437 ms, and the 3D interpolation only adds 10.4 ms to the processing time. Several recent reports have been published on GPU-based processing setups that capture, process and render multiple 3D OCT data blocks per second [11], without any of them performing geometric distortion correction in real time. As we assume that there is no a priori knowledge of the observed sample, the presented method cannot account for errors caused by refraction, absorption or scattering below the surface interface layer. Only the observed surface shape itself can be visualized aberration-free in real time. In this way, the method proves useful for various topographic applications such as corneal topography [12] and tympano-topography [13].

A three-dimensional recalibration procedure was described that effectively compensates for geometric distortion artefacts as they occur in setups where large scanning angles are required to obtain a large field of view. The quality of the proposed 3D recalibration scheme was confirmed by aligning the surface mesh of a highly curved tympanic membrane that was obtained with the described OCT setup with the surface model of the membrane, obtained using optical moiré profilometry. After application of the correction method, the two surfaces coincide within 15 µm. By implementing the proposed interpolation operations on the graphics card of a standard PC, the additional processing time required for distortion correcting large 3D OCT data sets is limited to nearly 10 ms, enabling its real-time implementation and rendering distortion-free surface information.

SVJ and JB acknowledge the support of the Research Foundation—Flanders (FWO). AB and AP acknowledge the support of the ERC grant COGATIMABIO 249889.