###### tags: `BioI-maging`

# Algebraic Reconstruction

## 一 Motivation: Fully discretize the problem.
Abstract the object into a matrix (2-D case). Each entry is represented by a pixel. Each pixel(entry) is set to a same default value (the first guess).
#### Question 1: What about 3-D cases? How can we store the information of the object?
3-D Matrix: tensor.
$a_{i,j}$ is the path length of the $i^{th}$ beam array in the pixel $j$. $u_j$ is the pixel value of the $j^{th}$ pixel grid.
Consider the radon transform:
$\int_{L_i} f(x) ds = - log \frac{I_i} {I_0} = b_i$
where $L_i$ is a line through the obejct. $I_0$ is the starting point of the beam. $I_i$ is the ending point of the beam. $I_i$ is defined by the line $L_i$. Beer-Lambert law tells that the line integral will be $- log \frac{I_i} {I_0}$ and we call it $b_i$.

We can approximate the same line integral using $a_{i,j}$ and $u_j$:
$$\sum_j a_{ij} u_j = b_i$$
In a large set of linear equation, this is:
$Au = b$
This is now a linear system question.

#### Question 2: How do we get matrix A in practice?
Sparse representation
kronecker product
Tensor algebra <> matrix algebra <> vector algebra
learn the difference between inner product , outer product, kronecker product, Cramér–Rao product.
They can be reduced to matrix operations.

However,
1. in reality, we often use non-square matrix to represent the object. Thus, the matrix is not invertible.
2. When the matrix goes big, it's expensive to solve the problem.

Solution for 1: Least Square!
Solution for 2: Interative Approach.

## 贰 Interative Approach
($i$) Motivation: We make a guess and compare it with the actual result, and we minimize the error. We can consider their difference as a function. Then we need to find the min value in an interval. This becomes a Gradient Descent problem.

($ii$) The Landweber Algorithm:
Iterate step:
$$u^{(k+1)} = u^{k} + t_k A^T (b - A u^{(k)})$$, where $k \in \Bbb N$

Remark:
1. This the gradient descent for the least square problem.
2. $t_k$ is the step size. It can be small so that our updated estimate won't go so far that it misses the $min$.

Interpretation:
1. $A u^{(k)}$ : Project current estimate $u^{(k)}$ to sinogram.
2. $b - A u^{(k)}$ : Substract from b (measured data) to obtain residual.
3. $A^T (b - A u^{(k)})$ backproject residual.
4. $u^{k} + t_k A^T (b - A u^{(k)})$ : update current estimate.

Warning: Least Square is not good for noisy data. It can magnify the noise. The solution is the iterate certain times and stop before the reconstructed object model becomes too noisy. At this point, the reconstructed model will have a general shape of the object and not too much noise, so it is a good time to stop.

The problem of the Landweber algorithm: to get the step size, we need to compute the largest eigenvalue of the matrix $A^TA$. This is too expensive. Thus, in practice, the Landweber algorithm is rarely used.

Other algorithm: ART, SIRT, CGLS, ...

($iii$) SIRT algorithm: try to overcome this problem.
$$u^{(k+1)} = u^{k} + t_k \color{red}{C} A^T \color{red}{R} (b - A u^{(k)})$$, where $k \in \Bbb N$
and for $A \in \Bbb R^{n \times n}$ ,
$C \in \Bbb R^{n \times n}$, and C is a diagonal matrix with $C_{j,j} = {1 \over \sum_i a_ij}$ .
$R \in \Bbb R^{m \times m}$, and R is a diagonal matrix with $R_{i,i} = {1 \over \sum_j a_ij}$ .

Why SIRT?
1. Faster than Landweber.
2. Easy to find $t_k$ since now, $0 < t_k < 2$.

Remark: Notice that both Landweber and SIRT only look at an iteration (projection result) back, not all the iterations we already have.

($iv$) CGLS: Krylov subspace method
Rro: Use information of all the previous iterates.
Con: Can not consider some special cosntraints e.g. non-negativity of the values.

($v$) Conclusion:
Pros and Cons of algebraic reconstruction methods:
Good:
1. Can be used in special geometry systems.
2. Can incorporate special constraints.

Bad:
1. More expensive than FBP(Filtered Back Projection)
2. When to stop: see ($ii$)Landweber (Warning).

## 叁 Regularized Method
Regularization: A way of enforcing additional information to obtain a nicer (more regular) solution. This helps us to handle the noise in the data.

Recall: In algebraic method, we sometimes stop the iteration early.
This is an example of enforcing regularization.

TV Regularization:
$$\underset{u}Min ||(b - A u||^2 + \lambda \cdot TV(u)$$
where $u$ is a particular object reconstruction, b is still the projection we obtained, $\lambda$ is something we try out. $TV(u)$ measure the overall variation of the data (the reconstruction $u$ we already have). When $TV(u)$ is small, the image (object) will be very smooth since data variation is limited, but the image can be blur at same time. So the prupose is to try out a good $\lambda$ that gives $TV(u)$ a appropreiate significance so that the reconstruction $u$ has the best balance between its correctness (how close is the projection of our current guess and the real projection data we get i.e. $b$) and smoothness (less noise).

##### The source till this point is
1. https://www.coursera.org/learn/cinemaxe/lecture/Puytk/iterative-reconstruction-jakob-sauer-jorgensen
2. https://www.coursera.org/learn/cinemaxe/lecture/M6vyB/regularization-jakob-sauer-jorgensen
3. Clinical Brain Imaging: Computerized Axial Tomography and Magnetic Resonance Imaging by William W. Orrison, Jr. John A. Sanders

## 肆 $L^2$-gradient flow
It seems this algorithm is an iterative regularization method.
#### Question 3: is it?
Answer: Yes. "In this paper we propose an $L^2$-gradient flow three-dimensional reconstruction algorithm of single-particle analysis for solving a variational model with TV regularization term." (Single-particle reconstruction using L2-gradient flow by Ming Li et al.)

This method is essentially a iterative method with regularization featuring special way of finding the function $f$ (denoted by $u$ in the previous sections).

Motivation: $f$ can be written as
$$f(x,y,z) = \sum_{ijk} f_{ijk}N_i^3(x)N_j^3(y)N_k^3(z)$$
where $N^3_a (a = i,j,k)$ is a cubic B-spline basis function and $f_{ijk}$ is the coefficient of tri-cubic B-spline tensor product $N_i^3(x), N_j^3(y), N_k^3(z)$. (Sort of like a decomposition).

With Schmidt process (Gram-Schmidt process?), we can obtain orthogonal basis. This makes the later computations even easier so that we don't solve big system of linear equations.
$$f(x,y,z) = \sum_{ijk} \tilde f_{ijk}\tilde N_i^3(x)\tilde N_j^3(y)\tilde N_k^3(z)$$
where $\tilde N^3_a (a = i,j,k)$ is an orthogonal basis of cubic B-spline basis function and $\tilde f_{ijk}$ is the coefficient of tri-cubic B-spline tensor product $\tilde N_i^3(x), \tilde N_j^3(y), \tilde N_k^3(z)$.

In this way, we can change $f$ to $\tilde f_{ijk}$, and $\tilde f_{ijk}$ back to $f$. An iteration scheme is specified to find $\tilde f_{ijk}^{(m+1)}$ according to $\tilde f_{ijk}^{(m)}$ (though I don't really understand what's it doing).

The algorithm goes like this (from the paper):

The reason to obtain $f^{(m+1)}$ everytime we have $\tilde f_{ijk}^{(m+1)}$ is we need $f^{(m+1)}$ in the iteration scheme for $\tilde f_{ijk}^{(m+2)}$.

##### Result: L2GF algorithm can obtain a resolution better than that of WBP, Fourier method, block ART and SIRT.

##### Source: Single-particle reconstruction using $L^2$-gradient flow by Ming Li et al.

Remark:
To understand the detail,
First learn discrete case, then continuous case.
L2GF flow.
L2GF. Mark the paper don't undertand