## Centimeter-sized Objects at Micrometer Resolution: Extending Field-of-View in Wavefront Marker X-ray Phase-Contrast Tomography

DOMINIK JOHN 1,3,4,5,‡,* , JUNAN CHEN 2,3,4,5,‡ , CHRISTOPH GASSNER 3,4,5 , SARA SAVATOVI ´ C 6,3,4,5 , LISA MARIE PETZOLD 3,4,5 , SAMI WIRTENSOHN 1,3,4,5 , MIRKO RIEDEL 1,2,3,4,5 , JÖRG U. HAMMEL 1 , JULIAN MOOSMANN 1 , FELIX BECKMANN 1 , MATTHIAS WIECZOREK 2 , AND JULIA HERZEN 3,4,5

1 Institute of Materials Physics, Helmholtz-Zentrum hereon, Max-Planck-Straße 1, 21502 Geesthacht 2 ImFusion GmbH, Agnes-Pockels-Bogen 1, 80992 München, Germany

4 Chair of Biomedical Physics, Department of Physics, TUM School of Natural Sciences, Technical University of Munich, 85748 Garching, Germany

3 Research Group Biomedical Imaging Physics, Department of Physics, TUM School of Natural Sciences, Technical University of Munich, 85748 Garching, Germany

5 Munich Institute of Biomedical Engineering, Technical University of Munich, 85748 Garching, Germany 6 Department of Physics and Astronomy "Galileo Galilei", University of Padua, 35131 Padua, Italy ‡ The authors contributed equally to this work.

* dominik.john@hereon.de

Abstract: Recent advancements in propagation-based phase-contrast imaging, such as hierarchical imaging, have enabled the visualization of internal structures in large biological specimens and material samples. However, wavefront marker-based techniques, which provide quantitative electron density information, face challenges when imaging larger objects due to stringent beam stability requirements and potential structural changes in objects during longer measurements. Extending the fields-of-view of these methods is crucial for obtaining comparable quantitative results across beamlines and adapting to the smaller beam profiles of fourth-generation synchrotron sources. We introduce a novel technique combining an adapted eigenflat optimization with deformable image registration to address the challenges and enable quantitative high-resolution scans of centimeter-sized objects with micrometre resolution. We demonstrate the potential of the method by obtaining an electron density map of a rat brain sample 15 mm in diameter using speckle-based imaging, despite the limited horizontal field-of-view of 6 mm of the beamline (PETRA III, P05, operated by Hereon at DESY). This showcases the ability of the technique to significantly widen the range of application of wavefront marker-based techniques in both biological and materials science research.

## 1. Introduction

The application of phase-contrast imaging with synchrotron sources has significantly enhanced the exploration of internal structures in materials and biological specimens, offering detailed insights for example into the mechanisms of disease and material structure [1-3]. A range of different approaches have been developed and successfully applied to extract the phase information imparted onto the X-ray wavefront by the sample, including e.g. propagation-based imaging [4,5], ptychographic methods [6], a Fokker-Planck based approach [7,8], and speckle-based imaging [9]. Currently, methods for imaging large objects at the highest possible resolution are an active area of research, because they can provide a better understanding of biological processes and material properties [10]. They are becoming even more relevant with the advent of fourth generation synchrotron sources, which typically feature very small angular source widths [11], leading to smaller fields of view at the same measurement distance. Recently, hierarchical propagation-based imaging has been used successfully to obtain morphological information of different human organs down to the micrometre level [12]. For this purpose, a large number of individual scans is acquired and stitched together to yield information on the complete organ.

To address this problem, approaches like a sample/diffusor drift correction for the case of UMPA [17] and an eigenflat approach have been developed. In the latter approach, a group of acquired flat-field images is decomposed into eigencomponents. A combination of these eigencomponents is then fitted to each projection separately to provide a synthetic flat-field image [18]. This method has also been successfully applied to phase-contrast scans using speckle-based imaging, but currently requires the presence of a sample-free region within the scan to accurately capture the displacement of the wavefront markers [19]. When acquiring multiple tomograms at different positions of the same object, only the tomograms capturing the edges of the object include such a sample-free region. For the inner tomograms, the entire field-of-view is covered by the object and the current eigenflat method is therefore not applicable.

Asimilar approach is currently difficult to replicate for wavefront-marker-based techniques, but would allow retrieving quantitative information on the electron density of the large object while also not requiring the single-material assumption (i.e. a fixed /u1D6FF / /u1D707 -ratio). In these techniques, a wavefront marker is placed in the beam path to create a reference pattern, which is then distorted by the object. Processing approaches such as cross correlation [13] or least-squares minimization (e.g. Unified Modulated Pattern Analysis (UMPA) [14]) then allow retrieving three different signals-attenuation, phase-contrast and dark-field-from multiple tomograms taken at different marker positions. For the purpose of accurate signal retrieval, the beam profile needs to remain constant for both the reference and the sample scan, i.e. the source must remain stable over the time horizon of the measurement-a condition assumed in conventional flat-field correction [15]. For most synchrotron radiation sources, however, the beam profile varies significantly over time (e.g. due to beamline component motion or the electron beam orbit control). These discrepancies between the reference image and the conditions during the measurement may cause the wavefront marker (e.g. a grid pattern [16]) to remain visible in the retrieved projections, compromising the image quality and quantitative accuracy of the imaging method.

An additional challenge for stitching multiple such tomograms are possible deformations experienced by the object due the prolonged exposure to intense radiation. Such deformations may not be accurately described by assuming rigid transformations (translation and rotation), as is typically the case for conventional stitching techniques used for panorama creation [20] or medical diagnosis [21]. For this problem, we therefore propose the use of deformable image registration, a technique applied successfully for example in image-guided surgery [22] and radiation therapy [23].

As an example application, the method is used to retrieve an electron density map of the full horizontal extent of a rat brain 15 mm in diameter at the Hereon-operated beamline P05 at PETRA III (DESY, Hamburg), despite the limited 6 mm horizontal FOV. This scan uses Talbot Array Illuminators as wavefront markers, though the method is also applied towards scans using sandpaper in the supplementary document.

In response to both of the challenges faced by wavefront-marker-based techniques, we present a novel approach combining an eigenflat method developed for this scenario with deformable image registration. In this way, we allow for an arbitrary expansion of the total field-of-view by combining multiple tomograms acquired at different positions to cover the entire sample. This means it becomes possible to obtain electron density maps for objects several centimeters in size with micrometer resolution without strict requirements on the stability of the beam profile or sample over time. Due to the quantitativeness of the method, the results are better comparable over time and at different beamlines.

Fig. 1. Illustration of the image acquisition strategy. The sample rotates around the same axis in every scan, but different parts of the sample are recorded in succession from left to right using the available FOV. While the conventional eigenflat optimization can be performed for the left-most and right-most scans due to the presence of a sample-free region, this is not the case for the inner scan(s). (A) shows the approach for /u1D45B overlapping scans, in which only scan 1 and scan /u1D45B contain a sample-free reference region. (B) illustrates the special case of /u1D45B = 3.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-a230de0c-67bc-4a46-abc3-c4d200976dd2.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-c2950ecf-7951-4eae-bef7-deee45c2ffcd.jpg)

## 2. Methodology

## 2.1. Scan Acquisition

The general scanning strategy for samples that do not fit into a single FOV is illustrated in Figure 1a. The approach is applicable to any number of recorded subregions, but is described for the case of /u1D45B = 3 scans as visualized by Figure 1b in the following.

At first, the left outer, inner, and right outer regions of the sample are scanned separately, resulting in three sets of computed tomography (CT) projections. The iso-centers of these scans are set fixed and the difference between the scans is equivalent to a change in the horizontal offset of the detector and source. The horizontal offsets of the detector and source are chosen so that there is a small overlap between projection images from different sample parts, indicated as left and right overlap respectively in Figure 1b.

## 2.2. Processing

The initial processing is shown in Figure 2a. The outer projection images are acquired such that there are sample-free areas at the borders. Then the eigenflat-field correction as first described by [18] is applied: Given a set of /u1D440 recorded flat-fields /u1D46D = ( /u1D487 1 , . . . , /u1D487 /u1D440 ) with /u1D441 pixels each, the mean of all flat-fields /u1D487 is first subtracted from each to obtain the centered flat-field matrix /u1D468 ∈ R /u1D441 × /u1D440 :

Fig. 2. Illustration of the scan processing pipeline. (A) For the left and right scan data, synthetic reference images are created by matching eigenflats to the sample-free region. For the middle scan, the averaged reference images are used. In all three cases, initial differential phase images are retrieved. (B) The differential phase images are used for deformable image registration, where the right image is matched to the middle image. A B-spline weighted free-form deformation model is defined by the control point grid, which is displayed as the yellow markers in the figure. The same procedure is applied to the left image together with the middle image (not illustrated). (C) A ROI is extracted from the overlapping region both from the middle and the right scan. The weights for the synthetic reference of the middle scan are adjusted iteratively until a satisfactory match between the differential phase images of the middle and right scan is achieved. (D) The calculated weights of the eigenflats for the middle scan are used to retrieve the final corrected middle phase image. This image is stitched with the previously deformed left and right phase images to yield the final projection image. The procedure is repeated for all projections.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-37e23203-b07c-48a6-97e6-ff2e9aa15369.jpg)

$$A = \left ( f _ { 1 } - \overline { f }, \dots,$$

$$\ g m s \text{s uvject to opumazuon} \\ \widehat$$

The eigenflat-fields /u1D496 /u1D458 then refer to the eigenvectors of the covariance matrix /u1D46A = /u1D468 /u1D468 /u1D447 . For each image, a matching flat-field is calculated by least-squares optimization of the difference between pixels in a sample-free region of the image and the corresponding region in the flat-field ̂ /u1D487 /u1D457 , where ˆ /u1D464 /u1D457 /u1D458 are the weights subject to optimization:

The optimized flats /u1D487 /u1D457 are used in the speckle-matching algorithm to obtain the differential phase images in /u1D465 -and /u1D466 -direction. In a first step, the inner projection images of the sample are also input into the speckle-matching algorithm together with an averaged reference image created by computing the average of all originally recorded flat-field images.

The data is then processed using an implementation of UMPA [24], yielding differential phase images in two orthogonal directions for all three parts of the image. The inner projections contain grid-shaped artifacts due to the time-varying nature of the beam, which does not allow to remove the wavefront marker pattern completely using an averaged flat approach. Before an adapted eigenflat correction may be applied to resolve this problem, the three image parts first need to be registered precisely.

Structural changes may be experienced by the samples due to the applied radiation dose. These include thermal expansion or changes in chemical composition due to radiolysis. This, and different beam profiles at different positions within the FOV, mean that the pair of images to be registered may not be perfectly aligned if the images are merely rigidly shifted based on the horizontal offsets. To cope with the potentially complicated deformation that may occur during the course of the scans, a deformable image registration algorithm ( ImFusion GmbH , Munich, Germany [25]) is applied. A free-form deformation (FFD) model [26] with B-Splines interpolation is used. The deformation model is modified by a control point grid. For each pixel inside the grid, its displacement is determined by the displacements of the nearest 4 control points. The resulting deformation field is continuous. Aiming at finding the optimal deformation model, the deformable registration is equivalent to an optimization problem. During the optimization, the displacements of the control points are the variables to be optimized, and the similarity metric between the images is to be maximized.

The deformable registration is demonstrated in Figure 2b. The inner image is assumed to be stationary and left undeformed, while the outer regions receive deformation models. This particular choice creates a computational advantage for the adapted eigenflat-approach introduced in the next subsection. The deformable image registration is then performed on either the differential phase images in the /u1D465 - or the /u1D466 -direction for all projections. We perform the registration consecutively for all the projections. To encourage smoothness between adjacent projections, as well as improve computational efficiency, we use the optimal control point displacements from the projection at the current angle as the initial guess of the projection at the next angle.

## 2.4. Adapted Dynamic Flat-field Correction

In a next step, eigenflats /u1D496 /u1D458 for the inner scan are calculated. However, the matching flat-field images cannot be obtained directly via the approach described in the previous section: The specific combination of weights ˆ /u1D464 /u1D457 /u1D458 (see Eq. (2)) cannot be computed because no sample-free region to be used for matching is present. A different optimization problem targeting the optimal weights for the eigenflats of the inner projection is therefore set up, as shown in Figure 2c.

$$I ( x, y ) \equiv \frac { \partial \Phi } { \partial x$$

In each optimization step, a flat-field image for the inner projection is calculated based on the current guess of the eigenflat weights ˆ /u1D464 /u1D457 /u1D458 and used to retrieve differential phase images /u1D715 Φ /u1D715 /u1D465 ( /u1D465 , /u1D466 ) and /u1D715 Φ /u1D715 /u1D466 ( /u1D465 , /u1D466 ) using the phase-retrieval algorithm (e.g. UMPA), where Φ is the phase of the measured object. Only one of the two signals is required for the optimization to follow, so we define the pixel value as:

for conciseness. We can then use the pixel values /u1D43C inner ( /u1D465 , /u1D466 ) and /u1D43C outer ( /u1D465 , /u1D466 ) from the inner and right outer scans respectively to create a difference image /u1D43C D ( /u1D465 , /u1D466 ) . Inside a region Ω , which we define as the overlap region(s) between the inner and outer scans, /u1D43C D ( /u1D465 , /u1D466 ) is approximated as follows where we have decomposed the difference between the two images into a residual speckle pattern /u1D43C SP ( /u1D465 , /u1D466 ) and a pixel-dependent noise term /u1D441 ( /u1D465 , /u1D466 ) . In an ideal case in which a perfectly matching flat-field has been used for the correction, the term /u1D43C SP ( /u1D465 , /u1D466 ) should equal 0, since the wavefront marker is successfully removed by the speckle-tracking algorithm. Neglecting the noise term, we may therefore use the homogeneity of the difference image /u1D43C D ( /u1D465 , /u1D466 ) as a measurement of the flat-field correction quality. We choose the total variation (TV) for this purpose, since it is a well

$$\begin{array} {$$

established metric for evaluating the spatial variation and the smoothness of an image [27]. For an image /u1D497 it is defined as:

where /u1D465 and /u1D466 denote the pixel indices in horizontal and vertical direction of the 2D image.

$$\text{it is defined as} \colon \\ & \text{TV} ( \nu )$$

The TV of the difference image is now computed and returned as the value of the cost function of the optimization. The cost function is defined as:

which is eventually minimized to find the optimal weights ˆ , /u1D464 as mentioned in Eq. 2.

$$\text{tion. The cost function is defined as} \quad \\ L = T V \left ( I _ { \text{inner} } ($$

Using the determined weights for the eigenflats, new flat-field images are generated and used to retrieve the final, corrected differential phase images of the middle scan. The left outer, inner and right outer differential phase-contrast images are then transformed and deformed using the information of the horizontal offsets and the deformation model previously determined by the deformable registration.

## 2.5. Stitching

Finally, the overlapping regions of the images are blended using a distance weighted sum

$$I _ { \text{blended} } ( x, y ) = \alpha ( x ) \, I _ { \text{inner} }$$

where

$$\alpha ( x ) = \frac { \Delta r - l ( x ) } { \Delta r }. \quad \quad \quad$$

Here, the intensity /u1D43C blended ( /u1D465 , /u1D466 ) of a pixel in the overlapping region is computed from the pixel intensity of the inner image /u1D43C inner ( /u1D465 , /u1D466 ) and the right outer image /u1D43C outer ( /u1D465 , /u1D466 ) . The weight /u1D6FC ( /u1D465 ) depends on width of the overlapping region Δ /u1D45F and the distance of the current pixel to the left boundary of the overlapping region /u1D459 ( /u1D465 ) .

The three projections are stitched using the same method to form one projection with enlarged FOV. The stitched differential phase images in the /u1D465 -and /u1D466 -direction are then integrated using a 2D Fourier integration approach [28] to retrieve the final integrated phase signal, which is used for the three-dimensional reconstruction, yielding an electron density map.

## 2.6. Cascading Approach for Arbitrary Number of Projections

It is to be noted that the proposed method is applicable to an arbitrary number of projections instead of three as shown in Figure 1b. The correction will then be performed consecutively from outer to inner projections. For a projection /u1D45B , where /u1D45B &gt; 2, the overlapping region with previously corrected projection /u1D45B -1 is considered. The cost function for the optimization on projection /u1D45B is defined as where Ω again refers to the overlapping region between the two scans.

$$\det \det & \det & \det & \det & \det & \det & \det & \det & \det & \det &$$

## 3. Experiments

## 3.1. Beamline Parameters

The imaging beamline P05 uses X-rays from an undulator located d = 85 m from the microtomography endstation. A double-crystal silicon monochromator provides a highly monochromatic beam

( Δ /u1D438 /u1D438 ≈ 10 -4 ) with an energy of 30 keV. The beam is approximately parallel (28 µ rad × 4 0 . µ rad at 10 keV [29]). A CMOS camera developed in collaboration with the Karlsruhe Institute of Technology (KIT) based on the CMOSIS CMV 20000 sensor is positioned 620 mm behind the sample. It has an array of 5120 × 3840 pixels with a physical pixel size of 6 4 . µ m, which leads to an effective pixel size of 1 28 . µ musing a 5x magnification lens. The objective is coupled to a 100 µ mCdWO4 scintillator.

The wavefront marker is a 2D phase shifting grating called a Talbot Array Illuminator [16] positioned 175 mm in front of the sample. The grating has a period of 10 µ m and produces a phase shift of Φ = 2 /u1D70B 3 using a hexagonal grating structure. The duty cycle is 1 3 and specifies the ratio of hole spacing to hole radius. The gratings were fabricated on 200 µ mthick silicon wafers using deep reactive ion etching to produce round holes of suitable depth. At fractional Talbot distances of 1 6 /u1D451 T, these TAIs exhibit a focusing effect with a theoretical compression ratio of 1:3 in each direction, allowing high visibility to be achieved. The gratings are stepped using a 2D piezoelectric stepper.

## 3.2. Sample

The sample investigated during this work is a whole rat brain. The animal housing and organ removal was carried out at Helmholtz-Zentrum Munich following the European Union guidelines 2010/63. The procedure was conducted in compliance with the ethical standards of the institution and approved by the responsible governmental body. Specifically, the experiments were performed under the license number ROB-55.2-2532.Vet\_02-21-133. After removal, the sample was fixated in Formaldehyde. Then, the sample was dehydrated using an ethanol series starting with 50%(vol/vol) ethanol per water up to 100%(vol/vol) ethanol in steps of 10% for one hour each. Afterwards, the sample was stained in an eosin stain on ethanol basis, with 0.6%(wt/vol) eosin on 100% ethanol. The ethanol-based staining was used due to the high lipid ratio in the brain tissue. The brain contains a MENX-associated mutation [30].

## 3.3. Measurement Protocol

The rat brain is imaged using a hanging-rotation-axis approach, as illustrated in Figure 3. The sample is placed inside a tube with a diameter of 15 mm and attached to a rotation stage at the top. To minimize variations in electron density, the tube is submerged in a bath filled with ethanol, which prevents the strong changes in phase that would occur if measured against air. This setup also maintains a consistent dynamic range for the phase retrieval method across all scans, which would not be the case if some regions contain air while the inner region is fully occupied by the sample. At the high radiation doses encountered at synchrotron-radiation sources, ethanol is preferred over water because water is prone to radiolysis. This leads to the formation of gas bubbles that can displace tissue during measurements and increase phase variations.

The scans are recorded at an energy of 30 keV using a continuous rotation mode; the position of the wavefront marker is moved after each 180 ◦ rotation, resulting in 9 steps per scanned FOV. At the beginning and the end of each scan a set of 70 flat-field images are taken. In total, 4000 projections are taken with an exposure time of 100 ms each.

## 3.4. Algorithm Parameters

The images are downsampled with a factor of 8 to increase the computational efficiency and mitigate the influence of the speckle background on the registration result. This mitigation is performed since the speckle background inside the overlapping region is not necessarily the same for the two scans and may therefore hinder the registration algorithm from matching the actual sample features. Considering a tradeoff between computational efficiency and deformation complicity, 16 control points on a 4 × 4 grid are assigned at the overlapping region. To emphasize

Fig. 3. Experimental setup used for the rat brain scan. The beam first passes the wavefront marker, a Talbot Array Illuminator, and then imparts on the sample. To allow for the presence of an ethanol bath, a hanging-rotation-axis approach is used in which the sample tube is supported from the top. The window into the bath consists of kapton foil. The passing X-rays are converted to visible light using a scintillator, magnified using a lens and recorded using a CMOS camera (not shown).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-36555443-92f8-456f-b8fc-1ddb8b600a62.jpg)

on the local motion, the local normalized cross-correlation [31] is used as the similarity metric for the registration.

For computational efficiency, only a subregion of the overlapping regions 180 × 240 pixels in size is used for the optimization. Since the target function wraps the phase images that are dynamically retrieved via the phase retrieval optimization, the target function is non-differentiable. We therefore use the derivative-free optimizer BOBYQA [32] from the NLopt library [33]. The optimization variables are the weights of eigenflats, they are therefore unitless. The boundary is empirically set to [-0 5 0 5 , the initial step is 0.1 and the optimization is terminated when . , . ] the relative change of the cost value falls below 1 × 10 -4 . We use the first three components (eigenflats) for the flat-field generation, which results in a total of 27 optimization variables (3 weights × 9 steps).

After the registration, the optimized deformation is applied to the image with the original resolution. The deformed image is then passed to the further steps. The geometrical error induced by beam divergence is neglected given the approximately parallel beam encountered at most synchrotron radiation sources.

## 4. Results and Discussion

Figure 4 shows a subregion of a projection taken during the inner scan. More specifically, the differential phase-contrast signal in the /u1D465 -direction is shown. It is once obtained by a conventional flat-field correction using a mean of the 70 recorded flat-field (Figure 4a) images and by our optimization approach (Figure 4b). While in Figure 4a, remnants of the 2D grating structure are superimposed onto the anatomical information, the effect is strongly reduced in Figure 4b by the correction. The same correction is repeated for all 4000 projection images of the middle part of the scan. The displacements of the FFD control points as determined by the registration are below 10 pixels. This marks the maximum amount by which image pixels at the overlapping region are shifted. The three subregions of the sample are then stitched as previously described

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-761e66ce-0d73-4279-adda-16b73c53e369.jpg)

Fig. 4. Comparison of averaged flat and adapted eigenflat approach for a region in the inner scan and the final stitched projection image. (A) Differential phase-contrast projection corrected with averaged flat-field images. The tissue is overlayed with artifacts in the form of regular grid patterns, which indicates that the speckle grid was at different positions during the flat and sample measurement. For this reason, it cannot be removed correctly using the conventional flat-field correction. (B) The same projection corrected with the eigenflats obtained from our proposed method. The appearance of the artifacts is strongly reduced. (C) Example stitched projection image with enlarged FOV showing the differential phase in the /u1D465 -direction.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-4127f9f5-0cfa-4992-93aa-654ae659a59d.jpg)

After the 2D phase retrieval using UMPA with a window size of 3 and reconstruction of the 2D Fourier integrated phase information using a filtered back projection, an electron density map displaying the entire horizontal extent of the rat brain is obtained. Figure 5 shows an axial slice of the brain. As becomes apparent, no obvious ring artifacts are visible in the scan. These would be created by grating structures in the projections. Moreover, no obvious artifacts due to the stitching of the three scans are present. At the same time, slight changes in electron density throughout the brain become discernible. The measured electron density is similar to the value of 305 previously reported for a human cerebellum [34].

The spatial resolution for the rat brain scan is calculated on a slice of the tomogram using a Fourier Shell Correlation (FRC) method [35, 36] with four image subsets. The resolution is determined to be 12 4 . µ m according to the full-bit and 10 9 . µ m according to the half-bit criterion. This resolution is most likely not the optimum achievable value since the sample was undersampled significantly in this proof-of-principle experiment. Given the approximately 5000 pixels (for 2-fold binning) inside the tube walls, 8000 projection angles would be required according to the Nyquist sampling theorem, but only 4000 projection angles were acquired. A higher number of angles should therefore be chosen for future experiments to obtain the highest possible resolution.

The overlap in this scan encompassed 1183 pixels in the horizontal direction, but may be reduced by approximately a factor of 6 in size in future scans. This is because only a subregion 180 pixels in horizontal extent of the overlap is used for reference-image matching. Such a reduction in overlap width is desirable to prevent redundancy in the acquired data leading to unnecessarily prolonged measurement times.

In the results of the proposed eigenflat optimization shown in Figure 4b, a residual is partially visible as diagonal patterns. This may be attributed to the fact that due to efficiency considerations

Fig. 5. (A) Axial slice displaying the entire horizontal extent (15 mm) of a rat brain inside a tube filled with ethanol. Small differences in electron density between the different tissue types become apparent. Despite being the result of three separate scans with a FOV of 6 mm each, the image gives a uniform impression. (B) Zoom-in region marked by the dotted border in (A).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-947772bd-5af7-4c73-a5cf-fd4c68cb0294.jpg)

we do not use all potentially available eigenflats, meaning not all the potentially available eigenflat-field information is considered for the optimization.

## 5. Conclusion and Outlook

In conclusion, we have presented a novel approach significantly expanding the application range of wavefront-marker-based imaging techniques. Our method enables the retrieval of small differences in electron density for centimetre-sized objects, given sufficient transmission, by cascading beam reference information inward from outer scans and combining subscans using deformable registration. This advancement is particularly relevant for upcoming fourth-generation synchrotron sources, which will feature smaller fields of view, and enables hierarchical scans yielding electron density maps comparable across beamlines, facilitating collaboration and validation of results. The ability to obtain quantitative electron density information for larger samples is a significant step forward in phase-contrast imaging, as it opens up new possibilities for studying complex biological structures and materials at multiple scales. This could be particularly interesting for fields such as biomedical imaging, materials science, and paleontology, where understanding the details of large objects is crucial.

The method was demonstrated using speckle-based imaging but is in principle applicable to any phase-contrast or attenuation-contrast method requiring accurate flat-field information, provided a suitable optimization metric is available. As described in the supplementary document, the method also allows for measurements of even larger objects with Talbot Array Illuminators using the cascading approach. However, it remains to be experimentally investigated whether a large number of stitchings ( /u1D45B &gt;&gt; 5 ) would lead to an accumulation of error when propagating information from the outermost scans inward. In such a case however, one may consider correcting the reference information locally for each FOV instead of using the overlap regions, as demonstrated for the ROI scan in the supplement.

While deformable registration can compensate for small sample shifts, it may not account for large-scale motions or strong movements between wavefront marker steps. Addressing faster changes likely requires a new type of deformation correction. Additionally, since a main advantage of speckle-based imaging is not relying on the single material assumption, future studies should examine large objects that strongly violate this assumption, such as metals in organic tissue [38], to further demonstrate the capabilities of the method.

Another drawback is the computational intensity of the optimization procedure for each pair of overlapping projections. Future work could explore adapting machine learning-based approaches for synthetic reference generation, such as those developed by [37], to alleviate this issue for the case of wavefront-marker based phase-contrast imaging as well.

Funding. Content in the funding section will be generated entirely from details submitted to Prism. Authors may add placeholder text in the manuscript to assess length, but any text added to this section in the manuscript will be replaced during production and will display official funder names along with any grant numbers provided.

Acknowledgments. The authors thank Fabio De Marco for fruitful discussions and Natalia S. Pellegata from Helmholtz Munich for providing the rat brain sample.

Disclosures. Junan Chen and Matthias Wieczorek are employees of ImFusion GmbH. The remaining authors declare no conflicts of interest.

Data availability. Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

- 1. M. Töpperwien, F. van der Meer, C. Stadelmann, and T. Salditt, 'Three-dimensional virtual histology of human cerebellum by x-ray phase-contrast tomography,' Proc. National Acad. Sci. 115 , 6940-6945 (2018).
- 3. K. Taphorn, M. Busse, J. Brantl, et al. , 'X-ray stain localization with near-field ptychographic computed tomography,' Adv. Sci. 9 , 2201723 (2022).
- 2. B. Müller, M. Stiefel, G. Rodgers, et al. , 'Three-dimensional imaging and analysis of annual layers in tree trunk and tooth cementum,' in Bioinspiration, Biomimetics, and Bioreplication XII, vol. 12041 (SPIE, 2022), pp. 98-113.
- 4. D. Paganin, S. C. Mayo, T. E. Gureyev, et al. , 'Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,' J. microscopy 206 , 33-40 (2002).
- 6. F. Pfeiffer, 'X-ray ptychography,' Nat. Photonics 12 , 9-17 (2018).
- 5. J. Reichmann, S. E. Verleden, M. Kühnel, et al. , 'Human lung virtual histology by multi-scale x-ray phase-contrast computed tomography,' Phys. Med. &amp; Biol. 68 , 115014 (2023).
- 7. D. M. Paganin and K. S. Morgan, 'X-ray fokker-planck equation for paraxial imaging,' Sci. reports 9 , 17537 (2019).

elements,' Opt. Express

8.

T. A. Leatham, D. M. Paganin, and K. S. Morgan, 'X-ray phase and dark-field computed tomography without optical

- 9. M.-C. Zdora, 'State of the art of x-ray speckle-based phase-contrast and dark-field imaging,' J. Imaging 4 , 60 (2018).

32

,

4588-4602 (2024).

- 10. M. Du, Z. Di, D. Gürsoy, et al. , 'Upscaling x-ray nanoimaging to macroscopic specimens,' J. applied crystallography 54 , 386-401 (2021).
- 12. C. Walsh, P. Tafforeau, W. Wagner, et al. , 'Imaging intact human organs with local resolution of cellular structures using hierarchical phase-contrast tomography,' Nat. methods 18 , 1532-1541 (2021).
- 11. H. N. Chapman, 'Fourth-generation light sources,' (2023).
- 13. S. Berujon, H. Wang, and K. Sawhney, 'X-ray multimodal imaging using a random-phase object,' Phys. Rev. A 86 , 063813 (2012).
- 15. J. A. Seibert, J. M. Boone, and K. K. Lindfors, 'Flat-field correction technique for digital detectors,' in Medical Imaging 1998: Physics of Medical Imaging, vol. 3336 (SPIE, 1998), pp. 348-354.
- 14. M. C. Zdora, P. Thibault, T. Zhou, et al. , 'X-ray phase-contrast imaging and metrology through unified modulated pattern analysis,' Phys. Rev. Lett. 118 (2017).
- 16. A. Gustschin, M. Riedel, K. Taphorn, et al. , 'High-resolution and sensitivity bi-directional x-ray phase contrast imaging using 2d talbot array illuminators,' Optica 8 , 1588 (2021).

tomography with the unified modulated pattern analysis (umpa) model,' J. Instrum.

17. S. Savatović,

F.

De Marco, M. Riedel, et

al.

,

'Helical sample-stepping for faster speckle-based multi-modal

- 18. V. V. Nieuwenhove, J. D. Beenhouwer, F. D. Carlo, et al. , 'Dynamic intensity normalization using eigen flat fields in x-ray imaging,' Opt. Express 23 , 27975 (2015).

18

, C11020 (2023).

- 19. M. P. Riedel, 'Setup development for high-resolution quantitative phase-contrast imaging at a synchrotron radiation source,' Ph.D. thesis, Technische Universität München (2023).
- 20. C. Herrmann, C. Wang, R. S. Bowen, et al. , 'Robust image stitching with multiple registrations,' in Computer Vision - ECCV 2018, V. Ferrari, M. Hebert, C. Sminchisescu, and Y. Weiss, eds. (Springer International Publishing, Cham, 2018), pp. 53-69.
- 22. P. Risholm, A. J. Golby, and W. Wells, 'Multimodal image registration for preoperative planning and image-guided neurosurgical procedures,' Neurosurg. Clin. North Am. 22 , 197-206 (2011). Functional Imaging.
- 21. S. Adwan, I. Alsaleh, and R. Majed, 'A new approach for image stitching technique using dynamic time warping (dtw) algorithm towards scoliosis x-ray diagnosis,' Measurement 84 , 32-46 (2016).
- 23. K. K. Brock, S. Mutic, T. R. McNutt, et al. , 'Use of image registration and fusion algorithms and techniques in radiotherapy: Report of the aapm radiation therapy committee task group no. 132,' Med. Phys. 44 , e43-e76 (2017).
- 25. 'ImFusion GmbH,' https://www.imfusion.com/ .
- 24. F. De Marco, S. Savatović, R. Smith, et al. , 'High-speed processing of x-ray wavefront marking data with the unified modulated pattern analysis (umpa) model,' Opt. Express 31 , 635-650 (2023).
- 26. D. Rueckert, L. Sonoda, C. Hayes, et al. , 'Nonrigid registration using free-form deformations: application to breast mr images,' IEEE Trans. on Med. Imaging 18 , 712-721 (1999).
- 28. C. Kottler, C. David, F. Pfeiffer, and O. Bunk, 'A two-directional approach for grating based differential phase contrast imaging using hard x-rays,' Opt. Express 15 , 1175-1181 (2007).
- 27. V. Caselles, A. Chambolle, and M. Novaga, Total Variation in Imaging (Springer New York, New York, NY, 2011), pp. 1016-1057.
- 29. F. Wilde, M. Ogurreck, I. Greving, et al. , 'Micro-ct at the imaging beamline p05 at petra iii,' in AIP conference Proceedings, vol. 1741 (AIP Publishing, 2016).
- 31. A. Baig, M. A. Chaudhry, and A. Mahmood, 'Local normalized cross correlation for geo-registration,' in Proceedings of 2012 9th International Bhurban Conference on Applied Sciences &amp; Technology (IBCAST), (2012), pp. 70-74.
- 30. N. S. Pellegata, 'Menx and men4,' Clinics 67 , 13-18 (2012).
- 32. M. J. D. Powell, 'The BOBYQA algorithm for bound constrained optimization without derivatives,' Tech. Rep. NA2009/06, Department of Applied Mathematics and Theoretical Physics, Cambridge University, Cambridge, UK (2009).
- 34. S. Savatović, M.-C. Zdora, F. D. Marco, et al. , 'Multi-resolution x-ray phase-contrast and dark-field tomography of human cerebellum with near-field speckles,' Biomed. Opt. Express 15 , 142-161 (2024).
- 33. S. G. Johnson, 'The NLopt nonlinear-optimization package,' https://github.com/stevengj/nlopt (2007).
- 35. W. O. Saxton and W. Baumeister, 'The correlation averaging of a regularly arranged bacterial cell envelope protein,' J. Microsc. 127 , 127-138 (1982).
- 37. A. Y. Grigorev and A. Buzmakov, 'Flat-field correction of x-ray tomographic images using deep convolutional neural networks,' Bull. Russ. Acad. Sci. Phys. 87 , 604-610 (2023).
- 36. M. van Heel and M. Schatz, 'Fourier shell correlation threshold criteria,' J. Struct. Biol. 151 , 250-262 (2005).
- 38. M. Riedel, K. Taphorn, A. Gustschin, et al. , 'Comparing x-ray phase-contrast imaging using a talbot array illuminator to propagation-based imaging for non-homogeneous biomedical samples,' Sci. Reports 13 (2023).

## Centimeter-sized Objects at Micrometer Resolution: Extending Field-of-View in Wavefront Marker X-ray Phase-Contrast Tomography: Supplementary Document

This supplemental document provides the information on additional applications of the proposed approach. The cascading reference approach is tested for other common scan cases. The performance of the proposed approach is tested on a dataset acquired with a different type of wavefront marker, namely P1000 sandpaper. Moreover, the feasibility of the adapted eigenflat-optimization is tested with region-of-interest scans.

## 1. CASCADING REFERENCE APPROACH

We perform a proof-of-concept experiment to verify the feasibility of the cascading approach (Eq. 9), in which the information from an outer flat-field image is propagated inward to the innermost projection. Note that the experimental measurement contains three scans as described in Figure 1a. Instead of correcting the left outer scan using the conventional approach, we assume it as a scan without sample-free reference area and utilize the reference from the left overlap region. The inner scan is first corrected as shown in Figure 4b. During each optimization step, an ROI from the left overlap region is extracted from the corrected inner phase image as the reference, and the corresponding ROI in the left outer phase image is also cropped out to compute the background, as demonstrated in Eq. 9. A comparison of the left outer phase image using the averaged flat-field and the optimized flat-field is shown in Figure S1. Similarly, the artifact caused by the wavefront marker is reduced after correcting the flat-field with the proposed approach. This demonstrates that the approach is also feasible for 3 inner regions, as would be encountered for n = 5 regions to be stitched.

Fig. S1. Comparison of the averaged flat and the cascading eigenflat approach for a region in the left outer scan. (A) Differential phase-contrast projection corrected with averaged flatfield images. (B) The same projection corrected with the eigenflats obtained from the cascading approach. In the optimization, the overlapping region with a previously corrected projection is used as the reference of the current projection.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-308bcdcc-56f7-47f4-8f67-59a43fec9f24.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-d2e411d9-fec2-4260-a484-f41ff718beb7.jpg)

## 2. SANDPAPER WAVEFRONT MARKER

Since the method is applied to Talbot Array Illuminators in the main section of this work, we want to demonstrate that its working principle is also applicable to sandpaper, a different type of commonly used wavefront marker in speckle-based imaging.

For this purpose, projection data previously acquired for [1] with an effective pixel size of 1.28 µm in a very similar setup is used. The examined object is a phantom sample consisting of polymethylmethacrylate (PMMA), polytetrafluorethylene (PTFE), polychlorotrifluoroethylene (PCTFE), and polyamide nylon-(6,6). It also includes a toothpick and a pipette tip containing polystyrene microspheres.

Fig. S2. Comparison of meaned flat and adapted eigen-flat approach for an ROI in a scan with sandpaper wavefront marker (A) Differential phase-contrast projection corrected with meaned flat-field images. (B) The same projection corrected with the eigen-flats obtained from our proposed method. The appearance of the artifacts is strongly reduced.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-60698125-ee90-4e85-ae42-91763a5adf56.jpg)

The wavefront marker used for this scan is P1000 sandpaper with a mean particle size of 5.8 µm. The same approach as mentioned in section 3 is applied to the acquired data. Figure S2 shows a ROI of the differential phase image retrieved with a meaned flat-field (a) and with an optimized flat-field from our method (b). In Figure S2b, a reduction in background noise is observed, which marks the removal of the random speckle pattern generated by the sandpaper wavefront marker. Moreover, the standard deviation of the differential phase image decreases from 0.81 to 0.56. The result visually and quantitatively confirms the feasibility of the proposed method on data acquired with a sandpaper wavefront marker.

## 3. REGION-OF-INTEREST SCANS

The reference region used in the proposed optimization is not always available in certain cases, e.g. ROI scans. However, a similar strategy is still feasible for the purpose of dynamic flatfield correction. One can adapt the proposed method by evaluating a sample region only. The optimization problem is then turned into minimizing the following cost function:

$$\text{ turned into munizing the following cost function: } \\ L = TV ( \frac { d \Phi } { d x } \Big$$

Figure S3a is the phase image acquired with a meaned flatfield. Figure S3b is the result using the method proposed in [2]. Note that this method is not applicable for ROI scans given that a sample-free reference region is required, so the result is only provided for comparison. S3c is the result from the same method as in Figure S3b, using a sample region instead. Figure S3d is the phase image retrieved using the flatfield generated from the proposed optimization S1. In analogy to Figure S3c, a sample region is used for the optimization. Our result in Figure S3d where Ω denotes a small region of interest inside the sample from the differential phase image. AROI with a relatively low variance is selected in order to mitigate the influence of the sample structures on the TV evaluation. Figure S3 provides a comparison of the result from the proposed eigenflat-optimization and the other approaches. Note that the scan itself is not originally a ROI scan since the edge of the sample is visible in the FOV, but the behavior of the eigenflat method is simulated by considering a ROI inside the sample region. At the same time, a comparison with approaches using the sample-free reference area becomes possible.

Fig. S3. Comparison of the impact on image quality of different approaches to flatfield generation. (A) Projection corrected with a meaned flat-field and without eigenflats. The grid pattern remains clearly visible. (B) Projection corrected with eigenflats (5 components) matched using a sample-free region. This approach would not be feasible for a region-of-interest scan, where no such region exists, but rather serves as a comparison. (C) Same approach as (B) but with a matching region inside the sample. This approach leads to strong artifacts. (D) Method proposed in the main section of the paper with a matching region inside the sample and 3 eigenflat components.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-37a8d380-8fe6-419f-8382-5719aa054087.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00482v1-figure-9aef0b30-98b7-48ee-95b3-18eabc36045a.jpg)

shows a comparable quality to the result in Figure S3b, however, our method does not require a sample-free reference region as in the other method. Moreover, less eigenflat components are needed for our method to generate an image with similar quality (3 vs. 5).

## REFERENCES

- 1. S. Savatovi´ c, F. De Marco, M. Riedel, et al. , 'Helical sample-stepping for faster speckle-based multi-modal tomography with the unified modulated pattern analysis (umpa) model,' J. Instrum. 18 , C11020 (2023).
- 2. M. P. Riedel, 'Setup development for high-resolution quantitative phase-contrast imaging at a synchrotron radiation source,' Ph.D. thesis, Technische Universität München (2023).
Dominik John Junan Chen Christoph Gaßner Sara Savatović Lisa Marie Petzold Sami Wirtensohn Mirko Riedel Jörg U. Hammel Julian Moosmann Felix Beckmann Matthias Wieczorek Julia Herzen
Centimeter-sized Objects at Micrometer Resolution: Extending Field-of-View in Wavefront Marker X-ray Phase-Contrast Tomography
Recent advancements in propagation-based phase-contrast imaging, such as
hierarchical imaging, have enabled the visualization of internal structures in
large biological specimens and material samples. However, wavefront
marker-based techniques, which provide quantitative electron density
information, face challenges when imaging larger objects due to stringent beam
stability requirements and potential structural changes in objects during
longer measurements. Extending the fields-of-view of these methods is crucial
for obtaining comparable quantitative results across beamlines and adapting to
the smaller beam profiles of fourth-generation synchrotron sources. We
introduce a novel technique combining an adapted eigenflat optimization with
deformable image registration to address the challenges and enable quantitative
high-resolution scans of centimeter-sized objects with micrometre resolution.
We demonstrate the potential of the method by obtaining an electron density map
of a rat brain sample 15 mm in diameter using speckle-based imaging, despite
the limited horizontal field-of-view of 6 mm of the beamline (PETRA III, P05,
operated by Hereon at DESY). This showcases the ability of the technique to
significantly widen the range of application of wavefront marker-based
techniques in both biological and materials science research.