TomoWarp2: A local digital volume correlation code

code


Motivation and significance
The primary interest in experimental mechanics is the understanding of the stress-strain responses of materials.In traditional mechanical tests, displacements and/or forces are measured at the boundaries of a specimen under load and a single scalar value for the associated strain is calculated making the hypothesis that strain is homogeneous.However, heterogeneous samples and nonhomogeneous deformation, especially in the form of localised deformation mechanisms, are exhibited by most materials.Therefore, to better understand material behaviour, full-field strain measurements that capture heterogeneous responses are needed.
Early approaches for full-field, local measurement of strain during mechanical testing include the use of stereophotogrammetry to compare photos of a specimen's surface, acquired at different stages of loading, and to map the displacement of points from one image to the next [1,2].The advent of digital cameras in the 1980s led to the development of computerised Digital Image Correlation (DIC) [3,4].DIC provides displacement vectors over a set of points in the images, from which local strains can be calculated to be used, for example, to constrain theoretical or numerical models.
DIC was initially based on photographs of a sample surface leading to the necessary assumption that the surface is representative of what happens in the volume.Non-destructive full-3D imaging techniques, such as X-ray tomography, allow 3D images to be acquired at different stages of deformation of a sample.Consequently, DIC can be extended to full 3D, i.e., Digital Volume Correlation (also known as volumetric DIC).Consequently, it is now possible to follow the evolution of 3D displacements throughout a sample and throughout a loading experiment (see, for example: [5][6][7]).This paper outlines a new open-source, Python/C-based code, TomoWarp2, for DVC.
TomoWarp2 has its origins in ''TomoWarp'' (a single-threaded C-code), originally developed to analyse 4D seismic images of hydrocarbon reservoirs [8] and later adapted to 4D X-ray tomography of deforming materials.The original codebase was used to study localised deformations in geomaterials, including  [6,[8][9][10][11][12].However, the desire for greater portability, parallel operation, flexibility in input image format and to facilitate further development prompted TomoWarp2.TomoWarp2 has been used in continuation of previous research topics in geomechanics, (for example: [13][14][15], and in more diverse research, including charge/discharge of batteries [16,17] and the failure of bones and bone-implant interfaces [18]. TomoWarp2 requires input of: (a) two (2D or 3D) images; (b) an estimated displacement range; (c) the size of the region to be tracked; (d) the spacing between the analysis points.The output measured displacements can be visualised and filtered through a post-processing routine and the resultant displacements can be used to calculate a tensor strain field.A number of advanced options are also available as described below.

Software description
TomoWarp2 belongs to the local variant of DIC techniques (as also implemented by Gates et al. [19] and in ''CMV3D'' by M. Bornert [20]).Local DIC approaches involve independent analysis of points within a reference volume (im1) using local subsets (''correlation windows'') about each point.The displacements of the points are determined by finding the best mapping of the corresponding subset in the reference volume with equivalent subsets in a subsequent volume (im2), by a search over a restricted ''search window'' defined by the user displacement range.The best mapping is determined by a correlation coefficient calculated over the subset.This approach provides displacements with a sensitivity of 1 pixel, but higher accuracy is usually required.Subpixel accuracy can be obtained by either interpolation of the local correlation coefficient field around its maximum in displacement space or by interpolating the image itself (which allows more degrees of freedom).The work of Pan and co-workers, [21], presents a comparison of sub-pixel approaches (in 2D).Such a local approach is distinct from global techniques, where a global mapping of one image into the other is sought (see: [22][23][24]).

Software architecture
TomoWarp2 is written in Python 2.7, making use of few external packages (principally numpy [25], scipy [26], tifffile [27]).The program is parallelised using the Python multiprocessing libraries with a specifically written SIMD parallelisation scheme that employs global queues with the necessary features (blocking and read/write semaphores) to manage synchronisation.To speed up the analysis, the pixel search has been written in C and connected to the main part using the Swig package [28].
The TomoWarp2 program is highly modular and consists of 3 principal processes that run concurrently: a single management process (DIC_setup), a single process for handling data loading and data communication (data_delivery_worker) and a number of identical processes performing the point-by-point correlation (DIC_worker).This architecture has been adopted for its simplicity and its amenability to running on massively parallel systems (e.g., easily portable to MPI).
DIC_setup calculates the extents of data required in im1 and im2 for each node and communicates the global extent of the required data to data_delivery_worker.The information about the nodes and their extents are individually placed into a global job queue to be read by the DIC_worker(s).The results are collected as they are returned by each DIC_worker.When all nodes have been processed, the output is written to a tab-separated-value (TSV) file containing the number, coordinates, displacement and rotation components, the correlation coefficient, and the potential error code for each node.
The management of the image data is handled by the data_delivery_worker.This task is responsible for (a) loading the necessary image data and (b) replying to individual requests for data from the DIC_worker(s) made on a global queue.Image data are fed back to each worker on a worker-specific queue.
Numerous DIC_worker(s) are initialised through a multiprocess call and all run concurrently (synchronised with input-job and output-results queues).Each DIC_worker initialises and blocks the job queue until the queue is filled with nodes by DIC_setup.The job package includes the node number and the extents of data required in im1 and im2, based on which a data request is then made to data_delivery_worker.If both data sets are received successfully, then the pixel search is executed, returning the integerdisplacement for the best correlation and its corresponding correlation coefficient.If requested, a sub-pixel optimisation of the result is performed.The displacement of the best match and correlation coefficient are placed into the results queue (which is read by DIC_setup and converted into displacements in the global reference frame).Each DIC_worker works on only local information, which simplifies book-keeping and, thus, the structure of these workers.

Software functionalities
The subpixel refinement of the displacements has two different implementations.The fastest (CC-interpolation) consists in finding the maximum of an function interpolating the correlation coefficient at the 27 points surrounding the best integer-displacement match.Another implemented option (image interpolation) numerically transforms the subset of the deformed image at new subpixel positions according to an optimisation procedure until an optimal mapping is determined.The latter is more computationally expensive, but allows rotational degrees of freedom to be taking into account.
A ''multi-scale'' calculation can be performed by interpolation of a prior estimate of the displacement field and performing the correlation over a finer grid with optimised search ranges; this significantly reduces the computational time.Prior fields are input as TSV files in the same format as the output from the program, i.e., they can be generated by a previous run of the code or by some external process.The correlation is generally performed over a regular rectangular grid of nodes.However, prior fields can be used to provide coordinates of arbitrarily distributed node points (in addition to providing the prior estimate of displacements).Reprocessing can also be performed for a subset of nodes that did not give satisfactory results in a previous correlation.The nodes to be reprocessed can be selected based on the output error or a threshold on the correlation coefficient in a given prior.
Post-processing of the results by median filtering and outlier removal is provided.In interactive mode, i.e., when the GUI is used, these filters can be applied several times and the result displayed and saved in a new TSV file.Tensor strain fields can be calculated using different frameworks, including with tetrahedral or cubic elements, and Geers theory [29] to calculate strains centred with respect to the node.

Illustrative example of effective use of TomoWarp2
TomoWarp2 has been used to measure the displacements and strain fields in a large range of different experiments.Here we select, as an illustrative example, a study based on X-ray tomography images of a Bentheim sandstone sample acquired before and after deformation in triaxial loading conditions [14].Fig. 1(a) shows horizontal and vertical slices through the X-ray tomography volumes of the sample acquired after deformation.The image  volumes (before and after deformation) are 1800 × 1800 × 2300 voxels 3 .
The displacements through the sample are expected to be very different just from the boundary conditions: the bottom was fixed during loading and the top displaced.Therefore, running a calculation with a fine grid of nodes covering the entire range of possible displacements in the three directions would be unnecessary and very computationally expensive.Consequently, the iterative multiscale facility of TomoWarp2, was utilised (see Fig. 2): (a) a first analysis on a coarse grid (node spacing, NS=400 px) and a large search window (SW=± 60 px) to find the minimum and maximum displacements in the three directions; (b) a second analysis with a less-coarse grid (NS=100 px) using the previous results as SW; (c) a fine grid calculation using the results from the previous as a prior (NS=40 and SW=± 2 px); the analysis is repeated for the points that did not correlate well, using (d) the original large SW and (e) again using the image-interpolation subxpixel refinement mode; (f) a median filter is applied to the displacement fields before strain calculation in the large strain framework (Fig. 1(b-c)).Where the localised deformation consists of a strong discontinuity, the high strain values can be used for fracture detection, e.g., the fracture appears in yellow, but the physical meaning of strain is lost.

Impact
One of the main focuses of this project was to obtain user friendly and flexible code.In this perspective, a big effort has been made to simplify the input setup through both a text input file and a GUI.Moreover, the number of required input parameters has been minimised and an automatic detection of image format has been implemented, while the possibility to specify a number of parameters allows the refinement and optimisation of the computation.The main output of the code is a TSV file, but all the results, including those from post-processing, can be saved in common image analysis software formats, which avoids the need for further conversion.The GUI also allows measured displacements (and associated histograms) to be displayed, which facilitates filtering and post-processing.
Another key aim is to enable continuous development of the code, e.g., to provide new features that can help users answer their specific research questions.This is facilitated by the architecture of the code, and the simplicity and independence of the single processes.Furthermore, the open nature allows users to, themselves, add-to or adapt the code to suite their needs.Moreover, the parallel computation and the possibility to run on server farms can significantly reduce computation times (and opens the door to more detailed analysis of larger image volumes).
TomoWarp2 is currently used by more than 15 research groups, situated in France, Sweden, UK, Spain, Chile, Norway, USA and Japan, in research ranging from geomechanics to biomechanics and Li-ion batteries.

Conclusions
TomoWarp2 is a Python based code for full-field vector displacement measurement between 2D or 3D image sets, based on local digital image correlation.Complete strain tensor fields can be computed from the displacements, e.g., for detailed characterisation of mechanical behaviour of a specimen in a deformation test.The simple structure of the input file, together with the GUI and the possibility to use different formats of input and output images, make the code straight forward and quick to use, while the code architecture enables new features to be implemented easily.Future developments will include discrete correlation in the style of [6] (to track displacements and rotation of grains) as well as an MPI implementation.

Fig. 1 .
Fig. 1.(a) horizontal and vertical slices from the X-ray tomography of the deformed sample.Corresponding slices from the calculated (b) maximum shear strain and (c) volumetric strain fields.

Fig. 2 .
Fig. 2. Vertical displacement of (a-c) multi-scale analysis from a very coarse grid to a fine grid, (d-e) Repeated analysis for the points that did not correlate well with larger search window and image interpolation subxpixel refinement mode, and (f) Application of a median filter.Positive values indicate downwards movement.