Sinogram Interpolation Method for Sparse-Angle Tomography

In sparse-angle X-ray tomography reconstruction, where only a small number of projection images are taken around the object, appropriate sinogram interpolation has a significant impact on image quality. A novel sinogram interpolation method is introduced for extreme sparse tomographic reconstruction where only nine measured projection images are available. The sinogram is interpolated by solving characteristics of the so-called warps, which can be considered as approximation sine waves in a limited region. The numerical evidence suggests that this approach gives superior results over standard interpolation methods when the tomographic data are extremely sparse and noisy.


Introduction
Three-dimensional (3D) X-ray imaging is typically needed for measuring the inner structures of a patient.For example, in implant dentistry where damaged or missing teeth are replaced by artificial teeth, a screw hole with accurate depth, angle and diameter is needed for the implant.The hole has to be deep enough to ensure firm attachment, but not too deep as nerves might then be damaged.This measurement cannot be based on two-dimensional (2D) images since they represent a distorted projection of the tissue.Another clinical application of 3D X-ray imaging is solving the superposition problem.Since a pixel of a 2D projection image is a sum of attenuation along the path of an X-ray, the overlapping low-contrast tissues are difficult to recognize.This can easily lead to misinterpretation and eventually cause a false diagnosis.However, in 3D imaging, the viewing angle can be chosen so that the boundaries between tissues can be accurately identified.
Lately, there has been an increasing interest in sparse-angle tomographic imaging, where fewer projection images are taken in order to reduce patient dose and scanning time.Despite that the sparse reconstruction does not offer high spatial resolution, it can offer a feasible overview on the imaging target, which extends the use of computed tomography (CT) imaging to less serious cases, like minor trauma studies or cosmetic operation planning.Moreover, sparse-angle imaging can also be used in pre-scans, where a low-resolution reconstruction is generated before the actual scan to verify patient positioning and to define optimal X-ray technique values (tube current, voltage, number of projection images) for the final scan.The pre-scan helps to optimize the patient dose and avoid re-exposures.
The challenge of sparse tomographic imaging is insufficient image quality caused by the inadequate number of projection images.This lack of data is illustrated in Figure 1 in the spatial, sinogram and frequency domains.In the final reconstruction, the missing information can be observed as artifacts in the 2D slice view as well as in the 3D volume view.In the slice view, the missing information generates streak artifacts (see Figures 2-4   especially when sharp edges are present [1].Numerical implementations of the tomographic reconstruction from sparse data are introduced for example in articles by Hyvönen et  (i.e.mechanical combination of the detector and the X-ray source) and the object.Finally, in the post-processing stage, the reconstruction artifacts are reduced and the clinically relevant information is emphasized.A practical and clinical overview of the reconstruction process can be found in the book by Hsieh (2009) and more theoretical approach in the book by Kak and Slaney (2001) [10,11].
A sinogram is a 2D representation of an X-ray CT scan, where each column represents a single row in the projection image arranged in increasing angular order.Typically, the horizontal axis represents the angle of the X-ray detector and the vertical axis, the distance from the rotation center in the detector row.Each volume element generates a sine wave in the sinogram when the X-rays are orthogonal to the detector (parallel beam geometry).Therefore, the sinogram consists of a number of overlapping sine functions, whose amplitude and phase depend on the location of the voxel and the grayvalue equals to the grayvalue of the volume element while the wavelength is constant ( 2π ).As a simple example of a two-point sinogram seen in Figure 5, Figure 1 illustrates the relation between spatial and sinogram domains, and the upper row in Figure 6 is an example of a sinogram arising from sparsely sampled data.See chapter 5.11.3 in the book by Gonzalez and Woods (2008) for more information about sinogram representation [12].
In this study, we focus on interpolating new sinogram columns without modifying the measured sinogram columns.This scenario is consistently called in this paper as sinogram interpolation.In the literature, there are also a number of other interpolation routines related to sinogram manipulation which should not be confused with the sinogram interpolation scenario.One of them is the projection image interpolation, where new values are generated between the data values to avoid Moire artifacts.This is typically handled by the standard nearestneighbor or linear interpolation [1].Another example is needed for an interpolation routine in helical CT device for gaining uniform sampling in the azimuthal direction or metal artifact removal processes [13][14][15][16][17][18].
Also, the sinogram extrapolation can be used for solving the truncation problem as indicated for example by Sze and Shum (1996), Gilland et al. (2000) or Baojun and Jiang (2007) to expand the projection data in limited angle tomography [19][20][21].Moreover, sinogram extrapolation can be used for expanding the truncated sinogram in a local tomography situation, where some of the attenuations take place outside the modeled volume (see Zamyatin and Nakanishi (2006) and references therein [22]).
Sinogram interpolation itself has been studied since the innovation of the medical CT system in the 1970's.For example, Brooks et al. (1978) studied both interpolating new data in projection images and interpolating new projection images to attenuate Moire effects [1].Also, Lahart (1981) studied a similar problem and used a least-squares approximation to interpolate projection images when the external shape of the object was known [23].However, in those days, the computational resources were very limited and therefore the spatial resolution as well as the image quality was not equivalent to today's de-facto standards.Nevertheless, these studies illuminate the limitations of sparse and limited angle tomographic imaging and therefore can be considered as fundamental studies of sinogram estimation (see also [24][25][26][27]).[28,29].In this study, about 60% of the sinogram columns were interpolated.Similarly, in the study by Penbel et al. (2005), only 50% of the sinogram columns were interpolated by using cumulant functions instead of linear interpolation [30].Moreover, Bertram et al. (2004 and2009) utilized tensors to gain directional interpolation based on sinogram data and achieved impressive results [31,32].
In this paper, we introduce a new approach, called Sinogram Interpolation Technique (SINT) for estimating new sinogram columns between known, measured sinogram columns in the extremely sparse situation of only nine known sinogram columns.
This article is organized as follows: In Section 1.2, the theoretical background of the SINT method is described.Secondly, in the Section 1.3, the capability and limitations of SINT by three different numerical examples with additional noise are demonstrated.Finally, in the Section 1.4, we discuss the computational and clinical aspects of SINT and its possible future applications.Also, we analyze the failure of standard interpolation methods in extremely sparse imaging situations.
Throughout this paper, capital letters are used for matrices.Moreover, for indicating a single element in a matrix, the standard [row][column]-order is applied.For example, , k h

S
indicates the element located at the k th row and h th column in the matrix S .Notwithstanding the aforesaid, a single subscript is used for referring to columns: h S denotes the h th column of the matrix S .Superscripts in brackets are used for additional information, like iteration round, special restrictions or relations.See Table 1 for the definition of each notation used in this work.We systematically use the term known sinogram column for the columns that are based on actual measured X-ray projection images.The interpolated sinogram columns are called estimated columns since they are a priori unknown.Furthermore, the known sinogram columns adjacent to an estimated column are called known neighbor columns.See also Figure 7 for an illustration of this terminology.
We use the term volume for the target object and image for the projection data, although in our examples the phantoms are two-dimensional and projection images are one-dimensional.Similarly, we use the terms voxel and pixel to describe the volume and image elements to be consistent with other tomographic literature and references.Finally, in this study, we are focusing on 2D parallel beam imaging situations, where all X-rays are perpendicular to the detector.Also, we are covering only cases where the whole object is covered in each projection image, and therefore the so-called local tomography problem is excluded from this study.

The SINT Method
In this Section the theoretical background of the SINT method is discussed.Firstly, in Section 2.1 two new concepts, called warp and weight factor, are introduced.Secondly, in the Section 2.2 a new method for defining the grayvalues in single sinogram column located between two a priori known neighbor columns is proposed.Finally, in Section 2.3 the method is summarized and extended to cover all sinogram columns in the sinogram domain.Weight factor

Weight Factor
In this study a sinogram is considered as a matrix

S S S S S S S S
and known projection angles corresponding to the columns of S as a vector [ ] where The purpose of the sinogram integration is to generate a feasible approximation for h S when other columns and vector θ are known.In this method the grayvalues of each unknown sinogram element 2 , i h S are approximated as weighted sums of nearby sinogram elements in the adjacent columns 1 h S − and 1 h S + .This can be determined in two ways, each producing the same result for h S ; either based on the column or based on the column where ( , whose solution is given by We end this chapter by highlighting two important identities related to weight factors.Firstly, consider a single sin wave that intersects the sinogram elements Secondly, because the total sum of each sinogram column is constant, the sum of all weight factors related to the same known sinogram component in columns 1 h − and 1 h + is always one, i.e. ( ) 1, 1 for each .

Warp
In extremely sparse-angle tomographic situation, the number of sine functions and their characteristic parameters (amplitude, grayvalue and phase) cannot be uniquely defined from the sinogram.For that reason, a new concept called warp is introduced.The warp can be considered as a sum of all the sine waves that intersects two sinogram elements 1 h S − and 1 h S + (see Figure 7).The most essential features of warps are: 1) Each warp has the shape of a sine wave and a wavelength of 2π .Therefore, the phase and amplitude of the warp are uniquely defined by two points.
2) Each warp connects two sinogram elements in known neighbors columns Our interpolation strategy is to define the estimated sinogram column h S by defining the amplitude, phase and so-called warp factor (will be explained in the Section 2.2) for each warp.We will show that by defining these values we can define the weight factors and thereby obtain an estimate for the sinogram column h S .Although the above warp-based approach does not completely characterize the underlying sine waves and provide an exact solution to the interpolation problem, it does give a feasible and stable approximation in limited-data situations as demonstrated in Section 1.3.

Definition of the Crucial Matrices
In this Section we will introduce three matrices to determine the connection between the warps and grayvalues of the estimated sinogram ( ) h S .These matrices are: 1) Matrix A to define the relation between the estimated and neighbor known sinogram elements based on the Equation (1.6) 2) Matrix G to fulfill the requirement in the Equation (1.7) 3) Matrix R to limit the weight factors related to the same known sinogram based on the Equation (1.8) We also build two temporary matrices Ψ and Φ to simplify the numerical implementation.All these matrices should be generated for each sinogram column separately as described in Section 2.3.

Warp Matrix Φ
The purpose of the warp matrix Since the warp has the shape of a sine wave with a frequency of one and an offset of ( ) , it can be defined as where the amplitude k ϑ and phase k φ depend on the known angles 1 h θ ± as follows: where ( ) ( ) 1 2 and 1 2 .
, then the warp is valid (see the requirement item 3 in the Section 2.1) and finally ( )  where [ ] ⋅ the rounding operator.Because this limitation, the height of the matrix Φ typically smaller than NN .The procedure for generating the warp matrix Φ is described in the Algorithm 1.

Relation Matrices A and Ψ
The purpose of this phase is to define two quite similar sparse matrices; Relation Grayvalue Matrix , where and Relation Index Matrix The relation matrix A can be constructed from the row vectors ( ) such that each row vector ( ) j A  consists of the grayvalues of known neighbor sinogram columns that are connected to an estimated sinogram element related to exactly one estimated sinogram element, only one non-zero value on each column is allowed.Therefore, the matrix A has the form where 0 are zero row vectors.
We also need the row indexes of the known points (indicated above as 1 i and 3 i ) ordered consistently with the structure of the matrix A .Those values will be stored into matrix . The procedure for generating the matrices A and Ψ is described in Algorithm 2.
At this point the matrices A and Ψ have been generated.The matrix A defines the relation between the estimated sinogram column and the weight factors as described in Equation (1.13).The matrix Ψ plays an important role in defining the eigenvector and rule matrices in the following Sections.

Eigenvector Matrix Λ
Having constructed the matrix A appearing in Equation (1.13), we now turn to specifying the weight vector w .From the Equation (1.7) can be seen that where and .solution.For additional information about matrices, eigenvectors and eigenvalues see Golub and Van Loan (1996) or Arfken and Weber (2001) [39,40].Next we will declare that a feasible solution for the Equation (1.15) can also be a linear combination of the eigenvectors.The Eigenvector Matrix

i h i h i i i i i h
is defined as where ( ) ( ) ( ) are all positive eigenvectors of the matrix Λ related to eigenvalue of one as discussed above.Then for every vector which shows that the linear combination of feasible eigenvectors fulfills the requirement (1.7).
The vector α is called as warp vector.Weight vector w can now be determined from the warp vector α by utilizing the Eigenvector Matrix The procedure for generating the matrix G is shown in the Algorithm 3. By combining the Equations (1.13) and (1.18) we get which indicates that if we can define the warp factors α , we can also determine the grayvalues of the estimated sinogram column h S .To calculate the warp factors α instead of weight factors w is a twofold benefit; the requirement (1.7) is now implicitly implemented and the dimension of the original problem is reduced from 2L to L .

Rule Matrix R
So far we have generated matrices to determine the relations between the warp factors and the grayvalues in the estimated column.However, the mission for defining the grayvalues for the column h S has not been accomplished since the warp factor vector α is still unknown.Therefore, in this Section we introduce a Rule Matrix ( 1 N and 2 N are number of non-zero elements in column where the 1 is a vertical vector of ones with a height of 1 2 ˆN N .The principle of generating the Rule Matrix R is following; for each known non-zero sinogram element a single row to the matrix R , say j R , is generated so that the value of each element

Summary of the Method
We summarize this Section by declaring that the grayvalues of the estimated column can be solved from the Equation Since the α cannot be uniquely defined from the Equation 1.22, the problem is ill-posed and it calls for regularization.The Tikhonov regularization algorithm was chosen, since it controls simultaneously residual 1 RGα − and the norm of the solution α (see chapter 2.3 in Kaipio and Somersalo (2005) for further information [41]).The Tikhonov solution for the Equation 1.22 has a form of ( ) ( ) ( ) is an identity matrix.The regularization factor can be based on the fact that in the non-local tomography the sum of each sinogram column is constant for throughout the sinogram.Therefore, the relaxation value β can be calculated from the Equation ( )
The complete Algorithm 5 for defining h S is described in details.To calculate The Equation 1.23 less memory consuming, the Singular Value Decomposition (SVD) can be implemented for the matrix RG such that T RG U V = Σ .Then the Equation 1.23 has a matrix-free representation 2 1, where n is the minimum dimension of the matrix Σ .See for example Theorem 2.5 in Kaipio and Somersalo (2005) for proof [41].So far we have described how to interpolate a single sinogram column.To interpolate all sinogram columns the procedure described in Sections from 2.2 to 2.2 are repeated for each unknown column.To numerically compute that, two nested loops are applied; an inner loop where the neighbor columns are fixed and the location of the estimated sinogram column varies and an outer loop where we change the neighbor columns.Before starting the actual interpolation process, we have to define how many new sinogram columns are needed for the isotropic resolution.If we generate too few new sinogram columns, we do not gain isotropic resolution.On the other hand, too many new sinogram columns just increases the computational burden without any improvements to the final reconstruction quality.
From Figure 1(c) we can see that ( ) where θ ∆ is the angular difference of the projection images, ω is the distance between two frequency components, and N is the number of frequency samples.For simplicity, N is assumed to be odd and θ ∆ is constant.
Based on the Fourier Slice Theorem, to gain uniform resolution in the frequency domain, the distance between frequency components within measurement ω should be equal to the maximum distance between adjacent frequency components d as discussed by Kak and Slaney (2001) and more deeply by Natterer (1986) [11,42].Therefore, by defining d ω = , we can define Equation for optimal angular difference where      is the upward rounding operator.The angle ( ) S can therefore be calculated from the Equation ( ) . The detailed procedure for whole sinogram interpolation is described in the Algorithm 6.

Results
We used the GNU Octave program (version 3.2.3)with the imaging toolbox for numerical implementation of OPEN ACCESS AM Algorithm 6.The SINT algorithm.
the SINT algorithm.The coding and executing environment was Ubuntu distribution version 10.04 including GNU/Linux operating system 2.6.32 and GNOME interface version 2.30.2.The computer that we employed for this study was a commercial 2.4 GHz dual-processor desktop computer with 8 GB of memory.
For the numerical implementation we adopted three synthetic phantoms with a size of 128 128× pixels; Shepp-Logan, Dental Arc and Boxes.The standard Shepp-Logan phantom was chosen since it is widely used for evaluating the reconstruction image quality.The Dental Arc was calculated from full CT reconstruction of an artificial head (i.e.real teeth and skull embedded in an acrylic head) and therefore it can be considered as a reference for the clinical capability of SINT.Thirdly, the Boxes phantom consists of two isolated homogeneous rectangles to demonstrate a simplest possible case.
Furthermore, these three phantoms were chosen to evaluate the capability of SINT because these phantoms generate essentially different kinds of sinogram representations, as shown in Figure 6.Since the Shepp-Logan phantom is a relatively isotropic object, it produces a smoothly varying sinogram.Secondly, the Dental Arc is a strongly anisotropic object, which produces a narrowing and widening sinogram with significant dynamic range in grayvalues.Finally, the Boxes phantom consists of two separate objects generating a branching and merging sinogram.
In all cases, we took nine synthetic projection images from 25 to 185 degree angles by using the built-in radon function.This function calculates the sum of pixels along a line with the given angles.The interpolation factor was set to 32 based on Equation (1.28), which gave about 0.6 degrees angular differentiation for the final interpolated sinogram.We also evaluated the robustness of this method against pixel noise by adding Gaussian noise with variation of 5% of the pixel grayvalue to each known sinogram element.
As reference interpolation methods, we used linear, cubic spline and nearest neighbor -interpolation methods, which are the most common interpolation methods.For further information about these interpolation methods, see Section 4.4 in the book by the Gonzalez and Woods (2008) [30].
As an interpolation quality metric, we calculated relative 2 L norms against the ground truth.That is  2.
We also reconstructed the object from the interpolated data for subjective comparison of reconstruction image quality.The reconstruction was executed by the built-in iradon -function, which executes Filtered Back-Projection (FBP) reconstruction with a Ram-Lak -filter.The FBP was chosen since since it exposes the artifacts caused by interpolation flaws better than iterative methods like Algebraic Reconstruction Technique (ART), which have been considered as more suitable algorithms for the ill-posed situations (see for example Kolehmainen et al. (2003) or Siltanen et al. (2003)) [6,5].For post-processing only linear grayvalue mapping was used to guarantee fair comparison.The authors prefer Kak and Slaney (2001) and Natterer (1986) for further information about the reconstruction algorithms and Gonzalez and Woods (2008) for post-processing [11,12,42].The reconstruction results can be seen in Figure 2.

Discussion
The results in this study are very preliminary.However, both qualitative and quantitative results suggest that SINT is significantly better than general interpolation routines when the angular data sampling is extremely sparse.This can be observed from the metrics as well as comparing the subjective image quality in Figures 2-4.
See also the profiles in Figure 8.
In this study, we did not compare the SINT method against any sinogram-dedicated or state-of-the-art interpolation method for a number of reasons.Firstly, none of them are dedicated to as sparse a situation as our proposed method.For example, in the studies by Zhe et al. (2003) and Köstler et al. (2006) about 60% of the sinogram columns were interpolated, while in our study over 92% of the sinogram columns were interpolated [28,29].Similarly, in the study by Penßel et al. (2005) only 50% of the sinogram columns were interpolated [30].Secondly, part of the methods published recently are not fully documented or they require the right parameter settings for optimal result.Therefore, they cannot be used as reference since we could not guarantee a fair comparison.Thirdly, some of the new interpolation methods are designed for specific cases like the studies by Constantino (2006 and 2009) and earlier by Tam et al. (1990) [43][44][45].These studies were limited to a situation where only the outer boundaries of isolated objects were modeled, while we concentrated on a more generic sinogram interpolation.Finally, we considered that a comparison with generally known interpolation algorithms gives a better understanding of the capability of this method than comparing to specific and dedicated interpolations.
There are two reasons why generic interpolation routines fail when they are applied to sparse sinograms.First, the generic interpolation methods do not take advantage of the characteristics of the sinogram.Second, the generic interpolation methods work only when there is a strong correlation between the neighbor sinogram values.However, when the distance from the known column increases the correlation between sinogram elements decreases.This vanishing correlation generates a whirlpool-shaped pattern which can be seen in the lower row of Figures 2-4.To verify these findings, we calculated the absolute error as a function of distance from the known sinogram column, i.e.: where i is the distance from the known column to the interpolated column (with index i satisfying ( ) , where κ is the interpolation factor as defined in Section 2.3).The influence of distance to the interpolation error can be seen in Figure 9, which shows that when the distance from the known sinogram column increases, the interpolation error of the SINT method increases significantly slower than the interpolation error of reference interpolation methods.The disadvantage of the SINT method is execution time, which in our examples was about 15 s per interpolation.The most time consuming task was to verify that all warps intersect non-zero elements in the known columns as described in Section 2.2.Since every combination of the left and right neighbor column elements has to be considered as a warp candidate, the warp tracking operation has to be repeated  algorithm, which is a well-known iterative method for defining the eigenvalues without the sorting procedure as discussed by Golub and Van Loan (1996) [39].Since the matrix Λ was very sparse, positive-definite and orthogonal with only eigenvalues of one, the sorting procedure was irrelevant and only few iterations were needed.
There are number of articles indicating that parallel computing radically decreases the execution time in the image processing [46][47][48].Despite the fact that we used a non-optimized high-level interpreted language in this study, the SINT method could easily be coded for a parallel processing environment, such as clustered computers or general purpose graphical processing unit (GP-GPU), by using OpenCL, Cuda or other parallel computing coding techniques.The parallelization of the SINT method is effective since the interpolated sinogram values are determined only from the original known sinogram columns, which are constant during the interpolation process and therefore they can be stored in fast read-only memory before executing the interpolation routine.
In this study, the method has been proved to work in parallel beam imaging geometry cases in 2D, where the X-ray beam is perpendicular to the detector, instead of applying it to fan-beam in 2D or cone-beam in 3D.This limitation has been implemented partly for simplifying the theory and keeping the focus on the interpolation method itself.However, our interest is to expand this study to fan-beam and cone-beam imaging geometries in the near future.Our basic strategy for the parallel to fan beam conversion is based on modifying the warp matrix Φ (see Section 2.2).As illustrated by Kak and Slaney (Chapters 3.4 and 3.5) [11], the projection angle in fan beam imaging geometry is the sum of the projection angle and the fan angle, i.e.: ( ) where i γ is the fan angle, which depends on the sinogram element index.The fan angle can be calculated from the Equation where c is the shortest distance from X-ray source point to detector and q is defined by Equation (1.12).This means that the angular value of each column h is no longer constant, but depends on the row index of the sinogram element.Despite that this modification makes generation of the warp matrix Φ more complex, the rest of the algorithm can be utilized per se.Finally, as shown by Kak and Slaney (Chapter 3.6), cone-beam geometry can be considered as a stack of fan beam sinograms and the 3D volume can be built by the interpolation process in the azimuthal direction.

Conclusion
Finally, very promising implication of the SINT method is in metal artifact reduction.Since metals are typically very radio-dense, they block a major part of the X-rays causing streak artifacts as indicated by Goldman and Fowlkes (2000) and Veldkamp et al. (2010) [49,50].Since metallic areas can easily be seen in the sinogram as white regions, they can be identified, and the interpolation process could be used across these regions in the

Figure 1 .
Figure 1.Representations of two projection images and their relations in (a) spatial, (b) sinogram and (c) frequency domain.The distance a is the pixel size, ω is the difference of the frequency components, N is the number of pixels, Δθ 0 is the angular difference between the projection images and d is the maximum distance between two frequency components in adjacent columns.The relation between the spatial and the frequency domain is based on well-known Fourier Slice Theorem, which states that one projection image defines frequencies in the 2D frequency domain within a single line.This line intersects origin, it has same angle than the imaging angle and it is limited by the Nyquist frequency theorem.Books by Kak (2001) and Gonzalez (2008) provide detailed theoretical background of the Fourier Slice Theorem.

Figure 2 .
Figure 2. The FBP reconstruction of the Shepp-Logan-phantom with 5% Gaussian noise.Upper row: (a) The original phantom, (b) reconstruction from original data and (c) SINT interpolation.Lower row: FBP reconstruction when reference interpolation methods were used; (d) linear; (e) nearest neighbor and (f) spline interpolation.

Figure 3 .Figure 4 .
Figure 3.The FBP reconstruction of the Dental Arc-phantom with 5% Gaussian noise.Upper row: (a) The original phantom, (b) Reconstruction from original data and (c) SINT interpolation.Lower row: FBP reconstruction when reference interpolation methods were used; (d) Linear; (e) Nearest neighbor and (f) Spline interpolation.

Figure 5 .
Figure 5. Example of sinogram presentation.Left is the original object with two non-zero points.Middle is the sinogram presentation of the same object when nine images are taken from 20 to 180 degree angles.Each point generates a sine shaped wave so that the amplitude is proportional to the distance between the point and the center of the object (i.e.rotation axis).Right is the reconstruction done from the projection data by using back projection algorithm.

Figure 6 .
Figure 6.The original sinograms (top row) and their SINT interpolations (bottom row).Left sinogram from Shepp-Logan-phantom, middle from Dental Arc-phantom and right from Boxes-phantom.Lately, there have appeared several studies on dedicated directional sinogram interpolation algorithms.For example, Zhe et al. (2003) successfully combined the Huesman algorithm with Principle Component Analysis in PET imaging, and Köstler et al. (2006) applied PDE-based methods for interpolation and the Neumann boundary condition for defining missing columns in the sinogram [28,29].In this study, about 60% of the sinogram columns were interpolated.Similarly, in the study by Penbel et al. (2005), only 50% of the sinogram columns were interpolated by using cumulant functions instead of linear interpolation [30].Moreover, Bertram et al. (2004 and 2009) utilized tensors to gain directional interpolation based on sinogram data and achieved impressive results [31,32].Also, interesting approaches for sinogram interpolation in ill-posed situations are the usage of Stackgrams by Happonen (2005) and the usage of object geometry estimation for sinogram interpolation by Prince et al. (1993) [33,34].The advanced general directional interpolation methods are also considered by Gerchberg (1974), Papoulis (1975) as well as La Riviere and Pan (1999) [35-37].Finally, Dong et al. (2013) introduced very promising results when combining the sinogram interpolation with the wavelet cost function [38].For more general information about different sinogram related extrapolations and interpolations, see for example Brooks et al. (1978), Lahart (1981) or Kak and Slaney (2001) and references within [1,11,23].In this paper, we introduce a new approach, called Sinogram Interpolation Technique (SINT) for estimating new sinogram columns between known, measured sinogram columns in the extremely sparse situation of only nine known sinogram columns.This article is organized as follows: In Section 1.2, the theoretical background of the SINT method is described.Secondly, in the Section 1.3, the capability and limitations of SINT by three different numerical examples with additional noise are demonstrated.Finally, in the Section 1.4, we discuss the computational and clinical aspects of SINT and its possible future applications.Also, we analyze the failure of standard interpolation methods in extremely sparse imaging situations.Throughout this paper, capital letters are used for matrices.Moreover, for indicating a single element in a matrix, the standard [row][column]-order is applied.For example,

Figure 7 .
Figure 7. Outline of the warps.A single warp (thick line) connects sinogram elements

∈
 are called as the weight factors.The first upper index of weight factors indicates the known column index and the second indicates the sinogram element index in the estimated column h .The Equations (1.3) and (1.4) clearly indicate that if the weight factors are known, then an estimate for the

S
+ to a single sinogram element in the estimated column h S by two strictly positive weight factors.3) Each warp intersects only strictly positive sinogram elements in the whole sinogram domain.4) All non-zero sinogram elements are associated with at least one warp and all weight factors are associated with exactly one warp.

S
, where L is the number of warps, is to simplify the creation of other matrices.The warp matrix Φ defines all warps related to the estimated sinogram column h S such that each row represents a single warp as follows: The first element in the row is the row index 1 i of the known left point 1 , 1 i h S − , the second element is the row index 2 i of the estimated sinogram column 2 , i h S and finally the row index 3 i is the right known point 3 , 1 i h S + .To define the value 3 i the warp function ( ) k W θ is defined for all angles θ based on two points in the sinogram ( ( ) + ) and the their angular values ( 1 h θ − and 1 h θ + ).
. The matrix A includes grayvalues of the known sinogram columns 1 h S ± and Ψ consists of the corresponding row indexes.First we define the Relation Grayvalue Matrix A , which specifies the relations between the weight factor vector w and the sinogram elements in the estimated column h S .The Equation (1.6) can be represented in a matrix form as , unlike in Equation (1.6), the vector w in Equation (1.13) consists only the weight factors that are included in warps and therefore the size of the vector w is 2L .The vector 2 L w + ∈  will be called in this study as weight vector.


From the Equation (1.15) can be observed that the weight vector w is an positive eigenvector for Λ related to an eigenvalue of one 1 .The number of feasible eigenvectors is exactly L based on the fact that for each positive eigenvector, say ( )2 p Lv ∈  , there is always another eigenvector ( ) has at least one negative component and it is not therefore feasible Algorithm 2. Generation of the matrices A and Ψ.
, the matrix Λ is positively semi-definite and or- thogonal and it has (only) eigenvalues of one.
to determine the actual values of the warp factors.As defined in Equation (1.8), the sum of weight factors connected to any known sinogram element , 1 i h S ± equals to one.To fulfill this requirement, we generate a Rule Matrix R which sums all weight factors that are connected to the same known sinogram element.Then 1 Rw = (1.20)

Algorithm 4 .
to the same known non-zero sinogram element, otherwise , 0 j h R = .The order of the rows in the matrix R is irrelevant.The pseudo-code for generating the matrix R is described in Algorithm 4. Generation of the matrix R.

Algorithm 5 .
Generation of the Sinogram column S h .Constant ε is a small true positive number.See also Equation (1.25) for alternative way to calculate( ) angular difference in sinogram.Then from the Equation (1.27) we get sinogram generated from the original phantoms with similar angular difference than the interpolated sinograms.The results are shown in the Table

Figure 8 .
Figure 8. Profiles of the middle row in Shepp-Logan (top), Dental Arc (mid.) and Boxes (bottom) phantoms.The solid line is the profile of the original phantom, dashed is reconstruction with the SINT interpolation and dotted is the reconstruction with linear interpolation.Nearest neighbor and spline were not plotted because they produced very similar result than the linear interpolation.
whole interpolation process, where N is height of the sinogram, H is the number of original sinogram columns and κ is the interpolation factor.In this study, this meant over 4 million warp tracking operations.Another time consuming operation was the calculation of the eigenvalues for the matrix G .In this study, we applied the Hessenberg and Schur decomposition based QR

Figure 9 .
Figure 9.The SINT method (solid line), linear (long dash), spline (short dash) and nearest neighbor (dotted line) interpolation errors as a function of the distance from the known column for all three cases; Shepp-Logan (left), dental-arc (mid) and boxes (right).The error is calculated as described in the Equation (1.31).

Table 1 . Parameters used in this article.
∈  Relation grayvalue matrix, see Equation (1.14) d d ∈  maximum distance between adjacent columns, see Figure 1 G 2 L L G × + ∈  Eigenvector matrix, see Equation (1.18) h 1 h H ≤ ≤ index of the sinogram matrix column or angle vector factor vector, see Equation(1.18)