1. Introduction

Spectral domain optical coherence tomography (OCT) is a volumetric imaging modality with broad clinical and research applications [1]. OCT images are built of A-Scans – one-dimensional reflectivity maps along the device’s optic axis. Formation of an image or volume requires sweeping the beam, acquiring A-Scans on a grid of points. As a result, large volumes built of many A-Scans take a substantial amount of time to acquire, and high SNR can only be achieved with averaging that further increases acquisition time. This leaves the imaging vulnerable to issues such as sample drift, or issues related to non-stationarity [2,3].

Compressed sensing imaging (CSI) has been developed to combat the acquisition time problem for images and volumes for many modalities by reducing the number of measurements required to achieve acceptable resolution [4–6]. Although CSI was originally formulated for acquisition in a transform domain incoherent to the imaging domain (e.g. Gaussian random matrices), “compressed sensing" has since been used more generally as a term encompassing sparsity-based signal recovery [7]. In this context, CSI has seen success in 3D volumetric OCT imaging [8,9] by employing either random or uniform image-domain subsampling schemes.

In addition to imaging, spectral domain OCT is also capable of vibrometry through spectral domain phase microscopy (SDPM) [10]. The A-Scan is generated by a Fourier transform of raw photodetector data, and its magnitude is proportional to the sample reflectivity (this is what is used for imaging), while the phase of the A-Scan is proportional to sub-pixel displacements. Taking many A-Scans at a single location over a period of time and tracking the phase, one can measure displacement responses with resolution on the order of an angstrom. This makes OCT particularly well-suited for the study of auditory mechanics, wherein structures in the inner ear all move at angstrom-to-nanometer scales to produce the sensation of hearing [11–17]. The application of compressed sensing to OCT-based displacement measurements has not yet been considered.

OCT vibrometry is significantly more time-consuming than imaging for two reasons: 1) measuring a displacement response requires acquiring many A-Scans at a single position over a period of time, and 2) displacements measured along an A-Scan are one-dimensional projections of the true three-dimensional motion onto the optical axis [10,14,18]. This second point means that complete description of three-dimensional motion can only be achieved if motion is measured at the same point at three or more angles, and then reconstructed from the resultant one-dimensional measurements [17,19,20]. This emphasizes the importance of developing compressed sensing vibrometry (CSVi).

Such a development is also timely; studies of recent interest in the field of cochlear mechanics use OCT to measure displacement responses at many locations in the cochlea over the course of a single experiment [14,18,19,21,22]. In in vivo experiments, time is restricted by the need to maintain animal anesthesia. This heavily limits the resolution and scale with which displacement maps can be acquired. CSVi could significantly reduce the amount of time required for such experiments, allowing for high resolution motion map acquisition, and opening the door to more complex experiments. For example, over-constrained 3-D motion reconstruction [19] or motion reconstruction in time-sensitive perturbation studies [13,21] are currently intractable, but could conceivably be performed in a reasonable amount of time using CSVi.

Recent years have seen a rise in deep-learning based image restoration and compressed sensing with state-of-the-art results [23,24]. These methods by and large rely on vast amounts of training data, with corresponding ground-truth acquisitions, and are often unable to perform reconstruction on samples outside of their training distribution (e.g. a change in instrument or subject) [25]. Optimization-based reconstruction methods are thereby favorable to learning based methods for the present application, as large in vivo OCT displacement map datasets are currently unfeasible to acquire. Optimization-based techniques remain especially popular in medical imaging contexts, wherein Total Generalized Variation regularization (see Section 4) and its derivatives comprise many state-of-the-art non-learned algorithms [26–28].

We provide a theoretical development of CSVi on OCT displacement scans, and provide evidence of its use in cochlear mechanics. We find that our method, based on a uniform subsampling pattern and total generalized variation optimization, can reduce acquisition time of dense areal motion maps in the cochlea by a factor of ten with less than 5% error. In addition, we argue that our method can be interpreted as a simultaneous denoiser by penalizing aphysical high-frequency components of displacement maps. This method is robust to stimulus frequency and level, as well as beam axis orientation with respect to the sample anatomy.

2. Background

2.1 Time-dependent OCT signals

Figure 1 shows the fundamental objects of interest for this study – OCT scans of the cochlea along with lineal and areal motion measurements. Panel A shows an anatomical drawing of a cross-section of the gerbil organ of Corti complex (OCC), the organ in the cochlea responsible for the transduction of fluid pressure stimuli into neurotransmitter.

The building block of OCT signals is the A-Scan – a complex-valued one-dimensional scan along the optic axis $\hat {\mathbf {z}}$ pointing into the sample. The magnitude of the A-Scan is related to the reflectivity of the sample at each depth, while the phase contains information about sub-pixel displacements. The nature of OCT is such that the values at all $z$-positions in an A-Scan are recorded simultaneously.

One can form a B-Scan by taking multiple A-Scans along a line segment. The magnitude of a B-Scan is a 2-D image, where pixel values encode sample reflectivity. Panel B in Fig. 1 shows a sample B-Scan of the OCC, along with a copy of this B-Scan with overlain markers of the anatomical structures derived from known anatomy (panel A). Shown also is the magnitude of a sample A-Scan, which makes up one vertical line in the B-Scan marked by a white dashed line.

The axis along which A-Scans are swept to form a B-Scan is called $\hat {\mathbf {y}}$, so that a B-Scan occupies the $yz$-plane. A C-Scan, or volumetric scan, is formed by taking parallel B-Scans along a third orthogonal axis, $\hat {\mathbf {x}}$. The coordinates $\hat {\mathbf {x}}$, $\hat {\mathbf {y}}$ and $\hat {\mathbf {z}}$ are referred to as optical coordinates, and they change dependent on the position and orientation of the scanner. We write positions in optical coordinates as $\mathbf {p}^\theta$, where $\theta$ is an index that describes the orientation of the scanner.

 

Fig. 1. Examples of OCT A-Scans, B-Scans and displacement measurements in the gerbil organ of Corti complex (OCC). A – Anatomical drawing of a cross-section of the OCC in gerbil, with important structures labelled. Axes at the top-right denote the longitudinal, radial and transverse axes of the cochlea. B – B-Scan of the gerbil cochlea taken in vivo through the round window membrane (RWM), 430 ${\mu}$m $\times$ 300 ${\mu}$m. The unlabeled B-Scan is recorded, and points are related to structures according to known anatomy (see panel A). The labeled B-Scan shows the RWM, basilar membrane (BM), outer hair cells (OHC), inner hair cells (IHC), pillar cells (PC), inner sulcus space (ISS), outer tunnel (OT), tectorial membrane (TM) and Reissner’s membrane (RM). Each vertical line in the B-Scan is an A-Scan magnitude. A white dashed line denotes the position of a sample A-Scan, shown to the right. C – OCT-measured time-domain displacement at the OHC in response to a Zwuis stimulus at 80 dB SPL, in vivo. D – Frequency-domain gain re ear canal (EC) pressure at the OHC in vivo, in response to Zwuis stimuli at 60, 70 and 80 dB SPL. Only data that are at least 2.5 standard deviations above the noise level are shown. E – An areal map of the 25 kHz component of the OHC displacement magnitude in vivo, in response to an 80 dB SPL Zwuis stimulus. Scan is 200 ${\mu}$m $\times$ 200 ${\mu}$m.

Download Full Size | PDF

Previous CSI work has focused on C-Scan magnitudes [8,9]. CSVi, on the other hand, concerns the phase. In SD-OCT these features are acquired simultaneously, but phase information can only be interpreted through acquisition of time-series of A-scans.

An M-Scan is formed by repeatedly taking A-Scans at a single location $(x,y)$ over a period of time. We write the continuous M-Scan in polar form as

The theory of SDPM gives a formula for the average sub-pixel displacement in terms of the phase of the M-Scan [10]. So long as the structures within a pixel are moving with displacements smaller than the pixel size (this assumption is always valid in cochlear mechanics, where displacements are usually on the order of nanometers), the displacements of the structures along the beam axis are given by

Mathematically, $\delta ^\theta$ is a scalar field mapping points in 3-D space and time to 1-D real-valued displacements:

This displacement is a 1-D projection of the true motion of the structure onto the optical $\hat {\mathbf {z}}$-axis. Generally speaking, a structure at a single point in space will have a 3-D displacement – that is, displacement is completely described by a vector field mapping points in space and time to a 3-D displacement vector:

It is prescient to discuss this object in terms of a global coordinate system rather than optical coordinates, as we will later consider the structure as being measured from more than one angle. Defining a general spatial coordinate $\mathbf {r}\in \mathbb{R}^3$, we can write the true displacement as

The measured displacement and true displacement are related by a projection map $\mathbf {P}^{\theta }$ that projects 3-vectors onto the beam axis. Similarly, the optical and global coordinates at orientation $k$ are related by a change of basis matrix $\mathbf {R}^{\theta }$ and a translation $\mathbf {t}^\theta$ (as any two stationary coordinate systems in $\mathbb{R}^3$ are related). We write

2.2 Frequency-dependent OCT signals

More often considered are the frequency domain versions of these displacement signals, which offer more clear information about the behavior of a system such as the cochlea. We define the temporal Fourier transforms of the true and measured displacement signals as

Under an assumption of stationarity, these signals describe a pattern of simultaneous motion within the sampled structure. This stationarity assumption is generally justified so long as no sample deterioration or drift has occurred over the sampling period. This further motivates the developing of CSVi, as accelerated sampling time better ensures the accuracy of this assumption.

The frequency axis is discretized by the character of the sound stimulus. Displacements within the cochlea are usually measured in response to either a sweep of pure tone stimuli, or a multitone “Zwuis" stimulus [29]. In either case, a finite number of frequency components are presented, and we observe the displacement responses only at those stimulus frequencies.

For our specific application, we must also note that the cochlea is a nonlinear system – displacement responses differ in character dependent on the stimulus magnitude. Magnitude is referred to as sound pressure level (SPL). Its units, dB SPL, are defined as

Figure 1 panel D shows an OCT-measured frequency-domain response of the outer hair cells (OHC) to a Zwuis stimulus at 60, 70 and 80 dB SPL. The data are plotted as gains with respect to the pressure at the ear canal (EC), which is where the input stimulus is applied – that is, this resembles the “transfer function" of the system (of course, nonlinear systems do not actually have well-defined transfer functions). As the gain varies as a function of SPL, it is clear that the system is nonlinear.

The signal we wish to reconstruct is $\mathbf {D}$. This must be done by measuring from $\Theta \geq 3$ orientations [17,19,20]. At each $\theta =1,2,\ldots \Theta$, we define $\hat {{\Delta} }^\theta$ to be the sampled version of the continuous complex scalar field ${\Delta} ^\theta$. Supposing that at each orientation we sample at $F$ frequencies at $S$ SPLs along an $M\times N$ grid, and supposing that A-Scans contain $L$ pixels along the $z$ axis, the sampled signal is

It should be noted that the process described in the present work is equally valid for areal scans (2-D scans where $N=1$). While we will develop the most general 3-D case, the focus of the results will be on areal scans, which have been the subject of recent inquiry in cochlear mechanics [14,21]. Figure 1 panel E shows an example of the magnitude of such an areal motion map from the gerbil cochlea at 25 kHz and 80 dB SPL.

2.3 Dimensions of sampling and required generality

Equation (9) shows that there are five axes on which we sample during cochlear mechanics experiments – the $x$, $y$ and $z$ spatial coordinates, frequency and SPL. Our goal is to reduce the number of samples taken at experiment time, but it is not necessary to subsample along all of these axes.

As for frequency, all $F$ frequency components are presented simultaneously if Zwuis multitone stimuli are used [29]. This means that sampling in frequency is not a time-limiting factor, so subsampling in frequency is unnecessary.

It is not standard in cochlear mechanics to present data as a function of SPL – generally, frequency-domain data is presented for a select hand-full of SPLs. SPL sampling is already relatively sparse, so its subsampling need not be considered.

As for the $z$ axis, the nature of OCT is that the displacements at all $z$ positions at a single $(x,y)$ location are recorded simultaneously. This means that $z$ sampling does not impact recording time, so we need not consider reducing the number of samples in $z$.

This leaves the $x$ and $y$ spatial domains in which to subsample, i.e. we look to compress the grid along which the M-Scans are taken. This means reducing the values of $M$ and $N$ in Eq. (9). At each SPL-orientation-frequency permutation, this is analogous to CSI for a complex volume.

Although orientation, frequency and SPL do not directly impact our subsampling scheme, they are still critical to consider in the validation of the CSVi method. Motion maps differ significantly between stimulus frequency, SPL and viewing angle [14,19,21,22]. Our method for CSVi must be sufficiently general to perform well for all frequency-orientation-SPL combinations.

3. Algorithms

We frame the compressed sensing vibrometry (CSVi) problem as complex-valued signal recovery problem. For notation, we consider our images in vectorized form, i.e. an $N_1 \times N_2$ complex-valued displacement map is a vector $\boldsymbol{x} \in \boldsymbol{C} ^N$, with $N_1 \times N_2 = N$. Our observations are modeled as $\boldsymbol{y} = \boldsymbol{M} \boldsymbol{x} + \boldsymbol{\nu} \in \boldsymbol{C} ^{N}$, where $\boldsymbol{M} = \mathbf{diag}(\boldsymbol{m} ) \in \{0, 1\}^{N \times N}$ is a mask operator, $\boldsymbol{m}$ is a boolean vector which denotes if a pixel in the displacement map was collected, and $\boldsymbol{\nu}$ is unknown measurement error. Note that for CSVi $\boldsymbol{m}$ is a boolean image with columns removed as we consider a subsampling of A-scans only.

In the following subsections, we present signal models and associated algorithms for the recovery of displacement map $\boldsymbol{x}$ from a masked observation $\boldsymbol{y}$ via use of sparsity-based priors/regularizers. We choose to use priors that have seen success in reconstruction problems involving anatomical images. While we are not reconstructing the anatomy but rather the displacement map, these maps are constrained by the morphology of the OCC. Moreover, observed and expected smoothness of these maps suggests that gradient- and wavelet-based sparsity are reasonable priors [14,21].

The following algorithms perform optimization of real-valued objective functions of complex variables by considering the optimization jointly over real and imaginary parts of their argument. For neatness of presentation, this joint optimization is written with complex arithmetic. For example, the simple linear program over two variables

3.1 Iterative shrinkage thresholding algorithm

In this section, we consider the signal prior of sparsity in a wavelet domain by performing a change of variables. Let $\boldsymbol{x} = \boldsymbol{W} \boldsymbol{z}$, where $\boldsymbol{W}$ is an (orthogonal) wavelet basis. This basis may be understood as representing a signal in terms of a single low-pass subband at a coarse image resolution, and several directional band-pass subbands at different resolutions. We then seek the wavelet coefficients $\boldsymbol{z}$ by using the sparsity-promoting $\ell _1$ norm as a regularizer, often referred to as the Basis Pursuit DeNoising problem (BPDN),

Here, $\boldsymbol{\lambda}$ is a vector of hyperparamters for determining tradeoff between data-fidelity and sparsity in each subband of the wavelet transform. The element of $\boldsymbol{\lambda}$ corresponding to the low-frequency component of the wavelet coefficients is set to zero as there is no reasonable assumption of sparsity.

A popular algorithm for solving BPDN is proximal gradient descent (PGD), which tackles objectives of the form $f(\boldsymbol{x} ) + g(\boldsymbol{x} )$ with $f$ smooth and $g$ possibly non-smooth, involving the proximal operator $g$, defined as,

When $g = {\cdot }_1$, the proximal operator has a closed form solution known as element-wise Soft-Thresholding,

Algorithm 1. Iterative Shrinkage Thresholding Algorithm (ISTA) for (10)

View Table | View all tables in this article

For our ISTA implementation, we chose to use a three-level Daubechies-7 wavelet transform. This choice was motivated by assessment of the sparsity of our data in various wavelet domains, as well as evaluation of the performance of ISTA when using 2-, 3- and 4-level wavelet transforms. More details can be found in the Supplemental Material.

3.2 Total variation

An alternative signal model is that of Total Variation (TV) regularization, in which the spatial gradient of the displacement map (in real and imaginary parts) is assumed to be sparse. We express this formulation as,

The $\ell _{1,2}$ norm, $ \|{\boldsymbol{D} \boldsymbol{x} }\|_{1,2} := \sum _{i,j} \|{\begin {bmatrix} (\boldsymbol{D} _v\boldsymbol{x} )[i,j] \\ (\boldsymbol{D} _h\boldsymbol{x} )[i,j] \end {bmatrix}}\|_2$, imposes an isotropic regularization by encouraging joint-sparsity in horizontal and vertical directions.

A efficient method for solving (13) can be found in Primal Dual Splitting (PDS, also known as the Chambolle-Pock Algorithm or the Primal Dual Hybrid Gradient Method) [34]. PDS tackles problems of the form

The PDS algorithm can then be written as taking alternating proximal gradient descent and proximal gradient ascent steps on the primal and dual variables, $\boldsymbol{x}$ and $\boldsymbol{z}$, respectively,

Returning to the CSVi TV regularization problem, our saddle-point problem is,

The proximal operator of the data-fidelity term can be found in closed form as,

The proximal operator of the indicator function of a convex set $\mathbb{C}$ is given by projection onto said set, $\Pi _{\mathbb{C}}$. Together, PDS for the TV-regularized CSVi problem is given in Algorithm 2.

Algorithm 2. Primal-Dual Splitting (PDS) for (22)

View Table | View all tables in this article

4. Total generalized variation

TV regularization assumes a piecewise constant signal model, which is ill-suited for describing the complex internal motions of the organ of Corti. A more sophisticated signal prior known as Total Generalized Variation (TGV) [34], in which the sparse spatial gradient prior is proposed in a hierarchical manner, offers a greater ability to model cochlear displacement. In the following, we consider the second order TGV prior (TGV-2). We express this formulation as,

The ratio of $\lambda _0$ to $\lambda _1$ play a role in balancing the importance placed on first and second spatial derivatives of the displacement map in the objective function, respectively.

A saddle point problem amenable to PDS can be formed as

The corresponding PDS algorithm for CSVi with a TGV-2 prior (26) is given in Algorithm 3. This algorithm is derived in an identical manner to the TV PDS algorithm (Alg. 2), wherein proximal gradient descent is performed on the primal variables ($\boldsymbol{x}, \boldsymbol{v}$) with dual variables ($\boldsymbol{z}, \boldsymbol{q}$) fixed, then proximal gradient ascent is performed on the duals with primals fixed.

Algorithm 3. Primal-Dual Splitting (PDS) for (26)

View Table | View all tables in this article

The presented algorithms were implemented in the Julia programming language and run on an Intel i5-8250U 3.4 GHz CPU with 8 threads. For a single reconstruction, ISTA and TV algorithms complete in roughly 45 seconds, and TGV completes in roughly 90 seconds.

5. Methods

5.1 Experimental preparation

The experiments were approved by the Columbia University Institutional Animal Care and Use Committee. Young adult gerbils were anesthetized, their scalps were removed and their heads were attached to a two-axis goniometer. Measurements were taken at the base of the cochlea through the left ear’s round window membrane using a Thorlabs Telesto 311C spectral domain OCT device, equipped with an LSM-03 lens. The device has a 1300 nm center wavelength, an axial resolution of 8 ${\mu}$m and a lateral resolution of 10 ${\mu}$m. Acquisition and initial processing of the OCT data was performed using the Thorlabs SpectralRadar software development kit. The temporal sampling rate of M-Scan recordings is approximately 100 kHz.

The cochlea was stimulated via one-second Zwuis multi-tone stimuli [29], generated by a Tucker Davis Technologies system and delivered closed-field to the animal’s EC using a plastic tube. The EC pressure, used as a reference for gain and phase responses, was measured using a Sokolich ultrasonic microphone with its probe tube stationed 1-2 mm from the tympanic membrane. Further details of the animal preparation and acoustical system can be found in Strimbu et al., 2020 [13].

The in vivo data set consists of 275 areal motion maps from three animals taken at very distinct orientations relative to the cochlea’s anatomy, together representing the three most common orientations considered in cochlear mechanics experiments. They are:

$\theta = 1$: Data taken near the 24 kHz region taken along the longitudinal axis of the cochlea through the OHC region (less technically, along the cochlea’s spiral orthogonal to the B-Scan in Fig. 1 B). The cochlea was stimulated with Zwuis stimuli containing 25 frequency components between 2 kHz and 30 kHz at 60, 70 and 80 dB SPL. That is, there are 75 maps from this animal. Each map is 270 ${\mu}$m $\times$ 300 ${\mu}$m, containing 100 rows and 200 columns.

$\theta = 2$: Data taken near the 24 kHz region along the radial axis at an oblong angle with respect to the transverse direction of the cochlea (less technically, a skewed version of the anatomical cross-section in Fig. 1 A). The cochlea was stimulated with Zwuis stimuli containing 25 frequency components between 2 kHz and 30 kHz at 50, 60, 70 and 80 dB SPL. That is, there are 100 maps from this animal. Each map is 270 ${\mu}$m $\times$ 300 ${\mu}$m, containing 100 rows and 200 columns.

$\theta = 3$: Data taken near the 40 kHz region in a transverse-radial plane (less technically, the exact cross-section represented in Fig. 1 A). The cochlea was stimulated with Zwuis stimuli containing 25 frequency components between 2 kHz and 50 kHz at 50, 60, 70 and 80 dB SPL. That is, there are 100 maps from this animal. Each map is 270 ${\mu}$m $\times$ 330 ${\mu}$m, containing 100 rows and 330 columns.

5.2 Pre-processing

Data were pre-processed on an M-Scan by M-Scan basis, as neighboring M-Scan information is not available in a sparsely sampled signal (i.e. areal maps are pre-processed column-by-column). The noise floors of the M-Scans were estimated, and M-Scan points were replaced with 0 if they were less than 2.5 standard deviations above the noise floor. The M-Scans were then median filtered (three-pixel kernel, where the median is selected by magnitude) in the $z$ dimension to remove outliers.

5.3 Evaluation

We look to compare the success of TV, ISTA and TGV using either uniform or random subsampling. We consider three subsampling rates, where the number of samples is reduced by a factor of $P=2$, $5$ or $10$. For random subsampling, we average over ten realizations of the sampling mask.

We begin by considering the performance of these six methods on a test set of 20 maps chosen from the full dataset at random. We then measure the performance of the best tested method on the entire dataset of 275 maps.

The reconstructions are compared on two metrics: normalized mean square error (NMSE) between the complex reconstruction and densely sampled signals, and structural similarity index (SSIM) between the magnitudes of the reconstruction and the densely sampled signals. For each tested displacement map, reconstruction method, subsampling paradigm and subsampling rate, a grid search was used to determine the parameters which give the optimal NMSE, and these values were used to compute the reported NMSE and SSIM.

5.4 Visualization of areal motion maps

The complex-valued areal motion maps in Figs. 3 and 5 are represented using cyclic colormaps with varying saturation. Colorwheels are present in the first panel of both figures.

The hue, representing the phase re EC pressure at each pixel, varies moving counter-clockwise around the colorwheel in the standard manner. That is, the hue at the right-hand horizontal represents a point moving in-phase with EC pressure, while a point at the top vertical represents a point +0.25 cycles out of phase with EC pressure.

The saturation represents gain magnitude at each pixel normalized to the maximum. It varies linearly from the center (black, no motion) to the outer edge of the colorwheel (most saturated, maximum gain magnitude across the areal map).

6. Results

6.1 Comparison of methods

We first compared TV, ISTA and TGV over a 20-map test set as described in Sec 5.3. The mean NMSEs and SSIMs, along with the respective standard deviations, across this test set are shown in Fig. 2.

The results show that across this test set, uniform sampling consistently provides better results than random sampling on both the NMSE and SSIM metrics for all tested methods and subsampling rates. Across the methods, TGV performs better than TV and ISTA. Within one standard deviation of the mean for this test set at $P=10$, uniform TGV gives an NMSE of less than 3.5% and the SSIM exceeds 0.95.

This is also true qualitatively, as shown in the representative example of Fig. 3 ($P=10$). The top row shows results for uniform subsampling. ISTA grossly reconstructs the signal well at its lowest-frequency (in a spatial sense) regions, but performs poorly near rapid changes between low- and high-amplitude positions with varying phases.

Physiologically, the bright region on the left-hand-side is at the junction between OHCs and Deiters cells sometimes referred to as the hotspot [14]. It is one of the regions of most interest in these motion maps, so poor performance here is not acceptable.

The TV reconstruction, which uses a piecewise-constant prior, is not as smooth as the densely sampled signal. The TGV signal, on the other hand, performs well at all positions in the OCC. The NMSE between this reconstruction and the densely sampled signal is below 1%.

As for the random sampling pattern, all reconstructions fail to resemble the densely sampled signal well. Random sampling patterns may include large breadths of columns that are not sampled, which we refer to as bands. Bands are results of the structured nature of OCT subsampling, where samples are removed column-by-column rather than pixel-by-pixel.

The reconstruction performs most poorly at these bands. ISTA reduces the signal value to zero at the centers of these bands, while TV performs nearly constant interpolation across them. TGV still outperforms ISTA and TV, but still shows significant artifacts at the largest bands. This further motivates the choice of TGV with uniform subsampling as the method of further evaluation.

 

Fig. 2. Normalized mean square error and magnitude SSIM for the three tested methods: TGV, TV and ISTA. Values displayed are sample means over a test set ($N=20$), and error bars represent one sample standard deviation from the mean. Results are compared for both uniform and random subsampling by a factor of $P=2$, 5 and 10. The inset in the NMSE panels shows the results for TGV on a smaller set of axes, as the error is far lower than those achieved using TV and ISTA. TGV with uniform subsampling is seen to significantly outperform all other methods – using this method at $P=10$, the mean NMSE is less than 2%.

Download Full Size | PDF

6.2 Uniform TGV performance

We have shown that on the test set, TGV with uniform subsampling outperforms TV and ISTA, as well as any method with random subsampling. We now evaluate this method on all 275 maps within our data set.

Figure 4 shows NMSE and SSIM values for uniform TGV with three subsampling rates. Average values and standard deviations are computed and presented according to animal. For all animals, at one standard deviation above the mean, NMSE does not exceed 4.1% and SSIM exceeds 0.9 at even $P=10$. At $P=5$, NMSE does not exceed 1%. This shows quantitatively that uniform TGV performs well with generality in frequency, SPL and orientation.

To show that this is qualitatively true as well, we show one map from each animal in Fig. 5 with a subsampling rate of $P=10$. The example maps are chosen to be at different amplitudes and frequencies so that the generality of the method can be assessed. For all three example maps, TGV performs qualitatively well. Reconstructions are slightly smoother than the original signal, but still maintain the same phase shifts and hotspot behavior as in the original map. Some amount of smoothing may be more physical than the densely sampled motion maps measured by OCT, which suffer from issues such as shot noise, signal competition and blurring by an anisotropic point spread function [1,35]. As such, the NMSE is an over-estimate for the error between our signal and the ground truth, as the signal to which we are referencing is inherently noisy. Further interpretation of CSVi as a denoiser is discussed in Sec 7.1.

7. Discussion and conclusions

We have developed a theoretical formulation for compressed sensing for OCT-measured displacement measurements. We have tested six possible methods for CSVi on in vivo cochlear mechanics motion maps, and have found that TGV with uniform subsampling is capable of reconstructing densely sampled maps from only 10% of the samples with less than 5% error. Our use of uniform A-scan subsampling is in contrast to Lebed et al. [8], which uses random subsampling of A-scans in their OCT volume sparse recovery. We note that, at high subsampling ratios, Lebed et al.’s reconstructions contain similar banding artifacts to those shown in Fig. 3 H, I, J. Our results suggest that uniform subsampling may also improve 3D volume recovery.

 

Fig. 3. Example reconstructions using 10% of the M-Scans from the gerbil OCC. Colormap is shown in A, with hue representing phase re EC and saturation representing gain normalized to the maximum (further described in Sec 5.4). Sample: $\theta = 2$, 80 dB SPL, 10 kHz component. Maps are 100 rows by 200 columns, or 270 ${\mu}$m $\times$ 300 ${\mu}$m. Data are further described in Sec 5.1. Top Row: Results for uniform subsampling. A – densely sampled motion map; B – map from A subsampled uniformly by a factor of 10; C – dense map reconstructed using ISTA (2.38% NMSE); D – dense map reconstructed using TV (5.84% NMSE); E – dense map reconstructed using TGV (0.95% NMSE). Bottom Row: Results for random subsampling. F – densely sampled motion map(identical to A); G – map from F subsampled randomly by a factor of 10; H – dense map reconstructed using ISTA (45.55% NMSE); I – dense map reconstructed using TV (35.59% NMSE); J – dense map reconstructed using TGV (10.79% NMSE). Both qualitatively and quantitatively, TGV with uniform subsampling is seen to outperform the other methods on this sample.

Download Full Size | PDF

 

Fig. 4. NMSE and magnitude SSIM for TGV with uniform subsampling by factors of $P=2$, 5 and 10 across the full dataset ($N=275$), organized by orientation. At $\theta =1$, $N=75$; at $\theta = 2$, $N = 100$; at $\theta = 3$, $N = 100$. Data are further described in Sec 5.1. Values displayed are sample means over the set at each orientation, and error bars represent one sample standard deviation from the mean. Inset shows the NMSE for subsampling by a factor of 2, as it is much smaller than that achieved when subsampling by factors of 5 or 10.

Download Full Size | PDF

 

Fig. 5. Representative sample of three reconstructions made using TGV with uniform subsampling. Colormap is shown in B, with hue representing phase re EC and saturation representing gain normalized by maxima. A – Cartoon of approximate anatomy at $\theta = 1$ with basilar membrane (BM), Deiters cells (DC) and outer hair cells (OHCs) labeled. B – Densely sampled motion map for $\theta =1$, 60 dB SPL, 15 kHz. Map is 100 rows by 200 columns, and 270 ${\mu}$m $\times$ 300 ${\mu}$m. C – Map in B uniformly subsampled by a factor of 10. D – Dense map reconstructed from the subsampled map in B. E-H – Same as A-D, but with $\theta = 2$, 80 dB SPL, 9 kHz. Maps are 100 rows by 200 columns, and 270 ${\mu}$m $\times$ 300 ${\mu}$m. I-L – Same as A-D, but with $\theta =3$, 70 dB SPL, 25 kHz. Maps are 100 rows by 300 columns, and 270 ${\mu}$m $\times$ 330 ${\mu}$m. All data are further described in Sec 5.1.

Download Full Size | PDF

The method is also shown to function well at various beam axis orientations, stimulus frequencies and stimulus amplitudes, showing that it is sufficiently general to be of use for any experimental setup. This is the first example of compressed sensing being applied to the phase of the OCT signal.

7.1 Denoising interpretation of CSVi

The objective functions chosen do not force the recovered signal to be equal to the observation at any point, instead rewarding nearness to the observation in an $\ell _2$ sense. TGV will also discourage high-frequency phenomena, as they will incur large second derivatives. Our method can be interpreted as simultaneously performing compressed sensing and denoising.

Reduction of phase noise in OCT is challenging in vivo. In theory, one could reduce noise by averaging recordings over a longer period of time, but this incurs two issues in practice: 1) acquisition time is a limiting factor in in vivo experiments, and 2) sample drift worsens over the acquisition period [2,3]. Other noise sources, such as signal competition, are not time-dependent and also cannot be solved by averaging [35]. CSVi not only aids in reduction of experiment time, but may aid in the reduction of noise as well. Although it is difficult to quantify noise when our “ground truth" displacement maps are inherently noisy, this denoising feature can be seen qualitatively in the more smooth (and thereby more likely physical) features apparent in recovered signals of Fig. 5.

Many important features, such as the important “hotspot" delineation (the difference between the higher-intensity motion at the OHC and the lower-intensity motion at the supporting cells and BM) appear high in spatial frequency [14,21]. The sharpness of this delineation is unlikely to be physical, as the OHCs are not mechanically isolated from the rest of the OCC (Fig. 1). It is more likely that this apparent sharpness is a result of noise, and is exaggerated by insufficient sampling. TGV optimization penalizes these aphysical high-frequency components.

It is useful to compare our uniform sampling method to the classical sampling theory of Nyquist. Within this framework, a bandlimited signal can be faithfully represented by fewer samples so long as the sampling frequency is greater than twice the bandwidth. High-frequency features will lead to aliasing when upsampling and interpolating discrete-space signals, and as high-frequency noise is unavoidable some noise aliasing will always be present. No denoising will take place and phantom edge phenomena will occur.

7.2 Experimental applications

CSVi is most useful in that it enables various experiments that would be otherwise intractable. OCT experiments for cochlear mechanics broadly fall into two categories – 1) comparison of the motion in healthy cochlea at different anatomical positions within the OCC, and 2) perturbation studies in which the OCC is modified by either a genetic mutation or administration of a drug (e.g. furosemide or salicylate) [11,13,14,18,21,36,37]. The former have been enhanced in recent years by the development of methods to probe 2-D or 3-D motion both quantitatively [17,19,20] and qualitatively [14,22]. All such methods require achieving measurement of a sample from multiple angles, which is extremely time-consuming.

This has made it challenging, if not impossible, perform recovery/perturbation studies while probing either 2-D or 3-D motion. For example, recovery from furosemide injection on timescales of 10 minutes are of interest [13], but dense reconstruction of 2-D motion currently takes at least 40 minutes (20 minutes at each measurement angle) [19]. A factor of 10 reduction in sample acquisition, as offered by the presented CSVi method, would accelerate experiments sufficiently for 2-D effects of such perturbations to be probed.

Even in healthy cochlea, densely sampled 3-D motion maps have not yet been achieved in vivo. Even methods which could theoretically reconstruct such measurements are far too time-limited. Any 3-D reconstruction would require at least three measurements from distinct angles, and the noise level is strongly multiplied by the fact that large angular differences are challenging to achieve. For this reason, even 2-D reconstructions have been limited to responses to 80 dB SPL stimuli for which the SNR is sufficiently high. Only over-constrained reconstructions involving measurements at many angles could achieve 2-D or 3-D reconstructions at lower SPL.

Densely sampled in vivo 3-D motion at many SPLs would offer the most complete description possible of the healthy cochlea. If acceleration by a factor of 10 is achieved, 20-minute motion measurements at each angle can be reduced to a 2-minute acquisition time.