Quantitative interior x-ray nanotomography by a hybrid imaging technique

Hierarchical structures appear often in life and materials sciences, and their characterization profits greatly from imaging methods that allow seamless probing of various length scales without sacrificing image quality. X-ray tomography is particularly adept at probing 3D structures; however, zooming in on a region of interest results in a loss of quantitativeness of image contrast and suffers from artifacts unless a priori knowledge or assumptions about the sample are used. Here, we demonstrate a hybrid technique that exploits a micrometer-resolution overview to realize ab initio nanoscale interior tomography with quantitative contrast. In a study of avian eggshell, a model for bionanoporous materials, our approach reveals a complex arrangement of vesicles with sizes ranging from hundred nanometers to a few micrometers. We anticipate that such an approach can be widely adopted and benefited from at synchrotron and laboratory sources, for instance, where such zooming capabilities are already present or can be readily realized. ©2015


INTRODUCTION
In technology as in nature, materials' properties are oftentimes optimized through structures whose characteristic length scales span many orders of magnitude. The need for studying such hierarchical structures has stimulated the development of x-ray tomography techniques with resolutions spanning a wide range of length scales. Sample penetration often necessitates the use of high-energy x rays, for which improved imaging contrast is obtained by techniques that directly use as contrast mechanism the x-ray phase advance upon propagation through the sample [1,2]. For samples on the size scale of tens of centimeters, concepts of Talbot interferometry are applied at synchrotron [3,4] and x-ray laboratory sources [5,6] and are being introduced for medical diagnosis applications [7]. On smaller length scales and finer resolution, Zernike phase contrast [8], holography [9][10][11][12][13], and techniques based on far-field coherent diffraction measurements [14][15][16] are applied. The interaction of x rays with matter frequently allows for a quantitative measurement of electron density, which is valuable for both biological and materials science applications. Furthermore, x rays lend themselves to in situ measurements under varying environmental conditions, and since the measurement is nondestructive, the method allows for multiple complementary studies on the same sample.
Many of these imaging techniques also offer the possibility of increasing resolution by reducing the field of view (FOV), analogous to increasing the magnification of an optical microscope. This capability can be used to obtain higher-resolution tomograms of an interior region in the sample, so-called interior, region-of-interest, or local tomography. While the direct application of standard tomography reconstruction usually yields the correct morphology and location of features, the reconstruction is inherently performed from incomplete data and the quantitative interpretation of contrast is foregone. Significant advances have been made to recover this quantitativeness if the sample can be assumed to be composed of a finite number of piecewise continuous sections or if the electron density for a portion of the interior FOV is known [6,17,18]. Thus, these techniques are able to reconstruct quantitative interior tomograms relying on educated assumptions or prior knowledge of the sample.
While the interior tomogram is often preceded by a coarse lower-resolution measurement used, for example, to identify regions for further investigation, most techniques seek the interior reconstruction independently. However, it has been shown that interior artifacts can be significantly reduced if a prior overview measurement can be brought to bear [19]. Here, we introduce a hybrid imaging and reconstruction technique combining a low-resolution overview with high-resolution projections with the same contrast mechanism and demonstrate that it yields an interior tomogram with nanoscale resolution and quantitative electron density contrast while placing the interior region in context with the entirety of the sample.

SAMPLE PREPARATION AND X-RAY MEASUREMENTS
In order to demonstrate the technique, we imaged a section of avian eggshell from Gallus gallus, with a combination of two scanning x-ray imaging methods. Avian eggshell is a complex structure that consists of biocalcite pillars growing on top of mammillary nodules rooted at the collagen fibers of the shell membrane [20][21][22]. These biomineralized structures grow unidirectionally outward and form the eggshell palisade tissue consisting of calcite interspaced with growth-directing proteins. The eggshell provides robust protection and serves as a calcium reservoir for the developing embryo [23]. At the same time, it allows gas exchange with the atmosphere through pore channels with diameters between 5 and 10 μm that begin between adjacent mammillae and grow towards the eggshell surface, kept open by an unknown mechanism. The eggshell biocalcite contains a large number of ellipsoidal vesicles, with diameters in the range of 100-400 nm, that originate from the uterine granular epithelial tissues [24]. These vesicles are arranged in layers that apparently follow the calcite crystal growth planes, as can be seen in ion-polished cross-sections by scanning electron microscopy backscattered electrons [25], and they play a central role in the deposition of crystal growthdirecting substances [24]. Physically, the insertion of vesicles in the calcite saves calcium imported via the bloodstream, while their arrangement, distribution, and connectivity substantially affect the biocalcite mechanical properties [26]. Studies in two dimensions provide limited insight into the vesicle distribution, as the cutting plane for observation plays a decisive role in the identification of their spatial arrangement. Hence, further study of the 3D morphology and spatial distribution of the palisade vesicles is desirable to elucidate their role in eggshell formation and function.
Sample preparation was carried out at the FOM Institute AMOLF, Amsterdam, The Netherlands. Eggs from chicken, Gallus gallus, were obtained from commercial sources. The shell was obtained after removal of the egg yolk and white, washing with deionized water, and drying at room temperature. An eggshell part of about 2 cm × 2 cm was glued with the outer shell down onto a silica substrate using two-component epoxy glue. After wetting with water the adhering shell membrane was removed with forceps as much as possible. Pillars with a square transverse section of 60 μm × 60 μm were cut in a grid pattern using an automated water-cooled diamond microsaw (Disco Corp. model DAC-2SP/86) along the full thickness of the eggshell such that the 350 μm height of the pillars contained the different layers related to the eggshell structure [20], and the full 360°angular range was available for tomographic measurements. One pillar was manually isolated, picked up with a fine brush, and glued with epoxy to the flattened tip of a steel needle mounted in a brass sample pin of 3 mm diameter with the eggshell's outer surface oriented downwards in contact with the needle.
X-ray experiments were carried out at the cSAXS beamline (X12SA) at the Swiss Light Source, Paul Scherrer Institut, Switzerland, with monochromatic radiation of 6.2 keV photon energy and using a 3D Piezo scanning system installed on top of a rotation stage [27].
The overview data set was measured to 1 μm resolution using x-ray beam deflection tomography [28]. A 1 μm diameter pinhole was used to define a pencil x-ray beam, and the sample was scanned continuously transversely to the direction of beam propagation, the z axis, while diffraction patterns were measured with a Pilatus detector [29] located 7.24 m from the sample. The scan had a field of view of 84 μm × 120 μm, raster scanned in 85 × 121 detector frames with 0.04 s exposure time, and was repeated for 60 different sample orientations covering an angular range of 180°and resulting in a total dose of 72 kGy.
For each scanning point, the x-ray beam deflection was measured from the diffraction patterns with subpixel accuracy using cross-correlation-based registration [30] of the pattern with respect to reference frames taken in air. Knowing the x-ray wavelength and experimental geometry, the beam angular deflection readily translates to the gradient of the phase advance which, measured at different sample orientations and using a modified filtered back projection (FBP) algorithm [28], yields a quantitative map of the 3D distribution of the sample's real part of the index of refraction and subsequently of electron density.
Figures 1(a) and 1(b) show renderings of the beam deflection tomogram. In Fig. 1(b), low density is shown as blue and high density is rendered semitransparent. Two gas exchange pore channels that traverse the eggshell section with diameters of about 4 and 12 μm, respectively, are shown as cyan isosurfaces. In red, we show an isolated cavity with a size of around 9 μm. Within this micrometer-resolution overview, we identified areas of lower density that are due to unresolved submicrometer porosity.
To investigate the submicrometer porous structure in the sample, a volume within the palisades that exhibited a layered structure was chosen for further inspection with highresolution interior tomography, as indicated in Fig. 1

(a) with
Research Article a dashed red line. Interior projections were measured via ptychography in transmission geometry [15,27,31,32], which is a coherent diffractive imaging technique that combines scanning and iterative phase-retrieval algorithms to reconstruct the sample amplitude and phase transmissivity with a resolution much finer than the size of the illumination. A coherent patch of the incident x-ray beam was selected with a 3 μm pinhole, and transverse scans with a FOV of 32 μm × 30 μm were taken in concentric circles with a radial step of 1.4 μm and 0.4 s exposure time per scanning point, totaling 390 detector frames. The measurement was repeated at 450 angular orientations of the sample spanning an angular range of 180°and imparting a total dose of 5.6 MGy for the whole tomogram, which is comparable to previous studies of ptychographic tomography with similar resolution [27,33].

RECONSTRUCTION
Two-dimensional ptychography reconstructions were performed using the difference-map algorithm [31] followed by maximum likelihood refinement [32,34], analyzing a 192 × 192 pixel area of the detector. The reconstructions had a pixel size of 43.8 nm from which the phase was unwrapped [35].
To correct for thermal drifts during the data acquisition, the phase projections were aligned prior to tomographic reconstruction. For alignment in the vertical direction, i.e., along the axis of rotation, we used the approach described in Ref. [36], in which the alignment is obtained by maximizing the correlation of sample mass fluctuations along the axis of rotation, taking advantage of the fact that, for tomographic projections, the vertical mass distribution is independent of the rotation angle if the entire extent of the sample is imaged transverse to the rotation axis. For interior projections, the exact vertical mass distribution is no longer available because only part of the sample is measured transverse to the axis of rotation. However, as noted in Ref. [36], this alignment approach is quite robust to deviations from the ideal scenario due to noise or incomplete measurements if low-order polynomials are removed from the vertical mass distribution. For alignment of the high-resolution interior projections, we removed polynomial functions up to third order.

A. Projection Alignment by Tomographic Consistency
Alignment of projections in the horizontal direction, i.e., transverse to the axis of rotation, was not possible using the center of mass approach [36]. For this task, we developed an iterative sinogram alignment technique based on tomographic consistency, in which we exploit the fact that the inverse Radon transform by FBP followed by a Radon transform only retrieves the original sinograms if they are consistent with a 3D representation, which is not the case for misaligned projections. We can thus use tomographic consistency to robustly refine the sinogram positions and align the projections. For alignment, a tomographic slice in the x; z plane is computed using FBP, and a synthetic sinogram is then computed from this slice using a Radon transform algorithm. The synthetic sinogram is subsequently used to refine the positions of the original projections by minimizing the mean squared error between them, and this process is repeated iteratively. The algorithm proves more robust by starting alignment with low-resolution sinograms, a coarse-to-fine resolution strategy that is often used in image registration [37]. Initially the sinogram resolution is reduced, in this case by a factor of 5, by means of a low-pass Hanning filter. Once convergence is achieved, typically in 10-20 iterations, the process is repeated with sinograms of finer resolution. Convergence can be improved by clipping from the intermediate tomographic slice nonphysical values that are caused by misalignment artifacts, for example by setting to zero negative electron density values. Note that we applied these constraints only during the alignment procedure and not for the final reconstruction. This alignment based on tomographic consistency has proven significantly more robust to noise and artifacts compared to the center of mass computation [36], even when applied to conventional projections, for which the entire sample is measured in the direction perpendicular to the axis of rotation. More details are provided in Supplement 1.
The algorithms for post-processing alignment used in [36], together with the alignment based on tomographic consistency introduced here, make interior ptychographic nanotomography possible with standard FBP methods. However, the quantitative measure of electron density of this technique [38] is foregone due to interior reconstruction artifacts. To exemplify this, Figs. 2(b) and 2(e) show the result of applying the FBP to the interior projections using a smooth sinogram extension, an approach in which the sinogram is padded using the values at its edges [39,40], and which is often used to reduce interior tomography artifacts. Alternatively, the same reduction in artifacts can be obtained from using differential data but the reconstruction remains nonquantitative [41,42]. Although the morphology is obtained with high resolution, an overall density offset is lost and, more importantly, the tomographic sections suffer from artifacts typically characterized by low-frequency variations which complicate threshold segmentation and even prevent relative density measurements between different structures [see Fig. 2(h)].

B. Quantitative Interior Tomogram
Our approach for quantitative reconstruction of the interior region builds conceptually upon previous findings that the interior tomography problem is solvable with very little a priori assumptions about the sample [6,17,18], and on previous simulations combining a low-resolution overview with interior highresolution data to reduce artifacts [19]. Here, we demonstrate that the overview tomogram with close to one order of magnitude worse resolution allows for a well-posed interior reconstruction that retains the quantitativeness of electron density contrast.
The quantitative interior tomogram shown in Figs. 2(c) and 2(f), with 43.8 nm voxel and 130 nm estimated resolution, is obtained from a numerical combination of the overview tomogram with the highly resolved ptychography interior projections. The interior projections were first aligned with respect to the low-resolution tomogram. Afterwards, for reconstruction of each section in the x; z plane, the overview tomogram is used as input to a gradient-based optimization algorithm to find a solution for the sample's real part of the refractive index, δx; y; z, that minimizes the error metric E min a n;m;θ n;m a n;m;θ p n;m x; y i 2 ; (1)  Fig. 1(b) and also seen in Fig. 1(c).
where Φ θ x; y are the ptychographic phase projections [36] with the sample oriented at an angle θ,Φ θ x; y are phase projections computed from the Radon transform of the current estimate of the index δx; y; z, W x; y is a binary mask that denotes the valid area of the ptychography reconstructions, p n;m x; y denotes 2D orthonormal Legendre polynomials of orders n and m, and a n;m;θ denotes the associated coefficients. We introduce p n;m x; y with n; m 0; 0, (0,1), and (1,0) to account for constant and linear offsets of the phase of the object, which are inherent degrees of freedom of ptychography when both the object and the illumination are reconstructed [36]. For conventional ptychographic tomography, these unconstrained parameters are found using the air around the sample as a reference [36], but in the case of interior tomography such a reference is not available. By explicitly minimizing the error metric with respect to the Legendre coefficients a n;m;θ , including constant and linear polynomials in x and y, we make the tomographic reconstruction insensitive to these unconstrained terms. Such an approach is used already in phase retrieval to provide reconstructions insensitive to additive or multiplicative terms in the data [34,43]. For each computation of the error metric, the coefficients a n;m;θ that minimize the error metric can be analytically obtained by a n;m;θ P x;y W x; yΦ θ x; y −Φ θ x; yp n;m x; y P x;y W x; y : (2)

RESULTS AND DISCUSSION
The electron density histogram for a 1 μm section of the sample along the y axis, shown in Fig. 3(d), is dominated by three peaks corresponding to the eggshell matrix (n e 0.79 e∕Å 3 ), sample mounting glue (n e 0.47 e∕Å 3 ), and air within pore channels and vesicles (n e 0.03 e∕Å 3 ). The recovered electron density of empty regions is expected to be zero and therefore provides an estimate of the measurement accuracy. The quantitativeness of the interior reconstruction is confirmed by comparing the obtained density for the calcite 2.62 g∕cm 3 , assuming the chemical composition CaCO 3 , versus its tabulated value of 2.71 g∕cm 3 [44]. The lower value obtained may be explained by errors introduced in the interpolation of the low-resolution tomogram and partial-volume effects of porosity below our imaging resolution. After median filtering with a 3 × 3 × 3 voxel window and threshold-based segmentation of the porous structure the morphology and distribution of vesicles can be quantified. The vesicles were labeled by identifying connected-component regions on a 6-voxel 3D neighborhood. The three orthogonal principal axes of each labeled structure were obtained, and the length of the vesicle was defined as the longest projected length of the vesicle along its principal axes. Figure 3(c) shows a histogram of vesicle size with a mean of 202 nm and standard deviation of 100 nm. The histogram's right tail hints at a secondary peak attributed to agglomerated vesicles, which was confirmed by a fit of the sum of two Gaussian distributions which resulted in peaks at 202 and 509 nm. A threedimensional rendering of the porous structure is shown in Fig. 3(a), where the distribution of vesicles is rendered according to the color code in Fig. 3(c). Figure 3(b) shows a bivariate histogram of vesicle 3D orientation described by spherical polar and azimuthal angles θ and ϕ, with the zenith along the axis of rotation. The individual orientation of the vesicles was obtained based on the direction of its principal axis along which the vesicle has the smallest diameter, as illustrated with arrows in the insets around Fig. 3(b). To reduce outliers from vesicles that are too small to reliably determine orientation, we restricted this analysis to vesicles with a volume larger than 125 voxels, which included about 26,000 of almost 55,000 structures that had been labeled. The mean orientation of vesicles was obtained by averaging their orientation unit vectors, and the spread shown corresponds to the contour lines that include 95% of a bivariate Gaussian probability density function, with the standard deviations and covariance obtained from the data points. The direction of eggshell growth, ϕ −90°and θ 10°, was determined from the overview tomogram. The mean vesicle orientation, ϕ −58.6°and θ 11.8°in the azimuthal and polar directions, respectively, is depicted with a yellow arrow and its spread indicated by an ellipse. Notably, the orientations of individual vesicles are not random, and they exhibit a tendency to align in the direction of eggshell growth.
The layered vesicle arrangement, inferred from the overview tomogram and in agreement with previous observations [20,25], was confirmed with the high-resolution interior tomogram, as is visible in Figs. 1(c) and 3(a). Based on the vesicle position and orientation, these layers were extracted using cluster analysis, and Fig. 3(b) shows their mean orientation and spread with a cyan arrow and ellipse, respectively. The layer orientation is mainly spread along one direction, and two layers in particular have an angle between them of 119.7°. These observations are consistent with the biocalcite crystal geometry and further indicate that the layered structure follows the crystal accretion direction [25]. The effect of the palisades vesicles on eggshell mechanical properties remains largely unknown. In particular, it is widely accepted that palisades vesicles are filled with liquid in fresh specimens, although direct observation of this fluid is complicated by difficulty of staining or loss of the fluid to sample preparation for electron microscopy [20,25]. Because of the quantitative density contrast in our measurement it is clear that the vesicles in our sample are empty. Future measurements with fresh specimens could yield the fraction of vesicles that are filled, which is relevant to assess the effect of this porous structure on the eggshell mechanical properties and may shed some light into statistical morphology and distribution of these vesicles. Such information may be used to model macroscopic mechanical properties, and detailed porosity maps can be used as input for finite-element simulations to elucidate their specific effects on eggshell robustness and crack propagation. Besides fundamental research of avian shell structure and function, such understanding is important for the development of techniques for fast production of robust mineralized materials by biologically inspired processes [45].

CONCLUSIONS
We introduce a hybrid measurement and reconstruction technique that allows tomography on an interior region of interest with quantitative contrast and demonstrate it by revealing the nanoscale porosity on a section of avian eggshell. The procedure described here is ab initio in the sense that it does not require any prior knowledge or assumption about the composition or complexity of the sample. Analogous to the way ptychography [15,31] and Fresnel coherent diffractive imaging (CDI) [16] increased dramatically the applicability and versatility of CDI by removing the need for the sample to be isolated in 2D [14], our approach extends the application of CDI tomography to samples that extend beyond the 3D FOV, thereby simplifying sample preparation, reducing the risk of damaging regions of interest, and expanding its scope of application while preserving the characteristics that make CDI attractive in the first place: quantitative electron density contrast without weak-object assumptions, resolution not limited by lens manufacturing, and nondestructive measurement. Quantitative x-ray phase imaging provides complementary information not easily obtained with other techniques at the nanoscale. For the particular case of ptychography, it finds important applications in the materials and life sciences [27,33,36,[46][47][48][49][50][51][52][53][54][55] and allows responses to different mechanical or environmental conditions on the same sample to be quantified [56,57]. Apart from preserving quantitative electron density contrast, our approach is inherently compatible with multiscale studies of samples down to nanometer-range resolution without the need to dissect the FOV of interest, and profits from high-accuracy dedicated 3D scanning instrumentation with recently demonstrated isotropic 3D resolution down to 16 nm [58].
Although exemplified here using x-ray phase nanotomography, the reconstruction technique is applicable to other quantitative techniques, wavelengths, and a wide gamut of length scales. Talbot interferometry [4] and holographic x-ray tomography [10], for instance, also provide a quantitative measure of electron density and have a natural zooming capability through increased optical magnification to the scintillator or by changing the sample-to-focus distance, respectively.