Analytical formula for two-dimensional ring artefact suppression

Titarenko, V.

doi:10.1107/S160057751601479X

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 23| Part 6| November 2016| Pages 1447-1461

https://doi.org/10.1107/S160057751601479X

Analytical formula for two-dimensional ring artefact suppression

Valeriy Titarenko ^a ^*

^aHenry Moseley X-ray Imaging Facility, Photon Science Institute, Alan Turing Building, The University of Manchester, Manchester M13 9PL, UK
^*Correspondence e-mail: valeriy.titarenko@manchester.ac.uk

Edited by M. Yamamoto, RIKEN SPring-8 Center, Japan (Received 30 March 2016; accepted 19 September 2016; online 17 October 2016)

Ring artefacts are the most disturbing artefacts when reconstructed volumes are segmented. A lot of effort has already been put into better X-ray optics, scintillators and detectors in order to minimize the appearance of these artefacts. However, additional processing is often required after standard flat-field correction. Several methods exist to suppress artefacts. One group of methods is based on minimization of the Tikhonov functional. An analytical formula for processing of a single sinogram was developed. In this paper a similar approach is used and a formula for processing two-dimensional projections is found. Thus suppression of ring artefacts is organized as a two-dimensional convolution of `averaged' projections with a given filter. Several approaches are discussed in order to find elements of the filter in a faster and accurate way. Examples of experimental datasets processed by the proposed method are considered.

Keywords: X-ray tomography; ring artefacts; ill-posed problem; convolution; filter; Tikhonov functional.

1. Introduction

A standard data acquisition procedure in X-ray tomography for synchrotron- or laboratory-based machines consists of recording three groups of images: dark-field images (when an X-ray beam is switched off), white/flat-field images (the beam is on but there is no sample in the beam) and projections (the beam is on and passes through a rotating sample). Several main scanning geometries are in use these days: parallel for synchrotron radiation facilities and nano-tomography laboratory machines, cone-beam for most medical and industrial applications, and helical/spiral cone-beam with moving and static gantries for medical and security applications (Natterer & Wübbeling, 2001 ; Kak & Slaney, 2001 ; Kalender, 2006 ; Warnett et al., 2016 ). To perform a reconstruction of a specimen one needs to find optical paths by the formula

$[P = \ln \left( {{I_{\rm{flat}} - I_{\rm{dark}}} \over {I - I_{\rm{dark}}}} \right), \eqno(1)]$

where I_flat, I_dark and I are flat-field, dark-field images and a projection, respectively. Note that the flat- and dark-field images should theoretically be taken at the same moment of time when the projection is recorded. Obviously, it is not possible in practical applications, so dark- and flat-field images are usually taken before and after all projections are acquired.

For many beamlines an X-ray beam may slightly move during an experiment. For instance, beam oscillations and linear movement in the vertical direction on beamline 05B1-1 at the Canadian Light Source are discussed by Hinebaugh et al. (2012 ). For tomography beamlines these fluctuations often lead to the flat-field image being changed during data acquisition. This effect can be explained in the following way. X-rays travel in a high-vacuum tube protected by a beryllium window just before entering an optical hutch where experiments are performed. Any small variation of the beam position and thermal instabilities at a monochromator force the intensity profile of the beam to move and stretch/shrink in the vertical direction. As a result, projection of the Be window onto a recording device [usually CCD (charge-coupled device) or CMOS (complementary metal-oxide semiconductor) cameras] behaves in a similar way. On the other hand, a scintillator, optical lenses and detector also have some dust/dirt and defects on their surfaces. All these lead to the fact that a background image has (at least) two components/profiles: one is stable and the other is slightly moving/stretching during data acquisition. Therefore the standard flat-field correction procedure defined by (1) leaves some artefacts on projections which in turn lead to so-called ring artefacts on reconstructed slices, i.e. concentric circles or arches. Some attempts to compensate intensity variations due to beam/monochromator instabilities have been given by Titarenko et al. (2010a ).

Unfortunately, defects on optical components behave differently. If there is an X-ray absorbing feature on a scintillator, then it usually affects those areas on recorded images which are projections of the feature onto the detector. Thus after the standard flat-field correction has been performed, the corresponding ring artefacts seen after reconstruction are regular, i.e. the circles/arches have the same intensity along the whole curve. However, if a scintillator has scratches or dust on its surface (see Fig. 1), then an additional visible light comes from these defects and depends on the total flux of X-rays incident on a wider area around these defects. One can change the position of slits in order to use a smaller shape of the X-ray beam and still be able to see these bright defects in areas with no X-rays (see Fig. 2). As a result, for dense materials or anisotropically attenuated specimens one can see the intensity of ring artefacts dependent on the angle of rotation. Some examples can be found in the paper by Titarenko et al. (2011 ). Of course, if defects on a scintillator are too strong and/or the exposure time to record an image is too long, then some pixels of the detector may become saturated, so we do not have any information from them.

Figure 1
A 620×360 region of a flat-field image acquired at beamline I12 of Diamond Light Source. Random high-intensity pixels and bright structures from defects on a scintillator are seen.

Figure 2
A flat-field image (a) and its dark part with increased contrast (b) acquired at beamline I12 of Diamond Light Source.

Several methods for ring artefact suppression have already been developed. In principle they can be placed into three main groups. The first group of methods uses a `hardware' solution by usually moving a detector with respect to the projections. Thus if there is a defect on a scintillator or Be window then it is either randomly moved to different pixels or blurred (Davis & Elliott, 1997 ; Davis et al., 2013 ). In this case no high-intensity rings can be seen on a reconstructed image. However, if a defect is strong, then one should expect to see a `wave' rather than a ring on slices. The second group is related to preprocessing sinograms before reconstruction is performed by applying various linear or non-linear filters. Methods from the third group process data after reconstruction and usually remap a slice in polar coordinates, so rings are transformed into stripes. Basically methods from these two groups are similar. A detailed description of some of these methods and their comparison for various data sets can be found in the literature (Sijbers & Postnov, 2004 ; Prell et al., 2009 ; Münch et al., 2009 ; Titarenko et al., 2010b ,c ; Sadi et al., 2010 ; Rashid et al., 2012 ; Brun et al., 2013 ; Wei et al., 2013 ; Altunbas et al., 2014 ; Miqueles et al., 2014 ; Paleo & Mirone, 2015 ; Wolkowski et al., 2015 ; Baek et al., 2015 ).

In this paper a new method for preprocessing sinograms/projections is presented. This method is a generalization of a method from Titarenko et al. (2010b), which has already been modified to be used for irregular ring artefacts (Titarenko et al., 2010c, 2011) and is a part of reconstruction software used at beamlines I12 and I13 of Diamond Light Source (Drakopoulos et al., 2015 ; Pešić et al., 2013 ). The original method was designed to process each sinogram independently, i.e. to be applied for one-dimensional (1D) projections, and assumes that `averaged' projections should not have sharp peaks. The same idea of smoothness is used in the proposed algorithm but applied to two-dimensional (2D) projections. Modern reconstruction algorithms usually exploit parallelization provided by graphics processing units (GPUs) and tend to process several slices at the same time. Thus it is more effective to process chunks of projections rather than individual rows. In the case of cone-beam geometry this approach is already used. Therefore, processing 2D projections may speed up other reconstruction steps by avoiding unnecessary movement or remapping of data. At the same time, additional prior information from neighbouring rows of projections should lead to better suppression of artefacts. Note that the solution provided by the original method of Titarenko et al. (2010b) is a vector obtained by multiplying a given matrix by an input vector. In the proposed method a solution (matrix) is a convolution of an input matrix with a 2D filter (matrix), so can be implemented in a more efficient way if a fast Fourier transform is used. The proposed 2D version of the method is a generalization of the 1D method published by the author (Titarenko, 2016 ). In the paper, performance advantages of the filtering approach versus the standard approach based on the inverse matrix from Titarenko et al. (2010b) are shown.

In the next section a brief mathematical description of steps to obtain the main formula is provided. Then examples of implementation of the proposed filter are presented followed by discussion. All proofs and strict mathematical formulas can be found in Appendices A, B and C.

2. Brief mathematical formulation

2.1. Tikhonov functional

Let there be an m ×n matrix A, a known m-vector u and an unknown n-vector z. We want to solve the following system of linear equations,

$[A z = u,\qquad{\rm{or}}\qquad \textstyle\sum\limits_{j\,=\,1}^n A_{ij} z_j = u_i. \eqno(2)]$

Instead of an exact matrix A and the right-hand side u we are given approximate A^h and $[u^\delta]$ , where h and δ are errors. For instance, $[\textstyle\sum_{i\,=\,1}^m\textstyle\sum_{j\,=\,1}^n |A^h_{ij} - A_{ij}|^2]$ $[\leq]$ h² , $[\textstyle\sum_{i\,=\,1}^m |u^\delta_{i} - u_i|^2]$ $[\leq]$ $[\delta^2]$ . A solution of A^h z = $[u^\delta]$ is defined as $[z^{h\delta}]$ . We aim to find an exact solution of (2) knowing A^h, $[u^\delta]$ , h and δ. Thus we have a classical inverse problem. The simplest way to estimate solution z is to find an approximate solution $[z^{h\delta}]$ . One may hope that, if errors h and δ tend to zero, then the corresponding solution $[z^{h\delta}]$ becomes closer to z. Unfortunately, for most inverse problems this is not the case as they are usually so-called ill-posed problems. This means that solutions z or $[z^{h\delta}]$ may not

(i) exist for the pairs (A, u) and/or $[(A^h, u^\delta)]$ ,

(ii) be unique,

(iii) be stable with respect to small variations.

To overcome these issues, Tikhonov proposed an approach based on minimization of the following smoothing functional,

$[M^{\,\alpha}[z] = \textstyle\sum\limits_{i\,=\,1}^m \left(\textstyle\sum\limits_{j\,=\,1}^n A_{ij}^h z_j - u^\delta_i\right)^2\, + \,\,\alpha \textstyle\sum\limits_{j\,=\,1}^n z^2_j, \eqno(3)]$

for a given parameter $[\alpha]$ > 0. Now $[M^{\,\alpha}[z]]$ is often called the Tikhonov functional. It was shown that the corresponding minimizers $[z^\alpha]$ of (3) tend to the exact solution z when errors h and δ approach zero. Several methods of how to choose α as a function of h and δ have been proposed. More details on how to solve ill-posed problems can be found in the literature (Tikhonov & Arsenin, 1977 ; Tikhonov et al., 1995 ; Engl et al., 1996 ).

2.2. Optimization problem

In order to better understand how the theory of inverse and ill-posed problems can be applied to correct ring artefacts, we consider the following continuous problem. Suppose there is a 2D rectangular sensor, so we are able to record an intensity of X-rays incident onto this sensor, i.e. for $[\xi\in[0,L_1]]$ and $[\eta\in[0,L_2]]$ , where ξ and η are the horizontal and vertical coordinates, respectively. As a sample is rotated around its vertical axis we may also record intensities as a function of angle θ (or time t). For a parallel beam we usually assume that θ is from 0 to 180° (or π radians). After the standard flat-field correction procedure (1) we obtain a continuous function $[\tilde p(\xi,\eta,\theta)]$ . If there are no errors in our data acquisition system the flat-field correction procedure should provide us with an exact function $[p(\xi,\eta,\theta)]$ . In the real world a function $[q(\xi,\eta,\theta)]$ = $[p(\xi,\eta,\theta) - \tilde p(\xi,\eta,\theta)]$ may not be zero. Of course, we try to organize data acquisition in such way that possible imperfections of our system do not affect the quality of data. Therefore we may always assume that function $[q(\xi, \eta, \theta)]$ is small with respect to $[p(\xi, \eta, \theta)]$ . So we want the following integral

$[\int_0^{\pi} \int_0^{L_2} \int_0^{L_1} q^2(\xi, \eta, \theta) \,{\rm{d}} \xi \,{\rm{d}}\eta \,{\rm{d}}\theta \eqno(4)]$

to tend to a small value (ideally zero). In the simplest case we may think that the error does not depend on the angle of rotation, i.e. we have a 2D function $[q(\xi, \eta)]$ and for a given recorded value $[\tilde p(\xi, \eta, \theta)]$ we should find an unknown value $[p(\xi,\eta, \theta)]$ such that

$[\int_0^{L_2} \int_0^{L_1} q^2(\xi,\eta,\theta)\,{\rm{d}}\xi\,{\rm{d}}\eta \eqno(5)]$

or equivalently

$[\int_0^{L_2} \int_0^{L_1} \big[\,\tilde p (\xi, \eta, \theta) - p(\xi,\eta,\theta)\big]^2\,{\rm{d}}\xi\,{\rm{d}}\eta \eqno(6)]$

are small. At the same time we have additional, so-called a priori, information that the unknown function $[p(\xi,\eta,\theta)]$ is smooth at least along the ξ and η coordinates (Titarenko et al., 2010d ). Thus we aim to find such a $[p(\xi,\eta,\theta)]$ that

$[\int_0^{L_2} \int_0^{L_1} \left\{\left[{{\partial p} \over {\partial \xi}} (\xi, \eta, \theta)\right]^2 \,+\,\, \left[{{\partial p} \over {\partial \eta}} (\xi, \eta, \theta)\right]^2 \right\}\,{\rm{d}}\xi\,{\rm{d}}\eta \eqno(7)]$

has a minimal value.

2.3. Original 1D problem

In the original method (see Titarenko et al., 2010b), a sinogram $[\tilde p_{k j}]$ is given, $[1\leq k \leq m]$ , $[1\leq j \leq n]$ , where m and n are the number of rows (projections) and columns (pixels in one row). It is assumed that the same unknown value q_j should be added to each element of the jth column in order to find the exact sinogram. For a regularization parameter γ, the Tikhonov functional can be written as

$[M^{\,\gamma} = \textstyle\sum\limits_{k\,=\,1}^m \textstyle\sum\limits_{j\,=\,1}^{n-1} \left(\,\tilde p_{kj} - \tilde p_{k, j+1} + q_j - q_{j+1}\right)^2 \,+\,\, \gamma m \textstyle\sum\limits_{j\,=\,1}^n q_j^2. \eqno(8)]$

Note that the Tikhonov functional can be written in several equivalent forms; for each of them different regularization parameters can be used. In order to separate results obtained for the 1D and 2D cases we use γ for 1D and α for 2D problems. If we average the sinogram along columns, i.e. use p_j = $[({{1}/{m}})\textstyle\sum_{k\,=\,1}^m \tilde p_{k j}]$ , then the above problem leads to minimization of

$[\sum\limits_{j\,=\,1}^{n-1} \left(\,p_j - p_{j+1} + q_j - q_{j+1}\right)^2\, +\,\, \gamma \sum\limits_{j\,=\,1}^n q_j^2, \eqno(9)]$

and an analytical formula to find q_j is presented. Titarenko et al. (2011) used weights w_k to find `averaged' values p_j = $[\textstyle\sum_{k\,=\,1}^m w_k \tilde p_{kj}]$ , so problems to correct for ring artifacts due to non-uniform stripes in sinograms can be solved.

2.4. Mathematical formulation of 2D problem

In the case of 2D projections we can write the Tikhonov functional in a similar way. Suppose we are given an `averaged' projection by adding up projections with some weights. Let this `averaged' projection be a rectangular matrix P with elements p_ij, i = $[1,\ldots,m]$ , j = $[1,\ldots,n]$ . This matrix contains some noisy data and we want to find an unknown exact matrix Z with elements z_ij such that the respective values of z_ij are close to p_ij. This approach is similar to the 1D case if we use q_ij = p_ij - z_ij; however, the choice to use z_ij instead of q_ij is more convenient for the 2D case due to the simpler description for matrices defined below. Our goal is to find matrix Z such that its values z_ij are close to values p_ij of the given averaged projection and use some additional information about Z. We suppose the neighbouring values of Z have similar values; it is a natural assumption (see Titarenko et al., 2010d). Then, for a regularization parameter $[\alpha]$ > 0, the Tikhonov functional can be written as

$[\eqalignno{ M^{\,\alpha} = {}& {{1}\over{2}} {\textstyle\sum\limits_{k\,=\,1}^n \textstyle\sum\limits_{l\,=\,1}^m} \left(\,p_{kl}-z_{kl}\right)^2 \,\,+\,\, {{\alpha}\over{2}} \Bigg\{ \textstyle\sum\limits_{k\,=\,1}^{n-1} \textstyle\sum\limits_{l\,=\,1}^m \left(z_{kl}-z_{k+1l}\right)^2 \cr& + \textstyle\sum\limits_{k\,=\,1}^n \textstyle\sum\limits_{l\,=\,1}^{m-1} \left(z_{kl}-z_{kl+1}\right)^2 \Bigg\}. & (10) }]$

To distinguish regularization parameters used in the original and proposed methods we use different symbols, γ and α. Comparing formulas (8) and (10) one can see $[\alpha]$ ∼ $[1/\gamma]$ . The necessary condition to attain the minimum is that all first derivatives of $[M^{\,\alpha}]$ equal zero, i.e.

$[{{\partial M^{\,\alpha}} \over {\partial z_{kl}}} = \sum\limits_{\hat k\,=\,1}^n \sum\limits_{\hat l\,=\,1}^m A_{kl, \hat k \hat l} z_{\hat k \hat l} - p_{kl} = 0, \eqno(11)]$

where matrix A is defined below.

For convenience, instead of (n×m) matrices P and Z we use nm vectors p and z by stitching rows. If we define i = (k-1)m+l, j = $[(\hat k-1)m +\hat l]$ , N = nm, then

$[{{\partial M^{\,\alpha}} \over {\partial z_i}} = \sum\limits_{j\,=\,1}^N A_{ij} z_j - p_i = 0. \eqno(12)]$

The elements of matrix A can be found in Appendix A. For example, we obtain (12×12) matrix A as

$[\eqalign{&\left( \matrix{ 1+2\alpha & -\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 \cr -\alpha & \hskip-8pt1+3\alpha & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 \cr 0 & \hskip-8pt-\alpha & \hskip-8pt1+3\alpha & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 \cr 0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt1 + 2\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 \cr -\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt1+3\alpha & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 \cr 0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt1+4\alpha & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt 0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 \cr 0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt1+4\alpha & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 \cr 0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt1+3\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha \cr 0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt1+2\alpha & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 \cr 0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt1+3\alpha & \hskip-8pt-\alpha & \hskip-8pt0 \cr 0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt1+3\alpha & \hskip-8pt-\alpha \cr 0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt0 & \hskip-8pt0 & \hskip-8pt-\alpha & \hskip-8pt1+2\alpha } \right)}]$

for n = 3 and m = 4. To find the solution of (12) we want to invert matrix A, so z = A^-1 p. These days images from X-ray tomography usually have a size of about 2000×2000 pixels. So N = 4000000 and it is very difficult to invert (N×N) matrix A even if it is a sparse matrix. It is shown in Appendix A that for n = m > 3 the matrix has (at least) two eigenvalues $[\lambda]$ = 1 and $[1+4\alpha]$ . Therefore the condition number $[\kappa(A)]$ = $[\lambda_{\rm{max}}/\lambda_{\rm{min}}]$ $[\geq]$ 1+ 4α, so for large values α the condition number is large and numerical matrix inversion becomes unstable with respect to small errors in data. Storing the inverse matrix may also be an issue due to its size. Therefore we try the following approximate approach. We find the analytical solution only for the central element of the matrix Z which will be the product of an (nm ×nm) matrix G by the nm vector p. This matrix G may be considered as a 2D filter, which can be used to find matrices Z from given images P when their sizes are greater than the size of G. As shown in the following sections, elements of filter G decay very fast as we go from the central element to the outer ones (several orders of magnitude). Therefore we may set them to zero from some distance or alternatively we may use a smaller filter G. For many practical datasets, (64×64) filters work very well. So for all elements of matrix Z inside a 32-element boundary we may always select a corresponding (64×64) part of matrix P and convolve it with the (64×64) filter G. Or, in an alternative approach, we just keep G and P the same power-of-two size, so a fast Fourier transform can be used.

We intend to find the values of the central element of the matrix Z; therefore we suppose n and m to be odd numbers. Let us define $[\beta]$ = 2 + 1/α, so $[\alpha]$ = $[1/(\beta-2)]$ , and define (m×m) matrix X as

$[X = \left( \matrix{ \beta & -1 & 0 & \cdots & 0 & 0\cr -1 & \beta+1 & -1 & \cdots & 0 & 0\cr 0 & -1 & \beta+1 & \cdots & 0 & 0\cr \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\cr 0 & 0 & 0 & \cdots & \beta+1 & -1\cr 0 & 0 & 0 & \cdots & -1 & \beta } \right), \eqno(14)]$

and denote the (m×m) identity matrix by I. Then the matrix A can be written as the following (n×n) block matrix,

$[A = {{{1}\over{\beta-2}}} \left( \matrix{ X & -I & 0 & \cdots & 0 & 0\cr -I & X+I & -I & \cdots & 0 & 0\cr 0 & -I & X+I & \cdots & 0 & 0\cr \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\cr 0 & 0 & 0 & \cdots & X+I & -I\cr 0 & 0 & 0 & \cdots & -I & X } \right). \eqno(15)]$

Let us introduce a function W_n(x), where the input parameter is a scalar x and the output is the following (n×n) matrix,

$[W_n(x) = \left( \matrix{ x & -1 & 0 & \cdots & 0 & 0\cr -1 & x + 1 & -1 & \cdots & 0 & 0\cr 0 & -1 & x + 1 & \cdots & 0 & 0\cr \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\cr 0 & 0 & 0 & \cdots & x + 1 & -1\cr 0 & 0 & 0 & \cdots & -1 & x } \right). \eqno(16)]$

Note that W_m(x) is an (m×m) matrix. Then matrix A can be written as

$[A = \left({{1}\over{\beta-2}}\right) W_n\left[W_m(\beta)\right], \eqno(17)]$

where $[W_m(\beta)]$ is a matrix and $[W_n[W_m(\beta)]]$ is a matrix function. Therefore the inverse of matrix A is

$[A^{-1} = (\beta-2)W_n^{-1} \left[W_m(\beta)\right], \eqno(18)]$

where W_n^-1(x) is the inverse of matrix W_n(x) (see Higham, 2008 ).

Symmetric matrix W_n(x) can be written in the form Y D Y^T, where D is a diagonal matrix and Y is an orthogonal matrix, i.e. Y^-1 = Y^T, where T stands for the transpose (Higham, 2008). Therefore the inverse matrix W^-1(x) can be written as Y D^-1 Y^T. In the Appendices we find explicit formulas for Y, D and consequently W^-1(x) and $[W^{\,-1}[W_n(\beta)]]$ . As we only need to find the value of the central element of matrix Z, then all matrices (Z, P and filter G) can be rewritten with respect to their central elements, so instead of m×n matrix Z we denote such a matrix as $[\check Z]$ with elements z_ij, i = $[-{({m-1})/{2}},\ldots,{({m-1})/{2}}]$ , j = $[-{({n-1})/{2}},\ldots,{({n-1})/{2}}]$ . Note that n and m are odd numbers. If these numbers are large, then the elements of the filter (matrix) $[\check G]$ can be written as

$[\check G_{jk} = {{\beta - 2} \over {4\pi^2}} \int_{-\pi}^\pi \int_{-\pi}^\pi {{\cos j\xi \cos k \eta} \over { \beta + 2\left(1 - \cos \xi - \cos \eta \right)}} \,{\rm{d}}\eta\,{\rm{d}}\xi]$

[see formula (59) in Appendix B]. If C^k_n $[\equiv]$ n!/[k!(n-k)!] is a binomial coefficient, then the elements of the filter $[\check G]$ can be written as the infinite sum

$[\check G_{jk} = (1-4\tau)\tau^{\,k+j} \textstyle\sum\limits_{q\,=\,0}^{\infty} C_{2q + j + k}^{\,q} \, C_{2q + j + k}^{\,q+j} \, \tau^{\,2 q}, \eqno(19)]$

where $[\tau]$ = $[\alpha/(1+4\alpha)]$ . The above formula [or formula (61)] is derived in Appendix B. In Appendix C it is shown how to find the infinite sum numerically with high precision. Some useful properties of matrix $[\check G]$ are derived in Appendix B. Examples of matrices $[\check G]$ are shown in Fig. 3. Matrix $[\check G]$ is symmetrical with respect to central horizontal, vertical and diagonal lines. All elements of the matrix are positive and for each row they monotonically decrease with the column index (see Fig. 4). The sum of all elements of matrix $[\check G]$ (of infinite width/height) is 1.

Figure 3
The central 51×51 section of matrix $[\check G]$ for (a) $[\alpha]$ = 1000, (b) $[\alpha]$ = 10000.

Figure 4
Values of matrix $[\check G]$ found for $[\alpha]$ = 300 and rows 0, 5, 10, 15, 20, 25 (blue, red, green, orange, black, violet, respectively).

Elements of matrix $[\check G]$ seem to be radial symmetric; however, there is a small error between the real values and those found from elements of the central row. If (i,j) is the index of an element of matrix $[\check G]$ and i,j > 0, then its distance from the central element is r = (i²+j²)^1/2, so we may find an index k such that k $[\leq]$ r $[\lt]$ k+1 and use elements $[\check G_{k-1,0}]$ , $[\check G_{k0}]$ , $[\check G_{k+1,0}]$ and $[\check G_{k+2,0}]$ to approximate value $[\check G_{ij}]$ by cubic interpolation. Relative errors of such interpolation with respect to the exact element $[\check G_{ij}]$ are shown in Fig. 5.

Figure 5
Relative error of cubic interpolation for (a) $[\alpha]$ = 100, (b) $[\alpha]$ = 2000.

If $[\alpha]$ = 0 then filtering should not change the input matrix, so all elements of $[\check G]$ are zeros except the central element, which is 1. The central element $[\check G_{00}]$ decreases when α increases; the other elements equal zero for $[\alpha]$ = 0, increase, have one maximum and decrease as in Fig. 6.

Figure 6
Values of elements (0,0), (0,1) and (0,2) of matrix $[\check G]$ as a function of α.

3. Example problems

The method proposed in this paper is applied to several datasets acquired at station 16.3 of the Synchrotron Radiation Source (SRS), Daresbury Laboratory (second-generation synchrotron source, now decommissioned) and beamline I12 of Diamond Light Source (DLS; third-generation synchrotron source) in 2008/2009. The choice to use these datasets is dictated by imperfections of an experimental setup in those days as these datasets were one of the first test datasets. For instance, in some cases the rotation axis was not perfectly vertical; a lot of high-intensity peaks were also present on projections as those shown in Fig. 1. A lot of work has been done at I12 to obtain images of sufficiently better quality (Drakopoulos et al., 2015). However, these datasets are still useful as similar issues can be seen at other beamlines worldwide. For all datasets a PCO 4000 14-bit CCD camera, monochromatic beam and CdWO₄ scintillator were used; reconstruction was performed using software developed by the author for beamline I12.

3.1. Graphite

A 5.0 mm × 2.5 mm piece of graphite with some metal particles was scanned at SRS (see Figs. 7 and 8). This specimen had produced almost ideal regular ring artefacts, even in the presence of these dense metal particles. The graphite was slightly tilted and as a result the outer regions of the sample are slightly blurred while the central part is sharp. The ring suppression method is applied for regularization parameter $[\alpha]$ = 1000.

Figure 7
A piece of graphite: (a) the whole section, (b) the central part. No artefact suppression.

Figure 8
Central part of the graphite sample after the ring suppression algorithm is applied, $[\alpha]$ = 1000.

3.2. TRISO particles

A cylindrical container with 1 mm tristructural-isotropic (TRISO) test particles was also scanned at SRS. A similar specimen was used by Titarenko et al. (2009 ) in order to show how a priori information (a known structure of these particles) can be used to remove rings. However, this prior information was applied in a manual way when corresponding areas of the same materials were selected. In this paper an automatic procedure is used.

Beamlines produce almost parallel X-rays, so reconstruction of a sample can be parallelized as all horizontal slices can be found independently of other slices. Therefore many ring artefact suppression algorithms can be used as they usually process 2D slices. However, if one has a cone-beam geometry or a sample is tilted then several neighbouring sinograms should be used to find one slice. In Fig. 9 the specimen is tilted (0.1° along the beam and 0.8° perpendicular to the beam). Therefore it is more efficient to process projections first and then reconstruct slices with a modified filtered backprojection algorithm for tilted samples.

Figure 9
Cylindrical container with TRISO particles, tilted axis of rotation. (a), (b) No ring artefact suppression or tilt correction have been applied, (c) only ring artefacts have been suppressed for $[\alpha]$ = 10, (d) both corrections have been performed.

3.3. Mouse skull

A mouse skull was scanned at beamline I12, DLS. Due to scattering, a lot of high-intensity photons reached the detector. Identification of the corresponding pixels is easier to perform on neighbouring projections as three-dimensional neighbouring pixels are used rather than 2D in the case of individual sinograms. A simple correction, e.g. based on linear interpolation of unaffected pixels, can be used (see Fig. 10). However, due to scratches on the scintillator (similar to those shown in Fig. 2), ring artefact suppression still leaves some rings present after correction.

Figure 10
Part of a mouse skull, with a lot of high-intensity pixels on the projections. (a) No correction for high-intensity pixels or ring artefacts, (b) high-intensity pixels are suppressed, (c) both corrections have been performed, $[\alpha]$ = 10.

4. Practical implementation

The main formula (19) is a series, i.e. a sum of infinite terms. Effective numerical implementation is discussed in Appendix C. A practical approach is to find the central row of matrix $[\check G]$ by the formula (64). Then the row above the central one can be found from

$[\check G_{1k} = {{1} \over {2}} \left({{\check G_{0k}} \over {\tau}} - \check G_{0 k -1} - \check G_{0 k+1}\right) \eqno(20)]$

for k > 0, and elements of other rows

$[\check G_{jk} = {{\check G_{j-1,k}} \over {\tau}} - \check G_{j-1, k -1} - \check G_{j-1, k+1} - \check G_{j-2, k} \eqno(21)]$

for k $[\geq]$ j. Suppose we want to find matrix $[\check G_{jk}]$ for [|j|,|k|] $[\leq]$ N. Then we find the values as shown in Fig. 11, i.e.

Figure 11
Use of the central row of matrix $[\check G]$ to find all other elements of the matrix.

(i) For k = $[0,\ldots,2N]$ , we use (64) to find the central row.

(ii) For k = $[1, \ldots, 2N-1]$ , we use (20).

(iii) For j = $[2, \ldots, N]$ and k = $[j,\ldots, 2N-j]$ , we apply (21).

(iv) For j = $[1,\ldots, N]$ , k = $[0,\ldots,j]$ : $[\check G_{j,-k}]$ = $[\check G_{-j, k}]$ = $[\check G_{-j, -k}]$ = $[\check G_{k, j}]$ = $[\check G_{k, -j}]$ = $[\check G_{-k, j}]$ = $[\check G_{-k, -j}]$ = $[\check G_{j, k}]$ .

5. Discussion

The 2D method presented in this paper is in some sense equivalent to the original 1D method presented in §2.3 or by Titarenko et al. (2010b). Thus they both should provide similar results in the case of regular ring artefacts when the rings' strength does not depend on the angle of rotation. For the data sets used as examples the author was able to find values of regularization parameter γ for the original 1D method which gives almost identical images as the proposed method. Of course, in the case of irregular rings the result of suppression may depend on a sample. The 1D method designed for irregular rings (Titarenko et al., 2011) can also be used with the 2D method as the original minimization problem leads to a set of independent problems, which can be solved with the new method.

The main advantage of the 2D method is that it can be considered as an image filtering procedure. This filtering operation is applied to an `averaged' optical path image and can be used after all projections have been collected. However, this `averaged' image is just a weighted sum of individual optical path images. Therefore this filtering operation can be applied to individual images just after they have been acquired. Taking into account that standard (non-iterative) reconstruction algorithms apply a 1D filter to each projection, then these 1D and 2D filters can be combined in one 2D filter, i.e. ring suppression can be embedded into a standard reconstruction procedure but with a modified filter.

Implementation of 2D filtering can be performed with fast Fourier transform operations, and high-performance libraries exist for both CPUs and GPUs.

APPENDIX A

A1. Two eigenvalues for matrix A

The matrix A is symmetric and all its non-zero elements are on five diagonals: A_ii, A_{i, i-1}, A_{i, i+1}, A_{i, i-m}, A_{i, i+m}. Consider each row of matrix A.

(1) Three non-zero elements:

(a) i = 1, $[A_{ii} = 1+2\alpha]$ , $[A_{i, i+1} = -\alpha]$ , $[A_{i, i+m} = -\alpha]$ ;

(b) i = n, $[A_{i, i-1} = -\alpha]$ , $[A_{ii} = 1+2\alpha]$ , $[A_{i, i+m} = -\alpha]$ ;

(d) i = n m, $[A_{i, i-m} = -\alpha]$ , $[A_{i, i-1} = -\alpha]$ , $[A_{ii} = 1+2\alpha]$ .

(2) Four non-zero elements:

(a) $[1\,\, \lt\,\, i \,\,\lt \,\,n]$ , $[A_{i, i-1} = -\alpha]$ , $[A_{ii} = 1+3\alpha]$ , $[A_{i, i+1} = -\alpha]$ , $[A_{i, i+m} = -\alpha]$ ;

(b) $[i = n + 1, n + m +1,\ldots, m (n - 2) + 1]$ , $[A_{i, i-m} = -\alpha]$ , $[A_{ii} = 1+3\alpha]$ , $[A_{i, i+1} = -\alpha]$ , $[A_{i, i+m} = -\alpha]$ ;

(c) $[i = n + m, n + 2 m,\ldots, m (n - 1)]$ , $[A_{i, i-m} = -\alpha]$ , A_{i, i-1} = $[-\alpha]$ , A_ii = $[1+3\alpha]$ , $[A_{i, i+m} = -\alpha]$ ;

(d) $[n (m - 1) + 1 \,\,\lt\,\, i \,\,\lt\,\, n m]$ , $[A_{i, i-m} = -\alpha]$ , $[A_{i, i-1} = -\alpha]$ , $[A_{ii} = 1+3\alpha]$ , $[A_{i, i+1} = -\alpha]$ .

(3) Five non-zero elements (the rest): $[A_{i, i-m} = -\alpha]$ , $[A_{i, i-1} = -\alpha]$ , $[A_{ii} = 1+4\alpha]$ , $[A_{i, i+1} = -\alpha]$ , $[A_{i, i+m} = -\alpha]$ .

It is clear that the sum of all elements of matrix A for each row is 1. So if (nm) vector $[v = \{1, 1, \ldots, 1\}]$ , then (A- I)v = 0 and thus $[\lambda]$ = 1 is the corresponding eigenvalue.

Now we consider case n = m $[\gt]$ 3 and check that $[\lambda]$ = 1 + 4α is another eigenvalue, and non-zero elements of one of the corresponding eigenvectors can be written as v_i(n-1) = v_{(i+1)(n-1) + 2} = (-1)ⁱ, i = $[1,\ldots, n-1]$ and v_n = 1, v_{n(n-1) + 1} = (-1)ⁿ. For convenience this vector can be represented as an (n×n) matrix by copying each n consequent elements to a row of the matrix, then it is an anti three-diagonal matrix, e.g. for n = m = 5,

$[v = \left(\matrix{ 0 & 0 & 0 & -1 & 1\cr 0 & 0 & 1 & 0 & -1\cr 0 & -1 & 0 & 1 & 0\cr 1 & 0 & -1 & 0 & 0\cr -1 & 1 & 0 & 0 & 0 } \right). \eqno(22)]$

Denote $[Q = (A - \lambda I)/\alpha]$ , then

(1) Three non-zero elements:

(a) i = 1, Q_ii = -2, Q_{i, i+1} = -1, Q_{i, i+n} = -1 and v_i = 0, v_i+1 = 0, v_{i +n} = 0;

(b) i = n, Q_{i, i-1} = -1, Q_ii = -2, Q_{i, i+n} = -1 and v_i-1 = -1, v_i = 1, v_i+n = -1;

(d) i = n n, Q_{i, i-n} = -1, Q_{i, i-1} = -1, Q_ii = -2 and v_i-n = 0, v_i-1 = 0, v_i = 0.

(2) Four non-zero elements:

(a) $[1\,\, \lt\,\, i\,\, \lt\,\, n]$ , Q_{i, i-1} = Q_ii = Q_{i, i+1} = Q_{i, i+n} = -1; v_i-1 = v_i = 0, v_i+1 = -1, v_i+n = 1 (for i = n-2), v_i-1 = 0, v_i = -1, v_i+1 = 1, v_i+n = 0 (for i = n-1) and v_i-1 = v_i = v_i+1 = v_i+n = 0 (otherwise);

(b) $[i = n + 1, 2n + 1,\ldots, n (n - 2) + 1]$ , Q_{i, i-n} = Q_ii = Q_{i, i+1} = Q_{i, i+n} = -1; v_i-n = v_i = 0, v_i+1 = (-1)ⁿ, v_i+n = -(-1)ⁿ [for i = n(n-3)+1], v_i-n = 0, v_i = -(-1)ⁿ, v_i+1 = 0, v_i+n = (-1)ⁿ [for i = n(n-2)+1], v_i-n = v_i = v_i+1 = v_i+n = 0 (otherwise);

(c) $[i = 2n, 3n,\ldots, n (n - 1)]$ , Q_{i, i-n} = Q_{i, i-1} = Q_ii = Q_{i, i+n} = -1; v_i-n = 1, v_i-1 = 0, v_i = -1, v_i+n = 0 (for i = 2n), v_i-n = -1, v_i-1 = 1, v_i = v_i+n = 0 (for i = 3n), v_i-n = v_i-1 = v_i = v_i+n = 0 (otherwise);

(d) $[n (n - 1) + 1\,\, \lt\,\, i\,\, \lt\,\, n n]$ , Q_{i, i-n} = -1, Q_{i, i-1} = Q_ii = Q_{i, i+1} = -1; v_i-n = 0, v_i-1 = (-1)ⁿ, v_i = -(-1)ⁿ, v_i+1 = 0 [for i = n (n-1)+2], v_i-n = (-1)ⁿ, v_i-1 = -(-1)ⁿ, v_i = v_i+1 = 0 [for i = n(n-1) +3], v_i-n = v_i-1 = v_i = v_i+1 = 0 (otherwise);

(e) the rest: Q_{i, i-n} = Q_{i, i-1} = Q_{i, i+1} = Q_{i, i+n} = -1; v_i-n = (-1)^j, v_i-1 = -(-1)^j, v_i+1 = (-1)^j, v_i+n = -(-1)^j (for i = jn + n - j), v_i-n = v_i-1 = 0, v_i+1 = -(-1)^j, v_i+n = (-1)^j (for i = j n - j -1), v_i-n = (-1)^j, v_i-1 = -(-1)^j, v_i+1 = v_i+n = 0 (for i = j n +2n - j +1), v_i-n = v_i-1 = v_i+1 = v_i+n = 0 (otherwise).

Now it is easy to check that $[\textstyle\sum_{j\,=\,1}^{n^2} Q_{ij} v_j]$ = 0.

APPENDIX B

B1. Formula for 2D filter

B1.1. Trigonometric formulas

From the formulas for the sum of sines/cosines we find

$[\sin(n\omega + \psi) + \sin[(n-2)\omega + \psi] = 2 \sin[(n-1)\omega + \psi]\cos\omega, \eqno(23)]$

$[\cos (n+1)x + \cos (n-1) x = 2 \cos n x \cos x, \eqno(24)]$

$[\eqalignno{&\cos(n&-1)x \cos my + \cos(n+1)x\cos my + \cos nx \cos (m-1)y \cr& + \cos nx \cos (m+1)y = 2(\cos x+\cos y)\cos n x \cos m y. & (25)}]$

The following formula is true (Gradshteyn & Ryzhik, 2015 ),

$[\sum\limits_{k\,=\,0}^{n-1} \cos(a + k b) = \cos\left\{a+\textstyle{[({n-1})/{2}]}b\right\} {{\sin (n b/2)} \over {\sin (b/2)}}. \eqno(26)]$

In the case of a = 0, equation (26) gives us

$[\sum_{k\,=\,0}^{n} \cos(k b) = {{1} \over {2}}\left\{1 + {{\sin[(2n + 1)b/2]} \over {\sin (b/2)}}\right\}.]$

Thus we find the Dirichlet sum

$[1 + 2 \sum_{k\,=\,1}^n \cos (k b) = {{\sin [(2 n + 1)b/2]} \over {\sin (b/2)}}. \eqno(27)]$

Some limits for trigonometric functions:

$[\lim\limits_{x \to 0}{{\sin n x} \over {\sin x}} = n, \eqno(28)]$

$[\lim\limits_{x\to \pi}{{\sin n x} \over {\sin x}} = (-1)^{n+1} n. \eqno(29)]$

B1.2. Combinatorial formulas

A factorial n! of a positive integer number n is defined as n! = $[n\cdot (n-1)\cdots 2 \cdot 1]$ , 0! $[\equiv]$ 1, a binomial coefficient C_n^k can be defined as

$[C_n^{\,k}\equiv {{n!} \over {k! (n-k)!}}.]$

From the definition C_n^k = C_n^n-k, the following formulas are true (Benjamin & Quinn, 2003 ):

$[{{1} \over {1-x}} = \sum_{n\,=\,0}^{\infty} x^n\qquad{\rm{for}}\,\,|x| \,\lt\, 1, \eqno(30)]$

$[(x + y)^n = \textstyle\sum\limits_{k\,=\,0}^n C_n^{\,k} x^k y^{n-k}, \eqno(31)]$

$[\textstyle\sum\limits_{j\,=\,0}^k C_m^{\,j} C_n^{\,k-j} = C_{m+n}^{\,k} \eqno(32)]$

(Vandermonde convolution/identity).

B1.3. Simple integrals

If n and m are non-negative integer numbers, then for $[n\geq m]$ ,

$[\int_0^{\pi/2} \sin^{2n +1} x \sin (2m+1)x \,{\rm{d}} x = {{(-1)^m \pi} \over {2^{2 n +2}}} \, C_{2n + 1}^{\,n-m}, \eqno(33)]$

$[\int_0^{\pi/2} \sin^{2n} x \cos 2m x\,{\rm{d}}x = {{(-1)^m \pi} \over {2^{2 n +1}}} \, C_{2n}^{\,n-m}, \eqno(34)]$

and the above integrals equal zero for n < m (see Gradshteyn & Ryzhik, 2015).

Let f(x) and h(x) be odd and even functions, respectively, i.e. f(-x) = -f(x) and h(-x) = h(x), then

$[\int_{-a}^{a} f(x)\,{\rm{d}}x = 0, \eqno(35)]$

$[\int_{-a}^{a} h(x)\,{\rm{d}}x = 2 \int_{0}^{a} h(x)\,{\rm{d}}x. \eqno(36)]$

Taking into account that $[\cos[x+({{\pi}/{2}})] = -\sin x]$ , we find for non-negative integer p,

$[\int_0^{\pi} \cos^{\,p} y\,{\rm{d}}y = (-1)^{\,p}\int_{-\pi/2}^{\pi/2} \sin^{\,p} x\,{\rm{d}}x.]$

If p is an odd number, then from (35)

$[\int_0^\pi \cos^{2n+1} x\,{\rm{d}}x = 0. \eqno(37)]$

In the case of an even p, equation (34) can be used and

$[\int_0^{\pi} \cos^{2n} x\,{\rm{d}}x = {{\pi} \over {2^{2 n}}} C_{2n}^{\,n}. \eqno(38)]$

In a similar way for $[n\geq m]$ we obtain

$[\eqalignno{\int_0^\pi \cos^{2n} y \cos 2m y\,{\rm{d}}y & = 2(-1)^m \int_0^{\pi/2} \sin^{2n} x \cos 2m x \,{\rm{d}}x \cr& = {{\pi} \over {2^{2n}}} \, C_{2n}^{\,n-m}, &(39)}]$

$[\int_0^\pi \cos^{2n+1} y \cos 2m y\,{\rm{d}}y = 0, \eqno(40)]$

$[\int_0^\pi \cos^{2n} y \cos (2m +1)y\,{\rm{d}}y = 0, \eqno(41)]$

$[\eqalignno{ \int_0^\pi \cos^{2n+1}& y \cos (2m +1)y\,{\rm{d}}y \cr& = 2(-1)^m\int_0^{\pi/2} \sin^{2n+1} x \sin (2m+1) x\,{\rm{d}}x \cr& = {{\pi} \over {2^{2n+1}}} C_{2n+1}^{\,n-m}. &(42)}]$

The above integrals are zeros in the case of n < m.

The Dirac δ-function can be written as the following sum,

$[\delta(x) = {{1} \over {\pi}}\,\lim_{n\to\infty} {{\sin nx} \over {x}}. \eqno(43)]$

B1.4. The first sequence

Suppose there is a number ρ, $[|\rho|]$ < 1/4, and for integer $[n\geq 0]$ we want to find

$[T_n = \int_0^\pi\int_0^\pi {{\cos nx} \over {1-2\rho (\cos x + \cos y)}}\,{\rm{d}}x\,{\rm{d}}y. \eqno(44)]$

As $[\cos x \leq 1]$ , $[\cos y \leq 1]$ , then $[1-2\rho (\cos x + \cos y)]$ > 0, so T_n is not a singular integral. From (30) we obtain

$[T_n = \sum_{k\,=\,0}^{\infty} (2\rho)^k U_{nk}, \eqno(45)]$

$[U_{nk}\equiv \int_0^\pi\int_0^\pi \cos nx (\cos x + \cos y)^k\,{\rm{d}}x\,{\rm{d}}y.]$

Taking into account the binomial theorem (31) we obtain

$[U_{n k} = \sum_{m\,=\,0}^k C_k^{\,q} \int_0^{\pi}\cos^m x \cos nx \,{\rm{d}} x \int_0^{\pi}\cos^{k-m} y\,{\rm{d}}y.]$

From (37) and (38) it follows that, for even k = $[2\tilde k]$ , we should consider only even m = 2s,

$[U_{n k} = \sum_{s\,=\,0}^{\tilde k} \,\,{{\pi} \over {2^{2(\tilde k - s)}}} \,\,C_{2\tilde k}^{\,2 s} C_{2(\tilde k - s)}^{\,\tilde k - s} \int_0^{\pi}\cos^{2 s} x \cos nx\,{\rm{d}}x]$

and for odd k = $[2\tilde k + 1]$ , we should consider only odd m = 2s + 1,

$[U_{n k} = \sum_{s\,=\,0}^{\tilde k} \,\,{{\pi} \over {2^{2(\tilde k - s)}}} \,\,C_{2\tilde k + 1}^{\,2 s+1} C_{2(\tilde k - s)}^{\,\tilde k - s} \int_0^{\pi}\cos^{2 s+1} x \cos nx\,{\rm{d}}x.]$

Due to (38)–(42) we see that U_kn is not zero only when both k and n are either even or odd numbers and $[2s\geq n]$ or $[2s + 1 \geq n]$ , respectively. So for even k and n, i.e. k = $[2\tilde k]$ , n = $[2\tilde n]$ , we find

$[\eqalignno{ U_{n k} & = {{\pi^2} \over {2^k}} \sum\limits_{s\,=\,\tilde n}^{\tilde k} C_{2\tilde k}^{\,2 s} C_{2(\tilde k - s)}^{\,\tilde k - s} C_{2s}^{\,s-\tilde n} = {{\pi^2} \over {2^k}} C_{2\tilde k}^{\,\tilde k -\tilde n} \sum_{s\,=\,\tilde n}^{\tilde k} C_{\tilde k -\tilde n}^{\,\tilde k - s} C_{\tilde k +\tilde n}^{\,\tilde k - s} \cr & = {{\pi^2} \over {2^k}} C_{2\tilde k}^{\,\tilde k -\tilde n} \sum_{j\,=\,0}^{\tilde k- \tilde n} C_{\tilde k -\tilde n}^{\,j} C_{\tilde k +\tilde n}^{\,\tilde k - \tilde n - j} = {{\pi^2} \over {2^k}} \left(C_{2\tilde k}^{\,\tilde k -\tilde n}\right)^2, & (46)} ]$

where the Vandermode identity (32) is used, and for odd k and n, i.e. k = $[2\tilde k + 1]$ , n = $[2\tilde n + 1]$ ,

$[U_{n k} = {{\pi^2} \over {2^k}} \sum_{s\,=\,\tilde n}^{\tilde k} C_{2\tilde k + 1}^{\,2 s+1} C_{2(\tilde k - s)}^{\,\tilde k - s} C_{2s+1}^{\,s-\tilde n} = {{\pi^2} \over {2^k}} \left(C_{2\tilde k+1}^{\,\tilde k -\tilde n}\right)^2. \eqno(47)]$

So instead of (46) and (47) we obtain one formula when n and k are both even or odd,

$[U_{nk} = {{\pi^2} \over {2^k}} \, \left[ C_k^{\,(k-n)/2} \right]^2.]$

U_{n k} = 0 for n < k or odd (k-n). We replace k = n+2q in (45),

$[T_n = \pi^2 \rho^n \sum\limits_{q\,=\,0}^{\infty} \left[ \rho^q C_{n+2q}^{\,q} \right]^2.]$

B1.5. The second sequence

The integral (44) depends only on one integer parameter n. Now we consider

$[S_{nm} = \int_0^\pi\int_0^\pi {{\cos nx\cos my} \over {1-2\rho (\cos x + \cos y)}} \,{\rm{d}} x \,{\rm{d}}y, \eqno(48)]$

which depends on two integer parameters n and m. Due to properties of the cosine function and symmetry of S_nm it is enough to consider only the case of $[0\leq m \leq n]$ as S_nm = S_mn = S_{-n m} = S_{n, -m}. Note that

$[{{(\cos x + \cos y)} \over {1-2\rho(\cos x + \cos y)}} = {{1} \over {2\rho}}\left[{{1} \over {1-2\rho(\cos x + \cos y)}} - 1\right].]$

Let n, $[m\ne 0]$ at the same time, then

$[\int_0^\pi\int_0^\pi \cos nx\cos my\,{\rm{d}}x\,{\rm{d}}y = 0,]$

$[\int_0^\pi\int_0^\pi {{(\cos x + \cos y)\cos nx\cos my} \over {1-2\rho (\cos x + \cos y)}}\,{\rm{d}}x\,{\rm{d}}y = {{1} \over {2\rho}} \, S_{nm}.]$

Therefore if we use identity (25) then

$[S_{n-1, m} + S_{n+1,m} + S_{n, m-1} + S_{n, m+1} = {{S_{nm}}\over{\rho}}. \eqno(49)]$

For n = m = 0 we obtain

$[S_{10} = {{1} \over {4}} \left(S_{-1 0} + S_{1 0} + S_{0, -1} + S_{0 1} \right) = {{S_{00} - \pi^2} \over {4\rho}}. \eqno(50)]$

Taking into account (49) and S_{n, -1} = S_{n, 1}, we obtain

$[2S_{n1} = {{S_{n,0}} \over {\rho}} - S_{n-1, 0} - S_{n+1,0} = {{T_n} \over {\rho}} - T_{n-1} - T_{n+1}.]$

As (C_n+2q^q)² - (C_n-1+2q^q)² = 0 for q = 0, we may write

$[\eqalign{{{T_n} \over {\rho}} - T_{n-1} & = \pi^2\rho^{n-1}\sum_{q\,=\,0}^\infty\rho^{2q}\left[(C_{n+2q}^{\,q})^2 - (C_{n-1+2q}^{\,q})^2\right]\cr {} & = \pi^2\rho^{n+1}\sum_{s\,=\,0}^\infty\rho^{2s}\left[(C_{n+2s+2}^{\,s+1})^2 - (C_{n+2s+1}^{\,s+1})^2\right],}]$

$[\eqalign{{{S_{n1}} \over {\pi^2 \rho^{n+1}}} & = \sum_{q\,=\,0}^{\infty}{{\rho^{2q}} \over {2}}\left[(C_{n+2q+2}^{\,q+1})^2 - (C_{n+2q+1}^{\,q+1})^2 - (C_{n+2q +1}^{\,q})^2\right]\cr {} & = \sum_{q\,=\,0}^{\infty}\rho^{2q} C_{n+2q+1}^{\,q} C_{n+2q+1}^{\,q+1}.}]$

Now check that the following formula is true for $[0\leq m \leq n]$ ,

$[S_{nm} = \pi^2 \rho^{n+m} \sum_{q\,=\,0}^{\infty} \rho^{2q} C_{2q+n + m}^{\,q}\, C_{2q+n+m}^{\,q+m}. \eqno(51)]$

We use mathematical induction. Above we have checked the formula for m = 0 and m = 1. Assume the formula is true for $[0\leq m ]$ < k. Then

$[{{S_{n-1, k}/\rho - S_{n-1, k-1} - S_{n-2,k}} \over {\pi^2 \rho^{n+k-2}}}]$

can be rewritten as a Taylor series of $[\rho^2]$ with coefficients

$[\eqalign{C_{2q+n + k - 1}^{\,q} & C_{2q+n+k-1}^{\,q+k} - C_{2q+n + k-2}^{\,q} C_{2q+n+k-2}^{\,q+k-1} \cr& - C_{2q+n + k-2}^{\,q} C_{2q+n+k-2}^{\,q+k} = C_{2q + n + k -2}^{\,q-1} C_{2q + n + k -1}^{\,q+k}}]$

for q > 0 and zero for q = 0. Therefore

$[\eqalign{{{S_{nk}} \over {\pi^2 \rho^{n+k}}} & = {{S_{n-1, k}/\rho - S_{n-1, k-1} - S_{n-2,k} - S_{n-1,k+1}} \over {\pi^2 \rho^{n+k-2}}} \cr& = \sum_{s\,=\,0}^{\infty} \rho^{2s}\left[C_{2s + n + k}^{\,s} C_{2s + n + k +1}^{\,s+k + 1} - C_{2s+n + k}^{\,s} C_{2s+n+k}^{\,s+k+1}\right] \cr& = \sum_{s\,=\,0}^{\infty} \rho^{2s} C_{2s + n + k}^{\,s} C_{2s + n + k}^{\,s+k}. }]$

B1.6. The third sequence

Let there be $[\varphi \in(0,{{\pi}/{2}})]$ and

$[V_n \equiv \int_0^{\pi} {{\cos n x} \over {1- \cos x \sin\varphi}}\,{\rm{d}}x.]$

Using formula (24) we obtain, for n > 0,

$[\eqalign{ V_{n+1} + V_{n-1} &= \int_0^{\pi} {{2\cos x \cos n x} \over {1- \cos x \sin\varphi}}\,{\rm{d}}x\cr {} & = {{2} \over {\sin\varphi}}\int_0^{\pi} \left({{1} \over {1- \cos x \sin\varphi}} - 1\right) \cos n x\,{\rm{d}}x\cr {} & = {{2} \over {\sin\varphi}}\, V_n,}]$

since $[\textstyle\int_0^{\pi} \cos n x\,{\rm{d}}x]$ = 0. We may find that

$[V_0 = {{\pi} \over {\cos \varphi}}, \qquad V_1 = {{\pi} \over {\cos\varphi}}\,\,{{\sin\varphi} \over {1+ \cos\varphi}}]$

and therefore use the recurrent formula

$[V_n = {{2} \over {\sin\varphi}} \,V_{n-1} - V_{n-2}.]$

Since

$[{{2} \over {\sin\varphi}}\,\, {{\sin\varphi} \over {1+ \cos\varphi}} - 1 = {{1-\cos\varphi} \over {1+\cos\varphi}} = {{\sin^2\varphi} \over {(1+\cos\varphi)^2}},]$

then V_n is a geometric sequence, i.e. V_n+1 = $[\gamma V_n]$ , where

$[\gamma = {{\sin\varphi} \over {1+ \cos\varphi}}.]$

So the values of sequence V_n can be written as

$[V_n = {{\pi} \over {\cos \varphi}} \left({{\sin\varphi} \over {1+ \cos\varphi}}\right)^{n}.]$

B1.7. Eigenvalues/eigenvectors

Let us consider an (n×n) matrix B_n,

$[B_n = \left( \matrix{ 2 \sigma - 1 & -1 & 0 & \cdots & 0 & 0\cr -1 & 2 \sigma & -1 & \cdots & 0 & 0\cr 0 & -1 & 2\sigma & \cdots & 0 & 0\cr \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\cr 0 & 0 & 0 & \cdots & 2\sigma & -1\cr 0 & 0 & 0 & \cdots & -1 & 2\sigma - 1 }\right). \eqno(52)]$

We aim to find its eigenvalues, therefore we want to solve equation det B_n = 0. For this purpose we replace the element (1,1) of B_n by $[2\sigma]$ and denote the corresponding matrix as F_n. If we use the first row of B_n to find its determinant, then det B_n = $[(2\sigma - 1) \det F_{n-1} - \det F_{n-2}]$ . In a similar way we obtain det F_n = $[2\sigma \det F_{n-1} - \det F_{n-2}]$ , thus det B_n = det F_n - det F_n-1, det F_n + det F_n-2 = $[2\sigma \det F_{n-1}]$ and

$[\det B_n - \det B_{n-2} = 2\sigma \det B_{n-1}.]$

Therefore we need to find a sequence {b_n} such that

$[b_{n} + b_{n-2} = 2\sigma b_{n-1}. \eqno(53)]$

Of course if one sequence {b_n} is found, then $[\{\tau b_n\}]$ is also a solution of the above equation. We want to find τ such that det B_n = b_n.

Let $[\sigma \in[-1,1]]$ , then taking into account (23) a solution of (53) is b_n = $[\tau \sin(n\omega + \psi)]$ , $[\sigma]$ = $[\cos\omega]$ . Note that matrices B_n are formally defined for $[n\geq 3]$ . One can find their determinants explicitly for n = 3 and 4: det B₃ = $[8 \sigma^3 - 8 \sigma^2 - 2 \sigma + 2]$ , det B₄ = $[16 \sigma^4 - 16 \sigma^3 - 8 \sigma^2 + 8\sigma]$ , so from (53) we obtain b_n = $[2\sigma b_{n+1} - b_{n+2}]$ and can extrapolate values b_n for n < 3: b₂ = $[4\sigma^2 - 4\sigma]$ , b₁ = $[2\sigma - 2]$ , b₀ = 0. Thus we obtain the following system of equations,

$[\eqalign{ b_0 & = \tau\sin(\psi) = 0,\cr b_1 & = \tau\sin(\omega + \psi) = 2\sigma - 2. }]$

The first equation is valid if either $[\tau]$ = 0 or $[\psi]$ = 0 or π. The case of $[\tau]$ = 0 is not possible, since the second equation would always be zero. Since $[\sin (\theta + \pi)]$ = $[- \sin\theta]$ , then we may only consider the case of $[\psi]$ = 0 and change the sign of τ to obtain the solution for $[\psi]$ = $[\pi]$ . Thus

$[\tau = {{2\sigma - 2} \over {\sin\omega}} = {{2(\cos\omega - 1)} \over {\sin\omega}},]$

$[\det B_n = 2(\cos\omega - 1) \, {{\sin n\omega} \over {\sin\omega}}.]$

Now we show that the above equation has n roots, so there is no reason to consider the case of $[|\sigma|]$ > 1.

Taking into account (28) we find the first root, $[\omega]$ = 0. The other roots can be found from $[\sin (n \omega)]$ = 0, so $[\omega]$ = $[k\pi/n]$ and k = $[0,\ldots,n-1]$ . Note that, due to (29), $[k\ne n]$ . As a result, det B_n = 0 if $[\sigma]$ = $[\cos ({{k \pi}/{n}})]$ , k = $[0,\ldots,n-1]$ . All eigenvalues of matrix B_n can be written as $[\lambda_k]$ = $[2(\sigma + \rho_k)]$ , where $[\rho_k]$ = $[\cos({{k\pi}/{n}})]$ , k = $[1,\ldots,n]$ .

For a given eigenvalue $[\lambda_k]$ we find an eigenvector y = $[\{y_1,\ldots,y_n\}]$ . In this paper we consider only the case of odd n. We need to solve the system of n linear equations

$[\left\{ \matrix{ (1+2\rho_k) y_1 + y_2 = 0_{\vphantom{\big|}},\hfill \cr y_{q-1} + 2\rho_k y_q + y_{q+1} = 0, \qquad q = 2,\ldots, n-1_{\vphantom{\big|}},\hfill \cr y_{n-1} + (1+2\rho_k) y_n = 0.\hfill } \right. \eqno(54)]$

Let us denote $[\theta]$ = $[k\pi/2n]$ and check that for k = n, $[\rho_k]$ = 1 and

$[y_q = 1 \eqno(55)]$

is the solution of (54); for k < n,

$[y_q = (-1)^q\sin(2q-1)\theta. \eqno(56)]$

If k is odd, then y_n+1-q = y_q, i.e. vector y is symmetric. If k is even, then y_n+1-q = -y_q, i.e. vector y is antisymmetric. Thus it is enough to check the case of q < n/2. We obtain y₁ = $[-\sin\theta]$ , y₂ = $[\sin 3\theta]$ = $[\sin\theta(3-4\sin^2\theta) = \sin\theta(2\cos 2\theta +1)]$ , $[2\rho_k + 1]$ = $[2\cos2\theta + 1]$ , so the first equation of (54) is true. As

$[\sin(2q-3)\theta + \sin(2q+1)\theta = 2\sin(2q-1)\theta\cos2\theta,]$

then the second equation in (54) is also true.

Now we find the norm of vector y. For k = n, it is n, and, for k < n,

$[\sum_{q = 1}^n y_q^2 = \sum_{q = 1}^n\sin^2(2q-1)\theta = {{1} \over {2}}\sum_{q = 1}^n\left[1-\cos2(2q-1)\theta\right].]$

Substituting a = $[2\theta]$ and b = $[4\theta]$ into (26) and taking into account that $[\sin 2n\theta]$ = $[\sin2k\pi]$ = 0, we find that the norm for k > 0 is n/2. Finally, we obtain the normalized vector y with elements

$[y_q = \left\{ \matrix{ {{(-1)^{q-1}} \over {\sqrt{n/2}_{\vphantom{\big|}}}} \sin{{(2q-1) k \pi} \over {2n}},\hfill & k = 1,\ldots,n-1,\hfill \cr {{1} / {\sqrt{n}}},\hfill & k = n.\hfill } \right.]$

B1.8. Inverse matrix B⁻¹

The original matrix $[B \equiv B_n]$ from (52) can be written as

$[B = Y D Y^{\,T},]$

where matrix Y is the orthogonal matrix, i.e. Y^-1 = Y^T, with elements

$[Y_{i j} = \left\{ \matrix{ {{(-1)^{i-1}} \over {\sqrt{n/2}_{\vphantom{\big|}}}} \sin{{(2i-1) j \pi} \over {2n}},\hfill & j = 1,\ldots,n-1,\hfill \cr {{1}/{\sqrt{n}}},\hfill & j = n.\hfill } \right.]$

and matrix D is the diagonal one with elements

$[D_{ii} = 2\left[ \sigma + \cos({{{i\pi}/{n}}}) \right],]$

(see Higham, 2008). The inverse of matrix B can then be written as

$[B^{\,-1} = Y D^{\,-1} Y^{\,T},]$

where D^-1 is a diagonal matrix with elements

$[D^{\,-1}_{ii} = {{1} \over {2[\sigma + \cos(i\pi/n)]}}.]$

Therefore the element (k,j) of matrix D^-1 Y^T,

$[D^{\,-1} Y^{\,T}|_{kj} = \sum_{q\,=\,1}^n D^{\,-1}_{k q} Y^{\,T}_{q j} = D^{\,-1}_{k k} Y^{\,T}_{k j} = D^{\,-1}_{k k} Y_{j k}]$

and the elements of B^-1 can be written as

$[B^{\,-1}_{ij} = \sum_{k\,=\,1}^n Y_{i k} D^{\,-1}_{k k} Y_{j k}.]$

Then

$[\eqalign{ &B^{\,-1}_{ij} = \cr& {{1}\over{n}} \left[\sum_{k\,=\,1}^{n-1}{{{(-1)^{i+j}} \over {\sigma + \cos (k\pi/n)}}} \sin{{{(2 i -1) k \pi} \over {2 n}}} \sin{{{(2 j -1) k \pi} \over {2 n}}} + {{{1} \over {2(\sigma - 1)}}}\right].}]$

Taking into account that $[\sin[(2i-1)(2n-k)\pi/2n]]$ = $[\sin[(2i-1)k\pi/2n]]$ and $[\cos(n\pi/n)]$ = −1, we find

$[B^{-1}_{ij} = {{1} \over {2n}} \sum_{k\,=\,1}^{2n-1}{{{(-1)^{i+j}}\over{\sigma+\cos({k\pi}/{n}) }}} \sin{{{(2 i -1) k \pi} \over {2 n}}} \sin{{{(2 j -1) k \pi} \over {2 n}}}.]$

All matrices mentioned above have indices starting from 1. However, our aim is to find a filter to be applied to a matrix. Therefore it may be more convenient to rewrite the same matrices as ones centered with respect to the central element. Since n and m are both odd, this central element always exists and is unique. If a standard matrix Q is given, then by $[\check Q]$ we denote a matrix centered with respect to its central element. So, instead of matrix B^-1_ij defined for i = $[1,\ldots,n]$ , j = $[1,\ldots,n]$ , we subtract n from k and $[(\check n + 1)]$ from the other indices of the above formula, where $[\check n]$ $[\equiv]$ (n-1)/2, and obtain

$[\eqalign{\check B^{-1}_{ij} = {}& {{1}\over{2n}} \sum_{k\,=\,-(n - 1)}^{n-1}\,\, {{{(-1)^{i+j}} \over {\sigma + \cos[(k+n)\pi/n] }}}\cr&\times \sin{{{[2 (i+\check n) +1] (k+ n) \pi} \over {2 n}}} \sin{{{[2 (j + \check n) +1] (k +n) \pi} \over {2 n}}},}]$

where i, j = $[-\check n, \ldots, \check n]$ . We may rewrite

$[\sin{{{[2 (i+\check n) + 1] (k+ n) \pi} \over {2 n}}} = (-1)^{i+\check n}\sin\left({{i k\pi} \over {n}} +{{k+1} \over {2}}\pi\right),]$

$[\cos{{{(k+n)\pi} \over {n}}} = -\cos {{{k \pi} \over {n}}},]$

$[\eqalign{& \check B^{\,-1}_{ij} = {{1}\over{2n}} \times \cr& \sum_{k\,=\,-(n - 1)}^{n-1} \!\!\!{{\sin\left\{({ik\pi}/{n}) +[({k-1})/{2}]\pi\right\} \sin\left\{({\,jk\pi}/{n}) +[({k-1})/{2}]\pi\right\}} \over {\sigma - \cos({k\pi}/{n})}}\cr {} & \!\qquad= {{1} \over {4n}} \sum_{k\,=\,-(n - 1)}^{n-1} {{\cos[{(i-j)k\pi}/{n}] + \cos[{(i+j+n)k\pi}/{n}]} \over {\sigma - \cos({k\pi}/{n}) }}. }]$

For the central row of $[\check B^{\,-1}]$ we obtain

$[\eqalign{\check B^{\,-1}_{0j} & = {{1} \over {4n}} \sum_{k\,=\,-(n - 1)}^{n-1} {{\cos({\,jk\pi}/{n})[1+(-1)^k]} \over {\sigma - \cos({k\pi}/{n})}}. }]$

For odd k the value of [1+ (-1)^k] = 0, so we need to consider only even k and use s = k/2,

$[\check B^{\,-1}_{0j} = {{1} \over {2n}} \sum_{s\,=\,-\check n}^{\check n} \,\,{{\cos({2js\pi}/{n})} \over {\sigma - \cos({2s\pi}/{n})}}. \eqno(57)]$

B1.9. Inversion of the matrix function

The matrix W_n(x) is defined in (16). According to (57) the central row of the function W_n^-1(x) can be written as

$[W^{\,-1}_n|_{0 j} = {{1} \over {n}} \sum_{s\,=\,- \check n}^{\check n} \cos {{2j s\pi} \over {n}} \left[x + 1 - 2\cos ({{2s\pi}/{n}}) \right]^{-1}. \eqno(58)]$

Now we need to find the matrix $[W_n^{\,-1}[W_m(\beta)]]$ . If x and v are two scalars, then W_m(x) + v = W_m(x+v). Therefore,

$[\left[W_m(x) + 1 - 2\cos { ({{ 2s \pi}/{n}})} \right]^{-1} = W^{\,-1}_m \left[x + 1 - 2\cos {({{ 2s \pi}/{n}})}\right].]$

In principle we only need the central row of (nm×nm) matrix $[W_n^{\,-1}[W_m(\beta)]]$ . Formula (58) can be used to find n (m×m) matrices which form the central (m×nm) block of $[W_n^{\,-1}[W_m(\beta)]]$ . So for the jth matrix $[\check R]$ we may write

$[\eqalign{ \check R_{0 k} & = {{1}\over{nm}} \sum_{s\,=\,-\check n}^{\check n} \sum_{r\,=\,-\check m}^{\check m} \,\,{{\cos({2js\pi}/{n}) \cos({2kr\pi}/{m}) }\over{ \beta+2\left[1-\cos {({2s\pi}/{n})} - \cos {({2r\pi}/{m})} \right] }}, }]$

where $[\check m]$ = (m-1)/2. Therefore the ( j m+k)th element of the central row of the matrix A^-1 defined in (18) can be written as

$[{{\beta-2}\over{nm}} \sum_{s\,=\,-\check n}^{\check n} \sum_{r\,=\,-\check m}^{\check m} \,\,{{\cos({2js\pi}/{n}) \cos({2kr\pi}/{m}) }\over{ \beta+2\left[1-\cos{({2s\pi}/{n})} - \cos{({2r\pi}/{m})} \right]}}.]$

If values n and m are relatively large, then instead of the sum we may consider a 2D integral. So define continuous variables ξ and η such that $[\xi,\eta \in[-\pi,\pi]]$ and s = $[{{\xi n}/{2\pi}}]$ , r = $[{{\eta m}/{2\pi}}]$ and introduce an integral

$[\check G_{jk} = {{\beta - 2} \over {4\pi^2}} \int_{-\pi}^\pi \int_{-\pi}^\pi \,\, {{\cos j\xi \cos k \eta} \over { \beta + 2\left(1 - \cos \xi - \cos \eta \right)}}\,{\rm{d}}\eta\,{\rm{d}}\xi,]$

then the elements of the central row of matrix A^-1 tend to $[\check G_{jk}]$ for large n and m. We denote $[\tau]$ = $[\alpha/(1+4\alpha)]$ and rewrite the last formula as

$[\check G_{jk} = {{1-4\tau} \over {\pi^2}} \int_{0}^\pi \int_{0}^\pi {{\cos j\xi \cos k \eta} \over { 1 - 2\tau \left(\cos \xi + \cos \eta \right)}}\,{\rm{d}}\eta\,{\rm{d}}\xi. \eqno(59)]$

Some useful properties of (59) can be found. As cosine is an even function, then $[\check G_{jk}]$ is symmetrical with respect to horizontal, vertical or diagonal flipping, i.e. $[\check G_{j k}]$ = $[\check G_{-j k}]$ = $[\check G_{j, -k}]$ = $[\check G_{k j}]$ . Thus it is enough to consider $[k\geq j\geq 0]$ . Formulas (49) and (50) give us

$[\check G_{j k} = \tau(\check G_{j-1 k} + \check G_{j+1 k} + \check G_{j k-1} + \check G_{j k+1}), \eqno(60)]$

$[\check G_{0 0} = \tau(4 \check G_{0 1} + 1).]$

In Appendic B1.5 it is shown how the integral in (59) can be calculated, so

$[\check G_{jk} = (1-4\tau)\tau^{\,k+j}\sum_{q\,=\,0}^{\infty} C_{2q + j + k}^{\,q} C_{2q + j + k}^{\,q+j} \tau^{\,2 q}. \eqno(61)]$

From the above formula we obtain other properties of $[\check G_{jk}]$ :

(i) $[\check G_{jk}]$ > 0, since $[\tau]$ < 1/4 and binomial coefficients are non-negative.

(ii) $[\check G_{j, k+1}]$ < $[\check G_{j, k}]$ , since

$[{{C_{2q + j+ k}^{\,q}} \over {C_{2q + j+ k +1}^{\,q}}} = {{q+j+k+1} \over {2 q + j + k + 1}}\,\,\gt\,\,{{1} \over {2}},]$

$[{{C_{2q + j+ k}^{\,q + j}} \over {C_{2q + j+ k +1}^{\,q+j}}} = {{q+k+1} \over {2 q + j + k + 1}}\,\,\gt\,\,{{1} \over {2}}]$

and therefore

$[\tau C_{2q + j + k +1}^{\,q} C_{2q + j + k +1}^{\,q+j} \,\,\lt\,\, C_{2q + j+ k}^{\,q} C_{2q + j+ k}^{\,q+j}.]$

(iii) The sum S_j of all elements of the jth row of matrix $[\check G]$ ,

$[\sum_{k\,=\,-\infty}^{\infty} \check G_{j k} = \sqrt{1 - 4\tau} \left({{1-2\tau - \sqrt{1 - 4\tau}} \over {2\tau}}\right)^{|j|}. \eqno(62)]$

Taking into account formula (27) and the property (43), we obtain

$[S_j = {{1-4\tau} \over {\pi}} \int_{0}^\pi {{\cos j\xi} \over { 1 - 2\tau \left(\cos \xi + 1\right)}}\,{\rm{d}}\xi.]$

If we introduce φ such that $[\sin\varphi]$ = $[2\tau/(1-2\tau)]$ , $[\cos\varphi]$ = $[\sqrt{1-4\tau}/(1-2\tau)]$ , then $[\varphi\in(0,{{\pi}/{2}})]$ , since $[0\,\lt\,\tau\,\lt\,{{1}/{4}}]$ , and the value of S_j follows from the results of Appendix B1.6.

(iv) The sum of all elements of matrix $[\check G]$ is 1,

$[\sum_{j\,=\,-\infty}^{\infty} S_j = 1 \qquad{\rm{or}}\qquad \sum_{j\,=\,-\infty}^{\infty} \sum_{k\,=\,-\infty}^{\infty} \check G_{jk} = 1,]$

since (62) and (30) for $[\gamma]$ = $[(1-2\tau - \sqrt{1 - 4\tau})/2\tau]$ give

$[\sum_{j\,=\,-\infty}^{\infty} \gamma^{|j|} = {{2} \over {1-\gamma}}-1 = {{1+\gamma} \over {1-\gamma}} = {{1} \over {\sqrt{1-4\tau}}}.]$

APPENDIX C

Numerical implementation

C1. High-precision calculations

In order to find elements $[\check G_{jk}]$ one can use a finite sum in (61). However, due to various rounding operations, there will be some additional error even if double precision numbers are used. Therefore it is better to try another approach. Formula (61) can be rewritten as $[\check G_{jk}]$ = $[(1-4\tau)\tau^{\,j+k}\chi]$ ,

$[\chi = \sum_{q\,=\,0}^{\infty} a_q \tau^{\,2q}, \qquad a_q = C_{2q + j + k}^{\,q} C_{2q + j + k}^{\,q+j}. \eqno(63)]$

or a sum $[\chi]$ = $[\chi_1 + \chi_2]$ of two terms

$[\chi_1 = \sum_{q\,=\,0}^{q^\star-1} a_q \tau^{\,2q}, \quad \chi_2 = \sum_{q\,=\,q^\star}^{\infty} a_q \tau^{\,2q},]$

where index $[q^\star]$ can be chosen later. The first term can be found using the recurrent formula

$[\eqalign{ c_{q^\star} & = 1,\cr c_q & = 1 + {{a_q} \over {a_{q-1}}} \,\tau^{\,2} c_{q+1}, \quad q = q^\star -1,\ldots, 1}]$

and since a₀ = C_{j + k}⁰ C_{j + k}^j = C_{j + k}^j, then $[\chi_1]$ = a₀c₁ = C_j+k^j c₁. Now we estimate from above the value of $[\chi_2]$ . We obtain

$[{{a_q} \over {a_{q-1}}} = {{(2q + j + k)^2(2q+ j + k - 1)^2} \over {q(q+j)(q+k)(q+j+k)}}.]$

Note that $[(q+j)(q+k) \geq q(q+j+k)]$ , so if we denote s = j+k we can estimate

$[{{a_q} \over {a_{q-1}}} \leq \left[{{(2q + s)(2q+ s - 1)} \over {q(q+s)}}\right]^2.]$

For $[q\geq s(s-1)/2]$ , we obatin

$[{{a_q} \over {a_{q-1}}} \leq 16]$

and $[({{a_q}/{a_{q-1}}})\tau^{\,2}]$ < 1, so we can estimate from above the value of $[\chi_2]$ ,

$[\chi_2 \leq a_{q^\star}\tau^{\,2q^\star} \sum_{s\,=\,0}^{\infty} (16\tau^{\,2})^s = {{a_{q^\star}\tau^{\,2q^\star}} \over {1- 16\tau^2}}.]$

As a result the value $[\check G_{jk}]$ can be found by the following formula,

$[\eqalign{ d_{q^\star} & = {{1} \over {1-16\tau^{\,2}}},\cr d_q & = 1 + d_{q+1} {{{(2q + j + k)^2(2q+ j + k - 1)^2} \over {q(q+j)(q+k)(q+j+k)}}},}]$

for $[ q = q^\star -1,\ldots, 1,]$

$[\check G_{jk} = (1-4\tau)\tau^{\,j+k} C_{j+k}^{\,j} d_1. \eqno(64)]$

C2. Faster calculations

Due to the symmetry of $[\check G_{jk}]$ we have $[\check G_{-1 k}]$ = $[\check G_{1 k}]$ , so from (60) we obtain

$[\check G_{1k} = {{1} \over {2}}\left({{\check G_{0k}} \over {\tau}} - \check G_{0 k -1} - \check G_{0 k+1}\right). \eqno(65)]$

The other rows can be found as

$[\check G_{jk} = {{\check G_{j-1, k}} \over {\tau}} - \check G_{j-1, k -1} - \check G_{j-1, k+1} - \check G_{j-2, k}. \eqno(66)]$

Acknowledgements

The author would like to thank Albrecht Kyrieleis, Philip Withers and beamline scientists of station 16.3 of Synchrotron Radiation Source (Daresbury Laboratory) and I12 beamline of the Diamond Light Source for the data provided.

References

Altunbas, C., Lai, C.-J., Zhong, Y. & Shaw, C. C. (2014). Med. Phys. 41, 091913. CrossRef PubMed Google Scholar
Baek, J., De Man, B., Harrison, D. & Pelc, N. J. (2015). Opt. Express, 23, 7514–7526. CrossRef PubMed Google Scholar
Benjamin, A. T. & Quinn, J. J. (2003). Proofs that Really Count: The Art of Combinatorial Proof. Washington: The Mathematical Association of America. Google Scholar
Brun, F., Accardo, A., Kourousias, G., Dreossi, D. & Pugliese, R. (2013). 8th International Symposium on Image and Signal Processing and Analysis (ISPA 2013), Trieste, Italy, 4–6 September 2013, pp. 672–676. Google Scholar
Davis, G. R. & Elliott, J. C. (1997). Nucl. Instrum. Methods Phys. Res. A, 394, 157–162. CrossRef CAS Web of Science Google Scholar
Davis, G. R., Evershed, A. N. Z. & Mills, D. (2013). J. Dent. 41, 475–482. CrossRef CAS PubMed Google Scholar
Drakopoulos, M., Connolley, T., Reinhard, C., Atwood, R., Magdysyuk, O., Vo, N., Hart, M., Connor, L., Humphreys, B., Howell, G., Davies, S., Hill, T., Wilkin, G., Pedersen, U., Foster, A., De Maio, N., Basham, M., Yuan, F. & Wanelik, K. (2015). J. Synchrotron Rad. 22, 828–838. Web of Science CrossRef CAS IUCr Journals Google Scholar
Engl, H. W., Hanke, M. & Neubauer, A. (1996). Regularization of Inverse Problems, Vol. 375 of Mathematics and its Applications. Dordrecht: Kluwer Academic Publishers Group. Google Scholar
Gradshteyn, I. S. & Ryzhik, I. M. (2015). Table of Integrals, Series and Products. Amsterdam: Elsevier/Academic Press. Google Scholar
Higham, N. J. (2008). Functions of Matrices: Theory and Computation. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics. Google Scholar
Hinebaugh, J., Challa, P. R. & Bazylak, A. (2012). J. Synchrotron Rad. 19, 994–1000. CrossRef CAS IUCr Journals Google Scholar
Kak, A. & Slaney, M. (2001). Principles of Computerized Tomographic Imaging. Philadelphia: Society for Industrial and Applied Mathematics. Google Scholar
Kalender, W. A. (2006). Phys. Med. Biol. 51, R29–R43. CrossRef PubMed Google Scholar
Miqueles, E. X., Rinkel, J., O'Dowd, F. & Bermúdez, J. S. V. (2014). J. Synchrotron Rad. 21, 1333–1346. Web of Science CrossRef CAS IUCr Journals Google Scholar
Münch, B., Trtik, P., Marone, F. & Stampanoni, M. (2009). Opt. Express, 17, 8567–8591. Web of Science PubMed Google Scholar
Natterer, F. & Wübbeling, F. (2001). Mathematical Methods in Image Reconstruction. Philadelphia: Society for Industrial and Applied Mathematics. Google Scholar
Paleo, P. & Mirone, A. (2015). J. Synchrotron Rad. 22, 1268–1278. CrossRef IUCr Journals Google Scholar
Pešić, Z. D., Fanis, A. D., Wagner, U. & Rau, C. (2013). J. Phys. Conf. Ser. 425, 182003. Google Scholar
Prell, D., Kyriakou, Y. & Kalender, W. A. (2009). Phys. Med. Biol. 54, 3881–3895. Web of Science CrossRef PubMed Google Scholar
Rashid, S., Lee, S. Y. & Hasan, M. K. (2012). Eurasip J. Adv. Signal. Process. 2012, 93. CrossRef Google Scholar
Sadi, F., Lee, S. Y. & Hasan, M. K. (2010). Comput. Biol. Med. 40, 109–118. Web of Science CrossRef PubMed Google Scholar
Sijbers, J. & Postnov, A. (2004). Phys. Med. Biol. 49, N247–N253. Web of Science CrossRef PubMed Google Scholar
Tikhonov, A. N. & Arsenin, V. Y. (1977). Solutions of Ill-Posed Problems. Washington: John Wiley and Sons. Google Scholar
Tikhonov, A. N., Goncharsky, A. V., Stepanov, V. V. & Yagola, A. G. (1995). Numerical Methods for the Solution of Ill-Posed Problems, Vol. 328 of Mathematics and Its Applications. Netherlands: Springer. Google Scholar
Titarenko, V. (2016). IEEE Signal Process. Lett. 23, 800–804. CrossRef Google Scholar
Titarenko, V., Bradley, R., Martin, C., Withers, P. J. & Titarenko, S. (2010c). Proc. SPIE, 7804, 78040Z. Google Scholar
Titarenko, S., Titarenko, V., Kyrieleis, A. & Withers, P. (2009). Appl. Phys. Lett. 95, 071113. CrossRef Google Scholar
Titarenko, S., Titarenko, V., Kyrieleis, A. & Withers, P. J. (2010d). J. Synchrotron Rad. 17, 540–549. CrossRef IUCr Journals Google Scholar
Titarenko, S., Titarenko, V., Kyrieleis, A., Withers, P. & De Carlo, F. (2011). J. Synchrotron Rad. 18, 427–435. Web of Science CrossRef IUCr Journals Google Scholar
Titarenko, V., Titarenko, S., Withers, P. J., De Carlo, F. & Xiao, X. (2010a). J. Synchrotron Rad. 17, 689–699. CrossRef CAS IUCr Journals Google Scholar
Titarenko, S., Withers, P. J. & Yagola, A. (2010b). Appl. Math. Lett. 23, 1489–1495. Web of Science CrossRef Google Scholar
Warnett, J. M., Titarenko, V., Kiraci, E., Attridge, A., Lionheart, W. R. B., Withers, P. J. & Williams, M. A. (2016). Meas. Sci. Technol. 27, 035401. CrossRef Google Scholar
Wei, Z., Wiebe, S. & Chapman, D. (2013). J. Instrum. 8, C06006. Google Scholar
Wolkowski, B., Snead, E., Wesolowski, M., Singh, J., Pettitt, M., Chibbar, R., Melli, S. & Montgomery, J. (2015). J. Synchrotron Rad. 22, 1130–1138. CrossRef CAS IUCr Journals Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 23| Part 6| November 2016| Pages 1447-1461

https://doi.org/10.1107/S160057751601479X

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Search IUCr Journals		doi		Advanced search
Author		volume	page

research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Analytical formula for two-dimensional ring artefact suppression

1. Introduction

2. Brief mathematical formulation

2.1. Tikhonov functional

2.2. Optimization problem

2.3. Original 1D problem

2.4. Mathematical formulation of 2D problem

3. Example problems

3.1. Graphite

3.2. TRISO particles

3.3. Mouse skull

4. Practical implementation

5. Discussion

APPENDIX A

A1. Two eigenvalues for matrix A

APPENDIX B

B1. Formula for 2D filter

B1.1. Trigonometric formulas

B1.2. Combinatorial formulas

B1.3. Simple integrals

B1.4. The first sequence

B1.5. The second sequence

B1.6. The third sequence

B1.7. Eigenvalues/eigenvectors

B1.8. Inverse matrix B−1

B1.9. Inversion of the matrix function

APPENDIX C

Numerical implementation

C1. High-precision calculations

C2. Faster calculations

Acknowledgements

References

research papers

B1.8. Inverse matrix B⁻¹