Full-Field OCT (FF-OCT) is able to image biological tissues in 3D with micrometer resolution. In this study we add elastographic contrast to the FF-OCT modality. By combining FF-OCT with elastography, we create a virtual palpation map at the micrometer scale. We present here a proof of concept on multi-layer phantoms and preliminary results on ex vivo biological samples such as porcine cornea, human breast tissues and rat heart. The 3D digital volume correlation that is used in connection with the 3D stack of images allows to access to the full 3D strain tensor and to reveal stiffness anisotropy.
© 2013 Optical Society of America
Organ structures, tissues and cells have distinctive intrinsic mechanical properties. Moreover, the mechanical properties of tissues and cells are related to their structure and function: changes in those properties can reflect healthy or pathological states . Adding a map of mechanical properties as a supplementary contrast mechanism to morphological images could help diagnosis. During the last decade intensive research has been done to combine elastographic contrast with several biomedical imaging modalities with valuable results, particularly in MRI  or in ultrasonic imaging  where the elasticity information contributes to determining the final diagnosis. In the field of OCT imaging this mapping of mechanical properties was introduced by Schmitt  in 1998. Since this pioneering work, a number of methods have been developed such as 2D, 3D, static or dynamic elastography [5–10]. Here, we would like to show that FF-OCT setups are well adapted to perform static elastography of tissue samples, with minimal implementation effort. To the best of our knowledge this is the first time that FF-OCT-elastography has been demonstrated. The main advantage of the technique over the various OCT approaches is the ability to perform 3D measurements (that are mandatory when a soft matter sample is deformed) with sub-cellular resolution.
More precisely in this study we present the theoretical background and the experimental setup leading to the static method we have developed in order to add an elastographic contrast to a high resolution FF-OCT images. By doing so, one should be able to recreate a virtual palpation map at the micrometer scale achieved by FF-OCT. In this method we register a volumetric image before and after mechanical solicitation of the sample. From these two sets of images we calculate the 3D strain maps inside the sample by using finite element digital volume correlation based algorithms.
2. Experimental setup
The FF-OCT system is based on a Linnik interferometer  (Fig. 1). The illumination that is provided by a halogen lamp (SCHOTT-KL 1500 compact) takes advantage of a Kohler illuminator to ensure homogeneous illumination of the sample. The spectrum of the halogen source is broad, and the light passes through two filters: a high pass filter at 610 nm and a low pass filter at 1000 nm in order to match the spectral response of the camera and minimize the energy incident on the sample. Taking into account the spectral response of the camera, the spectral bandwidth is approximately 150 nm, centered at 710 nm.
A pair of identical water-immersion microscope objectives (Olympus, 10X, NA = 0.3) are placed in the arms of the interferometer. A silicon mirror (reflectivity of about 20 %) is placed in the reference arm. Images are recorded with a CCD camera (Dalsa 1M60, 1024 1024 pixels, 50 Hz, 12 bits, full well capacity of 150,000). The interference signal is isolated from the background by a PZT driven modulation of the path difference.
As high numerical apertures are used, FF-OCT systems have usually a better lateral resolution than other OCT setups but the drawbacks is that using high numerical aperture makes the systems more sensitive to multiple scattering and to tissue induced aberrations. With this custom setup, a lateral resolution of 1.4 µm and an axial resolution of 1 µm is achieved with a sensitivity of 77 dB (for acquisition of a single image) and the typical penetration depth in biological tissues is about 200 µm.
By moving this very compact Linnik interferometer (15 × 15 × 5 cm3) vertically (along the sample arm axis), it is possible to record 3D images of the sample.
The second part of the setup is a custom sample holder developed at the company LLTech which commercializes FF-OCT systems (Fig. 2). The holder has a piston system that can induce precise displacements along the z-axis for accurate sample compression. Thus, mechanical stress can be applied to the sample in a controlled manner. The sample is placed between a 1 mm thick glass window (matched in the reference arm) and the piston. A small force sensor (in blue in Fig. 2, Honeywell FS micro force sensor) is included between the piston and the sample to ensure that the measurement is taken just after contact with a minimal preload. In order to remain within the region of linear elasticity for the tissues under examination, a mean strain of about 1 % is applied between imaged reference and loaded states.
In order to highlight the mechanical properties of soft tissues, a volumetric image of the sample before and after compression is captured using the above described sample holder. After compression, a waiting period of approximately 2 minutes is observed before taking the second FF-OCT volumetric image in order to allow relaxation of the sample (as can be observed on the force signal). Then, from these two sets of images, the 3D strain map is obtained by using the 3D digital volume correlation method described below.
Previous works have already discussed uncertainties of strain measurements in Optical Coherence Elastography (OCE)  and Digital Volume Correlation method (DVC) . A proper estimation of local uncertainty field is key to the quality of the DVC analysis, and a more thorough analysis will be reported in the future. Indeed, when using FF-OCT (or any other tomographic imaging instrument), one do not master the local contrast. Moreover, this displacement resolution is also spatially varying as very deep sections display a lower signal to noise ratio. Therefore, a full quantitative analysis has to rely on a pre-determination of the uncertainty, and noise evaluation for each specific set of images.
3.1. Estimation of the displacement map
In this paper, a global 3D DVC  is used to estimate the displacement field within the sample volume. Here, a brief overview of the algorithm used to perform this correlation is provided. This method is based on the hypothesis that the 3D texture of the sample (here the density of scatterers) is assumed to be passively advected by the deformation. This hypothesis allows one to relate the reference volumetric image, f(x), acquired before the application of any load, to the volumetric image after a chosen mechanical loading, g(x),
In order to estimate the displacement map u(x), DVC consists in minimizing the quadratic difference between the reference image, f(x) and the deformed image g(x) corrected with a known displacement map ũ. In practice this consists of minimizing the following functional Γ:
Under the assumption that β(x) is a Gaussian white noise, the above formulation is the most appropriate. Γ being strongly non-linear, its minimization requires that the trial displacement field is reasonably close to the actual solution. Close to its minimum a Newton-Raphson scheme is used, so that calling g̃(x) = g(x+ũ(x)), we estimate the incremental correction δu(x) to ũ(x) as the minimization of the linearized functional Γ0.
This linearization is a valid approximation for small incremental displacements. By iteratively correcting the displacement ũ(x) ← ũ(x) + δu(x) and image g̃ by the estimated displacement field, the algorithm is able to converge to a displacement field.
In the following, the displacement field ũ(x) is decomposed over a finite element functions basis, Ψk, namely piecewise tri-linear functions over 8-node cubic elements where a regular and uniform cubic mesh is used over the entire region of interest. The trial displacement field can thus be written as:
The resolution of this system enables calculation of the 3D local displacement field in the samples. It should be noted that the spatial extension of the finite element basis functions is an important parameter. Indeed, this parameter is comparable to the correlation window size in classical Digital Volume Correlation (DVC) .
In order to estimate the validity of the calculation, we define a quality parameter, R, as the square root of the functional Γ[ũ] at convergence, suitably normalized by the dynamic range of the reference image.
A low value of the latter quantity means that most of the image signal has been captured in the matching of both images. It has been shown that this specific formulation of DVC which ensures continuity of the displacement field offers better performance than classical DVC . Furthermore, this method has the advantage of being very flexible with regards to the choice of finite element shape functions and mesh geometry. With this method it is possible to estimate displacements of about 1/100 pixel . However, the spatial resolution of the displacement field is strongly dependent on the contrast and texture (via the pair correlation function) of the imaged tissue, on the image noise, but also on the chosen kinematic basis. In particular, it was shown in  that the mesh size, l, (size of the cubic elements expressed in voxel) plays a key role; A fine mesh (small l), allows for a good spatial resolution at the price of a high uncertainty on the measured displacements. On the contrary, a coarse mesh (large) will have a small uncertainty but will not be able to capture rapidly varying displacement fields. Hence a compromise has to be chosen depending on image quality and sought resolution. In all case, a small value of R is an essential feature to ensure that the analysis quantitatively accounts for the true kinematics.
3.2. From displacement map to stiffness
Most soft biological tissues can be considered as incompressible, and in particular the different test cases shown below. Moreover, only small strains are considered and hence we restrict our discussion to linear elastic incompressible media.
Incompressibility means that the displacement field is expected to obey everywhere in the tissue
There are two different ways of dealing with such a constraint:
- - Either DVC ignores this incompressibility constraint, and the result of the analysis will either prove this property, or alternatively div(u(x))2 can be used as an indicator of the quality of the displacement evaluation.
- - Or, one can restrict the minimization of Γ[ũ] to the sub-space of isochoric displacements.
The second route will be followed herein, through a penalization technique. Namely, in addition to Γ[ũ], a regularization functional is considered as the space integral of the quadratic norm of unbalanced stress with the bulk. This stress is estimated from linear elasticity of an incompressible solid, so that to dominant order it resumes to
A weight A is attributed to the regularization, and Γ[ũ] + AΓreg[ũ] is finally minimized. A is the penalty factor that ensures incompressibility (or uniform volume change is the boundary conditions reveal incompatible with an isochoric constraint). In addition, when the regularization functional is not built strictly from an incompressible solid (a Poisson ration of 0.49 was chosen in our case), it acts as a smoothing filter below a length scale which (after a suite normalization of the functionals) is equal to the fourth root of A. In the following, this length scale is chosen as twice the mesh size, l, so that its effect is quite modest apart from setting the incompressibility constraint.
Let us introduce the strain tensor, ε(x), as the symmetric part of the displacement gradient
Incompressibility implies that the strain tensor is traceless, tr(ε(x)) = 0. Within the framework of linear elasticity the strain magnitude can be captured through the equivalent von Mises strain which reads
Under the incompressibility assumption, a linear elastic law simply relates the stress tensor, σ(x), to the strain, ε(x), as
The full 3D components of the displacement field are to be properly evaluated in order to estimate relative variation of the shear modulus. Indeed as observed in the Lams equation, there is no way to determine absolute values of μ(x) if only a kinematic information is exploited.
4.1. Proof of principle on a phantom
The method was first applied to test optically homogeneous and mechanically heterogeneous samples. The sample was made of Polydimethylsiloxane (PDMS Sylgard 184) and Zinc oxide (ZnO – 205532 SIGMA-ALDRICH) particles. ZnO particles were selected to induce optical scattering and PDMS to control the stiffness of the sample because stiffness depends on the cross-linker concentration.
With a ZnO particle concentration of 0.2 mg/g and a stack of three 50 µm thick layers of differing stiffnesses manufactured using the spin-coating method: the first layer has 10 % of cross linkers corresponding to a Youngs modulus of about 2.6 MPa , the second layer has 3 % of cross linkers corresponding to a Youngs modulus of about 0.55 MPa , the third layer have 5 % of cross linkers corresponding to a Young modulus of about 1.03 MPa . The following results shown in Fig. 3 were obtained. Figure 3(a) shows a representation of the layered sample structure. A cross-section of the volumetric FF-OCT image is shown in Fig. 3(b). In Fig. 3(c) the cross-sectional image of the evaluated strain map corresponding to the FF-OCT image plotted. Figure 3(d) is the profile of the strain averaged over y.
With this 1D layered sample geometry, one can both clearly distinguish the three layers on the strain map (that cannot be distinguished on the FF-OCT section) and also have access to a relative value of the stiffness. The ratio between the strain in the second layer and the third layer is equal to 1.53 which corresponds well to the ratio of the Youngs moduli equal to 1.87. It is difficult to give a value for the first layer because the relative stiffness between this layer and the neighboring one is too high. Thus the first layer undergoes negligible strain that can hardly be measured.
In order to compute this strain map we used l=16 voxels. This size was chosen because 16 voxels is approximately the distance between two ZnO particles and there is not enough signal backscattered light between ZnO particles. With such an element size, the residual R of equation (2) is 0.1 % which is a remarkably low value.
4.2. Ex vivo biological tissues
The final objective of this study being to aid diagnosis by complementing the FF-OCT images, three preliminary results on biological samples were obtained. The first sample was fresh porcine cornea, the second was freshly excised human breast tissue and the last was fresh rat heart. The samples are cut using a standard surgeon scalpel. This procedure is enough to ensure a good contact with the windows with the minimum preload mentioned before in order to limit the artifacts induce by a non-uniform tissue loading.
4.2.1. Ex vivo porcine cornea
Human and animal cornea have previously been studied with FF-OCT . The current setup was tested on cornea because a strong mechanical contrast between the epithelium and the stroma is expected. The current study aimed at highlighting this heterogeneity.
Figure 4(a) is a cross-sectional slice of the volumetric FF-OCT image of an ex vivo porcine cornea. We can easily distinguish the epithelium with small round cells and underneath the stroma with undulated collagen sheets. Figure 4(b) shows a cross-sectional image of the evaluated strain map corresponding to the FF-OCT image of Fig. 4(a). The two layers clearly exhibit different stiffness. The strain of the epithelium layer is about two times larger than the deformation of the stroma. The element size was chosen to be l=8 voxels. The residual level R = 0.8 % is remarkably low in this example.
4.2.2. Human ex vivo breast tissue
Recent studies on fresh breast tissue have demonstrated good correspondence between histology slides and FF-OCT images . These studies show very good agreement between histopathological and FF-OCT based diagnoses by exploiting only the morphology of the tissue (sensitivity 93%, specificity 75%). In order to improve the sensitivity and specificity, the current study aims to combine FF-OCT and elasticity maps in the hope of mimicking the successful approach combining elastography with ultrasonic imaging  that was demonstrated in our laboratory. Here, this preliminary test looks at the mechanical contrast between an area of fatty tissue (adipocyte cells) surrounded by normal fibrous tissue. Figure 5(a) is a cross-sectional slice of the volumetric FF-OCT image of the sample. In this image, the adipose tissue (big black cells) is clearly distinguished from the surrounding fibrous tissue. Figure 5(b) is a cross-sectional image of the evaluated strain map corresponding to the FF-OCT image. Indeed, the far greater deformation of the fat cells in relation to the surrounding tissue can be appreciated. Due to the absence of signal from inside the adipocyte cells, an element size l=32 voxels was used, leading to a residual R of 1.1%.
4.2.3. Ex vivo rat heart
In this example the strain map was used to detect the muscle fiber direction in the ex vivo rat heart. Indeed, muscle fibers have anisotropic mechanical properties . The muscle fiber will not show the same strain along the fiber axis and perpendicular to the fiber axis under the same stress. Using this property of the muscle fibers and the fact that the FF-OCT elastography method provides access to a 3D strain map, it is possible to detect a change in the fiber orientation. For each pixel the eigenvector corresponding to the largest eigenvalue of the strain tensor is calculated. This vector ε⃗eigen gives the principal strain direction. The principal strain direction indirectly reflects the fiber direction.
In the heart, the muscle fibers are organized in a stack of layers of different orientations . In this example, the FF-OCT image was captured at the border between two layers. Figure 6(a) is a cross-sectional slice from the volumetric FF-OCT image of the sample. In the FF-OCT image it is difficult to evaluate the fiber direction, and the two different fiber orientations are not visible. However, if the angle of the projection of the vector ε⃗eigen on the (x, y) plane (Fig. 6(b)) is plotted, it is clear that there are two different fiber orientations.
This example shows that the local strain map provides access not only to stiffness information but also to other mechanical properties of the sample such as mechanical anisotropy or compressibility.
5. Discussion and conclusion
These preliminary results were obtained using an unmodified commercial setup. They show that it is possible to produce accurate 3D displacement maps with a residual below 1 % at the micrometer scale in biological tissues. From these displacement maps it is possible to calculate 3D local strain maps that reflect the local contrast in elastic properties. However, as the local stress map remains inaccessible, there is still a lack of quantitative information.
The stiffness information, although not quantitative, could aid medical diagnosis. The next step will be to validate that the strain information can detect or reveal stiff inclusions such as, for example, tumorous breast tissues that are known to be very stiff . Moreover, having access to a full 3D local strain tensor not only gives access to stiffness information but also open the way to other contrasts such as mechanical anisotropy as shown here on ex vivo rat heart example (Fig. 6) or compressibility.
In the current study, the stress applied to samples is taken into account only in order to have a idea of the preload applied to the samples. This measurement of the stress with a flat sensor is not well adapted to the current setup because if the sample is not flat the measurement is biased. In order to work with biological tissue, which is not flat and is mechanically non-linear , it is necessary to precisely control the stress being applied. In order to measure the applied stress, a custom force sensor is under development. With this new force sensor it should be possible to measure an accurate and local value of the applied force.
Furthermore, a FF-OCT front view endoscopic adaptation of this method is currently under development for in vivo and in situ imaging. The setup developed by A. Latrive et al.  is being used to capture the 3D FF-OCT image of the sample, while the sample is compressed by the tip of the endoscope. Preliminary results on ex vivo tissues are very encouraging .
The authors would like to thank Fabrice Harms, Anne Latrive, Franck Martins, Charles Brossollet, Kate Grieve and Eugenie Dalimier, members of LLTech, Paris (France), for technical support and advice, Karima Mouslim and Martine Antoine of Tenon Hospital, Paris (France) and Armelle Rancillac of the Neurobiology Laboratory, ESPCI-ParisTech, Paris (France), for providing the excised samples and advice. This research has been supported by the INSERM program DESP Physique Mathématiques Sciences ingénieurs et Cancer and also by the LABEX WIFI (Laboratory of Excellence within the French Program ”Investments for the Future”) under references ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*.
References and links
1. A. Sarvazyan, Shear acoustic properties of soft biological tissues in medical diagnostics, J. Acoust. Soc. Am. 93(4), 2329–2330 (1993). [CrossRef]
3. M. Tanter, J. Bercoff, A. Athanasiou, T. Deffieux, J.-L. Gennisson, G. Montaldo, M. Muller, A. Tardivon, and M. Fink, Quantitative assessment of breast lesion viscoelasticity: initial clinical results using supersonic shear imaging, Ultrasound Med. Biol. 34(9), 137386 (2008). [CrossRef]
4. J. Schmitt, OCT elastography: imaging microscopic deformation and strain of tissue, Opt. Express 3(6), 199211 (1998). [CrossRef]
5. C. Sun, B. Standish, and V. X. D. Yang, Optical coherence elastography: current status and future applications, J. Biomed. Optics 16(4), 043001 (2011). [CrossRef]
6. A. S. Khalil, R. C. Chan, A. H. Chau, B.E. Bouma, and M. R. K. Mofrad, Tissue elasticity estimation with optical coherence elastography: toward mechanical characterization of in vivo soft tissue, Ann. Biomed. Eng. 33(11), 1631–1639 (2005).
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. R. K. Wang, Z. Ma, and S. J. Kirkpatrick, Tissue Doppler optical coherence elastography for real time strain rate and strain mapping of soft tissue, Appl. Phys. Lett. 89(14), 144103 (2006). [CrossRef]
10. B. F. Kennedy, X. Liang, S. G. Adie, D. K. Gerstmann, B.C. Quirk, S.A. Boppart, and D.D. Sampson, In vivo three-dimensional optical coherence elastography, Opt. Express 19(7), 66236634 (2011). [CrossRef]
12. B. F. Kennedy, S. Howe Koh, R. A. McLaughlin, K. M. Kennedy, P. R. T. Munro, and D. D. Sampson, Strain estimation in phase-sensitive optical coherence elastography, Biomed. Opt. Express 3(8), 1865–1879 (2012). [CrossRef] [PubMed]
13. H. Leclerc, J.N. Pri, F. Hild, and S. Roux, Digital Volume Correlation: What are the limits to the spatial resolution?, Meca. Ind. 13, 361–371 (2013). [CrossRef]
14. S. Roux, F. Hild, P. Viot, and D. Bernard, Three dimensional image correlation from X-Ray computed tomography of solid foam, Compos. Part A-Appl. S. 39, 1253–1265 (2008). [CrossRef]
15. H Wang Zhixin, Polydimethylsiloxane Mechanical Properties Measured by Macroscopic Compression and NanoindentationTechniques(2011). Graduate School Theses and Dissertations.
16. K. Grieve, G. Moneron, A. Dubois, J.-F. Le Gargasson, and A. C. Boccara, Ultrahigh resolution ex vivo ocular imaging using ultrashort acquisition time en face optical coherence tomography, J. Opt. A-Pure Appl. Op. 7(8), 368373 (2005). [CrossRef]
17. O. Assayag, M. Antoine, B. Sigal-Zafrani, M. Riben, F. Harms, A. Burcheri, B. Le Conte de Poly, and A. C. Boccara, Large Field High Resolution Full Field Optical Coherence Tomography: A Pre-clinical study of human breast tissue and cancer assessment, Tech. Canc. Res. Treat. 1(1), 21–34 (2013).
18. J. Bercoff, M. Tanter, and M. Fink, Supersonic Shear Imaging : A New Technique for Soft Tissue Elasticity Mapping, IEEE Trans. Ultrason., Ferroelectr., Freq. Control 51(4), 396–409 (2004). [CrossRef]
19. J.-L. Gennisson, S. Catheline, S. Chaffai, and M. Fink, Transient elastography in anisotropic medium: Application to the measurement of slow and fast shear wave speeds in muscles, J. Acoust. Soc. Am. 114(1), 536–541 (2003). [CrossRef] [PubMed]
20. W.-N. Lee, M. Pernot, M. Couade, E. Messas, P. Bruneval, A. Bel, A.-A. Hagege, M. Fink, and M. Tanter, Mapping Myocardial Fiber Orientation Using Echocardiography-Based Shear Wave Imaging, IEEE Trans. Med. Imag. 31(3), 554–562 (2012). [CrossRef]
21. I. M. Ariel and J.B. Cleary, Breast Cancer Diagnosis and Treatment (New-York: McGraw-Hill, 1987).
22. P. Wellman, R. H. Howe, E. Dalton, and K. A. Kern, Breast tissue stiffness in compression is correlated to histological diagnosis, Harvard BioRobotics Laboratory Technical Report , (1999).
23. A. Latrive and A. C. Boccara, Full-field optical coherence tomography with a rigid endoscopic probe, Biomed. Opt. Express 2(10), 2897–2904 (2011). [CrossRef]
24. A. Latrive, A. Nahas, and A. C. Boccara, Multimodal endoscopic Full-Field OCT imaging and elasticity mapping with a needle-like probe, SPIE photonics WEST BIOS conference, paper 8571-25 of conference 8571 (2013).