Analysis of iterative deblurring: no explicit noise

Sinethemba Neliswa Mamba; Paweł Danielewicz

doi:10.1515/phys-2025-0252

Article Open Access

Analysis of iterative deblurring: no explicit noise

Sinethemba Neliswa Mamba and Paweł Danielewicz

Published/Copyright: December 24, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 23 Issue 1

Abstract

Iterative deblurring, notably the Richardson–Lucy algorithm with and without regularization, is analyzed in the context of nuclear and high-energy physics applications. In these applications, probability distributions may be discretized into a few bins, measurement statistics can be high, and instrument performance can be well understood. In such circumstances, it is essential to understand the deblurring first without any explicit noise considerations. We employ singular value decomposition for the blurring matrix in a low-count pixel system. A strong blurring may yield a null space for the blurring matrix. Yet, a nonnegativity constraint for images built into the deblurring may help restore null-space content in a high-contrast image with zero or low intensity for a sufficient number of pixels. For low-contrast images, control over null-space content can be achieved through regularization. When regularization is applied, the blurred image is, in practice, restored to one that is still blurred but less than the starting image.

Keywords: deblurring; deconvolution; Richardson–Lucy algorithm; heavy-ion collisions; collective flow

1 Introduction

Blurring is common in image acquisition, making the images less sharp and clear. When the blurring process is understood, it may be possible to remove, fully or partially, the blurring effects or to deblur the images. Here, we will analyze what a standard iterative deblurring methodology can and cannot do, with an eye toward improving the quantitative results of particle yield measurements in nuclear and high-energy physics. In the latter areas, the deblurring gained attention relatively recently [1], [2], [3], [4], [5]. The measurement statistics can be higher than in many optical applications, making non-nominal noise analysis relevant. At the same time, the number of bins in these areas, equivalent to pixels, can be modest.

A convenient framework in which the deblurring can be understood is that of the Singular Value Decomposition (SVD) for the blurring matrix, already used in the past in the literature [6], [7], [8]. With some advantages that the Richardson–Lucy (RL) [9], 10] and related iterative methods [2], 11], 12] have over other iterative methods for the restoration of brightness [13], [14], [15], we will concentrate here on the RL method. One novel aspect of our work will examine the contrast’s impact on the objects imaged in the deblurring process. Further, we will actively seek cases where deblurring fails, understand the failures, and predict when they may occur. Small or vanishing eigenvalues from SVD will play a role there. In the present work, we will focus on deblurring when no explicit noise is present in the problem, which is already complex. In future work, we shall consider practicalities around the noise in the problem, for which the present work will provide a base.

This paper is organized as follows. Section 2 discusses general blurring concepts and formulates a schematic model for explorations. We decompose the blurring matrix using SVD. Section 3 describes selected deblurring methods, including the RL and Landweber (LW) algorithms, and illustrates their practical operation. Section 4 explores the limits of the RL deblurring with and without regularization. Section 5 discusses how concurrent nuclear physics efforts connect to this paper’s explorations. We conclude in Section 6.

2 Blurring matrices and their singular value decomposition

The blurring relation can be stated in the form of the equation [9]

(1) g = T F .

Here, g ∈ R n is the vector of values for the blurred image, T ∈ R n × n is the blurring matrix or transfer matrix, and F ∈ R n is the vector for the blur-free image. Probability conservation for the blurring matrix implies that the sum of all the elements of each column (index i) is precisely equal to 1:

∑ j = 0 n − 1 T j i = 1 .

Without practical loss of generality, we shall consider that the indices correspond to locations in one dimension, representing pixels or other distribution bins. We shall adopt cyclic boundary conditions within that space, where the vector and matrix indices are computed modulo n. To fix the attention, we shall consider a couple of exemplary blurring (or point-spreading) functions, three- and five-bin, out of which we shall consider the five-bin one most often:

Three-bin

(2) T j i = 0.550 , for j = i , 0.225 , for j = ( i ± 1 ) mod n , 0 , otherwise.

Five-bin

(3) T j i = 0.250 , for j = i , 0.250 , for j = ( i ± 1 ) mod n , 0.125 , for j = ( i ± 2 ) mod n , 0 , otherwise.

In the grayscale convention, the pixel brightness values span the interval [0, 1], with 1 representing white and 0 representing black. In deblurring, it is important not to end up with values outside of that interval, as these would lack physical sense. We illustrate the three- and five-bin blurring functions in Figures 1(a) and 2(a), respectively, where we show the action of these functions on an original image where only pixel 5 has a value of 1, and the others have 0. The overall size of the space is n = 12.

Figure 1:

Illustration of the three-bin blurring function. (a) Image from blurring an original with pixel 5 at 1 and other at 0, for n = 12. (b) Singular values of the blurring function, ordered by their magnitude. Lines join the values to guide the eye.

Figure 2:

Illustration of the five-bin blurring function. (a) Image from blurring an original with pixel 5 at 1 and other at 0, for n = 12. (b) Singular values of the blurring function, ordered by their magnitude. Lines join the values to guide the eye.

Upon applying the transposed blurring matrix to both sides of Eq. (1), we find the norm equation

(4) T ⊤ g = T ⊤ T F ≡ M F .

If we apply SVD to T, we get

(5) T = ∑ α = 0 n − 1 σ α u α v α ⊤ ,

where u and v are, respectively, the separately orthonormal left and right singular vectors, and we adopt the convention of the singular values σ to be nonnegative and ordered in nonincreasing order. Then for the positive definite matrix M we arrive at the following spectral decomposition

(6) M = T ⊤ T = ∑ α , β σ α σ β v α u α ⊤ u β v β ⊤ = ∑ α σ α 2 v α v α ⊤ .

Generally, for original and deblurred images, we shall employ a decomposition in terms of the right singular vectors of the blurring matrix and for the blurred images – in terms of the left:

(7) F = ∑ α F α v α , g = ∑ α g α u α .

With (1), we then find

(8) g = T F = ∑ α σ α F α u α ,

or g _α = σ _α F _α.

The three- and five-bin blurring matrices, Eqs. (2) and (3), which we use here as examples, have some peculiar features. These features are often encountered in practice, possibly as approximate rather than exact features. Thus, these blurring matrices are symmetric, meaning their singular vectors are matrix eigenvectors. Further, with the matrix elements of T depending only on the absolute difference of indices, that dependence extends to the corresponding M matrices. This implies two types of invariance for the matrices T and M, which affect the singular vectors. One is the invariance under a shift of matrix indices by some integer. With this, a singular vector, corresponding to some singular value, remains a singular vector, corresponding to the same singular value, after a shift in its indices by an integer. The second invariance is under a reflection of the indices around any chosen index. With this, a singular vector, corresponding to some singular value, remains a singular vector, corresponding to the same singular value, after a reflection of its indices around any chosen index. With these two symmetries, any singular value is at least doubly degenerate unless the corresponding singular vector transforms onto itself, up to a factor of −1, under these transformations.

The spectra of singular values for the two exemplary blurring matrices are shown in Figures 1(b) and 2(b), and the singular vectors for the five-bin function are shown in Figure 3. Double degeneracy is seen for most of the singular values. For Figure 3, we choose the corresponding pairs of singular vectors in the degenerate subspaces so that some symmetry behind the degeneracies is evident, either a reflection or translational one, or both. In each of the spectra, in Figures 1(b) and 2(b), a nondegenerate singular value of 1 is seen. This value corresponds to the vector that is constant across the pixels. Under blurring, this vector transforms into itself. Also, this vector is invariant under the symmetries. The spectrum of singular values falls off faster for the five- than for the three-bin function, demonstrating a more significant deterioration of the information under the blurring with a broader function. In Figure 3, it is evident that a faster variation of the vectors over pixels is generally associated with lower singular values. The singular vector with the most rapid variation, alternating between the same positive and negative values from one pixel to another, is invariant under the symmetries. It is the other one that may be nondegenerate. Also, under blurring, it transforms into itself, up to a factor, because it is uniquely alternating and reflection-symmetric about every pixel, which will be preserved by any blurring that depends only on the difference of indices. With a broadening of the blurring function, this vector will likely be the first to be nullified. In Figure 1(b), the singular value for that vector, with index 11, approaches zero, and the vector is nullified when the blurring matrix evolves into T _jj = 0.5 and T _ji = 0.25, for j = (i ± 1) mod 12. In Figure 2(b), the maximally alternating vector is accompanied in the null space by two more vectors, cf., panels (a) and (f) in Figure 3.

Figure 3:

Singular vectors for the five-bin blurring function. Within degenerate spaces, pairs of vectors are chosen to transform onto the other or itself under reflection and/or translation. Panel (a) shows the two unique vectors that transform onto themselves under the transformations.

3 Deblurring methods

We next turn to deblurring. We illustrate the deblurring in Figure 4(a). We take a δ distribution where only one bin gets intensity 1, and the others get 0. We blur the original image with the 5-bin function (3) and deblur the result using the RL method [9], 10], details of which will be discussed. The blurred and restored images are shown in Figure 4(a) together with the original. The eye cannot distinguish between the RL-restored and original images, so they are represented with a single set of graphics. In addition, the intensity values for the images are shown in Table 1, and within the precision of the machine we employ, they cannot be distinguished. This may be perceived as astounding, as the singular vectors with the fastest bin-to-bin variation span the null space for the blurring matrix, and they most certainly contribute to the SVD of the δ distribution.

Figure 4:

Deblurring test with a single pixel with index 5 at the intensity 1 and all other pixels at the intensity of 0. (a) The original image F (circles), image g (triangles) from blurring F with the 5-pixel blurring function, and the restored images f from different deblurring methods are displayed. As the image restored by the RL method is indistinguishable from the original, the circles represent that image, too. The image restored using the LW method is shown as diamonds. Finally, the image restored with SVD, where the operational space was reduced (m = 6 in Eqs. (21) and (22)) compared to LW, is shown as circles. The symbols corresponding to each case are connected with lines. (b) Coefficients of decomposition in the basis of singular vectors versus singular value index for the original and blurred images and images restored with the RL and LW methods.

Table 1:

Bin intensities in deblurring the δ-distribution centered at index 5, illustrated in Figure 4(a). From left to right, the columns represent the bin index (i), and then the intensity for the original image (F _i), for the 5-bin blurred image (g _i), and the image restored with three different methods: RL f i RL , LW f i LW and SVD f i SVD .

i	F _i	g _i	f i RL	f i LW	f i SVD
0	0	0.000	0.000	0.083	−0.061
1	0	0.000	0.000	−0.250	−0.167
2	0	0.000	0.000	0.083	0.083
3	0	0.125	0.000	0.083	0.000
4	0	0.250	0.000	0.083	0.228
5	1	0.250	1.000	0.750	0.583
6	0	0.250	0.000	0.083	0.228
7	0	0.125	0.000	0.083	0.000
8	0	0.000	0.000	0.083	0.083
9	0	0.000	0.000	−0.250	−0.167
10	0	0.000	0.000	0.083	−0.061
11	0	0.000	0.000	0.083	0.250

In the following, we review exemplary deblurring methods and their potential utility for quantitative research analyses. After this section, we shall examine under what circumstances the deblurring can yield faithful results, such as in the example above, and we shall seek conditions and a manner under which the deblurring can fail. Iterative methods aim to restore a blurred image to its original form by iteratively improving estimates of the blur-free image. These methods take as input an initial guess of the blur-free image, the blurring matrix, and a criterion for the number of iterations. Typically, the initial guesses are the uniform or blurred image. Since the blur is to be suppressed at each iteration, naively, a more accurate output is expected with more iterations – in practice, this may not be the case.

While discussing the iterative methods below, we use k as an iteration index.

3.1 Richardson–Lucy method, classical and regularized

The classical RL algorithm iteratively upgrades the restored image f, given the blurred image g, according to [2], 9], 10]

(9) f i ( k + 1 ) = A i ( k ) f i ( k ) .

Here, the amplification factor A _i for the i’th pixel intensity is

(10) A i ( k ) = 1 q i ∑ j T j i g j g j ( k ) ,

i, j, m = 0, 1, 2, …, n − 1, g j ( k ) = ∑ m T j m f m ( k ) is the blurred-image prediction for the restored image f ^(k) at the k’th step of restoration, T is the before introduced point spread function or blurring matrix, and q _i = ∑_j T _ji serves the normalization. The normalization factors for the introduced blurring-function examples are unity, q _i ≡ 1.

The regularized RL algorithm modifies the iteration equation to

(11) f i ( k + 1 ) = A i ( k ) I i k f i ( k ) i , j , m = 0,1,2 , … , n − 1 ,

where the regularization factor is [2], 11]

(12) I i ( k ) = 1 1 − Λ ÷ ∇ f i ( k ) | ∇ f i ( k ) | = 1 1 − Λ , if f i ( k ) < f [ ( i − 1 ) mod n ] , [ ( i + 1 ) mod n ] ( k ) , 1 1 + Λ , if f i ( k ) > f [ ( i − 1 ) mod n ] , [ ( i + 1 ) mod n ] ( k ) , 1 , otherwise .

for i = 0, 1, 2, …, n − 1. Here, the r.h.s. is the one-dimensional version [2] of the center expression [11] in which ∇ is the gradient and | ∇ f i ( k ) | is ℓ ₁-norm; Λ is a small positive regularization coefficient. The role of the regularization factor is to suppress an uncontrolled rise, during iterations, in the contributions to the restored image from the fastest-varying singular vectors corresponding to vanishing or minimal singular values. The typical signature of such an uncontrolled rise is developing a bin-to-bin seesaw pattern in the restored image. As SVD might not be applied in parallel to the RL-restoration and the direct connection of the regularization factor I to SVD is not obvious, it is common to make Λ small enough that no pathological signs in the restored image appear.

The restoration of the δ-distribution of Figure 4(a) with the RL method, from the blurred image there, with or without regularization, yields a result that is not distinguishable by eye from the original, so they share the graphical representation in the figure. The numerical values for the restoration are provided in Table 1. Another test to which we subject the deblurring methods is presented in Figure 5. The original image F is a ramp with pixels 5–8 at intensity 1 and other pixels at intensity 0. The blurring, yielding g, is performed with the 5-bin blurring function. The results of RL deblurring, without or with regularization, are not distinguishable by eye from the original, but so are the results of deblurring with other methods that we will discuss next. The numerical values from the ramp tests are shown in Table 2.

Figure 5:

Deblurring test with a 4-pixel ramp. In the original image, the intensity is 1 for the pixels 5–8 and 0 for the others. The operational space for SVD is here the full row space of T.

Table 2:

The numerical values from the 4-pixel ramp deblurring test, with the RL, SVD, and LW methods, illustrated in Figure 5.

i	F _i	g _i	f i RL	f i LW	f i SVD
0	0	0.000	0.000	0.000	0.000
1	0	0.000	0.000	0.000	0.000
2	0	0.000	0.000	0.000	0.000
3	0	0.125	0.000	0.000	0.000
4	0	0.375	0.000	0.000	0.000
5	1	0.625	1.000	1.000	1.000
6	1	0.875	1.000	1.000	1.000
7	1	0.875	1.000	1.000	1.000
8	1	0.625	1.000	1.000	1.000
9	0	0.375	0.000	0.000	0.000
10	0	0.125	0.000	0.000	0.000
11	0	0.000	0.000	0.000	0.000

3.2 Landweber method

The Landweber (LW) method [14] can be viewed as the steepest-descent method for minimizing the square deviation between the brightness distribution for the generated and blurred images,

(13) g − T f 2 2 .

Iterative adjustments in f take the form

(14) f ( k + 1 ) = f ( k ) + β T T ( g − T f ( k ) ) ,

where 0 < β < 2/‖T ^T T‖ for convergence. With the T multiplication in (14), it is apparent that the iterations do not change the null-space content in f. The null-space content may be suppressed completely in f by replacing the minimized quantity (13) with

(15) g − T f 2 2 + ϵ f 2 2 ,

where 0 < ϵ ≪ ( σ α ) min 2 and ( σ α ) min is the smallest singular value of which the contribution to f should not be disturbed. The iterations are then modified to

(16) f ( k + 1 ) = ( 1 − β ϵ ) f ( k ) + β T T ( g − T f ( k ) ) .

The results from the deblurring of the blurred δ-distribution with the LW method are provided in Figure 4 and Table 1. The corresponding results from processing the 4-pixel ramp are provided in Figure 5 and Table 2. While the LW method struggles to restore the δ distribution, yielding, in particular, negative intensity values, it faithfully restores the 4-pixel ramp. This may be surprising, and we will return to this issue later in this Section.

3.3 Deblurring using singular value decomposition

The deblurring can be further carried out by employing SVD explicitly. In general, the SVD deblurring [6], [7], [8] seeks the best approximation to the blurred image in a reduced subspace of singular vectors, which in particular should exclude the null space for T, minimizing the square deviation

(17) g − T m f 2 2 ,

where

(18) T m = ∑ α = 0 m − 1 σ α u α v α ⊤ ,

and m ≤ n. The minimization of (17) is equivalent to the minimization of the square deviation within the reduced subspace only

(19) g m − T m f 2 2 ,

where

(20) g m = ∑ α = 0 m − 1 g α u α ,

with the minimum reached for

(21) f α = g α σ α , for α < m .

Outside of the considered subspace, the expansion coefficients for f are set to zero,

(22) f α = 0 , for m ≤ α < n .

In the tests of this Section, the SVD method yields the same results as LW when its operational space is maximal, i.e., identical to the row space of the blurring matrix T. Regarding the 4-pixel ramp test with results in Figure 5 and Table 2, we keep the operational space for SVD maximal. However, in the δ-function test shown in Figure 4 and Table 1, we reduce the operational space for SVD to m = 6. See Figure 3 for the corresponding singular vectors. With the reduction, the restoration’s quality significantly deteriorates. The location and strength of negative intensities in the restored image are similar for m = 8 (identical to LW within precision) and m = 6, cf. Eqs. (21) and (22).

3.4 Deblurring using deep learning

The above traditional deblurring methods rely on predefined blurring, which can be a limitation when the nature of blurring for the images is not fully understood. Deep Learning approaches can leverage large datasets and Deep Neural Networks (DNNs) to infer the mapping between blurred and original images [16]. Due to its data-driven nature, training the DNNs requires high-quality datasets containing pairs of blurred and blur-free images, often synthesized. The idea is to train DNNs to recover the original image F from blurred inputs g without explicitly modeling or knowing the blurring function T in Eq. (1).

3.5 Methods’ assessment

We call the original images, such as in Figures 4(a) and 5, where a substantial fraction of the pixels has an intensity zero or low compared to the maximal, as having high contrast. When such images are moderately blurred, the deblurring methods with the nonnegativity constraint for intensity built into the deblurring process, such as RL, generally perform far better than methods that lack the constraint, such as LW or SVD, see the case of the δ-distribution of Figure 4(a). In this case, the two latter methods yield negative intensity values in the restored images. Amazingly, while the δ-distribution contains contributions from the v ₁₀ and v ₁₁ vectors in the null space, the nonnegativity constraint allows the RL method to restore their contributions faithfully despite their singular values being zero, see Figure 4(b). Without the nonnegativity constraint, the content of the null space, in the case of LW, or excluded space, in the case of SVD, stays at zero.

One puzzling outcome of the tests so far is that the LW method and the SVD method, when its excluded space is identical to the null space, perform as well as RL in restoring the 4-pixel ramp in Figure 5, even though they struggle with the δ-distribution. The mystery is solved when looking at the null-space singular vectors in Figure 3 and the original ramp in Figure 5: a ramp with an even number of pixels at the same intensity lacks any null-space contribution for our blurring function of Eq. (3). When the ramp’s extension is changed to an odd number of pixels or the ramp is sloped, the LW and SVD methods can develop negative intensity values for the restored image, just like in the case of the δ distribution.

Given the benefits of the built-in nonnegativity constraint, we will further exclusively rely on the RL deblurring method. Other methods may perform similarly when supplemented with this type of constraint. In employing the RL method, the iteration progress is sometimes used as a tool. Figure 6 shows how the expansion coefficients in singular vectors for the restored image behave as a function of the RL iteration step. Here, the original image is the δ-distribution, and the starting image is the uniform distribution. One can observe that the higher the index of the singular vector, the generally slower the approach of the coefficient to its asymptotic value. This can be understood in terms of the weakening impact of expansion coefficients on the amplification factors (10) when the singular value decreases, as

(23) d A i ( k ) d f α ( k ) ∝ σ α ,

for g _i > 0. Null-space and low-singular-value coefficients will be affected by the positivity constraint and/or regularization, and in Figure 6 we observe a lumping in the pace of their approach to the asymptotic values.

Figure 6:

Restored-image expansion coefficients in singular vectors, as a function of the RL iteration step k, normalized to the coefficients for the original δ-distribution image, for even-index (a) and odd-index (b) singular vectors.

4 Exploring the limits of deblurring

In the context of deblurring applications in quantitative research in nuclear and particle physics, it is important to understand when and how deblurring can falsify results. We have observed the role that positive definiteness could play in restoring the null-space contributions to a distribution. However, intuitively, one might expect the positive definiteness to play a minor role when the original distribution is uniformly far from zero on the distribution’s general scale of variation. Below, we will compare the performance of the classic RL method in those different situations. When there is no positive definiteness to fall back on, regularization can improve the performance of a deblurring method. However, the regularization must be a compromise with some adverse effects, too, as may be evident from Eq. (15).

4.1 High versus low contrast

For high-contrast original images, such as in Figures 4(a) and 5, we have observed that the deblurring with the nonnegativity constraint can restore the null-space components from the blurred images. Notably, for the nonnegativity constraint to be impactful, the number of pixels at low intensity compared to the rest must be comparable to the dimension of the null space. We now confront that situation to the situation with low-contrast original images, where the intensity for most pixels varies within a small range relative to the maximal intensity.

An exemplary low-contrast case of a δ-function superimposed on a high uniform background is illustrated in Figure 7(a). When the classical RL method is applied to the blurred image, the method fails, on the scale of the δ-function norm, in a similar manner as the methods without the nonnegativity constraint in the high-contrast case of Figure 4(a). Again, we start the RL iterations in this illustration with a uniform image. Although the nonnegativity is present in the RL method, it does not affect restoration when the background is sufficiently elevated. Figure 7(b) complements the results of (a) with coefficients of decomposition in the basis of singular vectors, demonstrating that the coefficients for the null space are not restored for the classical RL method, unlike in Figure 4(b). For impactful constraints, we turn to the regularized RL method, in which restoration is modified relative to the classical method based on differences in the intensity of adjacent pixels, thereby suppressing local extrema. With a low Λ = 0.0001, the original image is satisfactorily restored in Figure 7(a). The Root Mean Squared Error (RMSE) values for the restorations of the image in Figure 7 and others, using different values of Λ, are provided in Table 3.

Figure 7:

Deblurring test with a single pixel with index 5 at intensity elevated by 0.25 over a high uniform background at 0.75. (a) Circles represent the original image F and triangles – the image g from blurring F with the 5-pixel blurring function. Moreover, the restored images f from different versions of the RL method are displayed. Wide diamonds represent the image restored using the classical RL method. Stars, slim diamonds, and rectangles represent the images restored with the regularized RL method at Λ = 0.0001, 0.001, and 0.01, respectively. The symbols corresponding to each case are connected with lines. (b) Coefficients of decomposition in the basis of singular vectors for the cases in (a) versus singular-value index.

Table 3:

RMSE values for the images of Figures 7 and 8, blurred (g) and restored (f _RL) at different regularization strengths Λ.

Figure	Case	RMSE
Figure 7	f _RL (Λ = 0)	0.032558
	f _RL (Λ = 0.0001)	0.002206
	f _RL (Λ = 0.001)	0.022123
	f _RL (Λ = 0.01)	0.059650
	g	0.061184
Figure 8(a)	f _RL (Λ = 0)	0.033004
	f _RL (Λ = 0.0001)	0.000377
	f _RL (Λ = 0.001)	0.003743
	f _RL (Λ = 0.01)	0.037838
	g	0.058463
Figure 8(c)	f _RL (Λ = 0)	0.035454
	f _RL (Λ = 0.0001)	0.000196
	f _RL (Λ = 0.001)	0.002004
	f _RL (Λ = 0.01)	0.019756
	g	0.057054
Figure 8(e)	f _RL (Λ = 0)	0.066372
	f _RL (Λ = 0.0001)	0.000171
	f _RL (Λ = 0.001)	0.001903
	f _RL (Λ = 0.01)	0.017943
	g	0.114109

The smallest RMSE value for each set of strengths is printed in bold.

Additional examples of image restoration at different contrast levels are shown in Figure 8.

Figure 8:

Deblurring tests such as in Figure 7 for other exemplary low-contrast images, with matching information for each of these presented in panels (a) and (b), (c) and (d), and (e) and (f), respectively.

4.2 Impact of Initialization

For a low-contrast original image F, it is possible to add a significant null-space component to that image, without violating intensity nonnegativity. When such a combination is fed to the classical RL method as the starting guess f ⁽⁰⁾, mathematically, the restoration iterations will never suppress the added spurious null-space component, returning that starting guess as the fully restored image f. Even components added from low-singular-value vectors may only be slowly modified over iterations.

To avoid excessive spurious null-space or low-singular-value admixtures in the restored image, it helps to start the classical RL iterations with an image that lacks such components, such as a uniform image. In Figure 9, we show what happens within the classical RL restoration when the restoration of a low-contrast image is started with a significant null-space component in the first guess for the image instead of the uniform guess. The null-space component’s strength persists at a significant level for Λ = 0 in the limit of many restoration steps. Notably, the null-space components are unlikely to converge to the values of the original image in the absence of an impact from the nonnegativity constraint, regardless of how the iterations are started, including with the uniform image. In Figures 7 and 8, we see that classical restorations of low-contrast images that start with uniform guesses tend to underestimate the magnitudes of the null-space components in the restored images. In particular, the panels (a) and (b) in Figure 8 are analogous to Figure 9 but differ importantly in initialization. As shown in Figures 7–9, the null-space components are brought under some control after a modest amount of regularization is applied.

Figure 9:

Sensitivity of an RL restoration to null-space content at iteration start and to the regularization strength. The RL iterations here start with the original image, with an α = 11 component added. Panels (a) and (b) show the images and singular-vector coefficients, respectively. The classical RL method (Λ = 0) leaves the null-space content intact. An RL method with regularization (Λ ≳ 10⁻⁴) suppresses the superimposed spurious null-space content, but excessive regularization, such as for Λ = 0.01, can chip away the null-space content in the restored image compared to the original. Panel (c) shows the RMSE values for the restored images as a function of Λ.

4.3 Impact of regularization

In examining the cases of restoration of interesting original images with few significant intensity jumps in Figures 7 and 8, we can see that the regularization can dramatically improve image restoration. The particular regularization of Eqs. (11) and (12) biases against multiple pronounced extrema and, with that, suppresses significant spurious null-space contributions that could be added to the restored low-contrast image for the classical RL method.

While a modest regularization within the RL method can help restore the null-space content of the original image, an excessive regularization can act to suppress not only that content but also contributions from singular vectors corresponding to low nonzero singular values. This is illustrated in Figures 7 and 8, which include results from RL restorations with progressively strengthening regularization. As the regularization strength increases, the features of the restored images evolve from those of the original image toward those of the blurred image. Figure 9(c) illustrates the dependence of RMSE for the restored image compared to the original image on the parameter Λ for the regularization strength. Excessively low values of Λ may not be sufficient to tame the null-space components intentionally superimposed on the image starting guess. The growth of RMSE with the growth of Λ for the higher Λ values can be seen, reflecting the gradual blurring in the features of the restored image. In effect, there is an optimal range of Λ in a restoration problem for low-contrast images. The choice of regularization amounts to trading a stronger blurring in the image before deblurring for a weaker blurring afterward.

4.4 More pixels

The case of just 12 pixels and a 5-pixel blurring function allows for an exploration of details in blurring and deblurring. For more pixels and wider blurring functions, the features need to be looked at in a more statistical manner. The separation between the null space and its complement may become less clear-cut, as the singular values may fall by many orders of magnitude with increasing index before hitting zero. Even with the need for some perspective shift, our experimentation reveals similarities in the deblurring across both higher and lower pixel counts. As an example, we take the case of 48 pixels and start with the 5-bin blurring function of Eq. (3). We construct a wider 13-bin blurring function, commensurate with the larger pixel space, by taking a third power of the 5-bin matrix, which is equivalent to a trifold blurring with the 5-bin function. With this, the singular vectors remain the same as for the 5-bin function, but now in the 48-pixel space, whereas the singular values are raised to the third power, causing the spectrum to fall off by orders of magnitude with the index before reaching the null space. In Figure 10, we show the case of blurring and deblurring of a low-contrast ramp similar to the one for 12 pixels in Figure 8(a). Again, the restoration with no regularization yields a much distorted image. A relatively wide range of regularization strengths yields satisfactory restored images, and only at λ ∼ 0.01 do they become excessively biased.

Figure 10:

Deblurring tests for 48 pixels and the original low-contrast image similar to those represented in panels (a) and (b) of Figure 8.

5 From experimentation to data analysis

We briefly discuss how this current backbone experimentation connects to the parallel nuclear data analysis. The blurring due to reaction-plane resolution in nuclear collisions [2], i.e., misassigning the azimuthal angle relative to the reaction plane to an emitted particle, nominally occurs in three dimensions. However, this blurring and deblurring occur only in one dimension of the azimuthal angle around the collision beam axis. The deblurring is repeated for separate velocity bins along and perpendicular to the beam. In other recent nuclear applications in nuclear physics, such as decay-energy and γ-ray spectra [3], 5] and the imaging of spatial sources from particle correlations [17], 18], deblurring is only one-dimensional. However, the applications may progress to more dimensions when, e.g., three-dimensional features of the sources are sought.

In the reaction-plane deblurring [2], the particle yield distributions are periodic, like here, as the azimuthal angle is used in the process. The transfer matrix for decay-energy spectrum deblurring is constructed by simulating the detection system’s performance and is independent of the distribution it operates on, as here. In the source-imaging from correlations, the transfer matrix is constructed theoretically and may be folded with detector resolution. It is independent of the source that is inferred. The case of reaction-plane deblurring is an interesting situation in which the intrinsic consistency of the deblurring procedures can be tested as they are applied in regions where detector efficiency nears 100 % and falls towards zero.

Elaborating on the reaction-plane case, the detectors perform well at detecting charged particles moving at high velocities, such as those close to the target velocity, but perform poorly at low velocities, approaching the target at rest in velocity space. When only the forward-velocity region is analyzed, exploiting the central limit theorem, the transfer matrix can be well estimated from a few moments of the detected distribution [19], making it independent of the deblurring outcome. With physically expected reflection symmetry relative to the reaction plane, the transfer matrix in the ideal-detector limit acquires a straightforward Fourier SVD form even ahead of any angular discretization:

(24) d N d ϕ ′ = ∫ d ϕ T ( ϕ ′ | ϕ ) d N d ϕ ≃ ∑ n σ n cos ⁡ n ϕ ′ ∫ d ϕ cos ⁡ n ϕ d N d ϕ .

Here, ϕ and ϕ′ are angles relative to the actual and estimated reaction plane, respectively, and N is particle yield at a particular longitudinal velocity and transverse velocity magnitude. The singular values are σ _n = ⟨cos nΔΦ⟩, where ΔΦ is the deviation in the angle of the estimated plane from the actual. However, the low-velocity transfer matrix results from a complex interplay of detector performance, especially when azimuthally asymmetric, and the inferred distribution. Still, when the projectile is the same or nearly the same as the target, the distributions from deblurring in the velocity space near the projectile and target should physically be near mirror reflections of each other. This can serve as the mentioned test of deblurring validity, as the procedures’ details and complexity differ across regions, with the singular values in the SVD decomposition considered here never approaching unity in the low-velocity areas. Notably, before deblurring, the collision data have been quantified in terms of Fourier coefficients, as evident in (24), primarily n = 1 and 2.

Deep Learning might be used to image nuclear collisions, but there are currently obstacles in generating training data for the networks, no matter how the training is envisioned. The theoretical transport models of collisions are phenomenological and may miss essential features of the original images one seeks to access. They are also computationally expensive, making it difficult to reach the statistics needed for training. The situation may change once the essential features of the images have been settled, allowing for simplifying the training.

6 Summary

We examined the iterative RL deblurring method and partially other methods in the context of the utility of deblurring in nuclear and high-energy physics. We ignore explicit noise in the images for the time. SVD of the blurring matrix is an essential tool in understanding both the blurring and deblurring processes. The singular vectors generally exhibit a faster variation with pixel position for decreasing singular values. Symmetries of the blurring matrix generally give rise to symmetries of the singular vectors and may make the left- and right-hand vectors coincide. Blurring strength may be characterized by the falloff of singular values with a singular-vector index. A sufficiently strong blurring can give rise to a null space in the SVD. A one-dimensional model system of 12 pixels with periodic boundary conditions has served to illustrate our points. For high-contrast images, where the number of pixels at zero or near-zero intensity relative to the image’s maximal intensity is comparable to the null-space dimension, the nonnegativity constraint built into a deblurring method can help restore the null-space content in the processed image. For lower-contrast images, a regularization term that penalizes multiple extrema in the restored image helps control the null-space content. However, overly strong regularization may start affecting the restored image, such as blurring, which one is trying to correct. Despite the simplicity of the primary model we studied here, it can come close to the data analyses carried out in parallel in various respects.

For a clearer picture, a similar analysis of deblurring performance in the presence of noise in the blurred images is underway.

Corresponding author: Paweł Danielewicz, African Institute for Mathematical Sciences, Kigali, Rwanda; and Facility for Rare Isotope Beams and Department of Physics and Astronomy, Michigan State University, East Lansing, MI, USA, E-mail: daniel@frib.msu.edu

Paweł Danielewicz was partly supported by the U.S. Department of Energy Office of Science Grant No. DE-SC0019209.

References

1. D’Agostini, G. A multidimensional unfolding method based on Bayes’ theorem. Nucl Instrum Methods Phys Res 1995;362:487–98. https://doi.org/10.1016/0168-9002(95)00274-x.Search in Google Scholar

2. Danielewicz, P, Kurata-Nishimura, M. Deblurring for nuclei: 3D characteristics of heavy-ion collisions. Phys Rev C 2022;105:034608. https://doi.org/10.1103/physrevc.105.034608.Search in Google Scholar

3. Nzabahimana, P, Redpath, T, Baumann, T, Danielewicz, P, Giuliani, P, Guèye, P. Deconvoluting experimental decay energy spectra: the 26O case. Phys Rev C 2023;107:064315. https://doi.org/10.1103/physrevc.107.064315.Search in Google Scholar

4. Vargas, J, Benlliure, J, Caamaño, M. Unfolding the response of a zero-degree magnetic spectrometer from measurements of the Δ resonance. Nucl Instrum Methods Phys Res 2013;707:16–25. https://doi.org/10.1016/j.nima.2012.12.087.Search in Google Scholar

5. Xu, J, Qin, Y, Qin, Z, Si, D, Zhang, B, Wang, Y, et al.. Reconstruction of Bremsstrahlung γ-rays spectrum in heavy ion reactions with Richardson–Lucy algorithm. Phys Lett B 2024;857:139009. https://doi.org/10.1016/j.physletb.2024.139009.Search in Google Scholar

6. Christian Hansen, P, Nagy, JG, O’Leary, DP. Deblurring images, fundamentals of algorithms. Society for Industrial and Applied Mathematics; 2006. https://epubs.siam.org/doi/book/10.1137/1.9780898718874 [Accessed 15 Jun 2023].10.1137/1.9780898718874Search in Google Scholar

7. Reichel, L. Introduction to numerical computing I. Lecture notes; 2013. https://www.math.kent.edu/∼reichel/courses/intr.num.comp.1/syllabus.html [Accessed 9 Jun 2024].Search in Google Scholar

8. Zuninga, CD. Singular value decomposition for imaging applications. SPIE; 2021. https://spie.org/Publications/Book/2611523 [Accessed 3 Nov 2023].10.1117/3.2611523.ch1Search in Google Scholar

9. Lucy, LB. An iterative technique for the rectification of observed distributions. Astron J 1974;79:745. https://doi.org/10.1086/111605.Search in Google Scholar

10. Richardson, WH. Bayesian-based iterative method of image restoration. J Opt Soc Am 1972;62:55–9. https://doi.org/10.1364/josa.62.000055.Search in Google Scholar

11. Nicolas, D, Blanc‐Feraud, L, Zimmer, C, Roux, P, Kam, Z, Olivo‐Marin, J, et al.. Richardson–Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution. Microsc Res Tech 2006;69:260–6. https://doi.org/10.1002/jemt.20294.Search in Google Scholar PubMed

12. Khan, MK. Iterative methods of Richardson–Lucy-type for image deblurring. NMTMA 2013;6:262–75. https://doi.org/10.4208/nmtma.2013.mssvm14.Search in Google Scholar

13. Biemond, J, Lagendijk, RL, Mersereau, RM. Iterative methods for image deblurring. Proc IEEE 1990;78:856–83. https://doi.org/10.1109/5.53403.Search in Google Scholar

14. Landweber, L. An iteration formula for fredholm integral equations of the first kind. Am J Math 1951;73:615–24. https://doi.org/10.2307/2372313.Search in Google Scholar

15. Vankawala, F, Ganatra, A, Patel, A. A survey on different image deblurring techniques. Int J Comput Appl 2015;116:15–18. https://doi.org/10.5120/20396-2697.Search in Google Scholar

16. Trippe, T, Genzel, M, Macdonald, J, März, M. Let’s enhance: a deep learning approach to extreme deblurring of text images. Inverse Probl Imag 2023;17:1041–68. https://doi.org/10.3934/ipi.2023019.Search in Google Scholar

17. Nzabahimana, P, Danielewicz, P. Source function from two-particle correlation through deblurring. Phys Lett B 2023;846:138247. https://doi.org/10.1016/j.physletb.2023.138247.Search in Google Scholar

18. Tam, CK, Chajęcki, Z, Danielewicz, P, Nzabahimana, P. Source function from two-particle correlation function through entropy-regularized Richardson–Lucy deblurring. Phys Rev C 2025;112:024613.10.1103/zfly-38pkSearch in Google Scholar

19. Voloshin, S, Zhang, Y. Flow study in relativistic nuclear collisions by Fourier expansion of azimuthal particle distributions, en. Z für Physik C Part Fields 1996;70:665–71. https://doi.org/10.1007/s002880050141.Search in Google Scholar

Received: 2025-05-27

Accepted: 2025-12-08

Published Online: 2025-12-24

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2025-0252

Keywords for this article

deblurring; deconvolution; Richardson–Lucy algorithm; heavy-ion collisions; collective flow

Creative Commons

BY 4.0