Image Reconstruction with Maclaurin Series Expansion

TL;DR

This paper explores a theoretical method for image reconstruction using a Maclaurin series expansion in the Fourier domain under idealized conditions.

Contribution

The novel approach uses a Maclaurin series to reconstruct images from minimal angular scanning data without prior knowledge or training.

Findings

The Maclaurin series expansion converges across the entire Fourier space under ideal conditions.

Computer simulations demonstrate successful 2D image reconstruction using a truncated Fourier-domain series.

The method requires minimal data and avoids the need for training or prior knowledge.

Abstract

This is a forward-thinking theoretical investigation and may not have practical values for current imaging systems. This investigation assumes that there is no noise in the measurements, the signals are continuous (not sampled), the computer has perfect precession, and there are no round-off errors. Under these unrealistic conditions, we can form a Maclaurin series expansion in the Fourier domain with measurements in a small scanning angular range. We show that this Maclaurin series expansion converges in the entire Fourier space. As a result, a complete data set is available for image reconstruction. The Fourier domain is complex; the expansion coefficients are most likely complex with real parts and imaginary parts. Computer simulations are performed to illustrate a 2D spatial-domain image can be obtained if a Fourier-domain truncated Maclaurin series expansion is available. Our goal is to use minimum data for trust-worthy reconstruction without any prior knowledge and training data.

Figures

Full text

Introduction

There are many data sufficiency conditions for various imaging modalities and imaging geometries [1–15]. For example, in two-dimensional (2D) imaging, the parallel-beam system requires a scanning angular range of 180°. The fan-beam system requires a scanning angular range of 180° plus the fan angle. If the scanning trajectories satisfy the data sufficiency conditions, we have stable image reconstruction algorithms, which can be analytical or iterative.

Even when the angular range satisfied the data sufficiency conditions, the sampling on the detector may not be adequate. If the detector is not large enough to cover the object to be imaged, the projection data is truncated, resulting in an under-sampling situation [16–33]. Another under-sampling situation is that the angular sampling is not dense enough, which is also known as few-view tomography [33–49].

When data is insufficient, some other assumes can make the inverse problem solvable. One of such situations is compressed sensing [50–60]. The compressed sensing methods consider the inverse problem solutions, which are sparse, that is, most of the elements are zero. In tomographic application, the solutions $[eqn]$ can be assumed as piecewise constant. The derivative of $[eqn]$ of the finite difference of $[eqn]$ is a sparse image. The compressed sensing theory suggests that a usable sparse solution can be obtained by minimization of the $[eqn]$ norm of the sparse solution $[eqn]$ . For a piecewise-constant solution, we minimize the $[eqn]$ norm of the finite difference of $[eqn]$ . The $[eqn]$ norm minimization is not an easy task, because the $[eqn]$ norm of an image is the total count of non-zero image pixels. The gradient of this total count with respect to each pixel does not exist. The gradient-based optimization algorithms do not work.

A practical work around is to use the $[eqn]$ norm to approximate the $[eqn]$ norm. For a piecewise-constant image, this remedy leads to the total variation (TV) minimization algorithms. In this paper, we do not assume that the images are piecewise constant. We do not treat the image reconstruction problem as a compressed sensing problem. We do not discretize the imaging problem. We assume the image and its projection measurements are continuous.

At the beginning of this section, we state that a data sufficiency condition for a 2D parallel-beam imaging problem is to acquire data over an angular range of 180°. We argue that this data sufficiency condition is not necessary. This condition is derived based on the assumption that the entire Fourier space must be completely measured. We will show in the next section that we do not have to measure the entire Fourier space to capture the complete information about the spatial-domain object. In fact, we only need to know the Fourier transform at one location, for example, at the center.

This paper assumes a perfect ideal world, where the detected signals are continuous and noiseless. As will be shown in the next section, perfect high-order mixed partial derivatives can be calculated by the measurements. All these assumptions guarantee that we can form a 2D Taylor series expansion in the Fourier domain of the object. This paper is forward-thinking, and we may not be able to implement the proposed system using today’s technology. It is likely that we are unable to do computer simulations for the proposed system due to the round-off errors and finite word-length limits in a practical computer. The errors may cause the Taylor series expansion unstable or even diverge.

For line-integral based imaging systems such as x-ray computer tomography (CT), positron emission tomography (PET), and single photon emission computed tomography (SPECT), it is challenging to estimate the derivatives in the Fourier domain. We will use the Central Slice Theorem to suggest some potential strategies in the next section.

Methods

In mathematics, a “compact support” refers to a function that is only non-zero within a bounded, closed set (a compact set), meaning it becomes zero outside of that specific region. In medical imaging, the patient’s body always has compact support.

If a function has a compact support, the Fourier transform of this function is an “entire function” (also known as holomorphic function and analytic function). An entire function has many desirable properties. An entire function can be differentiated with any order at any point in the complex plane and has no singularities. The Taylor series expansion of an entire function converges everywhere in the complex plane.

Let $[eqn]$ be a 2D real function that has finite support. Assume that $[eqn]$ vanishes for $[eqn]$ and $[eqn]$ According to the Paley-Wiener theorem [61], for a square-integrable function with a finite support, its Fourier transform is holomorphic.

The 2D Fourier transform of $[eqn]$ is given as

$$[eqn]$$

Interchanging of integration and summation is allowed because $[eqn]$ is supported on a finite region and because the series for the exponential function converges absolutely and uniformly [61].

From (1), we have

$$[eqn]$$

We also have

$$[eqn]$$ $$[eqn]$$

An alternative expression for this Maclaurin series expansion is

$$[eqn]$$

From the right-hand side of (3), we learn that if a mixed partial derivative is evaluated in the image domain, the whole image $[eqn]$ must be used in the calculation. Unfortunately, we do not know of any imaging system that can measure or calculate the double integral in (3).

Our innovative strategy is to use the directional derivatives to estimate the mixed partial derivatives. Let $[eqn]$ be a rotated version of $[eqn]$ with a rotation angle $[eqn]$ . The Radon transform (i.e., the parallel-beam projections) is given by

$$[eqn]$$

The Fourier transform of $[eqn]$ with respect to variable $[eqn]$ is

$$[eqn]$$

where $[eqn]$ is the frequency variable corresponding to $[eqn]$ . Let the 2D Fourier transform of $[eqn]$ be $[eqn]$ . According to the Central Slice Theorem [1], $[eqn]$ is a central slice of $[eqn]$ as illustrated in Figure 1.

The $[eqn]$ th-order directional derivative of $[eqn]$ can be calculated as the Fourier transform of the $[eqn]$ th moment

$$[eqn]$$

It can be shown that the $[eqn]$ th-order directional derivative of $[eqn]$ and the $[eqn]$ th-order mixed partial derivatives $[eqn]$ are related as

$$[eqn]$$

We can prove (9) using the mathematical induction method. The version of (9) for $[eqn]$ is

$$[eqn]$$

which can be immediately obtained by using the definition of the directional derivative when the partial derivatives $[eqn]$ and $[eqn]$ exist.

We now assume that (9) is valid for an integer $[eqn]$ . Then

$$[eqn]$$

In the derivation above, we used

$$[eqn]$$

which can be readily verified by definition.

To recap, Eq. (9) is the relationship between the $[eqn]$ th-order directional derivative and the $[eqn]$ th-order mixed partial derivatives.

Our goal is to form a 2D truncated Taylor series expansion. This expansion requires mixed partial derivatives to construct its coefficients.

A $[eqn]$ th-order directional derivative of the Fourier transform of Radon transform $[eqn]$ can be calculated by the Fourier transform of the $[eqn]$ th moment by Eq. (8).

There are $[eqn]$ mixed partial derivatives of order $[eqn]$ . Therefore, we need to measure $[eqn]$ at $[eqn]$ different angles: $[eqn]$ , and solve a system of linear equations:

$$[eqn]$$

where

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

with $[eqn]$ and $[eqn]$ .

The Taylor coefficients are then determined by the measurements at these $[eqn]$ angles: $[eqn]$ . These angles do not have to be uniformly distributed over 180°. They can, for example, be densely distributed in a very small angular range (say, 40°).

Thus, if an imaging system is able to measure all information of the object in the Fourier domain at (0, 0), including higher-order mixed partial derivatives at (0, 0), a Maclaurin series expansion can be formed; this expansion converges everywhere in the complex plane. Since this expansion converges, a truncated expansion (with finite number of terms) can be used to approximate the Fourier-domain function $[eqn]$ .

The next step is to evaluate the truncated Maclaurin expansion of $[eqn]$ at any location ( $[eqn]$ ). For example, we can evaluate the expansion at a regular grid of ( $[eqn]$ ). Finally, we perform the 2D inverse Fourier transform to obtain the reconstructed image $[eqn]$ .

The “entire function” idea presented above does not apply in the spatial domain. The spatial domain image cannot be an “entire function” because the spatial domain image has a finite support and thus the image must have discontinuities. Discontinuity prevents the spatial domain images from being differentiable everywhere.

In order to gain some intuitive about the feasibility whether a truncated Taylor series expansion is useful in image reconstruction, the next section will present a computer simulation example using an imperfect computer, which has limitations such as a finite-word length and round-off errors.

A square phantom $[eqn]$ has a known, close form 2D Fourier transform $[eqn]$ , which is a product of two sinc functions:

$$[eqn]$$

where the parameter $[eqn]$ determines the width of the square. It is known that the Maclaurin series expansion for the sine function

$$[eqn]$$

This immediately gives

$$[eqn]$$

Therefore, (18) becomes

$$[eqn]$$

A truncated version of (21) is

$$[eqn]$$

Results

A 1D Fourier transform pair is shown Figure 2, where the left diagram is in the 1D Fourier domain, and the right diagram is in the spatial domain. The left diagram contains two curves. The blue solid curve is a section of the sinc function. The orange broken curve is a truncated Maclaurin series expansion approximation with 50 terms of (20) and $[eqn]$ . The right diagram also contains two curves. The blue solid curve is the 1D inverse Fourier transform of the section of the sinc function shown in the left diagram. The broken orange curve is the 1D inverse Fourier transform of the truncated Maclaurin series expansion shown in the left diagram. The 1D invers Fourier transform was implemented in MATLAB numerically using 128 discrete samples over [−1.5, 1.5] for the variable $[eqn]$ . The MATLAB function ‘ifft’ was used in the computer simulation.

This 1D computer simulation demonstrates that the 50-term truncated Maclaurin series expansion is pretty close to the sinc function and that the discrete inverse Fourier transform can provide a pretty close approximation in the spatial domain. In other words, the “ifft” is not ill-conditioned.

A 2D Fourier transform pair is shown Figure 3, where the left diagram is $[eqn]$ in the 2D Fourier domain and the right diagram is its 2D inverse Fourier transform in the spatial domain. The object is the difference of two squares. The bigger square has a parameter $[eqn]$ and intensity of 1. The smaller square has a parameter $[eqn]$ and intensity of 0.7. Thus, $[eqn]$ is the difference of two 2D sinc functions and is calculated using MATLAB’s built-in sine function using discrete samples of $[eqn]$ and $[eqn]$ , in a $[eqn]$ regular grid. Taking MATLAB’s 2D inverse Fast Fourier Transform ‘ifft2,’ the left image $[eqn]$ becomes the right image $[eqn]$ .

Figure 4 is almost the same as Figure 3, except that the Fourier domain image on the left is a truncated version of the Taylor series expansion (22) with a maximum value $[eqn]$ being 50.

This 2D computer simulation demonstrates that the truncated Maclaurin series expansion is pretty close to the 2D sinc function and that the discrete 2D inverse Fourier transform can provide a pretty close approximation in the spatial domain. In other words, the “ifft2” is not ill-conditioned.

Discussion and Conclusions

In the spatial domain, the patient image $[eqn]$ is a square-integrable function with a finite support. Its 2D Fourier transform $[eqn]$ is an entire function on the 2D complex plane. A Taylor series expansion at $[eqn]$ converges everywhere in the complex Fourier space. The function $[eqn]$ can be evaluated with this Taylor series expansion anywhere in the complex Fourier space. The Taylor series expansion coefficients only depend on the local values of $[eqn]$ around $[eqn]$ . In other words, the spatial domain image $[eqn]$ can be reconstructed by the local information of $[eqn]$ in the Fourier domain. When $[eqn]$ , the expansion is referred to as the Maclaurin series expansion.

This paper assumes that the coefficients of the Taylor series expansion can be somehow exactly measured and calculated by a future hypothetical scanner. Then a truncated expansion can be formed and used for image reconstruction.

This paper considers image reconstruction only from the measurements of the object itself. Of course, prior information about the object can be useful to supplement unmeasured data. Machine learning is an excellent example of how prior information can make image reconstruction possible when measurements are incomplete. Some other prior information does not need neural network training, such as using the total variation (TV) norm minimization.

The robustness of the Taylor series expansion method is important but is beyond the scope of this paper. The singular value decomposition method has been used to study the ill-condition of limited angle tomography [62].