Highly resolved strain imaging during needle insertion_ Results with a novel biologically inspired device

Percutaneous needle insertions are a common part of minimally invasive surgery. However, the insertion process is necessarily disruptive to the substrate. Negative side effects are migration of deep-seated targets and trauma to the surrounding material. Mitigation of these effects is highly desirable, but relies on a detailed understanding of the needle–tissue interactions, which are difficult to capture at a sufficiently high resolution. Here, an adapted Digital Image Correlation (DIC) technique is used to quantify mechanical behaviour at the sliding interface, with resolution of measurement points which is better than 0.5mm, representing a marked improvement over the state of the art. A method for converting the Eulerian description of DIC output to Lagrangian displacements and strains is presented and the method is validated during the simple insertion of a symmetrical needle into a gelatine tissue phantom. The needle is comprised of four axially interlocked quadrants, each with a bevel tip. Tests are performed where the segments are inserted into the phantom simultaneously, or in a cyclic sequence taking inspiration from the unique insertion strategy associated to the ovipositor of certain wasps. Data from around the needle–tissue interface includes local strain variations, material dragged along the needle surface and relaxation of the phantom, which show that the cyclic actuation of individual needle segments is potentially able to mitigate tissue strain and could be used to reduce target migration. & 2013 The Authors. Published by Elsevier Ltd. All rights reserved. thors. Published by Elsevier Ltd. All rights reserved. under the terms of the Creative Commons Attribution License, which permits unrestricted edium, provided the original author and source are credited. (F. Rodriguez y Baena).


Introduction
Needle insertion is common practice in clinical intervention, with applications ranging from cyst drainage and tissue biopsy, to brachitherapy and drug delivery. However, while needles generally are minimally intrusive, soft tissue damage and target migration can still occur, which may cause complications for both the patient and the clinician. With a view to understanding the underlying physics of the insertion process and designing better and more effective tools, a number of studies have focussed on measuring and modelling tool-tissue interaction behaviour during needle insertion. In many instances, the inter-actions have been parameterised e.g. Majewicz et al. (2012), Misra et al. (2010), van Veen et al. (2012), or modelled numerically e.g. Abolhassani and Patel (2006), Dong et al. (2012), Glozman and Shoham (2007), and Kataoka et al. (2001). In previous work we have tried to tackle this challenge by developing high resolution finite element models that incorporate material deformation, material failure and contact mechanics (Oldfield et al., 2013a). Detailed numerical modelling of the interactions has proved challenging due to the number and complexity of competing phenomena at the needle-tissue interface.
To improve, tune and validate high resolution models that capture all aspects of needle-tissue interactions, experimental data examining behaviour around needle tips is required at a similar resolution. However, this data and the techniques for obtaining it are not currently available. In many instances, the consistency required of the experimental substrate has dictated that a soft tissue phantom is used. Despite the greater control that is offered by a phantom, to date, the deformation behaviour has only been captured on the surface, at low resolution or without revealing the deformed state of the substrate except at its boundaries. These limitations are addressed here.
We have recently demonstrated a laser-based technique for highly resolved planar deformation measurements in a tissue phantom (Kerl et al., 2012). The technique is based on Digital Image Correlation (DIC) and is applied in Fluid Dynamics to extract an Eulerian description of velocities in a spatially fixed grid, a technique known as Particle Image Velocimetry (PIV), see e.g. Raffel et al. (1998). In PIV, particles are suspended in the fluid phase and tracked over time. However, in Kerl et al. (2012), we embedded micrometre scale particles at random but spatially fixed locations in a tissue phantom to extract time resolved displacements (TRDs) and strains in an Eulerian description of behaviour. At present, this remains the only technique that obtains images at such a scale and resolution, is taken in a plane running through the contact interface and also reveals the internal deformations of the soft tissue phantom. However, in solid mechanics, a Lagrangian description of material behaviour is commonly used, for example when tracking speckle (O'Donnell et al., 1994) or surface definitions (Buckley et al., 2010(Buckley et al., , 2008Michalek et al., 2009). Consequently, a method is developed here, based on the internal, Eulerian description DIC measurements in Kerl et al. (2012), that provides output measurements in a Lagrangian description with an unprecedented level of detail at the tool-tissue interface. This method will subsequently be referred to as the Image Correlation Method (ICM).
Here, we use the ICM to evaluate measurements performed on a multi-segment needle which is being developed for soft tissue surgery (Frasson et al., 2011;Ko and Rodriguez y Baena, 2012) and is based on biomimetic concepts. As previously hypothesised (Parittotokkaporn et al., 2010;Aoyagi et al., 2012), we demonstrate that when each segment is actuated in turn (reciprocal motion), both the deformation and strain imparted on a gelatine phantom during insertion are reduced when compared to an equivalent, simultaneous push of the segments together. The discrepancy between mechanical and functional damage of brain tissue is highlighted by Hrapko et al. (2006) and functional damage due to strain is shown in the context of traumatic brain injury for white matter (Wright and Ramesh, 2012) and grey matter (Elkin and Morrison, 2007). As reduced strain at the tool-tissue boundary could be linked to reduced trauma, reciprocal insertion could thus be of benefit when interacting with soft and delicate tissue such as the brain, and particularly in applications where the needle is expected to deviate from a straight path to avoid obstacles e.g. Engh et al. (2010).
In this paper, illustration of the enhanced capability of the ICM is provided by a comparison of the strains developed during insertion of a four-part needle with all segments actuated simultaneously and in reciprocating motion. Following is a presentation of the theoretical construct upon which the Lagrangian output data is based. Experimental methodsincluding actuation strategiesand critical analytical parameters are presented in Section 3. A presentation and discussion of results, with implications as to the efficacy of the techniques described herein and the competing actuation mechanisms of a four part needle, are presented in Section 4. Conclusions and areas for further investigation are presented in Section 5.

Theory
In the following section, the theoretical basis for converting the Eulerian description of TRDs to a Lagrangian description of material deformation and strain is presented. Initially, the Lagrangian representation of displacement is defined. Subsequently, the interpolative steps required to convert the Eulerian data obtained from the DIC analysis to a Lagrangian description tracking material points are shown. This forms the basis of the ICM, and the procedure for obtaining material strains from the Lagrangian displacement data is finally given.

Undeformed and deformed states
The gelatine phantoms employed in this study were assumed to have an initially undeformed state. Strains associated with body forces were considered negligible compared to the strains resulting from needle insertion. In this undeformed state, the gelatine was defined using the material coordinates, X. During subsequent time steps, the gelatine was defined using the spatial coordinates, x, with both coordinate frames coincident. Points in the undeformed state X at time t¼ 0 are mapped into the deformed state at time, t: In the two-dimensional analysis performed here, x ¼ x y the undeformed and deformed positions are defined in a single Cartesian coordinate frame, in the plane of the applied laser light sheet, where X is aligned with the insertion axis. The material displacements, u, are defined in the same Cartesian frame as follows: 2.2. Lagrangian vs. Eulerian measurements By using DIC it is possible to extract TRDs at a fixed, stationary grid of measurement points for the gelatine phantom. In the Lagrangian description the temporal motion and mechanical properties of material points are tracked as they move through space. Data from the experiments described in Section 3.1 was obtained in an Eulerian description using DIC (Kerl et al., 2012) and, here, the conversion into the Lagrangian description is outlined. At time t ¼0 the gelatine was assumed to be at rest and in an unstrained state. The DIC measurement grid locations provided the material points that were subsequently traced during the tests that were performed. Rather than specific particles, these points in the material, tracked throughout the analysis, were only coincident with the fixed Eulerian measurement grid at the initial time step. To calculate the displacement of a material point from one time step to another, the Euler method was used: where n is the image frame number, _ x n represents the material point's TRD and Δt is the time increment between image frames n and n þ 1. To obtain _ x n at the material point itself, spline interpolation was performed from the TRD vector field at the fixed, Eulerian measurement grid. This was implemented through the interp2 function in MATLAB s at position x n .

Elemental strain calculation
Deformations observed in the gelatine were both finite and inhomogeneous. The finite and volumetric element methods offer a suitable framework for dealing with this type of response. The initial mesh of material can be constructed from the DIC measurement grid with subsequent deformations during needle insertions. Stresses are not required in the context of this paper and strains can be calculated directly from the deformed structure. The structure is based on triangular elements with the nodal points, or vertices, defined by the initial measurement points and their migrated positions (Fig. 1). Constant deformation and hence strain was assumed within the elements with its evaluation taking place at the centroid of triangular element e, X e c or x e c , defined by its three vertices, X e 1 ; X e 2 ; X e 3 or x e 1 ; x e 2 ; x e 3 respectively X e c ¼ 1 x Using the method of Teran et al. (2003), the lengths of two pairs of equivalent sides belonging to the material element in undeformed and deformed states, e.g. S 12 and s 12 respectively, lead to the deformation gradient: S e ¼ ½S e 12 S e 13 ð 7aÞ The deformation gradient for an element of material, F e , can then be calculated as follows: From the deformation gradient several different measures of finite strain can be evaluated. Here, the Green-Lagrange strain in the material element, G e , is used: where I is the identity matrix.

Experimental configuration
Needle insertions were performed in a tissue phantom to investigate the effectiveness of the ICM and different insertion mechanisms. The insertions were performed in a phantom using 6% by mass of gelatine powder (Kerl et al., 2012) and set in a 84 mm Â 58 mm Â 24.5 mm volume. The needle insertion axis ran parallel to the longest dimension and was at a height of 23 mm from the bottom. Aluminium Oxide (Al 2 O 3 ) powder was combined with the gelatine in its liquid state at a concentration of 1 g per litre to enable DIC. Images were collected using a CCD camera (LaVision Imager Intense, 1376 Â 1040 pixels, 290-1100 nm spectral range) in a single plane, which was aligned with the insertion axis of the four-part needle. A seven dot calibration target (LaVision) enabled pixel/mm conversion of camera images and corrected for the insignificant image distortion present. Projection errors were small due to the normal orientation of the camera to the laser light sheet and the single, transparent wall of the gelatine container. Any residual projection errors were eliminated by use of a dewarping algorithm. A field of view (FOV) of 39.4 mm Â 27.6 mm was achieved with a plane of laser light 100 μm deep. DIC was performed using LaVision Davis software (v7.2). Vectors of TRD were obtained on a 90 Â 63 grid with 0.439 mm horizonal and vertical spacing between measurement locations. These vectors had a 10 point median filter retrospectively applied to them. Due to the slow rate of Fig. 1 -Element, e, in undeformed material coordinates, X, and subsequently, deformed, in spatial coordinates, x. A description of the DIC process used to extract TRDs from the output images is provided in essence in Kerl et al. (2012) and in considerable depth in e.g. Raffel et al. (1998). Sub-regions of interest were compared between successive images with peak cross-correlation between these subregions yielding a displacement. Division by the image sampling rate provided TRDs at a fixed spatial grid, with each point centering the sub-regions of interest. In performing DIC a multi-pass operation was used with correlation windows decreasing in size from 128 Â 128 pixels to, finally, 32 Â 32 pixels (Kerl et al., 2012). A 50% overlap of the windows was specified which generated fixed grid positions for output with 16 pixel spacing.
A 6 mm diameter needle, made of four axially interlocked segments, was inserted into the gelatine block. Each segment tip was a single bevel and, when aligned, formed a prism. It was assumed that the needle segments were relatively rigid in comparison to the gelatine block. The four interlocking needle segments were manufactured from the proprietary material FullCure720 (Objet Geometries Ltd.). Each segment was actuated and controlled independently using DC motors (Maxon Amax22, Maxon Motors Ltd.) with built-in encoders. These motors were controlled using a CompactRio (National Instruments). Motion was transferred from the motors to the needle segments via individual lead screws and linear bearings. Insertion into the gelatine was made more consistent by a rigid trocar between the needle base and phantom block (Kerl et al., 2012).
Three different insertion mechanisms were used to compare the impact on the surrounding phantom. One approach actuated the needle segments independently in a 'reciprocal motion' insertion strategy. During reciprocating motion, working in a clockwise sense, the first, third, second and fourth segments were moved sequentially with a stroke length of 5 mm. There was no rest period between the end of the stroke of one needle segment and the start of the next segment in the cycle. Two different segment velocities, 4 mm/s (Rm Â 4) and 1 mm/s (Rm Â 1), were used for reciprocating motion. In an alternative mechanism, the four interlocking needle segments were simultaneously inserted (direct push) at a rate of 1 mm/s (Dp Â 1). For each insertion strategy, the relative progression of the four interlocked segments is shown in Fig. 3. Each insertion mechanism was tested three times, with individual tests identified by a subscript where necessary.

Image segmentation parameters
In order to segment the needle from the gelatine phantom, a series of image processing steps were implemented, as shown in Fig. 4. To create stronger edges (contrast) between the translucent needle material and the background before smoothing, a sigmoid function (Tanaka et al., 2007) was used initially. Correctly configured, this non-linearly mapped the needle edge pixels to higher intensities while the background pixel values were mapped to approximately the same intensity. A median filter (Brownrigg, 1984), 10 pixel square neighbourhood, removed small sharp peaks associated with particles of the order of 2 pixels in diameter. These peaks were caused by the particles suspended within the phantom and over saturation of the camera caused by strong scattering from the needle surface. A bilateral filter (Tomasi and Manduchi) then smoothed the image further whilst maintaining the strong edge between the needle and phantom.
To segment the smoothed image, a watershed algorithm (Roerdink and Meijster, 2000) with seed points was used  to divide the image into two regions. The seed points within the image were chosen by using the known geometry and position of the needle. A flood fill algorithm was, finally, performed to produce a binary mask (Fig. 4  (d)). This mask was then used to eliminate those TRD vectors associated with the needle features from the interpolation of TRD vectors at evolving material points.
Masking was also performed on material points that had moved beyond the PIV measurement grid.

Material tracking
During the insertion experiments, the maximum velocity of any needle segment was 4 mm/s. It is assumed that no particles in the gelatine block moved with a velocity greater than this. At a 10 Hz image sampling rate the largest conceivable displacement of a particle fell within the 50% overlap range of the smallest 32 Â 32 pixel correlation window. The image sampling was, therefore, deemed to be performed at a sufficiently high rate. A separate check was provided by a quality factor which related the highest correlation peak with its nearest alternative. Only 1% of TRDs failed to match the acceptability threshold and these were consequently rejected.
To evaluate strain in the Lagrangian frame, the initial and current position of the tracked material points must be known. For two representative tests (Dp Â 1 3 and Rm Â 1 3 ), the current position and path of material points are shown in Figs. 5 and 6. Even in these plots, showing only one in every five of the horizontal and vertical material points, a significant difference between the two actuation mechanisms is identifiable. For both mechanisms the material is displaced around the profile of the needle as it is inserted, firstly with increasing u in advance of the needle tip and then with an increasing magnitude of v as the needle nears. Trajectories for the direct push insertion are relatively smooth and display a sustained period of relaxation after the needle comes to rest (Fig. 5(b)). The phantom relaxation is characterised by a predominantly horizontal motion of decreasing u. During reciprocating motion, Fig. 6(a) and (b) demonstrates that the gelatine phantom will incrementally relax around the needle segments currently at rest. There is, however, no sustained, observable relaxation period once the overall needle motion has been arrested. It is also evident that different insertion mechanisms are most pronounced close to the insertion axis.

Detailed comparison of material motion
Figs. 7-9 show the displacements, u and v, once the needle tip has penetrated close to the (0,0) position. Displacements, u appear greater in the direct push test, particularly in the locality of the needle. Both the peak magnitudes and extent of large displacements are greater in the direct push insertions ( Fig. 7(a)). However, this effect varies during each cycle of reciprocal motion. Unsurprisingly, the vertical displacement magnitudes have a peak value similar to the needle radius in both direct push and reciprocating motion insertions. Displacements aligned with the insertion axis, u, demonstrate a more marked difference between the two insertion mechanisms.

Strain comparisons
Deflection of the gelatine tissue phantom has implications when considering migration of the target of a needle insertion. However, equally important is the strain generated during the interaction between the needle and its substrate, as strain is an indicator of the relative level of damage to surrounding tissue in clinical applications (Elkin and Morrison, 2007). Figs. 10-12 illustrate the direct x (G e 11 ), y (G e 22 ) and shear (G e 12 ) Green-Lagrange strains for tests Dp Â 1 3 , RmÂ 4 1 and Rm Â 1 3 respectively.
One feature of the strain plots is the relative lack of smoothness in comparison with the displacement plots. This is due to the local derivatives of the displacement data being used to generate strain output. It should be noted that this is a particular problem in the Rm Â 1 1,2,3 tests  There are areas of compressive x-strain ahead of the needle tip becoming a tensile strain when moving from the needle tip to the base of the shaft (Figs. 10(a), 11(a) 12(a)). Conversely, all tests (Figs. 10(b), 11(b) and 12(b)) show tensile y-strains around the needle tip as the phantom material is forced apart. This tensile strain becomes compressive as the separated material follows the widening profile of the needle and fills the reduced space between its edge and the walls of the gelatine container. In direct push, shearing (Fig. 10(c)) is larger than for reciprocating motion at equivalent points around the needle (Figs. 11(c) and 12(c))though, naturally, this effect varies in magnitude depending on the point in the reciprocating motion cycle when the measurements were made.
The relative merits of different actuation strategies are difficult to quantify from individual frames. Instead a metric that seeks to achieve this is proposed whereby the mean plus one standard deviation of the magnitudes of u and v in each frame are calculated (Figs. 13(a) and (b)). Although displacements are a less ideal metric than strain, the relative insensitivity to small errors prevents outliers and nonsmooth data from excessively distorting the results. Three features can be identified from these plots. Firstly, the response due to direct push insertion generally provides an upper limit to the behaviour of the reciprocating motion (Table 1). Secondly, the final relaxation component is more prominent in the axial motion. Thirdly, the final, mean u for all actuation mechanisms following phantom relaxation, are similar. The implications, combined with strain output in Figs. 10-12, are that reciprocating motion has the potential to be both a more stable insertion mechanism and reduce tissue strain. This has significant potential for reducing target migration. Deliberate calibration of the reciprocal actuation mechanism may also further reduce the amount of deformation induced by the insertion process. However, with three insertions by each mechanism, and the difficulty in extracting a reliable global strain metric, these findings are not presented as statistically significant. Gelatine can be inconsistent in its failure with planar cracks visibly forming in some instances and not in others. Experiments were performed sequentially on the same day under consistent environmental conditions to maintain a legitimate comparison between them. However, the small number of tests makes identifying mean behaviour and outliers difficult. Only a qualitative assessment has therefore been undertaken. In this assessment, the difference between direct push and reciprocal motion is more apparent than that between the two speeds of reciprocating motion. Subsequently, future experiments will have more tests to allow representative means of e.g. peak strains for similar actuation mechanisms and standard deviations indicating the consistency of the tissue phantom itself.
Individual Al 2 O 3 particles were not all present in the field of view throughout each test. One source of particle 'dropout' from frames is movement of the phantom material across the single laser light plane. This type of movement is assumed to be small but is reliant on careful alignment of the insertion axis with the laser light sheet. However, it has been shown that insertions into gelatine (Oldfield et al., 2013a), and more   generally (Shergold and Fleck, 2005), can result in planar cracks forming. Orientation of cracks relative to the laser light sheet may also result in asymmetric motion of the tissue phantom material across the laser light plane. A final source of lateral motion across the laser light plane is present in the reciprocal motion of the needlewhich is by definition asymmetric.

5.
Conclusions and future work    (2010) successfully quantify high resolution strain images of tissue but the approach is currently only applicable to arbitrarily thin samples capable of transmitting a polarised light source. Here, a method is presented for using a DIC technique to analyse needle-tissue interactions. The ICM effectively converts an Eulerian description of behaviour obtained using DIC to a Langrangian description. This enables displacement and strain data around the needle-phantom substrate interface to be obtained at a resolution at least an order of magnitude superior than the current state of the artfiducial marker-based strain mapping. Through a superimposed mesh, Green-Lagrange strains are generated. While the displacement output is relatively smooth, strain measurements are more susceptible to experimental uncertainty. Some information is also lost as material points move beyond the limits of the image FOV or are covered by the advancing needle. The output remains, though, ideal for optimising needle geometries and for local validation of contact and cutting models.
Tests on three needle actuation strategies reveal a significant difference in the behaviour of the tissue phantom. During direct push experiments, relaxation of the phantom only occurs when all segments come to rest. Conversely, reciprocal motion allows an incremental relaxation to take place around the stationary segments while another is still moving. Direct push experiments generally provide an upper boundary on the phantom displacement prior to relaxation. Reciprocal motion is thus demonstrated to be beneficial as a potentially tissue sparing process and in reducing target migration during insertion.
Output data may be used in the calibration of numerical models. In Oldfield et al. (2013a) for example, cohesive elements were tuned using fracture toughness to match a global force-displacement response. Here, output data is sufficiently refined to enable the fine-tuning of cohesive finite elements so that the global force-displacement and local displacement field around the needle tip could be matched during ongoing material failure. Concentrations of forces around the needle tip of the type shown in DiMaio and Salcudean (2002) may also be confirmed by the displacement and strain plots (Figs. 7-12) without extrapolation from low resolution or surface measurements.
Internal measurements in the gelatine tissue phantom are also of benefit in the calibration of contact interactions. Rather than matching gross behaviour (Oldfield et al., 2013a) or resistance at a particular location in numerical models through force measurements (Kataoka et al., 2002(Kataoka et al., , 2001, a more comprehensive validation may be achieved by comparing local displacements, stiction and drag shown at the needle-phantom interface. High resolution, temporally-evolving, local measurements are beneficial for this purpose (Asadian et al., 2011(Asadian et al., , 2012 and demonstrated here.
For predictive numerical models, rate dependency, cutting mechanics and contact forces must be correctly identified. The adaptation of DIC techniques has enabled displacement and strain to be examined around the needle tip at a high resolution. Additionally, the method is robust to disappearing Al 2 O 3 particles over the course of needle insertion, as displacements are inferred from TRD data. The Al 2 O 3 particles used here to seed the gelatine phantom are small enough to remain independent of the needle and gelatine's response. A typical resolution of particles in the undeformed state was 5635 particles in a FOV of 1091 mm 2 giving a resolution of more than 5 particles/mm 2 and these values could be improved with ease by changing particle concentration and optical configuration.
Irrespective of the quality and usefulness of the results achieved, the approach taken results in some uncertainties and errors. One particular benefit of the methods presented here is the ability to track interactions very close to the needle-tissue interface. This is also the region most prone to errors in TRD vectors due to confusion between the edge of the needle and Al 2 O 3 particles. As such, negative values of the determinant of the deformation gradient (det F) were used to check for unfeasible negative volumes. Occurrences were limited to a negligible number of material elements close to the needle-tissue interface. This suggests that masking of the needle could be developed further. However, high strain gradients were still well rendered around the needle tip despite some material point traces being covered by the masked needle. Some information was also lost due to masking at the right hand edge of measurement grid. Points passing beyond this limit were unable to show subsequent phantom relaxationthough sufficient data remained to understand the process. Subsequent investigations will focus on improving the post processing of TRD output. Two-dimensional temporal and spatial smoothing will be considered and, in the case of strain, techniques like those presented in Gao and Desai (2010); Meng et al. (2007) may prove beneficial. To improve image quality and ease needle segmentation, phosphorescent particles instead of purely scattering Al 2 O 3 offers significant potential (Fond et al., 2012). After excitation with ultraviolet light, these materials show phosphorescence emissions in the visible region which can be detected using conventional cameras. Ultraviolet light scattered from the needle or from the cuvette can be blocked using an ultraviolet filter while the particles themselves remain visible. Experiments that are planar in nature will also be implemented to remove uncertainty associated with material moving across the insertion axis. Measurements normal to the insertion axis, in conjunction with accurate force extraction e.g. Oldfield et al. (2013b), will also reveal the nature of the fracture mechanics that take place within gelatine tissue phantoms. Tests to establish axisymmetry and other threedimensional effects of the needle-tissue interactions will be performed by using concurrent measurements, along the insertion axis, in a dual plane arrangement (Pfadler et al., 2009).