# Optimal quadrature formulas for CT image reconstruction in the Sobolev space of non-periodic functions

In the present paper, optimal quadrature formulas in the sense of Sard are constructed for numerical integration of the integral ∫_a^b e^2π iω xφ(x)dx with ω∈R in the Sobolev space L_2^(m)[a,b] of complex-valued functions which are square integrable with m-th order derivative. Here, using the discrete analogue of the differential operator d^2m/d x^2m, the explicit formulas for optimal coefficients are obtained. The order of convergence of the obtained optimal quadrature formula is O(h^m). The optimal quadrature formulas of the cases m=2 and m=3 are applied for reconstruction of Computed Tomography (CT) images by approximating Fourier transforms in the filtered back-projection formula. In numerical experiments two kinds of phantoms are considered and numerical results are compared with the results of a built-in function of MATLAB 2019a, iradon. Numerical results show that the quality of the reconstructed images with optimal quadrature formula of the second and third orders is better than that of the results obtained by iradon.

• 3 publications
• 3 publications
• 6 publications
• 1 publication
07/30/2019

### On an optimal quadrature formula for approximation of Fourier integrals in the space L_2^(1)

This paper deals with the construction of an optimal quadrature formula ...
02/15/2021

### Optimal quadrature formulas for computing of Fourier integrals in a Hilbert space

In the present paper the optimal quadrature formulas in the sense of Sar...
02/15/2021

### On an optimal quadrature formula for approximation of Fourier integrals in the space W_2^(1,0)

The present paper is devoted to construction of an optimal quadrature fo...
07/31/2019

### Construction of optimal quadrature formulas exact for exponentional-trigonometric functions by Sobolev's method

The paper studies Sard's problem on construction of optimal quadrature f...
07/10/2022

### New Optimal Periodic Control Policy for the Optimal Periodic Performance of a Chemostat Using a Fourier-Gegenbauer-Based Predictor-Corrector Method

In its simplest form, the chemostat consists of microorganisms or cells ...
09/07/2022

### Numerical integration rules with improved accuracy close to singularities

Sometimes it is necessary to obtain a numerical integration using only d...
05/10/2020

### New conformal map for the trapezoidal formula for infinite integrals of unilateral rapidly decreasing functions

While the trapezoidal formula can attain exponential convergence when ap...

## 1 Introduction

It is known that when complete continuous X-ray data are available Computed Tomography (CT) image can be reconstructed exactly using the filtered back-projection formula (see, for instance, Buzug08 ; Feeman15 ; KakSlaney88 ). This formula gives interactions between the Radon transform, the Fourier transform and the back-projection transform. A description of the filtered back-projection formula along (KakSlaney88, , Chapter 3) is provided below.

In the Cartesian system with

-axes consider a unit vector

. Then the line perpendicular to this vector with the distance to the origin can be expressed as : . Assume the object is represented by a two variable function , which denotes the attenuation coefficient in X-ray CT applications. Then, the -view projection along the line can be expressed as

 P(t,θ)=∞∫−∞∞∫−∞μ(x,y)δ(xcosθ+ysinθ−t)dxdy,

where denotes the Dirac delta-function. The function is known as the Radon transform of . A projection is formed by combining a set of line integrals. The simplest projection is a collection of parallel ray integrals as is given by for a constant . This is known as a parallel beam projection. It should be noted that there are fan-beam in 2D and cone-beam in 3D projections Buzug08 ; Feeman15 ; KakSlaney88 .

The problem of CT is to reconstruct the function from its projections . There are analytic and iterative methods for CT reconstruction. One of the widely used analytic methods of CT reconstruction is the filtered back-projection method. It can be modeled by

 μ(x,y)=π∫0∞∫−∞S(ω,θ)|ω|e2πiω(xcosθ+ysinθ)dωdθ, (1.1)

where

 S(ω,θ)=∞∫−∞P(t,θ)e−2πiωtdt (1.2)

is the 1D Fourier transform of . The inner integral of (1.1),

 Q(t,θ)=∞∫−∞S(ω,θ)|ω|e2πiωtdω, (1.3)

is a 1D inverse Fourier transform of the product , which represents a projection filtered by a 1D filter whose frequency representation is . The outer integral performs back-projection. Therefore, the filtered back-projection consists of two steps: filtration and then back-projection.

Thus, the Fourier transforms play the main role in (1.1)-(1.3). But in practice, due to the fact that we have discrete values of the Radon transform, we have to approximately calculate the Fourier transforms in the filtered back-projection. For this purpose, it is necessary to consider the problem of approximate calculation of the integral

 I(φ)=b∫ae2πiωxφ(x)dx (1.4)

with . This type of integrals are called highly oscillating integrals. In most cases it is impossible to get the exact values of such integrals. Thus, they can be approximately calculated using the formulas of numerical integration. However, standard methods of numerical integration cannot be successfully applied for that. Therefore special effective methods should be developed for approximation of highly oscillating integrals. One of the first numerical integration formula for the integral (1.4) was obtained by Filon Filon28 in 1928 using a quadratic spline. Since then, for integrals of different types of highly oscillating functions many special effective methods have been developed, such as Filon-type method, Clenshaw-Curtis-Filon type method, Levin type methods, modified Clenshaw-Curtis method, generalized quadrature rule, and Gauss-Laguerre quadrature (see, for example, AvdMal89 ; IBab ; BabVitPrag69 ; BakhVas68 ; IserNor05 ; Mil98 ; NovUllWoz15 ; XuMilXiang ; ZhangNovak19 , for more review see, for instance, DeanoHuyIser18 ; MilStan14 ; Olver08 and references therein).

Recently, in BolHayShad16 ; BolHayShad17 ; BolHayMilShad17 , based on Sobolev’s method, the problem of construction of optimal quadrature formulas in the sense of Sard for numerical calculation of integrals (1.4) with integer was studied in Hilbert spaces and , respectively.

Here, we consider the Sobolev space of non-periodic, complex-valued functions defined on the interval which posses an absolute continuous st derivative on , and whose th order derivative is square integrable Sobolev74 ; SobVas . The space is a Hilbert space with the inner product

 ⟨φ,ψ⟩=b∫aφ(m)(x) ¯ψ(m)(x)dx, (1.5)

where is the -th order derivative of the function with respect to , is the complex conjugate function to the function and the norm of the function is correspondingly defined by the formula

 ∥φ∥L(m)2[a,b]=⟨φ,φ⟩1/2.

The aim of the present work is to construct optimal quadrature formulas in the sense of Sard in the Sobolev space for numerical integration of the integral (1.4) with real using the results of the work BolHayShad17 . Then to apply the obtained optimal quadrature formulas to the approximate reconstruction of CT image employing the filtered back-projection method (1.1). Recently we got the results when HayJeonLee19 .

The rest of the paper is organized as follows. Section 2 is devoted to construction of optimal quadrature formulas in the sense of Sard in the space for numerical calculation of Fourier integrals. There are obtained analytic formulas for optimal coefficients using the discrete analogue of the differential operator . In section 3, the obtained optimal quadrature formulas for the cases and are applied for CT image reconstruction by approximating Fourier transforms in the filtered back-projection formula.

## 2 Optimal quadrature formulas for Fourier integrals in the space L(m)2[a,b]

In space for approximation of the integral (1.4), we consider the following quadrature formula

 b∫ae2πiωxφ(x)dx≅N∑β=0Cβφ(xβ) (2.1)

with the error

 (ℓ,φ)=b∫ae2πiωxφ(x)dx−N∑β=0Cβφ(xβ), (2.2)

where is the value of the error functional at the function . Here the error functional has the form

 ℓ(x)=e2πiωxε[a,b](x)−N∑β=0Cβδ(x−xβ), (2.3)

are coefficients, () are nodes of the formula (2.1), , , , ,

is the characteristic function of the interval

, and is the Dirac delta-function. We mention that the coefficients depend on , and .

The error of the formula (2.1) defines a linear functional in , where is the conjugate space to the space . Since the functional (2.3) is defined on the space , the conditions

 (ℓ,xα)=0, α=0,1,...,m−1 (2.4)

should be fulfilled. The conditions (2.4) mean the exactness of the quadrature formula (2.1) for all algebraic polynomials of degree less than or equal to . Hence we get that for a function from the space the order of convergence of the optimal quadrature formula (2.1) is .

It should be mentioned that from equalities (2.4) one can get the condition for existence of the optimal quadrature formula of the form (2.1) in the space .

Sard’s optimization problem of numerical integration formulas of the form (2.1) in the space is the problem of finding the minimum of the norm of the error functional by coefficients , i.e., to find coefficients satisfying the equality

 ∥˚ℓ∥L(m)∗2[a,b]=infCβ∥ℓ∥L(m)∗2[a,b]. (2.5)

The coefficients satisfying the last equality are called optimal coefficients and they are denoted as . The quadrature formula with coefficients is called the optimal quadrature formula in the sense of Sard, and is the error functional corresponding to the optimal quadrature formula.

The solution of Sard’s problem gives the sharp upper bound for the error (2.2) of functions from the space as follows

 |(ℓ,φ)|≤∥φ∥L(m)2[a,b]∥˚ℓ∥L(m)∗2[a,b].

This problem, for the quadrature formulas of the form (2.1) with , was first studied by Sard Sard in the space for some . Since then, it was investigated by many authors (see, for instance, CatCom ; Com72a ; Com72b ; GhOs ; FLan ; MalOrl ; ShadHay11 ) using spline method, -function method, and Sobolev’s method. Finally, in the works Koh ; Shad83 ; Shad99 this problem was solved for any with equally spaced nodes and the explicit expressions for the optimal coefficients have been obtained (see Theorem 2.4).

It should be noted that the problem of construction of lattice optimal cubature formulas in the space of multi-variable functions was first stated and investigated by S.L. Sobolev Sobolev74 ; SobVas . Further, in this section, based on the results of the work BolHayShad17 , we solve Sard’s problem on construction of optimal quadrature formulas of the form (2.1) for with , first for the interval

and then using a linear transformation for the interval

. For this we use the following auxiliary results.

### 2.1 Preliminaries

We need the concept of discrete argument functions and operations on them. For this we give the definition for functions of discrete argument by following Sobolev74 ; SobVas .

Assume that the nodes are equally spaced, i.e., is a small positive parameter, and and are complex-valued functions defined on the real line or on an interval of . The function given on some set of integer values of is called a function of discrete argument. The inner product of two discrete argument functions and is defined by

 [φ(hβ),ψ(hβ)]=∞∑β=−∞φ(hβ)⋅¯ψ(hβ)

if the series on the right hand side of the last equality converges absolutely. The convolution of two functions and is the following inner product

 φ(hβ)∗ψ(hβ)=[φ(hγ),ψ(hβ−hγ)]=∞∑γ=−∞φ(hγ)⋅¯ψ(hβ−hγ).

We note that coefficients of optimal quadrature formulas and interpolation splines in the spaces

and depend on the roots of the Euler-Frobenius type polynomials (see, for instance, BabHay19 ; BolHayShad17 ; CabHayShad14 ; Koh ; Shad10 ; ShadHay11 ; ShadHay13 ; ShadHay14 ; ShadHayAkhm15 ; ShadHayNur13 ; ShadHayNur16 ; Sobolev06 ; SobVas ). The Euler-Frobenius polynomials , are defined as follows (see, for instance, Frob1910 ; SobVas ):

 Ek(x)=(1−x)k+2x(xddx)kx(1−x)2, k=0,1,2,.... (2.6)

The coefficients of the Euler-Frobenius polynomial of degree are expressed by the following formula which was obtained by Euler:

 as=s∑j=0(−1)j(k+2j)(s+1−j)k+1.

In Frob1910 it was shown that all roots , of the polynomial are real, negative and distinct, that is:

 q1

Furthermore, these roots satisfy the relation

 qj⋅qk+1−j=1, j=1,2,...,k.

For the Euler-Frobenius polynomials the following identity holds

 Ek(x)=xkEk(1x), (2.7)

and also the following is true.

###### Theorem 2.1

(Lemma 3 of Shad10 ). Polynomial which is defined by the formula

 Qk(x)=(x−1)k+1k+1∑i=0Δi0k+1(x−1)i (2.8)

is the Euler-Frobenius polynomial (2.6) of degree , i.e., , where .

We also use the formula

 n−1∑γ=0qγγk=11−qk∑i=0(q1−q)iΔi0k−qn1−qk∑i=0(q1−q)iΔiγk|γ=n, (2.9)

which is given in Ham62 , where is the finite difference of order of and is the ratio of a geometric progression.

In finding the analytic formulas for coefficients of optimal formulas in the space by Sobolev method the discrete analogue of the operator plays the main role. This discrete analogue satisfies the equality

 hDm(hβ)∗Gm(hβ)=δd(hβ), (2.10)

where is the discrete argument function for the function

 Gm(x)=|x|2m−12(2m−1)!, (2.11)

and is equal to 0 when and 1 when .

We note that the operator was introduced and studied by S.L. Sobolev Sobolev74 . In Shad85 the discrete function was constructed and the following was proved.

###### Theorem 2.2

The discrete analogue of the differential operator has the form

 Dm(hβ)=p⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩m−1∑k=1Akq|β|−1k for |β|≥2,1+m−1∑k=1Ak for |β|=1,C+m−1∑k=1Akqk for β=0, (2.12)

where

 p=(2m−1)!h2m,Ak=(1−qk)2m+1E2m−1(qk),C=−22m−1, (2.13)

is the Euler-Frobenius polynomial of degree , are the roots of the Euler-Frobenius polynomial , , and is a small positive parameter.

In addition, some properties of were studied in the works Shad85 ; Sobolev74 . Here we give the following.

###### Theorem 2.3

The discrete argument function and the monomials are related to each other as follows:

 ∞∑β=−∞Dm(hβ)(hβ)k={0\it for 0≤k≤2m−1,(2m)!\it for k=2m. (2.14)

Now we give the results of the works Shad83 ; Shad10 on the optimal quadrature formulas of the form (2.1) in the sense of Sard and on the norm of the optimal error functional corresponding to the case .

###### Theorem 2.4

Coefficients of the optimal quadrature formulas of the form

 1∫0φ(x)dx≅N∑β=0Cβφ(hβ) (2.15)

in the space have the form

 ˚C0 = h(12−m−1∑k=1dkqk−qNk1−qk), ˚Cβ = h(1+m−1∑k=1dk(qβk+qN−βk)), β=1,2,...,N−1, (2.16) ˚CN = h(12−m−1∑k=1dkqk−qNk1−qk),

where satisfy the system

 m−1∑k=1dkj∑i=1qk+(−1)i+1qN+ik(qk−1)i+1Δi0j=Bj+1j+1, j=1,2,...,m−1,

are roots of the Euler-Frobenius polynomial of degree with , and .

###### Theorem 2.5

The norm of the error functional of the optimal quadrature formula (2.15) on the space has the form

 ∥˚ℓ∥2L(m)∗2[0,1]=(−1)m+1(h2mB2m(2m)!+2h2m+1(2m)!m−1∑k=1dk2m∑i=1−qN+ik+(−1)iqk(1−qk)i+1Δi02m),

where is the Bernoulli number, and , and are given in Theorem 2.4.

### 2.2 Construction of optimal quadrature formulas for the interval [0,1]

Here we obtain optimal quadrature formulas of the form (2.1) for the interval when and . In the space , using the results of Sections 2, 3 and 5 of BolHayShad17 , for the coefficients of the optimal quadrature formulas in the sense of Sard of the form

 1∫0e2πiωxφ(x)dx≅N∑β=0Cβφ(hβ) (2.17)

for with , we get the following system of linear equations

 N∑γ=0CγGm(hβ−hγ)+Pm−1(hβ)=fm(hβ), β=0,1,...,N, (2.18) N∑γ=0Cγ(hγ)α=gα, α=0,1,...,m−1, (2.19)

where is a polynomial of degree with complex coefficients,

 fm(hβ) = 1∫0e2πiωxGm(x−hβ)dx = −2m−1∑α=0(hβ)2m−1−α(−1)αgα2α!(2m−1−α)!+e2πiωhβ(2πiω)2m−2m−1∑k=0(hβ)2m−1−k(2m−1−k)!(2πiω)k+1, gα = 1∫0e2πiωxxαdx = α−1∑k=0(−1)kα! e2πiω(α−k)!(2πiω)k+1+(−1)αα!(2πiω)α+1(e2πiω−1), α=0,1,...,

is defined by (2.11), and for .

In the system (2.18)-(2.19) unknowns are the optimal coefficients and , . We point out that when the system (2.18)-(2.19) has a unique solution. This solution satisfies conditions (2.4) and the equality (2.5). It should be noted that the existence and uniqueness of the solution for such type of systems were studied, for example, in ShadHayAkhm15 ; Sobolev74 ; SobVas .

We are interested in finding explicit formulas for the optimal coefficients , and unknown polynomial satisfying the system (2.18)-(2.19). The system (2.18)-(2.19) is solved similarly as the system (34)-(35) of BolHayShad17 by Sobolev’s method, using the discrete analogue of the differential operator .

We formulate the results of this section as the following two theorems.

###### Theorem 2.6

For real with , the coefficients of optimal quadrature formulas of the form (2.1) in the space when are expressed by formulas

 ˚C0 = h(e2πiωhKω,me2πiωh−1−12πiωh+m−1∑k=1(akqkqk−1+bkqNk1−qk)), (2.22) ˚Cβ = h(e2πiωhβKω,m+m−1∑k=1(akqβk+bkqN−βk)), β=1,2,...,N−1, (2.23) ˚CN = h(e2πiωKω,m1−e2πiωh+e2πiω2πiωh+m−1∑k=1(akqNk1−qk+bkqkqk−1)), (2.24)

where and are defined by the following system of linear equations

 m−1∑k=1ak[j∑t=1qkΔt0j(qk−1)t+1]+m−1∑k=1bk[j∑t=1qN+tkΔt0j(1−qk)t+1]=j!(2πiωh)j+1−j∑t=1e2πiωhKω,mΔt0j(e2πiωh−1)t+1, j=1,2,...,m−1,m−1∑k=1ak[j∑t=1qtkΔt0j(1−qk)t+1−j∑α=1hα−j(jα)α∑t=1qN+tkΔt0α(1−qk)t+1]+m−1∑k=1bk[j∑t=1qN+1kΔt0j(qk−1)t+1−j∑α=1hα−j(jα)α∑t=1qkΔt0α(qk−1)t+1]=(−1)j+1j!(2πiωh)j+1+j∑α=1hα−j(−1)αj!e2πiω(j−α)!(2πiωh)α+1−Kω,m1−e2πiωhj∑t=1(e2πiωh1−e2πiωh)tΔt0j+Kω,me2πiω1−e2πiωhj∑α=1hα−j(jα)α∑t=1(e2πiωh1−e2πiωh)tΔt0α,  j=1,2,...,m−1, (2.25)

are roots of the Euler-Frobenius polynomial with , and

 Kω,m=(sinπωhπωh)2m(2m−1)!2m−2∑α=0aαcos[2πωh(m−1−α)]+am−1. (2.26)

Here are the coefficients of the Euler-Frobenius polynomial of degree .

###### Theorem 2.7

For with , the coefficients of optimal quadrature formulas of the form (2.1) in the space when are expressed by formulas

 ˚C0 = h(−12πiωh+m−1∑k=1(akqkqk−1+bkqNk1−qk)), ˚Cβ = hm−1∑k=1(akqβk+bkqN−βk),    β=1,2,...,N−1, ˚CN = h(e2πiω2πiωh+m−1∑k=1(akqNk1−qk+bkqkqk−1)),

where and , , are defined by the following system of linear equations:

 m−1∑k=1ak[j∑t=1qkΔt0j(qk−1)t+1]+m−1∑k=1bk[j∑t=1qN+tkΔt0j(1−qk)t+1]=j!(2πiωh)j+1, j=1,2,...,m−1,m−1∑k=1ak[j∑t=1qtkΔt0j(1−qk)t+1−j∑α=1hα−j(jα)α∑t=1qN+tkΔt0α(1−qk)t+1]+m−1∑k=1bk[j∑t=1qN+1kΔt0j(qk−1)t+1−j∑α=1hα−j(jα)α∑t=1qkΔt0α(qk−1)t+1]=(−1)j+1j!(2πiωh)j+1+j∑α=1hα−j(−1)αj!e2πiω(j−α)!(2πiωh)α+1,  j=1,2,...,m−1.

Here are the roots of the Euler-Frobenius polynomial and .

We note that Theorem 2.6 is generalization of Theorem 6 in [6] for real with while Theorem 2.7 for with is the same with Theorem 7 of the work [6]. Theorem 6 is proved similarly as Theorem 6 of [6]. Therefore, it is sufficient to give a brief proof of Theorem 2.6.

The brief proof of Theorem 2.6. First, such as in Theorem 5 of the work BolHayShad17 , using the discrete function , for optimal coefficients , when is real and we get the following formula

 ˚Cβ=h(Kω,me2πiωhβ+m−1∑k=1(akqβk+bkqN−βk)),  β=1,2,...,N−1, (2.27)

where , , , are unknowns, and are roots of the Euler-Frobenius polynomial of degree with .

Next, it is sufficient to find , , , optimal coefficients , and unknown polynomial of degree . Now putting the form (2.27) of optimal coefficients , , into (2.18), using (2.7)-(2.9), after some simplifications, we get the following identity with respect to

 e2πiωhβ(h2mKω,m(2m−1)!2m−1∑i=0e2πiωhΔi02m−1(e2πiωh−1)i+1)+(hβ)2m−1(2m−1)!(˚C0−Kω,mhe2πiωhe2πiωh−1−hm−1∑k=1akqk−bkqNkqk−1)−2m−1∑j=1(hβ)2m−1−jj!(2m−1−j)!hj+1(j∑i=0Kω,me2πiωhΔi0j(e2πiωh−1)i+1+m−1∑k=1j∑i=0