3 D Strain and Elasticity Measurement of Layered Biomaterials by Optical Coherence Elastography based on Digital Volume Correlation and Virtual Fields Method

The three-dimensional (3D) mechanical property characterization of biological tissues is essential for physiological and pathological studies. A digital volume correlation (DVC) and virtual fields method (VFM) based 3D optical coherence elastography (OCE) method is developed to quantitatively measure the 3D full-field displacements, strains and elastic parameters of layered biomaterials assuming the isotropy and homogeneity of each layer. The integrated noise-insensitive DVC method can obtain the 3D strain tensor with an accuracy of 10%. Automatic segmentation of the layered materials is realized based on the full field strain and strain gradient. With the strain tensor as input, and in combination with the segmented geometry, the Young’s modulus and Poison’s ratio of each layer of a double-layered material and a pork specimen are obtained by the VFM. This study provides a powerful experimental method for the differentiation of various components of heterogeneous biomaterials, and for the measurement of biomechanics.


Introduction
Optical coherence elastography (OCE) uses optical coherence tomography (OCT) as an imaging modality to measure the mechanical properties of biological tissue [1,2].Thanks to the inherent high resolution and ~2 mm depth penetration of OCT, OCE has the potential for early disease diagnosis.To extract mechanical information such as strain and elastic modulus from OCT interference signals, phase-based [3,4] and speckle tracking based [5] methods can be used.Speckle tracking based OCE methods have an advantage over phase-based methods for simultaneous multi-dimensional static deformation measurement.
In 1998 Schmitt first proposed the concept of OCE or OCT elastography [5].He used a normalized cross-correlation (NCC) speckle tracking method to measure the displacements and strains of gelatin phantoms and tissue specimens including pork meat, and intact skin in two dimesnions (2D).This work demonstrated the feasibility of applying speckle correlation to extract mechanical properties from OCT images, which are fundamentally laser speckle images [6].In 2004, Rogowska et al. [7] improved the speckle tracking algorithm by using a zero-mean normalized cross-correlation (ZNCC) criterion.
They found that the noise in axial and lateral displacement maps could be reduced by increasing the size of the subset (cross correlation kernel).Chan et al. [8] and Khalil et al. [9] then proposed a variation framework to reconstruct the Young's modulus of arterial tissue from the speckle correlation displacement measurement, using the inverse problem solving.Rogowska et al. later applied the speckle correlation method to measure the deformation and elastic modulus of atherosclerotic arterial samples [10].However, only the average value of Young's modulus was obtained.Difficulties remain in measuring the Young's modulus of each layer of the aorta specimen.Sun et al. [11] conducted detailed experimental studies on the feasibility of combining a Newton-Raphson (NR) digital image correlation (DIC) algorithm [12] with OCT imaging for the strain measurement.They concluded that the speckle patterns were correlated well when the strain is less than 10%.Although measurement errors have not been well analysed and quantified from the above studies, speckle correlation methods have been proved to be applicable to integrate with OCT imaging for the deformation measurement.
Recently, with the development of digital volume correlation (DVC), 3D OCE techniques have been proposed.In 2013, Fu et al. [13] first combined OCT imaging with DVC for the 3D displacement and strain measurement by utilizing the commercial DVC software from LaVision (Göttingen, Germany).Later, they measured the full-field deformation of corneal considering the distortion induced by refraction in OCT reconstruction [14].They also proposed to use the virtual fields method (VFM) to calculate the Young's modulus and Poisson's ratio of the homogeneous phantoms of two shapes [13].The results were promising although there were fringe-like errors in the strain maps.In the same year, Nahas et al. [15] conducted proof-of-concept studies of combing full-field OCT with DVC for the strain tensor measurement of multi-layer phantoms and biological samples.The strains were calculated from the displacement gradient.This study demonstrated the feasibility of the method without a detailed evaluation of the accuracy of the results.Lately, Acosta Santamaría et al. measured the deformation of the porcine aorta immersed in tissue clearing agents under tensile loading by OCT and commercial DVC software [16].These studies demonstrated the feasibility of using a DVC algorithm to process OCT images.However, the data format is often limited by the commercial software.Girard et al. conducted a series of studies.They developed a DVC algorithm based on a differential evolution method for the 3D displacement measurement [17], then proposed a VFM method, which is insensitive to rigid body motion, to recover the constitutive parameters of the optic nerve head [18].This work provided a paradigm for combining OCT, DVC and VFM for the mechanical properties' measurement of heterogeneous tissue.However, using strains instead of displacements as the input for VFM calculation may be more accurate.Although there are not so many studies of combined OCT and DVC yet, it can be seen that it is completely feasible to combine them for the 3D displacement measurement.However, the noise insensitive method is required to calculate the strain tensors in dealing with the laser speckle noises in OCT images.In addition, when objects with complex structures are measured, a large amount of virtual displacements and large-scale linear equations maybe involved in VFM, which may result in unsolvable problems or wrong solutions.Thus, a method to simplify the VFM calculation is required when a nonhomogeneous sample is measured.
This study is to develop a DVC and VFM based OCE method for the quantitative characterization of mechanical properties of layered biological tissue.As the DVC algorithm has been advanced in all aspects since it was first proposed in 1999 [19], we picked out the steps that are insensitive to noise and integrated them into a 3D strain measurement method with OCT imaging.This included a zero-mean normalized sum of squared difference (ZNSSD) criterion [20], a 3D inverse compositional Gauss-Newton (IC-GN) matching method [21], and a local ternary quadratic polynomial fitting method [22] for strain calculation.We then chose a linear VFM method which is noise resistance [23].The automatic segmentation of the layers is proposed based on the strain gradient acquired by DVC and OCT, so that the Young's modulus and Poisson's ratio of each layer can be solved individually by VFM.Strains are used as the input for the VFM calculation of material properties of double-layer materials under a compressive load.This quantitative 3D OCE method can obtain the elastic modulus of each Appl.Sci.2019, 9, 1349 3 of 16 layer of the specimen simultaneously.As many biological tissues, such as the skin and the vascular wall, are layered structures, this study will have a wide application potential in experimental biomechanics.

3D Strain Measurement by OCT Imaging and DVC
A swept source OCT (SSOCT) system and a compressive loading device were developed in our laboratory.A schematic diagram of the system is shown in Figure 1a.Laser emitted from a swept laser source (HSL-20-100-B, Santec, Aichi, Japan) goes to a Mach-Zahnder interferometer.The interference signal was detected by a balanced photodetector (EBR370006-02, Exalos, Schlieren, Switzerland) then acquired by a data acquisition card in a computer.The central wavelength of the laser is 1315 nm; the bandwidth is 88 nm, which results in an axial resolution of 9.93 µm.The refractive index of biological tissue is around 1.4, so the axial resolution in tissue is 7.09 µm.Sample beams are scanned by two Galvo mirrors.The line scan rate is 100 kHz, so it takes 0.46 s to scan a bulk of 214 × 214 × 1024 voxels.A compressive loading device consisting of a base and an optical window was designed.A force sensor (FSH03221, FUTEK, Irvine, CA, USA) is set on the base to measure the applied load.An example of a 2D cross-sectional OCT image and a 3D volumetric image of a phantom are shown in Figure 1b,c.tissues, such as the skin and the vascular wall, are layered structures, this study will have a wide application potential in experimental biomechanics.

3D Strain Measurement by OCT Imaging and DVC
A swept source OCT (SSOCT) system and a compressive loading device were developed in our laboratory.A schematic diagram of the system is shown in Figure 1a.Laser emitted from a swept laser source (HSL-20-100-B, Santec, Aichi, Japan) goes to a Mach-Zahnder interferometer.The interference signal was detected by a balanced photodetector (EBR370006-02, Exalos, Schlieren, Switzerland) then acquired by a data acquisition card in a computer.The central wavelength of the laser is 1315 nm; the bandwidth is 88 nm, which results in an axial resolution of 9.93 μm.The refractive index of biological tissue is around 1.4, so the axial resolution in tissue is 7.09 μm.Sample beams are scanned by two Galvo mirrors.The line scan rate is 100 kHz, so it takes 0.46 s to scan a bulk of 214 × 214 × 1024 voxels.A compressive loading device consisting of a base and an optical window was designed.A force sensor (FSH03221, FUTEK, Irvine, USA) is set on the base to measure the applied load.An example of a 2D cross-sectional OCT image and a 3D volumetric image of a phantom are shown in Figure 1b,c.A DVC algorithm was developed by integrating a coarse search mechanism, an IC-GN based fine search algorithm and a local ternary quadratic polynomial fitting method for strain calculation.Firstly, because the ZNCC correlation function is not sensitive to noise, a ZNCC correlation function based coarse search was developed to find initial integer values of voxel shifts.Secondly, the IC-GN [21] based fine search algorithm was employed.The IC-GN algorithm changes the shape and location of target and reference sub-volumes simultaneously to find the maximum of the correlation function.Because the Hessian matrix does not need to update in each iteration, the computational speed can be ~1.7 times faster than the conventional DVC method [21].The computing speed is ~4 s per point calculating OCT images by DVC.The convergence criterion was set as either the increment of each displacement component was less than 0.001 voxels, or the maximum iteration of 20 was reached.A cubic spline interpolation scheme with not-a-knot end condition was utilized for the calculation of sub-voxel intensities.More details about the IC-GN algorithm can be found in paper [21].Through the coarse-fine search DVC calculation, the 3D displacements can be obtained with sub-voxel resolution.Thirdly, the local ternary quadratic polynomial fitting equations were used to fit 3 × 3 × 3 voxels of the displacement data with the least square method.Then, the strains were calculated using Cauchy's formulae [24].The displacement and strain resolution of this system are around 0.06 voxels (~1 μm) and 0.2%, respectively, tested by the stationary experiment as shown in A DVC algorithm was developed by integrating a coarse search mechanism, an IC-GN based fine search algorithm and a local ternary quadratic polynomial fitting method for strain calculation.Firstly, because the ZNCC correlation function is not sensitive to noise, a ZNCC correlation function based coarse search was developed to find initial integer values of voxel shifts.Secondly, the IC-GN [21] based fine search algorithm was employed.The IC-GN algorithm changes the shape and location of target and reference sub-volumes simultaneously to find the maximum of the correlation function.Because the Hessian matrix does not need to update in each iteration, the computational speed can be ~1.7 times faster than the conventional DVC method [21].The computing speed is ~4 s per point calculating OCT images by DVC.The convergence criterion was set as either the increment of each displacement component was less than 0.001 voxels, or the maximum iteration of 20 was reached.A cubic spline interpolation scheme with not-a-knot end condition was utilized for the calculation of sub-voxel intensities.More details about the IC-GN algorithm can be found in paper [21].Through the coarse-fine search DVC calculation, the 3D displacements can be obtained with sub-voxel resolution.Thirdly, the local ternary quadratic polynomial fitting equations were used to fit 3 × 3 × 3 voxels of the displacement data with the least square method.Then, the strains were calculated using Cauchy's formulae [24].The displacement and strain resolution of this system are around 0.06 voxels (~1 µm) and 0.2%, respectively, tested by the stationary experiment as shown in Figure 1d,e, where the region of interest (ROI) is marked out with a blue rectangle in (b) and a blue cube in (c), which may become worse in deformation measurements due to the speckle blinking [25,26].

VFM
The VFM is based on the principle of virtual work, which has where σ is the real stress tensor; ε* is the virtual strain tensor; u* is the virtual displacement vector which should meets the boundary conditions; T is the external force vector on the surface of S f ; b is the body force vector; a is the acceleration vector; ρ is the density; V is the volume; S is the surface; ":" means the contraction operation of the second order tensor, i.e., σ:ε* = ∑ i,jσi,jε * i,j , i, j = 1, 2, 3.In this paper, the body force b can be neglected because it is relatively small compared with the external force T; the samples are compressed statically, so the acceleration a can also be neglected.Hence, the Equation ( 1) can be simplified as For convenience, the stress tensor and strain tensor are written as vectors: According to elastic mechanics, the constitutive equations of 3D isotropic elastic materials are where the relationships between the two independent elastic parameters Q 11 , Q 12 and Young's modulus E and Poisson's ratio ν are Then, Equation (2) becomes where Q ij is elastic constant.If the material is homogenous, the Equation ( 7) can be expressed as Appl.Sci.2019, 9, 1349 5 of 16 The integration can be approximated to summation by where k represents the kth test point; M is the number of total points; v k is the volume of the kth point or mesh; ε* i,k is the virtual strain of the ith variable of the kth point; ε j,k is the real strain of the jth variable of the kth point.There are two unknown parameters Q 11 and Q 12 to solve, so two virtual fields are needed.Superscript (i) which means the ith virtual field and i = 1, 2 is used in order to build linearly independent equation: where the coefficient matrix A is a 2 × 2 square matrix; the constitutive parameter vector Q and the virtual work vector B are 2 × 1 vectors.The coefficient matrix A can be expressed as follows: The constitutive parameter vector Q is The virtual work vector B can be expressed as Hence, the two unknown elastic constants can be solved by Equations ( 14) and (15): For double-layered material, the principle of virtual work is still adaptable, while there will be two sets of linear equations, assuming each layer is homogeneous.Equation ( 2) can be written as: where n is the layer number; n A is the coefficient matrix; n Q is the constitutive parameters vector, and n B is the virtual work vector done by external forces.The coefficient matrix n A can be expressed as follows: where n represents the nth layered material and n = 1, 2 for a double-layered material; V n represents the volume of the nth layered material.The constitutive parameter vector n Q is where n Q 11 and n Q 12 represent the two independent elastic parameters of the nth layer.The virtual work vector n B can be expressed as where S n means the entire surface, and n T represents the external force vector of the nth layer.S fn is the surface where the external force applies.n T includes the forces on the interfaces of the layers, which becomes external forces after the layers are divided during calculation.u* (i) is the ith virtual displacement vector.A flow chart of the layered VFM is shown in Figure 2. where Sn means the entire surface, and n T represents the external force vector of the nth layer.Sfn is the surface where the external force applies.n T includes the forces on the interfaces of the layers, which becomes external forces after the layers are divided during calculation.u* (i) is the ith virtual displacement vector.A flow chart of the layered VFM is shown in Figure 2.

Strain and Young's Modulus Measurement of a Homogeneous Phantom
3D displacements and strains of a phantom under compression were measured by the DVC-based 3D OCE technique.A homogenous phantom was made by mixing food-grade translucent 45-degree silica gel and titanium dioxide (TiO2) scatterer with a diameter of 1 μm and density of 0.5 mg/mL.The phantom was cut into a lx = 4.24 mm, ly = 4.24 mm and lz = 11.60 mm cube shown as Figure 3a and was set on the surface of the force sensor of the loading device shown in Figure 1a and Figure 3b.The phantom was preloaded for about 5 min until the reading of the force sensor became stable at 1.072 N.Then, the first 3D OCT image was taken as shown in Figure 3d.The phantom was compressed by 0.046 N immediately after the first image was taken, and the second 3D OCT was taken.This procedure took less than 1 s, thus, tissue relaxation was neglected.The phantom was compressed for −41.0 μm reading from the translation of the optical window.The phantom was scanned 4 mm in both the width and length directions.The ROI was marked out with a blue cube shown in Figure 3d.It is 64.

Strain and Young's Modulus Measurement of a Homogeneous Phantom
3D displacements and strains of a phantom under compression were measured by the DVC-based 3D OCE technique.A homogenous phantom was made by mixing food-grade translucent 45-degree silica gel and titanium dioxide (TiO 2 ) scatterer with a diameter of 1 µm and density of 0.5 mg/mL.The phantom was cut into a l x = 4.24 mm, l y = 4.24 mm and l z = 11.60 mm cube shown as Figure 3a and was set on the surface of the force sensor of the loading device shown in Figures 1a and 3b.The phantom was preloaded for about 5 min until the reading of the force sensor became stable at 1.072 N.Then, the first 3D OCT image was taken as shown in Figure 3d.The phantom was compressed by 0.046 N immediately after the first image was taken, and the second 3D OCT was taken.This procedure took less than 1 s, thus, tissue relaxation was neglected.The phantom was compressed for −41.0 µm reading from the translation of the optical window.The phantom was scanned 4 mm in both the width and length directions.The ROI was marked out with a blue cube shown in Figure 3d.All components of displacements and strains obtained by the DVC are shown in Figure 4.It can be seen from Figure 4c that the value of displacement in the z-direction decreases from the top to the bottom, which corresponds to theory.The mean value of the displacement on the top surface is −42.7 μm.So, the relative error for the z-direction displacement measurement is 4.1% compared with the imposed displacement of −41.0 μm.The displacements in the x and y directions shown in Figure 4a,b change diagonally, probably because the sample was slightly tilted and the OCT scanning direction was not strictly parallel with the edge of the sample.Each of the strains shown in Figure 4d-i is supposed to be a constant, and therefore the variations of the results demonstrate measurement errors.Figure 4f shows that the values are fairly constant except a thin top layer.The bigger values on the top portion, are possibly induced by the friction between the glass window and the top surface of the sample.The mean value of the εz of the whole ROI is −0.39% as shown in Figure 4f.While the applied strain is −41.0 μm/11.60 mm = −0.35%,thus the relative measurement error is 11.4%.If the top portion with bigger values is excluded, the relative measurement error becomes 5.7%.The measurement error for shear strains shown in Figure 4g-i is within 0.2%.
The full-field strain ε obtained in Figure 4d-i and virtual strains ε* (i) were then utilized to form the coefficient matrix A of VFM in Equation (11).Considering the frictions on the top and bottom surface of the phantom, the boundary conditions were as follows:   The full-field strain ε obtained in Figure 4d-i and virtual strains ε* (i) were then utilized to form the coefficient matrix A of VFM in Equation (11).Considering the frictions on the top and bottom surface of the phantom, the boundary conditions were as follows: where T 0 was the applied pressure, i.e., T 0 = −0.046N/(4.24 × 10 −3 m × 4.24 × 10 −3 m) = −2558.7Pa; z t and z b were the locations of the top and bottom surfaces of the phantom where the in-plane displacements were zeros.z t = 908.8µm was calculated according to the location of the white line in Figure 3c.z b = −10,691.2µm calculated by subtracting l z from z t .
According to the boundary conditions, the virtual displacement fields were chosen as follows: Hence, the virtual strain fields were The coefficient matrix A in Equation (11) were written as where M is the number of the matched points in DVC and M = 405; v k is the volume of the kth mesh (point) and v k = LxLyLz/M.According to Equation ( 13) and static equilibrium, the virtual work vector B is Then, coefficient matrix A and the virtual work vector B were input to Equation ( 14) to solve the elastic constants.Q 11 and Q 12 were calculated to be 2.2046 × 10 6 Pa and 1.1070 × 10 6 Pa.Young's modulus and Poisson's ratio were then calculated to be 1464.5 kPa and 0.334 respectively, by inputting Q to Equation (15).The Young's modulus of the phantom was also measured by a uniaxial tensile test, which resulted in 1471.0 kPa.The relative discrepancy between the results obtained by the two methods is ~1%.The Poisson's ratio, calculated by dividing the average ε z from the lateral strain obtained by DVC, was 0.306, where the the quadratic mean of ε x and ε y was taken as the lateral strain.The Poisson's ratio calculated by VFM was 0.334, which is 8% different to 0.306.

Strain and Young's Modulus Measurement of a Double-Layer Phantom
A double-layered phantom was made by mixing 1 µm and 0.5 mg/mL TiO 2 scatterer with 15-degree silica gel for the bottom layer and 0-degree silica gel for the top layer.The dimensions of the double-layer phantom were l x = 4.20 mm, l y = 4.00 mm and l z = 1.90 mm.A compression experiment was conducted following the same procedure as the homogeneous phantom.The double-layer phantom was preloaded for 1.320 N, then compressed by 0.247 N as read from the force sensor.The phantom was scanned 4 mm along both the width and length directions, before and after deformation.A 3D and a 2D OCT image of the phantom are shown in Figure 5a, from which the interface of the two layers can be observed.The ROI is marked with a blue cube or rectangle, which is 99.3 µm below the top surface of the phantom as shown in Figure 5aii.The average of the correlation coefficient of the ROI was tested to be 0.581.The coordinate was built with the origin O at the bottom left corner of the ROI, and the three axes were along the three edges of the ROI as shown in Figure 5ai.The number of computed points was M = 9 × 9 × 5 = 405.The step lengths were 17 voxels in the x and y directions, and 8 voxels in the z-direction.The computed dimensions were Lx = Ly = 2861.1 µm and Lz = 397.2µm.The 3D displacements of the ROI with the defined coordinate system are shown in Figure 5b-d.The 2D images are the cross-sections along the dashed lines on the 3D images.The distributions of 3D strains are shown in Figure 6.From the 3D display of ε z in Figure 6c, the interface of the two layers can be observed.Strain-gradients were then calculated to segment the two layers automatically.Based on the abrupt change of the strain-gradient, the interface can be determined, which is drawn with dashed lines in Figure 6ci,cii.The strain-gradients of ε z in the z-direction of the cross-section in Figure 6ci,cii are plotted on their cross-sectional OCT images as shown in Figure 7.For accuracy of computation in the VFM, broken lines simulated the physical interface between two layers, the blue curves shown in Figure 7.Although the physical boundary is recognizable in the phantom, the interfaces of most biological tissues are invisible.Hence, the automatic estimation of the boundary is necessary.The location of the minimum of the strain-gradient is where the interface of the two layers, i.e., the purple dashed lines as shown in Figure 7a,b, shows a good agreement with the physical boundary.
After the full-field strain ε of the double-layer phantom was obtained, the constitutive parameters were calculated by the VFM.Because the errors of ε x are large on the left boundary as shown in Figure 6a, the values of the left two columns were dropped out during VFM calculation.Then Lx became: Lx = 2225.3µm.The number of computed points became M = (9 − 2) × 9 × 5 = 315 in the VFM.The boundary conditions were the same as that of the homogenous phantom test, where the displacements in the x and y directions on the top and bottom surfaces were approximated as zeros, considering the friction on the contact surfaces.The value of z t = 496.5 µm obtained from OCT imaging as shown in Figure 5aii  The full-field strain ε and virtual strains ε* (i) were input to Equation (17).The two coefficient matrixes n A were The full-field strain ε and virtual strains ε* (i) were input to Equation (17).The two coefficient matrixes n A were where the superscript n on the upper left corner represents the nth layer and n = 1, 2, where 1 means the bottom layer and 2 means the top layer; M n is the number of meshes in the nth layer; v k is the volume of the kth mesh and v k = LxLyLz/M.One key step for the calculation is to segment the two layers according to the strain-gradient as described earlier, so as to determine the value of M n .After the two layers were divided as shown in Figure 6c, the M n were calculated to be M 1 = 143 and M 2 = 172 respectively.According to Equation ( 19) and the static equilibrium, the virtual work vectors n B are where N x and N y are the number of points in the x and y directions respectively; N d is the dth maximum number in the z-direction of the bottom layer, where d = 1, 2, . . ., N x N y .Then the coefficient matrixes n A and the virtual work vectors n B were substituted into the linear equations, Equation (16).The elastic constants n Q 11 and n Q 12 of the two layers were then solved as: 1 Q = (1.1920× 10 6 Pa, 8.8362 × 10 5 Pa) T ; 2 Q = (3.2683× 10 5 Pa, 2.4166 × 10 5 Pa) T .The Young's moduli and Poisson's ratios were respectively E 1 = 439.7 kPa, ν 1 = 0.426; E 2 = 121.4kPa, ν 2 = 0.425, where the subscript 1 means the bottom layer and the subscript 2 means the top layer.The Young's moduli were also measured by tensile tests.The comparison of the results between the VFM and the tensile tests is listed in Table 1, which shows that the relative discrepancy between the two measurements is less than 5%, which is 5% better than a simple plane boundary as we used before.The Poisson's ratio obtained by DVC was 0.422.The relative difference between the method of DVC and VFM is 1%.
where Nx and Ny are the number of points in the x and y directions respectively; Nd is the dth maximum number in the z-direction of the bottom layer, where d = 1, 2, …, NxNy.Then the coefficient matrixes n A and the virtual work vectors n B were substituted into the linear equations, Equation ( 16).The elastic constants n Q11 and n Q12 of the two layers were then solved as: 1 Q = (1.1920× 10 6 Pa, 8.8362 × 10 5 Pa) T ; 2 Q = (3.2683× 10 5 Pa, 2.4166 × 10 5 Pa) T .The Young's moduli and Poisson's ratios were respectively E1 = 439.7 kPa, ν1 = 0.426; E2 = 121.4kPa, ν2 = 0.425, where the subscript 1 means the bottom layer and the subscript 2 means the top layer.The Young's moduli were also measured by tensile tests.The comparison of the results between the VFM and the tensile tests is listed in Table 1, which shows that the relative discrepancy between the two measurements is less than 5%, which is 5% better than a simple plane boundary as we used before.The Poisson's ratio obtained by DVC was 0.422.The relative difference between the method of DVC and VFM is 1%.

Strain and Young's Modulus Measurement of Pork
A piece of pork was bought from a grocery store.To test the DVC and VFM based method for the elasticity measurement of biological tissue, a fresh pork specimen was tested at room temperature.A small cube shown in Figure 8a was resected for measurement, which has a layer of pig's skin on the top.Dimensions of the pork sample were lx = 10.36 mm, ly = 10.30mm and lz = 7.50 mm.A reference volumetric OCT image was taken when 0.620 N preload was applied.Then another volumetric OCT image was taken when the specimen was compressed by 41.0 μm and the applied load was 0.082 N. 6 mm in both the width and length directions were imaged by the OCT.

Strain and Young's Modulus Measurement of Pork
A piece of pork was bought from a grocery store.To test the DVC and VFM based method for the elasticity measurement of biological tissue, a fresh pork specimen was tested at room temperature.A small cube shown in Figure 8a was resected for measurement, which has a layer of pig's skin on the top.Dimensions of the pork sample were l x = 10.36 mm, l y = 10.30mm and l z = 7.50 mm.A reference volumetric OCT image was taken when 0.620 N preload was applied.Then another volumetric OCT image was taken when the specimen was compressed by 41.0 µm and the applied load was 0.082 N. 6 mm in both the width and length directions were imaged by the OCT.A cross-sectional OCT image of the specimen with the ROI marked by a blue rectangle is shown in Figure 8b.A reference 3D volumetric OCT image is shown in Figure 8c with the ROI marked by a blue cube, which is 461.8 µm below the top surface of the specimen as shown in Figure 8b.The average of the correlation coefficient of the ROI was tested to be 0.391.The coordinate was built with the origin O at the bottom left corner of the ROI, and the three axes were along the three edges of the ROI as shown in Figure 8c.The number of computed points was M = 31 × 31 × 8 = 7688.The step lengths were 7 voxels in the x, y and z directions.The computed dimensions were Lx = Ly = 4057.9µm and Lz = 556.1 µm.The 3D displacements and strains of the ROI with the defined coordinate system are shown in Figure 9.The displacements in the x and y directions are less than 6 µm.The values of displacement and the normal strain in the z-direction shown in Figure 9c and f decrease with the increase of the depth.From Figure 9f the layered structure can be identified as the double-layered phantom.The other strain components are relatively small, mostly less than 1.0% but show heterogeneity of the tissue.

Discussion
The 3D OCE method based on DVC and VFM has been developed in this study.Experimental results demonstrate that this method can obtain the 3D displacements, strains and constitutive parameters including Young's modulus and Poisson's ratio of layered materials.

Discussion
The 3D OCE method based on DVC and VFM has been developed in this study.Experimental results demonstrate that this method can obtain the 3D displacements, strains and constitutive After the full-field strain ε of the pork specimen was obtained, the constitutive parameters were derived by the same procedure as the double-layer phantom using the VFM.Some of the geometric parameters in this experiment were z t = 1017.9µm as shown in Figure 8b; z b = z t − l z = −6482.1 µm.T 0 = −0.082N/(10.36 × 10 −3 m × 10.30 × 10 −3 m)= −768 Pa.By solving Equation ( 16), the Young's moduli and Poisson's ratios of the two layers were obtained.They were E 1 = 123.9kPa, ν 1 = 0.083 for the lower layer, and E 2 = 31.62kPa, ν 2 = 0.123 for the top layer.These values are in a similar range to the porcine skin's Young's modulus of 113 kPa measured by Yeung et al. [27].

Discussion
The 3D OCE method based on DVC and VFM has been developed in this study.Experimental results demonstrate that this method can obtain the 3D displacements, strains and constitutive parameters including Young's modulus and Poisson's ratio of layered materials.
The Young's modulus of the homogenous phantom obtained by the DVC and VFM based 3D OCE method is ~1% different from the value measured by a tensile test, which demonstrates the effectiveness of the new method.For the compressive loading experiments conducted in this paper, the theoretical strain distribution in the depth direction and the shear strains can be estimated.These distributions and analyses show the relative error is about 10%.Errors are also observable at the edges shown in Figure 4g-i.Reasons for the strain measurement errors include: (1) the influence of the glass window on the top of the specimen including the friction with the specimen and the slant compression; (2) the boundary effect of the polynomial fitting; (3) the influence of noise, which is unavoidable in laser interferometry.The measurement can be improved by improving the loading device, optimizing the virtual fields and applying data processing techniques to reduce noise.The Poisson's ratio can be easily calculated from the DVC results or by the VFM.The results obtained by the two methods are very close, which can be further verified by other experimental methods such as tensile testing.
A complete friction boundary condition was employed at both the top and bottom surfaces.This is reasonable, because the maximum ratio of shear stress to compressive normal stress on the top surface of the phantom is ~0.02 at the four corners, calculated by simulation.Generally, the static friction coefficient that is the ratio of the maximum static shear stress to normal stress between the glass and different rubbers is larger than 0.2 [28], which would be similar to the coefficient between the silica gel phantom and the glass.Because the shear stress is 10 times less than the maximum static shear stress, the complete friction boundary conditions can be applied in Equation (20).
The OCT noise will affect the gray value of speckle images and then affect the accuracy of DVC and VFM.Above all, the OCT noise will affect the correlation coefficient in DVC.The correlation coefficient decreased from 0.8 to 0.6 with the depth increasing in the experiment of the homogeneous phantom.Because the scattering light intensity decreases with the increase of depth, the information of the sample is reduced dramatically.The OCT noise will cause the results of DVC to no longer be reliable in deep regions.Hence, only the shallow regions were used.Furthermore, the OCT noise can mask the gray value of speckle images.It will result in matching errors even in the samples without any motion or deformation as the stationary experiment showed that the sensitivities of our system were 0.05 voxels and 0.2% for displacement and strain measurement, respectively.If the matching errors exist, the strains inputted into the VFM will lead to errors in the calculation of the constitutive parameters.By combining noise resistant methods, such as the ZNSSD criterion [20], the local ternary quadratic polynomial fitting method, and the linear VFM method [23], the experimental results will be insensitive to noise.
The 3D strain maps not only provide the necessary input, but also enable more automatic and accurate segmentation of the double-layer materials.This is important for the successful performance of the proposed VFM method.Although the interface was estimated as a series of broken lines, it is reliable because most of the purple dashed line is located at the blue curve in Figure 7b.The interface can also be described more accurately in the future work to improve the computing accuracy through finer meshing.The method of solving constitutive parameters of different components respectively has the merits of less virtual fields, large-scale linear equations and better solutions, because the solutions of the elastic constants acquired by the method of solving all components simultaneously had a negative value.From Figures 6c and 9f, it can be seen that the strain of the top layer is greater than the bottom layers, which indicates that the top layers are softer.Therefore, the image of strain distribution may be very useful for the diagnosis of diseases, as it shows the relative stiffness and boundaries of different components of soft tissue.
As a high-resolution laser interferometry imaging technique, the signal to noise ratio of OCT images decreases with the increases of depth as shown in Figure 8b, and the imaging is easy to be disturbed.In this experiment, only a thin portion of ~0.56 mm of the pork image was processed, as the correlation coefficient becomes too low anywhere below.As mentioned above, the OCT images are generated by laser interferometry.The scattering characteristic and anisotropy factor of the tissue influences the size and distribution of speckles of OCT images [29] and the DVC calculation.The size of OCT speckle was determined by OCT resolution, optical devices and microstructure of tissues [30,31].To improve the accuracy of the DVC based 3D OCE deformation measurement, the influence of OCT speckle properties needs to be studied further.Kurokawa et al. [32,33] and Wijesinghe et al. [34] proposed a more sensitive displacement and strain detection method using complex cross-correlation.The sensitivity can be within tens of nanometers and tens of micro-strains, which improves about tenfold the DVC method.In [32,33] out of plane motion was measured from 2D images; in [34], particularly, 3D displacements were estimated from 3D OCT volumes.The using of complex data for correlation calculation may be also helpful for improving DVC calculation.
The pork tissue tested in this paper is assumed to be layered isotropic elastic and homogenous based on the strain maps shown in Figure 9.However, soft tissue is usually anisotropic.DVC has the advantage of obtaining the displacement and strain tensors of anisotropic materials, although the algorithm may need to be further improved.The VFM may become very complicated though.Image segmentation methods may be used to differentiate between various components of heterogeneous material, so that the VFM can be simplified.If the constitutive model of the tested tissue is too complicated or even unknown, finite element method updating or other inverse problem solving may be a good choice.

Conclusions
In conclusion, a DVC and VFM based OCE is developed to quantitatively measure the 3D displacements, strains and elasticity of double-layered biological tissue.The integrated noise-insensitive DVC method can obtain the 3D strain tensor with an accuracy of 10%.With the full-field strain as input, the VFM can recover the elastic modulus of homogeneous material with the accuracy of ~1%.The automatic segmentation method proposed based on the strain gradient can simplify the VFM calculation of layered samples.There are less virtual displacements and large-scale linear equations involved in the VFM, once a double-layered material is divided into two parts.The full-field strain tensor provides excessive information for elastic modulus quantification of the double-layered material by the VFM.The relative measurement error of the Young's modulus of the double-layered phantom is ~5% in this study, which can be improved by optimizing the experimental setup and the virtual fields.The simultaneous elasticity quantification of each of the two layers of the pork specimen, shows great promise for the proposed method for mechanical properties measurement of biological tissues.In all, the noise insensitive DVC and VFM based OCE method developed in this study can obtain the 3D displacements, strains, Young's modulus and Poison's ratio of layered materials simultaneously.It can be a powerful tool for the differentiation of various components of heterogeneous biomaterials and for biomechanics measurement.experimental data.C.L. provided the experimental device and gave some valuable opinions for the methodologies.The programming codes were written by J.W. and F.M. J.C. and C.S. administrated and supervised this project.The original draft of this paper was written by F.M., and the revision of the paper was completed by C.S., J.W., X.Z. and F.M.

Figure 1 .
Figure 1.The SSOCT and compressive loading system.(a) A schematic diagram of the system; (b) a 2D OCT cross-sectional image (B-scan); (c) a 3D volumetric OCT image.Stationary experiment: (d) 3D displacements in the x, y and z directions of the ROI; (e) 3D strain tensors of the ROI.The ROI is marked out with a blue rectangle in (b) and with a blue cube in (c).

Figure 1 .
Figure 1.The SSOCT and compressive loading system.(a) A schematic diagram of the system; (b) a 2D OCT cross-sectional image (B-scan); (c) a 3D volumetric OCT image.Stationary experiment: (d) 3D displacements in the x, y and z directions of the ROI; (e) 3D strain tensors of the ROI.The ROI is marked out with a blue rectangle in (b) and with a blue cube in (c).

Figure 2 .
Figure 2. A flow chart of the layered VFM.
7 μm below the top surface of the phantom as shown in Figure 3c.The coordinate was built with the origin O at the bottom left corner of the ROI, and the three axes were along the three edges of the ROI as shown in Figure 3d.The number of DVC computed points was M = 9 × 9 × 5 = 405.The step lengths in the x, y and z directions were 17 voxels.By calibration, the dimensions of the ROI were Lx = Ly = 2861.1 μm, and Lz = 844.1 μm.The average of the correlation coefficient of the ROI was tested to be 0.659.

Figure 2 .
Figure 2. A flow chart of the layered VFM.
It is 64.7 µm below the top surface of the phantom as shown in Figure 3c.The coordinate was built with the origin O at the bottom left corner of the ROI, and the three axes were along the three edges of the ROI as shown in Figure 3d.The number of DVC computed points was M = 9 × 9 × 5 = 405.The step lengths in the x, y and z directions were 17 voxels.By calibration, the dimensions of the ROI were Lx = Ly = 2861.1 µm, and Lz = 844.1 µm.The average of the correlation coefficient of the ROI was tested to be 0.659.

Figure 3c .
Figure 3c.The coordinate was built with the origin O at the bottom left corner of the ROI, and the three axes were along the three edges of the ROI as shown in Figure 3d.The number of DVC computed points was M = 9 × 9 × 5 = 405.The step lengths in the x, y and z directions were 17 voxels.By calibration, the dimensions of the ROI were Lx = Ly = 2861.1 μm, and Lz = 844.1 μm.The average of the correlation coefficient of the ROI was tested to be 0.659.

Figure 3 .
Figure 3.A sample of homogeneous phantom.(a) A photo of the phantom; (b) a photo of the phantom under compression; (c) a 2D OCT image of the phantom, where the ROI is marked with a blue rectangle; (d) a 3D OCT image of the phantom, where the ROI is marked with a blue cube.All components of displacements and strains obtained by the DVC are shown in Figure4.It can be seen from Figure4cthat the value of displacement in the z-direction decreases from the top to the bottom, which corresponds to theory.The mean value of the displacement on the top surface is −42.7 µm.So, the relative error for the z-direction displacement measurement is 4.1% compared with the imposed displacement of −41.0 µm.The displacements in the x and y directions shown in Figure4a,b change diagonally, probably because the sample was slightly tilted and the OCT scanning direction was not strictly parallel with the edge of the sample.Each of the strains shown in Figure4d-i is supposed to be a constant, and therefore the variations of the results demonstrate measurement errors.Figure4fshows that the values are fairly constant except a thin top layer.The bigger values on the top portion, are possibly induced by the friction between the glass window and the top surface of the sample.The mean value of the ε z of the whole ROI is −0.39% as shown in Figure4f.While the applied strain is −41.0 µm/11.60 mm = −0.35%,thus the relative measurement error is 11.4%.If the top portion with bigger values is excluded, the relative measurement error becomes 5.7%.The measurement error for shear strains shown in Figure4g-i is within 0.2%.

Figure 3 .
Figure 3.A sample of homogeneous phantom.(a) A photo of the phantom; (b) a photo of the phantom under compression; (c) a 2D OCT image of the phantom, where the ROI is marked with a blue rectangle; (d) a 3D OCT image of the phantom, where the ROI is marked with a blue cube.
T0 was the applied pressure, i.e., T0 = −0.046N/(4.24 × 10 −3 m × 4.24 × 10 −3 m) = −2558.7Pa; zt and zb were the locations of the top and bottom surfaces of the phantom where the in-plane displacements were zeros.zt = 908.8μm was calculated according to the location of the white line in Figure 3c.zb = −10,691.2μm calculated by subtracting lz from zt.

Figure 4 .
Figure 4. 3D deformation measurement of the homogeneous phantom.(a-c) 3D displacements in the x, y and z directions respectively; (d-f) 3D normal strains in the x, y and z directions respectively; (g-i) 3D shear strains.

Figure 4 .
Figure 4. 3D deformation measurement of the homogeneous phantom.(a-c) 3D displacements in the x, y and z directions respectively; (d-f) 3D normal strains in the x, y and z directions respectively; (g-i) 3D shear strains.

Figure 5 .
Figure 5. 3D displacement measurement of a double-layer phantom.(a) OCT images: (a-i) a 3D OCT image with the ROI marked by a blue cube; (a-ii) a 2D OCT image with the ROI marked by a blue rectangle; (b-d) 3D displacements in the x, y and z directions respectively.Two displacement maps on the right side of each 3D image are the cross-sections along the dashed lines on (b-d).

Figure 5 .
Figure 5. 3D displacement measurement of a double-layer phantom.(a) OCT images: (a-i) a 3D OCT image with the ROI marked by a blue cube; (a-ii) a 2D OCT image with the ROI marked by a blue rectangle; (b-d) 3D displacements in the x, y and z directions respectively.Two displacement maps on the right side of each 3D image are the cross-sections along the dashed lines on (b-d).

Figure 6 .
Figure 6.3D strain fields of a double-layer phantom under compression.(a-c) Normal strains in the x, y and z directions respectively.The purple dashed lines in (c-i) and (c-ii) indicate the interface of the two layers.(d-f) Shear strains.Two strain maps shown underneath each of the 3D displays are the images of the cross-sections along the dashed lines.

Figure 6 .
Figure 6.3D strain fields of a double-layer phantom under compression.(a-c) Normal strains in the x, y and z directions respectively.The purple dashed lines in (c-i) and (c-ii) indicate the interface of the two layers.(d-f) Shear strains.Two strain maps shown underneath each of the 3D displays are the images of the cross-sections along the dashed lines.

Figure 7 .
Figure 7. Strain-gradient of εz in the z-direction.(a) and (b) are 2D strain-gradient maps of the cross-section (c-i) and (c-ii) in Figure 6 respectively plotted on their 2D OCT cross-sectional images.The purple dashed line and blue solid line indicate the estimated interface and physical interface of the two layers respectively.

Figure 7 .
Figure 7. Strain-gradient of ε z in the z-direction.(a,b) are 2D strain-gradient maps of the cross-section (c-i) and (c-ii) in Figure 6 respectively plotted on their 2D OCT cross-sectional images.The purple dashed line and blue solid line indicate the estimated interface and physical interface of the two layers respectively.

Figure 8 .
Figure 8. Pork sample.(a) A photo of the pork sample; (b) a central cross-sectional OCT image before deformation, where the ROI is marked with a blue rectangle; (c) a 3D OCT image before deformation, where the ROI is marked with a blue cube.

Figure 9 .
Figure 9. 3D deformation measurement of a pork specimen.(a-c) 3D displacements in the x, y and z directions respectively; (d-f) 3D normal strains in the x, y and z directions respectively; (g-i) 3D shear strains.

Figure 8 .
Figure 8. Pork sample.(a) A photo of the pork sample; (b) a central cross-sectional OCT image before deformation, where the ROI is marked with a blue rectangle; (c) a 3D OCT image before deformation, where the ROI is marked with a blue cube.

Figure 8 .
Figure 8. Pork sample.(a) A photo of the pork sample; (b) a central cross-sectional OCT image before deformation, where the ROI is marked with a blue rectangle; (c) a 3D OCT image before deformation, where the ROI is marked with a blue cube.

Figure 9 .
Figure 9. 3D deformation measurement of a pork specimen.(a-c) 3D displacements in the x, y and z directions respectively; (d-f) 3D normal strains in the x, y and z directions respectively; (g-i) 3D shear strains.

Figure 9 .
Figure 9. 3D deformation measurement of a pork specimen.(a-c) 3D displacements in the x, y and z directions respectively; (d-f) 3D normal strains in the x, y and z directions respectively; (g-i) 3D shear strains.

Funding:
This research was funded by [National Natural Science Foundation of China] grant number [11602166] and [Natural Science Foundation of Tianjin] grant number [16JCYBJC40500].The APC was funded by [National Natural Science Foundation of China] grant number [11602166].

Table 1 .
The comparison between the tensile test and calculated constitutive parameters of a double-layer phantom by VFM.

Table 1 .
The comparison between the tensile test and calculated constitutive parameters of a double-layer phantom by VFM.