Optical coherence elastography (OCE) provides deformation or the material property mapping of soft tissue.1,2 The addition of elastographic contrast may improve the inherent ability of optical coherence tomography (OCT) to differentiate the composition and structure of soft tissue.3 Moreover, the mechanical information extracted from OCE is important for analysis and identification of pathological changes in soft tissue. For example, OCE may provide high-resolution characterization of strains in arterial walls, which would be important complementary information for determining the stability of atherosclerotic lesions.4 There are two main categories of OCE techniques, phase-based methods5,6 and speckle tracking techniques,78.–9 which rely on the structure of the speckle pattern when it is fixed. In general, speckle-tracking based OCE can measure greater deformation than phase-based OCE methods as phase-based OCE is limited by the phase stability of OCT system10 and the phase wrapping induced by large physical deformations or high detected particle velocities within the imaging volume. Although phase unwrapping could extend the measurement range of the phase-based method, it is difficult to apply due to noise corruption or discontinuity of the wrapped phase maps in OCT imaging.6,11 The principle of speckle tracking techniques has been previously described by Schmitt.1 Briefly, the speckle can be temporally tracked by quantifying the displacement via cross correlation of the OCT images of prestressed and stressed tissue samples. However, the resolution for the displacement calculation was limited to 1 pixel and no strain elastograms were given, as the process of strain calculation by differentiating the displacements was very sensitive to noise. Kirkpatrick et al.12 demonstrated that a maximum likelihood speckle shift estimator is superior than cross correlation, when the tissue motion between frames is less than 0.8 pixels. However, in practice, it is difficult to estimate the pixel shift a priori. Moreover, if the deformation values have a wide range from subpixels to pixels, the maximum likelihood will not be effective. Another drawback of the existing speckle tracking methods is the use of numerical differentiation of displacements to obtain strains. This procedure is noise sensitive as any error in the displacement measurement will be amplified in its strain calculation.13 Due to these complications, the large majority of the present speckle-tracking based OCE techniques has not been verified for their measurement accuracy. Therefore, more advanced algorithms are required to improve the measurement resolution and accuracy. We aim to develop a robust speckle-tracking based OCE methodology with subpixel resolution and improved strain measurement accuracy.
Digital image correlation (DIC)1415.–16 is widely accepted for mechanical testing. The basic principle of DIC includes the tracking of the same points (or pixels) between the two images recorded before and after deformation. When DIC is applied to continuum mechanics, strains can be simultaneously obtained with displacements.17 Simultaneous computation of displacement and strain differentiates DIC from regular speckle tracking algorithms employed in image processing when compared to pattern recognition algorithms. DIC has been used to measure the material properties of biological tissues and biomaterials.18 This was accomplished by employing artificial speckle patterns on the surface of the specimen. DIC can be implemented if the speckle pattern deforms together with the specimen surface as a carrier of deformation. Under the condition of small perturbations due to time or mechanical loading where the changes of the speckle pattern between frames are small and free of discontinuities, we hypothesize that the DIC algorithm can be applied to OCT images, resulting in optical coherence elastograms. Such a hypothesis can be tested via static loading using low-frame rating imaging systems and when in practice, it will require high-frame rate devices to address dynamic loading conditions. Furthermore, DIC requires high-spatial-frequency information of the speckle to optimize the cross-correlation calculation. Therefore, image contrast and nonuniform distribution of speckle may affect the resolution and accuracy of the results.19,20 The contrast and brightness of OCT images decrease with imaging depth, which may also result in a decreasing correlation coefficient as a function of depth. The prerequisite for the effectiveness of this DIC-based OCE is that the speckle patterns are correlated. The correlation stability of the OCT speckle images free of load is evaluated. The focus of this article is to present this novel DIC-based OCE technique for deformation imaging of soft tissue. The technique was evaluated in terms of system error and differentiation of various biological composite materials. Advantages, limitations, and further development are also discussed.
Materials and Methods
DIC is fundamentally a cross-correlation or speckle tracking technique. As such it requires a reference image of the object before deformation and an image after deformation. At each pixel point, the DIC processing chooses a subset of pixels centered at the point and searches for a subset with a maximum correlation coefficient on the deformed image. Applying the theory of deformation of continuum mechanics, the differences in the positions of the reference subset center from the correlated subset center on the deformed image yield six deformation parameters including displacements and the displacement gradients. Theoretical derivation of DIC for simultaneous deformation and strain measurement can be found in Ref. 14. Here, we used a normalized cross-correlation criterion, which is insensitive to overall frame intensity fluctuations. The normalized cross-correlation equation is expressed as1) is usually converted to an optimization problem where the minimization of 1- provides the solution to the six unknowns. The six deformation parameters are independent variables of function , which avoids numerical derivation of displacements for strain calculations.
The maximum correlation searching process includes two steps. First, a simple search scheme is applied which yields integer pixel resolution. Second, a subpixel displacement registration is applied. In this article, bilinear interpolation is applied to obtain gray level information between pixels, and the Newton–Raphson algorithm is applied for high-accuracy subpixel registration.17,21,22 The results of the first step serve as initial approximation for the Newton–Raphson iteration, which continues in the subpixel domain about that integer pixel to obtain precision of 0.01 pixels or higher. When convergence is achieved, the six parameters are solved simultaneously. We set the downward displacement in the axial direction as positive in our software algorithm.
OCE Procedure and System Calibration
The OCE system consists of a custom polygon-based swept-source OCT (SSOCT) system and a loading structure. A diagram of the OCT system is shown in Fig. 1(a). The OCT system has a 36-kHz A-scan frequency, axial resolution in air, and lateral resolution determined by the focal spot size of the imaging lens (LSM03, Thorlabs, Newton, New Jersey). A photo of the loading structure under the scanning lens is shown in Fig. 1(b).
The loading structure shown in Fig. 1(b) is made of a plate fixed to a linear translation stage. A 20-mm-diameter hole was fabricated on the plate, and a 1.2-mm-thick cover slip was glued to the bottom of the plate to allow the transmission of light and application of the load. Static load was applied by adjusting the translation stage downwards. DIC calculations were performed on two consecutive frames of OCT images of a sample before and after deformation. A schematic diagram of OCE procedure is shown in Fig. 2. The sample was cut into pieces and covered with solution to match the refractive index if necessary. Preload was applied while the first OCT image was taken.8 Preloading was applied to ensure that the glass window was in contact with tissue and decrease the correlation noise.23 It was applied by moving the compression plate down a certain distance, usually in the range of a few microns depending on the shape of the sample, beyond the point of first contact. The two frames of OCT images taken under preload and compressive load were then computed by DIC to obtain displacement distribution maps and elastograms of normal strains and shear strain simultaneously.
We first evaluated the system error by processing two consecutive images of a phantom without applying a load. The phantom was thick and made of RTV silicone (ELASTOSIL RT 601 A/B, WACKER, Germany) with as an optical scatter. The concentration of was , and the mixing ratio of A/B was 15:1, resulting in elastic modulus of .24 The two frames of OCT image of the phantom were taken at an interval of . One of the images is displayed in Fig. 3(a). The displacements and strains in this sample were expected to zero as no load was applied. Therefore, the resultant OCE images indicate the inherent measurement error of the system. Figure 3(b) demonstrates the correlation coefficient in a region of interest (ROI). Correlation coefficients were observed throughout the ROI indicating a high correlation in the speckle pattern. The displacement in the axial and lateral direction is shown in Figs. 3(c) and 3(d), respectively, from which it can be seen that the system error for displacement measurement was pixels in axial direction and pixels in lateral direction. The measurement error is slightly greater in the lateral direction, possibly due to jitter in the Galvo scanning mirror, subject to feedback mechanisms in its control loop. The normal strains in the axial direction and lateral direction are shown in Figs. 3(e) and 3(f), which demonstrate a measurement error of . We also evaluated the correlation changes of the speckle patterns with time. The phantom free of load was imaged with various time intervals of 0 to 15 min. The average and standard deviation of the correlation coefficients within the ROI was plotted with time in Fig. 3(g). The error bars demonstrated that the variations of the coefficient within the ROI are not significant. Although the correlation coefficient varies with time in general, it seems to be of a random nature. This random change could be explained by the random characteristics of speckles. The speckle pattern is sensitive to the system unstableness, such as minor vibrations in the optical table or jitter in the Galvo scanning mirror. The lowest correlation coefficient was 0.75 within the 15-min timeframe. A was defined as strong correlation as per previous publication.25 Thus, results obtained when were considered to be poor data as decorrelation between the speckle patterns existed.
Due to light absorption of tissue, the intensity of the image decreases with depth. In this experiment, the total imaged depth was representing 300 pixels in the axial direction. The changes of the root mean square (RMS) contrast26 and correlation coefficient with depth were evaluated and plotted in Figs. 3(h) and 3(i). As estimated, the correlation becomes weaker with the decreasing of image contrast.
To calibrate the relationship between physical deformation and pixel shift in the OCE image domain, a homogeneous phantom made of candle gel and , with a thickness of , was tested. The resultant relationship was termed as the ratio. A series of compressions ranging from 1 to 150 μm was applied to deform the phantom. The compressions were applied by adjusting the translation stage. The measurement results are shown in Fig. 4. One OCT image is shown in Fig. 4(a). The axial displacements of the sample undergoing 10-μm compression are overlaid on a structural image and shown in Fig. 4(b). The width of the ROI was , assuming that the refractive index of the sample was 1.4. The ROI includes 800 and 240 pixels in the horizontal and vertical directions, respectively. The 300 pixels evenly distributed within the ROI with 8-pixel interval were calculated by Eq. (1). The values of the pixels that are not directly calculated were obtained by interpolation. The displacement inflicted via the translation stage was 10 μm on the top surface of the phantom, which gradually decreased with depth. The measurement results in Fig. 4(b) do not show a gradual decreasing pattern in the axial direction as the variations of the displacement in the ROI were too small to be resolved by the DIC-based OCE algorithm. Figure 4(c) shows that the pixel shift in images linearly increases with the increasing of the compression. The horizontal coordinate of the “” is the average displacement applied in the ROI, and the vertical coordinate is the average value in the ROI measured by OCE. For example, at 150-μm compression, the displacement distribution of the phantom is from 150 μm on the top surface linearly dropping to 0 at the bottom of the phantom. Thus, the average value within the ROI can be calculated byFig. 4(c) for the calibration will improve the calibration accuracy. The ratio was calculated to be 5.8 based on the linear curve fitting. More specifically, one pixel shift between the two OCT images was caused by 5.8 μm deformation of the sample. The correlation coefficient decreases with the increasing of the displacement as shown in Fig. 4(d). The correlation coefficient was about 0.6 at 150 μm, which was considered to be the maximum measurable displacement. Maps of displacement in lateral and axial directions of the phantom under 150-μm compression are shown in Figs. 4(e) and 4(f). In the lateral direction, the scanned length was 5 mm including 1024 pixels and resulted in . Negative values demonstrate that the ROI is on the left side of the center of the phantom. The axial displacements in Fig. 4(f) show a decreasing trend with depth as expected.
Evaluation of OCE for Material Differentiation
A candle gel and phantom including a small stiff part made of silicon and were tested by the DIC-based OCE technique. To estimate the elastic modulus of the candle gel and the silicon inclusion, compression tests were conducted. Samples of thick were resected from the candle gel and silicon material, respectively. A load cell (FSH01045, FUTEK, Irvine, California) was fixed at the bottom of the loading structure as shown in Fig. 1(b). The 100-μm compression was applied to each sample by the translation stage while the load was recorded by the load cell. Then, the elastic modulus was estimated by assuming linear elastic material, where is the Young’s modulus, is the stress, and is the strain in the equation. The Young’s modulus () was experimentally derived to be 8.98 MPa for the small silicon inclusion and 2.99 MPa for the surrounding candle gel material. To make the candle gel phantom with stiff silicon inclusion, a few random sized pieces of the silicon phantom were cut off and embedded to the candle gel before it solidified. A photo of the phantom is shown in Fig. 5(a). A small piece enclosed by a rectangle on Fig. 5(a) was cut off for OCE experiment. The stiffer piece is vaguely visible in the OCT images of Figs. 5(b) and 5(c), which are taken before and after compression. The phantom was compressed by . Displacements and strains in the axial direction of an ROI obtained by the DIC-based OCE are shown in Figs. 5(d) and 5(e). It can be seen that the stiff block deforms less than the surrounding region. From the strain elastogram displayed in Fig. 5(e), the boundary of the stiff block is clearly visible. The strains generated on the block are significantly smaller than the strain on the surrounding material. A strain distribution map obtained by numerical derivation of the displacement is shown in Fig. 5(f), which is severely corrupted by noise. As it is difficult to estimate the strain distribution of the phantom which has such a nonregular shape inclusion, we conducted a finite element simulation using COMSOL4.3 to verify the OCE results. Isotropic linear elastic models were used in the simulation. Poisson’s ratio for both materials is 0.499 and Young’s modulus is 8.98 and 2.99 MPa, respectively, as per the experimental measurements. Under 20-μm compression, the normal strain distribution in the axial direction was displayed in Fig. 5(g). The pattern of the strain distribution matches well with the experimental results in Fig. 5(e). A large strain region above the inclusion is shown in green–blue in Fig. 5(e) and blue in Fig. 5(g). Discrepancy of the values between Figs. 5(g) and 5(e) is due to a few factors including: (1) shape: the geometry of the phantom in the simulation is not identical to the actual phantom, because the actual phantom, especially the stiff inclusion, does not have perfectly straight edges; (2) material of the phantom may not be homogeneous, as suggested by Fig. 5(d); (3) two-dimensional simulation was conducted while the actual phantom is a three-dimensional object; and (4) the material properties used in the simulation may differ from the actual values.
Biological Tissue Imaging
The DIC-based OCE technique was applied to image biological tissues. Unlike the phantoms, where the fundamental shape could be reasonably controlled, the biological tissue did not have regular shapes. Thus, the load applied may not be uniformly distributed, resulting in complicated stress analysis. Therefore, only the strain distributions of these samples through the DIC-based OCE were investigated without detailed conformational analysis of stress distributions.
A chicken breast sample was imaged by the DIC-based OCE algorithm. A small piece containing a layer of fat on the top, a soft membrane in the middle, and a muscle layer at the bottom was imaged. The sample was compressed by 10 μm. Imaging results are shown in Fig. 6. Figure 6(a) is an OCT image of the sample. The three layers are labeled as F for fat, Me for membrane, and Mu for muscle, based on estimation, as clear boundaries are not visible on the OCT image. Due to the nonregular shape of the sample, the cover slip was only in contact with the right part of the top surface. Figure 6(b) shows the map of displacement in the axial direction. The displacement distribution varies from 9 to 10 μm. Figure 6(c) is the elastogram of normal strain in axial direction, which shows that the softer membrane has higher strain of comparing with almost zero strain on the top and bottom layers.
A cerebral aneurysm is a blood-filled dilation (balloon-like bulge) in the wall of the cerebral artery. Rupture of the cerebral aneurysm can cause hemorrhage to the brain with potentially fatal consequences. It is essential to understand the mechanical conditions and consequent stress distribution in arterial walls that involved in aneurysmal rupture. Local variations in the anisotropy and inhomogeneity of arterial tissues have not been examined in detail due to limitations of imaging and experimental approaches. Walraevens et al.27 suggest that a compression test could be able to discriminate healthy from calcified aortic vascular wall tissue. In this work, a cerebral aneurysm sample using the DIC-based OCE was imaged as a feasibility study. The aneurysm sample was obtained under discarded human tissue exemption and had been fixed in formalin for several months before the OCE test. A photo of an aneurysm is shown in Fig. 7(a). The upper right corner of the sample was resected for the OCE test, which is shown in Fig. 7(b). The aiming beam indicates the imaged cross-section. One OCT image of the cross-section of the sample is shown in Fig. 7(c). Figure 7(d) is the histology of the cross-section correlating to Fig. 7(c). A dashed rectangular region is zoomed in Fig. 7(f). Figure 7(e) is the elastogram of normal strain in the axial direction overlaid on the OCT image. The ROI indicated by white arrows in Fig. 7(e) is correlated to the region indicated by black arrows in Fig. 7(d) where the adventitia and media have different structures. OCE image clearly shows the boundary of the adventitia layer, which is not clearly identified in OCT structural image in Fig. 7(c). The region indicated by a red asterisk is not processed because the speckle patterns are decorrelated. The sample was also imaged from the bottom side of Fig. 7(d). The elastogram of axial normal strains overlaid on an OCT image is shown in Fig. 7(g). A high strain region in blue indicated by a black arrow correlates to the local heterogeneity in the fibrotic aneurismal wall shown in the histological image in Fig. 7(f), as indicated by a black arrow. Although microcalcification is shown in Fig. 7(f) indicated by a red arrow, it is difficult to correlate it to the strain elastogram in Fig. 7(g). The small red regions in Fig. 7(g) pointed by a red arrow might be induced by noise.
While the error sources in the DIC measurement include OCT imaging noise, out-of-plane displacement, and error in the subpixel registration process, the most troublesome factor affecting OCE was speckle decorrelation. Speckle patterns may become decorrelated with time especially speckle images of a sample with a fluid component. The lowest correlation coefficient observed was 0.75 within the 15-min test shown in Fig. 3(g), which demonstrated that the decorrelation induced by random speckle changes was not an issue within a reasonably long (15-min) time interval. However, correlation stability varies with various stiffness of tissue under loading condition.28 In all experiments performed during this work, the correlation coefficient was in the ROI to ensure robust results. Subset size directly determines the area of the subset being used to track the displacements between the reference and target subsets, which are critical to the accuracy of the displacement measurement. The subset size should be chosen based on images to ensure sufficient distinctive intensity patterns are contained in the subsets.29 A -pixel subset was chosen for most of the experiments conducted in this article. The subset size determined the top surface of ROI, which has to be at least half of the width of the subset below the top of the image to perform correlation calculation. In this case, the upmost boundary of the ROIs is row 25 (in pixels). For phantom experiments shown in Figs. 3–5, the ROIs are below 25 pixels. For random-shaped samples in Figs. 6 and 7, a structural OCT image-based mask was used to create ROIs. From the system evaluation and calibration in Figs. 3 and 4, the resolution for axial displacement was measured to be 0.1 pixels, which correlated to . However, the calibration was dependent on the medium refractive index, as the axial resolution of the OCT image depends on the medium refractive index. Thus, if the refractive index is unknown, each sample should be calibrated to improve measurement accuracy. The resolution for lateral displacement measurement was about 0.2 pixels, equivalent to as demonstrated in Fig. 3(d). The accuracy for axial displacement measurement is , as demonstrated in Fig. 4(f). More homogeneous phantoms and well-designed experiments will be conducted to further evaluate the accuracy of the system for displacement and strain measurements in two dimensions.
The maximum OCE measurable deformation tested in Fig. 4 was before speckle decorrelation occurred. Factors for achieving large deformation measurements include (1) a larger subset of 69 pixels is used for deformation (for uniform distribution of speckle, larger subsets will provide higher accuracy); (2) the phantom material is relatively stiff; and (3) a phantom with thickness of was tested. At 150-μm deformation, only 1.5% of strain is generated. Similar amounts of compression applied to a thin sample will generate larger strains. The correlation between two speckle patterns of an object, before and after deformation, is more affected by the strain than the displacement. However, when imaging thick samples, the strain may be too small to be within the resolution of the DIC-based OCE. Under this situation, numerical derivation will be required to calculate strains. The phase-sensitive technique has the ability of measuring deformation in nanometer scale.6 It would be advantageous to small strain measurement.
Soft-tissue images rarely have uniform speckle distribution where speckle decorrelates under small compressions. The compression that was applied in the soft-tissue experiments was usually . In addition, as the tissue samples had irregular shapes, it was almost impossible to have ideal contact between the top surface of the sample and the compression plate. For example, in both Figs. 6 and 7, the glass plates were in contact with only a part of the top surface. In both experiments, speckle in the region under direct compression decorrelated, whereas on the other region, correlation remained. Therefore, OCE could be applied to the region without direct compression. This may provide some information for future experimental designs of loading schemes, although the stress analysis under this type of nonuniform compression may be quite complicated. In the case of Figs. 6 and 7, the regions in contact with the compression glass are either at the edge of the sample or too small [a point () contact]. Those regions were not in the interest of our analysis. Therefore, greater compressions (10 to 20 μm) were applied. However, it is completely possible to apply a small amount of compression and analyze the strain distribution in the region right under such compression. As a glass window has to be used to transmit light and apply load, the artifact due to reflection of the glass window, as shown in Fig. 6(c), can be destructive for image correlation. Such artifact should be eliminated by setting the imaging beam on optimal imaging angle of to 75 deg.
The experiment shown in Fig. 5 demonstrated the capability of the DIC-based OCE for identifying structural features and strains of various compositions of soft tissue. The three-layer structure could be differentiated in Fig. 6, even though a nonuniform compression was applied. However, to interpret the results, care should be given to the geometrical shape of the sample and the loading condition. For example, there is always strain concentration in the vicinity of structural discontinuities, where materials can not be differentiated simply based on strain maps. As the tissue samples and various composites included did not have regular shapes, interpretation of the strain elastograms must be careful. Strain elastogram in the axial direction alone is usually enough for the purpose of differentiation of various components of a sample. Normal strains in both axial and lateral directions and shear strain have all been obtained simultaneously through the DIC-based OCE. These data will be used for quantitative analysis of material properties in the future where constitutive models of the sample may have to be built and inverse problem solving be applied.30
In general, the speckle-based OCE methods are limited to measure large deformations ( for this DIC-based method) and affected by speckle decorrelation. However, the mechanical setup, control of the loading system, and data acquisition are generally simpler than phase-sensitive methods, which usually employ mechanical or acoustic waves to deform a tissue sample with potential advantage in terms of dynamic range.6,11,31,32 For speckle-tracking based methods, a loading plate with an optical window allowing the optical beam to pass is sufficient. Speckle/echo-tracking based ultrasound elastography has been successfully used in clinical applications by applying a slight pressure through the ultrasound probe.33 Similarly, this DIC-based speckle tracking algorithm has the potential to be integrated with a forward-viewing OCT probe34 to perform elastography.
OCE experimental results in Fig. 7 showing the displacements and strains of local components of the aneurysm wall had the potential to provide additional contrast for aneurysm characterization, despite the fact that formalin fixation could alter the aneurysm samples’ mechanical property. Although the stress condition of the samples was not completely analyzed, the strain elastograms were exciting. The stress failures of aneurysms will be studied next, which may provide crucial information for understanding the pathophysiology in the future.
The DIC-based OCE included interpolation and Newton–Raphson algorithms in the computation to solve 2-D displacements and strains simultaneously. This process avoided the numerical derivation of displacements for strain calculation. It greatly improved the measurement resolution and accuracy than normally used cross-correlation techniques in image processing. The DIC-based OCE showed promise for tissue deformation and strain measurement. The resolution of this DIC algorithm can reach 0.01 pixels, under ideal conditions. With our SSOCT system and subset size, the resolution for displacement measurement was (0.1 pixels) in the axial direction and 1 μm (0.2 pixels) in the lateral direction. The resolution for strain measurement was 0.5% in both the axial and lateral directions. The maximum measurable value was determined by the correlation status of speckle patterns. Displacements as high as 150 μm could be measured by the OCE system. More sophisticated data analysis algorithms will be developed by including the displacements and strains in two dimensions obtained from the DIC-based OCE. Information about the load applied will also be included to study the material properties, such as Young’s modulus, quantitatively in the future. Interesting chicken breast tissue layers are revealed by the OCE tests. Elastograms of aneurismal wall show strong contrast corresponding to features in histology including local heterogeneity and the layers of media and adventitia. More vascular samples will be measured and correlated with histology to verify the OCE findings. Benefiting from our technique’s inherent high resolution, the DIC-based OCE approach has a potential for characterization of atherosclerotic plaques and aneurysms,35 as well as other lesions.
This work was supported by Canadian Institutes of Health Research Natural Sciences and Engineering Research Council, and MITACS Elevate Postdoctoral Fellowship Program. We thank Androu Abdalmalak for helping with the COMSOL simulation.