Interior reconstruction

1) Projections of Sheep-Logan phantoms 2) Truncated projections (zero extrapolation) 3) Constant, 4) exponential and 5) quadratical extrapolations 6) Mixed extrapolation of 4 and 5

${\displaystyle {\begin{bmatrix}f\\g\end{bmatrix}}={\begin{bmatrix}A&B\\C&D\end{bmatrix}}{\begin{bmatrix}x\\y\end{bmatrix}}}$

${\displaystyle A}$, ${\displaystyle B}$, ${\displaystyle C}$, ${\displaystyle D}$ are known matrices; ${\displaystyle x}$ and ${\displaystyle y}$ are unknown vectors; ${\displaystyle f}$ is a known vector, and ${\displaystyle g}$ is an unknown vector. We need to know the vector ${\displaystyle x}$. ${\displaystyle x}$ and ${\displaystyle y}$ are the original image, while ${\displaystyle f}$ and ${\displaystyle g}$ are measurements of responses. Vector ${\displaystyle x}$ is inside the region of interest ${\displaystyle X}$, (${\displaystyle x\in X}$). Vector ${\displaystyle y}$ is outside the region ${\displaystyle X}$. The outside region is called ${\displaystyle Y}$, (${\displaystyle y\in Y}$) and ${\displaystyle f}$ is inside a region in the measurement corresponding to ${\displaystyle X}$. This region is denoted ${\displaystyle F}$, (${\displaystyle f\in F}$). The region of vector ${\displaystyle g}$ (outside the region ${\displaystyle F}$) also corresponds to ${\displaystyle Y}$ and is denoted as ${\displaystyle G}$, (${\displaystyle g\in G}$). In CT image reconstruction, it has

${\displaystyle C=0}$

To simplify the concept of interior reconstruction, the matrices ${\displaystyle A}$, ${\displaystyle B}$, ${\displaystyle C}$, ${\displaystyle D}$ are applied to image reconstruction instead of a complex operator.

The response in the outside region can be a guess ${\displaystyle G}$; for example, assume it is ${\displaystyle g_{ex}}$

${\displaystyle {\begin{bmatrix}x_{0}\\y_{0}\end{bmatrix}}={\begin{bmatrix}A&B\\C&D\end{bmatrix}}^{-1}{\begin{bmatrix}f\\g_{ex}\end{bmatrix}}}$

a) Shepp-Logan head phantom b) Crop of the phantom c) Reconstruction without extrapolation d) Reconstruction with constant, (e) quadratic and (f) mixed extrapolation

A solution of ${\displaystyle x}$ is written as ${\displaystyle x_{0}}$, and is known as the extrapolation method. The result depends on how good the extrapolation function ${\displaystyle g_{ex}}$ is. A frequent choice is

${\displaystyle g_{ex}|_{\partial G}=f|_{\partial F}}$

at the boundary of the two regions.[1][2][3][4] The extrapolation method is often combined with a priori knowledge,[5][6] and an extrapolation method which reduces calculation time is shown below.

Assume a rough solution, ${\displaystyle x_{0}}$ and ${\displaystyle y_{0}}$, is obtained from the extrapolation method described above. The response in the outside region ${\displaystyle g_{1}}$ can be calculated as follows:

${\displaystyle g_{1}=Cx_{0}+Dy_{0}}$

The reconstructed image can be calculated as follows:

${\displaystyle {\begin{bmatrix}x_{1}\\y_{1}\end{bmatrix}}={\begin{bmatrix}A&B\\C&D\end{bmatrix}}^{-1}{\begin{bmatrix}f\\g_{1}+g_{1ex}\end{bmatrix}}}$

It is assumed that

${\displaystyle f|_{\partial F}=(g_{1}+g_{1ex})|_{\partial G}}$

at the boundary of the interior region; ${\displaystyle x_{1}}$ solves the problem, and is known as the adaptive extrapolation method. ${\displaystyle g_{1ex}}$ is the adaptive extrapolation function.[7][8][9][10][5]

It is assumed that a rough solution, ${\displaystyle x_{0}}$ and ${\displaystyle y_{0}}$, is obtained from the extrapolation method described below:

${\displaystyle {\begin{bmatrix}f_{1}\\g_{1}\end{bmatrix}}={\begin{bmatrix}A&B\\C&D\end{bmatrix}}{\begin{bmatrix}0\\y_{0}\end{bmatrix}}}$

or

${\displaystyle f_{1}=By_{0}}$

The reconstruction can be obtained as

${\displaystyle {\begin{bmatrix}x_{1}\\y_{1}\end{bmatrix}}={\begin{bmatrix}A&B\\C&D\end{bmatrix}}^{-1}{\begin{bmatrix}f-f_{1}\\g_{ex}\end{bmatrix}}}$

Here ${\displaystyle g_{1ex}}$ is an extrapolation function, and it is assumed that

${\displaystyle (f-f_{1})|_{\partial F}=g_{1ex}|_{\partial G}}$

${\displaystyle x_{1}}$ is one solution of this problem.[11]

Local tomography, with a very short filter, is also known as lambda tomography.[12][13]

The local inverse method extends the concept of local tomography. The response in the outside region can be calculated as follows:

${\displaystyle f=Ax+By}$

Consider the generalized inverse ${\displaystyle B^{+}}$ satisfying

${\displaystyle BB^{+}B=B}$

Define

${\displaystyle Q=[I-BB^{+}]}$

so that

${\displaystyle QB=0}$

Hence,

${\displaystyle Qf=QAx}$

The above equation can be solved as

${\displaystyle x_{1}=A^{+}Q^{+}Qf}$,

considering that

${\displaystyle QQ=Q}$

${\displaystyle QQQ=Q}$

${\displaystyle Q}$ is the generalized inverse of ${\displaystyle Q}$, i.e.

${\displaystyle Q^{+}=Q}$

The solution can be simplified as

${\displaystyle x_{1}=A^{+}Qf}$.

The matrix ${\displaystyle A^{+}Q=A^{+}[I-BB^{+}]}$ is known as the local inverse of matrix ${\displaystyle {\begin{bmatrix}A&B\\C&D\\\end{bmatrix}}}$, corresponding to ${\displaystyle A}$. This is known as the local inverse method.[11]

Here a goal function is defined, and this method iteratively achieves the goal. If the goal function can be some kind of normal, this is known as the minimal norm method.

${\displaystyle \min(R\|x\|+S\|y\|+T\|g\|)}$,

subject to

${\displaystyle {\begin{bmatrix}x\\y\end{bmatrix}}={\begin{bmatrix}A&B\\C&D\end{bmatrix}}^{-1}{\begin{bmatrix}f\\g\end{bmatrix}}}$

and ${\displaystyle f}$ is known,

where ${\displaystyle R}$, ${\displaystyle S}$ and ${\displaystyle T}$ are weighting constants of the minimization and ${\displaystyle \|\cdot \|}$ is some kind of norm. Often-used norms are ${\displaystyle L_{0}}$, ${\displaystyle L_{1}}$, ${\displaystyle L_{2}}$, ${\displaystyle L_{+\infty }}$ total variation (TV) norm or a combination of the above norms. An example of this method is the projection onto convex sets (POCS) method.[14] [15]

In special situations, the interior reconstruction can be obtained as an analytical solution; the solution of ${\displaystyle x}$ is exact in such cases.[16][17][18]

Extrapolated data often convolutes to a kernel function. After data is extrapolated its size is increased N times, where N = 2 ~ 3. If the data needs to be convoluted to a known kernel function, the numerical calculations will increase log(N)·N times, even with the fast Fourier transform (FFT). An algorithm exists, analytically calculating the contribution from part of the extrapolated data. The calculation time can be omitted, compared to the original convolution calculation; with this algorithm, the calculation of a convolution using the extrapolated data is not noticeably increased. This is known as fast extrapolation.[19]