Open Access
Add to BibSonomy
Abstract
A digital volume correlation (DVC)-based optical coherence elastography (OCE) method with inverse compositional Gauss-Newton (IC-GN) algorithm and second-order shape function is presented in this study. The systematic measurement errors of displacement and strain from our OCE method were less than 0.2 voxel and 4 × 10−4, respectively. Second-order shape function could better match complex deformation and decrease speckle rigidity-induced error. Compared to conventional methods, our OCE method could track a larger strain range up to 0.095 and reduce relative error by 30-50%. This OCE method has the potential to become an effective tool in characterising mechanical properties of biological tissue.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
The mechanical properties of tissues are related to their physiological and pathological states. For example, fibroadenoma and carcinomas breast tissue are four and eight times stiffer than normal tissue, respectively [1,2]. Various elastographic technologies have been developed to supply additional biomechanical contrast to morphological images for clinical diagnosis. Elastography based on ultrasound and magnetic resonance imaging (MRI) were proposed by Ophir in 1991 [3] and Muthupillai in 1995 [4], respectively. In 1998, Schmitt first reported the optical coherence elastography (OCE), an elastography based on optical coherence tomography (OCT) [5]. Benefit from OCT’s the ultra-high resolution, the OCE technology has been quickly developed in the last decade, and can be applied for precisely evaluating the microscale mechanical properties of tissue [6].
To date, two technique routes of OCE have been developed based on different principles, i.e. speckle tracking and phase-sensitive OCE [6]. Thereinto, the phase-sensitive OCE captures the interferometric phase from complex OCT signals. The sensitivity of the retrieved tissue displacement is in nano-scale, which is expedient to extract nonlinear mechanical parameters of tissue [7–9]. Nevertheless, phase-sensitive OCE measures the displacement along the OCT scanning direction, which is often locked into the phase wrapping when tracking large deformation. A robust phase unwrapping algorithm is required to unwrap the phase map [10]. The unwrapping algorithm is sensitive to the noises and artifacts in the raw OCT images, and the induced errors hinder the translation of phase-sensitive OCE to clinical practice.
On the other hand, Schmidt firstly proposed the two-dimensional (2D) speckle tracking OCE method in 1998 [5]. In this study, the resolution of displacement measurement was in pixel-level, and subsequently produced amplified error in the differential-based strain mapping. Considering anisotropic characteristics of biological tissue, 2D speckle tracking method is incapable of detecting out-plane displacement and thus three-dimensional (3D) full-field deformation measurement is essential in biomechanical research [11]. Fu et al. [12] firstly introduced digital volume correlation (DVC) technology to OCE. They reconstructed the 3D deformation map of silicone phantoms from OCT imaging by using a commercial DVC program (La Vision, Göttingen, Germany). The performance of their protocol reached the best size of sub-volume as $36\times 36\times 36$ voxels and the maximum applicable strain level as 0.0168. The DVC algorithm used in this study was based on fast Fourier transform (FFT), which may generate artificial fringes in deformation field and lead to confused interpretation of the calculation results. Meng et al. [13] put forward a static 3D OCE method by integrating inverse compositional Gauss-Newton (IC-GN) DVC method with swept source OCT system. This method enabled displacement measurement and strain calculation at a higher precision of 2 µm and 0.003, respectively. Besides, speckle tracking OCE method has been successfully applied to investigate the mechanical properties of porcine muscle, rat heart, porcine aorta and cornea [13–16]. However, the key challenge of further translating the speckle tracking based OCE technology into clinical practice is that the applicable strain range in conventional methods is too small to meet the actual measurement demand [17]. Therefore, a DVC algorithm that can measure a larger range of strain from OCT imaging is needed.
Since proposed, DVC algorithm has been optimized and improved in many aspects such as shape function, subset optimization algorithm and interpolation method [18]. For example, subset optimization algorithm evolved from gradient method, and Newton-Raphson (NR) method into Gauss-Newton(GN) method to improve the accuracy and speed of calculation [19]. Shape function was also developed from zero-order to second-order shape function. Zero-order shape function only describes the constant translational motions of sub-volumes. First-order shape function further takes the gradient components into account and meets the demand of describing simplified sub-volume deformation. However, biological tissues have complex micro-structure, which induces anisotropic and non-linear mechanical behavior [20]. The deformation in sub-volume scale is too complicated to be tracked using conventional first-order shape function. Besides, scattering particles in phantom, which are commonly used in OCT imaging as speckle information carrier for the correlation matching are stiffer than phantom materials. The particles are treated as non-deformable body and only rigid replacement(movement) are tracked to present the sample’s deformation. This hypothesis induces a calculation error named speckle rigidity induced error (SRI error). Theoretically, second-order shape function depicts high-order complex deformation in sub-volume could reduce the SRI error. Lan et al. [21] proposed an improved second-order shape function-based IC-GN DVC method and applied to confocal microscope imaging data, which achieved a local error less than 0.1 voxels. However, the application of second-order shape function-based IC-GN DVC method in OCE has not been reported and its performance has not been investigated yet.
In this study, numerical simulation of pure translation, linear displacement and non-uniform deformation were conducted with the IC-GN DVC method based on both first- and second-order shape functions. Phantom compression experiments were processed to determine the limit of strain range and assess the performance of using the DVC methods based on first- and second-order shape functions. In the following context, IC-GN DVC method based on first- and second-order shape functions are abbreviated as IC-GN$^{1}$ and IC-GN$^{2}$ for convenience.
2. Methods
2.1 DVC method
DVC is a 3D speckle tracking method, which is the extension of 2D digital image correlation (DIC). DVC tracks the point of interest (POI) between the reference volume and the deformed volume and calculates their relative displacement. Before computation, the whole volume is pre-divided into sub-volumes. The calculating process of DVC algorithm (Fig. 1) includes two main steps: integer voxel search and sub-voxel iteration. The coarse integer voxel search uses 3D fast Fourier transform-based cross-correlation (3D FFT-CC) to obtain the integer voxel displacement of sub-volumes in all directions. The displacement map is subsequently inputted to the Gauss-Newton sub-voxel iteration. In this study, the pipeline of DVC analysis was designed in MATLAB (R2022a, MathWorks, Inc., Natick, MA, US) environment, based on a published inverse compositional Gauss-Newton (IC-GN) algorithm [22], which could incorporate the shape and location of both reference and target sub-volumes simultaneously and thus avoided redundant computations.
Zero-mean normalized sum of squared difference (ZNSSD) criterion is used as correlation coefficient in IC-GN to evaluate the similarity between the reference and the deformed volumes because it is insensitive to the noise and the overall intensity fluctuation of sub-volume,
The relationship of voxels in the reference and the deformed sub-volumes are approximately described by shape functions. Considering the point $P(x, y, z)$ in the reference sub-volume and the point $P'(x', y', z')$ in the deformed sub-volume. First-order shape function can be described as:
However, linear deformation mapping in first-order shape function can only describe translation, rotation, scaling. Comparatively, second-order shape function accounts for the effect of second-order deformation parameters in sub-volume, therefore, could describe the complex deformation fields more accurately:
After obtaining the displacement components, strain can be subsequently calculated through point-wise least square fitting (LSF) approach. A local displacement field centered at an interrogated point is fitted using first-order polynomials:
The polynomial coefficients are determined by LSF, and each component in the strain tensor is calculated by the Cauchy’s formulas as:
A detailed procedure of LSF for calculating the strain tensor can be found in Supplement 1.
2.2 Simulated experiments
In order to eliminate potential errors from out-plane deformation, image acquisition quality and volume distortion during loading process, the precision of IC-GN$^{1}$ and IC-GN$^{2}$ were preliminarily investigated using simulated speckle patterns. The simulated speckle patterns were generated using an in-house built MATLAB code. In this algorithm, a set of Gaussian-shaped speckles were generated. The brightness and location of each speckle elements in the reference and deformed volumetric image were numerically defined as:
Diverse deformation fields were imposed to generate deformed volume. Firstly, a 20-frame series of volumetric images with sub-voxel translation were produced, ranging from 0 to 1 voxel with an incremental step of 0.05 voxel in x direction by Eq. (6). To simulate the imaging noise in reality, random Gaussian noise with the mean value of 0 and different standard variance of 0%, 3% and 5% of the full 8-bit grayscale were added to all of images. Then, each pair of the reference and the deformed volumetric images were computed by the aforesaid two DVC methods. In the calculation, the size of sub-volume was selected as $31\times 31\times 31$ voxels and grid step between the neighboring points was 4 voxels. Region of interest (ROI) was $500\times 500\times 216$ voxels which resulted in total 6.4$\times$10$^{5}$ calculated points. After calculation, the difference between the calculated displacement values and the prescribed displacement field was treated as computational errors, presented as two components, i.e., mean bias error and standard deviation (SD) error. The mean bias error of the displacement is defined as:
where $u_m$ is computed displacement of the mth calculated points and $u_{imp}$ denotes the actual imposed value of sub-pixel displacement; N is the total number of computation points. The SD error, which indicates the mean displacement that a set of numbers deviates from their mean variation of the measured displacement relative to the mean value, is defined as follows:In the second simulated experiment, a linear-changed displacement field in x direction was applied in deformed images. It was defined as $u=\varepsilon x$, where $\varepsilon$ was the strain value. In this study, a total of ten volumetric test images were produced with the strain values from 0.02 to 0.2 and an incremental step of 0.02 voxel. Subsequently the process of the first test was repeated here, i.e. adding noise and DVC calculation. Strain components were calculated from displacements by Eq. (5). Mean bias error and SD error of calculated strain were analyzed by Eqs. (7),(8).
To further quantitatively compare the spatial resolution and limitations of the two DVC methods, in the third part of simulated experiment, three superimposed nonuniform Gaussian-shaped deformation fields with different widths were given as Equation S7. In addition, the performance of IC-GN$^{1}$ and IC-GN$^{2}$ under different sub-volume sizes and noise levels were investigated. Random Gaussian white noise with different standard variance of 0%, 3% and 5% of the full 8-bit grayscale was added to the reference and deformed volumetric image. The root-mean-square-error (RMSE) is defined to represent the total error:
2.3 Mechanical experiment
2.3.1 OCT imaging system
The OCE experimental set-up in this study consisted of a Spectral Domain OCT system (TEL220PS, Thorlabs Inc.), a mechanical clamp and a translation stage (Fig. 2). The OCT system uses a broadband source with central wavelength of 1300 nm and the bandwidth of 100 nm. The maximum scan rate of source is 76 kHz. The axial resolution of OCT imaging is 3.5 µm (in air). Further illustration of OCT imaging principle is included in Supplemental Document.
2.3.2 Experimental procedure of OCE
The experimental flow chart of OCE is shown in Fig. 3. The Polydimethylsiloxane(PDMS) phantom (The process of preparing PDMS phantom refers to Supplemental Document) was given compression loading, while the OCT recorded the grayscale volumetric images of the phantom before and after deformation. Before DVC calculation, the OCT volumetric images were pre-processed using Gaussian filter to reduce the noise level. Then the aforementioned DVC algorithm was employed to track the volumetric internal grayscale speckle pattern by means of searching sub-volume in ROI of reference and deformed images. After getting the full displacement field, strain components were mapped from displacement by LSF.
During the OCT imaging acquisition, a camera in-line mounted with OCT scanner simultaneously captured the top surface images of PDMS sample before and after compression. The captured images were subsequently used for validation by using an open-source 2D DIC algorithm Ncorr [23]. Theoretically, the captured displacement and strain on the parallel layer snipped from 3D OCT-DVC method should be equal to these captured from 2D photography-DIC method. As the 2D photography-DIC technique has been well-developed and widely adopted in mechanical experiment, the 2D DIC results were used as a baseline to assess the performance of IC-GN$^{1}$ and IC-GN$^{2}$ methods.
2.3.3 System error evaluation
To evaluate the system errors of OCE system, the phantom was firstly captured OCT imaging at its static status, two successive frames of volumetric images were captured at an interval of 5s. The phantom was applied a horizontally translation by 0.5 mm and two volumetric images of phantom were captured before and after translation. The cubic computation ROI in DVC was $512\times 512\times 216$ voxels in x, y, z directions with a calculation grid step of 5 voxels and a sub-volume size of $41\times 41\times 41$ voxels. After the displacement field was obtained from DVC calculation, the point-wise LSF method with a window of $11\times 11\times 11$ voxels were implemented to calculate the internal strain field. In these two non-loading conditions, the displacement for static test and strain for translation test were expected to be zero. The inherent measurement errors of OCE system can be characterized by the calculated displacement and strain field, referred as virtual translation error and virtual strain error. These two type of errors represented the uncertainties of the developed OCE system, which indicated the measurable least displacement and strain value.
2.3.4 Linear compression test
To quantify the upper bound of strain that our proposed OCE system can track effectively, a series of compression experiments were conducted. The phantom was compressed on the clamp by the strain range of 0~0.1, with an interval of 0.01. Before test, an initial pre-compression was applied to the phantom to eliminate the effect of residual strain and eccentric compression. Followed the aforementioned protocol in Fig. 3, the OCT volumetric images were acquired and the displacement field were tracked using DVC algorithm. When calculating the strain, the displacement filed was cropped to avoid the boundary effect of the LSF. The 2D strain map processed using 2D DIC algorithm from captured phantom surface images was treated as the ground truth of applied strain and used for evaluating the 3D DVC strain calculated by IC-GN$^{1}$ and IC-GN$^{2}$. The relative mean errors were defined as the strain difference (accumulation of all computation points) between DVC and DIC.
2.3.5 Non-uniform deformation experiment
To investigate the performance of IC-GN$^{1}$ and IC-GN$^{2}$ on nonuniform deformation field, the phantom prefabricated with a hole was uniaxially compressed. Before compression, an initial pre-compression was also applied and then the phantom was compressed of 2.2mm. Two 3D volumetric and 2D surface images were captured before and after compression. After experiment, the 3D volumetric images were computed by the two DVC method and 2D surface images were calculated by DIC as ground truth as previous experiments. Calculation points inside the hole of volumetric images were discarded in computation. After computation, the relative displacement errors (accumulation of all computation points) between DVC and DIC were obtained.
3. Results
3.1 Simulated experiment
The simulated imaging dataset was given translation, linear and nonlinear deformation successively. Figure 4(b&c) presents the mean bias error of the computed displacement for IC-GN$^{1}$ and IC-GN$^{2}$, compared to the preset translation. The curves of mean bias (y-axis) are all in the form of a sinusoidal tendency against to the fractional voxel position (x-axis). This phenomenon is attributed to positional intensity interpolation error from sub-voxel reconstruction. The mean bias of IC-GN$^{1}$ was lower than that of IC-GN$^{2}$ by 10%-20%. Benefit by the IC-GN method’s anti-noise feature, with noise level increasing, only a small increase of bias error was observed. The maximum bias error of the two methods was less than 0.005 voxel. Figure 4(d&e) shows the SD errors of two methods. The SD errors of each method increased slightly following the noise level increasing and were nearly constant against different pre-assigned sub-voxel translations. With the noise level increasing from 0% to 5%, the maximum SD errors of IC-GN$^{1}$ raised from 0.0013 to 0.0017 while that of IC-GN$^{2}$ raised from 0.0013 to 0.0035.
Fig. 4. Schematic diagram of translation numerical experiments. The computed mean bias and standard deviation (SD) errors calculated in a $31\times 31\times 31$ voxels sub-volume: (a) the simulated volumetric image was given the translational movement ranging from 0 to 1 voxel with an incremental step of 0.05 voxel; (b) Mean bias errors of IC-GN$^{1}$; (c) Mean bias errors of IC-GN$^{2}$; (d) SD errors of IC-GN$^{1}$; (e) SD errors of IC-GN$^{2}$.
The simulated linear tension test is shown as Fig. 5. The error bar plot in Fig. 5(b) presents the difference between the preset strain and computed strain using IC-GN$^{1}$ and IC-GN$^{2}$. The rhombus denotes the average strain errors as in Eq. (7) and the range of bar denotes the SD as in Eq. (8). With the preset strain increasing, the SD of errors between preset and computed strain using IC-GN$^{1}$ and IC-GN$^{2}$ firstly constrained within 0.002 until the preset strain reached to 0.14. Thereafter, the SD of errors in IC-GN$^{1}$ significantly transcended that in IC-GN$^{2}$. Similarly, the average strain errors in IC-GN$^{1}$ and IC-GN$^{2}$ were initially close to each other, and after the preset strain reached to 0.14, the IC-GN$^{1}$ presented larger average strain errors.
Fig. 5. (a) Schematic diagram of tension numerical experiments with the cube representing the simulated volumetric image. (b) The error bar plots between computed and preset strain with IC-GN$^{1}$ (blue) and IC-GN$^{2}$ (red). The rhombus denotes the average strain errors and bar denotes the standard deviations. (c) The 3D maps of computed normal strain in tension direction of IC-GN$^{1}$ and IC-GN$^{2}$ with prescribed strain of 0.02 and 0.2.
The 3D maps of computed strain of IC-GN$^{1}$ and IC-GN$^{2}$ with prescribed strain of 0.02 and 0.2 are presented in Fig. 5(c). When the prescribed strain equaled to 0.02, the metrological performance of IC-GN$^{1}$ outperformed IC-GN$^{2}$. It is noted that the maximum errors of IC-GN$^{2}$ was 0.002 while that of IC-GN$^{1}$ was only 0.001. However, with the prescribed strain as 0.2, a large bias region, where error exceeded the prescribed strain by 0.01 was observed in IC-GN$^{1}$. On the contrary, the standard deviations of relative error of IC-GN$^{2}$ were less than that of IC-GN$^{1}$ when tracking this large strain.
Three analytical Gaussian-shaped non-uniform deformation field with different widths were shown in Fig. 6(a). Deformation field with diverse spatial frequencies can effectively examine the spatial resolution of IC-GN$^{1}$ and IC-GN$^{2}$-based DVC algorithms. Figure 6(b) is the plot of calculated displacement from IC-GN$^{1}$ and IC-GN$^{2}$ along the dotted line in Fig. 6(a), comparing with the prescribed analytical deformation field as reference. In all the three deformation fields, IC-GN$^{2}$ presented closer results to the analytical reference than IC-GN$^{1}$. Another finding was that both IC-GN methods underestimated the amplitude of the prescribed displacement, and with the spatial frequency region getting smaller, the peak gap between the calculated and the analytical curves was getting larger. In the worst condition, IC-GN$^{1}$ obtained nearly 50% of error of the peak strain compared to the prescribed counterpart at the spatial frequency region of 16 voxels. This was caused by the intrinsic low-pass smoothing feature of cross-correlation operation as a sliding average window.
Fig. 6. (a) Analytical Gaussian displacement field adopted to compare spatial resolution of IC-GN$^{1}$ and IC-GN$^{2}$. The widths of the three Gaussian functions are 16, 22, 32, respectively. The displacement field was defined using Equation S7. (b) A line plot across the red dotted line in (a) compares the computed displacement of IC-GN$^{1}$ and IC-GN$^{2}$ with the prescribed analytical displacement field.
As shown in Fig. 6(c), the RMSE of IC-GN$^{1}$ raised rapidly and continuously with the increasing of sub-volume size, while that of IC-GN$^{2}$ firstly dropped and then increased when sub-volume size raised from 11 to 41. The performance of IC-GN$^{2}$ was better than that of IC-GN$^{1}$, and its leading-edge enlarged with increasing of sub-volume size. Due to the noise robustness characteristic of IC-GN method, no significant difference was observed between different noise levels.
3.2 OCE experiment
The system error of the designed OCE system was evaluated by applying static status and translational displacement field to the phantom. Figure 7 shows the box-plot of the error of virtual displacement (captured displacement at static status) and virtual normal strain (calculated strain when applying only translation) in three directions calculated by IC-GN$^{1}$ and IC-GN$^{2}$. The range of virtual displacement and strain error from IC-GN$^{2}$ was larger than that from IC-GN$^{1}$ in all three directions. Measurement error in axial direction (z axis) was slightly smaller than lateral directions (x&y axis), due to the smaller resolution in axial direction which was determined by OCT source bandwidth rather than focal length of objective lens in lateral direction. Correlation coefficient of all computed points were above 0.9 in the ROI. The absolute value of virtual displacement was less than 0.2 voxel and absolute value of virtual normal strain was less than $4\times 10^{-4}$.
Fig. 7. The OCE system errors evaluation: (a) Virtual displacement error in the static situation (no translation); and (b) virtual strain error with translation only (no strain).
When applying linear displacement onto PDMS phantom, the average correlation coefficient and the mean strain error of all computed points compared with DIC (treated as ground truth) were plotted in Fig. 8(b). With the increasing of compressed strain, the average correlation coefficient dropped and the mean strain error increased. In our DVC algorithm, the accepted minimum of correlation coefficient was 0.8, where the corresponding strain was 0.095. The maximum mean strain error of all applied strain value was less than 15%, demonstrating the robustness of the proposed DVC algorithm. In addition, IC-GN$^{1}$ outperformed IC-GN$^{2}$ when the applied strain was less than 0.57; nevertheless, while prescribed strain increased to 0.57, IC-GN$^{2}$ presented better performance than IC-GN$^{1}$, where the mean strain error can be reduced by 30%-50%.
Fig. 8. The compression experiment with linear displacement: (a) A 2D OCT imaging slice sample and the calculated displacement field (locality). (b) The average correlation coefficient and mean strain errors of all computed points relative to DIC are plotted. (c) The normal strain in tension direction of IC-GN$^{1}$ and IC-GN$^{2}$ with prescribed strain of 0.02 and 0.07.
The computed normal strain of IC-GN$^{1}$ and IC-GN$^{2}$ with applied strain of 0.02 and 0.07 are presented in Fig. 8(c). When a small prescribed strain ($\varepsilon = 0.02$) was applied, the maximum and mean strain errors of IC-GN$^{2}$ was 0.002 while that of IC-GN$^{1}$ was only 0.0015. For the prescribed strain equaled to 0.07, there was a corner region where the computed strain value was underestimated to 0.04 in IC-GN$^{1}$, which was a large gap to the applied strain 0.07. In contrast, the calculated strain value of IC-GN$^{2}$ was close to the applied value and its maximum absolute error and relative error value were only 0.005 and 8%, respectively.
The phantom with a hole was used to generate the nonlinear strain field under the compression loading (Fig. 9(a)). The displacement error with different sub-volume size between the DIC ground truth and IC-GN$^{1}$ (blue) and IC-GN$^{2}$ (red) were plotted in Fig. 9(b). The rhombus denotes the average displacement errors and the range of bar denotes the SD error. When the size of sub-volume was less than 21, there was no adequate speckle information for IC-GN method to converge in iteration. With the increasing of the sub-volume size, the average displacement error of IC-GN$^{2}$ oscillated and slowly increased while that of IC-GN$^{1}$ increased rapidly. Additionally, the displacement SD error of IC-GN$^{2}$ also raised tardily with the increase of sub-volume size. Nevertheless, displacement SD error of IC-GN$^{1}$ kept increasing.
Fig. 9. The compression experiment with no-linear strain: (a) A 2D OCT imaging slice sample and the calculated displacement field around the hole (locality). (b) The error bar plots between computed and preset displacement with IC-GN$^{1}$ (blue) and IC-GN$^{2}$ (red). The rhombus denotes the average displacement errors and bar denotes the standard deviations.
4. Discussion
In this study, a robust OCT imaging and DVC algorithm-based OCE method was investigated. This method integrated the algorithms of 3D-FFT coarse search, IC-GN fine iteration, second-order shape function and point-wise least square fitting. The performance of IC-GN$^{1}$ and IC-GN$^{2}$ were analytically compared based on both numerical simulated dataset and OCT-imaged PDMS phantom experiments.
4.1 Advantage of IC-GN algorithm
The IC-GN algorithm used in this study is approx. two times faster and approx. 60% more accurate than the traditional NR method. Optimized by Pan et al. [24], the calculation procedure of IC-GN became more simple benefiting from the Hessian matrix pre-calculation during iteration and the elimination of the reconstruction of intensity gradients at sub-voxel locations within the deformed volume image. Experimental results indicated that the optimized DVC method could further reduce the iteration number by nearly 60%, consequently enhanced the computational speed by approximately 45-50 times. Because IC-GN and LSF methods are both noise insensitive, the integrated OCE method therefore remained to be noise robust when processing real OCT images with noise. Combined with parallel computing acceleration, OCE method presented in this study may be suitable for real-time clinical diagnosis, which adds additional biomechanical contrast to the conventional morphological images.
4.2 IC-GN$^{1}$ versus IC-GN$^{2}$
Compared with first-order shape function, the second-order shape function bring into the additional second-order displacement gradient terms. Thus, it is considered to have a higher accuracy when dealing non-uniform deformation. However, redundant parameters in IC-GN$^{2}$ may induce over-matched error between the low-order real displacement field and high-order shape function in the scenario with simple displacements, such as rigid translation. In this study, this phenomenon has been noticed in both the simulated and the mechanical experiments when applying only translation.
Instead of using sample’s natural speckle, adulterated scatters particles are usually used as deformation carriers to maintain the correlated characters of sub-volume in OCE method. Non-deformable character of scatters tends to induce negligible SRI error in mechanical experiments. In the simulated experiment, first-order shape function could correctly track particles of sub-volume in the initial phase. But when strain exceeded 0.12, SRI errors increased significantly and the mean errors and SD errors of strain of IC-GN$^{1}$ became larger than that of IC-GN$^{2}$. Correspondingly, in the OCE experiment, another error source affecting the DVC accuracy was speckle decorrelation, which was due to the blinking/boiling effect of speckle with an increase in applied strain. On this occasion, IC-GN$^{2}$ was found to have better performance than IC-GN$^{1}$ when preset strain value exceeded 0.06, which was smaller than simulation. In addition, a high correlation coefficient (0.8) was maintained until the prescribed strain increased to around 0.09. In this study, the upper limit of strain in OCE method was found to be superior to the previous results [12–14], which showed the robust performance of our proposed OCE method.
4.3 Application potentials
In general, nonuniform deformation is more common in practice. First-order shape function can hardly match prominent nonuniform deformation field. On the contrary, second-order shape function can describe complex displacement field with strain gradient. In vivo soft tissues have complex non-linear deformation and its information carrier structure have disparate mechanical properties in real world. Therefore, IC-GN$^{2}$ will reveal its advantage in the application on the biomechanical analysis of biological tissues. In addition, OCT is particularly suitable for aqueous translucent tissue and its micron-level resolution and millimeter penetration depth are beneficial for imaging micro-structure of biological tissue and for consequent DVC calculation.
4.4 Limitations and future works
In the simulated experiments, the Gaussian speckle pattern only ideally simulated the location and brightness of each speckle. The simulated speckle pattern has not been able to mimic the real OCT imaging process, hence some of the characteristics of OCT imaging, such as evolution of speckles during deformation, and the sub-wavelength speckle-induced multiplicative noise could not be included in the simulation study yet. The simulation study focused on the speckle pattern generated from OCT imaging simulation program will be studied in the future research.
To improve the performance of OCE method, the Gauss filtering was applied to reduce the noise level of captured OCT images. While other imaging pre-processing methods, such as polarization diversity, spatial compounding and frequency compounding [25], had not been investigated in this study, which will be included in the following studies.
Practically, it is critical but difficult in DVC algorithm to determine the suitable size for sub-volume to obtain precise displacement field. Although small sub-volume size brings low RMSE in the simulated imaging data, excessively small sub-volume size is not recommended for IC-GN$^{2}$ because there is no adequate speckle information in real image to fit unknown parameters of shape function. Otherwise, distribution and size of speckles may not be as uniform as the idealized simulated images and the noise level of real OCT image could also be higher. Thus, the subset size should be precisely overall considered based on image resolution, level of noise, image quality and strain level. The development of self-adaptive sized sub-volume for DVC method [26,27] is necessary for an accurate OCE measurement.
In addition, before it can be translated to clinical application, the proposed OCE technology has to be further developed and validated on ex vivo biological tissue samples, which is planed as future study.
5. Conclusion
In this paper, an OCE method based on IC-GN$^{2}$ was presented to solve 3D deformation field of OCT imaging of PDMS phantom. Noise-robustness characteristic of this proposed method was analyzed numerically. Benefit from the inherent high resolution of DVC method, the resolutions of this OCE method for displacement measurement were < 0.1 voxel and < 0.2 voxel in axial and lateral direction, respectively. Corresponding resolutions for strain measurement were about $3\times 10^{-4}$ and $4\times 10^{-4}$ in axial and lateral direction, respectively. The DVC methods based on first and second-order shape functions were compared in this study. According to the compression experiment, the maximum measurable strain based on IC-GN$^{2}$ method reached as high as 0.095. And IC-GN$^{2}$ could reduce the SRI error caused by the non-deformation of scatters. IC-GN$^{2}$ had better performance than IC-GN$^{1}$ when strain exceeded 0.12 and 0.06 in theoretical simulation and in real OCE experiment, respectively. The IC-GN$^{2}$ was found to reduce relative error by 30%-50% compared with IC-GN$^{1}$. Under nonuniform deformation field, IC-GN$^{2}$ method performed its robustness in both simulated images and experimental dataset for tracking the displacement with large strain gradient. Overall, the 3D OCE method based on IC-GN$^{2}$ presented a robust performance in the calculation of displacement and strain field. It has great potential to be further developed and translated to the characterization of biomechanical properties of biological tissue.
Funding
National Natural Science Foundation of China (11972118, 12172089, 61821002); Australian Research Council (DP200101970, DP200103492); Centre for Biomedical Technologies, Queensland University of Technology (Early Career Researcher Grant).
Acknowledgments
As the receiver of Roland Bishop Award, Dr. Jiaqiu Wang would like to express his appreciation to the Bishop family for their generous support of Biomedical Engineering Research.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. S. Chen, M. Fatemi, and J. F. Greenleaf, “Quantifying elasticity and viscosity from measurement of shear wave speed dispersion,” J. Acoust. Soc. Am. 115(6), 2781–2785 (2004). [CrossRef]
2. J. F. Greenleaf, M. Fatemi, and M. Insana, “Selected methods for imaging elastic properties of biological tissues,” Annu. Rev. Biomed. Eng. 5(1), 57–78 (2003). [CrossRef]
3. J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastography: a quantitative method for imaging the elasticity of biological tissues,” Ultrason. imaging 13(2), 111–134 (1991). [CrossRef]
4. R. Muthupillai, D. Lomas, P. Rossman, J. F. Greenleaf, A. Manduca, and R. L. Ehman, “Magnetic resonance elastography by direct visualization of propagating acoustic strain waves,” Science 269(5232), 1854–1857 (1995). [CrossRef]
5. J. M. Schmitt, “Oct elastography: imaging microscopic deformation and strain of tissue,” Opt. Express 3(6), 199–211 (1998). [CrossRef]
6. B. F. Kennedy, K. M. Kennedy, and D. D. Sampson, “A review of optical coherence elastography: fundamentals, techniques and prospects,” IEEE J. Sel. Top. Quantum Electron. 20(2), 272–288 (2014). [CrossRef]
7. R. K. Wang, S. Kirkpatrick, and M. Hinds, “Phase-sensitive optical coherence elastography for mapping tissue microstrains in real time,” Appl. Phys. Lett. 90(16), 164105 (2007). [CrossRef]
8. J. Chen, Y. Hu, X. Lin, H. Liu, and C. Sun, “Optical coherence elastography of bilayer soft tissue based on harmonic surface wave spectroscopy,” Opt. Lasers Eng. 145, 106667 (2021). [CrossRef]
9. C. Yang, Z. Xiang, Z. Li, N. Nan, and X. Wang, “Optical coherence elastography to evaluate depth-resolved elasticity of tissue,” Opt. Express 30(6), 8709–8722 (2022). [CrossRef]
10. C. Sun, B. A. Standish, and V. X. Yang, “Optical coherence elastography: current status and future applications,” J. Biomed. Opt. 16(4), 043001 (2011). [CrossRef]
11. F. Meng, X. Zhang, J. Wang, C. Li, J. Chen, and C. Sun, “3d strain and elasticity measurement of layered biomaterials by optical coherence elastography based on digital volume correlation and virtual fields method,” Appl. Sci. 9(7), 1349 (2019). [CrossRef]
12. J. Fu, F. Pierron, and P. D. Ruiz, “Elastic stiffness characterization using three-dimensional full-field deformation obtained with optical coherence tomography and digital volume correlation,” J. Biomed. Opt. 18(12), 121512 (2013). [CrossRef]
13. F. Meng, C. Chen, S. Hui, J. Wang, Y. Feng, and C. Sun, “Three-dimensional static optical coherence elastography based on inverse compositional gauss-newton digital volume correlation,” J. Biophotonics 12(9), e201800422 (2019). [CrossRef]
14. A. Nahas, M. Bauer, S. Roux, and A. C. Boccara, “3d static elastography at the micrometer scale using full field oct,” Biomed. Opt. Express 4(10), 2138–2149 (2013). [CrossRef]
15. V. A. Acosta Santamaría, M. Flechas García, J. Molimard, and S. Avril, “Three-dimensional full-field strain measurements across a whole porcine aorta subjected to tensile loading using optical coherence tomography–digital volume correlation,” Front. Mech. Eng. 4, 3 (2018). [CrossRef]
16. J. Fu, M. Haghighi-Abayneh, F. Pierron, and P. Ruiz, “Assessment of corneal deformation using optical coherence tomography and digital volume correlation,” in Mechanics of Biological Systems and Materials, Volume 5, (Springer, 2013), pp. 155–160.
17. K. V. Larin and D. D. Sampson, “Optical coherence elastography–oct at work in tissue biomechanics,” Biomed. Opt. Express 8(2), 1172–1202 (2017). [CrossRef]
18. B. K. Bay, T. S. Smith, D. P. Fyhrie, and M. Saad, “Digital volume correlation: three-dimensional strain mapping using x-ray tomography,” Exp. Mech. 39(3), 217–226 (1999). [CrossRef]
19. B. Pan and B. Wang, “Some recent advances in digital volume correlation,” Opt. Lasers Eng. 135, 106189 (2020). [CrossRef]
20. Y. Gao, T. Cheng, Y. Su, X. Xu, Y. Zhang, and Q. Zhang, “High-efficiency and high-accuracy digital image correlation for three-dimensional measurement,” Opt. Lasers Eng. 65, 73–80 (2015). [CrossRef]
21. S. Lan, Y. Gao, X. Xu, Y. Su, Y. Liu, C. Bai, S. Wu, and Q. Zhang, “Error analysis of surface-distribution and non-deformation of fluorescent beads for the ic-gn2 dvc algorithm,” Opt. Lasers Eng. 140, 106541 (2021). [CrossRef]
22. C. Li and R. Shu, “Accurate and simple digital volume correlation using pre-interpolation,” Meas. Sci. Technol. 31(9), 095201 (2020). [CrossRef]
23. J. Blaber, B. Adair, and A. Antoniou, “Ncorr: open-source 2d digital image correlation matlab software,” Exp. Mech. 55(6), 1105–1122 (2015). [CrossRef]
24. B. Pan, B. Wang, D. Wu, and G. Lubineau, “An efficient and accurate 3d displacements tracking strategy for digital volume correlation,” Opt. Lasers Eng. 58, 126–135 (2014). [CrossRef]
25. J. M. Schmitt, S. Xiang, and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). [CrossRef]
26. M. Gates, J. Gonzalez, J. Lambros, and M. Heath, “Subset refinement for digital volume correlation: numerical and experimental applications,” Exp. Mech. 55(1), 245–259 (2015). [CrossRef]
27. B. Wang and B. Pan, “Self-adaptive digital volume correlation for unknown deformation fields,” Exp. Mech. 59(2), 149–162 (2019). [CrossRef]
