Local inverse

From HandWiki

The local inverse is a kind of inverse function or matrix inverse used in image and signal processing, as well as in other general areas of mathematics. The concept of a local inverse came from interior reconstruction of CT[clarification needed] images. One interior reconstruction method first approximately reconstructs the image outside the ROI (region of interest), and then subtracts the re-projection data of the image outside the ROI from the original projection data; then this data is used to make a new reconstruction. This idea can be widened to a full inverse. Instead of directly making an inverse, the unknowns outside of the local region can be inverted first. Recalculate the data from these unknowns (outside the local region), subtract this recalculated data from the original, and then take the inverse inside the local region using this newly produced data for the outside region.

This concept is a direct extension of the local tomography, generalized inverse and iterative refinement methods. It is used to solve the inverse problem with incomplete input data, similarly to local tomography. However this concept of local inverse can also be applied to complete input data.

Local inverse for full field of view system or over-determined system

[math]\displaystyle{ \begin{bmatrix} f \\ g \end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix}. }[/math]

Assume there are [math]\displaystyle{ E }[/math], [math]\displaystyle{ F }[/math], [math]\displaystyle{ G }[/math] and [math]\displaystyle{ H }[/math] that satisfy

[math]\displaystyle{ \begin{bmatrix} E & F \\ G & H \end{bmatrix} \begin{bmatrix} A & B \\ C & D \end{bmatrix} = J. }[/math]

Here [math]\displaystyle{ J }[/math] is not equal to [math]\displaystyle{ I }[/math], but [math]\displaystyle{ J }[/math] is close to [math]\displaystyle{ I }[/math], where [math]\displaystyle{ I }[/math] is the identity matrix. Examples of matrices of the type [math]\displaystyle{ \begin{bmatrix} E & F \\ G & H \end{bmatrix} }[/math] are the filtered back-projection method for image reconstruction and the inverse with regularization. In this case the following is an approximate solution:

[math]\displaystyle{ \begin{bmatrix} x_0 \\ y_0 \end{bmatrix} = \begin{bmatrix} E & F \\ G & H \end{bmatrix} \begin{bmatrix} f \\ g \end{bmatrix}. }[/math]

A better solution for [math]\displaystyle{ x_1 }[/math] can be found as follows:

[math]\displaystyle{ \begin{bmatrix} x_1 \\ y_1 \end{bmatrix} = \begin{bmatrix} E & F \\ G & H \end{bmatrix} \begin{bmatrix} f - B y_0\\ g - D y_0 \end{bmatrix}. }[/math]

In the above formula [math]\displaystyle{ y_1 }[/math] is useless, hence

[math]\displaystyle{ x_1=E(f-B y_0)+F(g - D y_0). }[/math]

In the same way, there is

[math]\displaystyle{ y_1=G(f-A x_0)+H(g - C x_0). }[/math]

In the above the solution is divided into two parts, [math]\displaystyle{ x }[/math] inside the ROI and [math]\displaystyle{ y }[/math] is outside the ROI, f inside of FOV(field of view) and g outside the FOV.

The two parts can be extended to many parts, in which case the extended method is referred to as the sub-region iterative refinement method [1]

Local inverse for limited-field-of-view system or under-determined system

[math]\displaystyle{ \begin{bmatrix} f \\ g \end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix}. }[/math]

Assume [math]\displaystyle{ A }[/math], [math]\displaystyle{ B }[/math], [math]\displaystyle{ C }[/math], and [math]\displaystyle{ D }[/math] are known matrices; [math]\displaystyle{ x }[/math] and [math]\displaystyle{ y }[/math] are unknown vectors; [math]\displaystyle{ f }[/math] is a known vector; [math]\displaystyle{ g }[/math] is an unknown vector. We are interested in determining x. What is a good solution?

Assume the inverse of the above matrix exists

[math]\displaystyle{ \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} E & F \\ G & H \end{bmatrix} =J . }[/math]

Here [math]\displaystyle{ J }[/math] is or is close to [math]\displaystyle{ I }[/math]. The local inverse algorithm is as follows:

(1) [math]\displaystyle{ g_{ex} }[/math]. An extrapolated version of [math]\displaystyle{ g }[/math] is obtained by

[math]\displaystyle{ g_{ex}|_{\partial G}=f|_{\partial F}. }[/math]

(2) [math]\displaystyle{ y_0 }[/math]. An approximate version of [math]\displaystyle{ y }[/math] is calculated by

[math]\displaystyle{ y_0=H g_{ex}. }[/math]

(3) [math]\displaystyle{ y' }[/math]. A correction for [math]\displaystyle{ y }[/math] is done by

[math]\displaystyle{ y'=y_0+y_\text{co}. }[/math]

(4) [math]\displaystyle{ f' }[/math]. A corrected function for [math]\displaystyle{ f }[/math] is calculated by

[math]\displaystyle{ f' = f-B y' . }[/math]

(5) [math]\displaystyle{ g_{1ex} }[/math]. An extrapolated function for [math]\displaystyle{ g }[/math] is obtained by

[math]\displaystyle{ g_{1ex}|_{\partial G}=f'|_{\partial F}. }[/math]

(6) [math]\displaystyle{ x_1 }[/math]. A local inverse solution is obtained

[math]\displaystyle{ x_1 = E f' +F g_{1ex}. }[/math]

In the above algorithm, there are two time extrapolations for [math]\displaystyle{ g }[/math] which are used to overcome the data truncation problem. There is a correction for [math]\displaystyle{ y }[/math]. This correction can be a constant correction to correct the DC values of [math]\displaystyle{ y }[/math] or a linear correction according to prior knowledge about [math]\displaystyle{ y }[/math]. This algorithm can be found in the following reference:.[2]

In the example of the reference,[3] it is found that [math]\displaystyle{ y'=y_0+y_\text{co}=k y_0 }[/math], here [math]\displaystyle{ k=1.04 }[/math] the constant correction is made. A more complicated correction can be made, for example a linear correction, which might achieve better results.

A^+ B is close to 0

Shuang-ren Zhao defined a Local inverse[2] to solve the above problem. First consider the simplest solution

[math]\displaystyle{ f = A x + B y, }[/math]

or

[math]\displaystyle{ A x = f - B y = f'. }[/math]

Here [math]\displaystyle{ f'=f - B y }[/math] is the correct data in which there is no influence of the outside object function. From this data it is easy to get the correct solution,

[math]\displaystyle{ x' = A^{-1} f'. }[/math]

Here [math]\displaystyle{ x' }[/math] is a correct(or exact) solution for the unknown [math]\displaystyle{ x }[/math], which means [math]\displaystyle{ x' = x }[/math]. In case that [math]\displaystyle{ A }[/math] is not a square matrix or has no inverse, the generalized inverse can applied,

[math]\displaystyle{ x' = A^{+} (f - B y)= A^{+} f'. }[/math]

Since [math]\displaystyle{ y }[/math] is unknown, if it is set to [math]\displaystyle{ 0 }[/math], an approximate solution is obtained.

[math]\displaystyle{ x_0 = A^{+} f. }[/math]

In the above solution the result [math]\displaystyle{ x_0 }[/math] is related to the unknown vector [math]\displaystyle{ y }[/math]. Since [math]\displaystyle{ y }[/math] can have any value the result [math]\displaystyle{ x_0 }[/math] has very strong artifacts, namely

[math]\displaystyle{ \mathrm{error}_0=|x_0-x'|=|A^{+} B y| }[/math].

These kind of artifacts are referred to as truncation artifacts in the field of CT image reconstruction. In order to minimize the above artifacts in the solution, a special matrix [math]\displaystyle{ Q }[/math] is considered, which satisfies

[math]\displaystyle{ QB = 0, }[/math]

and thus satisfies

[math]\displaystyle{ QA x = Qf - QB y = Qf. }[/math]

Solving the above equation with Generalized inverse gives

[math]\displaystyle{ x_1 =[QA]^{+} Qf = [A]^{+} Q^{+} Qf. }[/math]

Here [math]\displaystyle{ Q^{+} }[/math] is the generalized inverse of [math]\displaystyle{ Q }[/math], and [math]\displaystyle{ x_1 }[/math] is a solution for [math]\displaystyle{ x }[/math]. It is easy to find a matrix Q which satisfies [math]\displaystyle{ QB=0 }[/math], specifically [math]\displaystyle{ Q }[/math] can be written as the following:

[math]\displaystyle{ Q = I-BB^{+}. }[/math]

This matrix [math]\displaystyle{ Q }[/math] is referred as the transverse projection of [math]\displaystyle{ B }[/math], and [math]\displaystyle{ B^{+} }[/math] is the generalized inverse of [math]\displaystyle{ B }[/math]. The matrix [math]\displaystyle{ B^{+} }[/math] satisfies

[math]\displaystyle{ BB^{+}B=B, }[/math]

from which it follows that

[math]\displaystyle{ QB=[I-BB^{+}]B=B- BB^{+}B=B-B=0. }[/math]

It is easy to prove that [math]\displaystyle{ QQ=Q }[/math]:

[math]\displaystyle{ \begin{align} QQ & =[I-BB^+][I-BB^+] = I-2 BB^+ + BB^+ BB^+ \\ & = I-2 BB^+ + BB^+ = I-BB^+ = Q, \end{align} }[/math]

and hence

[math]\displaystyle{ QQQ=(QQ)Q=QQ=Q. }[/math]

Hence Q is also the generalized inverse of Q

That means

[math]\displaystyle{ Q^{+}Q=QQ=Q. }[/math]

Hence,

[math]\displaystyle{ x_1 =A^{+}[Q]^{+} Qf =A^{+} Q f }[/math]

or

[math]\displaystyle{ x_1 =[A]^{+} [I-BB^{+}] f. }[/math]

The matrix

[math]\displaystyle{ A^L=[A]^{+} [I-BB^{+}] }[/math]

is referred to as the local inverse of the matrix [math]\displaystyle{ \begin{bmatrix} A & B \\ C & D \end{bmatrix}. }[/math] Using the local inverse instead of the generalized inverse or the inverse can avoid artifacts from unknown input data. Considering,

[math]\displaystyle{ [A]^{+} [I-BB^{+}]f'=[A]^{+} [I-BB^{+}](f-B y) = [A]^{+} [I-BB^{+}]f. }[/math]

Hence there is

[math]\displaystyle{ x_1 =[A]^{+} [I-BB^{+}] f'. }[/math]

Hence [math]\displaystyle{ x_1 }[/math] is only related to the correct data [math]\displaystyle{ f' }[/math]. This kind error can be calculated as

[math]\displaystyle{ \mathrm{error}_1 = |x_1-x'| =| [A]^{+}BB^{+} f' |. }[/math]

This kind error are called the bowl effect. The bowl effect is not related to the unknown object [math]\displaystyle{ y }[/math], it is only related to the correct data [math]\displaystyle{ f' }[/math].

In case the contribution of [math]\displaystyle{ [A]^{+}BB^{+}f' }[/math] to [math]\displaystyle{ x }[/math] is smaller than that of [math]\displaystyle{ [A]^{+} B y }[/math], or

[math]\displaystyle{ \mathrm{error}_1\ll \mathrm{error}_0, }[/math]

the local inverse solution [math]\displaystyle{ x_1 }[/math] is better than [math]\displaystyle{ x_0 }[/math] for this kind of inverse problem. Using [math]\displaystyle{ x_1 }[/math] instead of [math]\displaystyle{ x_0 }[/math], the truncation artifacts are replaced by the bowl effect. This result is the same as in local tomography, hence the local inverse is a direct extension of the concept of local tomography.

It is well known that the solution of the generalized inverse is a minimal L2 norm method. From the above derivation it is clear that the solution of the local inverse is a minimal L2 norm method subject to the condition that the influence of the unknown object [math]\displaystyle{ y }[/math] is [math]\displaystyle{ 0 }[/math]. Hence the local inverse is also a direct extension of the concept of the generalized inverse.

See also

References

  1. Shuangren Zhao, Xintie Yang, Iterative reconstruction in all sub-regions, SCIENCEPAPER ONLINE. 2006; 1(4): page 301–308, http://www.paper.edu.cn/uploads/journal/2007/42/1673-7180(2006)04-0301-08.pdf
  2. 2.0 2.1 Shuangren Zhao, Kang Yang, Dazong Jiang, Xintie Yang, Interior reconstruction using local inverse, J Xray Sci Technol. 2011; 19(1): 69–90
  3. S. Zhao, D Jaffray, Iterative reconstruction and reprojection for truncated projections, AAPM 2004, Abstract in Medical Physics 2004, Volume 31, P1719, http://imrecons.com/wp-content/uploads/2013/02/iterative_extro.pdf