FoilTrack: a package to increase strain-resolution by improved X-radiographic image processing

ABSTRACT In high pressure multi-anvil experiments X-radiography is used to ascertain strain in deforming samples because the tooling prevents optical or other direct observations of the sample. The processing of these X-radiographic images to determine bulk sample strain is one of the limiting factors to making measurements closer to the strains and strain-rates that occur during mantle convection or the passage of seismic waves. Typically, sample deformation in these experiments is tracked by the displacement of high-contrast marker foils in X-radiographs. X-radiographs are treated individually or pairwise in a multi-step process that tracks the displacement of marker foils during experiments. Here I develop a new algorithm, FoilTrack, that treats all the X-radiographic observations in a single-step process, resulting in improved accuracy and consistency of length changes determined from X-radiographic images, as well as providing more realistic parameter uncertainty. The improvements are demonstrated using data from small-strain sinusoidal deformation experiments.


Introduction
The tooling and cells that form a multi-anvil high pressure experimental apparatus completely surround the sample [1,2].This prevents the use of optical imaging and Digital Image Correction (DIC) for deformation in the sample.The cells are mechanically 'soft' to prevent huge deformation of the sample while generating pressure.This is not a problem for hydrostatic multi-anvil experiments where only stable pressures are required but in deformation multi-anvil apparatus [e.g.[3][4][5][6] it means that displacement of the differential actuators does not directly correspond to strain in a sample.In-situ sample deformation during these experiments is therefore observed using synchrotron X-rays [e.g.[7][8][9][10][11][12].X-rays are able to penetrate the cells and X-radiographs are collected using visible-light cameras focused on phosphorescent YAG scintillators.Although it is possible to map full deformation fields in such high pressure experiments [3], it is much more common to track just the ends of the sample.The deformation of highly absorbing samples can be monitored in X-radiographs without the need to mark the ends of the sample [e.g.9].In general, though, silicate samples have similar X-ray absorptions to the cell assemblies and deformation is therefore tracked by placing high-Z metal foils (e.g.gold, platinum) at the ends of the sample (Figure 1).During experiments, the position and displacement of these foils in X-radiographs is used to determine the sample length, strain and strain-rates in the sample.The processing of X-radiographs to determine sample strain is based on finding the position, or change in position, of the shadows of the high-Z marker foils.X-radiographs are processed in a number of ways: manually [13], with image-J or with custom-written code [e.g.7,11,14,15]; many manuscripts though exclude any details of radiograph processing.These radiograph processing pipelines treat each X-radiograph as a single independent observation and use a multi-stage process to reduce the multi-pixel images to a set of positions or displacements, that are then further processed although they were errorfree displacement transducer measurements (e.g.Section 2.1 and references therein).Formal errors are not propagated between all steps resulting in formal errors that underestimate the true uncertainty in the measurements.The signal-to-noise ratio in the Xradiographs propagates through into the final positions, strain or strain-rate, providing a lower limit to the possible resolution or measurement precision.This becomes a problem when high pressure deformation experiments are aiming to get as close to Earth-like strain-rates and deformations as possible.Constant strain-rate deformation experiments are typically conducted at 10 −5 − 10 −7 s −1 [e.g.7,8,10-12,16], which are significantly faster than the typical strain-rate of , 10 −12 s −1 occurring in the bulk of the Earth's mantle.Slower strain-rates are in principle achievable in the apparatus but the allocation of synchrotron X-rays in short blocks (usually less than 2-3 days per allocation period) puts lower bounds on what is measurable in the time available.Seismic frequency, sinusoidal deformation is either studied using other techniques, small-strain amplitudes and insignificant pressure [strains c. 10 −6 ; pressures c. 200 MPa 17,18] or with high pressure synchrotron-radiography at higher pressures and larger amplitudes [19].But these applied strains are greater than the , 10 −7 strain that occurs during the passage of seismic waves.
Substantial improvements in resolution are needed to get closer to Earth-like strain-rates.These improvements are unlikely to come from significant advances in X-radiograph resolution because the images are already optimised within the geometrical constraints of the high pressure apparatus, the use of X-ray scintillators and visible-light cameras.Instead, improvements can be made by utilising the functional form of the deformation applied to the samples.The displacement-time profile applied to the samples is generally either linear or sinusoidal; at worst it can be approximated by a low-order polynomial.In this study, I outline the package FoilTrack which uses knowledge of the deformation profile to improve the precision and coherence of measured sample length changes (Section 2.2).Data fitting is performed as a single-stage process which propagates the errors inherent in the measurement.The improved algorithm is validated using experimental data from small-strain sinusoidal deformation experiments (Section 3).

X-Radiographic image processing
Multi-anvil apparatus has been positioned on synchrotron beamlines for more than 30 years.Sample lengths, from which strain is calculable, were initially measured manually in X-radiographs [13].Automated processing of radiographs followed soon thereafter [10,14] and has since become common place.There are two basic approaches to tracking the change in sample lengths: feature finding [14] and image-correlation techniques [10].Feature finding involves locating an identifiable feature (usually the minimum intensity of the maker foil shadow) in each radiograph; this approach is used in both manual [13] and automated processing of radiographs [14].

Prior algorithm
Image correlation and sub-pixel tracking of marker foil displacements in X-radiographs was first implemented for X-radiographs by Li et al. [10].They based their method on the image processing algorithms of Pratt [20] and Trucco and Verri [21].The basis of the algorithm is finding the minimum of the Sum Squared Differences (SSD) in pixel intensity between regions of interest in a reference radiograph (I r {X, Y} z ) and a search region in the comparison radiograph (I c ) where {X, Y} z is the array of pixel coordinates for the region of interest z.The SSD of the pixel intensities between the two radiographs is calculated for a range of offsets, o, as: For small total offsets, this function is approximately parabolic in o; an example set of D z r,c (o) values are plotted as black crosses in Figure 2(a).The sub-pixel displacement (d z r,c ) of the region of interest is found by finding the minimum of a cubic spline interpolated between the values of D z r,c (o).The offset range, o, used to determine the displacement, is generally +5 pixels (≈ 10 m) vertically in the radiograph (y-direction in Figure 1).Li et al.  5)).The small displacement (here <0.1 pixels) is reflected in the sinusoidal variation in SSD values at constant offset.The combination of displacements above and below a region gives the length changes (Equation ( 3)) which are plotted in Figure 3.
[10] were primarily interested in tracking comparatively large total displacements (. 10 pixels) during deformation experiments and tracked the displacement between consecutive pairs of radiographs (d r=1,c=2 , d r=2,c=3 ,…d r=N−1,c=N , where N is the total number of radiographs).In this implementation, the position of the region of interest changes as the marker foil moves and the total displacement of each region of interest is the sum of all prior displacements.The position of the marker foil in radiograph n r is: where P 0 is the initial position of the region of interest determined from the mean pixel Y values in that region.The sample length, l, is the difference in position of two regions of interest above and below the sample: Time stamps convert the length as a function of radiograph number, n r , into length as a function of time, t.Sample length (or strain) is fit to these values for the deformation profile applied in the experiment.This algorithm has been utilised to measure deformation in high-strain experiments [e.g.7,10,22] and small sinusoidally varying displacements during thermal conductivity measurements [23,24].Constant strain-rate deformation experiments use the change in length with time to determine the strain-rate by fitting a degree-1 polynomial to the  1 and 2. Blue dots are minima of the SSD polynomial for independent calculations of the displacement between each radiograph and the reference radiograph (Equation ( 3)).Red lines are fit to the minima calculated using the method of [25,26] (Section 2.1).Black lines are the length changes calculated from the surface fits to all the SSD data (Equation ( 5)) and thus not fit to the points, but ideally they should reproduce them.The black dashed lines are the secular length changes of the samples calculated from the SSD surface without the sinusoidal displacements.
strain.Data from small-strain sinusoidal deformation experiments are fit by detrending the data with a spline or polynomial function and then fitting the remaining sinusoidally varying length changes.The detrending of the data takes no account of the phase or amplitude of the sinusoidal deformation.
Hunt et al. [25,26] improved the precision in later thermal conductivity measurements by using the middle radiograph as a single reference (d r=N/2,c=1 , d r=N/2,c=2 , . . .d r=N/2,c=N ; N .100) and degree-6 polynomials rather than splines to find the SSD minima.The thermal conductivity studies were primarily interested in phase differences and not specifically concerned with the amplitude of the sinusoidal displacements.These studies pushed the limits of resolution with this image processing algorithm and could resolve the phase differences because of the multiple cycles of the driving wave that were observed.Hunt et al. [16] used this improved algorithm to analyse highstrain experiments.
Subsequently, we have attempted experiments in which both the phase and amplitude of the strain are critical.Neither the original algorithm of Li et al. [10] nor the improved versions discussed above were sufficiently precise or coherent for the requirements of the subsequent studies.In addition, the multi-step data processing routine does not propagate errors between steps, resulting in underestimation of the formal errors.However, any errors are likely to be dominated by the error in absolute sample length which Li et al. [10] stated to be +5.We have therefore further refined the algorithm to return more consistent and accurate displacements of the marker foils.

Refined algorithm
To gain precision in the image processing, rather than treating each pair of images in isolation and then interrogating the displacements, a new algorithm was developed that describes all the SSD for a sequence of images, D(o, t), as a single polynomial surface, S(O, t) Figure 2(b).The surface has degree p in offset (o), q in time (t) and the total degree of the polynomial surface is the greater of p and q [27].For p>q, the surface is: where C uv are the polynomial coefficients and O is a modified form of the pixel offset (o) that accounts for the sinusoidal displacement of the foils: where a is the amplitude, ϕ the phase and f the frequency of the displacement.If the amplitude, a, of the sinusoidal displacement is zero then O = o and Equations (4a) and (4b) are no longer sinusoidal.The surface S z (O, t) is fit to D z (o, t) by ordinary least-squares minimisation.The fit is performed simultaneously for all the regions of interest; the only common parameter between the regions of interest is the period of the driving force (f ).During testing, it was found that lower-degree surfaces do not fully capture the shape of the D(o, t) surface while higher-degree surfaces have artefacts in the fit and unrealistically low standard errors.
Figure 2(a) shows the D(o, t) values for the first (top) region of interest from an example dataset collected with a 100 s period.This clearly shows the significant increase in D(o, t) with offset as the intensity mismatch between the regions of interest in the reference and comparison images get larger.The sinusoidal displacements of the marker foil with time are reflected in the sinusoidal variation of D with time at constant offset.The fitted surface, S(O, t), is shown in Figure 2(b) and is a good description of the data.The residuals of the fit Figure 2(c) are small, not systematic with time and form a bow-tie shape in offset Figure 2 (c).The bow-tie shape is because at small offsets the intensity differences are small and the addition of noise to small differences squared has less effect than at large offsets where the intensity differences are larger.For example, the addition of 1 arbitrary unit of noise to an intensity difference of 3 increases the difference squared from 9 to 16 (difference of 7) whereas 1 arbitrary unit added to an intensity difference of 10 increases the difference squared from 100 to 121 (difference of 21).
The displacement of each region of interest, d z (t), is the minimum of the surface with respect to time: which is found by differentiation of the polynomial surface, S z (o, t), with respect to unmodified offset (o) .Because the reference image is now a single image, the position of each marker foil is now: rather than the sum it was previously (Equation ( 2)).The sample length can be calculated as previously (Equation ( 3)) but because P z (t) is an algebraic expression, length is instead calculated by combining the polynomial and sinusoidal components of the displacements algebraically, to give both a secular length change over the dataset and a single sinusoidal amplitude (A) and phase angle (Φ).
The largest source of error, in this calculation, is the absolute length of the sample [10].To minimise both the absolute and relative length error between data acquisitions, the regions of interest in the reference radiographs are positioned automatically, using features in the data (e.g.maximum gradient, minimum intensity).The setting of the region of interest is a feature-finding routine equivalent to the algorithm of Higo et al. [14].

Example data acquisition
The data used to validate the improved algorithm was collected from a small-strain sinusoidal deformation experiment, undertaken in the D-DIA [3,6] on beamline X17B2 at the NSLS, Brookhaven National Laboratory, New York.The high pressure assembly was the same as that used by Dobson et al. [22] and the samples consisted of high-purity zinc wire (1 mm long, 1 mm diameter; 99.9985 % metal basis, Puratronic from Alfa Aesar) bracketed by two corundum pistons (, 0.9 mm long, 1 mm diameter solid rods of Alsint-23 corundum, from Alfa Aesar, with ends polished flat and parallel).At the ends and between each of the samples was a 25m Platinum disk.
After compression of the experiment, over ∼2 hours, and heating to the target temperature, the samples were strained sinusoidally by driving the D-DIA's deformation pumps, during which, X-radiographs were acquired using a yttrium aluminium garnet scintillator and a visible-light camera, for 10 nominal periods.Between 20 and 40 X-radiographs were collected per driving period, with an exposure time of 0.3 s.Data were collected multiple end loads (pressures), temperatures and driving periods.Data were acquired, at each pressure and temperature condition, for all the required periods before the conditions changed and the cycle repeated.A typical radiograph is shown in Figure 1.The amplitude of the deformation was the minimum needed to observe sinusoidal strains in both the sample and elastic standard using the algorithm available when the experiments were performed (the prior algorithm, see Section 2.1).
The data were fit with both the prior (Section 2.1) and refined (Section 2.2) algorithms using the same regions of interest.The regions of interest were positioned automatically to minimise systematic errors in sample length.Horizontally, the regions of interest were centred in the bright part of the image and ended close to but not overlapping with the anvil shadows.The regions of interest not adjacent to the zinc sample (Figure 1, bottom red box) were centred over the minimum in a spline interpolation of the intensity profile; the width and depth of which remained very similar throughout the experiment.The regions of interest adjacent to the zinc sample became broader throughout the experiment as the platinum marker foil diffused into the zinc.To account for this, the regions of interest were centred over the maximum gradient (as interpolated by a spline) on the side of the foil away from the sample.
Various degrees were tested for the surface coefficients and p = 6 and q = 2 were deemed to be the best.It was found that lower-degree surfaces do not fully capture the shape of the D(o, t) surface while higher-degree surfaces have artefacts in the fit and unrealistically low standard errors.

Results
The data were initially processed using the prior algorithm to produce sample length changes with time (blue dots, Figure 3, Equation ( 3)).The samples start each measurement by shortening, this is consistent with the deformation imparted by the pumpswhich start by advancing, acting to shorten the samples.The sinusoidal length changes of the central zinc sample (Figure 3(b)) are larger than those in the alumina, which is to be expected because of its smaller Young's modulus.The length change in the bottom corundum samples is significantly noisier than those of the other samples.This is because of the much less defined marker foil in the region of interest (Figure 1).Better fits can be obtained by increasing the number of terms used detrending the data.However, the use of higher order background functions reduces the apparent amplitude of the sinusoidal deformation.As well as containing more random noise the bottom displacement also contains apparently systematic excursions from the expected sinusoidal signal (e.g. between 900 and 940 s).Such excursions have a strong influence on the phase and amplitude returned by the prior algorithm.They are possibly caused by lowintensity non-random features in the image [e.g.fixed pattern noise in the CCD, 28,29] that have a greater influence on the SSD values (Equation ( 1)) in low-contrast regions of the radiographs.

Sinusoidal deformation
To determine the amplitude and phase of the deformation the data were reprocessed using both the original and improved data methods; the results of which are shown in Figure 3.Both algorithms produce length changes that match the point-by-point sample length changes reasonably well.The prior algorithm (red lines) directly fits the length change data and the largest residuals are at the beginning and end of the time series.These arise from how the secular length change is dealt with by the prior algorithm [see Section 2.1, 25,26, for details].The refined algorithm (black lines) however, shows a much better correspondence with the length changes despite not directly fitting these points.The length change of the samples determined from the fits (Figure 3) is formally a degree-3 polynomial, because q = 3 (Equation (4a)), but often appears to be approximately a degree-1 polynomial because the higher order terms are small.
The sinusoidal deformation was driven with an exact period and for both fits the period of the fit was a free parameter.The refined algorithm recovers the known period of the data significantly better than the prior algorithm (Table 1, Figure 4).The normalised periods returned by the prior algorithm strongly depend on the driving period.In contrast, there is significantly less scatter in the normalised periods for the Refined algorithm and negligible period dependency (Figure 4).
The amplitude of the length changes in the zinc sample (dots, Figure 3 (b)) are not as well fit by the prior algorithm (red line) as by the refined algorithm (black line).The prior algorithm systematically underestimates the amplitude of the deformation, relative to that of the refined algorithm (Figure 5 (a)).The underestimation is more significant in the noisy bottom corundum sample where the low contrast between the foil and the background (fourth box, Figure 1) results in a poorer and very noisy displacement signal (dots, Figure 3 (c)).Consequently the frame-by-frame displacements (D(o), Equation ( 1)) and associated length changes are noisy and the fitted sinusoidal length change is poorly constrained using the prior algorithm (red line, Figure 3 (c)).The refined algorithm, on the other hand, returns a fit with a much larger amplitude from the noisy data (Figures 3(c) and 5(a)) because the SSD (D) values at the highest offsets (o) constrain the fit (Figure 2).
The phase values are the least well correlated of the fitted parameters Figure 5 (b).Most of the scatter is again in the bottom corundum sample and is caused by the low contrast of the marker foil in the fourth region of the X-radiograph.There is also a small period dependency in the plot of phases returned by both algorithms which is likely caused by the poor period fitting in the prior algorithm.The expected phase angle in the sample is less well known than the other parameters because the sample is not well coupled to the deformation pumps.Dissipation in the 3-4 m length of small-bore hydraulic pipe between the deformation pumps and the actuators is significant and is the likely cause of the separation of the phases by period Figure 5 (b).Although the deformation pumps were driven with the greater amplitudes in the short period data the amplitudes were lower.However, the improved period and amplitudes from the refined algorithm leads confidence to a corresponding improvement in the phases.
The errors returned by the prior algorithm are unrealistically small.The typical phase error for the top sample is less than 0.1 • and smaller than those of the middle sample.This is despite the larger amplitude in the middle sample providing visually a better constraint to the phase.Formal errors for the phase and amplitude of the middle and bottom Prior samples are very similar (Table 1) despite the obviously very different quality of the data (Figure 3).In contrast, errors for the phase and amplitude returned by the refined algorithm are significantly larger for the bottom sample than the other two samples.The refined algorithm also returns a typical phase error for the middle sample that is smaller than those of the top sample.Overall, refined algorithm shows remarkable improvement over the previous algorithm.The larger errors in the new algorithm are a more realistic representation of the errors inherent in the data than those from the prior algorithm where the multi-stage fitting process significantly underestimates the uncertainties.

Pseudo-constant strain-rate deformation
Although the new algorithm was designed for analysing sinusoidal deformation in experiments it is equally applicable to more traditional, constant strain-rate experiments.The instantaneous length-change-rate was calculated from the point-by-point lengths using a moving window of 11 images.The length-change rate is the slope of a first-degree polynomial and these are plotted in Figure 6 as dots.
These instantaneous length-change rates are closely comparable to the length-change rates determined by differentiation of the refined algorithm's fit to the data (black lines, Figure 6).This provides a validation of the refined algorithm but closer inspection of the two rate calculations shows that the length-rates determined from the point-by-point calculation underestimate the maximum, minimum and range of rates relative to that of the refined algorithm.It is unclear why the refined algorithm has a greater maximum strainrate when the comparison between the length changes appears to be reasonable (Figure 3).But as the refined algorithm produces more consistent amplitudes, phases and periods is it reasonable to infer that the strain-rates determined from the refined algorithm are also an improvement on the instantaneous strain-rates.Although unreported here, the strain-rate uncertainties determined by the refined algorithm are more realistic than those of the prior algorithm.

Distribution
FoilTrack is an open-source program licenced under the GPL-3.0licence.The scripts are available for download from https://github.com/S-Hunt/Foil-Track.The scripts are written for MATLAB 2020 or later and require the 'Image processing', 'Parallel Computing' and 'Statistics and machine learning' toolboxes.Documentation and working examples of the processing pipeline for the scripts are included as part of the GitHub repository.The scripts do not require any special computational hardware and have been tested on Mac OS, Windows and Linux versions of MATLAB.
The GitHub repository is 23 MB (2MB without example data files).The data included with the example is a cut-down version of the data discussed here.The images have been reduced in size (clipped), density (every 4th image is included) and duration (only images in the first 200 of 400 images are included).The example is therefore unable to reproduce the values presented in Table 1, but output values and files from processing just the example data are contained within the repository.

Conclusion
This study shows that significant increases in the coherence and accuracy of fitted parameters can be gained by treating X-radiographic images as a continuous data series rather than as individual observations.By treating the data sets in a single-step process more realistic uncertainties can be estimated for each parameter.Treating the data set as a continuous data series requires applying the known (or assumed) functional form of the deformation to the data at the earliest possible stage of processing.Applying this knowledge improves both the constraints on each individual fit and the coherence of the overall data set.The high degree of constraint provided to the fits by the SSD surface (Figure 2) will enable experiments to be performed at lower total strains while still accurately measuring the response of the sample.
Although sinusoidal or linear deformation is assumed here, it is equally possible to apply and process other deformation functions by the use of higher-degree polynomials.This new algorithm will allow measurements to be resolved at lower total strains, allowing experiments to be performed closer to the conditions found in the bulk of the deep Earth.
In general, keeping the data in its rawest possible form for as long as possible should be the guiding principle of data processing.The principle, employed here, of treating multiple images as a single series, is generally applicable to Digital Image Correlation techniques.The improvement in coherence of the data indicates that corresponding improvements in resolution and accuracy should be achievable in generic DIC datasets by assuming continuous behaviour across multiple images.Strain-rates calculated from the point-by-point data with a moving window of 11 images (blue dots) compared to the strain-rate from the refined algorithm's fit (black line).Note that point-by-point strain-rates underestimate the maximum and minimum strain-rates.

Figure 1 .
Figure 1.Example X-radiograph from the high pressure synchrotron experiment used to validate the software algorithms in this study.The dark areas at either side of the image are the shadows of the tungsten carbide anvils and the samples are observed in the bright central stripe.The red boxes are the positions of the regions of interest tracked between images.The scale of the image is 2 m/pixel.

Figure 2 .
Figure 2. Example of Sum Squared Differences (SSD) data and its fit with our improved algorithm.(a) SSD data, D(o, t), for a region of interest.(b) The best-fit surface, S(o, t), to the data and (c) the residuals (D(o, t) − S(o, t)).In (a) the first set of values, at t = 0, are highlighted by black crosses.The values of D(o, t) and S(o, t) are in arbitrary units squared.The marker foil displacement is the minimum of the SSD at each time step (Equation (5)).The small displacement (here <0.1 pixels) is reflected in the sinusoidal variation in SSD values at constant offset.The combination of displacements above and below a region gives the length changes (Equation (3)) which are plotted in Figure3.

Figure 3 .
Figure 3. Length changes in (a, c) corundum and (b) zinc samples calculated for the same data shown in Figures1 and 2. Blue dots are minima of the SSD polynomial for independent calculations of the displacement between each radiograph and the reference radiograph (Equation (3)).Red lines are fit to the minima calculated using the method of[25,26] (Section 2.1).Black lines are the length changes calculated from the surface fits to all the SSD data (Equation (5)) and thus not fit to the points, but ideally they should reproduce them.The black dashed lines are the secular length changes of the samples calculated from the SSD surface without the sinusoidal displacements.

Figure 4 .
Figure 4. Comparison of normalised periods (measured period/driving period) returned by the prior and refined algorithms.If both algorithms were identical and recovered the driving period the points would lie on the point (1,1).Standard errors are smaller than the symbols.

Figure 5 .
Figure 5.Comparison of the (a) amplitude and (b) phases of the sinusoidal deformation returned by the prior and refined algorithms.If both algorithms were identical the points would lie on the solid black line.Standard errors have been omitted for clarity but are listed in Table1.

Figure 6 .
Figure 6.Strain-rates calculated from the point-by-point data with a moving window of 11 images (blue dots) compared to the strain-rate from the refined algorithm's fit (black line).Note that point-by-point strain-rates underestimate the maximum and minimum strain-rates.
I r {X, Y} Intensity in refernce radiograph region denoted by {X, Y} z l

Table 1 .
Periods, phases and amplitudes of the sinusoidal deformation returned by both the original and refined algorithms.The precision of the given values is required to show the errors returned by the fitting algorithm.