## 1.

## Introduction

Minimally invasive surgery has revolutionized surgical intervention. It enables a quicker recovery for the patient related to reduced postoperatory stress. Moreover, it offers some aesthetic advantages in particular for exposed parts of the human body, e.g., the face. However, these advantages come at the expense of losing the sense of touch, as palpation is an important tool for tumor localization in open surgery. The localization process is based on the increased stiffness of tumorous tissue compared to healthy tissue. Preoperational data acquired from computer tomography (CT), magnetic resonance tomography (MRT), and ultrasound imaging can reveal information about the size of the tumor and location of the transition between malignant and healthy regions within the organ, but do not account for any changes in the location or shape of the organ due to increased deformability of soft tissue in comparison to hard tissue, different hydration levels of the body compared to the time when the preoperatory data were collected or changed position of the body (e.g., supine versus lateral recumbent). Furthermore, the application of pressure, which is needed for expanding a body cavity (e.g., thorax) in order to provide the surgeon with more space for carrying out the operation, results in a deformed shape and location of the organ under investigation. The goal of the minimally invasive intervention is to remove the tumor under best possible preservation of important functional tissue such as blood vessels and nerves.

Hence, there exists an unmet need for the development of sensor systems, which supports the surgeon during minimally invasive surgery in tumor localization.

Elastography is an imaging technique that enables the recovery of quantitative elastic parameters such as strain, stress, elastic modulus (linear elastic), or shear modulus (nonlinear elastic) from displacement measurements. Elastography has proven to be a successful complementary technique with increased sensitivity and specificity for cancer detection. This is based on the fact that tumorous tissue is 7 to 14 times stiffer in comparison to healthy tissue, according to Ref. 1. An increase of specificity from 61% to 79% was reported in Ref. 2 for the combination of elastography with ultrasound imaging. Specificity expresses the ratio of true negatives in comparison to the sum of true negatives and false positives.

The goal of this paper therefore is the development of an elastographic imaging technique in combination with finite element (FE) modeling to recover quantitative elastic parameters of tissue, which can be applied in minimally invasive surgery to support the surgeon during the three-dimensional (3-D) localization and discrimination process.

The paper is structured in the following manner, after presenting the background, the methodology and measurement principle will be explained in more detail. The FE model used will be presented, followed by the experimental setup based on rolling indentation and the results achieved for the investigation of the silicone phantom. The experimentally obtained and simulated two-dimensional (2-D) displacement fields and strain maps are compared to demonstrate that a realistic FE model has been generated, which can be used to solve the direct problem and hence the recovery of elastic parameters, which can directly not be measured. It will be shown how the strain contrast can be improved and strain ambiguities be reduced when taking into account many different strain maps, which are generated in the rolling indentation process. Afterwards, the 3-D distribution for different stress tensor elements and locations of the indenter will be represented. Moreover, the successful application of the image correlation technique to a biological sample will be shown.

## 2.

## Background

The working principle of elastography relies on the comparison of at least two different states of the object under investigation (mechanically loaded and unloaded). The loaded state can be obtained when applying an external force or body internal forces, which are mostly associated with the cardiovascular system, such as blood pressure change.^{3} The latter is limited to the investigation of blood vessels or tissue in close proximity only. Under load, tumorous tissue deforms less than healthy soft tissue due to its increased stiffness. By comparison of the two data sets, loaded and unloaded, the elastic properties of the tissue can be obtained and tissue abnormalities be localized. In that manner small tumors and/or enlarged lymph nodes hidden underneath a layer of tissue such as peritoneum, which may not have been registered in the preoperatory data, can be found, resulting in an increased accuracy of the diagnosis. Elastography has successfully been combined with some of the most established imaging modalities such as ultrasound imaging,^{2} OCT imaging,^{4} MRI,^{5} and CT.^{6} Digital image correlation for the dynamic analysis of biomechanical systems such as bones, tendons, and muscles was implemented in Ref. 7. The application of an endoscopic electron speckle pattern interferometry) system for deformation measurement of porcine kidney is described in Ref. 8. Digital holographic distal endoscopy for 3-D shape and deformation measurement has been reported in Ref. 9. Recently, the application of photoacoustic tomography for tissue investigation and discrimination was demonstrated in Ref. 10. In Ref. 11, it could be demonstrated that recording speed and optical resolution can be increased when employing digital double exposure holography for the investigation of the deformations caused by acoustic waves. In Refs. 12 and 13, digital image correlation elastography has been used.

In addition to the different imaging modes used at different resolution levels (ultrasound, CT, MRT, and visible light-based optical techniques), elastography is not only restricted to macroscopic organ and tissue level but can furthermore be applied at a microscopic cell level. A recent technique for the determination of elastic behavior of cells is Brillouin microscopy, which is based on inelastic scattering similar to Raman spectroscopy. In this case, the application of an external force results in a spectral shift. The amount of spectral shift can be related to stiffness properties.^{14}^{,}^{15}

Furthermore, atomic force microscopy has been applied for the investigation of individual cells, whereas the tip of the scanning diamond needle penetrates the specimen vertically while recording the force and displacement.^{16} It is worth noticing that, in contrast to macroscopic properties of tissues and organs, at microscopic level cancerous cells have a smaller Young’s modulus than healthy cells.^{17}^{,}^{18} The different elastic behavior of tissue and organs can be explained by a more densely meshed extracellular matrix of cancerous tissue in comparison to healthy tissue, resulting in an increased stiffness.^{19}^{,}^{20} The reduced stiffness of cancerous cells on the other hand may facilitate migration and invasion of cancerous cells in comparison to healthy cells. The increased stiffness of the extracellular matrix likewise promotes cell migration due to a feedback between forces.

In elastography, different types of loading can be characterized, according to Ref. 21, as quasistatic, harmonic, and transient. The quasistatic loading is the most widely used and has been applied initially in ultrasound elastography and was later extended to other elastographic imaging modalities. Harmonic loading consists of a low frequency acoustic wave transmitted within the tissue using a sinusoidal vibration source in contact mode. The local shear modulus can be obtained directly from the velocity of the propagating shear wave and the density of the specimen, as discussed in Ref. 22 for homogeneous tissue. However, shear waves attenuate rapidly in biological soft tissue, which limits the penetration depths. The problem can be overcome using transient loading, which is based on an acoustic radiation force that is locally focused within the specimen. The loading of the transient source is performed in a noncontact manner.

## 3.

## Methodology

The requirements imposed on the elastographic measurement system are recovery of cancer relevant elastic parameters, clinically applicable (no harm to patient), an endoscope implementable measurement principle, acquisition of data in real time, optical performance parameters as encountered in endoscopy (field-of-view typically $50\times 50\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ and sub-mm optical resolution^{23}) and robustness against disturbing environmental conditions (air flow and vibrations). Digital image correlation elastography satisfies all these requirements. Furthermore, it can be implemented using a visible light camera, which not only reduces setup costs but also enables the extraction of other morphological features relevant to cancer, such as hypervascularization.^{24} The hypervascularization effect is characterized by an increased density of blood vessels with a randomly distributed structure and irregular branching.^{24}

## 3.1.

### Digital Image Correlation

Digital image correlation is based on the comparison of two images at different states (loaded and unloaded). It is a well-known technique, widely applied in digital particle imaging velocimetry,^{25} in order to study fluid flow^{26} and air flow.^{27} More recently it has also been applied in biomedical imaging in order to obtain displacement vector fields.^{28}^{,}^{29}

Each image is divided in overlapping subimages, also known as interrogation areas. A displacement vector is obtained by comparing the corresponding subimages of the two states using crosscorrelation, as shown in Fig. 1. Extraction of the displacement vector via crosscorrelation can computationally best be accomplished via calculation of a 2-D FFT of both corresponding interrogation areas, computation of cross product of the corresponding FFTs of the first interrogation area, the resulting conjugate of the second interrogation area, and determination of mean displacement vector from location of maximum with respect to the center.

Design rules that have to be taken into account for digital image correlation are:

– The corresponding size of the interrogation area should be larger than the diameter of the markers that are tracked (e.g., speckles and particles).

– Adjacent interrogation areas need to overlap in order to ensure that marker points that are positioned close to the edge of the interrogation area can still be tracked.

– Homogeneously distributed fine structural details with a density of at least seven marker points per interrogation area need to be available in order to obtain a detection probability of 95%.

^{30}– The in-plane motion difference between individual marker points should be smaller than the diameter of a marker point.

^{30}– The out of plane motion of individual marker points should be smaller than quarter the depth of field.

^{30}– The maximum in-plane displacement of less than a quarter of the interrogation area employed.

^{30}– The maximum in-plane rotation should not exceed 10 deg.

^{31}– The maximum marker point displacement of half the size of a marker point during the camera integration time.

^{32}

The spatial resolution in digital image correlation matches the pixel size but can be improved successively, via the use of interpolation routines, to 0.001 pixels, as discussed in Ref. 28. The spatial separation of displacement vectors is determined by the size of the interrogation area and the overlap between consecutive interrogation areas, e.g., a $20\times 20\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ interrogation area with 50% overlap results in a density of 10 pixels (for every 10 pixels a new displacement vector is calculated). From the displacement, the strain that represents a quantitative elasticity parameter can be calculated. The strain is a tensor that is defined as the partial derivative of the displacement.

## (1)

$$\u03f5=\left[\begin{array}{ccc}{\u03f5}_{xx}& {\u03f5}_{xy}& {\u03f5}_{xz}\\ {\u03f5}_{yx}& {\u03f5}_{yy}& {\u03f5}_{yz}\\ {\u03f5}_{zx}& {\u03f5}_{zy}& {\u03f5}_{zz}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial {u}_{x}}{\partial x}& \frac{1}{2}(\frac{\partial {u}_{x}}{\partial y}+\frac{\partial {u}_{y}}{\partial x})& \frac{1}{2}(\frac{\partial {u}_{x}}{\partial z}+\frac{\partial {u}_{z}}{\partial x})\\ \frac{1}{2}(\frac{\partial {u}_{y}}{\partial x}+\frac{\partial {u}_{x}}{\partial y})& \frac{\partial {u}_{y}}{\partial y}& \frac{1}{2}(\frac{\partial {u}_{y}}{\partial z}+\frac{\partial {u}_{z}}{\partial y})\\ \frac{1}{2}(\frac{\partial {u}_{z}}{\partial x}+\frac{\partial {u}_{x}}{\partial z})& \frac{1}{2}(\frac{\partial {u}_{z}}{\partial y}+\frac{\partial {u}_{y}}{\partial z})& \frac{\partial {u}_{z}}{\partial z}\end{array}\right].$$^{34}The unit length is usually represented by the pixel-size in the object space. In our case, two parameters of the strain tensor or strain vector can be calculated, which can be derived from the displacements in the $x$ and $y$ directions. In our case, the indenter is actuated in the $y$ direction; therefore, the largest displacement occurred along this direction. Therefore, the corresponding elastic parameters along the $y$ direction will predominantly be considered in the following.

## 3.2.

### Inverse Problem and Three-Dimensional Finite Element Modeling-Based Solution

In addition to the 2-D strain map, which can be recovered from the displacement field, other elasticity parameters can be obtained only at the interface between indenter and specimen. Knowledge of the indenter geometry, the force applied, and the displacement obtained at the indenter position enables the recovery of elastic parameters such as stress and shear modulus for only this point. However, the goal is to retrieve the elastic parameters for all imaged object points and to show their corresponding 3-D distribution. The 3-D information is essential in order to locate the tumor in minimally invasive surgery and hence to support the navigation process. Hence, a solution to an under-defined system needs to be found in order to recover the spatial distribution of elastic parameters. An inverse problem has to be solved, which can be accomplished via solving a direct problem. The underlying physical laws applicable to both, inverse and direct problem, are the same, but inverted, and relate the elastic properties to a measurable mechanical response. The direct problem is well defined with respect to the geometrical (dimensions of tissue and foreign body, indenter geometry, and force) and elastic input parameters (elastic modulus for linear elastic, shear modulus, and locking stretch for nonlinear elastic). It starts with the desired elastic parameters as input parameters, whereas the experimentally obtained parameters (displacement field and strain map) represent the output, which can be compared with the measured values, see Fig. 2. Hence, it is an inverted process of the inverse problem and is therefore referred to as a forward problem.^{21} In that manner, the desired unknown elastic parameters, which experimentally cannot be retrieved, can be obtained from a realistically well-matched FE model that enables generating the computational solution of the direct problem. Solving the direct problem can be implemented using a noniterative and iterative approach. The noniterative approach requires a good estimate of the input parameters that can be generated from *a priori* knowledge based on stress–strain measurements about the different materials/tissues and the geometry of the indenter and object under investigation. The generation of relevant data to apply the noniterative method for medical applications can be obtained from previous operations on patients of similar conditions (age, gender, and life style) or preoperational data. Moreover, the geometry/morphology of inner organs such as kidney and the locations of abnormalities (cancer) can be obtained from preoperational data such as CT and MRT data.

The noniterative approach enables a quicker recovery of elastic parameters in combination with computational methods such as the FE method. The iterative approach on the other hand requires less input parameters; usually only the displacement field and the general elastic behavior (nonlinear hyperelastic behavior is very common among soft tissue^{35}) of the specimen are provided. The loading and elastic parameters are iteratively adjusted to result in a good match for the modeled and measured displacement field. Under the assumption that the same constitutive model is applied, the iterative approach is more time consuming, e.g., 6000 iterations have been employed in Ref. 36 to obtain a good match.

Many FE models of different dimensionality and elastic properties have been used for the implementation of the direct problem, among others a 2-D hyperelastic Arruda Boyce FE model,^{1} a 3-D axisymmetric FE model of linear elastic behavior,^{34} and a four-dimensional (4-D) hyperelastic model taking into account the temporal component.^{37} A detailed discussion of the many different existing models used in elastography to implement the direct problem can be found in Ref. 21. Since our measurement principle is based on rolling indentation, a 4-D hyperelastic FE model has been generated, taking into account the three spatial dimensions and the changing location of the indenter over time. The chosen constitutive FE model to simulate the hyperelastic behavior of soft tissues is the Arruda–Boyce model.^{38} Besides the excellent modeling of the nonlinear elasticity behavior of soft tissue, it furthermore offers the advantage of requiring only two parameters, shear modulus G and locking stretch ${\lambda}_{m}$. The shear modulus is defined as the ratio between shear stress and shear strain whereas the locking stretch corresponds to the value of the chain stretch at which the chain length reaches its fully extended state. A more detailed discussion can be found in Ref. 38.

## (2)

$${\sigma}_{0}=G[({\u03f5}^{2}-{\u03f5}^{-1})+\frac{1}{5{\lambda}_{m}^{2}}({\u03f5}^{4}-{\u03f5}^{-2})+\frac{11}{175{\lambda}_{m}^{4}}({\u03f5}^{6}-{\u03f5}^{-3})+\frac{19}{875{\lambda}_{m}^{6}}({\u03f5}^{8}-{\u03f5}^{-4})+\frac{519}{\mathrm{67,375}{\lambda}_{m}^{8}}({\u03f5}^{10}-{\u03f5}^{-5})].$$## 3.3.

### Silicone Phantoms and Corresponding FE Models

As a proof of concept, before examining the morphologically highly complex kidney, a silicone sample of similar dimension and hyperelastic behavior has been created, which is common practice in optical elastography.^{1}^{,}^{36}^{,}^{34} In order to generate a realistic silicone phantom, the elastic behavior of porcine kidney has first been studied in an uniaxial compression test, carried out on a commercial microindenter Zwicki line 500 (Roell GmbH, Ulm, Germany). The hyperelastic behavior mentioned in Ref. 35 could be confirmed, see Fig. 3(a). Afterward, various materials such as, Agar-Agar and different types of silicone ZA-OF, ZA-OO (O ShA) (ordered from Polymerschmiede GmbH), and PDMS Dow Corning 184 (ordered from Biesterfeld Spezialchemie GmbH) have been investigated to mimic the biological tissue. It was found that ZA-OO (O ShA) was best suited for mimicking the healthy soft tissue and PDMS Dow Corning 184 for mimicking the pathological tissue. By comparison with the force indentation curve for a porcine kidney shown in Fig. 3(a), a similar hyperelastic behavior and viscoelastic behavior of the silicone sample, shown in Figs. 3(b) and 3(c) could be confirmed. Moreover, the silicone samples also enable good reproducibility and temporal stability (not the case with Agar-Agar). The stress–strain curves obtained, as shown in Fig. 3(c), could be used as *a priori* information in order to generate a realistic FE model.

In Fig. 3(d), the dimensions of the phantom with a foreign body are shown. In correspondence to the phantom, an FE-model with a centrally located cylindrical inclusion of increased shear modulus (10 mm diameter) has been generated. The same outer dimensions have been used to produce another homogeneous silicone phantom to model an entirely healthy organ, which will be used to improve the strain contrast, which will be discussed in more detail at a later stage in this paper. The coordinates of indentation with respect to the edge of the phantom (8 mm), and the geometry of the indenter (two miniature ball bearings with 6 mm diameter and 3 mm width), which introduces the loading, have been taken into account for the generation of a realistic 3-D FE model. The entire FE model is composed of 110,565 elements. The size of the FEs ranged from $0.3\times 0.3\times 0.3\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{3}$ to $5\times 5\times 5\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{3}$. Small mesh sizes were employed in contact areas, i.e., in the contact area between the indenter and the tissue model and in the contact area between the homogeneous tissue and the foreign body (see Fig. 4).

In order to obtain a realistic modeling of the elastic behavior of malignant and normal tissues, different combination ratios between hard (PDMS 184) and soft-silicone (ZA-OO) have been generated and the force-indentation curves have been obtained from uniaxial compression tests applied to cylindrical silicone samples, as shown in Fig. 3(b). The stress–strain curves for soft and hard silicone were used as input parameters for the FE model. A 100% soft-silicone (ZA-OO) was chosen for mimicking normal tissue and a 100% hard-silicone (PDMS Dow Corning 184) for malignant tissue. The shear modulus G and locking stretch ${\lambda}_{m}$ for soft silicone ($G=0.73\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{kPa}$ and ${\lambda}_{m}=7.7$) and hard silicone ($G=7.73\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{kPa}$ and ${\lambda}_{m}=7$) obtained from the stress–strain curve [generated by uniaxial compression test, see Fig. 3(b)] have been used as input parameters for the FE model. These values are within the naturally occurring range according to Ref. 1 (malign tissues are 7 to 14 times stiffer in comparison with normal tissue).

Figure 4 shows the 3-D displacement field for a centrally located indenter. It took 2.5 h to compute the FE model, which was run on the bwUniCluster (one node and seven processors).

## 4.

## Experimental Results

## 4.1.

### Elastographic Measurements

In order to apply a well-known force and indentation over the entire length of the organ within a minimum amount of time, a two-axis force rolling indenter was manufactured in the institute engineering workshop, see Fig. 5(a).

The rolling movement is well suited for minimal invasive surgery, based on the fact that large areas can be scanned within a short time, in comparison to a pointwise measurement mode. Rolling indentation is related to the “quasi static loading,” as described in Ref. 21, albeit with temporally changing local coordinates of the indenter.

The force rolling indenter supports two measurement modes, namely constant indentation and constant force. Both are made possible via a feedback loop between the force-meter (Alluris FMI-B30BI, measurement range: 0 to 10 N, measurement uncertainty: 15 mN), and the stepper motor driving the $z$-axis (range: 19 mm, positioning uncertainty: $10\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$). The rolling indenter can be displaced laterally by 156 mm with a positioning uncertainty of $50\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$.

Furthermore, the device enables quick change of indenters with different geometries and sizes.

The aim is to use a small indenter in order to reduce the influence of disturbing artifacts, such as inhomogeneities that lie outside the imaged area and of an uneven measurement surface, on the measurement result.^{39} However, the size of the indenter cannot be too small since it may otherwise pierce into the sample instead of rolling across it. Moreover, in order to carry out the digital image correlation technique, the indenter size needs to match the density of the structural details available at the specimen. An indenter too small introduces strongly localized deformation, which at a small density of structural details, may not be detectable via digital image correlation, as pointed out in the design rules for digital image correlation in Sec. 3.1 based on Ref. 30 (seven markers per interrogation area).

At the next stage, a silicone phantom with a cylindrical inhomogeneity mimicking malign tissue, in accordance with the FE model has been generated. For mimicking normal tissue 100% soft-silicone (ZA-OO) was chosen and for malign tissue 100% hard-silicone (PDMS Dow Corning 184). The front surface of the silicone phantoms was covered with markers consisting of ash particles with different sizes ($300\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ to 1.3 mm) and shades (light gray to black), to enable a dense distribution of tracers while ensuring a good discrimination between adjacent markers. A strainer was used to restrict the particle size to 1.3 mm. Prior to the application of the ash particles the silicon phantom has been cleaned by water, resulting in a stronger adhesion of the ash particles. The ash particles have then been dispensed on the phantom. Afterward the front surface of the phantom was oriented toward the floor and shacked several times to remove those particles, which are not well fixed on the surface of silicon phantom. These particles would otherwise during the rolling indentation process fall off and hence result in measurement errors.

The silicone sample was illuminated with a white light source. Two miniature ball bearings with 6 mm diameter and 3 mm width have been used for the indenter, as shown in Fig. 5(a) enlarged image on right-hand side. The indentation depth employed was 3 mm. The indenter was placed at a distance of 8 mm to the front surface of the silicone phantom. The rolling indenter was moved at a speed of $2.4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}/\mathrm{s}$ across the specimen. A digital camera (Zyla 5.5, $2560\times 2160\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$, $6.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ pixel size) in combination with a Tamron 08831 camera objective was used to record the images. The resolution of the system is in our case not determined by the camera objective, which can deliver an NA of 0.03 ($20\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ resolution) at 200 mm recording distance, but by the magnification of the camera objective of 0.115 and the sensor’s pixel-size, resulting in resolvable features of $56\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$. In that manner renal blood vessels (interlobular arteries and veins) with perimeter size ranging from $\text{150\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ to $200\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$^{40} can well be resolved and imaged. The field of view recorded by the camera was $144\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 121\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$.

In order to save memory and more importantly postprocessing time, the images have been cropped to $1800\times 520\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$, which corresponds to a field of view of $101\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 29\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$. Two-hundred images have been recorded with an exposure time of 0.04 s (25 fps). This exposure time was well chosen in order to ensure that motion artifacts of the rolling indenter do not influence the measurement. According to Ref. 32, the lateral movement in the image plane is limited to half the size of the smallest imaging detail. In our case, the rolling indenter moved roughly $100\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ during the exposure time, which is perfectly suitable for the particles employed in our experiments. Hence, two different object states can be recorded within a time frame of less than $100\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{s}$.

In addition to the task of applying a load, the indenter can also be used to measure the force along the direction of indentation, as shown in Fig. 5(b). In this manner, the indenter fulfills the functionality of an elastographic sensor, which in combination with the 2-D displacement field and its derived elastic parameters, can be used for 3-D location of the elastic inhomogeneities.

The recorded images were analyzed using PIVLab.^{41} An interrogation area of $65\times 65\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$ was used with an overlap of 50%. Hence, displaced marker points that are positioned at the edge of the interrogation window could be tracked. The interrogation area was 3.6 mm while the largest marker size was 1.3 mm, which means that more than seven marker points fall within an interrogation area. In a preprocessing step, the contrast of the images was increased using histogram equalization, by which detection probability for valid vectors in images could be increased fivefold.^{42} Moreover, it was possible to increase the spatial resolution of the displacement vectors by iteratively changing the size of the interrogation window, whereby an increase of resolution down to 0.04 pixels was reported in Ref. 43. Moreover, a $2\times 3$ Gaussian interpolation was performed on the cross-correlation peak, which reportedly increases the resolution, as discussed in Ref. 44. Combination of both the effects results in a resolution for the displacement vectors of 0.02 pixels. The displacement vector field and its contour map are displayed in Fig. 6.

## 4.2.

### Comparison with FE-Model

Subsequently, the displacement fields and cross-section plots of the FE simulated data were compared with the experimental data obtained for the silicone phantom with a foreign body, see Fig. 7(a). The difference map between simulated and experimentally obtained displacement fields is displayed in Fig. 7(b). A good overall match was obtained with a mean measurement error of 0.0604 mm and a standard deviation of 0.0654 mm, which have both been calculated from the difference map displayed in Fig. 7(b).

High-displacement values are preserved when the pressure wave passes through the foreign body due to their increased stiffness in comparison to the soft silicone, which already indicates the location of the foreign body. However, the displacement does not represent an elastic parameter.

Therefore, the strain ${\u03f5}_{yy}$ along the direction of indentation ($y$-direction) has been calculated, which in addition to representing an elastic parameter, enables the visualization of elastic inhomogeneities with high contrast. In addition to the high positive strain values of the experimental strain map located at the upper edge of the silicone phantom, a similar strain distribution between the FE model and the experimentally obtained strain map was observed; see Fig. 7(c). The mean error for strain difference displayed in Fig. 7(d) was 0.0598 with a standard deviation of 0.0562. A way to overcome the identification ambiguity arising from the high strain values at the edge of the phantom will be discussed in Sec. 4.3 based on strain map averaging.

In conclusion, a good agreement between FE modeling and experiments could be achieved for the displacement and strain along the indentation direction.

Therefore, the FE model represents a realistic model that can be employed to reveal other parameters, which can only be obtained solving the inverse problem. This, in particular, refers to the stress distribution. FE modeling enables the recovery of the strain and stress tensors in 3-D, as shown in Figs. 8(a)–8(d) for ${\sigma}_{xx}$, ${\sigma}_{yy}$, and ${\sigma}_{zz}$. Moreover, the stress tensor at different indenter positions can be recovered. To illustrate the changed indenter location, the 3-D stress distributions ${\sigma}_{zz}$ obtained from FE modelling are displayed in Fig. 8(d) for the indenter location right above the foreign body, and Fig. 8(e) for the indenter location at 10 mm distance from the foreign body.

## 4.3.

### Strain Contrast Enhancement Techniques

Two different techniques to enhance the strain map will be discussed. The first is based on the additional recording of displacements fields from an elastically homogeneous specimen. The second approach is based on the averaging of multiple strain maps, which have become available due to the application of the rolling indenter.

The first approach was implemented subtracting the displacement field obtained from a silicone phantom from the one obtained from the silicone phantom with a foreign body. The position of the foreign body can already be seen in the difference displacement map, see Fig. 9. Finally, a difference strain map $\mathrm{\Delta}{\u03f5}_{yy}$ was calculated from the difference displacement map and an improved contrast could be achieved, as highlighted in the cross-section plot in Fig. 9. For calculating the strain contrast, the smallest maximum and highest minimum have been taken into account and the strain cross-sections normalized for better comparison. The significant values for the normal strain map [shown in Fig. 7(c)] are indicated by the red dots, whereas black dots have been used for indicating the difference strain map, see Fig. 9. A contrast of 0.20 was obtained for the normal strain map, whereas an improvement to a contrast value of 0.64 for the difference strain map could be obtained.

The second approach is based on averaging, which is made possible by the large amount of displacement fields and consequently strain maps produced by rolling indentation. This approach does not only enable the improvement of strain contrast in order to support the identification of the foreign body but also minimizes the detection ambiguity caused by high strain values in close proximity to the indenter position. In Fig. 10, the averaging results obtained from two different amounts of strain maps [61 strain maps in Fig. 10(a) and 247 strain maps in Fig. 10(b)] are shown. It becomes obvious that the positive high strain values at the top edge of the strain map are caused by the indenter, which when moved across the specimen results in an upward directed displacement within its close proximity. These artifacts, which for a single strain map lead to a false interpretation [see Fig. 7(c)], can therefore be associated to the indenter. They can hence be neglected for the identification of the foreign body. To evaluate the detection improvement, the contrast within close proximity of the foreign body has been calculated. For calculating the contrast, the smallest maximum and highest minimum have been taken into account, as indicated by the red (individual) and black (averaged) dots in Fig. 10. For the individual strain map a contrast of 0.20, for averaging 61 strain maps a contrast of 0.92 and for averaging 247 strain maps a contrast of 0.71 could be obtained. Hence, the contrast improvement is restricted to rolling indentation, which is located in close proximity of the foreign body. The averaging approach highlights the advantage of the rolling indentation approach, which enables a strain contrast improvement after the data have been selected. It can furthermore be fine-tuned with respect to the chosen amount of data to result in the best possible contrast.

The implementation in an FEM environment for the difference strain map method should not be too problematic, whereas the averaging method will be more challenging due to the large amount of data that will have to be created.

## 4.4.

### Porcine Kidney

Conventionally, digital image correlation on biological objects is implemented with incoherent light, such as white light, while the surface of the object under investigation is spray coated in order to generate surface-localized markers that can be used to calculate the displacement field.^{29}^{,}^{45} Spray coating is not possible in a surgical environment. To complicate things further, no structural details are available under visible white light conditions on the porcine kidney surface. Commonly, this problem can be overcome via the application of coherent light to produce speckles. However, speckles strongly fluctuate due to the interaction of the coherent light with the biological tissue (absorption, dehydration, and multiple interference from different tissue layers). In fact, the correlation time of speckles is only a few milliseconds.^{46} However, the renal capsule hosts a network of capillaries called glomerulus, at which the blood filtration takes place. Applying a spectrally adjusted illumination, which is sensitive to strong absorption of hemoglobin, can help to visualize these capillaries, so that they can be used for digital image correlation. Therefore, the spectral properties of kidney tissue had been obtained first. For realistic measurements, fresh samples from a nearby slaughter house (Huber Fleischzentrum OHG) have been shock frozen using liquid nitrogen at 77 K in order to preserve the morphology and to prevent the appearance of ice crystals. Afterward, they were slowly warmed up to $-25\xb0\mathrm{C}$ so that the tissue does not split when cutting it. A microtom (Leica CM1900) was employed for cutting the tissue in $60\text{-}\mu \mathrm{m}$-thick sections, which were then embedded between a microscopy slide and a cover slip. Two spectral measurements (spectrometer from MUT GmbH, 450 to 750 nm wavelength range and 0.9 nm spectral resolution), one with and one without the sample, have been taken.

The recording parameters that represent multiplicative factors, such as integration time, recording area, quantum efficiency, and electron-to–digital count conversion, have not been changed during the two measurements. Hence, a conversion from counts into radiant flux $\varphi $ was not necessary to calculate the degree of transmission $T$, which enables the recovery of the absorption coefficient according to Beer–Lambert law^{47}

## (3)

$$T=\frac{{\varphi}_{\text{out}}}{{\varphi}_{\text{in}}}=\frac{{\text{Count}}_{\text{out}}}{{\text{Count}}_{\text{in}}}={e}^{-\alpha \xb7d},$$The approach of tracing the displacement map using the capillaries offers the additional benefit of enabling further tissue discrimination based on the hypervascularization effect in malignant tissue.

## 5.

## Discussion and Conclusion

A 4-D FE-model has been generated taking into account the 3-D spatial coordinates and the changing indenter position over time. The underlying constitutive FE-model applied is the Arruda–Boyce model, which is well suited for the nonlinear elastic behavior of soft tissue since it is hyperelastic. Furthermore, it requires only two elastic input parameters, locking stretch and shear modulus. A first guess of these parameters has been obtained from a uniaxial compression test on well-defined cylindrical silicone samples and the corresponding resulting stress–strain curves, which were then used as input parameters for the FE model. A good match for the 2-D displacement field and strain maps generated from the 3-D FE-model and the experimental 2-D digital image correlation technique using the rolling indenter could be obtained.

It could be demonstrated that a realistic FE model enables the recovery of elastic parameters of the object under investigation at different indenter positions. In this manner, the inverse problem could be solved. It was furthermore demonstrated that the strain contrast could be increased via computing the displacement difference map, based on the recording of the displacement map of a homogeneous sample and via the averaging approach, which furthermore enables the reduction of strain ambiguities that result as a disturbing artifact from the indentation process. In reality, the homogeneous sample required for the first strain contrast improvement method could be a healthy sample or region of the sample. An instrument enabling the application of well-defined experimental parameters has been developed. Future research will be focused on the miniaturization of the setup, to enable its implementation in a minimally invasive tool. A possible solution to implement rolling indentation has been demonstrated in Ref. 49, albeit without the application of digital image correlation.

The issue of missing structural details for the application of digital image correlation for the investigation of kidney could be overcome via a spectrally adjusted illumination source in the green wavelength range.

At this wavelength, a good absorption contrast of hemoglobin in comparison to kidney tissue could be obtained. Further future work will be focused on the extension to kidney sample with a foreign body, its FE modeling implementation and the measurement of 3-D displacement fields.

## Acknowledgments

This project was jointly founded by the state of Baden-Württemberg, Aesculap AG and the Universities of Tübingen and Stuttgart under the scope of the Industry on Campus initiative IoC105. We would like to thank Mr. Schierbaum from the Institute of Applied Physics at the University of Tübingen for carrying out the uni-axial compression tests.

## References

K. Sangpradit et al., “Finite-element modeling of soft tissue rolling indentation,” IEEE Trans. Biomed. Eng. 58, 3319–3327 (2011).IEBEAXIEBEAX0018-9294http://dx.doi.org/10.1109/TBME.2011.2106783Google Scholar

W. A. Berg et al., “Shear-wave elastography improves the specificity of breast US: the BE1 multinational study of 939 masses,” Radiology 262(2), 435–449 (2012).RADLAXRADLAX0033-8419http://dx.doi.org/10.1148/radiol.11110640Google Scholar

W. Derwich et al., “High resolution strain analysis comparing aorta and abdominal aortic aneurysm with real time three dimensional speckle tracking ultrasound,” Eur. J. Vasc. Endovasc. Surg. 51, 187–193 (2016).http://dx.doi.org/10.1016/j.ejvs.2015.07.042Google Scholar

B. F. Kennedy, K. M. Kennedy and D. D. Sampson, “A review of optical coherence elastography: fundamentals, techniques and prospects,” IEEE J. Sel. Top. Quantum Electron. 20(2), 272–288 (2014).IJSQENIJSQEN1077-260Xhttp://dx.doi.org/10.1109/JSTQE.2013.2291445Google Scholar

E. A. Stewart et al., “Magnetic resonance elastography of uterine leiomyomas: a feasibility study,” Fertil. Steril. 95(1), 281–284 (2011).Google Scholar

T. Sasaki, M. Haruta and S. Omata, “CT elastography: a pilot study via a new endoscopic tactile sensor,” Open J. Biophys. 4, 22–28 (2014).http://dx.doi.org/10.4236/ojbiphy.2014.41004Google Scholar

T. Schmidt, J. Tyson and K. Galanuli, “Full-field dynamic displacement and strain measurement using, advanced 3D image correlation photogrammetry,” Exp. Tech. 27(3), 47–50 (2003).EXPTD2EXPTD20732-8818http://dx.doi.org/10.1111/ext.2003.27.issue-3Google Scholar

B. Kemper et al., “Quantitative determination of out-of plane displacements by endoscopic electronic-speckle-pattern interferometry,” Opt. Commun. 194, 75–82 (2001).OPCOB8OPCOB80030-4018http://dx.doi.org/10.1016/S0030-4018(01)01272-XGoogle Scholar

E. Kolenovic et al., “Miniaturized digital holography sensor for distal three-dimensional endoscopy,” Appl. Opt. 42(25), 5167–5172 (2003).APOPAIAPOPAI0003-6935http://dx.doi.org/10.1364/AO.42.005167Google Scholar

P. Hai et al., “Photoacoustic elastography,” Opt. Lett. 41(4), 725–728 (2016).OPLEDPOPLEDP0146-9592http://dx.doi.org/10.1364/OL.41.000725Google Scholar

C. Buj et al., “Speckle-based off-axis holographic detection for photoacoustic tomography,” Curr. Dir. Biomed. Eng. 1, 356–360 (2015).http://dx.doi.org/10.1515/cdbme-2015-0088Google Scholar

K. M. Moerman et al., “Digital image correlation and finite element modelling as a method to determine mechanical properties of human soft tissue in vivo,” J. Biomech. 42, 1150–1153 (2009).JBMCB5JBMCB50021-9290http://dx.doi.org/10.1016/j.jbiomech.2009.02.016Google Scholar

S. J. Kirkpatrick and D. D. Duncan, “Laser speckle strain measurements in soft tissue,” in SEM X Int. Congress and Exposition on Experimental and Applied Mechanics (2004).Google Scholar

G. Antonacci et al., “Quantification of plaque stiffness by Brillouin microscopy in experimental thin cap fibroatheroma,” J. R. Soc. Interface 12, 20150843 (2015).1742-5689http://dx.doi.org/10.1098/rsif.2015.0843Google Scholar

P. So, “Brillouin bioimaging,” Nat. Photonics 2, 13–14 (2008).NPAHBYNPAHBY1749-4885http://dx.doi.org/10.1038/nphoton.2007.260Google Scholar

R. Jacobson, K. Radmacher and M. Matzke, “Direct, high-resolution measurement of furrow stiffening during division of adherent cells,” Nat. Cell Biol. 3, 607–610 (2001).http://dx.doi.org/10.1038/35078583Google Scholar

M. Laidler et al., “Elasticity of normal and cancerous human bladder cells studied by scanning force microscopy,” Eur. Biophys. J. 28, 312–316 (1999).http://dx.doi.org/10.1007/s002490050213Google Scholar

S. E. Cross et al., “Nanomechanical analysis of cells from cancer patients,” Nat. Nanotechnol. 2, 780–783 (2007).NNAABXNNAABX1748-3387http://dx.doi.org/10.1038/nnano.2007.388Google Scholar

O. Chaudhuri et al., “Extracellular matrix stiffness and composition jointly regulate the induction of malignant phenotypes in mammary epithelium,” Nat. Mater. 13, 970–978 (2014).NMAACRNMAACR1476-1122http://dx.doi.org/10.1038/nmat4009Google Scholar

J. R. Lange and B. Fabry, “Cell and tissue mechanics in cell migration,” Exp. Cell Res. 319(16), 2418–2423 (2013).http://10.1016/j.yexcr.2013.04.023Google Scholar

M. M. Doyley, “Model-based elastography: a survey of approaches to the inverse elasticity problem,” Phys. Med. Biol. 57, R35–R73 (2012).http://10.1088/0031-9155/57/3/R35Google Scholar

A. Manduca et al., “Magnetic resonance elastography: non-invasive mapping of tissue elasticity,” Med. Image Anal. 5(4), 237–254 (2001).http://dx.doi.org/10.1016/S1361-8415(00)00039-6Google Scholar

K. C. Maitland and T. D. Wang, “Endoscopic image resolution,” in Biomedical Technology and Devices, 2nd ed., p. 226, Taylor and Francis Group, Boca Raton, Florida (2010).Google Scholar

J. Didier and C. Roy, “Angiography as a diagnostic tool in renal cell carcinoma,” in Renal Cell Cancer Diagnosis and Therapy, pp. 127–130, Springer, London (2008).Google Scholar

C. E. Willert and M. Gharib, “Digital particle image velocimetry,” Exp. Fluids 10(4), 181–193 (1991).EXFLDUEXFLDU0723-4864http://dx.doi.org/10.1007/BF00190388Google Scholar

B. L. Smith and D. R. Douglas, “Particle image velocimetry,” in Handbook of Fluid Dynamics, pp. 1–48, CRC Press, Boca Raton, Florida (2016).Google Scholar

P. J. Bryanston-Cross and A. Epstien, “The application of sub-micron particle visualisation for PIV (Particle Image Velocimetry) at transonic and supersonic speeds,” Prog. Aerosp. Sci. 27, 237–265 (1990).http://dx.doi.org/10.1016/0376-0421(90)90008-8Google Scholar

D. Zhang and D. D. Arola, “Application of digital image correlation to biological tissues,” J. Biomed. Opt. 9(4), 691–699 (2004).http://dx.doi.org/10.1117/1.1753270Google Scholar

T. Luyckx et al., “Digital image correlation as a tool for three-dimensional strain analysis in human tendon tissue,” J. Exp. Orthop. 1(7), 1–9 (2014).http://10.1186/s40634-014-0007-8Google Scholar

R. D. Keane and R. J. Adrian, “Theory of cross-correlation analysis of PIV images,” Appl. Sci. Res. 49, 191–215 (1992).ASRHAUASRHAU0003-6994http://dx.doi.org/10.1007/BF00384623Google Scholar

D. Zhang, X. Zhang and G. Cheng, “Compression strain measurement by digital speckle correlation,” Exp. Mech. 39(1), 62–65 (1999).EXMCAZEXMCAZ0014-4851http://dx.doi.org/10.1007/BF02329302Google Scholar

T. Kreis, Handbook of Holographic Interferometry: Optical and Digital Methods, Springer, Weinheim (2005).Google Scholar

W. Voigt, Lehrbuch der Kristallphysik, 11928th ed., Teubner, Stuttgart (1966).Google Scholar

K. M. Kennedy et al., “Analysis of mechanical contrast in optical coherence elastography,” J. Biomed. Opt. 18(12), 121508 (2013).http://dx.doi.org/10.1117/1.JBO.18.12.121508Google Scholar

M. Farshad et al., “Material characterization of the pig kidney in relation with the biomechanical analysis of renal trauma,” J. Biomech. 32, 417–425 (1999).JBMCB5JBMCB50021-9290http://dx.doi.org/10.1016/S0021-9290(98)00180-8Google Scholar

L. Dong et al., “Quantitative compression optical coherence elastography as an inverse elasticity problem,” IEEE J. Sel. Top. Quantum Electron. 22(3), 277–287 (2016).IJSQENIJSQEN1077-260Xhttp://dx.doi.org/10.1109/JSTQE.2015.2512597Google Scholar

A. Wittek et al., “A finite element updating approach for identification of the anisotropic hyperelastic properties of normal and diseased aortic walls from 4D ultrasound strain imaging,” J. Mech. Behav. Biomed. Mater. 58, 122–138 (2015).http://dx.doi.org/10.1016/j.jmbbm.2015.09.022Google Scholar

E. M. Arruda and M. C. Boyce, “A three-dimensional constitutive model for the large stretch behaviour of rubber elastic materials,” J. Mech. Phys. Solids 41, 389–412 (1993).JMPSA8JMPSA80022-5096http://dx.doi.org/10.1016/0022-5096(93)90013-6Google Scholar

P. Wijesinghe, B. F. Kennedy and D. D. Sampson, “Localized loading to reduce artifacts in compression optical coherence elastography,” in SPIE Biomedical Optics Conf. (BiOS), San Francisco, pp. 9327–9325 (2015).Google Scholar

F. Sozio et al., “Morphometric analysis of intralobular, interlobular and pleural lymphatics in normal human lung,” J. Anat. 220(4), 396–404 (2012).JOANAYJOANAY0021-8782http://dx.doi.org/10.1111/joa.2012.220.issue-4Google Scholar

W. Thielicke and E. J. Stamhuis, “PIVlab—towards user-friendly, affordable and accurate digital particle image velocimetry in MATLAB,” J. Open Res. Softw. 2(1) (2014).http://doi.org/10.5334/jors.blGoogle Scholar

U. Shavit, R. J. Lowe and J. V. Steinbuck, “Intensity capping: a simple method to improve cross-correlation,” Exp. Fluids 42, 225–240 (2007).http://dx.doi.org/10.1007/s00348-006-0233-7Google Scholar

F. Scarano and M. L. Riethmuller, “Iterative multigrid approach in PIV image processing with discrete window offset,” Exp. Fluids 26, 513–523 (1999).EXFLDUEXFLDU0723-4864http://dx.doi.org/10.1007/s003480050318Google Scholar

T. Roesgen, “Optimal subpixel interpolation in particle image velocimetry,” Exp. Fluids 35, 252–256 (2003).http://10.1007/s00348-003-0627-8Google Scholar

D. Zhang and D. D. Arola, “Applications of digital image correlation to biological tissue,” J. Biomed. Opt. 9(4), 691–699 (2004).http://dx.doi.org/10.1117/1.1753270Google Scholar

X. L. Dean-Ben et al., “Wavefront shaping based on three-dimensional optoacoustic feedback,” in European Conf. on Biomedical Optics, Munich (2015).Google Scholar

H. Naumann and G. Schröder, Bauelementer der Optik, 6th ed., Carl Hanser, München (1992).Google Scholar

S. Prahl, “Optical absorption of hemoglobin,” http://omlc.org/spectra/hemoglobin/ (15 December 1999).Google Scholar

P. Puangmali et al., “Miniature 3-axis distal force sensor for minimally invasive surgical palpation,” IEEE/ASME Trans. Mechatron. 17(4), 646–656 (2012).IATEFWIATEFW1083-4435http://dx.doi.org/10.1109/TMECH.2011.2116033Google Scholar

## Biography

**Daniel Claus** received his MScEng from Technical University of Ilmenau in 2006 and his PhD from University of Warwick, Great Britain, in 2010. He joined the Institut fuer Technische Optik at the University of Stuttgart in 2013. His research areas include digital holography, ptychography, phase retrieval, light field imaging, shape measurement, optical testing, optical elastography, biomedical imaging, and opto-mechanical design.

**Marijo Mlikota** received his master’s degree in mechanical engineering from the University of Zagreb in 2010. Currently, he is working on a doctoral degree at the University of Stuttgart in the field of metal fatigue and water droplet erosion modelling. Other research areas include modelling of soft tissue behavior and fracture behavior of polymer materials.

**Jonathan Geibel** received his BS degree in mechanical engineering in 2015. After that he continued pursuing his master’s degree in mechanical engineering in Stuttgart. He worked at the Institut für Technische Optik, University of Stuttgart during his bachelor's thesis from September 2015 to April 2016.

**Giancarlo Pedrini** received his MS degree in physics from the Swiss Federal Institute of Technology, ETH-Zurich, in 1982 and his PhD in optical sciences from the University of Neuchatel, Switzerland, in 1990. He joined the Institut fuer Technische Optik at the University of Stuttgart in 1991. His research areas include digital holography, vibration analysis, shape measurement, optical testing, measurement of the elastic parameters of biological samples, endoscopy, and phase retrieval.

**Siegfried Schmauder** received his diploma in mathematics in 1981 and his PhD in materials science in 1988. Since 1994, he has been a professor in strength of materials and materials techniques in mechanical engineering at the University of Stuttgart. His research is focused on micromechanical analyses and multiscale materials modelling (MMM) to predict mechanical and functional properties of different materials. In these fields he published several books and about 500 articles in scientific journals.

**Wolfgang Osten** received his MSc in physics from the Friedrich-Schiller-University Jena in 1979 and his PhD from the Martin-Luther-University Halle-Wittenberg in 1983. Since 2002, he has been a full professor at the University of Stuttgart. His research work focuses on new concepts for industrial inspection, combining modern principles of optical metrology, sensor technology and image processing. Special attention is directed to the development of resolution-enhanced technologies for the investigation of micro and nanostructures.