Propagation-based phase-contrast tomography of a guinea pig inner ear with cochlear implant using a model-based iterative reconstruction algorithm.

Propagation-based phase-contrast computed tomography has become a valuable tool for visualization of three-dimensional biological samples, due to its high contrast between materials with similar attenuation properties. However, one of the most-widely used phase-retrieval algorithms imposes a homogeneity assumption onto the sample, which leads to artifacts for numerous applications where this assumption is violated. Prominent examples are biological samples with highly-absorbing implants. Using synchrotron radiation, we demonstrate by the example of a guinea pig inner ear with a cochlear implant electrode, how a recently developed model-based iterative algorithm for propagation-based phase-contrast computed tomography yields distinct benefits for such a task. We find that the model-based approach improves the overall image quality, removes the detrimental influence of the implant and accurately visualizes the cochlea.


Introduction
Phase-contrast x-ray imaging (PCI) is a method that is sensitive to the refraction of x-rays in matter, yielding distinct advantages for the visualization of details with similar attenuation properties like those often encountered in biology and medicine. By extending PCI to computed tomography (CT), it has become a widely used imaging method in laboratory and preclinical studies as it allows for high-contrast three-dimensional visualization at high spatial resolution [1].
Several techniques have emerged to exploit and visualize phase contrast [2,3]. The experimentally most intuitive technique is referred to as propagation-based phase-contrast imaging. It solely relies on a sufficiently coherent x-ray wavefield transversing a sample and propagating sufficiently far to the detector. Thereby, the phase shifts induced by the sample lead to distinct interference effects in the measured intensities [4,5]. From these interference effects, the phase-shifting properties of the sample can be reconstructed from a single image if certain assumptions are imposed onto the sample [6][7][8].
The two most common assumptions are either that the absorption is constant and thus can be neglected [7] or that absorption and phase are coupled, which is usually referred to as the homogeneity assumption and results in the single-material phase-retrieval algorithm of Paganin et al. [8]. However, for the presented sample, the attenuation has a crucial influence and cannot be neglected, therefore we restrict ourselves to the second approach.
Due to the benefits of PCI and the ease of implementation, propagation-based PCI methods have yielded distinct benefits for numerous applications [1]. One of the most prominent applications is the investigation of biological samples [9]. However, additional materials that violate the homogeneity assumption degrade the reconstruction quality significantly, which remains the most considerable drawback of this imaging technique. Different approaches have been developed to alleviate this problem: A generalization of Paganin's algorithm is able to accurately reconstruct known materials that are embedded in another medium. This method requires multiple reconstructions, which are obtained for each pair of materials sharing an interface, to be spliced together manually [10][11][12]. In terms of the streak artifacts, which play a predominant role in the presented application, this generalization would thus not change the overall behaviour and is therefore not applied in this case. Another approach relies on multiple images being acquired at multiple distances for every angle. This enables the reconstruction of heterogeneous samples, but makes the experimental setup more cumbersome [13,14].
In the following, we discuss the reconstruction of a guinea pig cochlea with implant as depicted in Fig. 1(a) as an example of the examination of a biological sample with foreign objects that violate the homogeneity assumption. We show how recently developed model-based iterative reconstruction algorithms for this imaging technique [15] can improve image quality of such tasks and account for the detrimental influence of the implant compared to the conventional two-step reconstruction approach consisting of a phase retrieval prior to tomographic reconstruction [6].
The study of animal cochleae with implants is particularly interesting because the insertion of the electrode array into the narrow turns of the cochlea often causes subtle damage, which can result in wide variations in experimental results. Sometimes, the electrode touches or penetrates the basilar membrane and disrupts the organ of Corti, which lies on top. In chronic experiments, trauma caused by inserting the electrode often triggers soft tissue growth. In all these cases it would be extremely valuable to assess the state of the implanted cochleae as precisely as possible. However, metal parts of the electrodes usually cause large imaging artifacts. These  artifacts mask the delicate structures of the inner ear, especially close to the electrode, where damage is most likely to occur and where the highest resolution would be required. In theory, the cochlear implant electrode could be retracted from the cochlea, but this procedure is likely to cause additional damage. Therefore, in-situ imaging of the cochlea with the inserted implant is highly desirable [16][17][18].

Methods
In the following, the underlying assumptions around the image formation process are briefly addressed. Next, the conventional two-step reconstruction approach is described. Afterwards, the model-based iterative reconstruction algorithm is formulated. Fig. 2 illustrates the different approaches.

Image formation and assumptions
In general, for small propagation distances, the phase and attenuation properties of the sample cannot be recovered independently from a single measurement. Therefore, we employ the widely used homogeneity assumption [8], which couples the intensity to the phase of the wavefield. If the complex refractive index of the sample is known, the measured x-ray wavefield behind the sample can provide the sample thickness, which is usually referred to as the trace of the sample.
The propagation of the x-ray wavefield exiting the sample to the detector is modelled here by the transport-of-intensity equation (TIE) [19]. In addition, one assumes small propagation distances, such that the intensity evolves linearly with the propagation distance [8,20].

Conventional two-step reconstruction
The conventional approach as illustrated by the blue arrows in Fig. 2 consists of two steps. The first step is the single-material phase-retrieval filter by Paganin et al. (PAG) [8], which recovers the trace t of the objects from the flat-field corrected intensities y under every angle independently. It is derived from TIE [8] with the assumptions described previously. This phase-retrieval step can be written as   [8]. The volume is then reconstructed using the filtered back-projection (FBP) algorithm. In comparison, we use a model-based iterative reconstruction algorithm (MBIR) [15], which recovers the volume directly from the measured intensities. Thereby, the whole image formation is modeled, including the statistical properties of the x-rays, attenuation, phase-shifts, propagation and prior knowledge about the sample.
where the propagation distance is given by z, the linear attenuation coefficient of the material is denoted by µ, and δ is the real part of the deviation of the complex refractive index from unity. In addition, F ⊥ denotes the two-dimensional Fourier transform along the spatial dimensions perpendicular to the propagation distance and k ⊥ are the corresponding Fourier frequencies. This filter is a regularized low-pass filter, which removes the edge-enhancement effects on the measurements and increases the contrast [8]. Due to the design of the filter, noise is reduced as it is predominantly present in the high frequency components.
The second step of the conventional approach is a filtered back-projection (FBP) of the recovered traces, which can be denoted as and yields the three-dimensional density distribution x of the sample. Thereby, the conversion from intensities to line integrals is already performed within the phase-retrieval step.

Model-based iterative reconstruction
Various model-based iterative reconstruction algorithms have been reported for conventional absorption tomography that directly recover the three-dimensional density distribution of the linear attenuation coefficients from the measured intensities, yielding distinct advantages in terms of image quality [21,22]. In these, one models the entire image formation process and finds the optimal solution based on a statistical optimization problem. Built upon the success of these methods, recently this approach was extended to propagation-based computed tomography with the assumptions presented above [15]. Fig. 2 illustrates in orange how this approach aligns with the previously discussed two-step approach. Thereby, the forward model, which estimates the expected intensitiesȳ directly from the three-dimensional density distribution of the sample x, can be written as where D denotes a diagonal matrix with the vector in brackets on its diagonal. The exponential function is to be understood element-wise. The matrix A denotes the forward projection operation and the matrix L denotes the Laplacian operator using a five-point stencil finite-difference method. Consequently, the first part of this forward model, denoted by the diagonal matrix, accounts for the attenuation of the x-rays inside the sample and coincides with the models for conventional absorption imaging. The second part in brackets models the interference effects due to the phase-shifts induced by the sample, which is the dominant contribution for the contrast in propagation-based PCI.
To find the optimal density distributionx, which coincides best with the acquired measurements y according to their statistical properties as well as prior knowledge about the sample, an objective function is established and optimized according tô where the first summand denotes the data-fidelity term and the second summand denotes the regularization term. The data-fidelity term introduces two weighting functions. The first term K y = D[y] is the covariance matrix [23], which accounts for the reliability of the acquired measurements with respect to noise and arises from the generation of x-rays following a Poisson distribution and individual x-rays not being correlated. Its inverse is denoted by K −1 y . The second term D[w] is a measure for the reliability of the established forward model excluding therefore rays, where the homogeneity assumption is violated. Lastly, the regularization term favors smooth solutions while preserving edges inside the volume. It consists of a coupling parameter β, which has to be chosen accordingly and a Huber potential R γ given by where γ is the transition parameter between the quadratic and linear penalty [24]. In this formulation the indices refer to specific entries of the vectors, thereby i runs over all voxels and n over the next neighbors of i. The neighborhood N i therefore holds 26 voxels. The different distances to the neighbors are taken into account by the additional factors ∆ i,n ∈ {1, 2 1/2 , 3 1/2 }. The objective function given by Eq. (4) cannot be solved analytically and a full representation of the operators (especially the projection matrices) is not feasible due to memory limitations. Therefore, gradient-based iterative optimization routines are employed. For the current estimate of x, the gradient of the objective function with respect to x is calculated to provide information about the subsequent search direction. Employing a line-search routine that provides the step size along the search direction, the new estimate of x is defined and the procedure is repeated iteratively.

Preparation, scan and reconstruction setup
To prepare the cochlea for the measurement, the sample was boiled out and bleached with hydrogen peroxide. Afterwards, the implant was deployed into the sample. The custom made guinea pig electrodes from MED-EL GmbH are made from Platinum and the wire material consists of Platinum-Iridium embedded in Silicon. The implant is with the exception of its size identical to cochlear implant electrodes used in patients.
The measurements of the cochlea sample were acquired at the micro-tomography end-station of the imaging beam-line P05 at PETRA III at DESY (Hamburg, Germany) operated by the Helmholtz-Zentrum Geesthacht. The experimental setup is described in [25]. The provided CCD camera system with 12.0 µm pixel size was used with a magnification factor of approximately 4.96 for the visible light optics, giving an effective pixel size of around 2.42 µm. The energy of the incoming monochromatic photons was set to 35 keV. Due to the limited field-of-view of 7.2 mm in horizontal direction, an off-axis tomography was performed. As the illumination in the vertical direction is limited to approximately 2 mm, two measurements were performed at different height levels. The propagation distance between the sample and the camera was set to be 1 m. The exposure time was set to 0.1 s.
Due to fluctuations of the monochromator, the intensity distribution of the incident beam changed over time. For each projection, the corresponding flat-field image for the empty-beam correction was therefore chosen by minimizing the variance to all available flat-fields in a sample-free region. To correct for a slight rotation of the detector around the beam-axis with respect to the tomographic axis, a cross-correlation of overlapping data in the off-axis tomography was performed. The information obtained by this procedure was also used for calculating the number of overlapping pixels for the stitching procedure which allowed to create a virtual half-scan of the whole sample. The stitching of the two tomographies at different height levels was performed by maximizing the variance of sum of the overlapping regions. Finally, the data was cropped and binned by a factor of four. The resulting dimensions were 1434 × 450 pixels for the 1500 virtual projections with a pixel size of 9.68 µm. An example of such a preprocessed virtual projection image is depicted in Fig. 1(b).
When selecting the qualitative parameters for sharp phase retrieval, we set the ratio δ/µ = 10 −9 , which defines the strength of the phase retrieval filter. To avoid any scaling, we thereby set µ = 1 m −1 .
All algorithms were executed on a heterogeneous compute server with two Intel Xeon CPUs (E5-2690 v4), 1 TB of DDR3 RAM and four Nvidia Titan Xp (12 GB GDDR5).

Results and discussion
b a 500 µm 500 µm For the conventional reconstruction approach, the traces were recovered with the single-material algorithm in Eq. (1). As an example, the recovered trace of the corrected intensity measurement shown in Fig. 1(b) is depicted in Fig. 1(c). The subsequent tomographic reconstruction was performed with the filtered back-projection algorithm according to Eq. (2). A part of a reconstructed slice is depicted in Fig. 3(a), with its position marked by a red rectangle in the overview slice at the upper right corner. One finds strong streak-like artifacts arising from the implant. In addition, both the single-material phase-retrieval as well as the analytic tomographic reconstruction introduce additional blurring, which reduces the spatial resolution significantly. One possible answer to the loss in resolution are post-processing techniques. For instance, qualitative information from a reconstruction without phase-retrieval can be manually merged to improve the resolution at the edges again [26].
The model-based iterative reconstruction approach omits the intermediate step of recovering the traces from the intensity measurements. This step is directly performed as part of the tomographic reconstruction inside a statistical optimization problem given by Eq. (4). Since our model in Eq. (3) does not describe the metal implant, we chose the weights D[w] such that pixels in the projection images, where the metal implant is present, do not contribute to the reconstruction by simple thresholding. For this, the projection images are directly reconstructed by a filtered backprojection. The thresholding is then performed on the volume and the obtained binary mask is dilated to also account for the edge-enhancement. Subsequently, the binary weights are obtained by selecting all values which are zero after forward projecting the binary mask of the volume. The reason for applying the thresholding to the volume and not directly on the projection measurements is that this approach is less prone to noise and the mask is then tomographically consistent. The transition parameter γ = 5 of the regularizer was estimated from the background noise of the conventional reconstruction approach. The strength of the regularization was set to β = 5 · 10 −5 . As an optimization routine, an unmodified L-BFGS algorithm [27] was used with 2000 iterations ensuring converged estimates. As the conventional reconstruction is corrupted by streak artifacts, it is not suitable as an initial guess for the MBIR algorithm. Therefore, we perform additional 250 iterations of the MBIR algorithm initialized with zeros and a very strong regularization parameter of β = 10 −3 to obtain the low-frequency components. The same reconstructed slice as for the conventional approach is depicted in Fig. 3(b). Due to the masking of the weights, the implant is removed in the tomographic reconstruction, therefore omitting any artifacts from the implant. In addition to the weighting of remaining data points according to their noise properties, the regularization suppresses noise while maintaining the resolution at the edges. Finally, by directly integrating the propagation of the x-rays into the tomographic reconstruction, the phase retrieval of different projections is coupled by the regularization in the volume and weighted according to the noise properties.
In Fig. 4, renderings created by using Avizo Fire 8.1 (Thermo Fisher Scientific, USA) from the two reconstruction approaches are depicted. As the implant gives strong attenuation contrast, as seen on the measured intensities in Fig. 1(b), its position within the volume can be segmented by thresholding a conventional attenuation-based FBP reconstruction. This information is depicted in red in Fig. 4. Due to the streak artifacts arising from the implant in the conventional reconstruction approach, the cochlea itself cannot be accurately recovered. Furthermore, the position of the implant within the cochlea cannot be visualized precisely as seen in Fig. 4(a) and indicated by the insert and the blue arrows. In contrast, when taking the model-based reconstruction as a basis for the rendering of the cochlea, an artifact-free rendering can be obtained and the position of the implant can be accurately depicted within the volume, as seen in Fig. 4(b). Especially, the regions around the cochlear implant can be accurately visualized.

Conclusion
We examined an animal cochlea with a strongly absorbing implant, using propagation-based phase-contrast tomography at a synchrotron facility, discussing the challenges arising from such implants that violate the homogeneity assumption. Thereby, we compared the conventional two-step reconstruction approach to a novel model-based iterative reconstruction algorithm, which is applied to synchrotron radiation for the first time. The image quality was improved overall in terms of resolution and the artifacts (streaks etc.) arising from the metal implant were addressed and removed by the iterative approach. For the presented cochlea sample, this is for instance beneficial to evaluate whether the electrode array has penetrated the basilar membrane.
Our approach can be transferred to similar imaging tasks that involve biological samples with  additional features that violate the homogeneity assumption, such as metal implants or bones. Thereby, the high phase contrast between materials with similar attenuation properties can be utilized without being diminished in the vicinity of highly absorbing features. The metal streak artifacts examined in this study are predominantly related to the violation of the homogeneity assumption and thus are not present in a direct reconstruction of the measured intensities without performing any phase retrieval. Also, in purely attenuation-based tomography, there are numerous sources that lead to similar metal streak artifacts, which potentially influence the artifacts. The most prominent underlying effects are beam hardening, noise from low photon counts (photon starvation), scattering, motion artifacts and nonlinear partial volume effects (NLPV, also known as exponential edge-gradient effects). Accounting for these effects is especially crucial in medical applications [28][29][30][31]. Due to the monochromaticity of the source, beam hardening should be negligible, although some contamination with higher harmonics can occur. With sufficient transmission seen even behind the implants, the artifacts are also not attributed to beam starvation, although the presented approach could naturally account for beam starvation in the same way it accounts for the violation of the homogeneity assumption [15]. The relatively large sample-to-detector distance in propagation-based imaging should also reduce the effect of scattering sufficiently. The inanimate sample should also not introduce motion artifacts. Lastly, NLPV can indeed influence the image quality for our application. Traditional model-based iterative reconstruction approaches may alleviate these artifacts by a finer sampling of the image volume, but this increases computational cost. However, recent model-based iterative reconstruction algorithms have been developed that tackle this problem by including prior knowledge of the implant into the tomographic reconstruction [31]. Closely related are under-sampling artifacts (aliasing), if the Nyquist theorem for tomography [32] N projections ≥ π 2 × N columns (6) is not fulfilled. Due to regularization techniques, iterative reconstruction can however account for this sparsity extremely well [33,34]. An overview of additional sources for artifacts in synchrotron microtomography, not necessarily restricted to metal streak artifacts, can be found in [35].
In general, there are several components of our MBIR algorithm that improve the image quality compared to the conventional reconstruction approach. Firstly, the introduction of the weights D[w], which we introduced in this paper as a measure for the reliability of our forward model, enables us to exclude all data that is not described by our forward model. Secondly, apart from being able to handle different types of artifacts, iterative methods in general, which employ a statistical noise model and regularization techniques, have proven to increase the overall image quality in terms of resolution and noise and especially regarding dose requirements [36,37]. Although the weights D[w] are the crucial components for the removal of the streak artifacts arising otherwise from the violation of the homogeneity assumption of the forward model, the regularization alone would still mitigate the influence of these artifacts.
With current computer hardware, iterative approaches can be applied even to large data sets. There are two main components that increase the total computational time. Firstly, the most time-consuming operations are the projection operations. In each iteration, two projections are performed: one forward projection in the forward model given by Eq. (3) and one backprojection in the evaluation of the gradient of the objective function given by Eq. (4) with respect to x. However, using state-of-the-art graphic cards, individual projection and backprojection operations for this data set took only 3 − 5 s. Implementation details of the projection operations can be found in [38]. The projection operations and thus the computational time scales roughly linearly with the number of voxels times the number of acquired projections. Secondly, the total computational time is directly proportional to the number of iterations. Solving the objective function in Eq. (4) in fewer iterations is challenging, due to the non-linear forward model of Eq. (3). However, some progress has recently been made in this respect for the optimization of non-linear models for attenuation-based computed tomography [22]. Furthermore, if a sufficiently accurate analytical reconstruction is available as a starting guess, the number of iterations and thus the computational time decreases significantly and the MBIR approach can be thought of as a refinement step. The versatility of model-based approaches allows for arbitrary extensions [37]. As an example, models for the source or the detector may be introduced in the forward model [22] or more advanced regularization techniques might be used, which have proven beneficial for multi-distance propagation-based PCI applications [14].