Paper deep dive
Consecutive Posterior Fusion for Diffusive Recovery of Unobservable Image Structures
Elena Morotti, Davide Evangelista, Elena Loli Piccolomini
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Solving severely ill-posed imaging inverse problems requires recovering image structures that are unobservable or weakly constrained by the measurements. Diffusion models provide expressive learned priors for inferring such missing information, while posterior sampling incorporates measurement consistency along the reverse process. Standard diffusion posterior samplers, however, rely on instantaneous measurement-aware estimates, without explicitly exploiting information carried by previous posterior corrections. We introduce Consecutive Posterior Fusion Denoising Diffusion Null-Space Models (CPF-DDNM), an inference-time strategy that fuses consecutive measurement-aware estimates to improve the diffusive recovery of unobservable image structures, without requiring retraining or additional denoiser evaluations. We instantiate this principle within DDNM, whose range/null-space decomposition reveals that consecutive fusion preserves the measurement-determined component while acting exclusively on the prior-driven null-space estimate. We thus provide a geometric interpretation of CPF-DDNM and a local error analysis that characterizes the optimal time-dependent fusion coefficient, including the extrapolative regime. Experiments on sparse-view and simulated low-dose computed tomography, as well as medical image super-resolution, show consistent improvements over DDNM and competitive performance against diffusion-based inverse solvers.
Tags
Links
- Source: https://arxiv.org/abs/2610.03261v1
- Canonical: https://arxiv.org/abs/2610.03261v1
Trouble viewing inline? Open PDF directly →
Full Text
79,312 characters extracted from source content.
Expand or collapse full text
Consecutive Posterior Fusion for Diffusive Recovery of Unobservable Image Structures Elena Morotti Affiliation: Department of Computer Science and Engineerin Affiliation: University of Bologna Affiliation: Bologna, 40126, Italy Email: elena.morotti4@unibo.it Davide Evangelista Affiliation: Department of Computer Science and Engineering Affiliation: University of Bologna Affiliation: Bologna, 40126, Italy Email: davide.evangelista5@unibo.it Elena Loli Piccolomini Affiliation: Department of Computer Science and Engineering Affiliation: University of Bologna Affiliation: Bologna, 40126, Italy Email: elena.loli@unibo.it Abstract Solving severely ill-posed imaging inverse problems requires recovering image structures that are unobservable or weakly constrained by the measurements. Diffusion models provide expressive learned priors for inferring such missing information, while posterior sampling incorporates measurement consistency along the reverse process. Standard diffusion posterior samplers, however, rely on instantaneous measurement-aware estimates, without explicitly exploiting information carried by previous posterior corrections. We introduce Consecutive Posterior Fusion Denoising Diffusion Null-Space Models (CPF-DDNM), an inference-time strategy that fuses consecutive measurement-aware estimates to improve the diffusive recovery of unobservable image structures, without requiring retraining or additional denoiser evaluations. We instantiate this principle within DDNM, whose range/null-space decomposition reveals that consecutive fusion preserves the measurement-determined component while acting exclusively on the prior-driven null-space estimate. We thus provide a geometric interpretation of CPF-DDNM and a local error analysis that characterizes the optimal time-dependent fusion coefficient, including the extrapolative regime. Experiments on sparse-view and simulated low-dose computed tomography, as well as medical image super-resolution, show consistent improvements over DDNM and competitive performance against diffusion-based inverse solvers. Keywords D iffusion models, medical image reconstruction, inverse problems, null-space methods, posterior fusion 1 Introduction Imaging inverse problems aim to recover an unknown image from indirect measurements related to it through a physical acquisition process. In many practical imaging settings, these measurements may be incomplete, undersampled, or affected by noise, making the reconstruction problem severely ill-posed [2]. Such situations arise, for example, in computed tomography (CT), magnetic resonance imaging (MRI), inpainting, and super-resolution. In these settings, the acquisition process may not uniquely determine the image: undersampling leaves components that are unobservable from the data, while measurement noise further increases the uncertainty along weakly observable directions of the forward operator. The quality of the reconstruction therefore depends critically on the prior used to infer the missing or poorly constrained information. Consider the noiseless linear imaging model = y= K x, where K is the forward operator modelling the imaging system, y contains the acquired data and x is the image to be reconstructed. Whenever K has a non-trivial null-space Ker()Ker( K) and ∈Range() y ( K), the set of data-consistent solutions is the affine space =†+,∈Ker(), x= K y+ z, z ( K), (1) where † K denotes the pseudo-inverse of the forward operator, † K y is the minimum-norm solution determined by the measurements, and z is invisible to the acquisition operator. Infinitely many data-consistent reconstructions can thus be obtained by varying z. In the presence of noise, this ambiguity remains, together with the additional uncertainty associated with small singular values of K. Classical variational methods resolve this ambiguity through handcrafted regularizers that promote properties such as sparsity, smoothness, or structural consistency [6, 11, 12]. Deep learning provides an alternative route by learning image statistics directly from data. In medical reconstruction, learned components have been employed as post-processing operators [13, 23], embedded in unrolled optimization schemes [10, 20], or combined with explicit model-based procedures [24]. More recently, diffusion generative models have emerged as expressive learned priors for inverse problems. Among the broad family of approaches employing diffusion models to drive the image reconstruction problem, in this work we focus on diffusion posterior sampling methods. These methods typically combine a pretrained generative prior with measurement-aware corrections throughout the reverse trajectory, thereby selecting plausible solutions that remain compatible with the observed data. This is done by adjusting the clean-image prediction ^0|t x_0|t obtained by the diffusion model at timestep t using a measurement-aware term, producing a correction t D_t, which is then inserted into the reverse diffusion process in place of the original ^0|t x_0|t. An important representative of the class of diffusion posterior samplers is the Denoising Diffusion Null-Space Model (DDNM) [27], which is particularly relevant because it makes the separation between measured and prior-inferred information explicit by defining its correction as t=†+(−†)^0|t, D_t= K y+( I- K K) x_0|t, (2) where † K y is the component determined by the measurements and (−†)^0|t( I- K K) x_0|t is the null-space completion inferred by the learned prior. This structure is especially relevant for severely undersampled medical imaging problems, such as sparse-view CT and high-factor super-resolution, whose solutions depend strongly on the recovery of unobservable image components. However, a limitation of classical diffusion posterior sampling methods is that, particularly during the early timesteps, when noise dominates the reverse diffusion process, the correction terms t D_t are volatile, sometimes causing the process to diverge towards an out-of-distribution solution. In our recent Lagged Multistep Posterior (LAMP) method [8], we showed that consecutive measurement-aware estimates are highly correlated and may therefore be combined to reduce their variance and stabilize the reconstruction trajectory. While LAMP provides a general mechanism, based on a multistep approximation of an underlying differential equation, for exploiting this interaction, it does not investigate how consecutive posterior corrections interact with the specific structure of the underlying inverse problem. In this work, we address this question within DDNM, whose explicit range/null-space decomposition makes it possible to identify where the use of consecutive posteriors acts. To this aim, we define ~t:=t+βt(t+Δt−t), D_t:= D_t+ _t ( D_t+ t- D_t ), (3) where t+Δt D_t+ t is the correction produced at the preceding reverse iteration and βt _t is the time-dependent fusion coefficient. We refer to this general corrective strategy as Consecutive Posterior Fusion (CPF) and we denote its application to DDNM as CPF-DDNM. In the following, we show that consecutive DDNM corrections t+Δt D_t+ t and t D_t share the same measurement-determined component † K y. Therefore, the CPF correction ~t D_t in Eq. (3) preserves the component determined by the measurements and modifies only the prior-driven recovery of the unobservable image component. The range/null-space structure of DDNM, in addition, provides a direct geometric interpretation of the fusion mechanism and allows the role of the fusion coefficient to be analyzed specifically in terms of null-space recovery, including both interpolative (0<βt<10< _t<1) and extrapolative (βt>1 _t>1) choices. In particular, we investigate time-dependent coefficients that can move beyond the segment connecting consecutive null-space estimates while leaving the measurement-determined component unaffected. The main contributions of this work are summarized as follows: • We specialize consecutive posterior fusion to the range/null-space decomposition of DDNM. We show that the fusion preserves the measurement-determined component and modifies exclusively the completion inferred by the diffusion prior. • We study the extrapolative regime, interpreting this as an over-relaxation within the null-space, and derive a local optimality condition showing when an extrapolative coefficient can minimize the expected completion error. • We introduce time-dependent schedules for βtt∈[0,T]\ _t\_t∈[0,T], and evaluate CPF-DDNM on simulated sparse-view CT, low-dose CT, and medical image super-resolution. CPF-DDNM improves the original DDNM baseline across all experiments, without retraining the diffusion prior or adding denoiser evaluations. The results also demonstrate competitive or superior performance with respect to other state-of-the-art diffusion-based inverse solvers. The remainder of the paper is organized as follows. Section 2 reviews diffusion-based methods for imaging inverse problems. Section 3 introduces CPF-DDNM and presents its geometric interpretation and theoretical analysis. Section 4 describes the experimental protocol, and Section 5 presents and discusses the results. Finally, Section 6 summarizes the main findings and outlines directions for future work. 2 Related Works Diffusion models have been widely adopted as learned priors for imaging inverse problems. Existing approaches differ primarily in how measurement information is incorporated into the reverse process. DDRM [16] derives posterior updates for linear degradation operators through their spectral decomposition. DPS [5] introduces approximate likelihood gradients and can address noisy and nonlinear inverse problems. DDNM [27] enforces measurement consistency through an explicit range/null-space decomposition, whereas DiffPIR [35] embeds a pretrained diffusion denoiser within a plug-and-play half-quadratic splitting scheme. Although these approaches implement data consistency differently, they all generate a sequence of measurement-aware estimates along the reverse trajectory. To speed up the generation process or improve the stability of the diffusion trajectory, several works have recently investigated multistep discretizations of the reverse diffusion dynamics [18, 21, 33]. These approaches exploit information from multiple diffusion steps to construct higher-order approximations of the generative trajectory. By comparison, the temporal relation between consecutive measurement-aware estimates in diffusion posterior samplers remains comparatively underexplored. A recent step in this direction is LAMP [8], which introduces a lagged temporal correction derived from a second-order discretization of the posterior dynamics. LAMP provides initial evidence that consecutive measurement-aware estimates contain useful information that is discarded by standard posterior-sampling updates. Its formulation is backbone-independent and was validated on DiffPIR and DDRM using a linear combination of consecutive estimates. However, due to the second-order discretization formulation, it always resulted in a convex combination of successive measurement-aware terms with coefficients βt∈[0,1] _t∈[0,1]. On the other side, CPF-DDNM focuses on the extrapolative regime, where consecutive estimates define an over-relaxation direction within the null-space. Thus, the two methods exploit the same source of temporal information in different regimes and for different purposes. Diffusion models have also attracted growing interest in medical image reconstruction, where they provide flexible learned priors for highly ill-posed acquisition problems [17, 29]. Applications include sparse-view and low-dose CT, accelerated MRI, PET reconstruction, and multimodal imaging. In CT, diffusion-based methods have been proposed for sparse-view reconstruction and artifact reduction, including approaches that incorporate the acquisition model directly into the generative process [30, 32]. Related posterior-sampling strategies have been developed for accelerated MRI and joint PET–MRI reconstruction [1, 31]. Most of these works focus on modality-specific conditioning mechanisms, data-consistency operators, or posterior-sampling rules. CPF-DDNM is complementary to these developments: it modifies how an existing DDNM solver uses consecutive measurement-aware estimates, without retraining or changing the underlying diffusion prior. 3 Method In this section, we present the proposed posterior fusion strategy for diffusion-based medical image reconstruction. The proposed CPF-DDNM method builds on the null-space formulation of DDNM and replaces the posterior estimate computed at the current reverse iteration with a fusion of two consecutive posterior estimates. 3.1 Diffusion models for image generation Diffusion models constitute one of the most effective paradigms for deep learning-based image generation. In the context of diffusion posterior sampling, we focus on Denoising Diffusion Implicit Models (DDIMs) [25], whose sampling formulation provides the basis for the reconstruction methods considered in this work. We briefly review their generative mechanism to establish the notation used throughout the paper. Let 0 x_0 be a sample from the distribution p0()p_0( x) of clean images. Given two sequences αtt=1T\ _t\_t=1^T and σtt=1T\ _t\_t=1^T satisfying αt2+σt2=1 _t^2+ _t^2=1, respectively referred to as the signal-rate and noise-rate schedules, define the stochastic process tt=1T\ x_t\_t=1^T as t=αt0+σtϵt,ϵt∼(,). x_t= _t x_0+ _t ε_t, ε_t ( 0, I). (4) When αT≈0 _T≈ 0, the terminal variable approximately follows a standard Gaussian distribution, i.e. T∼(,) x_T ( 0, I). A neural network-based denoiser ϵΘ(t,t) ε_ ( x_t,t) is trained to predict the noise component ϵt ε_t from the noisy sample t x_t by minimizing ℒ(Θ)=0∼p0,t∼[0,T],ϵt∼(,)[‖ϵΘ(t,t)−ϵt‖22] ( )=E_ x_0 p_0,t [0,T], ε_t ( 0, I) [|| ε_ ( x_t,t)- ε_t||_2^2 ] (5) Once trained, the denoiser can be used to progressively generate samples from p0()p_0( x) by reversing the process in Eq. (4). In particular, the clean sample can be approximated through the Tweedie estimate ^0|t:=t−σtϵΘ(t,t)αt, x_0|t:= x_t- _t ε_ ( x_t,t) _t, (6) which follows by isolating 0 x_0 in Eq. (4) and replacing the unknown noise ϵt ε_t with its neural-network estimate. The quantity ^0|t x_0|t therefore represents an estimate of the sample 0 x_0 based on the information available at timestep t. As the reverse process progresses from T towards 00, this estimate typically becomes increasingly detailed and coherent. Starting from Eq. (6), DDIM constructs the sample at the next reverse timestep by reintroducing noise according to the lower noise level σt−Δt _t- t: t−Δt=αt−Δt^0|t+σt−ΔtϵΘ(t,t), x_t- t= _t- t x_0|t+ _t- t ε_ ( x_t,t), (7) where Δt t denotes the temporal step size of the reverse process. By repeatedly applying Eq. (6) and Eq. (7), starting from T∼(,) x_T ( 0, I) and proceeding backwards in time until 0 x_0, DDIM generates a new sample from p0()p_0( x). The parameter Δt>0 t>0 determines the number of denoising iterations, since generating 0 x_0 from T x_T requires N=TΔtN= T t reverse steps. A larger number of steps generally increases the computational cost while providing a finer discretization of the reverse trajectory and, potentially, a more accurate generated sample. 3.2 Diffusion posterior samplers The generative process described in Section 3.1 can produce high-quality samples, although it does not directly control which image is generated. Solving an inverse problem instead requires the generated image to remain compatible with the observed measurements. Diffusion posterior samplers address this requirement by introducing a measurement-aware correction at each reverse timestep, immediately after computing the clean estimate ^0|t x_0|t. This correction progressively guides the reverse trajectory towards consistency with the measurements, while the diffusion prior complements the available information by estimating the unobservable component ∈Ker() z ( K) introduced in Eq. (1). Consider a pre-trained diffusion model ϵΘ(t,t) ε_ ( x_t,t), trained in an unsupervised manner on a dataset representative of the images to be processed. Given an initial iterate T∼(,) x_T ( 0, I), the clean estimate ^0|t x_0|t is computed at each timestep according to Eq. (6). Then, a measurement-aware correction operator (^0|t,t,,)C( x_0|t,t, K, y) is introduced, whose formulation depends on the selected posterior-sampling method. We call the corrected image t D_t, i.e. t:=(^0|t,t,,) D_t:=C( x_0|t,t, K, y). Replacing ^0|t x_0|t with t D_t in Eq. (7) yields the diffusion posterior-sampling update t−Δt=αt−Δtt+σt−ΔtϵΘ(t,t), x_t- t= _t- t D_t+ _t- t ε_ ( x_t,t), (8) through which the diffusion trajectory is progressively conditioned on the acquired measurements. Iterating Eq. (8) until t=Δt= t yields the final reconstruction 0 x_0. Several posterior samplers can be interpreted within this measurement-correction framework, and the specific choice of (^0|t,t,,)C( x_0|t,t, K, y) determines the diffusion posterior-sampling method. In addition, the precise point at which the correction is introduced further differs across methods. Focusing on the DDNM method proposed in [27], it adopts the correction described in Eq. (2). Denoting the projection operator onto the kernel of K as Ker:=−†P_Ker:= I- K K [6], the DDNM correction reads: t=†+Ker^0|t. D_t= K y+P_Ker x_0|t. (9) For the noiseless setting, this expression provides a clear geometric interpretation of the correction of ^0|t x_0|t: the DDNM preserves only the estimated null-space component Ker^0|tP_Ker x_0|t, while it replaces the range-space component of ^0|t x_0|t with the measurement-determined term † K y. When † K is the exact Moore-Penrose pseudo-inverse and the measurements belong to the range of K, this construction enforces exact data consistency at every corrected estimate. Additionally, this explicit separation also makes the role of the diffusion model particularly transparent, since the learned prior affects only the image components that cannot be recovered directly from the measurements. Applying the same correction directly to noisy measurements may, however, amplify noise in the reconstruction. In an ill-conditioned problem, the pseudo-inverse strongly amplifies measurement components associated with small singular values of K, and the term † K y may consequently become unstable. 3.3 The CPF-DDNM method Although the reverse diffusion process in Eq. (8) generates a sequence of strongly correlated estimates of the same unknown image, standard DDNM constructs its correction using exclusively the denoised estimate produced at the current timestep. The estimate obtained at the preceding reverse iteration is therefore discarded, even though it provides a closely related prediction of the same unobservable image component. Motivated by this observation, we propose to combine the DDNM corrections obtained at two consecutive reverse iterations. As anticipated in Eq. (3), the resulting fused correction is defined as ~t=t+βt(t+Δt−t) D_t= D_t+ _t( D_t+ t- D_t), where βt∈ℝ _t controls the contribution of the preceding correction. By expanding the two DDNM terms, this expression becomes ~t=†+(−†)(^0|t+βt(^0|t+Δt−^0|t))=†+Ker(^0|t+βt(^0|t+Δt−^0|t)), split D_t&= K y+ ( I- K K ) ( x_0|t+ _t ( x_0 t+ t- x_0|t ) )\\ &= K y+P_Ker ( x_0|t+ _t ( x_0 t+ t- x_0|t ) ), split (10) where Ker:=−†P_Ker:= I- K K. The proposed reverse step is then obtained by replacing the standard DDNM correction t D_t with its fused counterpart ~t D_t in Eq. (8), yielding t−Δt=αt−Δt~t+σt−Δtϵθ(t,t). x_t- t= _t- t D_t+ _t- t _θ( x_t,t). (11) This extension of DDNM, which we refer to as Consecutive Posterior Fusion DDNM (CPF-DDNM), requires only the storage of the correction obtained at the preceding reverse iteration. It introduces no additional denoiser evaluations and requires no modification or retraining of the pretrained diffusion model. Equation (10) also shows that CPF-DDNM retains the range/null-space structure of the original DDNM correction. To make this property explicit, let t:=Ker^0|t z_t:=P_Ker x_0|t (12) denote the null-space completion predicted at reverse timestep t. The standard DDNM correction can then be written as t=†+t D_t= K y+ z_t. Similarly, the correction obtained at the preceding reverse iteration is t+Δt=†+t+Δt D_t+ t= K y+ z_t+ t, and the CPF-DDNM correction admits the decomposition ~t=†+~t,~t:=t+βt(t+Δt−t). D_t= K y+ z_t, z_t:= z_t+ _t ( z_t+ t- z_t ). (13) Hence, the measurement-determined component † K y remains unchanged, while the fusion acts exclusively on the prior-driven null-space completion. The following proposition formalizes this property and establishes the corresponding data-consistency result. Proposition 1 (Null-space preservation and data consistency). Assume that † K is the Moore-Penrose pseudo-inverse of K, so that Ker=−†P_Ker= I- K K is the orthogonal projector onto Ker()Ker( K). Then, for any βt∈ℝ _t , the CPF-DDNM correction satisfies ~t=†+~t,~t∈Ker(), D_t= K y+ z_t, z_t ( K), (14) and consequently ~t=†. K D_t= K K y. (15) If ∈Range() y ( K), then ~t= K D_t= y, and the fused correction is exactly data-consistent. Proof. Substituting t=†+t D_t= K y+ z_t and t+Δt=†+t+Δt D_t+ t= K y+ z_t+ t into Eq. (3) gives ~t=†+t+βt(t+Δt−t)=†+~t. D_t= K y+ z_t+ _t ( z_t+ t- z_t )= K y+ z_t. (16) Since t,t+Δt∈Ker() z_t, z_t+ t ( K) and Ker()Ker( K) is a linear subspace, their affine combination ~t=t+βt(t+Δt−t) z_t= z_t+ _t( z_t+ t- z_t) also belongs to Ker()Ker( K). It follows that ~t= K z_t= 0 and therefore ~t=† K D_t= K K y. When ∈Range() y ( K), the Moore-Penrose identity †= K K y= y establishes exact data consistency. ∎ Proposition 1 assumes access to the exact Moore-Penrose pseudo-inverse. In the CT experiments considered in this work, the required pseudo-inverse operations are approximated through a finite number of CGLS iterations. Consequently, the corresponding null-space projection and data-consistency properties hold up to the numerical error introduced by this approximation. 3.4 Insights on the fusion parameter The role of βt _t can now be interpreted geometrically within the null-space. We define the displacement between two consecutive null-space estimates as t:=t+Δt−t. d_t:= z_t+ t- z_t. (17) Using this notation, Eq. (13) becomes ~t=t+βtt. z_t= z_t+ _t d_t. (18) Hence, the pair (t,t+Δt)( z_t, z_t+ t) determines a local direction t d_t within the null-space, while βt _t determines the position of the fused estimate along the affine line passing through these two estimates. Since the measurement-dependent component † K y is unchanged, this displacement affects exclusively the component that is not determined by the measurements. For 0≤βt≤10≤ _t≤ 1, the fused estimate lies on the segment connecting t z_t and t+Δt z_t+ t, and the fusion is therefore interpolative. In particular, βt=0 _t=0 recovers the current estimate t z_t, whereas βt=1 _t=1 selects the estimate t+Δt z_t+ t obtained at the preceding reverse iteration. When βt>1 _t>1, the fused estimate lies beyond t+Δt z_t+ t along the direction t d_t, and the update becomes extrapolative. In this regime, Eq. (18) can equivalently be written as ~t=t+Δt+t, z_t= z_t+ t+ m_t, (19) where t:=(βt−1)t=(βt−1)(t+Δt−t). m_t:=( _t-1) d_t=( _t-1)( z_t+ t- z_t). (20) The quantity t m_t represents the additional displacement introduced beyond t+Δt z_t+ t and can therefore be interpreted as an over-relaxation term acting exclusively within Ker()Ker( K). It also has a momentum-like algebraic form, since it extends the direction identified by the difference between two consecutive posterior estimates. The geometric interpretation of this extrapolative regime is illustrated in Figure 1. Figure 1: Geometric interpretation of the CPF-DDNM fusion in the extrapolative regime βt>1 _t>1. The current and preceding null-space estimates Ker^0|tP_Ker x_0|t and Ker^0|t+ΔtP_Ker x_0|t+ t determine the direction t:=Ker^0|t+Δt−Ker^0|t d_t:=P_Ker x_0|t+ t-P_Ker x_0|t. CPF-DDNM extends the estimate beyond Ker^0|t+ΔtP_Ker x_0|t+ t through the additional contribution t:=(βt−1)t m_t:=( _t-1) d_t, while the measurement-determined component † K y remains unchanged. The geometric representation above describes how the fused estimate moves as βt _t varies, while the potential advantage of this operation can be analyzed through a local uncertainty model. Let ∗ z^* denote the ideal null-space component associated with the unknown image, and represent two consecutive estimates as t=∗+t,t+Δt=∗+t+Δt, z_t= z^*+ η_t, z_t+ t= z^*+ η_t+ t, (21) where t η_t and t+Δt η_t+ t denote the errors associated with the corresponding prior-driven estimates. In the CT experiments, these terms may additionally include perturbations introduced by the numerical approximation of the pseudo-inverse. Substituting Eq. (21) into Eq. (13) gives ~t−∗=(1−βt)t+βtt+Δt. z_t- z^*=(1- _t) η_t+ _t η_t+ t. (22) CPF-DDNM can therefore be interpreted as combining two correlated estimates of the same unknown null-space component. In particular, the error of the fused estimate is itself a linear combination of the errors associated with the two consecutive posterior estimates, suggesting that an appropriate choice of βt _t can reduce the expected reconstruction error. For a fixed reverse timestep t, we formalize this observation by defining Jt(β):=[‖(1−β)t+βt+Δt‖22],J_t(β):=E [ \|(1-β) η_t+β η_t+ t \|_2^2 ], (23) where the expectation is taken over the randomness inducing the reconstruction errors at the considered timestep. Let At:=[∥t∥22],Bt:=[∥t+Δt∥22],Ct:=[⟨t,t+Δt⟩]. gatheredA_t:=E [\| η_t\|_2^2 ], B_t:=E [\| η_t+ t\|_2^2 ],\\ C_t:=E [ η_t, η_t+ t ]. gathered (24) The following result characterizes the fusion coefficient minimizing the expected error under this local model. Proposition 2 (Optimal fusion coefficient). Assume that the consecutive null-space estimates satisfy Eq. (21). If At+Bt−2Ct>0A_t+B_t-2C_t>0, then Jt(β)J_t(β) is strictly convex and admits the unique minimizer βt∗=At−CtAt+Bt−2Ct. _t^*= A_t-C_tA_t+B_t-2C_t. (25) Proof. Expanding Eq. (23) gives Jt(β)=(1−β)2At+β2Bt+2β(1−β)Ct.J_t(β)=(1-β)^2A_t+β^2B_t+2β(1-β)C_t. (26) Differentiating with respect to β yields ∂Jt∂β=2(Ct−At)+2β(At+Bt−2Ct). ∂ J_t∂β=2(C_t-A_t)+2β(A_t+B_t-2C_t). (27) Setting the derivative equal to zero gives Eq. (25). Moreover, ∂2Jt∂β2=2(At+Bt−2Ct)>0, ∂^2J_t∂β^2=2(A_t+B_t-2C_t)>0, (28) which proves that βt∗ _t^* is the unique minimizer. ∎ Proposition 2 shows that the optimal fusion coefficient depends on both the error magnitude of the consecutive estimates and their correlation. Its value is therefore not restricted to the interpolation interval [0,1][0,1]. In particular, under the condition of Proposition 2, βt∗>1⟺Ct>Bt, _t^*>1 C_t>B_t, (29) showing that an extrapolative update may minimize the expected error when the cross-correlation between the current and preceding errors is sufficiently large relative to the error energy of the preceding estimate. This result provides a statistical counterpart to the geometric interpretation in Figure 1, while avoiding the need to introduce a separate realization-wise optimization criterion. The statistics AtA_t, BtB_t, and CtC_t depend on the unknown reconstruction errors and are therefore unavailable during inference. Rather than attempting to evaluate Eq. (25) explicitly, we parametrize βt _t through a task-dependent schedule of the form βt=c1+c21−tT, _t=c_1+c_2 1- tT, (30) where c1c_1 and c2c_2 control the evolution of the fusion strength along the reverse trajectory. In particular, βT=c1,β0=c1+c2. _T=c_1, _0=c_1+c_2. (31) Different choices of c1c_1 and c2c_2 therefore allow the fusion strength to increase or decrease during sampling, to vanish as the reverse process approaches the clean-image regime, or to remain active throughout the trajectory, including in the extrapolative regime. The schedules considered in this work and their effect on reconstruction performance are investigated experimentally in Section 5.2. 4 Experimental setup In this section, we describe the implementation details of the proposed method and the medical imaging applications whose results are presented and discussed in the following section. All experiments were implemented in Python 3.12.12 using PyTorch v2.7.1, MONAI v1.5.1, and ASTRA Toolbox v2.4.1. Training and inference were performed on an NVIDIA A4000 GPU with 24 GB of memory. Unless otherwise specified, all experiments used the random seed 42. The source code, configuration files, pretrained checkpoint, and patient-level data partitions will be made available at code-released-upon-acceptance. 4.1 Dataset and preprocessing The dataset used is the Mayo Clinic dataset [22], consisting of 3631 images with 512×512512× 512 pixels of human chests from 11 patients, which we have downscaled to 256×256256× 256 pixels and normalized in the range [0,1][0,1]. In particular, we selected 3305 images from 10 patients for training and the remaining 326 from a single patient for testing. Figure 2 shows representative CT images from the test set, selected to illustrate different anatomical regions and structural characteristics. The corresponding enlarged regions displayed in the second row highlight fine anatomical details, including soft-tissue boundaries, pulmonary structures, and bone interfaces. Images belonging to the held-out test set are used as ground-truth references. Figure 2: Representative axial CT images from the Mayo Clinic dataset used as ground-truth references. The second row shows zoomed regions of interest. 4.2 Diffusion model and training As a diffusion prior for all the experiments, we consider the one released in [7], which was trained from scratch on the Mayo Clinic images. In particular, we employed a two-dimensional DiffusionModelUNet architecture provided by the MONAI Generative Models framework [3]. The network gets as input a batch of grey-scale 256×256256× 256 images and processes it via a four-level UNet architecture with 6464, 128128, 256256, and 512512 feature channels, respectively. Each resolution level contains two residual blocks. As is common practice, we employed self-attention only at the two deepest resolution levels, where the spatial dimensionality is smaller and long-range dependencies can be modeled at a reduced computational cost. The attention layers use 3232 channels per attention head. The network receives the noisy image t x_t together with the corresponding diffusion timestep t, embedded through a sinusoidal embedding layer, and predicts a single-channel estimate of the noise ϵΘ(t,t) ε_ ( x_t,t) appearing in Eq. (5). Training was performed using a discrete DDPM forward process with T=1000T=1000 diffusion timesteps. We adopted a linear variance schedule, with β1=10−4 _1=10^-4 and βT=2×10−2 _T=2× 10^-2. At every training iteration, a clean image 0 x_0 and a timestep t were sampled, with t drawn uniformly from the 10001000 available diffusion steps. Before applying the forward diffusion process, the image intensities, originally normalized to [0,1][0,1], were linearly mapped to [−1,1][-1,1]. For any timestep t, a Gaussian noise realization ϵt∼(,) ε_t ( 0, I) was then generated and added to the clean image according to the noise level associated with the sampled timestep. The resulting noisy image was provided to the U-Net, which was trained to recover the injected noise by minimizing the mean squared error ℒDDPM(Θ)=0,t,ϵt[‖ϵΘ(t,t)−ϵt‖22],L_DDPM( )=E_ x_0,t, ε_t [|| ε_ ( x_t,t)- ε_t||_2^2 ], (32) consistently with the noise-prediction objective introduced in Eq. (5). The model parameters were optimized using Adam with a learning rate of 10−410^-4 and a batch size of 1616, for a total of 100100 epochs. 4.3 Medical applications and imaging tasks We evaluate CPF-DDNM on two medical imaging inverse problems: computed tomography reconstruction and single-image super-resolution. For CT, we consider both sparse-view and low-dose acquisition settings. Although these problems differ in the way the measurements are degraded, they all contain image components that are either unobservable or only weakly constrained by the acquired data, and therefore require prior information to obtain a plausible reconstruction. 4.3.1 CT image reconstruction Computed tomography reconstructs a cross-sectional image from projection measurements acquired at different viewing angles. In the two-dimensional parallel-beam geometry considered in this work, each projection angle defines a family of parallel X-ray paths through the object, while the detector records the corresponding line-integral measurements. The collection of all projections acquired over the angular range forms the sinogram. Let ∈ℝM2 x ^M^2 denote a ground-truth image and let K be the discretized parallel-beam projection operator. For each image, the simulated measurements are generated according to =+δ‖‖∞,∼(,), y= K x+δ \| K x \|_∞ e, e ( 0, I), (33) where δ controls the relative noise level. Both K and ⊤ K were implemented using the AstraToolbox [26]. For an M×M× M image acquired using NθN_θ projection angles and a detector with NdN_d elements, the corresponding discretized operator satisfies ∈ℝNθ⋅Nd×M2 K ^N_θ· N_d× M^2. In all our CT experiments, we used M=256M=256 and Nd=512N_d=512. We considered two acquisition settings, corresponding to sparse-view and low-dose CT. In the sparse-view experiments, the number of acquired measurements is smaller than the number of unknown image pixels. In particular, for all the geometries considered in this work, NθNd<M2N_θN_d<M^2, which implies that the system is underdetermined and that K necessarily has a non-trivial null-space. As the number of projection angles decreases, the amount of available information is further reduced, and the unresolved image subspace generally becomes larger. We test three acquisition geometries, denoted by G60G_60, G45G_45, and G30G_30, corresponding respectively to Nθ∈60,45,30N_θ∈\60,45,30\ projection angles uniformly distributed over [0,180∘)[0,180 ). A relative noise level of δ=0.01δ=0.01 was used in all sparse-view experiments. For the low-dose experiments, the parallel-beam geometry was instead fixed to Nθ=180N_θ=180 projection angles uniformly distributed over [0,180∘)[0,180 ), so that the number and spatial arrangement of the measurements remained unchanged across the experiments. In actual low-dose CT acquisitions, the reduced photon count increases measurement uncertainty and decreases the signal-to-noise ratio of the projection data. We approximate this effect using the additive post-log Gaussian noise model in Eq. (33), considering relative noise levels δ∈0.1,0.3,0.5δ∈\0.1,0.3,0.5\. This model should be interpreted as a simplified approximation of the quantum and electronic noise affecting CT measurements after logarithmic preprocessing, rather than as a complete physical simulation of the photon-counting acquisition process. Although the system matrix in this setting has more rows than columns, the nominal number of measurements NθNdN_θN_d overestimates the amount of useful information available to the reconstruction, since detector bins corresponding to rays that do not intersect the image support produce zero measurements. For a square image with side length M, the projected width of the support at angle θ is approximately M(|cosθ|+|sinθ|)M(| θ|+| θ|). Consequently, for uniformly distributed projection angles, the effective number of nonzero measurements can be approximated as meff≈M∑k=0Nθ−1(|cosθk|+|sinθk|)≈4πMNθ.m_eff≈ M _k=0^N_θ-1 (| _k|+| _k| )≈ 4πMN_θ. (34) A complete proof of this result can be found in [14]. Since rank()≤meffrank( K)≤ m_eff, the condition meff<M2m_eff<M^2 is sufficient to imply rank()<M2rank( K)<M^2 and therefore Ker()≠Ker( K)≠\ 0\. Using the approximation in Eq. (34), this condition reduces to Nθ<π4MN_θ< π4M. In our experiments, M=256M=256 and Nθ=180N_θ=180, while π4M≈201.1 π4M≈ 201.1. According to this support-based dimensional estimate, the discretized projection operator considered in the low-dose experiments is therefore expected to have a non-trivial kernel. Independently of this exact null-space component, the added measurement noise also introduces substantial uncertainty along poorly conditioned and weakly observable directions of the CT operator. 4.3.2 Medical image super-resolution Single-image super-resolution aims to reconstruct a high-resolution image from a spatially downsampled observation. In the setting considered in this work, each low-resolution pixel is obtained by averaging the intensities over a non-overlapping block of high-resolution pixels. This operation removes fine-scale anatomical information, including small structures and sharp tissue boundaries, which cannot be uniquely recovered from the observed image alone. The acquisition model is written as =+ y= K x+ e, where ∈ℝM2 x ^M^2 is the unknown high-resolution image, K denotes the block-averaging downsampling operator, and e represents additive measurement noise defined as in Eq. (33). We considered both the noiseless case and a mildly noisy setting, with δ∈0, 0.01δ∈\0,\ 0.01\. For an integer scale factor s>1s>1, the observed image contains (M/s)2(M/s)^2 pixels and the corresponding forward operator satisfies ∈ℝ(M/s)2×M2 K ^(M/s)^2× M^2. Since (M/s)2<M2(M/s)^2<M^2, the resulting inverse problem is necessarily underdetermined and K has a non-trivial null-space. More precisely, dimKer()=M2−rank()≥M2−(Ms)2. ( K)=M^2-rank( K)≥ M^2- ( Ms )^2. (35) Multiple high-resolution images can therefore produce the same low-resolution observation, and the fine-scale image components removed by the downsampling operation must be inferred using prior information. This null-space ambiguity makes super-resolution a natural application for DDNM-based reconstruction and for the proposed CPF-DDNM strategy. In all super-resolution experiments, we used a scale factor s=4s=4. 4.4 Pseudo-inverse computation and algorithmic implementation The proposed method requires the application of the pseudo-inverse † K as in (10). The way † K is evaluated depends on the structure of the forward operator, and this distinction is important for the inverse problems considered in this work. For super-resolution, the degradation operator is an average downsampling operator acting blockwise. In this case, the Singular Value Decomposition (SVD) of K is available in closed form, since each low-resolution pixel corresponds to the average of a non-overlapping high-resolution block. Therefore, the pseudo-inverse can be computed exactly by applying the corresponding blockwise SVD. This allows us to implement the DDNM correction with an exact † K in the super-resolution experiments. For computed tomography, the forward operator K maps an image into its corresponding sinogram. At the image sizes considered in this work, the matrix representation of K is too large to be explicitly formed and factorized, and therefore the action of the Moore-Penrose pseudoinverse cannot be computed through an explicit SVD. Therefore, rather than exactly constructing † K , we approximate its action numerically using the conjugate-gradient least-squares (CGLS) algorithm initialized from the zero image. CGLS requires only applications of the forward operator K and its adjoint ⊤ K , and is therefore well suited to large-scale tomographic reconstruction problems. For a fixed number J of iterations, we denote the resulting approximate inverse action by CGLSJ(,)≃~J†,CGLS_J( K, w) K _J w, where w is a generic measurement-space vector. Since to compute ~t D_t we only need † K y and (−†)^0|t ( I- K K ) x_0|t, the DDNM correction can be implemented without explicitly constructing the projector † K . In the CT experiments, we replace the exact pseudo-inverse action with its finite-iteration CGLS approximation and use t=CGLSJ(,)+CGLSJ(,^0|t),J=16. D_t= CGLS_J( K, y)+ CGLS_J( K, K x_0|t), J=16. (36) Hence, the corrected estimates may be affected by two distinct sources of approximation: the error in the prior-driven prediction produced by the diffusion model and the numerical error associated with the finite-iteration approximation of † K . Unlike the exact-pseudoinverse setting considered in the theoretical analysis, data consistency and the range/null-space decomposition are therefore satisfied only up to the accuracy of the CGLS approximation. We observe that finite-iteration CGLS introduces an implicit regularization effect, which can be relevant in the presence of measurement noise. The complete reconstruction procedure is summarized in Algorithm 1. Algorithm 1 CPF-DDNM posterior sampling 0: Measurement y, forward operator K, pretrained denoiser ϵθ ε_θ, fusion schedule β(t)β(t), number of CGLS iterations J∈ℕJ 1: Initialize T∼(,) x_T ( 0, I) 2: Compute †=CGLSJ(,) K y= CGLS_J( K, y) 3: Set prev=None D_prev= None 4: for t=T,T−Δt,…,Δt=T,T- t,…, t do 5: Compute ^0|t=(t−σtϵθ(t,t))/αt x_0|t= ( x_t- _t ε_θ( x_t,t) )/ _t 6: Compute †^0|t=CGLSJ(,^0|t) K K x_0|t= CGLS_J\! ( K, K x_0|t ) 7: Compute t=†+^0|t−†^0|t D_t= K y+ x_0|t- K K x_0|t as in (9) 8: if prev≠None D_prev≠ None then 9: Set βt=β(t) _t=β(t) 10: Compute ~t=t+βt(prev−t) D_t= D_t+ _t ( D_prev- D_t ) as in (3) 11: else 12: Set ~t=t D_t= D_t 13: end if 14: Update t−Δt=αt−Δt~t+σt−Δtϵθ(t,t) x_t- t= _t- t D_t+ _t- t ε_θ( x_t,t) 15: Set prev=t D_prev= D_t 16: end for 17: return 0 x_0 4.5 Comparison methods, parameter settings, and evaluation metrics The proposed CPF-DDNM method was compared with DDNM [27], DiffPIR [35], Diffusion Posterior Sampling (DPS) [5], and its Posterior Sampling Plus variant (PS+) [5]. DDNM represents the main baseline, since CPF-DDNM directly modifies its null-space completion strategy. For the sparse-view and low-dose CT experiments, we additionally considered the analytical filtered backprojection (FBP) method [15] and two model-based iterative reconstruction methods: the Simultaneous Iterative Reconstruction Technique (SIRT) [9] and an ℓ2 _2-TV regularized least-squares method solved using the Chambolle-Pock algorithm [4, 19]. The parameters of the comparison methods were tuned on the training set and then kept fixed across all test images within the same task and acquisition setting. For DDNM and CPF-DDNM, we set Δt=50 t=50 in all experiments which, together with T=1000T=1000, results in N=20N=20 reverse diffusion steps. In the CT experiments, the action of the pseudo-inverse † K was approximated using CGLS with J=16J=16 iterations, a setting that provided stable performance across the considered acquisition conditions. Note that limiting the number of CGLS iterations introduces an implicit regularization effect, since early termination reduces the amplification of measurement noise along poorly determined components of the inverse problem. This truncated CGLS approximation was therefore used for both DDNM and CPF-DDNM as an alternative to the explicit noise-handling strategy introduced by DDNM+. The choice of the fusion coefficient βt _t is a crucial aspect of CPF-DDNM. Accordingly, we conducted an ablation study on sparse-view CT, reported in Section 5.2, considering both constant coefficients and time-dependent schedules. The experimental analysis indicates that the CT reconstruction tasks benefit from maintaining βt>1 _t>1 throughout the reverse trajectory, thereby preserving the extrapolative effect of the fusion. For all CT experiments, the fusion coefficient was defined according to the time-dependent schedule βt=1.8−0.61−tT, _t=1.8-0.6 1- tT, (37) which keeps βt _t in the extrapolative regime throughout the reverse diffusion process, with βt∈[1.2,1.8] _t∈[1.2,1.8]. For super-resolution, the most effective schedule instead gradually decreases toward zero during the reverse process and is defined as βt=1−1−tT. _t=1- 1- tT. (38) Reconstruction quality was quantitatively evaluated using three complementary metrics: SSIM [28], which measures structural similarity; PSNR, which quantifies pixel-wise reconstruction fidelity; and LPIPS [34], which measures perceptual similarity in a learned feature space. Each reconstruction was repeated using five random initializations. For each test image, the metrics were first averaged across the five runs, after which the reported mean and standard deviation were computed across the test images. 5 Results and discussion In this Section we present and discuss the results obtained on the experiments described in Section 4. 5.1 Results on sparse-view CT reconstruction The first numerical experiment addresses sparse-view CT reconstruction. We evaluate the considered methods under increasingly severe angular undersampling by using 60, 45, and 30 projection angles. The quantitative results are summarized in Table 1. Table 1: Quantitative evaluation metrics computed over the entire test set across five runs with different initial noise realizations, for acquisition geometries with 30, 45, and 60 projection angles. In bold, the best results for each experiment. Method SSIM (↑ ) LPIPS (↓ ) PSNR (↑ ) G30G_30 FBP 0.407±0.0140.407± 0.014 0.442±0.0160.442± 0.016 22.83±0.4622.83± 0.46 SIRT 0.636±0.0210.636± 0.021 0.384±0.0140.384± 0.014 26.52±0.5526.52± 0.55 ℓ2 _2-TV 0.828±0.0090.828± 0.009 0.132±0.0150.132± 0.015 30.61±0.4830.61± 0.48 DPS 0.756±0.0180.756± 0.018 0.136±0.0240.136± 0.024 26.62±0.7226.62± 0.72 PS+ 0.524±0.0180.524± 0.018 0.231±0.0130.231± 0.013 19.96±0.5219.96± 0.52 DiffPIR 0.851±0.0120.851± 0.012 0.085±0.0110.085± 0.011 32.12±0.6232.12± 0.62 DDNM 0.784±0.0180.784± 0.018 0.132±0.0110.132± 0.011 31.06±0.6231.06± 0.62 CPF-DDNM 0.851±0.009 0.851± 0.009 0.083±0.015 0.083± 0.015 32.40±0.52 32.40± 0.52 G45G_45 FBP 0.518±0.0180.518± 0.018 0.364±0.0120.364± 0.012 25.84±0.4925.84± 0.49 SIRT 0.705±0.0160.705± 0.016 0.322±0.0110.322± 0.011 28.19±0.5328.19± 0.53 ℓ2 _2-TV 0.860±0.0060.860± 0.006 0.092±0.0130.092± 0.013 32.62±0.4532.62± 0.45 DPS 0.759±0.0170.759± 0.017 0.135±0.0220.135± 0.022 26.74±0.6626.74± 0.66 PS+ 0.525±0.0180.525± 0.018 0.231±0.0130.231± 0.013 19.98±0.5319.98± 0.53 DiffPIR 0.888±0.0110.888± 0.011 0.060±0.0090.060± 0.009 34.10±0.6334.10± 0.63 DDNM 0.824±0.0160.824± 0.016 0.106±0.0110.106± 0.011 32.67±0.6132.67± 0.61 CPF-DDNM 0.913±0.006 0.913± 0.006 0.036±0.003 0.036± 0.003 34.86±0.49 34.86± 0.49 G60G_60 FBP 0.602±0.0190.602± 0.019 0.317±0.0120.317± 0.012 27.63±0.5027.63± 0.50 SIRT 0.758±0.0150.758± 0.015 0.264±0.0100.264± 0.010 29.55±0.5429.55± 0.54 ℓ2 _2-TV 0.880±0.0060.880± 0.006 0.076±0.0100.076± 0.010 33.72±0.4533.72± 0.45 DPS 0.759±0.0160.759± 0.016 0.135±0.0210.135± 0.021 26.74±0.6126.74± 0.61 PS+ 0.525±0.0180.525± 0.018 0.230±0.0130.230± 0.013 19.99±0.5319.99± 0.53 DiffPIR 0.912±0.0090.912± 0.009 0.043±0.0070.043± 0.007 35.61±0.6335.61± 0.63 DDNM 0.853±0.0140.853± 0.014 0.084±0.0100.084± 0.010 33.89±0.6233.89± 0.62 CPF-DDNM 0.938±0.006 0.938± 0.006 0.021±0.003 0.021± 0.003 37.27±0.66 37.27± 0.66 As expected, the performance of most reconstruction methods improves as the number of acquired views increases. CPF-DDNM achieves the best mean values for all three evaluation metrics and for all the considered acquisition geometries, yielding the highest SSIM and PSNR values and the lowest LPIPS scores. For the most challenging 30-view geometry, the CPF-DDNM performance is very close to that of DiffPIR. Nevertheless, CPF-DDNM provides a clear improvement over the original DDNM baseline, increasing the SSIM from 0.78370.7837 to 0.85110.8511 and the PSNR from 31.0631.06 to 32.40dB32.40\,dB, while reducing the LPIPS from 0.13220.1322 to 0.08290.0829. For the 45- and 60-view geometries, the advantage of CPF-DDNM over the competing diffusion-based methods becomes more pronounced. With 45 views, CPF-DDNM improves the PSNR by approximately 0.77dB0.77\,dB over DiffPIR and reduces the LPIPS from 0.0601 to 0.0356, while also achieving a higher SSIM (0.9125 versus 0.8877). The improvement is even more evident with 60 views, where CPF-DDNM reaches a PSNR of 37.27dB37.27\,dB, approximately 1.66dB1.66\,dB higher than DiffPIR, together with an LPIPS reduction from 0.0433 to 0.0208. Compared with the original DDNM baseline, the proposed fusion consistently improves all three metrics in both acquisition settings. Overall, the increasing performance gap observed from 45 to 60 views is consistent with the hypothesis that consecutive posterior fusion can more effectively exploit the additional measurement information available at denser angular samplings. FBP SIRT ℓ2 _2-TV DPS PS+ DiffPIR DDNM CPF-DDNM Figure 3: Qualitative comparison of sparse-view CT reconstructions obtained with the different methods. From top to bottom, the results correspond to 30, 45, and 60 projection angles. For each acquisition setting, the second row shows magnified details of the corresponding reconstructions. Figure 3 provides a qualitative comparison of the reconstruction methods for the 30-, 45-, and 60-view geometries, shown from top to bottom, respectively. As expected, the reconstruction quality progressively improves as the number of projection angles increases, with better-defined anatomical structures and a reduction of streaking and blurring artifacts. For the most challenging 30-view geometry, FBP and SIRT are strongly affected by reconstruction artifacts, whereas ℓ2 _2-TV suppresses most of the streaks but produces an overly smoothed appearance. The diffusion-based methods provide more anatomically plausible reconstructions, although DPS and PS+ show limited structural fidelity. DiffPIR, DDNM, and CPF-DDNM recover the main anatomical structures more accurately, with CPF-DDNM providing slightly sharper boundaries and better preservation of fine details. The differences among the methods become clearer for the 45- and 60-view geometries. As more projection information becomes available, CPF-DDNM more effectively combines measurement consistency with the learned prior. These qualitative observations are consistent with the quantitative results reported in Table 1. Overall, these results confirm the effectiveness of consecutive posterior fusion in improving the reconstruction quality of the original DDNM baseline. FBP ℓ2 _2-TV DiffPIR DDNM CPF-DDNM Figure 4: Quantitative and qualitative comparison for the 45-view sparse-CT geometry. The boxplots show the distributions of SSIM and LPIPS values computed over the entire test set across five runs with different initial noise realizations. Representative reconstructions and enlarged regions of interest are shown for visual comparison. Figure 4 provides a more detailed analysis of the 45-view geometry and is divided into two complementary parts. The upper part reports the boxplots of the SSIM and LPIPS values computed over the entire test set and across five runs with different initial noise realizations. CPF-DDNM achieves the best overall distributions, with higher SSIM values and lower LPIPS values than the competing methods. In particular, the improvement over DDNM is clearly visible for both metrics. The lower part presents a further qualitative comparison on a test image different from the one shown in Figure 3. The corresponding enlarged regions of interest highlight the differences in the recovery of the low-contrast area by the considered methods. 5.2 Ablation Study on the Fusion Coefficient We investigated the influence of the fusion coefficient βt _t on the reconstruction dynamics by considering both constant coefficients and time-dependent schedules. The analysis was performed on sparse-view CT reconstruction with 4545 projection angles. For each choice of βt _t, we monitored the evolution of SSIM, LPIPS, PSNR, and the squared data residual along the reverse diffusion trajectory averaged over the training set. Following the geometric interpretation introduced in Section 3.3, we first considered the constant values βt=β,β∈0,0.5,1,1.5. _t=β, β∈\0,0.5,1,1.5\. Here, β=0β=0 corresponds to the non-fused DDNM approach, values in (0,1](0,1] reflects an interpolative regimes, whereas β=1.5β=1.5 produces an extrapolative update. We additionally evaluated three time-dependent schedules. Let st=1−tT∈[0,1] s_t=1- tT∈[0,1] (39) denote the normalized progress of the reverse diffusion process. Thus, st=0s_t=0 at the beginning of sampling, when t=Tt=T, and st=1s_t=1 at the end of the reverse trajectory, when t=0t=0. The three schedules are defined as Schedule A:βtA A: _t^A =1−st, =1- s_t, Schedule B:βtB B: _t^B =1.8(1−st), =1.8 (1- s_t ), Schedule C:βtC C: _t^C =1.8−0.6st. =1.8-0.6 s_t. All three schedules monotonically decrease the fusion coefficient as the reverse process approaches the clean image. However, they differ in the range of values attained by the fusion coefficient. Schedule A satisfies 0≤βtA≤1, 0≤ _t^A≤ 1, and therefore remains interpolative throughout the reverse process. Schedule B starts from the extrapolative value βTB=1.8 _T^B=1.8, crosses the threshold βt=1 _t=1, and approaches zero toward the end of the reverse trajectory. Schedule C instead satisfies 1.2≤βtC≤1.8, 1.2≤ _t^C≤ 1.8, and therefore preserves the extrapolative regime throughout the entire reverse trajectory while progressively reducing its strength. Figure 5: Evolution of SSIM, LPIPS, PSNR, and the squared data residual along the reverse diffusion trajectory for different constant coefficients and time-dependent schedules of βt _t. Results refer to sparse-view CT reconstruction with 4545 projection angles. Figure 5 compares the reconstruction dynamics obtained with the different choices of βt _t. During the initial part of the reverse trajectory, the metric curves are very similar. At these early stages of the reverse process, characterized by high diffusion noise, the consecutive denoised estimates still contain incomplete structural information, and their fusion has only a limited effect. The differences among the configurations become more evident during the final reverse steps, when the posterior estimates are more informative and their displacement defines a more meaningful reconstruction direction. The constant choices βt≤1 _t≤ 1 generally produce smoother trajectories than the extrapolative configurations. Since these choices do not introduce over-relaxation, their trajectories remain closer to those obtained with the standard DDNM update. A different behavior is observed for the constant coefficient βt=1.5 _t=1.5 and for Schedule C, both of which remain extrapolative throughout the reverse process. In this regime, the additional term t=(βt−1)(t+Δt−t) m_t=( _t-1) ( D_t+ t- D_t ) moves the reconstruction beyond the preceding posterior estimate along the direction identified by two consecutive reverse iterations. This over-relaxed update produces visible oscillations in SSIM, LPIPS, and PSNR during the final part of the trajectory. Nevertheless, the favorable phases of these oscillations reach higher SSIM and PSNR values and lower LPIPS values than the interpolative configurations. The comparison between the constant choice βt=1.5 _t=1.5 and Schedule C indicates that maintaining βt>1 _t>1 is not sufficient by itself. Although both choices produce an oscillatory behavior, the constant coefficient may end the reverse process in an unfavorable phase, resulting in a degradation of the final reconstruction. Schedule C progressively reduces the extrapolation strength from 1.81.8 to 1.21.2 and attains the best final image-quality metrics among the considered configurations. This result motivates the use of Schedule C in all the CT experiments reported in this work. The squared data residual exhibits an overall decreasing trend for all configurations, although the extrapolative choices introduce noticeable oscillations during the final reverse steps. For some extrapolative configurations, however, the last reverse iteration may occur during an unfavorable phase of the oscillatory dynamics, leading to worse image-quality metrics than those obtained at the preceding iteration. Interestingly, the validation experiments showed that the iteration with the smallest data residual did not necessarily correspond to the best image-quality metrics. This observation suggests that the residual trajectory may contain useful information for identifying favorable terminal iterates, although the design of a systematic stopping criterion is left for future work. Overall, the ablation study suggests that a properly modulated extrapolative schedule can exploit the directional information contained in consecutive posterior estimates more effectively than either a purely interpolative fusion or a constant over-relaxation coefficient. 5.3 Low-dose CT reconstruction In this section, we present the results obtained in the experiments simulating low-dose acquisitions as described in Section 4.3.1. Figure 6 evaluates the robustness of the reconstruction methods to different noise levels. The bar plots report the mean SSIM and LPIPS values obtained over the entire test set and across five runs with different initial noise realizations. The whiskers represent the corresponding standard deviations. As expected, reconstruction quality generally deteriorates as the noise intensity increases, as indicated by decreasing SSIM values and increasing LPIPS values. CPF-DDNM achieves the highest SSIM values across all the considered noise levels. Regarding LPIPS, DiffPIR provides the best performance at the lowest noise level, δ=0.1δ=0.1, whereas CPF-DDNM achieves the lowest values for δ=0.3δ=0.3 and δ=0.5δ=0.5. Interestingly, the relative advantage of CPF-DDNM becomes more pronounced as the measurement noise increases, particularly at the highest noise level. These results indicate that consecutive posterior fusion improves the robustness of the original DDNM formulation and remains effective under increasingly severe measurement noise. The qualitative results confirm that CPF-DDNM yields the best reconstruction, as further highlighted by the magnified view of a low-contrast region. In particular, noise is effectively suppressed while preserving sharp and well-defined object boundaries. FBP ℓ2 _2-TV DiffPIR DDNM CPF-DDNM Figure 6: Robustness to measurement noise in low-dose CT reconstruction. The upper part reports the mean SSIM and LPIPS values over the complete test set for noise levels δ=0.1δ=0.1, 0.30.3, and 0.50.5, with whiskers indicating the corresponding standard deviations. The lower part presents representative reconstructions and enlarged regions of interest for δ=0.3δ=0.3, comparing the preservation of anatomical structures and fine details across the considered methods. 5.4 Super-resolution of CT images Table 2: Quantitative comparison for ×4× 4 super-resolution in the noiseless setting and with noise level δ=0.01δ=0.01. Results are reported as mean ± standard deviation over the test images, after averaging each metric over five random initializations per image. Bold values indicate the best result for each setting. Noise Method SSIM ↑ LPIPS ↓ PSNR ↑ δ=0δ=0 DPS 0.497±0.0550.497± 0.055 0.288±0.0170.288± 0.017 19.60±1.2819.60± 1.28 PS+ 0.774±0.0490.774± 0.049 0.117±0.0150.117± 0.015 28.40±1.6628.40± 1.66 DiffPIR 0.863±0.0280.863± 0.028 0.093±0.0140.093± 0.014 30.18±1.3230.18± 1.32 DDNM 0.856±0.0310.856± 0.031 0.076±0.014 0.076± 0.014 31.40±1.60 31.40± 1.60 CPF-DDNM 0.872±0.028 0.872± 0.028 0.101±0.0140.101± 0.014 31.35±1.6031.35± 1.60 δ=0.01δ=0.01 DPS 0.495±0.0550.495± 0.055 0.288±0.0170.288± 0.017 19.60±1.2819.60± 1.28 PS+ 0.771±0.0550.771± 0.055 0.118±0.0140.118± 0.014 28.23±1.6028.23± 1.60 DiffPIR 0.839±0.028 0.839± 0.028 0.111±0.019 0.111± 0.019 29.59±1.27 29.59± 1.27 DDNM 0.808±0.0280.808± 0.028 0.131±0.0170.131± 0.017 29.03±1.3029.03± 1.30 CPF-DDNM 0.828±0.0270.828± 0.027 0.142±0.0140.142± 0.014 29.19±1.3929.19± 1.39 The final numerical experiment evaluates the considered methods for the ×4× 4 image super-resolution setting described in Section 4. Table 2 reports the mean SSIM, LPIPS, and PSNR values computed over the complete test set, while Figure 7 shows two representative reconstructions for the noisy case. In the noiseless setting, CPF-DDNM achieves the highest SSIM value and a PSNR essentially comparable to that of DDNM, indicating an effective preservation of structural information. A different behavior is observed when measurement noise is introduced. According to the test-set averages reported in Table 2, DiffPIR achieves the best values for all three metrics. CPF-DDNM nevertheless improves upon the original DDNM baseline in terms of SSIM and PSNR, although it obtains a higher LPIPS value. Thus, from a purely quantitative perspective, DiffPIR provides the best overall performance in the noisy setting. The qualitative comparison in Figure 7, however, reveals a marked discrepancy with this quantitative ranking. In both displayed examples, the DiffPIR reconstructions appear noticeably smoother, with some low-contrast tissues being less clearly resolved in the enlarged regions of interest. CPF-DDNM instead produces sharper tissue boundaries and preserves several localized details that are less visible in the DiffPIR reconstructions. Among the methods shown, CPF-DDNM therefore provides the most convincing visual results for these two noisy examples, despite the lower average metric scores reported in Table 2. This discrepancy illustrates that global image-quality metrics such as LPIPS may also be insensitive to localized differences that are clinically or anatomically relevant. DPS PS+ DiffPIR DDNM CPF-DDNM Figure 7: Qualitative comparison of the different methods for ×4× 4 super-resolution in the presence of measurement noise with δ=0.01δ=0.01. Two representative images from the test set are shown together with the corresponding reconstructions and enlarged regions of interest. 6 Conclusions In this work, we introduced CPF-DDNM, a lightweight inference-time strategy for improving the prior-driven recovery of unobservable image components in DDNM-based inverse solvers. The method exploits information already available along the reverse diffusion trajectory by combining consecutive measurement-aware estimates, without modifying or retraining the diffusion prior. The range/null-space decomposition of DDNM provides a direct geometric interpretation of the proposed fusion. With an exact Moore–Penrose pseudo-inverse, consecutive DDNM corrections share the same measurement-determined component and differ only in the null-space of the forward operator. CPF-DDNM therefore modifies exclusively the prior-driven component while preserving the information determined by the measurements. Thus, for noiseless data, CPF-DDNM preserves exact data consistency. Our local error analysis further provides conditions under which the extrapolative optimal fusion coefficient can yield a small null-space recovery error. Experiments demonstrate that the benefits of consecutive posterior fusion depend on both the forward operator and the fusion schedule. The largest gains were observed in sparse-view CT, where CPF-DDNM consistently improved upon DDNM across acquisition geometries and noise levels, and outperformed the considered diffusion-based inverse solvers. In super-resolution, CPF-DDNM improved structural similarity over DDNM, although DiffPIR achieved better average performance for some metrics in the noisy setting. The method was effective with both numerically approximated pseudo-inverses in CT and an analytically available pseudo-inverse in super-resolution. CPF-DDNM requires no additional denoiser evaluations or pseudo-inverse operations compared with DDNM, adding only the storage of the preceding correction and a linear combination. The main remaining limitation is the selection of the fusion schedule, which is currently determined on validation data and depends on the forward operator. Future work will investigate adaptive inference-time strategies based on quantities such as consecutive correction variations and data residuals, as well as plugging the CPF strategy to different posterior sampling methods. References [1] Yamin Arefeen, Brett Levac, Zach Stoebner, and Jonathan I. Tamir. INFusion: Diffusion regularized implicit neural representations for 2d and 3d accelerated mri reconstruction. In 2024 58th Asilomar Conference on Signals, Systems, and Computers, pages 1886–1890. IEEE, 2024. [2] Mario Bertero, Patrizia Boccacci, and Christine De Mol. Introduction to inverse problems in imaging. CRC press, 2021. [3] M Jorge Cardoso, Wenqi Li, Richard Brown, Nic Ma, Eric Kerfoot, Yiheng Wang, Benjamin Murrey, Andriy Myronenko, Can Zhao, Dong Yang, et al. Monai: An open-source framework for deep learning in healthcare. arXiv preprint arXiv:2211.02701, 2022. [4] 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(1):120–145, 2011. [5] Hyungjin Chung, Jeongsol Kim, Michael Thompson Mccann, Marc Louis Klasky, and Jong Chul Ye. Diffusion Posterior Sampling for General Noisy Inverse Problems. In The Eleventh International Conference on Learning Representations, pages 1–30, 2023. [6] Heinz Werner Engl, Martin Hanke, and Andreas Neubauer. Regularization of inverse problems, volume 375. Springer Science & Business Media, 1996. [7] Davide Evangelista. Physics-guided flow matching for ct image reconstruction, 2026. [8] Davide Evangelista, Elena Morotti, Francesco Pivi, and Maurizio Gabbrielli. Improving diffusion posterior samplers with lagged temporal corrections for image restoration. arXiv preprint arXiv:2605.12573, 2026. - accepted at NeurIPS 2026. [9] Peter F. C. Gilbert. Iterative methods for the three-dimensional reconstruction of an object from projections. Journal of Theoretical Biology, 36(1):105–117, 1972. [10] Harshit Gupta, Kyong Hwan Jin, Ha Q Nguyen, Michael T McCann, and Michael Unser. Cnn-based projected gradient descent for consistent ct image reconstruction. IEEE transactions on medical imaging, 37(6):1440–1453, 2018. [11] Per Christian Hansen. Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion. SIAM, 1998. [12] Per Christian Hansen, James G Nagy, and Dianne P O’leary. Deblurring Images: Matrices, Spectra, and Filtering. SIAM, 2006. [13] 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. [14] Jakob S Jorgensen, Emil Y Sidky, and Xiaochuan Pan. Quantifying admissible undersampling for sparsity-exploiting iterative image reconstruction in x-ray ct. IEEE Transactions on Medical Imaging, 32(2):460–473, 2012. [15] Avinash C Kak and Malcolm Slaney. Principles of computerized tomographic imaging. SIAM, 2001. [16] Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song. Denoising diffusion restoration models. Advances in Neural Information Processing Systems, 35:23593–23606, 2022. [17] Amirhossein Kazerouni, Ehsan Khodapanah Aghdam, Moein Heidari, Reza Azad, Mohsen Fayyaz, Ilker Hacihaliloglu, and Dorit Merhof. Diffusion models in medical imaging: A comprehensive survey. Medical Image Analysis, 88:102846, 2023. [18] Luping Liu, Yi Ren, Zhijie Lin, and Zhou Zhao. Pseudo numerical methods for diffusion models on manifolds. arXiv preprint arXiv:2202.09778, 2022. [19] Elena Loli Piccolomini and Elena Morotti. A model-based optimization framework for iterative digital breast tomosynthesis image reconstruction. Journal of imaging, 7(2):36, 2021. [20] Elena Loli Piccolomini, Marco Prato, Margherita Scipione, and Andrea Sebastiani. Ctprintnet: an accurate and stable deep unfolding approach for few-view ct reconstruction. Algorithms, 16(6):270, 2023. [21] Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. DPM-Solver++: Fast solver for guided sampling of diffusion probabilistic models. Machine Intelligence Research, 22(4):730–751, 2025. [22] C McCollough. Tu-fg-207a-04: Overview of the low dose ct grand challenge. Medical physics, 43(6Part35):3759–3760, 2016. [23] Elena Morotti, Davide Evangelista, and Elena Loli Piccolomini. A green prospective for learned post-processing in sparse-view tomographic reconstruction. Journal of Imaging, 7(8):139, 2021. [24] Elena Morotti, Davide Evangelista, Andrea Sebastiani, and Elena Loli Piccolomini. Adaptive weighted total variation boosted by learning techniques in few-view tomographic imaging. Journal of Scientific Computing, 106(3):74, 2026. [25] Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021, pages 1–20, 2021. [26] 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. [27] Yinhuai Wang, Jiwen Yu, and Jian Zhang. Zero-shot image restoration using denoising diffusion null-space model. In The Eleventh International Conference on Learning Representations, pages 1–31, 2023. [28] Zhou Wang, Eero P Simoncelli, and Alan C Bovik. Multiscale structural similarity for image quality assessment. In The Thirty-Seventh Asilomar Conference on Signals, Systems & Computers, 2003, volume 2, pages 1398–1402. Ieee, 2003. [29] George Webber and Andrew J. Reader. Diffusion models for medical image reconstruction. BJR—Artificial Intelligence, 1(1):ubae013, 2024. [30] Wenjun Xia, Wenxiang Cong, and Ge Wang. Patch-based denoising diffusion probabilistic model for sparse-view ct reconstruction. arXiv preprint arXiv:2211.10388, 2022. [31] Taofeng Xie, Zhuo-Xu Cui, Chen Luo, Huayu Wang, Congcong Liu, Yuanzhi Zhang, Xuemei Wang, Yanjie Zhu, Guoqing Chen, Dong Liang, Qiyu Jin, Yihang Zhou, and Haifeng Wang. Joint diffusion: Mutual consistency-driven diffusion model for PET-MRI co-reconstruction. Physics in Medicine & Biology, 69(15):155019, 2024. [32] Liutao Yang, Jiahao Huang, Guang Yang, and Daoqiang Zhang. CT-SDM: A sampling diffusion model for sparse-view ct reconstruction across various sampling rates. IEEE Transactions on Medical Imaging, 44(6):2581–2593, 2025. [33] Qinsheng Zhang and Yongxin Chen. Fast sampling of diffusion models with exponential integrator. In International Conference on Learning Representations, pages 1–14, 2023. [34] Richard Zhang, Phillip Isola, Alexei A Efros, Eli Shechtman, and Oliver Wang. The Unreasonable Effectiveness of Deep Features as a Perceptual Metric. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 586–595, 2018. [35] Yuanzhi Zhu, Kai Zhang, Jingyun Liang, Jiezhang Cao, Bihan Wen, Radu Timofte, and Luc Van Gool. Denoising Diffusion Models for Plug-and-Play Image Restoration. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1219–1229, 2023.