arXiv:2606.26991v1 [eess.IV] 25 Jun 2026
Enabling self-supervised learned primal dual with Noise2Inverse Antti Sällinen1 , Siiri Rautio3 , Santeri Kaupinmäki1 , and Andreas Hauptmann1,2 1
Research Unit of Mathematical Sciences, University of Oulu, Finland 2 Department of Computer Science, University College London, United Kingdom 3 Department of Mathematics and Information Science, Josai University, Japan June 26, 2026 Abstract X-ray computed tomography reconstruction is an ill-posed inverse problem, particularly in low-dose and sparse-angle settings where measurements are noisy and incomplete. While learned reconstruction methods such as the Learned Primal-Dual algorithm achieve strong performance, they typically rely on supervised training with access to ground-truth data, which is often unavailable in practice. In this work, we propose a self-supervised reconstruction method by extending the Noise2Inverse framework to the Learned Primal-Dual algorithm. The resulting approach, called Noise2Inverse Learned Primal-Dual (N2I-LPD), enables training of a learned iterative reconstruction operator without ground-truth images by exploiting the statistical independence of noise in distinct measurements with respect to angular rotation of the CT-scan. We compare the proposed method with classical reconstruction methods, as well as neural network–based approaches such as a U-Net trained within the same N2I framework. The results demonstrate that N2I-LPD achieves improved reconstruction quality, highlighting the potential of combining learned reconstruction operators with self-supervised training strategies for practical CT imaging scenarios where ground-truth data is unavailable.
1
Introduction
X-ray computed tomography (CT) is a widely used imaging modality in medical and industrial applications due to its ability to provide high-resolution crosssectional reconstructions of objects. Mathematically, CT reconstruction can be formulated as an inverse problem where the goal is to recover an unknown object 1
from its line integral measurements [17, 23]. The inverse problem is ill-posed, particularly when the number of projections is limited, and becomes increasingly challenging in low-dose settings due to high noise levels [7, 10]. In scenarios where projection data is densely sampled over a full angular range and acquired with sufficiently high radiation dose, analytical reconstruction methods such as filtered backprojection (FBP) [23] produce reconstructions of sufficient quality. However, if these conditions are not met, the reconstruction task becomes significantly more challenging. In low-dose settings, the measurement is corrupted by increased noise, which is directly propagated to the reconstruction by FBP, resulting in degraded image quality. In sparse-angle settings, where only a limited number of projections are available, the data becomes incomplete and the inverse problem becomes more ill-posed. As a result, reconstructions suffer from artifacts. To address these challenges, reconstruction methods based on variational regularization and iterative optimization have been widely studied. In classical approaches, the reconstruction is formulated as a regularized optimization problem, for instance, using Total Variation (TV) [26], and solved with iterative methods such as the Primal-Dual Hybrid Gradient (PDHG) algorithm [6, 29]. In the past decade, machine learning approaches have emerged as powerful tools for CT reconstruction. These data-driven methods leverage neural networks to improve reconstruction quality and are now considered state-of-the-art in many settings [4, 32, 33]. One of the simplest ways to incorporate neural networks into CT reconstructions is through post-processing [4, 33]. In this approach, some knowledgedriven reconstruction is first applied to the measured data to obtain an initial reconstruction, which is then processed by a neural network. The network aims to remove noise and artifacts from the initial reconstruction. These methods are relatively easy to train and implement, but they typically require large amounts of training data. Convolutional Neural Networks (CNNs), particularly U-Net–based architectures, have been shown to perform well in this setting [16, 18]. An alternative approach is to learn the reconstruction operator as an unrolled iterative scheme [2, 3]. In these methods, neural networks are embedded into iterative reconstruction schemes, allowing the entire reconstruction process to be learned from data. A prominent example is the Learned Primal-Dual (LPD) algorithm [3], which unrolls a primal-dual optimization method into a trainable architecture. Most data-driven reconstruction methods are trained in a supervised setting, where paired data consisting of inputs (measurements or preliminary reconstructions) and corresponding ground-truth images are required. In CT imaging, this would require both a low-dose and a high-dose, densely sampled scan, where the latter serves as data for a reference reconstruction. However, in clinical practice, acquiring such data is infeasible, as it would expose patients to unnecessary levels of ionizing radiation by multiple and high-dose scans. This lack of groundtruth data constitutes a major limitation for supervised learning approaches in medical imaging. This limitation of supervised learning approaches motivates the use of selfsupervised learning methods, which enable training without ground-truth images by exploiting the statistical structure of the measurement noise. In recent years, several such methods have been proposed for imaging applications in2
cluding Noise2Noise [21], Noise2Self [5], Noise2Void [20], and specifically for CT imaging as Noise2Inverse [15]. Even more recent approaches exist, including Equivariance2Inverse [28] and Noisier2Inverse [11]. These approaches typically rely on the assumption that noise is mean-zero and statistically independent across measurements. In this work, we extend the self-supervised Noise2Inverse (N2I) framework, which utilizes the FBP, to the Learned Primal-Dual reconstruction method. This results in a self-supervised unrolled learned reconstruction operator that does not require ground-truth images during training. We evaluate the proposed method on simulated data as well as experimental data to test generalizability. For comparison, we train U-Net models both in a supervised manner and in a self-supervised manner using the Noise2Inverse framework. Additionally, we investigate the effect of the rectified linear unit (ReLU) activation function on training and reconstruction quality. The paper is organized as follows. In section 2, we introduce the inverse problem of X-ray tomography and explain how CT scanning is modeled mathematically. Section 3 recalls the Noise2Inverse framework, the Learned primal-dual method, and introduces how to combine them to train the learned reconstruction operator. In section 4, the experimental setup is discussed, including introducing the data sets and explaining the practical neural network implementation details. The results are presented in section 5. The section 6 is for discussion, where we lay out the implications of our findings and potential future research. Finally we have the section 7 for conclusions, where we shortly state the final remarks of this work.
2
X-ray tomography
The task of CT reconstruction can be framed as an inverse problem [23, 17]: g = A f + ε,
(2.1)
where g denotes the measured X-ray data (sinogram), A : X → Y is the X-ray transform describing the measurement process and geometry, and ε represents measurement noise which is assumed to be zero-mean Gaussian. The inverse problem (2.1) is ill-posed and due to real-world complications such as incomplete sampling and measurement noise, regularization is required to obtain stable reconstructions [9]. Here, we will use the X-ray transform to model a fan-beam geometry. It is an integral transform closely related to the Radon transform and integrates the attenuation coefficient along lines through the object [23]. The measurement data in CT are obtained by collecting line integrals of the object from multiple projection angles. The X-ray transform is defined as follows. Let f ∈ X = L2 (Rn ), we define lines L ⊂ Rn with directions defined by angles ϕ ∈ S n−1 , incident at detector locations s ∈ Rn , such that the line coordinates x ∈ L(ϕ, s) are parametrized by t ∈ R as x(t) = s − tϕ. We then have that the X-ray transform A : X → Y is defined by the line integral Z Z ∞ A f (ϕ, s) = f (x) dx = f (s − tϕ) dt. (2.2) −∞
x∈L(ϕ,s)
3
In a CT scanner, X-ray beams are emitted with initial intensity I0 , they pass through the object, and are detected with reduced intensity I1 . According to the Beer–Lambert law [23], the logarithmic attenuation measured by the detector corresponds to the line integral of the attenuation coefficient: Z ∞ I0 ln = f (s − tϕ) dt. (2.3) I1 −∞ This relationship forms the basis for modeling the measurement process and analyzing measurement noise. In practice, CT scanners acquire only a finite number of such measurements from different projection angles. Consequently, in practice the forward operator used in reconstruction algorithms corresponds to a discretized version of the X-ray transform. In this work, we consider the fan-beam geometry, which is a two-dimensional scanning setup where a point source emits X-rays that propagate through the object and are detected by a detector array. During acquisition, the sourcedetector system rotates around the object, collecting projections at multiple angular positions and producing the sinogram y. More precisely, this setting is defined with two variables: angle increments and the amount of X-rays. The angle increments give one dimension of the sinogram. Those increments tell us how many angular samples the sinogram has. The second dimension is given by the amount of X-rays, and it informs the density of the detected information. To add, angle increments act as an independent variables, which is necessary to know for Noise2Inverse later in this paper. The goal of the inverse problem is to recover x from the measurements in (2.1). There are many approaches to solving this problem, ranging from analytical reconstruction methods to variational regularization and modern machine learning techniques. One classical approach is the Filtered Back Projection (FBP) algorithm, which approximates a filtered inverse of the X-ray transform in (2.2) [23]. The filtering operation is defined as Z ∞ 1 (Fg)(x)|r|eirt dr. Gg(ϕ, s) = 2π −∞ The FBP algorithm reconstructs the object by first filtering the projection data in the frequency domain as above and then backprojecting the filtered projections over the image domain. The reconstruction formula is then given by f ≈ A† g := A∗ Gg = A∗ G A f,
(2.4)
where A∗ is the adjoint of the forward operator A in certain geometry. The FBP method provides very fast reconstructions, but assumes densely sampled projection data covering the full angular range. When these assumptions are violated, for example in sparse-angle settings, the resulting reconstructions contain artifacts. Furthermore, in the presence of noisy measurements, such as in low-dose CT imaging, the reconstructed images are typically noisy since the FBP algorithm does not incorporate explicit denoising or regularization. An alternative classical approach is variational regularization, which incorporates prior information to stabilize the reconstruction. This is particularly 4
useful in ill-posed situations, such as when the data are noisy or incomplete. The inverse problem (2.1) can be formulated as the minimization problem arg min Eα (f ) = arg min || A f − g||22 + α|||∇f |||1 , f
(2.5)
f
where ||·||22 is the data discrepancy term and ||·||1 is the regularization functional with regularization parameter α controlling the trade-off between data fidelity and regularization. In this work, we consider Total-Variation (TV) regularization [26], as written above. Many modern reconstruction algorithms can be interpreted as iterative schemes for solving variational problems of the general form of (2.5). In recent years, learning-based approaches have been proposed that replace parts of these iterative algorithms with trainable neural network components. This leads to learned reconstruction schemes that combine the structure of classical optimization methods with the expressive power of deep learning. Such learned reconstruction operators can be designed so that they are able to solve the inverse problem 2.1. Let us consider a family of parameterized mappings Rθ : Y → X, called learned reconstruction operators [13], where θ ∈ Θ is a specific parameter configuration obtained from the parameter space Θ via a training procedure. The goal is to provide suitable training data, both in quality and amount, in order to train the operator Rθ , solving the inverse problem as f = Rθ (g). Post-processing methods, which pair FBP reconstruction with a neural network, first apply an analytical reconstruction operator A† to obtain an initial reconstruction, which is then processed by a neural network Λθ . The network takes the initial reconstruction as input and outputs a denoised or artifactreduced image. This kind of learned reconstruction operator can be expressed as Rθ = Λθ ◦ A† . Such post-processing methods are relatively simple to train and implement in practice. On the other hand, unrolled learned reconstruction operators aim to learn the reconstruction process itself. The goal is to learn a reconstruction operator that maps the measured data to high-quality image reconstructions. These methods incorporate the forward model into the network architecture and learn the reconstruction mapping directly from data. One prominent example is the Learned Primal-Dual (LPD) algorithm [3], which unrolls a primal-dual optimization scheme into a trainable neural network architecture.
3
Methods
We extend the self-supervised Noise2Inverse framework [15] to the Learned Primal-Dual reconstruction method [3]. This yields a learned reconstruction operator that can be trained without access to ground-truth images, in contrast to supervised approaches such as U-Net–based denoisers [12]. In the following, we first describe the Noise2Inverse and Learned Primal-Dual methods separately, and then present their combination within a unified reconstruction framework. 5
3.1
Noise2Inverse framework
Noise2Inverse [15] is a self-supervised framework, which utilizes independent subsets of projection angles in the training of the neural network, giving the possibility to train the networks without ground-truth images. Considering the full sinogram data, it can be interpreted as a collection of sparsely sampled sub-sinograms. For example, if one has every even angle increment in one sub-sinogram, and every odd in the other, then one can combine these subsets into a more densely sampled one. Under the assumption that the measurement noise at different projection angles is mean-zero and statistically independent, the resulting sub-sinograms also contain independent noise realizations. Analytical reconstruction methods can then be applied separately to each sub-sinogram, producing reconstructions with independent noise components. Noise2Inverse utilizes these reconstructions as input–reference pairs for neural network training. Since the noise realizations vary independently between the input and target images, the network cannot consistently learn the noise and instead learns the underlying image structures. This splitting strategy can be implemented in different splits, which will be discussed later. More precisely, let us consider a set of sinogram data {g̃i }N i=1 , where g̃i denotes a noisy measurement g̃i = A(fi ) + ε. Here, the measurement operator A has dimensions Nα × M , where Nα denotes the number of projection angles and M the number of detector elements (beams). Thus, each measurement g̃i is obtained from Nα projection angles α1 , α2 , . . . , αNα . The Noise2Inverse method proceeds as follows. First, a hyperparameter K is chosen, which defines an angle splitting αj , αj+K , αj+2K , . . .. The sinograms are then split into sub-sinograms g̃i,1 , . . . , g̃i,K , where each g̃i,j contains measurements corresponding to every Kth projection angle. After splitting the data, sub-reconstructions are obtained as f˜i,j = Rj (g̃i,j ),
j = 1, . . . , K.
Here Rj is a reconstruction operator, and in this work, we use the filtered backprojection operator. Next, the measurements are partitioned into a collection J of index sets J ⊂ {0, 1, . . . , Nα }. For each J ∈ J , we define the mean reconstruction 1 X 1 X˜ fi,j = R(g̃i,j ). f˜i,J = |J| |J| j∈J
j∈J
In this way, input–target pairs for training are constructed as (f˜i,J C , f˜i,J ), lead-
6
ing to the training problem θ̂ = arg min θ
N 1 XX ∥Λθ (f˜i,J C ) − f˜i,J ∥22 |J | i=1
(3.1)
J∈J
N 1 XX = arg min ∥Λθ (R(g̃i,J C )) − f˜i,J ∥22 |J | i=1 θ
(3.2)
N 1 XX = arg min ∥(Rθ (g̃i,J C ) − f˜i,J ∥22 |J | θ i=1
(3.3)
J∈J
J∈J
Note that here, in contrast to the original Noise2Inverse publication [15], we consider a learned reconstruction operator Rθ , since the aim is to utilize measurement data in the training. ∗ To obtain a clean reconstruction fi,out , the network is evaluated as ∗ fi,out =
1 X Λθ̂ (f˜i,J C ), |J |
(3.4)
J∈J
that is, predictions from different splits are averaged. The key idea of Noise2Inverse is that the input and target images contain independent realizations of noise. As a result, the network learns to reconstruct the underlying signal while suppressing noise components. Intuitively, since the noise varies across the splits, it cannot be consistently learned by the network. Let g and f denote random variables in the measurement space Y and the image space X, respectively. The noisy measurements are given by g̃ = A(f̃) + ε. The corresponding sub-reconstructions are f̃J C = RJ C (f̃J C ), and f̃J = RJ (f̃J ). The trained network Λθ̂ from equation 3.1 approximates the regression Λ†θ = argmin θ
1 X Ef,ε ||Λθ̂ (f̃J C ) − f̃J ||22 , |J |
(3.5)
J∈J
which minimizes the expected prediction error. To generalize this, let us randomize J, denoting it by J. Now J takes uniformly random indices from the collection J . This generalizes the equation 3.5 such that it becomes Λ†θ = Eµ ||Λθ̂ (f̃JC ) − f̃J ||22 , (3.6) where µ is a joint measure of f, ε, and J. Before stating a proposition, let us define a clean sub-reconstruction of a clean measurement as fJ∗ = RJ (gJ ). Proposition 3.1. Let f̃J , f̃JC , fJ∗ , and µ be defined as above. Let ε be mean-zero and element-wise independent. Also, let RJ be a linear operator for all J ∈ J . Then, for all learned reconstruction operators Rθ : Y → X, we have Eµ ||Rθ (f̃JC ) − f̃J ||22 = Eµ ||Rθ (f̃JC ) − fJ∗ ||22 + Eµ ||fJ∗ − f̃J ||22 . (3.7) 7
The proof for the same case with regression functions h can be found in [15]. The proposition shows that the expected prediction error decomposes into a supervised prediction error term and a noise variance term. In practice, the splits are implemented by selecting subsets of projection angles. For experimental sinogram data, this can be done by indexing the sinogram array along the angular dimension. For simulated data, the splits can be implemented more accurately by defining corresponding forward operators for each subset.
3.2
Effect of ReLU on Noise2Inverse training
In this section, we analyze how the ReLU activation function influences training when used to enforce non-negativity in the Noise2Inverse setting. We use ReLU in all training, testing, and validation scenes, where in training we found that the loss is higher when ReLU is used, compared to the case where ReLU is absent. The proposition 3.1 does not include this case, but only considers the testing. Next, we derive the reasoning behind that. Let us consider the L2 loss functional L(θ) =
1 X ||Λθ (f˜J C ) − f˜J ||22 . |J | J∈J
Let the output of the network be a reconstructed image with an additive noise component εs . We denote this as Λθ (f˜JC ) = f˜J′ C + εs , where εs is a random variable representing residual noise. Then, using the definition of the ReLU function ReLU(x + y) = max(x + y, 0) =
x + y + |x| + |y| x + y + |x + y| ≤ , 2 2
we obtain f˜′ C + εs + f˜J′ C + εs f˜J′ C + εs = J 2 f˜J′ C + εs + |f˜J′ C + εs | ≤ 2 = ReLU(f˜J′ C + εs ) f˜′ C + εs + |f˜J′ C | + |εs | . ≤ J 2
8
With the inequality above, we can derive that 1 X ||Λθ (f˜J C ) − f˜J ||22 |J | J∈J 1 X ≤ ||ReLU(Λθ (f˜J C )) − f˜J ||22 = LReLU (θ) |J |
L(θ) =
J∈J
1 X |J |
f˜J′ C + εs + |f˜J′ C + εs | − f˜J 2
1 X ≤ |J |
f˜J′ C |f˜′ C | |εs | εs + + J + − f˜J 2 2 2 2
2
f˜J′ C
|εs | 2
=
J∈J
J∈J
≤
1 X |J |
2
J∈J
+
εs + 2
|f˜J′ C | 2
2 2
2
− f˜J
+ 2
2 2
. 2
If we assume that the reconstruction is non-negative, then |f˜J′ C | = f˜J′ C , and the equation above reduces to L(θ) ≤
1 X ˜′ εs − f˜J fJ C + |J | 2 J∈J
2
+ 2
|εs | 2
2
.
(3.8)
2
From equation 3.8, we observe that an additional constant term |ε| 2 is introduced into the loss. Furthermore, in the comparison between the output and the target reconstruction, the noise is effectively reduced by a factor of 12 .
3.3
Learned Primal-Dual method
Learned Primal-Dual (LPD) algorithm [3] is a data-driven approach for solving image reconstruction problems. It belongs to the class of learned iterative reconstruction algorithms, which are inspired by classical iterative optimization methods [13]. The key idea is to use machine learning to learn the update rules at each iteration based on the current iterate. The LPD method is based on adapting convolutional neural networks (CNNs) within the Primal-Dual Hybrid Gradient (PDHG) algorithm [6]. In contrast to the classical PDHG algorithm, the proximal operators are replaced by learned proximals, that are neural networks. Algorithm 1 Learned Primal-Dual 1: Initialize f0 ∈ X, h0 ∈ U 2: for i = 1, . . . , I do 3: hi ← Γθd hi−1 , Afi−1 , g i
4: fi ← Λθip fi−1 , A∗ hi 5: end for (1) 6: return fI
The LPD method is presented in Algorithm 1. The algorithm consists of 2I learnable parameter spaces; Γθid , corresponding to the dual updates, and Λθip , corresponding to the primal updates. These components are implemented as 9
neural networks and are separately optimized during training. In this work, the operator A denotes the ray transform and A∗ its adjoint operator. The same operator A is also used as the forward model to generate the measurement data from the images in the first place. For the network architecture, we employ residual networks [14]. The effect of different parameterizations is presented in Section 5.
3.4
Learned Primal-Dual in the Noise2Inverse framework
To combine the LPD reconstruction method with the N2I framework, several modifications are required, as the LPD architecture must be adapted to operate on split measurement data. First, due to the data splitting in N2I, the forward operator A and its adjoint A∗ must be adapted to the corresponding subsets of projection angles. This leads to a collection of split operators Ak and A∗k , where k = 1, . . . , K denotes the split index. Second, since the LPD algorithm operates on full sinograms in the dual space, the split sinograms must be recombined into a single sinogram. This is achieved using the sinogram combining function in Algorithm 3, which arranges the split sinograms according to their angular ordering, as illustrated in Figure 1. The overall workflow of the proposed method is described in Algorithm 2. We initialize the primal and dual variables f0 , g ′ , and h0 , where f0 is given by an FBP reconstruction, g ′ is the full sinogram obtained from the splits, and h0 is initialized as a zero tensor of the same size as g ′ . Within each unrolled iteration, the split sinograms are first computed by applying the operators Ak to the current primal iterate, resulting in sub-sinograms stored in gf . These are then combined into a full sinogram g ∗ using the sinogram extension function. The dual variable hi is updated by passing hi−1 , g ∗ , and g ′ through the dual network. This gives the full sinogram as an output. After the dual network, an analytical reconstruction operator A†k , (for example, FBP) is applied to obtain K sub-reconstructions fh . These are averaged to form a single reconstruction, which is then used to update the primal variable fi via the primal network, using fi−1 and fh . After I unrolled iterations, the network produces an approximation of the reconstruction operator, and the final output is used for training and evaluation. This defines the learned reconstruction operator Rθ for the N2I-LPD method as outlined in Algorithm 2. The code for the implementation of this algorithm can be found from GitHub [27]
4
Experiments
In this section, we describe the CT data used in our experiments, how the corresponding measurements are simulated, and details of the neural network implementation.
10
Algorithm 2 Noise2Inverse Learned Primal-Dual 1: Initialize f0 ∈ X j×x1 ×x2 , g ′ ∈ Y y1 ×y2 , h0 ∈ Y y1 ×y2 2: for i = 1, . . . , I do 3: {gfk }jk=1 ← Ak (fi−1 ) for k = 1, . . . , j 4: 5:
g ∗ ← combine_sinogram({gfk }jk=1 ) hi ← hi−1 + Γθid (hi−1 , g ∗ , g ′ )
{hki }jk=1 ← split(hi ) 7: {fhk }jk=1 ← A†k (hki ) for k = 1, . . . , j Pj 8: fh ← 1j i=1 fhi 9: fi ← fi−1 + Λθip (fi−1 , fh ) 10: end for 11: Rθ (g) ← ReLU(fI ) 6:
Algorithm 3 combine_sinogram 1: Initialize g1 ∈ Y y1 ×y2 , g2 ∈ Y k×y1 ×y2 2: for k = 0, . . . , n do
for j = 0, . . . , y1 do k+|k|·j,y2 4: g1 ← g2k,j,y2 5: end for 6: end for 7: return g1 3:
4.1
Walnut data set
The data set consists of reconstructions of 42 different walnuts, acquired using a high-dose cone-beam CT setup [8]. For each walnut, the data set provides a 5013 tensor of reconstructed cross-sectional images. These reconstructions are effectively noise-free and are therefore suitable for supervised learning, where they can serve as ground-truth images. For more information about the data set, we refer the reader to [8]. Since we consider a self-supervised setting in this work, we simulate the measurement data from these noise-free reconstructions. To do this, we use ODL [1], a Python package. Specifically, we employ the FanBeamGeometry and RayTransform operators to generate sinograms from the reconstruction data. The fan-beam geometry used for each reconstruction is parameterized as follows: 512 projection angles uniformly distributed over [0, 2π], and 496 detector elements (beams), with a source-to-object distance of 2.0 and an object-to-detector distance of 1.0. After generating the sinograms, we add 5% zero-mean Gaussian noise to simulate a realistic measurement setting. Reconstructions are then obtained using the ODL function fbp_op. Finally, we note that some slices in the data set contain mostly air. These slices are excluded from our experiments, as they are not relevant for the reconstruction task considered in this work.
11
Figure 1: Visual explanation of the sinogram combining function (Algorithm 3). This example is for four splits.
Model U-Net max U-Net min
# Parameters 7696706 25522
First channel 64 16
Depth 4 2
Table 1: Architectures of the U-Net models considered in this work. We include U-Net max with 7,696,706 parameters and U-Net min with 25,522 parameters. Here, First channel denotes the number of channels in the first layer of the U-Net architecture, and Depth refers to the total number of layers.
Model LPD max LPD min
# Parameters 578260 28730
# Channels 32 16
# Layers 4 2
# Unrolled iterations 10 5
Table 2: Architectures of the Learned Primal-Dual models considered in this work. We include LPD max with 578,260 parameters and LPD min with 28,730 parameters.
12
4.2
Network architecture and training
We use the U-Net [25] and ResNet [14] architectures in our experiments, following the setups in Noise2Inverse [15] and Learned Primal-Dual [3]. For the ResNet, we used the PReLU activation function between convolutional layers, following the original LPD implementation. In addition, we applied a final ReLU activation to both networks to enforce non-negativity of the output. For the U-Nets, we used skip connections and 3 × 3 kernels for all 2D convolutions. For max-pooling and transposed convolutions, we used 2 × 2 kernels. Table 1 shows how the number of channels in the first convolutional layer and the network depth affect the total number of parameters. For the ResNets, we used 3 × 3 kernels for the 2D convolutions. The number of parameters was varied by changing the number of layers in the network, as well as the number of unrolled iterations in the algorithm. The corresponding configurations are summarized in Table 2. All neural networks were implemented using PyTorch [24]. Training was performed for 50,000 iterations using the ADAM optimizer [19] with a learning rate of 0.001, together with a cosine annealing learning rate schedule. We also applied gradient clipping with a maximum gradient norm of 1, with respect to the ℓ2 -norm. For the loss function, we used the standard L2 -loss in the training and evaluation, implemented using PyTorch’s MSELoss.
4.3
Experimental setup
In the Noise2Inverse framework, the data are partitioned into different bins (disjoint subsets). Following [15], we use a 3:1 split, where three subsets are used as input and one subset as comparison data. In our experiments, we also observed that a 1:1 split is not suitable for the N2I-LPD method, although the exact reason for this remains unclear. This splitting strategy effectively leads to sparse-angle reconstruction, since increasing the number of splits reduces the angular sampling in each sub-reconstruction. For the network inputs, we use sinograms in methods that operate in the measurement domain, and corresponding FBP reconstructions for methods that operate in the image domain. At test time, the final reconstruction is obtained by evaluating the trained network on all four splits and averaging the resulting outputs.
4.4
Comparison methods
We consider several other reconstruction methods for comparison. The walnut data set contains ground truth images that were computed using FDK and accelerated gradient decent algorithm. For more information, see the original publication of the data set [8]. FBP reconstructions are computed from the full noisy sinogram, and a classical Primal-Dual reconstruction [6] with TV regularization. For the FBP reconstruction, we use the ODL function fbp_op, which provides an approximate inverse of the RayTransform. The parameters used are padding=1, filter_type=’Ram-Lak’, and frequency_scaling=1.0. The Primal-dual TV reconstruction algorithm is implemented as outlined in [29], using isotropic total variation. The forward operator is constructed us-
13
(a) Ground truth
(b) Low-dose FBP
(c) PD TV-I
(d) PD TV-U
(e) N2I-Unet min
(f) N2I-Unet max
(g) N2I-LPD min
(h) N2I-LPD max
Figure 2: Comparison of results for the Walnut data set. ing the ASTRA Toolbox [30, 31] along with the Spot linear operator toolbox [22]. Standard values for the primal-dual algorithm parameters were chosen: θ = 1, σ = τ = 1/L, where L = ||(A, ∇)T ||2 , and the forward operator A follows the fan-beam geometry outlined in 4.1. In order to translate the ODL geometry parameters into a format suitable for ASTRA, a specific geometry dependent scaling factor1 must be applied to the detector pixel size, the sourceto-object distance, and the object-to-detector distance. A scaling factor value of 1.57356942 was empirically determined to produce the best image registration between the ASTRA and ODL reconstructions. The finite differences in the total variation calculations were scaled by an empirically determined step size of 1/h = 85/(scaling factor) ≈ 54, which was found to best balance the magnitude of the data discrepancy and total variation terms in the primal-dual iterations.
5
Results
In this section, we present the results of the proposed Noise2Inverse Learned Primal-Dual method. Both quantitative and qualitative evaluations are considered. For comparison, we include filtered backprojection, classical Primal-Dual methods with TV regularization, and N2I U-Net reconstructions. See Figure 2. For the classical primal-dual method with TV regularization, reconstructions were generated for all of the 50 sinograms using 31 different regularization parameter values chosen uniformly from the interval [0.0025, 0.019]. The PSNR with respect to the respective ground truth image was computed for each of these 1 url: https://odlgroup.github.io/odl/_modules/odl/tomo/backends/astra_cuda.html 2 The
theoretical scaling factor is given by (source-to-object distance + object-to-detector distance)/source-to-object distance, which in our geometry equates to 1.5. However, this value leaves a few pixels of registration error between the ASTRA and ODL reconstructions.
14
35
30 28
Mean PSNR
PSNR
30
25
20
26 24 22
15
20 0.005
0.01
0.015
0.005
Regularization parameter
0.01
0.015
Regularization parameter
(a) PD TV-I: individualized regularization (b) PD TV-U: uniform regularization paparameter selection. rameter selection.
Figure 3: Regularization parameter choice for classical primal-dual with isotropic TV regularization. PSNR is computed for each choice of regularization parameter with respect to the ground truth image, for all 50 reconstructions. (a) PSNR versus regularization parameter curves for all reconstructions. The maximum of each curve is indicated by the solid circles. (b) Mean PSNR value for a uniform choice of regularization parameter for all reconstructions. Dashed vertical line indicates the regularization parameter producing the maximal PSNR. reconstructions. Two different methods for choosing the regularization parameter were considered: one where we choose the regularization parameters which maximize the PSNR for each individual reconstruction, and another where we choose a uniform regularization parameter which maximizes the mean PSNR of all the reconstructions, as seen in Figure 3. The individualized regularization parameter choice in Figure 3a produces a best-case scenario for the average PSNR, whereas the uniform regularization parameter choice in Figure 3b presents a more realistic value for the average PSNR of classical primal-dual total variation reconstructions, as ground truth images are not available in practice. Table 3 summarizes the quantitative results for all methods. The results show that the N2I-LPD with the largest parameter configuration achieves the best performance, with an average PSNR of 29.62 dB. The N2I-LPD with the smallest parameter configuration and the TV-regularized PD method both achieve an average PSNR of approximately 29.00 dB. Both N2I U-Net models reach an average PSNR of 28.67 dB. In comparison, FBP performs significantly worse, with an average PSNR of 19.01 dB and SSIM of 0.20. The evaluation was performed on 50 slices that were not used during training. For SSIM, the Primal-Dual TV method achieves the highest value (0.88), while the learned methods yield values between 0.68 and 0.76. All quantitative metrics were computed using the PSNR and SSIM implementations from PyTorch Ignite. We verified the PSNR and SSIM values for the TV-regularized reconstructions using both PyTorch Ignite and MATLAB, obtaining consistent results. Next, we discuss the qualitative results presented in Figure 2. Although the TV-regularized reconstructions achieve higher SSIM values, the visual differ-
15
ences in fine structural details compared to the N2I U-Net reconstructions are relatively small. The N2I U-Net results tend to preserve slightly sharper fine details, while the classical PD method produces stronger denoising effects. For the N2I-LPD reconstructions, both configurations exhibit improved preservation of fine details and smoother noise suppression compared to the other methods. In particular, they avoid the piecewise-constant appearance typically associated with TV regularization. Overall, the N2I-LPD methods achieve the best performance in terms of PSNR and provide visually improved reconstructions compared to the other methods.
5.1
Generalizability
We also investigate the generalizability of the trained networks. To this end, we consider the cone-beam CT (CBCT) projections provided in the walnut data set, and construct a corresponding 2D-sinogram from them. This is done by using the projections acquired when the X-ray source is in the central position, as described in [8]. For each projection, the center column is taken and assigned to the sinogram one column at a time. After constructing the full sinogram, the corresponding ODL geometry and operators are defined, after which the trained networks are evaluated on this previously unseen dataset. Figure 4 presents the reconstructions obtained by the aforementioned networks on previously unseen data. We note that the reference reconstructions provided in the Walnut dataset are computed using the FDK algorithm [8], and therefore do not exactly correspond to the FBP-based reconstruction setup considered in this work. Consequently, obtaining a perfectly matching ground-truth reconstruction for comparison is not possible in this setting. However, since the considered data are noise-free, the FBP reconstruction can still be used as a reliable qualitative reference. The ground truth reference reconstruction from the dataset is still included as the closest available approximation to the ground truth. While most structures are consistent across the reconstructions, Figure 4 highlights regions where visible differences occur, indicated by the blue square. All networks produce qualitatively similar reconstructions. However, in the highlighted region shown by the red square in Figure 4, N2I-LPD max provides the visually most accurate reconstruction. We also provide the difference map (see Figure 5) between the FBP reconstruction and the N2I-LPD max reconstruction, illustrating that no significant structural details are lost between the two reconstructions. Table 4 shows that Method FBP PD TV-I PD TV-U N2I U-Net min N2I U-Net max N2I-LPD min N2I-LPD max
PSNR 19.01 29.00 28.31 28.67 28.67 29.39 29.62
SSIM 0.20 0.88 0.89 0.68 0.70 0.75 0.76
Table 3: Results for the walnut data set for all different methods used. Best values for PSNR and SSIM are highlighted in bold.
16
(a) Ground truth
(b) Low-dose FBP
(c) N2I-Unet min
(e) N2I-Unet max
(f) N2I-LPD min
(g) N2I-LPD max
Figure 4: Comparison of results for the generalized setting. all neural networks evaluated in this generalized setting produce quantitatively similar results in terms of PSNR. These observations suggest that the networks generalize without introducing hallucinated structures into the reconstructions. Method N2I-LPD max N2I-LPD min N2I U-Net max
N2I U-Net min 36.73 39.27 41.54
N2I U-Net max 36.81 39.38
N2I-LPD min 38.26
Table 4: PSNR-values for all the combinations of the different network results in generalized setting used in this study.
6
Discussion
The results indicate that the proposed Noise2Inverse Learned Primal-Dual method achieves superior performance in terms of PSNR and visual reconstruction quality, while the TV-regularized Primal-Dual method attains the highest SSIM values. These findings suggest that incorporating learned reconstruction operators within self-supervised frameworks is a promising direction for low-dose and sparse-angle CT. These results further motivate several open questions regarding the proposed approach and its implementation. First, in Algorithm 2, the dual network operates on a full sinogram constructed from the split measurements. An alternative design would be to assign 17
0.100
0.075
0.050
0.025
0.000
0.025
0.050
0.075
0.100
Figure 5: The difference map between the FBP reconstruction and the N2I-LPD max reconstruction in the generalized case. separate dual-space networks to each split. Such a modification could potentially allow the model to better exploit the statistical independence and zero-mean assumptions of the noise across the splits. A similar extension could be considered in the primal space, where instead of averaging back-projections prior to the network, separate networks could be applied to each reconstruction and their outputs combined. Second, the choice of split ratio in the Noise2Inverse framework remains an open question. In this work, we follow [15] and use a 1 : 4 split, which was shown empirically to perform well. They divided the splitting strategy to be either 1 : X or X : 1. However, more general splitting strategies, such as Y : X or X : Y configurations, may affect both training stability and reconstruction quality. A systematic study of different split ratios could provide further insight into the behavior of the method. Finally, an interesting question remains of how sparse of a setting could be used in the training to still achieve reasonable results. Increasing the number of splits effectively reduces the angular sampling in each subset, leading to more illposed sub-problems. Understanding how far this sparsification can be pushed while maintaining reconstruction quality is an important question for future work. These observations highlight several directions for improving and extending the proposed N2I-LPD framework. In particular, they suggest that the design of the splitting strategy, the role of learned components in the primal and dual spaces, and the trade-off between sparsity and reconstruction quality are key factors that merit further investigation. Overall, the results demonstrate the potential of combining learned iterative reconstruction methods with self-supervised training strategies for challenging CT reconstruction settings.
7
Conclusion
We introduce a novel reconstruction method, Noise2Inverse Learned PrimalDual, which combines the N2I framework with the LPD algorithm, enabling training without ground-truth data. The proposed approach extends learned it-
18
erative reconstruction methods to a self-supervised setting. Our results demonstrate that the N2I-LPD method outperforms classical TV-regularized PrimalDual reconstruction and N2I U-Net in terms of PSNR and visual reconstruction quality, while remaining competitive in SSIM. These findings highlight the potential of combining unrolled learned reconstruction operators with selfsupervised training strategies for low-dose and sparse-angle CT imaging.
Acknowledgments This work has been supported in parts by the Research council of Finland (Project No. 359186, Flagship of Advanced Mathematics for Sensing Imaging and Modelling; Project No. 338408, Academy Research Fellow (AI-SOL); Project No. 370528, Academy project (AequiLoFi)) and by the Finnish Ministry of Education and Culture’s Pilot for Doctoral Programmes (Pilot project Mathematics of Sensing, Imaging and Modelling). S. Rautio also acknowledges support from the Väisälä Fund through the Finnish Academy of Science and Letters.
19
References [1] Jonas Adler, Holger Kohr, and Ozan Öktem. Operator discretization library (odl). Zenodo, 2017. https://github.com/odlgroup/odl. [2] Jonas Adler and Ozan Öktem. Solving ill-posed inverse problems using iterative deep neural networks. Inverse Problems, 33(12):124007, 2017. [3] Jonas Adler and Ozan Öktem. Learned primal-dual reconstruction. IEEE transactions on medical imaging, 37(6):1322–1332, 2018. [4] Simon Arridge, Peter Maass, Ozan Öktem, and Carola-Bibiane Schönlieb. Solving inverse problems using data-driven models. Acta numerica, 28:1– 174, 2019. [5] Joshua Batson and Loic Royer. Noise2self: Blind denoising by selfsupervision. In International Conference on Machine Learning, pages 524– 533. PMLR, 2019. [6] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of mathematical imaging and vision, 40:120–145, 2011. [7] Hu Chen, Yi Zhang, Mannudeep K Kalra, Feng Lin, Yang Chen, Peixi Liao, Jiliu Zhou, and Ge Wang. Low-dose ct with a residual encoderdecoder convolutional neural network. IEEE transactions on medical imaging, 36(12):2524–2535, 2017. [8] Henri Der Sarkissian, Felix Lucka, Maureen van Eijnatten, Giulia Colacicco, Sophia Bethany Coban, and Kees Joost Batenburg. A cone-beam x-ray computed tomography data collection designed for machine learning. Scientific data, 6(1):215, 2019. [9] Heinz Werner Engl, Martin Hanke, and Andreas Neubauer. Regularization of inverse problems, volume 375. Springer Science & Business Media, 1996. [10] Jeffrey A Fessler, Milan Sonka, and J Michael Fitzpatrick. Statistical image reconstruction methods for transmission tomography. Handbook of medical imaging, 2:1–70, 2000. [11] Nadja Gruber, Johannes Schwab, Markus Haltmeier, Ander Biguri, Clemens Dlaska, and Gyeongha Hwang. Noisier2inverse: Self-supervised learning for image reconstruction with correlated noise. IEEE Access, 2025. [12] Javier Gurrola-Ramos, Oscar Dalmau, and Teresa E Alarcón. A residual dense u-net neural network for image denoising. IEEE Access, 9:31742– 31754, 2021. [13] Andreas Hauptmann and Ozan Öktem. Learned iterative networks: An operator learning perspective, 2025. [14] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
20
[15] Allard Adriaan Hendriksen, Daniël Maria Pelt, and K Joost Batenburg. Noise2inverse: Self-supervised deep convolutional denoising for tomography. IEEE Transactions on Computational Imaging, 6:1320–1335, 2020. [16] Kyong Hwan Jin, Michael T McCann, Emmanuel Froustey, and Michael Unser. Deep convolutional neural network for inverse problems in imaging. IEEE transactions on image processing, 26(9):4509–4522, 2017. [17] Avinash C. Kak and Malcolm Slaney. Principles of Computerized Tomographic Imaging. IEEE Press, 1988. [18] Eunhee Kang, Junhong Min, and Jong Chul Ye. A deep convolutional neural network using directional wavelets for low-dose x-ray ct reconstruction. Medical physics, 44(10):e360–e375, 2017. [19] Diederik P Kingma. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014. [20] Alexander Krull, Tim-Oliver Buchholz, and Florian Jug. Noise2voidlearning denoising from single noisy images. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 2129–2137, 2019. [21] Jaakko Lehtinen, Jacob Munkberg, Jon Hasselgren, Samuli Laine, Tero Karras, Miika Aittala, and Timo Aila. Noise2noise: Learning image restoration without clean data. arXiv preprint arXiv:1803.04189, 2018. [22] Ewout van den Berg Michael P. Friedlander. Spot: A linear-operator toolbox for Matlab. https://github.com/mpf/spot, 2013. [Online; accessed 09-April-2026]. [23] Frank Natterer. The Mathematics of Computerized Tomography. SIAM, 1986. [24] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32, 2019. [25] Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015. [26] Leonid I Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: nonlinear phenomena, 60(14):259–268, 1992. [27] Antti Sällinen. Noise2Inverse Learned Primal-Dual. https://github.com/ asalline/N2I-LPD, June 2026. [28] Dirk Elias Schut, Adriaan Graas, Robert van Liere, and Tristan van Leeuwen. Equivariance2inverse: A practical self-supervised ct reconstruction method benchmarked on real, limited-angle, and blurred data. arXiv preprint arXiv:2510.23317, 2025. 21
[29] Emil Y Sidky, Jakob H Jørgensen, and Xiaochuan Pan. Convex optimization problem prototyping for image reconstruction in computed tomography with the chambolle–pock algorithm. Physics in medicine and biology, 57(10):3065–3091, 2012. [30] Wim Van Aarle, Willem Jan Palenstijn, Jeroen Cant, Eline Janssens, Folkert Bleichrodt, Andrei Dabravolski, Jan De Beenhouwer, K Joost Batenburg, and Jan Sijbers. Fast and flexible x-ray tomography using the astra toolbox. Optics express, 24(22):25129–25147, 2016. [31] Wim Van Aarle, Willem Jan Palenstijn, Jan De Beenhouwer, Thomas Altantzis, Sara Bals, K Joost Batenburg, and Jan Sijbers. The astra toolbox: A platform for advanced algorithm development in electron tomography. Ultramicroscopy, 157:35–47, 2015. [32] Ge Wang, Jong Chu Ye, Klaus Mueller, and Jeffrey A Fessler. Image reconstruction is a new frontier of machine learning. IEEE transactions on medical imaging, 37(6):1289–1296, 2018. [33] Ge Wang, Jong Chul Ye, and Bruno De Man. Deep learning for tomographic image reconstruction. Nature machine intelligence, 2(12):737–748, 2020.
22