Correction for Mechanical Inaccuracies in a Scanning Talbot-Lau Interferometer

Grating-based X-ray phase-contrast and in particular dark-field radiography are promising new imaging modalities for medical applications. Currently, the potential advantage of dark-field imaging in early-stage diagnosis of pulmonary diseases in humans is being investigated. These studies make use of a comparatively large scanning interferometer at short acquisition times, which comes at the expense of a significantly reduced mechanical stability as compared to tabletop laboratory setups. Vibrations create random fluctuations of the grating alignment, causing artifacts in the resulting images. Here, we describe a novel maximum likelihood method for estimating this motion, thereby preventing these artifacts. It is tailored to scanning setups and does not require any sample-free areas. Unlike any previously described method, it accounts for motion in between as well as during exposures.

but for higher sensitivities, the Talbot effect can be utilized [1], [2], which creates an intensity pattern resembling the grating at distinct propagation distances.Resorting to the Talbot effect also has the advantage that it can be induced by a diffraction grating, which preserves more of the source's flux.The pitch of this pattern is often much smaller than the detector's resolution and thus cannot be resolved directly.In this case, an additional absorption grating is placed in front of the detector, translating this intensity pattern to an intensity modulation depending on the relative position of the gratings [3].
Standard X-ray tubes providing sufficient flux for medical imaging, however, do not feature the required spatial coherence; the intensity pattern would vanish due to focal spot blur.This limitation is circumvented by placing a third grating close to the source [4], effectively splitting it into an array of line sources.Each aperture is sufficiently small to create an interference pattern.Even though these are mutually incoherent, the distance of these apertures is chosen so that they all create the same intensity pattern, which makes the technique accessible for clinical applications.
In the past decade, this idea has been developed from proof of concept, first demonstrations of X-ray dark-field imaging [5] to various trials on animal models, human samples, and even patients [6].Imaging of the lung has become one of the most promising candidates for future clinical application, since there is no non-invasive imaging modality available so far that can probe the lung's microstructure, and because its alveolar structure causes a large amount of coherent ultra-small angle scattering, which amounts to a strong contrast in the darkfield image.After encouraging results in in vivo pigs [7] and human cadavers [8], a clinical dark-field chest radiography demonstrator system has been conceived and built [9], which is also subject of this paper.
In order to achieve a field of view that is larger than the gratings currently available, the interferometer has to be moved with respect to the sample.By detuning the interferometer, sampling the interference pattern and scanning the sample is done simultaneously, moving either the sample [10] or the entire interferometer including the detector [11].This comes at a cost: the alignment of the gratings is harder to maintain during acquisition for a moving interferometer.For thoracic imaging of patients, the scan has to be quick to avoid motion artifacts from breathing, which makes it even harder to maintain the grating alignment.
Several methods have been proposed to remove artifacts caused by random alignment inaccuracies.A recently published method based on deep learning [12] requires repeated training whenever the flat field changes, and is built upon the assumption of a weakly scattering sample, which is not valid for our setup.Methods based on principal component analysis (PCA) [13], [14] require a fixed position of the sample with respect to the interferometer.In [15], the authors provide expressions for the artifacts resulting from a small deviation from ideal stepping, and propose a removal procedure based on the assumption that image and artifacts are uncorrelated.Because sampling in a scanning acquisition scheme is highly irregular, this approach is not compatible with scanning acquisition either.The algorithm described in [16] does not have that limitation: grating positions are estimated assuming that the correct ones yield the most smooth attenuation and darkfield images.Again, this requires an image that is entirely uncorrelated with the artifacts.In [17], this issue becomes apparent: describing an improved variant of the aforementioned algorithm, the authors note that estimated grating positions differ on whether differential-phase or dark-field images are assumed to be smooth.
For scanning setups (as opposed to phase stepping), only one procedure has been discussed so far [18]: rotation, offset and curvature of the gratings are estimated by assuming that the true grating position maximizes the correlation between reference and sample raw image.This works best in samplefree areas on the detector, which are not always available in a clinical setting due to careful collimation.
A completely different approach solely relies on a leastsquares fit, extending the cost function of per-pixel phase retrieval by per-exposure parameters.While being less robust, this has the advantage that there are no additional assumptions on the image.The most popular algorithm for that is the one described in [19], which has successfully been applied to grating-based X-ray imaging [20], [21].Unfortunately, their approach of alternating optimization for images and per-shot parameters shows very slow convergence for data acquired by scanning.An extension to more mechanical degrees of freedom has been demonstrated [22], still relying on some kind of alternating optimization.Only recently, a superior optimization scheme for the very same cost function has been suggested [23] for linear displacements in static interferometers.
In this paper, we present a novel way of mitigating vibration artifacts in a large, quickly moving scanning grating interferometer.Formulated as a maximum likelihood problem, it incorporates for the first time both motion in between as well as during exposures, which is shown to be impossible for stepping setups.Alternating optimization is replaced by a gradient descent scheme, which is both faster for scanned data sets and opens up further possibilities for optimization, like quasi-Newton methods or preconditioning.

II. METHODS
Fig. 1.Setup and origin of vibration artifacts.The setup is a Talbot-Lau interferometer, operated in scanning mode; the active area sweeps across the detector as the interferometer is rotated upward around the focal spot of the X-ray source.Without loss of generality, the third grating G 2 is subject to vibrations.If the grating itself is assumed to be rigid, two degrees of freedom are particularly notable: a "tip" movement φ s 0j and a "tilt" movement φ s 1j .

A. Image formation model
For a Talbot-Lau interferometer, the expectation value of the intensity measured by a single detector pixel can be approximated as [4]: Here, the interferometer itself is characterized by the socalled flat field, i.e. the mean intensity t f , the visibility of the reference pattern v f , and its position φ f with respect to the analyzer grating in units of p 2 /2π, where p 2 is the pitch of the analyzer grating.The influence of the sample is described by its transmittance t as well as a visibility reduction v and a shift φ of the reference pattern perpendicular to the grating lines.If one neglects propagation effects within the sample, t, v, and φ are related to projections of the sample's attenuation coefficient µ, linear diffusion (dark-field) coefficient ϵ, and refractive index decrement δ along the beam path as follows: where n is the unit vector from source to detector pixel, c v > 0 and c φ > 0 depend on the grating pitch and the setup geometry [24]- [26], and ∂ g denotes the partial derivative taken perpendicular to the grating lines.

B. Acquisition scheme and basic image reconstruction
In the following, we will assume a full-field detector whose position is fixed with respect to the sample 1 .In order to be able to reconstruct t, v, and φ, one has to take a series of intensity measurements I j with varying reference phases φ f j .In a fringe scanning scheme [10], [27], this is done by slightly detuning the interferometer such that there are moiré fringes   I. perpendicular to the scanning direction.By taking acquisitions while moving the interferometer across the sample, we obtain the necessary variation of the reference phase (i.e.local grating alignment) without moving the gratings with respect to each other.
We obtain transmission t, differential phase φ, and darkfield v by solving the following least-squares problem: t, φ, v = arg min t,φ,v j l I j Īj .
Here, we use a data-weighted least-squares model l to approximate the log-likelihood that arises from photon shot noise (variance σ 2 j ), neglecting electronic noise.The optimization problem ( 7) can be solved directly; for an X-ray grating interferometer, this has been described first in [28].By a change of variables, we can reformulate (5) as a linear equation, The resulting linear least squares problem, can be solved directly by standard techniques from linear algebra, such as LDL ⊺ Cholesky decomposition.Then, the visibility reduction and the differential phase are given by:

C. A model for mechanical vibrations
Experimentally, we observe "tip" and "tilt" vibrations of the interferometer, whereas the gratings themselves are rigid, see Fig. 1.Equivalently, these vibrations can be described as a random displacement of the third grating G 2 : For the jth exposure, the grating is shifted perpendicular to the grating lines by φ s 0j as well as rotated around the optical axis by φ s 1j .For small angles and taking the translational symmetry of the gratings into account, this rotation is equivalent to a local displacement of the third grating, which is proportional to the location's distance to the center of rotation, measured along the grating lines.Accounting only for motion between exposures results in the following model, x i being the x coordinate of the ith pixel.The coordinate system is chosen so that x i of a pixel from the left/right edge of the detector is ±1.
Furthermore, motion may also happen during one exposure lasting over a time T .We obtain where Taking only uniform movement into account, this leads to a visibility reduction, where sinc x := (sin x)/x.If the grating displacement αij T happening during one exposure is sufficiently small, we can make the following approximation: where v s 0 , v s 1 , v s 2 denote the motion-induced visibility reduction, grouped by its degree in x: We would like to point out that this approach is valid for higher-order vibration modes as well: vibrations whose displacements can be described by nth degree polynomials of x i result in 2nth degree visibility polynomials of x i .
For our setup, we found "tilt" vibrations to be slow enough that a constant visibility drop v s j := v s 0j is sufficient, yielding the following model: In order to obtain images t, v, φ and vibration parameters vs , φs 0 , φs 1 , we fit this model to the measured intensities, t, v, φ, vs , φs 0 , φs Notably, the solution is not unique.For any choice of λ ∈ R 3 , (21), since all occurrences of λ 0 , λ 1 , λ 2 cancel out in (19).We select one solution by requiring median ( φs 0 ) = 0, median ( φs assuming that there are no long-term drifts between reference and sample measurements, and that the motion in the sample scans follows the same random distribution as the one observed in the reference scans.

D. Gradient-based optimization
In order to find the most likely images t, v, φ and grating motion parameters φs , vs , we utilize that calculating the images while keeping φ s , v s constant is computationally cheap.We define a new cost function L where this part of the optimization has already been done, As per the multivariable chain rule, its gradient is given by (ξ s being v s , φ s 0 or φ Conveniently, this corresponds to first calculating the image and then taking the partial derivatives of the log-likelihood with respect to vibration parameters.L is minimized using the L-BFGS algorithm [29], [30]. 1) Subsampling: In order to increase performance, we may take only a subset of all pixels into account.Using only every nth detector column is approx.n times faster.This is similar to assuming a detector that is just 1/n as wide as the actual one.We found that e.g.n = 20 does not visibly alter the resulting images.
There are two limitations to subsampling: First, aliasing needs to be avoided: the value of the subsampled cost function does not change when substituting one of the tilt displacements φ s 1j by φ s 1j + w/n • 2π, w being the number of pixels in one detector row, which is not the case for the full cost function.Second, the estimation becomes more noisy, since it is fitted to less measured data.
2) Preconditioning: For faster convergence, we aim to rescale the problem so that the diagonal of its Hessian is approximately one.The diagonal of the Hessian of L is given by We will neglect everything but the leading term on the right side.The accuracy of this approximation can be checked empirically for a specific setup, but may be motivated by the fact that in the limit of many exposures contributing to each pixel, the dependency of the images on one single exposure's parameter ∂ ζi /∂ξ s j vanishes.Furthermore, we approximate this expression by its statistical expectation value ⟨•⟩ (w.r.t. the photon statistics) to obtain an expression akin to the Fisher information in standard phase retrieval [31], since   To obtain an estimate that depends less on the initial guess, we take the average over all possible grating positions, neglect the influence of motion during the exposure (v s ≈ 1) and obtain Expanding this expression and similar ones for the remaining variables up to the vanishing third order in v f ij vi yields: This can be used to rescale the problem so that the diagonal of its Hessian is nearly constant.I.e., instead of L, we optimize where ⊘ and • 1 2 denote the element-wise division and square root, respectively.

E. Dealing with sample motion and non-convexity
If the sample or a part of it is moving during acquisition, which is inevitably the case for a human patient and the relatively long acquisition time of a scanning setup, the algorithm described above falsely attributes that inconsistent transmission data to grating motion, often resulting in strong additional artifacts.This is mitigated by excluding pixels subject to motion, which can be identified by examining the "intensity-normalized" residual εi for each detector pixel, εi = min ti,vi,φi j before starting the optimization: if the transmission in said pixel is inconsistent, the stepping curve is no longer a good model for the measured intensities and εi becomes large.Thus, pixels where εi exceeds some empirically established threshold are excluded.We note that this choice is not crucial: discarding some pixels that are not actually subject to motion will not have a significant influence on the outcome, as long as enough pixels of each shot are assumed to be valid.We take advantage of the same redundancy in Sec.II-D when taking only a subset of the detector's columns into account.After determining the vibration parameters, the final images are retrieved taking all pixels into account.
The optimization problem stated in ( 21) is non-convex.We found a combination of two measures suitable to avoid local minima: First, we replace the phase image by its circular mean during the first iterations: This is similar to setting the transmission to unity in [32].Second, we start by optimizing φ s , and add v s to the optimization problem when φ s has already mostly converged.

F. Inapplicability of the algorithm to a stepping acquisition scheme
For a phase stepping acquisition scheme, there is no unique solution when optimizing for per-exposure phase and visibility.Suppose we found the correct solution locally minimizing the log-likelihood.Then, the model intensities Îij fitting the data optimally are given by Îij = ti 1 + vs j vi cos φs j + φi .
The key difference to (19) is that the flat field does not change from exposure to exposure.Hence, unlike for inhomogeneous gratings used in a scanning acquisition scheme, distinguishing between sample and flat field is not necessary.Just like in (8), per-exposure and per-pixel quantities can be separated:  By a change of variables, ĉs j := vs j cos φs j , ŝs j := vs j sin φs j , ĉi := ti vi cos φi , ŝi := ti vi sin φi , we obtain an equivalent linear equation, We can construct many different solutions • ′ that fit the data just as well, e.g.these ones, parameterized by λ ∈ R: To illustrate these spurious solutions, we obtain expressions for the resulting images by applying (10): Thus, this kind of optimization scheme is not suitable for stepped data when grating motion during exposures is taken into account.This is also the case for flat fields in scanning setups, which are measured in a stepping acquisition scheme.
In theory, this situation also occurs in a perfectly uniform scanning interferometer, when t f ij and v f ij are independent of j, and the flat-field phase is separable, The latter is the case if the flat-field phase gradient in scanning direction is constant.In practice, this is unlikely and may easily be avoided by detuning the interferometer slightly.

A. Simulation study
We conceived a simulation study for a scanning X-ray Talbot-Lau interferometer affected by mechanical vibrations roughly resembling the clinical prototype.The simulation is monochromatic at an energy of 45 keV.The distance from source to the patient centre is 1.5 m, the distance from source to detector is 1.8 m.We assume a detector having 1 mm × 1 mm pixels, with a total size of 45 cm × 51 cm, resulting in a field of view of 37.5 cm × 42.5 cm.
The interferometer features horizontal grating lines, and has characteristic parameters c v = 1, c φ = 40 µm and a spatially constant flat-field visibility v f = 0.3 and illumination t f = 2•10 4 photons.The active area of the interferometer is 7.5 cm high and as wide as the detector.There are two moiré fringes in the center, and 2.5 moiré fringes at the left and the edge of the grating.The active area moves 0.25 cm in between two exposures and traverses the entire detector, resulting in 30 exposures of each detector pixel, and 209 exposures in total. .Retrieved attenuation (a-d), dark-field (e-h) and differential phase images before (first column) and after (second column) estimation of grating vibrations.The third and the fourth column show zoom-ins of the uncorrected and corrected images, respectively.Zoom-ins of the dashed regions are shown in Fig. 7.For the logarithmic dark-field images, a window from − ln 1.05 to − ln 0.55 was chosen, for the differential phase images, it ranges from -0.2 to 0.2.The zoom-ins feature narrower windows.The scan was performed at 70 kVp, with a total scan duration of 7 s, and an estimated effective dose [33] of 43 µSv.The field of view is 37 cm×37 cm.
The phantom is the FORBILD thorax phantom [34], extended by a refractive index decrement and a dark-field material property, i.e. linear diffusion coefficient.Fig. 2 shows the phantom data.Assumed material parameters are given in Table I, which were calculated from [35] except ϵ, whose choice is motivated by the values typically found in images of patients.The flat-field phase and three exemplary raw images are displayed in Fig. 3.
Simulating mechanical vibrations, φ s 0j and φ s 1j were drawn from a normal distribution with µ φ = 0, σ φ = 2π/16.The flat fields are assumed to be known exactly, i.e. contain neither noise nor vibrations.The visibility reduction v s due to motion was calculated by |1 − |u||, where u was drawn from a normal distribution with µ u = 0, σ u = 0.2.This choice corresponds to vibrations that are approximately twice as large as what we observe experimentally.The probability distribution of v s was chosen to avoid negative visibilities, which are physically impossible, as well as visibilities larger than one: due to the choice of perfect flat fields, vibrations in the sample scan will always result in a visibility reduction.In the reconstruction, we assume max (v s ) = 1 accordingly.
Fig. 4 shows the resulting images before and after applying the proposed algorithm.Even without correction, there are no visible artifacts in the transmission image for these simulation parameters.The uncorrected dark-field image is not only impaired by fringe-like artifacts, it is also quantitatively wrong, since the vibration-induced loss of fringe visibility is falsely interpreted as a dark-field signal.All visible artifacts from vibrations are successfully removed, both in the dark-field and in the differential phase image.A slight global ramp in x direction is visible in the corrected differential phase image, due to the ambiguity mentioned in Sec.II-C.

B. Experimental validation
After a first validation in simulations, we applied the algorithm to a real scan of a patient.The setup is the one described in [9], a scanning-type Talbot-Lau interferometer with horizontal grating lines and an active area of 42 cm × 6.5 cm, which is moved in vertical direction to reach a total field of view of 42 cm×42 cm measured in the detector plane.One scan takes 7 s, while acquiring 195 partially overlapping frames, resulting in at least 24 exposures of each detector pixel.Raw images are recorded using a flat-panel detector; the rotating anode X-ray source is operated at 70 kVp.Regardless of whether the vibration correction has been applied, the images are corrected for the effects of Compton scatter, beamhardening, and dark-field bias [36] in regions with particularly low transmission.Adapted Unique processing (Unique X-ray image processing software, Philips Healthcare, The Netherlands) was applied to attenuation images to obtain an appearance akin to regular chest X-ray images [37].The darkfield image has been iteratively denoised to achieve spatially constant noise variance [38].The resulting images with and without vibration correction are shown in Fig. 5. Unlike in the simulation study, even the attenuation image exhibits some faint, yet visible fringe artifacts.Most artifacts in all three images are successfully removed; we attribute residual artifacts Comparison of convergence speed of our novel algorithm with the alternating optimization one used in [18], [20].
close to the heart, which feature a different spatial frequency, to cardiac motion.In [39], we provide a detailed description of their formation and an algorithm for their reduction.Some of the residual artefacts at the right shoulder are presumably caused by patient motion as well, as this area features large gradients of the transmission.Besides, residual Compton scatter would impair this particular area strongly as well, since the illuminated area of the sample is changing particularly quickly.
There are phase ramps with discontinuities at the boundaries of the G 2 grating's tiles, which are present in the uncorrected as well as in the corrected images.These are caused by a single random vertical offset of the whole interferometer during the sample scan with respect to the reference scans.Because the the vertical gradient of the flat-field phase is discontinuous at the grating boundaries, the effect is most striking there.Removing these artefacts is beyond the scope of this paper.

C. Evaluation of the convergence rate
Here, we demonstrate that our method outperforms alternating minimization of the same loglikelihood [20], [21] in scanning setups, and that the preconditioning described in Sec.II-D.2 is effective.
In order to enable a performance-only comparison with the previous algorithms as originally presented, we simulated and reconstructed only linear displacements in between exposures (φ s 1 = 0, v s = 1).All other parameters of the simulated measurement were the same as before.For our method, we omitted the subsampling described in Sec.II-D as well as the measures to avoid local minima sketched in Sec.II-E.Further, we included an evaluation of the effect of rescaling the pershot parameters in our model.
Since both methods involve a similar computational effort for each iteration, we compared the convergence rate by number of iterations, see Fig. 6.Except for very early iterations, our method converges faster than alternating optimization approaches, and rescaling the problem clearly accelerates convergence.
The algorithm presented in [22] comprises an extension of the one in [20] by rotational vibrations at the expense of replacing the formerly accurate optimization step of the per-shot parameters by an approximation.Thus, we expect a behaviour similar to [20] at best.We are not aware of a quantitative analysis, we analyzed noise power spectra of darkfield images.Ideally, these would be taken from sample-free areas of patient images.However, these are very small due to careful collimination, and often overexposed.Instead, we utilized empty scans.These scans were acquired after each of the patient scans from the study in [9], respectively; here, we take the first 50 of these into account.By separating the whole image into its spatial frequencies, we can focus on the ones associated with vibrations -namely the fringe frequency and its first harmonic [15].These noise power spectra are displayed in Fig. 11.As can be seen, our method works reliably and removes a much larger fraction of fringe artefacts than the one presented in [18].The artefacts' spatial frequency is consistent with our fringe size of about 14 px.

IV. DISCUSSION AND CONCLUSION
We propose an algorithm for correcting artifacts due to grating motion occurring both during and between exposures.We first demonstrated its functionality on simulated data, and then successfully applied it to data from a state-of-the-art scanning Talbot-Lau interferometer.While this kind of setup is generally more prone to vibrations than static setups, it turns out that motion both in between and during exposures can be accounted for.Because the algorithm is formulated as an extension of standard least-squares phase retrieval, implementation is straightforward.These findings may lead to more relaxed stability constraints, and thus, easier construction of future interferometers.While processing becomes much more robust, this approach still has some requirements on setup stability, namely uniform motion of small amplitude during exposures.It has to be noted that while motion-induced bias can be avoided, random grating positions and randomly reduced visibility still result in an increased noise level that is not even spatially homogeneous.
Most recently, an alternative large animal scanning system was presented [40] where the pig is moved through a stationary interferometer.While this avoids vibrations, it suffers from additional artifacts due to sample motion.More importantly, this concept seems impractical for patient imaging in upright position, which is the preferred acquisition setup for chest radiography.
Because our method builds solely upon a stochastic model for the measured intensities, effects that are not included in the model must be corrected for beforehand.In our case, these are motion of the sample and Compton scatter.Compton scatter can be estimated; regions affected by motion can be excluded from the optimization problem.Unlike in this work, the source spot is not fixed with respect to the sample in most scanning setups.In that case, the beam paths (2)-( 4) for the same pixel change slightly during the scanning process, so there is additional inconsistency in the measured intensities.Future research will show to what extent the presented algorithm is compatible with that approximation.

Fig. 2 .
Fig. 2. Setup simulation study.Simulated (a) attenuation, (b) dark-field, (c) differential phase images, generated by a parallel-beam forward projection of an extended FORBILD phantom.The overlay in (a) indicates some of the outlines of the four distinct materials (1) lung, (2) soft tissue, (3) bone marrow, and (4) bone, whose assumed properties are listed in TableI.
This article has been accepted for publication in IEEE Transactions on Medical Imaging.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TMI.2023.3288358This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/

Fig. 3 .
Fig. 3. Simulation study.The flat-field phase of one exposure is shown in (a).Raw measurements of the phantom were simulated, of which (bd) depict three exemplary ones (positions in scan 52, 104, and 156).

Fig. 4 .
Fig. 4. Results of the simulation study.Retrieved attenuation (a, d), dark-field (b, e), and differential phase (c, f) images before (a-c) and after (d-f) estimation of grating vibrations.The grayscales are as in Fig. 2.

Fig. 5 .
Fig.5.Correction of vibration artefacts in images of a patient's thorax (77 year old male).Retrieved attenuation (a-d), dark-field (e-h) and differential phase images before (first column) and after (second column) estimation of grating vibrations.The third and the fourth column show zoom-ins of the uncorrected and corrected images, respectively.Zoom-ins of the dashed regions are shown in Fig.7.For the logarithmic dark-field images, a window from − ln 1.05 to − ln 0.55 was chosen, for the differential phase images, it ranges from -0.2 to 0.2.The zoom-ins feature narrower windows.The scan was performed at 70 kVp, with a total scan duration of 7 s, and an estimated effective dose[33] of 43 µSv.The field of view is 37 cm×37 cm.
Fig.6.Comparison of convergence speed of our novel algorithm with the alternating optimization one used in[18],[20].

Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .
Fig. 7.The effect on the dark-field image of modelling motion during exposures (c) compared with correcting only motion between exposures (b) and no correction at all (a).The region for these zoom-ins is marked in Fig. 5.
This article has been accepted for publication in IEEE Transactions on Medical Imaging.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TMI.2023.3288358This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/ This article has been accepted for publication in IEEE Transactions on Medical Imaging.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TMI.2023.3288358 This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/