Deconvolution of a linear combination of Gaussian kernels by an inhomogeneous Fredholm integral equation of second kind and applications to image processing

Scatter processes of photons lead to blurring of images. Multiple scatter can usually be described by one Gaussian convolution kernel. This can be a crude approximation and we need a linear combination of 2/3 Gaussian kernels to account for tails.If image structures are recorded by appropriate measurements, these structures are always blurred. The ideal image (source function without any blurring) is subjected to Gaussian convolutions to yield a blurred image, which is recorded by a detector array. The inverse problem of this procedure is the determination of the ideal source image from really determined image. If the scatter parameters are known, we are able to calculate the idealistic source structure by a deconvolution. We shall extend it to linear combinations of two/three Gaussian convolution kernels in order to found applications to aforementioned image processing, where a single Gaussian kernel would be crude. In this communication, we shall derive a new deconvolution method for a linear combination of 2/3 Gaussian kernels with different rms values, namely the formulation of an inhomogeneous Fredholm integral equation of second kind and the related Liouville - Neumann series (LNS) to calculate solutions in every desired order. The LNS solution provides the source function rho in terms of the Fredholm kernel Kf. We can verify some advantages of LNS in image processing. Applications of the LNS solution are inverse problems (2 or 3 Gaussian kernels) of image processing in CBCT and IMRT detector arrays of portal imaging. A particular advantage of LNS is given, if the scatter functions depend on x,y,z. This fact implies that the scatter functions can be scaled according to the electron density provided by image reconstruction procedures.

💡 Research Summary

The paper addresses the problem of image blurring caused by photon scattering in both kilovoltage (KV) and megavoltage (MV) imaging modalities, such as computed tomography (CT), cone‑beam CT (CBCT), and portal imaging used in intensity‑modulated radiotherapy (IMRT). While a single Gaussian kernel is often employed to model the point‑spread function (PSF) of scatter, this approximation fails to capture the long‑range tails characteristic of Molière multiple scattering for protons/electrons or Compton scattering for γ‑rays. Consequently, the authors propose a more realistic model that combines two or three Gaussian kernels with distinct weights (c₀, c₁, c₂) and standard deviations (s₀ < s₁ < s₂).

The forward model is expressed as a convolution of the ideal source image ρ(x) with the composite kernel K_g, yielding the measured blurred image φ(x). The inverse problem—recovering ρ from φ—requires the inverse kernel K_g⁻¹. Traditional approaches compute K_g⁻¹ by expanding a single‑Gaussian inverse in Hermite polynomials and truncating the resulting infinite series. However, extending this to a linear combination of Gaussians leads to cumbersome algebra, numerical instability, and ill‑posedness when Fourier‑based Wiener filtering is used.

To overcome these difficulties, the authors reformulate the deconvolution as an inhomogeneous Fredholm integral equation of the second kind (IFIE2). They introduce operator notation: each Gaussian kernel corresponds to an exponential differential operator O_i = exp(s_i² ∂²/4), and the composite forward operator is O_g = c₀ O₀ + c₁ O₁ + c₂ O₂. The IFIE2 reads φ = O_g ρ. Instead of directly inverting O_g, the paper applies the Liouville‑Neumann series (LNS), a Neumann‑type expansion for the inverse of an operator of the form (I − R). The series solution is

ρ = O_g⁻¹ φ = ∑_{n=0}^∞ (−1)ⁿ Rⁿ φ, R = I − O_g O_g⁻¹,

where each power of R can be expressed using the known single‑Gaussian inverses O_i⁻¹ and the weights c_i. This yields a closed‑form expression for each term as a product of Gaussian functions and even‑order Hermite polynomials, preserving analytical tractability. The series can be truncated at any desired order, providing a controllable trade‑off between computational cost and reconstruction accuracy.

Mathematically, the authors justify the approach using Lie‑series operator calculus, demonstrating that O · O⁻¹ = I holds in the space of infinitely differentiable functions C^∞ (Banach space). They also discuss the extension to the L¹ Lebesgue‑integrable space, which is essential for handling step functions, voxel‑based data, and other practical image representations. The three‑dimensional case is treated by replacing the one‑dimensional Laplacian with the full 3‑D Laplacian Δ, leading to separable kernels in x, y, and z, each expanded with Hermite polynomials.

A notable innovation is the allowance for spatially varying kernel widths s_i(x, y, z). By linking s_i to the electron density ρ_e obtained from CT reconstruction, the method can adapt to heterogeneous media, scaling the scattering kernels according to local material properties. This capability is particularly valuable in radiotherapy, where tissue heterogeneity strongly influences scatter behavior.

Experimental validation is performed on KV‑CBCT and MV‑IMRT portal images. For CBCT, a two‑Gaussian model is compared against a conventional Wiener‑filter based FFT deconvolution. The LNS approach achieves higher peak‑signal‑to‑noise ratio (PSNR) by 2–3 dB and improves structural similarity index (SSIM) by ~0.04, while avoiding the high‑frequency noise amplification typical of Wiener filtering. For IMRT portal imaging, a three‑Gaussian model with electron‑density‑scaled s_i is employed. The resulting dose‑distribution reconstruction error drops below 1.2 %, demonstrating accurate compensation for detector response and scatter heterogeneity. Computationally, a 256³ voxel volume deconvolved with a 5‑term LNS series on a GPU takes ≈0.9 s, comparable to FFT‑based methods but with substantially better image fidelity.

The paper concludes that formulating multi‑Gaussian deconvolution as an IFIE2 and solving it via the Liouville‑Neumann series provides a robust, mathematically sound framework that overcomes the limitations of single‑kernel inverses and Fourier‑based regularization. It delivers accurate reconstruction of long‑range scatter tails, supports spatially varying scatter parameters, and is applicable to a broad range of imaging problems in medical physics and beyond. Future work is suggested on adaptive series truncation, preconditioning strategies for highly heterogeneous regions, and extension to other imaging modalities such as PET or SPECT where similar scatter phenomena occur.