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,²⁵ in order to study fluid flow²⁶ and air flow.²⁷ 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.

Fig. 1

Schematic representation of the working principle for digital image correlation.

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

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{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.

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.²¹ 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.

Fig. 2

Parameter-identification based on the comparison between the outcome of the FE-model (solution of direct problem) and experimentally measured values (displacement field and strain map).

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³⁵) 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,¹ a 3-D axisymmetric FE model of linear elastic behavior,³⁴ and a four-dimensional (4-D) hyperelastic model taking into account the temporal component.³⁷ 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.³⁸ 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.

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.

Fig. 3

Uniaxial compression test results: (a) force indentation diagram for porcine kidney highlighting its hyperelastic and viscoelastic behavior, (b) force-indentation diagram obtained from cylindrical silicone samples of hardsilicone and softsilicone (dimensions in mm), (c) corresponding stress–strain diagram, and (d) dimensions of silicone phantom.

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{\mathrm{mm}}^3$ to $5\times 5\times 5{\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).

Fig. 4

3-D FE model showing the displacement along $y$-direction of silicone phantom and indenter, sliced at central location of foreign body (parameters applied 3 mm indentation, indenter located on top of foreign body).

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\mathrm{kPa}$ and ${\lambda }_m=7.7$) and hard silicone ($G=7.73\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.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).

Fig. 5

(a) Setup: roll-indenter with silicone phantoms, (b) graphical representation of the force applied along the indenter’s rolling direction with and without foreign body.

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\mu \mathrm{m}$). The rolling indenter can be displaced laterally by 156 mm with a positioning uncertainty of $50\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.³⁹ 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\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\mathrm{mm}/\mathrm{s}$ across the specimen. A digital camera (Zyla 5.5, $2560\times 2160\text{pixels}$, $6.5\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\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\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\mu \mathrm{m}$⁴⁰ can well be resolved and imaged. The field of view recorded by the camera was $144\mathrm{mm}\times 121\mathrm{mm}$.

In order to save memory and more importantly postprocessing time, the images have been cropped to $1800\times 520\text{pixels}$, which corresponds to a field of view of $101\mathrm{mm}\times 29\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\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\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.⁴¹ An interrogation area of $65\times 65\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.⁴² 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.

Fig. 6

Experimental results: displacement maps obtained: (a) displacement vector field and (b) contour map representation (absolute value of displacement vector).

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).

Fig. 7

Comparison between FE simulation and experimentally obtained results (top FE Simulation, center experimentally obtained displacement field, cross-section plot, highlighted by dark line in experimental data) (a) displacement field, (b) strain map, (c) difference map between simulated and measured displacement field, and (d) difference map for simulated and calculated strain.

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 ${ϵ}_{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.

Fig. 8

Flow chart and results for obtained displacements fields with and without foreign body, the resulting difference displacement field, the calculate strain field, and corresponding cross-section plot, which compares the difference approach strain field with the conventional results shown in Fig. 7(b).

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 }{ϵ}_{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.

Fig. 9

Stress distributions obtained using the hyperelastic Arruda–Boyce model for Abaqus (a) ${\sigma }_{xx}$ at the front surface, (b, c, d) 3-D stress distribution along (b) the $x$-direction ${\sigma }_{xx}$, (c) $y$-direction ${\sigma }_{yy}$, (d) and $z$-direction ${\sigma }_{zz}$ with indenter positioned on top of inhomogeneity, (e) 3-D stress distribution along $z$-direction ${\sigma }_{zz}$ with indenter positioned in 10 mm distance to inhomogeneity.

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.

Fig. 10

Averaging the strain maps obtained from: (a) 61 recorded indenter position, (b) 247 different indenter positions, which demonstrates that the high strain values at top edge of the strain map arise from the indenter. The cross-section plot at the bottom of each subfigure highlights the contrast improvement for identifying the foreign body in comparison to the individual strain map obtained from the centrally position indenter shown in Fig. 7(c).

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.⁴⁶ 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°\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⁴⁷

Fig. 11

(a) Absorption coefficient hemoglobin and kidney, (b) section of kidney illuminated with LED 505-nm peak wavelength, and (c) same section of kidney illuminated with white light.

Fig. 12

(a) Reference shape measurement using a constant force of 0.1 N, (b) displacement field obtained via digital image correlation from a porcine kidney, and (c) strain map for kidney.

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.