New GPU Software Maps 3D Strain Inside Materials at Voxel-Scale Resolution

When engineers want to know exactly how a material deforms on the inside—how a crack threads through sandstone, or how stress concentrates around glass fibres in a polymer—they turn to X-ray computed tomography. The technique delivers stunning three-dimensional images of a specimen’s internal structure, but turning those images into quantitative maps of displacement and strain has long been a computational bottleneck. Now a researcher at the University of Freiburg has released an open-source tool that promises to dissolve that bottleneck, running the entire analysis on a single commodity graphics card while resolving deformation at a scale of just a few voxels.

The software, called voxelDVC, was developed by Georg Ganzenmüller and published in the journal SoftwareX. It performs digital volume correlation, or DVC, a technique that compares tomographic scans of a specimen taken before and after deformation. By registering the two volumes against each other, DVC recovers the full three-dimensional displacement field: a point-by-point account of how every region of the material moved. The result is essentially a four-dimensional picture of mechanical behaviour, with strains computed from the gradients of the recovered displacements.

DVC comes in two flavours. Local methods chop the image into small independent subvolumes and correlate each one separately, which is fast and simple but produces noisy, discontinuous results. Global methods, by contrast, represent the displacement field as a continuous function expanded in finite-element shape functions and determine all the unknowns in a single coupled minimization. Continuity is guaranteed by construction, and the formulation admits regularization—mathematical penalties that suppress implausibly rough displacement fields, ranging from simple smoothness constraints to schemes grounded in mechanical equilibrium. Regularization is what makes very fine meshes usable, allowing elements as small as a single voxel, as earlier work by Leclerc and colleagues demonstrated on small volumes.

The trouble with fine global meshes has always been the linear algebra. The Gauss-Newton Hessian that drives the solve is sparse, with roughly 81 nonzero entries per row for trilinear hexahedral elements, but assembling it touches every voxel with a dense 24-by-24 element update, and its storage grows with the number of degrees of freedom—which scales as the inverse cube of the element edge length. Worse, a sparse direct factorization of a three-dimensional grid problem requires operations growing with the square of the degrees of freedom and memory growing as the degrees of freedom raised to the four-thirds power, even with an optimal nested-dissection ordering. For a modern micro-CT scan containing billions of voxels, that memory bill quickly becomes prohibitive, forcing researchers to choose between meshes fine enough to resolve strain concentrations around fibres, pores or cracks and meshes that actually fit in a computer’s memory.

voxelDVC takes a fundamentally different route: it keeps the full finite-element space but eliminates the assembled matrices altogether. Both the correlation Hessian and the regularization operator are applied matrix-free, element by element, an approach well established in high-performance finite-element computing and, on structured voxel grids, in giga-voxel topology optimization. Because the solver works on a structured Hex8 mesh aligned with the voxel grid, the mesh connectivity follows from simple integer index arithmetic and never needs to be stored. The voxels themselves serve as quadrature points, so the data term is a direct sum over voxels rather than an approximation of an integral, and the reference image is never interpolated. Memory then scales with the number of image voxels, independently of the element size, and the whole solve runs on a single GPU.

The numerical machinery is elegant in its economy. The Hessian is a sum of rank-one terms, one per voxel, and its action on a vector is computed without ever storing it: one gather of nodal values, two small tensor contractions per voxel, and one scatter-add, all vectorized over chunks of elements on the GPU. The normal equations are solved with Jacobi-preconditioned conjugate gradients, which requires only this operator application. A comparison with explicitly assembled finite-element matrices on a small mesh confirmed that the matrix-free operators, their diagonals and the quadrature weights agree with the assembled versions to within round-off error. The deformed image is sampled at non-integer positions using cubic B-spline interpolation, a choice that matters: the ground-truth tests showed that conventional trilinear sampling introduces a systematic sub-voxel bias of up to 0.03 voxels, whereas cubic B-spline sampling keeps the mean error at or below 0.003 voxels.

Large displacements, which would defeat a single-shot correlation, are handled by a two-stage strategy. A preprocessing step first estimates an affine map—rotation, per-axis scaling and translation—between the two volumes and crops them to their valid overlap. The main solver then works through a coarse-to-fine image pyramid, typically three levels at which the images are subsampled by factors of four, two and one, with the mesh refined accordingly at each level. The solution is passed upward simply by rescaling nodal displacements, without any interpolation between meshes. The regularization length is scaled with the voxel size at each level so that the regularization weight is recomputed consistently, and the Gauss-Newton loop on each level stops when the relative update falls below a threshold or the iterations stagnate at the float32 round-off floor.

Accuracy was verified with a built-in ground-truth framework that imposes a known displacement field on a real reference image—synthesized with cubic B-splines so that the test does not share the solver’s interpolation assumptions—and then checks how well the pipeline recovers it. On a CT dataset of short-glass-fibre-reinforced polyamide, imposing a 5 percent uniaxial compression with local perturbations caused by the stiffness contrast between glass fibres and matrix, the solver recovered the interior displacement field with root-mean-square errors as low as 0.041 voxels once the regularization length was properly tuned. The tests also revealed a subtle systematic bias from the H1 regularizer, whose natural boundary condition pulls normal strains toward zero near unconstrained domain faces over a distance on the order of the regularization length—exactly the kind of effect the bundled parameter-sweep tool is designed to expose before a user trusts the solver on data without independent ground truth.

Real-data performance is striking. On a 753-by-401-by-353-voxel scan of a notched sandstone column fractured under compression, voxelDVC captured the inclined fracture band and a secondary crack in 34 seconds on an NVIDIA L40 S GPU, with results qualitatively matching the established SPAM code. Halving the element size to resolve the fracture more sharply took 82.7 seconds—changed with a single command-line option. Benchmarks on a laptop-class RTX 4050 GPU showed peak memory essentially independent of element size: a fixed overhead of roughly 1.05 to 1.3 gigabytes plus 44 to 54 megabytes per million reference voxels. That translates to cubic volumes of about 500 cubed voxels on an 8-gigabyte card and more than 900 cubed on a 48-gigabyte card, all without a cluster or a large-memory workstation. For a fixed element size, runtime grows linearly with volume on large datasets.

The implications reach well beyond one laboratory. Experimentalists working with in-situ CT on fibre-reinforced composites or geomaterials have often had to fall back on general-purpose image-registration tools for lack of an accessible fine-mesh global DVC alternative. Because voxelDVC’s memory footprint depends only on the number of voxels, the element size can now be chosen according to the physics of the problem and the texture of the image rather than the limits of the hardware. Just as importantly, the low cost per solve makes the bundled ground-truth and sweep tools genuinely practical: a multi-run calibration against a synthetic known-displacement problem built from the user’s own reference image offers a principled way to choose solver settings for any dataset, rather than relying on defaults tuned for a different material. The software is released under the MIT license on GitHub, with an archived snapshot on Zenodo, and it forms the basis of a forthcoming systematic assessment of global DVC at voxel-scale resolution. Cracks, it should be noted, still appear as localized strain bands rather than true displacement jumps, since the continuous finite-element field cannot represent discontinuities—but for the growing community mapping deformation inside opaque materials, the resolution ceiling has just been lifted dramatically.

Subject of Research: A matrix-free, GPU-native global digital volume correlation solver for measuring 3D displacement and strain fields from high-resolution computed tomography

Article Title: voxelDVC: a matrix-free, GPU-native global digital volume correlation solver for high-resolution computed tomography

Article References: Ganzenmüller, G. C. (2026). voxelDVC: a matrix-free, GPU-native global digital volume correlation solver for high-resolution computed tomography. SoftwareX, 36, Article 103111. https://doi.org/10.1016/j.softx.2026.103111

Image Credits: AI Generated

DOI: Not provided

Keywords: digital volume correlation, computed tomography, GPU computing, finite elements, strain mapping, open-source software, matrix-free solvers, image registration, materials testing, fracture, composites, geomaterials