Learning-Based Approaches for Reconstructions With Inexact Operators in nanoCT Applications

Imaging problems such as the one in nanoCT require the solution of an inverse problem, where it is often taken for granted that the forward operator, i.e., the underlying physical model, is properly known. In the present work we address the problem where the forward model is inexact due to stochastic or deterministic deviations during the measurement process. We particularly investigate the performance of non-learned iterative reconstruction methods dealing with inexactness and learned reconstruction schemes, which are based on U-Nets and conditional invertible neural networks. The latter also provide the opportunity for uncertainty quantification. A synthetic large data set in line with a typical nanoCT setting is provided and extensive numerical experiments are conducted evaluating the proposed methods.


I. INTRODUCTION
T OMOGRAPHIC X-ray imaging on small scales such as on the nano scale -in short: nanoCT -is an important imaging technique that allows a visualization of the inner structure of small objects in the micro-or nanometer range with a suitable resolution.In contrast to computerized X-ray tomography on larger scales, for example in medical imaging, small vibrations that are transferred from the environment as well as manufacturing tolerances make a significant impact on the data: the resulting relative motion between object and scanner lead to motion artifacts when classical reconstruction techniques such as the filtered backprojection are applied.In this sense, nanoCT is an intrinsically dynamic inverse problem, i.e., an inverse problem that is heavily underdetermined since only one projection per state of the object is available.Hence, the unknown dynamic, T. Lütjen, F. Schönfeld, J. Leuschner, M. Schmidt, and T. Kluth are with the Center for Industrial Mathematics, University of Bremen, Germany (e-mail: {jleuschn,maximilian.schmidt,tkluth}@uni-bremen.de).
A. Wald is with the Institute for Numerical and Applied Mathematics, University of Göttingen (e-mail: a.wald@math.uni-goettingen.de).resp.inaccuracy in the forward map, needs to be taken into account.
Solving inverse problems typically relies on a proper formulation of the forward operator.But often and as mentioned before the forward map/model is not known with sufficient accuracy.Typical examples can be found in several imaging applications such as magnetic particle imaging [1] and also in CT applications [2], [3] as addressed in the present manuscript.Different directions taking into account the model's inexactness explicitly and/or implicitly have been considered in the past.An explicit treatment is, for example, the total least squares approach [4] explicitly taking into account a deviation on the operator when solving the original linear problem.This has also been extended and investigated with respect to regularization properties for general bilinear operators [5], [6] and also for a combination of model-and databased prior information [6].The original total least squares approach can also be equivalently formulated by a problem taking the inexactness implicitly into account.While this approach relies on the minimization of a certain type of functional, there also exists a class of iterative regularization methods taking into account the inexactness implicitly.In [7], an iterative reconstruction technique is proposed that relies on sequentially projecting iterates onto suitably defined subsets of the solution space, which are designed in such a way that the solution set is included in each of these sets and measurement noise as well as modeling errors are reflected in the design of the subsets.In the case of moving objects in CT, the unknown motion is interpreted as an inexactness with respect to a simpler model such as the (standard) Radon transform.This approach has been shown to yield promising results, at least for simulated data.More recently, a learning-based approach has been proposed to explicitly correct an inexact operator for the purpose of parameter reconstruction [8], which still relies on the knowledge of the exact forward operator.In general, several supervised learning-based approaches such as, for example, unrolled iteration schemes [9]- [11] can be interpreted as methods taking into account operator inexactness implicitly.Typically, the modeled and potentially inexact operator information is somehow incorporated into the structure of the trainable reconstruction method and the training is then performed with data pairs obtained from the true operator.
In the present work we address the problem of image reconstruction with an inexact operator in the framework of nanoCT applications.On the algorithmic side we consider sequential algorithms which take the inexactness of the operator into account and which still suffer from artifacts due to the inexactness.These methods are then extended by a learned post-processing scheme to obtain improved reconstructions less prone to the operator inexactness.In this context we use classical post-processing schemes via a UNet and we further develop new approaches based on conditional invertible neural networks, which allow for additional feature extraction for the purpose of uncertainty quantification.Furthermore, in order to evaluate the proposed methods a large synthetic dataset adapted to the nanoCT problem is generated.
The work is structured as follows: In Sec.II the general setting and the problem is specified.Sec.III then provides a detailed description of the non-learned sequential algorithms taking into account the operator inexactness, the learned post-processing schemes via neural networks with a particular focus on conditional invertible neural networks, and the data set generation.The methodological part is followed by extensive numerical experiments in Sec.IV and we conclude with a discussion in Sec.V.

II. PROBLEM DESCRIPTION AND SETTING
The small scale of the object in the low micrometer or upper nanometer range means that even small disturbances during the measuring process have a relatively large impact on the scanning geometry and can lead to inconsistent data.This means that the position of the object undergoes changes during data acquisition due to the resulting unavoidable relative motion between object and tomograph.For example, the X-radiation used in the scanner is generated by the interaction of an electron beam with a tungsten needle.Both the electron beam and the position of the needle may vary slightly.However, the largest impact on the relative motion is caused by the environment, when vibrations are transmitted to the tomograph, or by manufacturing tolerances in the tomograph itself, such as a thermal drift or an imprecise motion of the rotating stage on which the object is placed [12], [13].Additionally, the rotation axis may be tilted and the detector is skewed.In order to prevent reconstructions of being impaired by strong artifacts due to the relative motion and geometry deviations, this inconsistency has to be compensated in the reconstruction process.
We use the following notation to outline our methods and results: The unknown density function of the tested object (i.e., the ground truth) is denoted by x ∈ X , the measurable data by y ∈ Y, where X and Y are suitable vector spaces, e.g., scalar functions X ⊂ {x : Ω ⊂ R 2 → R} or X = R n×n representing an image.The model describing the relation between density and data is denoted by A, and we refer to it as the forward operator.Hence, we consider the underlying inverse problem as an operator equation Noisy data y δ := y + ξ with an additive noise vector ξ is supposed to have a noise level δ, such that y δ − y Y ≤ δ.If we use an inexact forward operator A η , we assume that the respective error in the model is given by η i.e., A 0 = A. The main goal is now to reconstruct an approximation x of x from (II.1) for unknown A and y but for given inexact A η and potentially noisy y δ , i.e., from the perspective of the given A η we have access to the corrupted measurement only, where y η would be a measurement obtained from A η .Those desired general reconstruction methods for the inverse problem (II.1) are denoted by T : Y → X in the remainder.Here, classical methods are combined with learning-based components in terms of neural network schemes F θ : X → X with network parameters θ.

III. MATERIALS AND METHODS
Classical reconstruction methods for CT include analytical inversion formulas, such as the standard filtered back-projection (FBP), as well as iterative reconstruction algorithms [14].However, these methods do not take into account potential inexactness of the forward operator as occurring for nanoCT.One classical iterative scheme is given by the Kaczmarz algorithm [15], which in the context of CT is known as the algebraic reconstruction technique (ART) [16].The Kaczmarz algorithm forms the basis for both non learning-based algorithms described in the following Sec.III-A, which extend it in different ways in order to handle operator inexactness.Then we consider the learning-based extension in Sec.III-B.
A. Sequential algorithms explicitly taking into account inexactness 1) Dremel: An approach for explicit geometry correction during the Kaczmarz algorithm was described in [12], which we will refer to as the Dremel method.It interleaves the iterations of Kaczmarz with correction steps that adapt the previously assumed geometry by maximizing the cross-correlation between the measured data and the forward-projection of the intermediate reconstruction over the plane of possible detector shifts.We focus on the case in which shifts are identified for each angle independently.While the correction step in the Dremel method naturally determines shifts in the detector plane, it can also be used to estimate object shifts (or source shifts) as they can be approximately translated into each other as described in [12].We use the Dremel method to estimate the object shifts while reconstructing from the perturbed data as we consider data sets featuring random object shifts (Sec.III-C).The full algorithm for the Dremel method is specified in the Supplementary Material A (Algorithm 1).
2) RESESOP-Kaczmarz: The combination of regularizing sequential subspace optimization (RESESOP) as introduced in [17], and Kaczmarz' method, in short RESESOP-Kaczmarz, was first introduced in [7] for linear semi-discrete inverse problems.For a mathematical analysis of this method we refer to [7].Within the present manuscript we write RESESOP to refer to RESESOP-Kaczmarz.The general idea of the method is outlined in the following.The respective algorithm in the version used for this work can be found in the Supplementary Material B (Algorithm 2).
We need to preface a few necessary concepts.Let us first define a set of linear mappings A k,l : X → Y k,l , where k ∈ K := {0, . . ., K − 1} and l ∈ L := {0, . . ., L − 1} are the index sets of scanner angles and detector points of a CT image.A k,l with k ∈ K and l ∈ L describes an X-ray of one scanner angle k and one sensor l on the detector plane.By A * we denote the adjoint of an operator A. The semi-discrete inverse problem reads and the solution set can be defined as and a stripe is defined by The algorithm then includes multiple concepts: SESOP component Sequential subspace optimization (SESOP) is an iterative method which can apply several search directions in each iteration step.
The full iteration step for SESOP is given by with a chosen finite index set I n , w n,i ∈ Y, search directions A * w n,i , and tn,i , i ∈ I n , minimizing the function This procedure is equivalent to calculating the metric projection P H (x n ) =: x n+1 of the current iterate x n onto the intersection of hyperplanes H := i∈In H(u i , α i ) with u i = A * w n,i and α i = ⟨w n,i , y⟩.Kaczmarz component SESOP alone does not take into account the nature of a semi-discrete problem A q x = y q , q = 0, ..., Q−1, (e.g., here q = (l, k), Q = K •L) since it is not defined for different realizations A q of the forward mapping A. Kaczmarz' method, however, requires such a discretized setting, and it works by iteratively projecting onto the solution set of one realization A q of A. The iteration step can then be written as where [i] = i mod Q. Kaczmarz' method thus makes it possible to merge subproblems in RESESOP-Kaczmarz.
Regularization The advantage of combining the SESOP algorithm with Kaczmarz' method is that the iteration can be adapted easily to regularizing semidiscrete inverse problems with modeling errors and noise levels depending on the subproblem.For inverse problems with inexact forward operator A η and noisy data y δ , it is important to introduce a form of regularization.As has been described in [7], SESOP-Kaczmarz can be regularized by replacing the hyperplanes H by stripes with a width defined by η and δ.Such a stripe is defined by with ρ > 0 and for a direction w ∈ Y.In the semi-discretized setting the method then terminates if a discrete version of the discrepancy principle is satisfied, i.e., for all k, l with constants τ k,l > 1 and noise and inexactness levels δ k,l , η k,l given for any semi-discrete subproblem.The full algorithm for the RESESOP method is specified in the Supplementary Material B (Algorithm 2).

B. Post-processing via neural networks
Neural networks F θ are finding increasing applications in the field of CT [18], [19].These include both detection and reconstruction tasks [20]- [22].A common approach for the latter is to post-process CT reconstructions of established methods T , such as filtered backprojection, using neural networks [23]-[27], e.g., where F θ : X → X .The idea is to correct artifacts and noise in the image using learned techniques.The postprocessing networks are trained in a supervised way with (simulated) ground truth data from idealized scanning conditions, e.g., with a higher dose, more scan angles, or less motion of the scanned object.These scan settings result in superior image quality but may not be feasible in practical use due to time or security constraints.The post-processing approaches have the advantage of working mainly in the space of images X .Thus, convolutional neural networks, like the U-Net [28] architecture are suitable for this application.In addition, both the classic and the post-processed reconstruction are available to the user.This can increase the acceptance of data-driven approaches among users who already have years of experience evaluating images from a particular method.Another network type, so called invertible neural networks, also allows for image postprocessing and it additionally allows for an uncertainty quantification.Compared to the U-Net architecture, the concept of conditional invertible neural networks, which are considered in this work, is less standard such that we highlight the important differences in architecture design, training, and reconstruction in the following.
1) Conditional invertible neural networks: For our application in nanoCT, the uncertainties due to the inaccuracy of the forward model must be taken into account.To this end, it is helpful to switch to the statistical view of inverse problems [29], [30].Here, instead of creating a single reconstruction, we are interested in recovering the whole conditional distribution p x|y (x|y δ ), respectively p x|x (x|x) with x = T (y δ ) in the post-processing setting of the present work.The idea is to approximate this distribution via a generative method, called conditional invertible neural network (iNN) [31]- [33] which have been applied to CT reconstruction [34], [35] from exact operators so far.
Architectures: Conditional invertible network architectures are typically given as mappings G θ : X × X → Z with latent space Z.The conditioning is encoded in the second input variable and invertibility is then given with respect to the first variable, i.e., Gθ,x := G θ (•, x) : X → Z is an invertible mapping for any x ∈ X .In the remainder we consider two specific architectures: • CiNN: A typical choice is a multi-scale architecture, based on NICE [36] and RealNVP [37].We use additive coupling blocks, the learned invertible downsampling proposed in [38], and a ResNet conditioning network processing the conditioning input before feeding it into the coupling blocks.The conditioning ResNet does not need to be invertible.• CiUNet: We also consider a conditional variant of the invertible U-Net proposed in [38], which has been studied in [35].We again use additive coupling blocks, and the conditioning network is implemented as a (non-invertible) U-Net, which is connected to the coupling blocks of the invertible U-Net at each respective scale.
Further details and illustrations of the architectures are provided in the Supplementary Material C.
Training: The goal of the conditional iNN during training is to minimize the Kullback-Leibler (KL) divergence to the conditional distribution.Equivalently, for any x this can be done using the negative log-likelihood where q θ is parameterized via the conditional iNN G θ and an assumed distribution p z , which allows for easy sampling, for a variable z defined on Z.Those are linked by the relation x|x = G−1 θ,x (z)|x, i.e, where z|x ∼ p z for any x.Joint training for all x in terms of expectation minimization and applying change of variables thus results in the considered training loss for a given data set of tuples (x (i) , x(i) ) including ground truth and reconstruction from a predefined method T .In the remainder we distinguish two variants of training for the conditional iNNs.First, we consider ground truth images x as the input variable or alternatively we consider the residual ∆x = x − x as input variable for the invertible path of the network.The latter case is denoted by the suffix Res in the present work (i.e., CiNNRes, CiUNetRes) as for given latent variable z the reconstruction can be interpreted as a ResNet with respect to x, i.e., (III.12) Reconstruction and uncertainty quantification: After training the conditional iNNs we can derive reconstructions and standard deviations from samples drawn from p z , i.e., given M samples z (m) ∈ Z and a x = T (y δ ), we obtain the final reconstruction x reco via which defines our learned post-processing routine F θ .In addition this approach allows for an uncertainty quantification, e.g., in terms of pixel-wise standard deviation, i.e., . (III.14) The Res case (here, G−1 θ,x (z (m) ) = ∆x (k) ) works analogously with minor adaptation, i.e., . (III.16)

C. Data set generation
In order to train the methods described in Sec.III, it is necessary to simulate tomographic imaging data including inexactness in the operator.The generated dataset is composed of 32 095 samples split into 30 490 (95%) for training, 1 284 (4%) for validation and 321 (1%) for testing.Phantoms x (i) were generated for a field of view of 255 × 255 pixels corresponding to spatial positions r ∈ R 2 and are constructed from randomly generated rectangles and ellipses which overlay each other and have different levels of density ranging from 0 to 1, 0 being no density.One phantom has a main shape and up to 3 subshapes within (see Fig. 1), which is inspired by common experimental settings in nanoCT [12], [39], [40].
The objective in the considered problem setting is to simulate distortions in the measurements caused by small vibrations near the CT scanner during the scanning process.As a result, this is incorporated as movements of the entire phantom with respect to the scanner position/scanning angle ϕ, where an example of randomly  generated movements is illustrated in Fig. 2. The vibrations have been mimicked by overlapping dampened sinus waves (38 for parallel and 9 for fan beam) with different starting points during the scanning process ϵ(ϕ).The starting points and parameters for the sinus waves are randomly generated for each sample, but limited to a maximum r 1 -and r 2 -direction shift distance from the center.
In addition to the sinus noise ϵ(ϕ), which is a deterministic function with respect to ϕ randomly differing for each phantom x, there is a small random noise ξ r−dir := (ξ r1−dir , ξ r2−dir ) t and ξ ϕ−dir taken normal distribution with randomly drawn standard deviation for each angle and applied to r 1 -, r 2 -direction shifts and rotations, each drawn for any angle ϕ.Here, the standard deviations are drawn from N (0.127, 0.0254 2 ).In summary, this can be formalized as considering the randomly moved phantom x non−perturbed (r) for the particular scan angle ϕ as (r ∈ R 2 ) where R is the corresponding rotation matrix in R 2 .
The generated phantoms and perturbations were then used to calculate the sinograms y η and y δ for the nonperturbed and the perturbed phantoms (see Fig. 3) for a parallel-and a fan-beam geometry.The sinograms were created using Operator Discretization Library (ODL) [41] and ASTRA Toolbox [42] and the respective geometry parameters from Tab. I.Note that the perturbed phantom x perturbed (ϕ), which is represented as a sequence of images with respect to ϕ, is assigned to the exact but unknown operator A case and that the non-perturbed phantom x non−perturbed is assigned to the inexact but known operator A η case.Here we use A η = A † which is chosen as the operator without phantom perturbation in any case, i.e., measurement from exact unknown operator are y with y(ϕ, s) = (A † x perturbed (ϕ))(ϕ, s) (s being the detector position) and hypothetical measurement from inexact known operator y η = A † x non−perturbed .In the present work emphasis is given on the operator inexactness such that the additive noise ξ = 0 is chosen in any case, i.e., y δ = y and δ = 0.The data set includes the tuples (x (i) , y (i) ) where non−perturbed and y (i) = y δ,(i) .In addition y η,(i) is also available in the data set.The dataset is made available under https: //doi.org/10.5281/zenodo.8123498.

IV. NUMERICAL EXPERIMENTS AND RESULTS
We conducted numerical experiments on the simulated phantom data set described in Sec.III-C.For this we proceed as follows: First, we evaluate each non-trained reconstruction method T on the respective data sets {(x (i) , y (i) )} i for parallel and fan-beam geometry to obtain the reconstructions x(i) T from reconstruction method T using A η = A † with A † being either a standard parallel or fan beam setting.Here, for T we used FBP, Dremel (see Sec. III-A1), and RESESOP (see Sec. III-A2).For the RESESOP method an estimate of η needs to be obtained prior the reconstruction.Here, we assume that a conservative estimate of η k,l for any semi-discrete subproblem is given in terms of the maximum deviation over all detector positions for each angle, i.e, η k,l = max ℓ∈L (∥y k,ℓ − y η k,ℓ ∥) for any l ∈ L. In the numerical experiments we also include the cases of ±20% overand underestimation.The Dremel method includes some algorithmic options, which were selected by evaluation on validation samples (see Supplementary Materials A).
Second, we train a learned post-processing routine F θ using the data tuples (x (i) , x(i) T ).For the UNet architecture F θ , the loss is given by the Euclidean norm, i.e., minimize i ∥F θ (x (i) T ) − x (i) ∥ 2 with respect to network parameters θ.The invertible network architectures are trained according to the descriptions in Sec.III-B1 using a standard Gaussian density p z .Here, we use a CiNN and a CiUNet as outlined in Sec.III-B.All trainings are performed using the Adam optimizer, and the checkpoint with minimum validation loss is selected.
Third (for iNNs only), the iNNs serve as an uncertainty estimator in terms of an estimate for pixel-wise standard deviation.
Finally, the resulting reconstruction schemes are then evaluated quantitatively (PSNR, SSIM) and qualitatively on the test set of the respective data set for parallel-and fan-beam geometry.
Further algorithmic details are provided in the Supplementary Materials A, B, and C. The code for the learned post-processing schemes via UNet and conditional iNNs is made available under https://gitlab.informatik.uni-bremen.de/inn4ip/cond-inn4nanoct.

A. Architecture specifications
The UNet architecture used in the present work is an adapted standard U-Net architecture with 5 scales and 32 to 128 channels using strided convolution downsampling, bilinear upsampling, skip connections and a sigmoid output activation.A detailed illustration of the used UNet can be found in the Supplementary Material C (Fig. 8).
The CiNN architecture has 6 scales and uses learned invertible downsampling, splitting, additive conditional coupling blocks, and a ResNet conditioning network.In each downsampling the number of channels naturally is multiplied by 4, after which half of the channels are split off to form a part of the latent output z (except after the first downsampling).The other half of the channels is further processed by a conditional coupling block, followed by the next downsampling, until the coarsest scale with 128 channels is reached.A small part of this signal is finally processed by random permutation and a conditional coupling block, before concatenating all parts of the latent output z.See the Supplementary Material C (Fig. 9) for an illustration.
The CiUNet architecture is based on the iUNet architecture [38] using 5 scales, learned invertible down-and upsampling, additive conditional coupling blocks, and skip connections forwarding half of the channels at each scale from the encoder to the decoder.A non-invertible U-Net is used for the conditioning, connected to the conditional coupling blocks at the respective scales in reversed order, i.e., U-Net's decoder activations are connected to CiUNet's encoder coupling blocks.After the downsampling in the encoder to the coarsest resolution, at which 32 channels are used, the features are upsampled back to the original resolution by the decoder, while concatenating with the forwarded channels from the skip connections.The decoder output then forms the latent output z (via reshaping).A sequence of four additional layers, namely activation normalization, downsampling, conditional coupling and upsampling, is inserted at the beginning before the encoder.An illustration and further details can be found in the Supplementary Material C (Fig. 10).
Reconstructions and standard deviations are computed with M = 100 samples from the latent space.

B. Non-perturbed phantom case
As a sanity check and also as a reference, we also trained and evaluated selected methods on unperturbed measurements, i.e., those measurements generated by unperturbed phantoms and A † = A η which is illustrated in Tab.II.All reconstruction methods T without postprocessing deliver accurate reconstructions as expected.The post-processing via UNet and the residual iNNs (CiNNRes, CiUNetRes) result in the largest performance for the phantoms in Fig. 1.

C. Perturbed phantom case
Quantitative results: The quantitative results on the perturbed phantom test data set for parallel beam are sum-  Qualitative results: Qualitatively, we can also observe different characteristics in the remaining artifacts in the reconstruction.The outcome of the non-learned reconstruction methods T which also serve as the input to the learned post-processing schemes is illustrated in Fig. 4. All methods suffer from streaking artifacts emerging mainly on edges and corners in the phantom.
The intensity of distortions decreases from FBP over Dremel to RESESOP.For FBP and RESESOP we can also observe more severe distortions within the objects.The post-processing via the UNet is illustrated in Fig. 5. Qualitatively, the reconstructions in Fig. 5 are superior to those provided solely by the methods T .In all cases the streaking artifacts are removed, also in the objects interior for FBP and RESESOP.Qualitative differences which cause quantitative differences become apparent in the difference images in Fig. 5. FBP+UNet and Dremel+UNet have more severe deviations particularly at phantom edges.For RESESOP these distortions are less distinct but in some cases other artifacts appear in the reconstruction, e.g., in the reconstruction at the bottom right in Fig. 5.
Qualitative results cond.iNNs: The CiUNetRes showed the best performance among the iNNs such that we restrict the presented qualitative results to this particular network case.Further image reconstructions of other methods and also for the fan beam setting are provided in the Supplementary Materials E and F. The post-processing via the CiUNetRes is illustrated in Fig. 6.Here, we can observe that the CiUNetRes tends to provide smoother reconstructions.In the FBP case we observe severe smoothing which is most likely the reason for the quantitative performance drop in PSNR.For Dremel and RESESOP we observe much less smoothing and in the RESESOP case the artifact in the interior of the right phantom disappeared.The extent of smoothing becomes more apparent in the difference images in Fig. 6.Here, the rectangular shapes show more severe differences in the corners where these effects decrease from FBP to RESESOP again.In any T case the CiUNetRes is able to remove the streaking artifacts but the more severe these artifacts have been the more smoothing is imposed which results in oversmoothed image reconstructions with varying intensity.Phantoms with smooth contours are likely to be advantageous for the CiUNetRes due to the observed effect in the corners.In sum, the CiUNetRes is able improve the reconstructions from Dremel and RESESOP but the quantitative improvement is smaller when compared with the UNet approach which is likely to be caused by the imposed invertibility in the network and the sampling mean which can also contribute to the observed smoothing effect.
Uncertainty estimation cond.iNNs: But the invertibility and the sampling also opens the door for extracting additional features.For this we illustrate the pixel-wise standard deviation in Fig. 7.In the FBP reconstruction the outer contour is only extracted which gives limited information about the actual uncertainty, e.g., when compared with the differences to ground truth illustrated in Fig. 6.In case of Dremel and RESESOP we observe Reconstructions Differences Fig. 5. Image reconstructions xreco and differences to ground truth using the non-trained reconstruction methods T combined with UNet post-processing on perturbed parallel beam data for the phantoms in Fig. 1. a larger similarity to the actual differences to ground truth.The standard deviations indicate an uncertainty in the reconstructed location of the edges.But similar to the actual reconstruction the uncertainty in non-smooth contours like the edges is not properly predicted.But on the other side the standard deviation can also indicate differences in the amount of uncertainty, e.g., being Reconstructions Differences Fig. 6.Image reconstructions xreco differences to ground truth using the non-trained reconstruction methods T combined with CiUNetRes post-processing on perturbed parallel beam data for the phantoms in Fig. 1. apparent in the middle phantom for RESESOP where the left and right inclusions in the bigger rectangle are more severely mis-located than the middle one.This can also be observed in the standard deviation.On the circular phantom at the left the standard deviation also predicts less uncertainty in the horizontal edges compared to the vertical edges which is in line with the difference to ground truth.

V. DISCUSSION AND CONCLUSION
In the present work we investigated different classical and learning-based computational strategies how to deal with inexactness in the underlying forward operator motivated by nanoCT applications where the inexactness has a severe probabilistic characteristic due to environmental influences.On the methodological side sequential algorithms already taking into account that the forward operator is inexact were combined with learning-based components such as classic U-Nets and more novel network architectures such as conditional invertible neural networks.The latter class also allowed for additional feature extraction for the purpose of uncertainty quantification.Overall, the non-learned reconstructions can be further improved by a learned post-processing network in most cases.In the proposed computational pipeline the performance of the learned component relies on the quality of the preliminary reconstruction, i.e, the postprocessing method only has access to the information from the output of the previous process and not from the actual measurement.Consequently, an initial reconstruction with too little or too distorted information can only be improved to a limited extent.Future works include an extension taking into account a conditioning depending on both actual measurement and preliminary reconstruction.
One particular focus of the work was the investigation of conditional iNN approaches as learned postprocessing tools for reconstructions from inexact operators as they also provide the opportunity to extract features to quantify the uncertainty.Here, the kind of training has a larger impact on the result, i.e., training with respect to the residual is superior.This might be related to the complexity of the learned feature distributions.In the residual case the error pattern of the preliminary reconstruction needs to be learned solely, while in the full reconstruction case the error pattern as well as the image class needs to be encoded in the network somehow.
The overall reconstruction performance in particular for the CiUNetRes has been only slightly worse compared to the classic U-Net approach.We observed that this might be related to an inherent smoothing effect in the iNNs' reconstructions along contours.Here, one needs to take into account that the underlying architectures have a different structure in particular within the invertible path of the iNN.The particular choice of invertible network architecture is likely to be one of the origins of the observed smoothing effect.In the present work we use coupling layers which tend to result in universal approximators of diffeomorphisms [43], [44] which are compared to classic U-Nets being composed of solely continuous and not necessarily differentiable layers.Here, future works may relax this property in the invertible path of the post-processing network using invertible residual networks [45] which tend to approximate homeomorphisms only [46].In summary, the conditional invertible networks for post-processing can provide a powerful tool to further improve image reconstructions from inexact forward operators and to simultaneously provide insights into the remaining operator uncertainty.
Besides the previously mentioned possible extensions, future research directions include the investigation of estimated higher order statistics, the application to measured nanoCT data, and also the transfer to other imaging modalities suffering from operator inexactness, e.g., such as magnetic particle imaging [1].SUPPLEMENTARY MATERIAL E ADDITIONAL FIGURES -PERTURBED PARALLEL BEAM Fig. 15.Difference of ground truth phantoms (Fig. 1) and T reconstructions (Fig. 4).

Fig. 1 .
Fig. 1.Three sample phantoms with 255 × 255 pixels included in the generated test data set.

Fig. 2 .
Fig. 2. Randomly generated perturbations for the data set generation with respect to scanner positions ϕ.

Fig. 3 .
Fig. 3. Non-perturbed (top) and perturbed (bottom) sinograms of the phantoms in Fig. 1 resulting from the perturbations like illustrated in Fig. 2. Scanner positions ϕ are on the horizontal axis and detector positions are on the vertical axis.

Fig. 4 .
Fig. 4. Image reconstructions using the non-trained reconstruction methods T on perturbed parallel beam data for the phantoms in Fig. 1.

Fig. 9 .Fig. 10 .Fig. 13 .Fig. 14 .
Fig. 9. CiNN architecture.Part of the information is factored out at different scales, and concatenated to form the latent vector z.The other part is continued to be transformed, using additive coupling blocks, which receive conditioning inputs from a ResNet, and invertible down-sampling.

Fig. 20 .Fig. 21 .
Fig. 17.Image reconstructions xreco and differences to ground truth using the non-trained reconstruction methods T combined with CiNNRes post-processing on perturbed parallel beam data for the phantoms in Fig. 1.

TABLE I PARALLEL
AND FAN BEAM GEOMETRY.S DENOTES THE NUMBER OF RECONSTRUCTION PIXELS, K DEFINES THE NUMBER OF SCANNING ANGLES (TIME STEPS), L DEFINES THE NUMBER OF DETECTOR PIXELS, AND D SPECIFIES THE SOURCE RADIUS IN RECONSTRUCTION PIXEL SIZE UNIT.THE DETECTOR EXTENT IS SUCH THAT THE PROJECTION BEAMS COVER THE WHOLE IMAGE.

TABLE III QUANTITATIVE
EVALUATION ON THE PARALLEL BEAM TEST DATA FOR MEASUREMENTS FROM PERTURBED PHANTOMS.OVERALL BEST PERFORMANCE IS HIGHLIGHTED IN BOLD FOR ANY T , BEST PERFORMANCE UNDER THE INN POST-PROCESSING NETWORKS IS HIGHLIGHTED IN COLOR, AND OVERALL BEST PERFORMANCE IS UNDERLINED.

TABLE IV QUANTITATIVE
EVALUATION ON THE FAN BEAM TEST DATA FOR MEASUREMENTS FROM PERTURBED PHANTOMS.OVERALL BEST PERFORMANCE IS HIGHLIGHTED IN BOLD FOR ANY T , BEST PERFORMANCE UNDER THE INN POST-PROCESSING NETWORKS IS HIGHLIGHTED IN COLOR, AND OVERALL BEST PERFORMANCE IS UNDERLINED.