Iterative phase contrast CT reconstruction with novel tomographic operator and data-driven prior.

Breast cancer remains the most prevalent malignancy in women in many countries around the world, thus calling for better imaging technologies to improve screening and diagnosis. Grating interferometry (GI)-based phase contrast X-ray CT is a promising technique which could make the transition to clinical practice and improve breast cancer diagnosis by combining the high three-dimensional resolution of conventional CT with higher soft-tissue contrast. Unfortunately though, obtaining high-quality images is challenging. Grating fabrication defects and photon starvation lead to high noise amplitudes in the measured data. Moreover, the highly ill-conditioned differential nature of the GI-CT forward operator renders the inversion from corrupted data even more cumbersome. In this paper, we propose a novel regularized iterative reconstruction algorithm with an improved tomographic operator and a powerful data-driven regularizer to tackle this challenging inverse problem. Our algorithm combines the L-BFGS optimization scheme with a data-driven prior parameterized by a deep neural network. Importantly, we propose a novel regularization strategy to ensure that the trained network is non-expansive, which is critical for the convergence and stability analysis we provide. We empirically show that the proposed method achieves high quality images, both on simulated data as well as on real measurements.

Introduction

Despite ever-improving diagnostics and therapies, breast cancer remains the most common malignancy in women [1]. To fight this public health burden, better imaging techniques are needed to improve early detection of the disease and increase survival rates. In fact, none of the currently used breast imaging techniques (mammography, breast ultrasound, breast MRI and absorption-based breast CT and tomosynthesis [2, 3]) is able to deliver fully three-dimensional images with sufficiently high isotropic resolution and soft-tissue contrast necessary to identify critical breast cancer imaging biomarkers [4]. X-ray phase contrast CT could potentially offer a solution by combining higher soft tissue contrast with the high three-dimensional resolution which characterizes CT [5].

When X-ray waves interact with matter, their amplitude and phase are modified according to the refractive index of the material they interact with. In particular, the refractive index of a given material is described as n = 1 − δ + iβ. The real part of the index of refraction δ determines the change in the wavefront’s phase Φ as where l indicates the direction of the X-ray beam. From this, the refraction angle α can then be calculated using where λ is the X-ray’s wavelength. The imaginary part of the index of refraction β determines the attenuation coefficient μ via μ = 4πβ/λ, which can then be used to calculate the beam’s attenuation by using the Beer-Lambert law.

It is well known that different soft tissues have similar β’s [5], which leads to little contrast between different tissue types in conventional CT. On the contrary, these differences are larger in δ’s [5]. Consequently, this can theoretically lead to higher soft tissue contrast in the reconstructed CT volumes [6].

X-ray detectors cannot directly measure the X-ray’s phase. Many techniques have been proposed over the years to solve this problem: propagation-based [7], crystal interferometry [8], analyzer-based [9], edge-illumination [10] and grating interferometry [11-13]. The latter is a promising approach for in-vivo imaging since it satisfies the prerequisites for clinical compatibility: it has non-restrictive requirements in terms of temporal and spatial coherence of the X-ray beam, it can be operated at large fields-of-view (FOV) and it has a comparably high mechanical robustness [11].

Talbot-Lau grating interferometry detects the X-ray’s refraction angle α by exploiting a particular interference pattern called Talbot carpet [14]. When an X-ray beam is refracted, this results in a lateral shift in the interference pattern. Therefore, by measuring this shift, the wavefront’s change in phase can be easily obtained by integrating Eq (2). To obtain the Talbot carpet, three gratings are placed between the source and the detector [12]. The first grating (source grating or G0) is made of highly absorbing material and is placed right in front of the X-ray source to increase the beam’s coherence. The second grating (phase grating or G1), which is not designed to absorb photons, imposes a significant phase shift to the X-ray beam and creates the interference pattern. Typically, the detector has a resolution in the range of hundreds of micrometers and cannot resolve the pattern of a few micrometer period. Therefore, a highly absorbing third grating (analyzer grating or G2) is placed in front of the detector. By moving one of the gratings with respect to the others in the x-direction (see (2)), it is then possible to obtain an interferogram called the phase stepping curve [15], from which it is in turn possible to compute the lateral shift of the interference pattern.

The phase stepping curve can be modeled as

The transmission signal T, the dark-field signal D, and the differential phase signal are obtained with:

Here, g2 is the pitch of the G2 grating and d2 denotes the distance between the origin and G2. I0, V0 and Φ0 are the corresponding flat-field intensity, visibility and phase maps, respectively, and k is the k-th phase step.

By combining the phase stepping curves of a flat-field measurement with the actual scan, it is possible to obtain the interference pattern’s shift φ, i.e. the differential phase contrast (DPC) signal, with simple a Fourier analysis. While with grating interferometry (GI) it is also possible to obtain the absorption signal μ, which is related to the average intensity of the curve, and the dark-field signal ϵ, which is related to the curve’s amplitude, in this paper we focus exclusively on the phase δ. If combined with a CT acquisition protocol, GI naturally extends to GI-CT. Note that since grating interferometry detects the first derivative of the phase (see (6)), during reconstruction an integration step is required if phase contrast tomograms are to be reconstructed.

In an attempt to translate this promising technology into clinical practice, our group has embarked in a long term effort to build a first-of-its-kind GI breast CT (GI-BCT) prototype. As such, our scanning protocol must respect clinical constraints in terms of deposited radiation dose, scanning times and patient comfort.

It is very challenging to obtain high quality images with GI-CT when operated in low dose regimes and with fast acquisition protocols. In fact, to date, the successful use of GI phase contrast CT has been limited to synchrotron beamlines [16] and laboratory setups [17, 18] where high image quality is achieved by a high X-ray flux or by long scanning times, respectively.

Since conventional absorption-based CT has been used in the clinics for decades now, a natural question that arises is: “Why is it more challenging to perform fast low-dose GI-based phase contrast CT as compared to fast low-dose absorption contrast CT?”

This question has two main answers, both of which are related to the fact that the phase cannot be measured directly. First of all, there is an intrinsic noise amplification during the signal retrieval of the differential phase in the projection domain [19]; secondly because the differential nature of the phase contrast forward operator causes the inverse problem to be more severely ill-conditioned as compared to conventional CT.

The amount by which the noise gets amplified during signal retrieval depends on the quality of the grating interferometry setup. The dominant characteristic here is the visibility, i.e. the amplitude of the interference pattern in a flat-field scan. The higher the visibility, the more precisely one can compute the interference pattern’s lateral shift and, consequently, the X-ray beam’s refraction. High visibilities are notoriously challenging to achieve and the problem of noise amplification must thus be addressed by using a powerful denoising (or regularization) strategy, which is the first challenge this article aims to address. We would like to highlight that while for attenuation denoising is a relatively easy task because of the nature of its noise power spectrum (NPS), for phase it is more challenging since the noise increases with decreasing spatial frequency [6].

The second challenge is more difficult to tackle since the ill-conditioning of the problem is intrinsically coupled to the physics of the signal acquisition in grating interferometry. Nonetheless, the implementation of the forward-backward operator pair has an impact on the stability of the reconstruction process. This article proposes new operators which are less ill-conditioned and thus less prone to instabilities as compared to the mainstream implementation [20].

In light of the aforementioned two challenges, it is important to reconstruct the tomograms with a stable and powerful inversion algorithm. A pseudo-inverse such as the filtered backprojection (FPB) algorithm could be applied in conjunction with the Hilbert filter [21]. However, it is widely accepted in the CT community that iterative reconstruction algorithms are better suited to deal with highly ill-conditioned problems.

In iterative reconstruction a variational loss function is typically defined which comprises a data-fidelity term and a regularization functional that incorporates prior knowledge about the expected reconstruction. Minimization of this loss with an optimizer of choice then allows to reconstruct an image which is simultaneously consistent with the measurements as well as with the prior knowledge.

A widely used prior in image reconstruction is the total variation (TV) prior [22] which promotes sparse edges in the solution by assuming piece-wise constant signals. While TV is still considered to be a powerful baseline algorithm to regularize ill-conditioned inverse problems, recent years have witnessed the rise of data-driven algorithms which outperform traditional methods.

In addition to pure black-box models, two main data-driven approaches that draw inspiration from classical variational optimization schemes have been proposed. The first approach comprises end-to-end methods which unroll iterative schemes, thereby transforming each iteration of the iterative reconstruction algorithm into a distinct layer of a neural network [23, 24]. Since the imaging physics is embedded into the network, these models are generally believed to be more robust to noise and adversarial perturbations compared to pure black-box neural networks and more data efficient too. In the second approach, the idea is to learn the regularizer a-priori on a representative training set, and to then use the trained regularizer in conjunction with the data-fidelity term in a classical variational optimization framework [25-31].

The first type of approach tends to yield superior results and delivers faster reconstructions [24]. However, it has the significant disadvantage that no convergence guarantees can be derived and that it needs explicit supervision. Algorithms in the second category tend to be slower and potentially yield slightly inferior results. However, they have several advantages. The algorithms can be trained independently of the forward operator, they are more data-efficient (and possibly unsupervised), and they are amenable to stability and convergence analysis [27]. Since it is of utmost importance in the medical field to reliably reconstruct the tomograms, we propose a novel algorithm that fits into the second category.

Contributions

To address the two major challenges in GI-BCT outlined above and thereby be able to reconstruct phase contrast tomograms from noisy measurement data, we propose an iterative reconstruction algorithm which 1) uses a less ill-conditioned forward-backward tomographic operator pair, and 2) leverages the power of deep learning to regularize the (still) highly ill-conditioned tomographic inversion problem. In particular, we propose an algorithm that alternates between data updates governed by the L-BFGS algorithm [32] and regularization steps performed with a deep neural network in a Plug-and-Play fashion [33]. Furthermore, we propose an efficient spectral norm regularization strategy which enforces the network to be globally non-expansive.

We apply the proposed approach to both simulated data and real measurements and show that it achieves excellent results. Importantly, we show that it is possible to train the regularizer on simulated data and later apply it to real measurements. This is an extremely useful characteristic since training data is always scare in a medical context.

We would like to mention that while to the best of our knowledge this paper is the first one to combine phase contrast CT reconstruction with a data-driven prior, many articles have been published in the past on iterative phase contrast CT reconstruction [34-37]. Likewise, many works have investigated the synergy of deep learning and iterative methods. Romano et al. introduced the concept of regularization by denoising (RED) in which a denoiser is used in alternation with data updates [25]. This concept is similar to the idea of Plug-and-Play regularization [33] and has then later been implemented with neural networks by [28]. A different approach has been taken by Lunz et al. who proposed a data-driven adversarial regularization approach in which the prior is trained on unpaired data to distinguish between “good”, artefact- and noise-free images, and “bad” images containing artefacts and/or noise [26], while Mukherjee et al. later proposed a convex variant of the former approach [27]. On the other hand, also generative approaches towards learning a prior have been investigated. Two notable ones include score-based models in which learn the the gradient of the log-likelihood of the distribution [31], and generative adversarial networks which learn the prior distribution by letting a discriminator and a generator network compete with each other [38].

Our work adds an important contribution to the field since, to the best of our knowledge, it is the first one to combine the quasi-Newton optimizer L-BFGS with a data-driven non-expansive Plug-and-Play prior. Our method is especially useful for highly ill-conditioned inverse problems where a quasi-Newton method greatly accelerates convergence compared to simpler first-order methods.

Materials and methods

GI-BCT prototype and scanned sample

Our GI-BCT prototype is a 5th Talbot order 3-gratings symmetric interferometer with a pitch of 4.2 μm. A Comet MXR-225HP/11 tube operated at 70kVp and 10mA was used together with a CdTe photon-counting Dectris prototype detector with a pixel size of 75 μm operated at 10Hz.

To demonstrate the effectivenss of the proposed method on real data, we scanned a formalin-fixed mastectomy specimen from an autopsy without any gross pathological changes obtained at the Department of Pathology and Molecular Pathology, University Hospital Zürich. The data has been obtained with written informed consent under the ethical approval KEK-2012_554 obtained from the Cantonal Ethics Commission of canton Zürich. We acquired 600 projections under continuous circular rotation. After a full rotation, we phase-stepped G0. We combined 5 rotations to construct the sinograms sampling the phase stepping curve. A total of 10 scans have been averaged to partially compensate for the low visibility (thus high noise amplitudes) we are currently working with. Future improvements in grating fabrication are expected to make the latter step obsolete.

Simulated data

We used in-silico breast phantoms as described in [39]. DPC sinograms have then been simulated by forward projecting clean phantom data, adding realistic Poisson noise and finally retrieving the differential phase contrast sinogram with Fourier analysis [39]. We simulated 120000 photons leaving the source to get as close as possible to the best data quality we can achieve on our prototype. These DPC sinograms then served a double purpose: 1) to generate training and validation data to train the denoising network, and 2) to quantify the reconstruction performance of the proposed data-driven iterative algorithm.

To generate training and validation data for the denoising network, we ran a few iterations of the (unregularized) optimization scheme by starting from the clean phantom as the initial guess. Thereby, realistic phase wrapping artifacts and noise are introduced in the images, as shown in Fig 1. The clean phantoms and the generated images were then used for supervised training of the network in a 2D slice-by-slice manner. A total of 560 images were used for training and validation.

Fig 1

Training data for the denoising network.

On the left: a noisy iterate obtained after applying 15 gradient updates based on noisy sinogram data; in the middle: a clean breast phantom; on the right: the difference image.

DPC forward and backward tomographic operators

The optimization problem we aim to solve is where δ is the tomographic volume of the real part of the index of refraction belonging to the set of clean phase contrast CT slices, φ is the retrieved DPC sinogram and A is a linear forward operator modeling (6). We would like to note that our loss function assumes a constant variance within the retrieved DPC data. While it is possible to compute the expected variance in the data based on the noise model developed in [19], this model requires the true attenuation and dark-field signals as input. In particular, the variance depends quadratically on the (noiseless) dark-field signal induced by the sample. Since it is difficult to accurately estimate the dark-field signal from highly noisy projection data, this uncertainty in the noise model input would lead to large errors in the DPC variance estimation. Therefore, we decided not to incorporate this into our reconstruction pipeline.

In order to iteratively optimize (7) with a gradient-based method, we need to model the forward operator A as well as its adjoint, the backward operator A⊤. As shown in (6), the forward operator in grating interferometry’s phase contrast channel includes a differential term. While this sets it apart from attenuation-based CT, it remains a linear operator, which is an important prerequisite for analyzing convergence of the proposed algorithm as it ensures convexity of the variational loss.

Two main approaches have been proposed to compute this derivative. In both cases, the differentiation step in the direction perpendicular to the gratings follows the computation of the line integral. The approaches differ in how they compute this derivative. One approach is to numerically approximate it, e.g. by computing a finite difference approximation [40] or by using splines [37]. The other approach parameterizes the image pixels with spherically symmetric expansion functions which allow to analytically obtain the derivative [35]. By computing this derivative a-priori and storing it in a lookup table, it is then possible to sample from this function in the detector plane, thereby circumventing the need for numerical differentiation.

The second approach is more appealing since analytical differentiation is known to be more stable than numerical differentiation, and thus a higher robustness to noisy measurement data can be expected. However, the method proposed in [35] has the disadvantage that it uses a voxel-driven forward model. These operators are known to be non-parallel, and thus less scalable since it is not possible to compute every detector pixel value independently. Ray-driven forward operators on the other hand, are parallel and are thus better suited for large-scale imaging problems. For the inverse operator the opposite is true: voxel-driven methods are easier to parallelize than ray-driven methods.

In this paper, we adapt the method proposed in [35] in order to obtain a ray-driven forward operator and a voxel-driven backward operator, thereby achieving higher scalability. Both operators have been implemented in CUDA for 3D cone beam geometry.

Forward operator

What sets the proposed operator apart from the one in [35] is that we sample from the analytical derivative before computing the line integrals. In other words, we first calculate how much each pixel refracts the X-rays, before summing up the individual contributions to obtain the total refraction. Besides offering advantages in terms of computational speed and problem conditioning, the proposed implementation more closely resembles the physical process of refraction compared to a somewhat abstract differentiation step in the projection plane.

The rationale of our forward operator is displayed in Fig 2 and described in the following. Let be our phase contrast image and our sinogram with α projections, m detector columns and n detector rows. Before applying the forward operator we compute the analytical derivative of the following spherically symmetric expansion function parameterized by the Kaiser-Bessel window [35]: with where we chose α = 3 and w twice the image resolution as this empirically led to superior results. Once the lookup table is filled, we first compute the rays connecting the X-ray source S to detector element D. Then, we determine the minimum distance d between the ray and the voxel centers p. Based on d we then sample the lateral shift imposed by voxel p from the derivative of the spherically symmetric expansion function stored in a lookup table. Finally, a sum over all non-zero r yields the detector count φ[α, m, n]. Given the highly parallel nature of this approach, the sinogram φ is obtained by computing all its entries in parallel.

Fig 2

Ray-driven forward operator for phase contrast CT.

Please note that only a 2D view is provided for simplicity. Image voxel centers p are displayed as orange dots, the ray as a green line and the distances d of the voxel centers p to the ray are shown in grey. Note that only those pixels are considered which lie within a distance of w (the width of the spherically symmetric expansion function) from ray . The spherically symmetric expansion function is plotted as a dotted line, its derivative as a continuous line.

Note that the derivative information is considered solely in the horizontal direction. In the vertical direction we used classical linear interpolation to increase computational speed.

Backward operator

For the backward operator we follow the reverse logic. Here, every image voxel can be computed independently to yield an highly parallel operator. The concept is illustrated in Fig 3 and explained below.

Fig 3

Voxel-driven backward operator for phase contrast CT.

Please note that only a 2D view is provided for simplicity. One image voxel center p is displayed as an orange dot, the line as a green line. Detector centers D are plotted as black dots. The intersection q of the line with the detector plane is shown as a blue dot. The spherically symmetric expansion function is plotted as a dotted line, its derivative as a continuous line.

First, we compute the lines connecting the X-ray source S and the voxel centers p. We then calculate the intersection of line with the detector plane, thereby obtaining q. Computing the distance between detector element centers D and q, we obtain t, which allows us to sample from the lookup table to get the lateral shift

Note that, for cone beam geometry, the lookup table used in the backward operator, which is acting in the detector plane, must be scaled to account for the geometric magnification. Finally, a summation over all angles and detector elements yields the voxel δ[i, j, k].

In order to save computations, both for the forward and the backward operator, only those voxel centers, respectively detector pixel centers, are considered which lie within a distance w from the traced rays.

Problem (ill-)conditioning

In Fig 4 we show the spectrum of the singular values of the normal matrix A⊤A, i.e. the Hessian of (7). We observe that the spectrum of the proposed operator decays at a slower rate compared to the finite difference-based operator. Especially for the larger singular values, our operator is superior. This is somewhat expected since the image has been parameterized by circular expansion functions covering more than one ‘classical pixel’, which is the parameterization used by the finite difference operator. Therefore, each element in the solution is influenced by a larger number of elements, which increases the curvature of the Hessian in its principal directions, and makes the proposed operator more robust to noise.

Fig 4

Spectrum of singular values of the normal matrix .

The spectrum of the proposed operator is shown as a dotted line, the spectrum for the finite difference-based operator is displayed as a continuous line. The spectrum of the latter method decreases more rapidly compared to the spectrum of the proposed method. The singular values have been normalized.

Unfortunately, (7) remains a highly ill-conditioned problem even with the new operator. In other words, given the quickly decaying spectrum in Fig 4, the loss in (7) has a very heterogeneous curvature, with some directions being steep, while others being almost flat. The latter makes it difficult to optimize (7) and leads to slow convergence. Moreover, analytically solving (7) yields

Given that A⊤A has some very small singular values (), a multiplication with its inverse will yield a highly unstable solution. Since an unregularized optimization of (7) will converge to this fixed point, a powerful regularization strategy is necessary to stabilize the reconstruction.

Data-driven regularizer

As mentioned in the introduction, data-driven regularizers are pushing aside classical regularization strategies such as total variation (TV) [22]. We here propose to learn a regularization network which is able to remove both noise and artifacts from the image iterates in a Plug-and-Play fashion as they converge to the final reconstruction.

Given noisy and clean images as shown in Fig 1, we aim to train a biasless network with ReLU-nonlinearities R, which maps the corrupted image x to its clean counterpart y, because for a biasless net there exists an input x = 0 s.t. f(x) = 0, and hence 0 lies in . This fact will be used to theoretically motivate our algorithm in the Theoretical motivation for PnP-BFGS section. Furthermore, since in the Theoretical motivation for PnP-BFGS section we consider our network to be a projection operator, we aim at ensuring the non-expansiveness of f, i.e. f should be Lipschitz continuous with Lipschitz constant L ≤ 1 i.e.

Since explicitly constraining the Lipschitz constant of f to be upper bounded by 1 is an NP-hard problem [41], we penalize the product of the Lipschitz constants of the layers W of the network if it is larger than 1. To achieve this we propose the following regularization term: with ϵ = 10−8, and where L(W) denotes the Lipschitz constant of layer l. Since a composition of L-Lipschitz functions yields at most an L-Lipschitz function, i.e. any minimizer of the regularizer ensures that the network f is non-expansive.

The Lipschitz constant of a layer l can be obtained by computing the spectral norm of W. Since the spectral norm of W corresponds to the largest singular value σ of W, we can compute it using the power method [42].

Our training loss then combines the classical denoising loss with the proposed regularizer : where N is the number of training samples, the clean image δ* is sampled from the manifold , and the noisy image δ is simulated as explained above. We set λ = 10−4 and empirically, a (quasi-)minimum of the regularization term, i.e. , is always reached during optimization. Note that, while the regularizer leads to a network which is 1-Lipschitz, the individual layers may not be non-expansive (see Fig 5). In fact, especially the first layer has a fairly high Lipschitz constant. We would like to mention that allowing some layers to be expansive led to significantly improved denoising performance compared to explicitly constraining every layer to be non-expansive. This seems to suggest that it is important for the network to have the freedom to learn some expansive layers. On the other hand, we observed that without spectral norm regularization the product of the spectral norms of the layers diverged.

Fig 5

Lipschitz constants of the network’s layers.

The product of the layers’ constants is 0.99. This sets an upper bound on the overall network’s Lipschitz constant, thus making it non-expansive.

We parameterized f with a 7-million-parameter U-net [43]. The model has been trained with Adam [44] until the validation performance stopped improving. Importantly, we removed all biases because it has been shown that 1) this leads to higher model interpretability, and 2) that biasless networks generalize better to different noise amplitudes [45].

The higher interpretability comes from the fact that, in the absence of biases, we can regard the denoising process as being locally linear [45]. Therefore, the Jacobian of the trained network shows how the pixel neighborhood is used to denoise a particular pixel (see Fig 6).

Fig 6

Model interpretability.

Due to the local linearity of a bias-free network, it is possible to visualize how the denoising happens by computing the Jacobian of the network with respect to the input image [45]. A small image patch is shown for computational reasons. From left to right, we plot the noisy input patch, its denoised counterpart with two pixels highlighted in red. The third and fourth image show how strongly each pixel is weighted to denoise the two pixels highlighted in red.

The high robustness against different amounts of noise is important in our case since 1) the noise encountered during image reconstruction might slightly differ from iteration to iteration, and 2) it makes it easier to train the network as we do not have to find the perfect amount of noise to train on.

Crucially, before inputting the training images into the network, we randomly moved the mean of both the noisy and clean images. This ensures that the network learns to remove artifacts and noise independently of the mean of the image. This is essential in our case because, given the differential nature of the problem, pixels far away from an edge take longer to converge compared to pixels close to an edge. Therefore the local mean of the iterates varies within the images and changes over iterations.

Iterative phase contrast CT reconstruction with data-driven regularizer

To optimize (7), we propose an alternating optimization scheme which alternates between data updates and denoising steps performed by the data-driven non-expansive denoiser. In particular, the proposed algorithm performs k steps to optimize the data-fidelity term in (7), followed by a denoising step. Given the high ill-conditioned nature of our optimization problem, and the inability of first-order methods to effectively optimize (7), we chose to use the L-BFGS quasi-Newton method [32] to perform the data updates, because of its ability to approximate the Hessian of the loss function, while being very memory efficient.

By removing noise and artifacts, the denoising step with a non-expansive network can be interpreted as an (approximate) projection of the current image iterate x onto the data manifold of clean images. Iteratively alternating between data-fidelity optimization and projection to the data manifold thereby allows to approach the measured data, while simultaneously staying within the manifold, and thus yielding a physically relevant solution.

The algorithm starts with an all-zero image and stops when the loss is sufficiently small, i.e. when where ϵ is given by the amount of noise in our data. Since the L-BFGS optimizer restarts with an empty inverse Hessian estimate after each denoising step, our algorithm can be interpreted as a meta algorithm which iteratively applies the L-BFGS method with ever improving starting points.

Both the L-BFGS optimizer and the denoising network have been implemented in Tensorflow 2.0 [46]. The forward and backward operators implemented in CUDA have been registered as Tensorflow functions. The complete algorithm is summarized in Algorithm 1.

Algorithm 1: L-BFGS optimization with data-driven Plug-and-Play denoiser

input: i = 0; δ0 = 0; kmax = 15; ∊;

while do

k = 0; y = 0; s = 0;

if i > 0 then

δ = δreg;

end

while k < kmax do

;

;

if k > 0 then

s[k − 1] = δ − δ;

;

end

;

k = k + 1;

end

i = i + 1;

;

end

output: δreg

Theoretical motivation for PnP-BFGS

The following motivating analysis of the proposed reconstruction algorithm holds for BFGS, i.e. the non-memory efficient version of L-BFGS. The key reason for choosing the L-BFGS algorithm with a PnP denoiser is its potential to recover the underlying true image with fewer iterations as compared to first-order methods (e.g., gradient-descent (GD) with PnP denoising). In this section, we formally show that any quasi-Newton (QN) method that builds a progressively accurate approximation of the inverse Hessian can potentially achieve a better convergence rate, thereby leading to significant saving in computations.

Our objective is to recover the image by solving (7). To motivate the use of PnP-BFGS, we consider the following algorithm that alternates between updating the current iterate along a descent direction , followed by a projection step onto the image manifold :

δ = P(z).

Here, g = ∇f(δ) = A⊤(Aδ − φ) is the gradient of the quadratic loss in (7) evaluated at the kth iterate δ, and B denotes a positive-definite matrix that is updated recursively using a suitable QN strategy (which is L-BFGS in our case). To make the analysis simpler, we make the assumption that the data is noise-free, i.e., φ = Aδ*. With this simplifying assumption, we have g = A⊤A(δ − δ*) = H(δ − δ*), where H is used as a shorthand for the Hessian matrix A⊤A.

We make the assumption that the image manifold is a closed and non-empty set containing the zero image, i.e., . The pre-trained PnP denoiser f can be thought of as an (approximate) projection onto , since it projects noisy images (that are outside ) onto . Since the convolutional layers in f have no bias term, f(0) = 0, which is consistent with the assumption .

The cone at the true solution δ* is defined as

Suppose, the null-space of H satisfies , i.e., contains no vector in ker(H) apart from the zero vector. This is akin to the standard set-restricted eigenvalue condition considered in [47].

Under the assumptions made above, if the QN algorithm builds a progressively accurate estimate of the left inverse of H over , then the proposed PnP-BFGS algorithm can potentially achieve a super-linear convergence rate. More precisely, a super-linear rate of convergence can be achieved if the following holds, as we argue later in this section:

At this point, we state the following facts (see [47] for proofs):

Fact-1: If is a closed set, .

Fact-2: Let be a closed cone and let . Then, it holds that , where is the unit ball in .

Fact-3: Let be a closed and nonempty set that contains 0. Let be a nonempty and closed cone containing , i.e., . Then for all , it holds that , where κ = 1 if is convex and κ = 2 if is non-convex. (Lemma 6.4 in [47]).

Using the results above, the error in the (k+ 1)th iteration can be bounded as where the second equality is a consequence of Fact-1. Now, using Fact-2 and Fact-3 above and setting η = 1, we have

Now, dividing both sides of (19) by ‖x − x*‖2 and taking the limit as k → ∞, we notice that as a consequence of (17). Therefore, under reasonable assumptions on the image manifold , (20) indicates that PnP-BFGS can achieve super-linear convergence, leading to significantly faster convergence that GD. Nevertheless, it remains to be shown rigorously that (17) indeed holds for L-BFGS-based updates, and, therefore, the analysis presented here only serves as a motivation behind choosing PnP-LBFGS for reconstruction.

Results

To benchmark the performance of the proposed method, we compared it to three alternative methods: 1) a classical iterative phase contrast reconstruction algorithm which uses a finite difference-based numerical differentiation strategy for the forward and backward operators along with the well-known TV-based regularization [20]; 2) a filtered backprojection; and 3) a filtered backprojection followed by a deep learning-based post-processing step. For the comparison to be fair, we used the same network for the post-processing as we used within the iterative reconstruction. Please note the latter comparison is easily applicable only when an analytical solution exists (and thus not applicable to cases like [48]). If not, two problems might arise: 1) it could be computationally demanding to create a dataset for training a post-processing net; and 2) in highly ill-posed scenarios, an iterative algorithm might not even be able to converge to a sufficiently high-quality image, which would consequently make the denoising nearly impossible.

To quantitatively evaluate the performance of our model, we reconstructed a 32 × 1280 × 1280 voxel volume based on a 32 × 600 × 1280 pixel noisy sinogram randomly selected from our simulated dataset. As seen in the upper image of Fig 7, this data is corrupted by both heteroscedastic noise (see red arrows) as well as by phase wrapping artefacts (see yellow arrows). Both factors need to be compensated for by the regularization functional.

Fig 7

Simulated sinogram (top) and real sinogram (bottom) used for reconstructing the images in Figs 9 and 10.

The scale bar measures 25 mm. The yellow arrows indicate phase wrapping artefacts, the red arrows indicate localized high amplitude noise.

To validate the performance of our algorithm on real measurements we scanned a formalin-fixed mastectomy which yielded the DPC sinogram in the lower part of Fig 7. We reconstructed a 32 × 1280 × 1280 voxel volume based on a 32 × 600 × 1280 pixel sinogram. This data is again corrupted by the heteroscedastic high-amplitude noise highlighted by the red arrows. Phase wrapping artifacts are less pronounced here since there are no straight edges in the scanned mastectomy (see Fig 10).

The computation time for a single data update step of the proposed method was 39s on average. The regularization step performed by the trained denoiser had a negligible computational overhead of roughly 100μs. Until convergence the proposed algorithm needed 12 iterations, each composed of 15 data updates combined with one regularization step. Therefore, the whole reconstruction took 117 minutes, i.e. roughly 2 hours. We would like to mention that further optimization of the CUDA implementation of the forward and backward operators will likely enable faster reconstruction times. The baseline algorithm based on finite-difference-based operators and TV regularization had a total reconstruction time of 175 minutes. Here, the long reconstruction time is caused by the regularization step performed with the Chambolle algorithm [49], which is significantly slower compared to a deep learning-based denoiser. On the contrary, the finite-difference-based operator was ca. 8 times faster compared to the proposed operator, which is to be expected since in the proposed operator the voxels are parameterized by basis functions which span multiple classical cubic voxels.

Effect of the proposed operators

To compare the effect of the proposed operator and the finite difference-based one within the iterative reconstruction scheme, we performed a simple experiment. By starting from the clean ground truth image, we performed a single gradient update based on noisy sinogram data. Since this operation adds noise to the image by propagating it from the sinogram space to the image space, it can be visualized as the displacement of the image iterate away from the manifold of clean images . Since the task of the denoiser is to project the iterates back to the very same manifold, it is evident that a smaller displacement will result in an easier task for the denoiser.

The result in Fig 8 shows that the proposed operator introduces less noise after one gradient step compared to the finite difference-based operator. While this may be hard to capture by eye, the gradient update of the proposed operator leads to a peak-signal-to-noise-ratio (PSNR) of 11.08 dB, while the baseline method achieves only 9.60 dB. Since these reconstruction algorithms are iterative in nature, this effect gets amplified as the algorithm proceeds, i.e. as the image iterate moves farther away from the manifold .

Fig 8

The effect of one gradient update step.

The proposed operator is shown on the left, the finite difference-based operator in the middle. The clean image used as a starting point is shown for reference (right).

We can thus confirm that the proposed operator has an important advantage over the baseline operator since it eases the task of the denoising network, thereby enabling the algorithm to handle larger amounts of noise.

Effect of the proposed data-driven regularizer

Simulated data

Now that we have shown the advantages of the new operator, we will turn our attention to the benefits of using a data-driven prior over a classical one. Fig 9 shows the reconstruction results of simulated data we obtained with the proposed algorithm, along with the comparisons to the TV-based iterative reconstruction [20], filtered backprojection (FBP) and deep learning-based post-processing. We can see that the proposed algorithm achieves excellent results and significantly outperforms the TV-based reconstruction, both qualitatively as well as in terms of PSNR and SSIM (see Table 1). Moreover, we observe that the proposed method is better at removing phase wrapping artifacts compared to TV-based reconstruction, where some persist at the phantom edge.

Fig 9

Reconstruction results on simulated data.

First row: clean phantom; second row: proposed algorithm; third row: TV-regularized algorithm; fourth row: filtered backprojection. On the left a full slice is shown and on the right a zoomed-in section given by the red rectangle is displayed. The white scale bar is 15 mm. The two small red rectangles have been used to compute SNR and CNR values.

Table 1

Quantitative results for the reconstructions in Figs 9 and 10.

Simulated data	PSNR	SSIM	SNR	CNR
Proposed algorithm	30.50	0.96	41.12	3.90
TV	24.73	0.91	45.35	5.71
FBP	23.05	0.77	22.30	2.09
DL post-processing	35.32	0.97	52.57	4.42
Real data	PSNR	SSIM	SNR	CNR
Proposed algorithm	-	-	94.80	18.38
TV	-	-	51.78	8.70
FBP	-	-	26.78	3.51
DL post-processing	-	-	77.84	9.71

A closer look at the TV-based reconstruction reveals the typical staircase artefacts. Interestingly, TV performs better in terms of SNR and CNR. This is attributable to the fact that the latter method promotes piecewise constant signals which translate into low standard deviations and thus higher SNR and CNR values.

A further observation that should be made is that the difference in performance between the TV-based reconstruction and the filtered backprojection is not very large. We believe this is because of the smoothing effect of the Hilbert filter, which suppresses some of the noise, and which is not used when iteratively reconstructing the data.

Finally, a comparison to the deep learning-based post-processing reveals that the latter achieves higher image quality metrics as well as slightly superior visual image quality. This is to be expected since in this case the network has been trained to solve exactly the task it is evaluated on, and it acts on pseudo-inverse solutions of the inverse problem, i.e. FBP with the Hilbert filter. On the contrary, within the iterative reconstruction algorithm, the network acts on intermediate iterates and was not trained to maximize the performance of post-processing, i.e. the network does not get applied to a pseudo-inverse solution of the inverse problem.

Real data

We applied exactly the same algorithms as for the simulated data. Importantly, we used the neural networks that have been trained on simulated data, both for the proposed algorithm as well as for the post-processing.

A look at Fig 10 reveals that the proposed algorithm offers superior results compared to TV-based reconstruction also on real measurements, both qualitatively, as well as in terms of SNR and CNR (see Table 1). PSNR and SSIM values could not be calculated since no ground truth image was available. The fact that the denoising network generalizes well to real data, even though it was trained solely on simulations, indicates that our simulations match well with the real data.

Fig 10

Reconstruction results on real data.

First row: proposed algorithm; second row: TV-regularized algorithm; third row: filtered backprojection. On the left a full slice is shown and on the right a zoomed-in section given by the red rectangle is displayed. The white scale bar is 15 mm. The two small red rectangles have been used to compute SNR and CNR values.

The results in Table 1 reveal that, on real data, the post-processing approach achieves an inferior performance compared to the proposed algorithm both in terms of SNR and CNR. A closer look at the images in Fig 10 reveals that there are some artefacts in the post-processed image (see red arrow) which are not present in the proposed method. Given the superior performance of post-processing on simulated data, these results suggest that the proposed method might offer higher generalizability. This phenomenon could be explained by the fact that the latter is truly a hybrid approach that takes into consideration the imaging physics and the image prior in a principled way (inspired by quasi-Newton algorithms), whereas a post-processing method does not incorporate the imaging physics.

Discussion

In this work we have proposed an iterative reconstruction algorithm with novel tomographic operators and a data-driven regularization strategy, along with an efficient strategy to train a non-expansive network.

We could show that the spectrum of the proposed operators, in which the derivative is computed analytically, decays less rapidly compared to the finite difference-based operators, in which the derivative is computed numerically, thus facilitating data-fidelity optimization. Moreover, the new operators propagate less noise, thus easing the task of the denoising network.

The data-driven regularization strategy proved to be significantly superior to TV-based regularization and perhaps more robust in terms of generalization compared to a deep learning-based post-processing step. Importantly, the training loss based on spectral norm regularization allows to efficiently learn a non-expansive network. This approach is applicable to many different fields where a bound on the Lipschitz constant of a network is desired.

In conclusion, we proposed an algorithm which can handle larger amounts of noise compared to the baseline algorithm. Given the noisy data acquired on our setup, the results in this article suggest that a data-driven prior might be indispensable for reconstructing high-quality phase contrast CT images.

27 Jun 2022

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Jude Hemanth

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Please amend your Methods section to include the ethics approval number included in your Ethics Statement.

3. In your Data Availability statement, you have not specified where the minimal data set underlying the results described in your manuscript can be found. PLOS defines a study's minimal data set as the underlying data used to reach the conclusions drawn in the manuscript and any additional data required to replicate the reported study findings in their entirety. All PLOS journals require that the minimal data set be made fully available. For more information about our data policy, please see http://journals.plos.org/plosone/s/data-availability.

Upon re-submitting your revised manuscript, please upload your study’s minimal underlying data set as either Supporting Information files or to a stable, public repository and include the relevant URLs, DOIs, or accession numbers within your revised cover letter. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories. Any potentially identifying patient information must be fully anonymized.

Important: If there are ethical or legal restrictions to sharing your data publicly, please explain these restrictions in detail. Please see our guidelines for more information on what we consider unacceptable restrictions to publicly sharing data: http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions. Note that it is not acceptable for the authors to be the sole named individuals responsible for ensuring data access.

We will update your Data Availability statement to reflect the information you provide in your cover letter.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: N/A

Reviewer #2: N/A

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: In this work the authors presented an iterative phase contrast CT reconstruction with novel tomographic operator and data-driven prior. The combination of conventional methods and deep learning methods is an important topic. The presented work combines L-BFGS optimization scheme with a deep learning parameterized prior for iterative reconstruction for breast cancer imaging. The theoretical analysis provides the motivation of proposing PnP-LBFGS. The experimental results on simulated and real datasets confirm the advantages of the proposed method.

This reviewer thinks that this paper could be improved in the following ways.

1. The compared baseline methods only include conventional methods. The literature has shown that the deep learning has shown better/competitive performance than iterative algorithms in some applications. Therefore, deep learning-based methods, such as post-processing FBP reconstructed images, should be considered as a competitive baseline.

2. I am wondering if the authors could provide the comparison of reconstruction times, which is a factor for clinical use.

3. In literature there are many works synergizing deep learning and iterative methods. A brief discussion on key references may be important to further highlight the novelty of the proposed methods.

Reviewer #2: The manuscript “Iterative phase contrast CT reconstruction with novel tomographic operator and data-driven prior” discusses an iterative procedure for phase-contrast tomography that make use of a neural network denoising step.

The manuscript is generally well-written and clear, even if many aspects of the proposed strategy are not properly discussed. In addition, often, the authors make statements that are not properly justified or that would need at least some references. About these remarks from my side, more detailed information can be found in the attached pdf.

In general:

1) I am not convinced the proposed approach can be called data-driven as the denoising part is associated with the model space (the images), not the data space;

2) In both the synthetic and the experimental tests, the proposed approach tends to remove many details (that are spatially consistent, so, probably not related to random noise, and that are connecting different subdomains in the reconstruction). My impression is that similar results could be obtained, for example, with a much simpler spatial filter or a more appropriate choice of the regularization within the framework of more "standard" approaches. This, of course, does not mean that the proposed approach is not interesting, but a more accurate and fair comparison with the “mainstream” approaches would make the paper more relevant to the readers.

3) A more in-depth discussion of the noise in the data would be important. Moreover, why no data covariance matrix is included in the inversion algorithm? It seems to me that the proposed strategy might be affected, for example, by severe modeling error that is not taken into account and that might lead to dangerous over-interpretation of the results. I do not expect the authors to modify their algorithm accordingly (at least for this manuscript), but I feel a short discussion about that would be important.

I hope these comments might be helpful in further improving the quality of the paper.

Best

P.S.

Authors must improve the quality of the images.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

**********

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

22 Jul 2022

Response to Reviewers:

Please find below our rebuttals to the editor’s and reviewers’ comments. We hope we could satisfactorily address the concerns. We highlighted the changes in red in the text.

Editor’s comments:

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

- We changed the figure labelling from Fig. to Fig and multiple figure labelling to Figs … and …

- We justified the text

- We changed the author affiliations

- We added the corresponding author’s initials after the corresponding email address

- We changed the last section’s name from conclusions to discussion

Please amend your Methods section to include the ethics approval number included in your Ethics Statement.

- We added the ethics statement

In your Data Availability statement, you have not specified where the minimal data set underlying the results described in your manuscript can be found. PLOS defines a study's minimal data set as the underlying data used to reach the conclusions drawn in the manuscript and any additional data required to replicate the reported study findings in their entirety. All PLOS journals require that the minimal data set be made fully available. For more information about our data policy, please see http://journals.plos.org/plosone/s/data-availability. Upon re-submitting your revised manuscript, please upload your study’s minimal underlying data set as either Supporting Information files or to a stable, public repository and include the relevant URLs, DOIs, or accession numbers within your revised cover letter. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories. Any potentially identifying patient information must be fully anonymized.

Important: If there are ethical or legal restrictions to sharing your data publicly, please explain these restrictions in detail. Please see our guidelines for more information on what we consider unacceptable restrictions to publicly sharing data: http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions. Note that it is not acceptable for the authors to be the sole named individuals responsible for ensuring data access.

We will update your Data Availability statement to reflect the information you provide in your cover letter.

- Due to legal restrictions, we are not allowed to share the patient data obtained with written informed consent under the ethical approval KEK-2012 554 granted by the Cantonal Ethics Commission of canton Zürich. We can however share our in-silico breast phantoms and the corresponding sinogram data. They will be uploaded to the ETH research collection archive as supporting material to our paper: https://www.research-collection.ethz.ch. The exact URL will follow.

Reviewers’ comments:

The compared baseline methods only include conventional methods. The literature has shown that the deep learning has shown better/competitive performance than iterative algorithms in some applications. Therefore, deep learning-based methods, such as post-processing FBP reconstructed images, should be considered as a competitive baseline.

- We added the comparison to a deep learning-based post-processing step in the “Results” section and in the “Effect of the proposed data-driven regularizer” section. For the comparison to be fair, we used the same network as within the iterative reconstruction. We would like to stress that this comparison only makes sense when it is possible to analytically reconstruct the phase contrast volumes (which is not always the case, e.g. in Teuffenbach et al., 2017).

- We accordingly edited Figs 9 and 10.

I am wondering if the authors could provide the comparison of reconstruction times, which is a factor for clinical use.

- We added a paragraph in the “Results” section describing the overall reconstruction times of the two iterative methods, as well as a more detailed description of the computational times of the single steps of the algorithms.

In literature there are many works synergizing deep learning and iterative methods. A brief discussion on key references may be important to further highlight the novelty of the proposed methods.

- We added a paragraph in the “Contributions” section that discusses prior work on combining deep learning with iterative methods and explained where the novelty of our method lies and why it is relevant to the scientific community.

I am not convinced the proposed approach can be called data-driven as the denoising part is associated with the model space (the images), not the data space;

- In the machine learning community, the term “data-driven” refers to algorithms that are fitted on training data, to differentiate them from algorithms that are solely designed by human engineering, and which do not contain trainable parameters. The fact that the mapping is applied in image space does not preclude the term data-driven in our opinion.

In both the synthetic and the experimental tests, the proposed approach tends to remove many details (that are spatially consistent, so, probably not related to random noise, and that are connecting different subdomains in the reconstruction). My impression is that similar results could be obtained, for example, with a much simpler spatial filter or a more appropriate choice of the regularization within the framework of more "standard" approaches. This, of course, does not mean that the proposed approach is not interesting, but a more accurate and fair comparison with the “mainstream” approaches would make the paper more relevant to the readers.

- The article already contains a comparison to the most celebrated and successful classical regularization scheme for CT, i.e. TV regularization. As it can be seen from the results in Figs 9 and 10 and in Table 1, TV regularization is not able to achieve comparable performance to the data-driven PnP strategy.

- Moreover, the new comparison with FBP-denoising uses a CNN that is composed of spatial filters, which can thus be regarded as a spatial filter.

- It is true that some small details are lost during the denoising. However, given the high amount of noise in the data and the very small size of the structures, we believe it is unrealistic to hope for those small features to be recovered.

A more in-depth discussion of the noise in the data would be important. Moreover, why no data covariance matrix is included in the inversion algorithm? It seems to me that the proposed strategy might be affected, for example, by severe modeling error that is not taken into account and that might lead to dangerous over-interpretation of the results. I do not expect the authors to modify their algorithm accordingly (at least for this manuscript), but I feel a short discussion about that would be important.

- We added a paragraph in the “DPC forward and backward tomographic operators” section in which we argue why we did assume constant variance in the data. The reason is that to accurately estimate the DPC variance, one needs to know the dark-field signal. Since the dark-field signal is hard to accurately compute based on highly noisy data, we decided not to include this into our model. We are currently working on a new algorithm which explicitly models the variance in the intensity data instead of the retrieved DPC data.

P.S. Authors must improve the quality of the images.

- All images were very high-resolution and all passed the PACE system test. Please let us know in case we must edit them.

Submitted filename: Response_to_Reviewers.pdf

Click here for additional data file.

1 Aug 2022

Iterative phase contrast CT reconstruction with novel tomographic operator and data-driven prior

PONE-D-22-15849R1

Dear Dr. Gogh

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Jude Hemanth

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: (No Response)

Reviewer #2: Partly

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: (No Response)

Reviewer #2: N/A

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: (No Response)

Reviewer #2: No

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: (No Response)

Reviewer #2: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

Reviewer #2: (No Response)

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

**********

16 Aug 2022

PONE-D-22-15849R1

Iterative phase contrast CT reconstruction with novel tomographic operator and data-driven prior

Dear Dr. van Gogh:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Jude Hemanth

Academic Editor

PLOS ONE