A multi-scale geometric flow method for molecular structure reconstruction

Given an even integer $m\geqslant 4$, suppose the domain ${{\Omega }_{c}}={{\left[ -\frac{n}{2}-1,\frac{n}{2}+1 \right]}^{3}}$ is uniformly partitioned with grid points

$$\begin{eqnarray*}&&{{{\bf x}}_{{\bf i}}}={\bf i}h,\ -\frac{m}{2}\leqslant {\bf i}\leqslant \frac{m}{2},\ {\bf i}=({{i}_{1}},{{i}_{2}},{{i}_{3}}),\end{eqnarray*}$$

where $h=\frac{n+2}{m}$. The function f is represented as

$$\begin{eqnarray}&&f({\bf x})=\mathop{\sum }\limits_{-\frac{m}{2}+2\leqslant {\bf i}\leqslant \frac{m}{2}-2}{{f}_{{\bf i}}}{{\phi }_{{\bf i}}}({\bf x},h),\ \ \ {\bf x}={{[{{x}_{1}},{{x}_{2}},{{x}_{3}}]}^{T}}\in {{\Omega }_{c}},\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray*}&&{{\phi }_{{\bf i}}}({\bf x},h)=N(\parallel {\bf x}-{\bf i}h\parallel /h),\end{eqnarray*}$$

and N(s) is a cubic B-spline basis function. The cubic B-spline basis function, defined on the uniform knots $-2,-1,0,1,2$, is

$$\begin{eqnarray}N(s)=\left\{ \begin{array}{ccccccccccccccc} \frac{2}{3}-{{s}^{2}}+\frac{|s{{|}^{3}}}{2}, & 0\leqslant |s|\lt 1, \\ \frac{1}{6}{{\left( 2-|s| \right)}^{3}}, & 1\leqslant |s|\lt 2, \\ 0, & 2\leqslant |s|. \\ \end{array} \right.\end{eqnarray} \tag{ 9 }$$

The support of $N(s/h)$ is the interval $(-2h,2h)$. Hence, the support of ${{\phi }_{{\bf i}}}({\bf x},h)$ is $\{{\bf x}\in {{\mathbb{R}}^{3}}:\parallel {\bf x}-{\bf i}h\parallel \lt 2h\}$. It is not difficult to show that

Theorem 4.1. The functions in the set ${{\{{{\phi }_{{\bf i}}}\}}_{-\frac{m}{2}+2\leqslant {\bf i}\leqslant \frac{m}{2}-2}}$ are linearly independent.

Therefore, they form a basis functional space. Since N(s) is a C²-continuous function, ${{\phi }_{{\bf i}}}({\bf x},h)$ is also C² continuous and

$$\begin{eqnarray*}\begin{array}{rcl} \frac{\partial {{\phi }_{{\bf i}}}({\bf x},h)}{\partial {{x}_{j}}} & = & \left\{ \begin{array}{ccccccccccccccc} \frac{{{y}_{j}}N^{\prime} (\parallel {\bf y}\parallel /h)}{h\parallel {\bf y}\parallel },\ & \parallel {\bf y}\parallel \ne 0, \\ 0,\ & \parallel {\bf y}\parallel =0, \\ \end{array} \right. \\ j & = & 1,2,3, \\ \frac{{{\partial }^{2}}{{\phi }_{{\bf i}}}({\bf x},h)}{\partial {{x}_{j}}\partial {{x}_{k}}} & = & \left\{ \begin{array}{ccccccccccccccc} \frac{{{y}_{j}}{{y}_{k}}N^{\prime\prime} (\parallel {\bf y}\parallel /h)}{{{h}^{2}}\parallel {\bf y}{{\parallel }^{2}}}-\frac{{{y}_{j}}{{y}_{k}}N^{\prime} (\parallel {\bf y}\parallel /h)}{h\parallel {\bf y}{{\parallel }^{3}}},\ & \parallel {\bf y}\parallel \ne 0, \\ 0,\ & \parallel {\bf y}\parallel =0, \\ \end{array} \right. \\ j,k & = & 1,2,3,\ j\ne k, \\ \frac{{{\partial }^{2}}{{\phi }_{{\bf i}}}({\bf x},h)}{\partial x_{j}^{2}} & = & \left\{ \begin{array}{ccccccccccccccc} \frac{y_{j}^{2}N^{\prime\prime} (\parallel {\bf y}\parallel /h)}{{{h}^{2}}\parallel {\bf y}{{\parallel }^{2}}}-\frac{(y_{j}^{2}-\parallel {\bf y}{{\parallel }^{2}})N^{\prime} (\parallel {\bf y}\parallel /h)}{h\parallel {\bf y}{{\parallel }^{3}}}, & \parallel {\bf y}\parallel \ne 0, \\ -\frac{2}{{{h}^{2}}}, & \parallel {\bf y}\parallel =0, \\ \end{array} \right. \\ j & = & 1,2,3, \\ \end{array}\end{eqnarray*}$$

where ${\bf y}={\bf x}-{\bf i}h={{[{{y}_{1}},{{y}_{2}},{{y}_{3}}]}^{T}}$.

The volumetric electron density maps of molecules are often approximated by Gaussian maps in the literature (see [2, 10, 26]). In such approximations, each atom is simulated by a sphere. Using radial basis functions, the spherical property of atoms can be ideally approximated. The left figure of figure 1 shows the iso-contours of a bi-cubic B-spline basis function. When the iso-value approaches 0, the iso-contours differ greatly from the circles. The right figure shows the iso-contours of the cubic B-spline radial basis function. Furthermore, the projections of the cubic B-spline basis functions have no closed form, while the projections of the cubic radial B-spline basis functions can be exactly evaluated from their closed forms. This increases greatly the computational efficiency.

Zoom In Zoom Out Reset image size

In our bi-gradient method, we need to compute the partial derivatives of J(f) with respect to the coefficients of f. Hence, we consider first the computation of these partial derivatives. It is easy to see that

$$\begin{eqnarray}&&\frac{\partial (J(f))}{\partial {{f}_{{\bf i}}}}=\frac{2}{p}\mathop{\sum }\limits_{l=1}^{p}\int _{-\frac{n}{2}}^{\frac{n}{2}}\int _{-\frac{n}{2}}^{\frac{n}{2}}[{{X}_{{{{\bf d}}_{l}}}}f-{{g}_{{\bf d}}}]{{X}_{{{{\bf d}}_{l}}}}{{\phi }_{{\bf i}}}{\rm d}u{\rm d}v.\end{eqnarray} \tag{ 10 }$$

Now let us explain how each of the terms in (10) is efficiently computed.

Computing ${{X}_{{{{\bf d}}_{l}}}}{{\phi }_{{\bf i}}}$.

Let ${\bf d}\in {{S}^{2}}$ be a given direction. Then the projection of ${{\phi }_{{\bf i}}}$ in the direction ${\bf d}$ is a 2D function, defined as

$$\begin{eqnarray*}\begin{array}{rcl} ({{X}_{{\bf d}}}{{\phi }_{{\bf i}}})(u,v) & = & \int _{-\infty }^{\infty }{{\phi }_{{\bf i}}}(u{\bf e}_{{\bf d}}^{(1)}+v{\bf e}_{{\bf d}}^{(2)}+t{\bf d},h)\ {\rm d}t \\ {} & = & \int _{-\infty }^{\infty }N(\parallel u{\bf e}_{{\bf d}}^{(1)}+v{\bf e}_{{\bf d}}^{(2)}+t{\bf d}-{\bf i}h\parallel /h)\ {\rm d}t \\ {} & = & \int _{-\infty }^{\infty }N\left( \left\|\left( u/h-{{{\bf i}}^{T}}{\bf e}_{{\bf d}}^{(1)} \right){\bf e}_{{\bf d}}^{(1)}+\left( v/h-{{{\bf i}}^{T}}{\bf e}_{{\bf d}}^{(2)} \right){\bf e}_{{\bf d}}^{(2)}+t{\bf d}/h\right \| \right)\ {\rm d}t \\ {} & = & h\int _{-\infty }^{\infty }N\left( \left\|{{{\bf a}}_{{\bf i}}}(u,v)+t{\bf d}\right \| \right)\ {\rm d}t \\ {} & = & h\int _{-\infty }^{\infty }N\left( \sqrt{{{a}^{2}}+{{t}^{2}}} \right)\ {\rm d}t. \\ \end{array}\end{eqnarray*}$$

For the cubic B-spline radial basis function ${{\phi }_{i}}$, we have

$$\begin{eqnarray*}&&\begin{array}{l} ({{X}_{{\bf d}}}{{\phi }_{{\bf i}}})(u,v) \\ \qquad =\left\{ \begin{array}{ccccccccccccccc} 2h\int _{0}^{\sqrt{1-{{a}^{2}}}}N\left( \sqrt{{{a}^{2}}+{{t}^{2}}} \right)\ {\rm d}t+2h\int _{\sqrt{1-{{a}^{2}}}}^{\sqrt{4-{{a}^{2}}}}N\left( \sqrt{{{a}^{2}}+{{t}^{2}}} \right)\ {\rm d}t, & 0\leqslant a\leqslant 1, \\ 2h\int _{0}^{\sqrt{4-{{a}^{2}}}}N\left( \sqrt{{{a}^{2}}+{{t}^{2}}} \right)\ {\rm d}t, & 1\lt a\lt 2, \\ 0, & 2\leqslant a\lt \infty , \\ \end{array} \right. \\ \end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl} a & = & \parallel {{{\bf a}}_{{\bf i}}}(u,v)\parallel , \\ {{{\bf a}}_{{\bf i}}}(u,v) & = & \left( u/h-{{{\bf i}}^{T}}{\bf e}_{{\bf d}}^{(1)} \right){\bf e}_{{\bf d}}^{(1)}+\left( v/h-{{{\bf i}}^{T}}{\bf e}_{{\bf d}}^{(2)} \right){\bf e}_{{\bf d}}^{(2)}, \\ \end{array}\end{eqnarray*}$$

${\bf e}_{{\bf d}}^{(1)}$ and ${\bf e}_{{\bf d}}^{(2)}$ are two directions defined by (1), which span the (u,v)-plane in space ${{\mathbb{R}}^{3}}$. The integrations above can be exactly computed using the expression (9) and the following integral formulas

$$\begin{eqnarray*}\begin{array}{rcl} \int {{({{t}^{2}}+{{a}^{2}})}^{\frac{1}{2}}}{\rm d}t & = & \frac{1}{2}\left[ t{{({{t}^{2}}+{{a}^{2}})}^{\frac{1}{2}}}+{{a}^{2}}{\rm log} \left( t+{{({{t}^{2}}+{{a}^{2}})}^{\frac{1}{2}}} \right) \right]+C, \\ \int {{({{t}^{2}}+{{a}^{2}})}^{\frac{3}{2}}}{\rm d}t & = & \frac{1}{16}\left[ 2t(5{{a}^{2}}+2{{t}^{2}}){{({{t}^{2}}+{{a}^{2}})}^{\frac{1}{2}}}+8{{a}^{4}}{\rm log} \left( t+{{({{t}^{2}}+{{a}^{2}})}^{\frac{1}{2}}} \right) \right. \\ {} & {} & -\left. {{a}^{4}}{\rm log} \left( {{t}^{2}}+\frac{{{a}^{2}}}{2}+t{{({{a}^{2}}+{{t}^{2}})}^{\frac{1}{2}}} \right) \right]+C. \\ \end{array}\end{eqnarray*}$$

Computing ${{X}_{{{{\bf d}}_{l}}}}f$. Here we propose an efficient approach for computing ${{X}_{{{{\bf d}}_{l}}}}f$.

$$\begin{eqnarray*}\begin{array}{rcl} ({{X}_{{{{\bf d}}_{l}}}}f)(u,v) & = & {{X}_{{{{\bf d}}_{l}}}}\left( \mathop{\sum }\limits_{{\bf i}}{{f}_{{\bf i}}}{{\phi }_{{\bf i}}} \right)(u,v) \\ {} & = & \mathop{\sum }\limits_{{\bf i}}{{f}_{{\bf i}}}({{X}_{{{{\bf d}}_{l}}}}{{\phi }_{{\bf i}}})(u,v), \\ \end{array}\end{eqnarray*}$$

where $({{X}_{{{{\bf d}}_{l}}}}{{\phi }_{{\bf i}}})(u,v)$ is computed as above. Notice that N(s) is locally supported. The cost for computing ${{X}_{{{{\bf d}}_{l}}}}f$ is $O({{m}^{3}})$. The total cost for the projection is $O(p{{m}^{3}})$, where p denotes the total number of projections. Thus, compared with using fast Fourier transform, the cost of this method is one order higher; however, its performance is much better. The computation can be accelerated by removing small coefficients.

$$\begin{eqnarray*}&&({{X}_{{{{\bf d}}_{l}}}}f)(u,v)\approx \mathop{\sum }\limits_{\{{{f}_{{\bf i}}}\gt \epsilon \}}{{f}_{{\bf i}}}({{X}_{{{{\bf d}}_{l}}}}{{\phi }_{{\bf i}}})(u,v),\end{eqnarray*}$$

where $\epsilon \gt 0$ is a given small number.

Let

$$\begin{eqnarray}&&f({\bf x})=\mathop{\sum }\limits_{-\frac{m}{2}+2\leqslant {\bf i}\leqslant \frac{m}{2}-2}{{f}_{{\bf i}}}{{\phi }_{{\bf i}}}({\bf x},h),\ \ \ {\bf x}={{[{{x}_{1}},{{x}_{2}},{{x}_{3}}]}^{T}}\in {{\Omega }_{c}}.\end{eqnarray} \tag{ 11 }$$

In the following, the vector consisting of the coefficient of $f({\bf x})$ is denoted as ${\bf f}$. Now we refine $f({\bf x})$ by replacing h and m with $\frac{h}{2}$ and $2m$, respectively. To obtain a representation in the form

$$\begin{eqnarray*}&&f({\bf x})=\mathop{\sum }\limits_{-m+2\leqslant {\bf i}\leqslant m-2}{{f}_{{\bf i}}}{{\phi }_{{\bf i}}}({\bf x},h/2),\ \ \ {\bf x}={{[{{x}_{1}},{{x}_{2}},{{x}_{3}}]}^{T}}\in {{\Omega }_{c}},\end{eqnarray*}$$

from (11), we need to refine formulas for $N(s,h)$. It is easy to derive that

$$\begin{eqnarray*}\begin{array}{rcl} N(s,h) & = & \frac{1}{8}N(s-h,\frac{h}{2})+\frac{1}{2}N(s-\frac{h}{2},\frac{h}{2}) \\ {} & {} & +\frac{3}{4}N(s,\frac{h}{2})+\frac{1}{2}N(s+\frac{h}{2},\frac{h}{2})+\frac{1}{8}N(s+h,\frac{h}{2}). \\ \end{array}\end{eqnarray*}$$

Hence

$$\begin{eqnarray*}&&{{\phi }_{{\bf 0}}}({\bf x},h)=N(\parallel {\bf x}\parallel ,\;h)\approx \mathop{\sum }\limits_{-2\leqslant {\bf i}\leqslant 2}{{n}_{{\bf i}}}{{\phi }_{{\bf i}}}({\bf x},h/2),\end{eqnarray*}$$

where values for ${{n}_{{\bf i}}}$ are computed by interpolation. There are a total of 125 basis functions involved, and the coefficients could be determined by interpolating 125 points ${{[-h,-h/2,0,h/2,h]}^{3}}$.

One method to minimize the energy J(f) is to use the following L2GF

$$\begin{eqnarray*}&&{{\int }_{{{\mathbb{R}}^{3}}}}\frac{\partial f}{\partial t}{{\phi }_{{\bf i}}}\ {\rm d}{\bf x}=-\frac{\partial (J(f))}{\partial {{f}_{{\bf i}}}},\ -\frac{m}{2}\leqslant {\bf i}\leqslant \frac{m}{2}.\end{eqnarray*}$$

In matrix form, the flow can be written as

$$\begin{eqnarray*}&&M\frac{\partial {\bf f}}{\partial t}=-\nabla J(f),\end{eqnarray*}$$

where $M=\left[ {{\int }_{{{\mathbb{R}}^{3}}}}{{\phi }_{{\bf i}}}{{\phi }_{{\bf j}}}{\rm d}{\bf x} \right]$ is a sparse matrix and $\nabla J(f)=R\ {\bf f}-G$ with R and G are defined by the right-hand side of (10). We call the vector ${{M}^{-1}}\nabla J(f)$ as L²-gradient of J(f). To compute ${{M}^{-1}}\nabla J(f)$ efficiently, matrix elements of M are approximated by

$$\begin{eqnarray*}&&{{\int }_{{{\mathbb{R}}^{3}}}}{{\phi }_{{\bf i}}}{{\phi }_{{\bf j}}}{\rm d}{\bf x}\approx \mathop{\prod }\limits_{k=1}^{3}\int _{\mathbb{R}}^{{}}N\left( \frac{{{x}_{k}}}{h}-{{i}_{k}} \right)N\left( \frac{{{x}_{k}}}{h}-{{j}_{k}} \right){\rm d}{{x}_{k}},\ \ {\bf i}=({{i}_{1}},{{i}_{2}},{{i}_{3}}),\ {\bf j}=({{j}_{1}},{{j}_{2}},{{j}_{3}}).\end{eqnarray*}$$

Using this approximation, ${{M}^{-1}}$ can be approximately computed quickly.

We present the main steps of the bi-gradient method. The method depends on a parameter $\beta \in [0,1]$.

Algorithm 4.1. Bi-gradient method.

• (i)  
  Compute a threshold ${{\varepsilon }_{{\rm stop}}}\gt 0$ of J(f) for stopping the iteration.
• (ii)  
  Given an initial value ${{f}^{(0)}}=0$. Set k = 0.
• (iii)  
  Compute ${{{\bf r}}_{k}}:=-\nabla J({{f}^{(k)}})$ and then ${{{\bf h}}_{k}}:={{M}^{-1}}{{{\bf r}}_{k}}$.
• (iv)  
  Compute ${{\alpha }_{k}}$ and ${{\beta }_{k}}$ such that

  $$\begin{eqnarray}&&J({{f}^{(k)}}+{{\alpha }_{k}}{{r}^{(k)}}+{{\beta }_{k}}{{h}^{(k)}})={\rm min} ,\end{eqnarray} \tag{ 12 }$$

  where ${{r}^{(k)}}$ and ${{h}^{(k)}}$ are spline functions with ${{{\bf r}}_{k}}$ and ${{{\bf h}}_{k}}$ as their coefficient vectors, respectively. The real numbers ${{\alpha }_{k}}$ and ${{\beta }_{k}}$ are obtained by solving a 2 × 2 linear system derived from (12).
• (v)  
  Set ${{d}_{k}}=({{\alpha }_{k}}/{{\beta }_{k}}){{r}^{(k)}}+{{h}^{(k)}}$. Then for the given $\beta =\frac{\alpha }{1+\alpha }\in [0,1]$, determine a ${{\tau }_{k}}$, such that

  $$\begin{eqnarray*}&&J\left( {{f}^{(k)}}+{{\tau }_{k}}((1-\beta ){{r}^{(k)}}+\beta {{d}_{k}}) \right)={\rm min} ,\end{eqnarray*}$$
• (vi)  
  Compute

  $$\begin{eqnarray*}&&{{f}^{(k+1)}}={{f}^{(k)}}+{{\tau }_{k}}((1-\beta ){{r}^{(k)}}+\beta {{d}_{k}}).\end{eqnarray*}$$
• (vii)  
  If the following stopping condition

  $$\begin{eqnarray}&&J({{f}^{(k+1)}})\leqslant {{\varepsilon }_{{\rm stop}}}\ \ {\rm or}\ \ k=M-1\end{eqnarray} \tag{ 13 }$$

  is satisfied, stop the iteration. Then ${{f}_{1}}:={{f}^{(k+1)}}$ is the required result. Otherwise, set k as $k+1$ and then go back to step 3. The integer $M\gt 1$ in (13) is a given bound for the iteration.

Computing stopping threshold. Let f^* be the exact function to be reconstructed, which is unknown. Then the measured image ${{g}_{{{{\bf d}}_{l}}}}$ can be represented as

$$\begin{eqnarray*}&&{{g}_{{{{\bf d}}_{l}}}}({\bf u})=({{X}_{{{{\bf d}}_{l}}}}{{f}^{*}})({\bf u})+{{n}_{{{{\bf d}}_{l}}}}({\bf u}),\end{eqnarray*}$$

where ${{n}_{{{{\bf d}}_{l}}}}({\bf u})$ are additive noise images. Therefore

$$\begin{eqnarray}&&J({{f}^{*}})=\frac{1}{p}\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{\mathbb{R}}^{2}}}}{{\left[ ({{X}_{{{{\bf d}}_{l}}}}{{f}^{*}})({\bf u})-{{g}_{{{{\bf d}}_{l}}}}({\bf u}) \right]}^{2}}{\rm d}{\bf u}=\frac{1}{p}\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{\mathbb{R}}^{2}}}}{{\left[ {{n}_{{{{\bf d}}_{l}}}}({\bf u}) \right]}^{2}}{\rm d}{\bf u}.\end{eqnarray} \tag{ 14 }$$

Hence we can take

$$\begin{eqnarray*}&&{{\varepsilon }_{{\rm stop}}}=\frac{1}{p}\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{\mathbb{R}}^{2}}}}{{\left[ {{n}_{{{{\bf d}}_{l}}}}({\bf u}) \right]}^{2}}{\rm d}{\bf u}.\end{eqnarray*}$$

To compute ${{\varepsilon }_{{\rm stop}}}$, we need to estimate ${{n}_{{{{\bf d}}_{l}}}}({\bf u})$. Though ${{n}_{{{{\bf d}}_{l}}}}({\bf u})$ are unknown in general, some parts of ${{n}_{{{{\bf d}}_{l}}}}({\bf u})$ are presented in the image ${{g}_{{{{\bf d}}_{l}}}}({\bf u})$. We can use these known parts to estimate ${{\varepsilon }_{{\rm stop}}}$. In this paper, we regard the part ${{[-R,R]}^{2}}\setminus \{{\bf u}:\parallel {\bf u}\parallel \leqslant R\}$ in ${{g}_{{{{\bf d}}_{l}}}}({\bf u})$ as noise. We therefore compute ${{\varepsilon }_{{\rm stop}}}$ as follows

$$\begin{eqnarray*}\begin{array}{rcl} {{\varepsilon }_{{\rm stop}}} & \approx & \frac{4}{p(4-\pi )}\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{[-R,R]}^{2}}\setminus \{{\bf u}:\parallel {\bf u}\parallel \leqslant R\}}}{{\left[ {{g}_{{{{\bf d}}_{l}}}}({\bf u}) \right]}^{2}}{\rm d}{\bf u} \\ {} & \approx & \frac{4}{p(4-\pi )}\mathop{\sum }\limits_{l=1}^{p}\mathop{\sum }\limits_{{{i}^{2}}+{{j}^{2}}\geqslant \frac{n}{2}}{{\left[ {{g}_{{{{\bf d}}_{l}}}}(i,j) \right]}^{2}}, \\ \end{array}\end{eqnarray*}$$

where ${{[-R,R]}^{2}}$ is the domain of the measured images. The integrations above are approximated by summations of the image values over the integer grid points.

The second stage for reconstructing f is as the following steps.

Algorithm 4.2. Combining geometric flow.

• (i)  
  Given an initial value ${{f}^{(0)}}={{f}_{1}}$, where f₁ is the output of algorithm 4.1 Set k = 0.
• (ii)  
  Compute $\nabla J({{f}^{(k)}})$ and $\nabla G({{f}^{(k)}})$, where J is defined by (2) and G is defined by the geometric flow in section 3.
• (iii)  
  Compute combined direction ${{{\bf h}}_{k}}=-{{\alpha }_{k}}\nabla J({{f}^{(k)}})-{{\beta }_{k}}\nabla G({{f}^{(k)}})$ (${{\lambda }_{k}}=\frac{{{\beta }_{k}}}{{{\alpha }_{k}}}$, see section 4.6.1 for detail).
• (iv)  
  Compute ${{\tau }_{k}}$ (see (21) and section 4.6.1 for detail) and then compute

  $$\begin{eqnarray}&&{{f}^{(k+1)}}={{f}^{(k)}}+{{\tau }_{k}}{{h}^{(k)}},\end{eqnarray} \tag{ 15 }$$

  where ${{h}^{(k)}}$ is a linear combination of the B-spline radial basis function with ${{{\bf h}}_{k}}$ as its coefficient vector.
• (v)  
  If the following stopping condition

  $$\begin{eqnarray}&&\frac{\parallel \nabla G({{f}^{(k)}})\parallel }{\parallel \nabla G({{f}^{(0)}})\parallel }\lt \epsilon \ \ {\rm or}\ \ k=N-1\end{eqnarray} \tag{ 16 }$$

  is satisfied, stop the iteration. Then $f:={{f}^{(k+1)}}$ is the final result. Otherwise, set $k=k+1$ and then go back to step 2. In (16), is a small number: we take it as 10⁻⁵. The integer $N\gt 1$ in (16) is a given bound for the iteration.

4.6.1. Combined direction

Let

$$\begin{eqnarray*}&&{{{\bf r}}_{k}}:=-\frac{\nabla J({{f}^{(k)}})}{\parallel \nabla J({{f}^{(k)}})\parallel },\ \ {{{\bf g}}_{k}}:=-\frac{\nabla G({{f}^{(k)}})}{\parallel \nabla G({{f}^{(k)}})\parallel }.\end{eqnarray*}$$

If ${\bf r}_{k}^{T}{{{\bf g}}_{k}}=1$ (${{{\bf g}}_{k}}={{{\bf r}}_{k}}$), then the combine direction is taken as ${{{\bf r}}_{k}}$. If ${\bf g}_{k}^{T}{{{\bf r}}_{k}}=-1$ (${{{\bf g}}_{k}}=-{{{\bf r}}_{k}}$), then the spanned space by ${{{\bf g}}_{k}}$ and ${{{\bf r}}_{k}}$ is one-dimensional (1D). Then we replace ${{{\bf g}}_{k}}$ by ${{{\bf h}}_{k-1}}$ and restart the determination of the combined direction. In the following, we assume $|{\bf r}_{k}^{T}{{{\bf g}}_{k}}|\lt 1$. Let

$$\begin{eqnarray*}&&{{\tilde{{\bf r}}}_{k}}=\frac{{{{\bf g}}_{k}}-{{\alpha }_{k}}{{{\bf r}}_{k}}}{\sqrt{1-\alpha _{k}^{2}}},\ {{\tilde{{\bf g}}}_{k}}=\frac{{{{\bf r}}_{k}}-{{\alpha }_{k}}{{{\bf g}}_{k}}}{\sqrt{1-\alpha _{k}^{2}}},\end{eqnarray*}$$

where ${{\alpha }_{k}}={\bf r}_{k}^{T}{{{\bf g}}_{k}}$. Then it is easy to see that

$$\begin{eqnarray*}&&\parallel {{\tilde{{\bf r}}}_{k}}\parallel =\parallel {{\tilde{{\bf g}}}_{k}}\parallel =1,\ \tilde{{\bf g}}_{k}^{T}{{{\bf g}}_{k}}=0,\ \tilde{{\bf r}}_{k}^{T}{{{\bf r}}_{k}}=0,\ \tilde{{\bf g}}_{k}^{T}{{{\bf r}}_{k}}=\tilde{{\bf r}}_{k}^{T}{{{\bf g}}_{k}}=\sqrt{1-\alpha _{k}^{2}}\geqslant 0.\end{eqnarray*}$$

Let

$$\begin{eqnarray*}&&{{\theta }_{0}}=-{\rm sign}({{\alpha }_{k}}){{{\rm cos} }^{-1}}\left( {\bf r}_{k}^{T}{{\tilde{{\bf g}}}_{k}} \right).\end{eqnarray*}$$

Then

$$\begin{eqnarray*}&&|{{\theta }_{0}}|\lt \pi /2.\end{eqnarray*}$$

In the plane spanned by the vectors ${{{\bf g}}_{k}}$ and ${{{\bf r}}_{k}}$, we first introduce an angle θ measured from vector ${{{\bf r}}_{k}}$ to vector ${{{\bf g}}_{k}}$. Then we define the two step-size curves

$$\begin{eqnarray}&&{{\tau }_{r}}(\theta )=\frac{2{\rm cos} (\theta ){{a}_{r}}}{{{{\rm cos} }^{2}}(\theta ){{b}_{r}}+2{\rm cos} (\theta ){\rm sin} (\theta ){{c}_{r}}+{{{\rm sin} }^{2}}(\theta ){{d}_{r}}},\ \ \ \theta \in (-\pi /2,\pi /2),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{{\tau }_{g}}(\theta )=\frac{{\rm sin} (\vartheta ){{a}_{g}}}{{{{\rm sin} }^{2}}(\vartheta ){{b}_{g}}+2{\rm cos} (\vartheta ){\rm sin} (\vartheta ){{c}_{g}}+{{{\rm cos} }^{2}}(\vartheta ){{d}_{g}}},\ \ \vartheta =\theta -{{\theta }_{0}}\in (0,\pi ),\end{eqnarray} \tag{ 18 }$$

for reducing the energies J(f) and G(f) in the directions

$$\begin{eqnarray*}\begin{array}{rcl} {\bf r}(\theta ) & = & {\rm cos} (\theta ){{{\bf r}}_{k}}+{\rm sin} (\theta ){{\tilde{{\bf r}}}_{k}},\ \ \ \theta \in (-\pi /2,\pi /2), \\ {\bf g}(\theta ) & = & {\rm sin} (\vartheta ){{{\bf g}}_{k}}+{\rm cos} (\vartheta ){{\tilde{{\bf g}}}_{k}},\ \ \vartheta =\theta -{{\theta }_{0}}\in (0,\pi ), \\ \end{array}\end{eqnarray*}$$

respectively. The numbers a_r, b_r, c_r and d_r are given by (24)–(27). Constants a_g, b_g, c_g and d_g are given by (29)–(31). Note that

$$\begin{eqnarray*}&&{\bf r}(\theta )={\bf g}(\theta ),\ \theta \in [{{\theta }_{0}},\pi /2].\end{eqnarray*}$$

Let

$$\begin{eqnarray*}&&\tau (\theta )={\rm min} \left\{ {{\tau }_{r}}(\theta ),{{\tau }_{g}}(\theta ) \right\},\ \theta \in [{{\theta }_{0}},\pi /2].\end{eqnarray*}$$

Then the curve $\tau (\theta )$ is the right step-size for the direction ${\bf r}(\theta )$, which makes J(f) not increase and G(f) decrease. We want to find a ${{\theta }^{*}}\in [{{\theta }_{0}},\pi /2]$, that leads to the maximal reduction of the energy G(f). Since

$$\begin{eqnarray*}\begin{array}{rcl} G(f+\tau (\theta ){\bf g}(\theta )) & = & G(f)+\tau (\theta )\nabla G{{(f)}^{T}}{\bf g}(\theta )+O({{\tau }^{2}}(\theta )) \\ {} & = & G(f)+{\rm sin} (\theta -{{\theta }_{0}})\tau (\theta )\nabla G{{(f)}^{T}}{{{\bf g}}_{k}}+O({{\tau }^{2}}(\theta )) \\ {} & = & G(f)-{\rm sin} (\theta -{{\theta }_{0}})\tau (\theta )\parallel \nabla G(f)\parallel +O({{\tau }^{2}}(\theta )). \\ \end{array}\end{eqnarray*}$$

Omitting the higher order term $O({{\tau }^{2}}(\theta ))$, we determine our best ${{\theta }^{*}}$ as follows

$$\begin{eqnarray*}&&{{\theta }^{*}}={\rm arg} \mathop{{\rm max} }\limits_{\theta \in [{{\theta }_{0}},\pi /2]}{\rm sin} (\theta -{{\theta }_{0}})\tau (\theta ).\end{eqnarray*}$$

Then the combined direction is given as

$$\begin{eqnarray}&&{{{\bf h}}_{k}}={\rm sin} ({{\theta }^{*}}-{{\theta }_{0}}){{{\bf g}}_{k}}+{\rm cos} ({{\theta }^{*}}-{{\theta }_{0}}){{\tilde{{\bf g}}}_{k}}={\rm cos} ({{\theta }^{*}}){{{\bf r}}_{k}}+{\rm sin} ({{\theta }^{*}}){{\tilde{{\bf r}}}_{k}}.\end{eqnarray} \tag{ 19 }$$

Remark 4.2. Since $\tilde{{\bf r}}_{k}^{T}{{{\bf r}}_{k}}=0$, we know from (19) that

$$\begin{eqnarray}&&{\bf h}_{k}^{T}{{{\bf r}}_{k}}={\rm cos} ({{\theta }^{*}})\gt 0.\end{eqnarray} \tag{ 20 }$$

Hence ${{{\bf h}}_{k}}$ is a descent direction of J(f).

Having ${{\theta }^{*}}$ and the combined direction ${{{\bf h}}_{k}}$, we compute the step-size in the direction ${{{\bf h}}_{k}}$ as follows

$$\begin{eqnarray}&&{{\tau }_{k}}=\tau ({{\theta }^{*}}).\end{eqnarray} \tag{ 21 }$$

4.6.2. Step-size curves

We first consider the step-size curve ${{\tau }_{r}}(\theta )$. Given a descent direction h, from minimizing $J({{f}^{(k)}}+\tau h)$, we have

$$\begin{eqnarray*}&&J^{\prime} ({{f}^{(k)}}+\tau h)=0,\end{eqnarray*}$$

we can determine τ as

$$\begin{eqnarray*}&&\tau =-\frac{N1}{{{D}_{1}}},\end{eqnarray*}$$

with

$$\begin{eqnarray}&&{{N}_{1}}=\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{\mathbb{R}}^{2}}}}\left[ {{X}_{{{{\bf d}}_{l}}}}{{f}^{(k)}}-{{g}_{{{{\bf d}}_{l}}}}({\bf u}) \right]{{X}_{{{{\bf d}}_{l}}}}h\ {\rm d}{\bf u},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{D}_{1}}=\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{\mathbb{R}}^{2}}}}{{\left[ {{X}_{{{{\bf d}}_{l}}}}h \right]}^{2}}\ {\rm d}{\bf u},\end{eqnarray} \tag{ 23 }$$

To make J(f) not increase, we redefine

$$\begin{eqnarray*}&&\tau =-\frac{2{{N}_{1}}}{{{D}_{1}}}.\end{eqnarray*}$$

Taking

$$\begin{eqnarray*}&&{\bf h}={\rm cos} (\theta ){{{\bf r}}_{k}}+{\rm sin} (\theta ){{\tilde{{\bf r}}}_{k}},\end{eqnarray*}$$

then we obtain (17) with

$$\begin{eqnarray}&&{{a}_{r}}=-\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{\mathbb{R}}^{2}}}}\left[ {{X}_{{{{\bf d}}_{l}}}}{{f}^{(k)}}-{{g}_{{{{\bf d}}_{l}}}}({\bf u}) \right]{{X}_{{{{\bf d}}_{l}}}}{{r}_{k}}\ {\rm d}{\bf u},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{{b}_{r}}=\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{\mathbb{R}}^{2}}}}{{\left[ {{X}_{{{{\bf d}}_{l}}}}{{r}_{k}} \right]}^{2}}\ {\rm d}{\bf u},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{{c}_{r}}=\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{\mathbb{R}}^{2}}}}{{X}_{{{{\bf d}}_{l}}}}{{r}_{k}}{{X}_{{{{\bf d}}_{l}}}}{{\tilde{r}}_{k}}\ {\rm d}{\bf u},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{{d}_{r}}=\mathop{\sum }\limits_{l=1}^{p}{{\int }_{{{\mathbb{R}}^{2}}}}{{\left[ {{X}_{{{{\bf d}}_{l}}}}{{\tilde{r}}_{k}} \right]}^{2}}\ {\rm d}{\bf u}.\end{eqnarray} \tag{ 27 }$$

Now we consider step-size curve ${{\tau }_{g}}(\theta )$. Since

$$\begin{eqnarray*}\begin{array}{rcl} G(f+\tau (g(\theta ))) & = & G(f)+\tau \nabla G{{(f)}^{T}}{\bf g}(\theta ) \\ {} & {} & +\frac{1}{2}{{\tau }^{2}}{\bf g}{{(\theta )}^{T}}{{\nabla }^{2}}G(f){\bf g}(\theta )+O({{\tau }^{3}}), \\ \end{array}\end{eqnarray*}$$

then from $G^{\prime} (f+\tau (g(\theta )))=0$ and omitting the higher order term $O({{\tau }^{2}})$, we obtain

$$\begin{eqnarray}&&\tau =-\frac{\nabla G{{(f)}^{T}}{\bf g}(\theta )}{{\bf g}{{(\theta )}^{T}}{{\nabla }^{2}}G(f){\bf g}(\theta )}.\end{eqnarray} \tag{ 28 }$$

Substituting

$$\begin{eqnarray*}&&{\bf g}(\theta )={\rm sin} (\vartheta ){{{\bf g}}_{k}}+{\rm cos} (\vartheta ){{\tilde{{\bf g}}}_{k}},\end{eqnarray*}$$

into (28), we obtain (18), with

$$\begin{eqnarray}&&{{a}_{g}}=-{{[\nabla G({{f}^{(k)}})]}^{T}}{{{\bf g}}_{k}}=-{{G}_{{{g}_{k}}}}({{f}^{(k)}}),\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{{b}_{g}}={\bf g}_{k}^{T}{{\nabla }^{2}}G({{f}^{(k)}}){{{\bf g}}_{k}}={{G}_{{{g}_{k}}{{g}_{k}}}}({{f}^{(k)}}),\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{{c}_{g}}={\bf g}_{k}^{T}{{\nabla }^{2}}G({{f}^{(k)}}){{\tilde{{\bf g}}}_{k}}={{G}_{{{g}_{k}}\tilde{{{g}_{k}}}}}({{f}^{(k)}}),\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{{d}_{g}}=\tilde{{\bf g}}_{k}^{T}{{\nabla }^{2}}G({{f}^{(k)}}){{\tilde{{\bf g}}}_{k}}={{G}_{\tilde{{{g}_{k}}}\tilde{{{g}_{k}}}}}({{f}^{(k)}}),\end{eqnarray} \tag{ 32 }$$

where ${{G}_{{{g}_{k}}}}({{f}^{(k)}})$ and ${{G}_{{{g}_{k}}{{g}_{k}}}}({{f}^{(k)}})$ are the first- and second-order variations of G with respect to g_k. Appendix B gives the computational details.

4.6.3. Discussions

Now we show that the f produced by algorithm 4.2 satisfies the condition (4). To illustrate this, let us consider a general form functional

$$\begin{eqnarray*}&&\mathcal{E}(f)={{\int }_{{{\mathbb{R}}^{3}}}}E(f({\bf x}))\ {\rm d}{\bf x}.\end{eqnarray*}$$

Then the partial derivative with the coefficient of the spline function f is

$$\begin{eqnarray*}&&\frac{\mathcal{E}(f)}{\partial {{f}_{{\bf i}}}}={{r}_{{\bf i}}}\ \ {\rm with}\ \ {{r}_{{\bf i}}}=-{{\int }_{{{\mathbb{R}}^{3}}}}E^{\prime} (f){{\phi }_{{\bf i}}}\ {\rm d}{\bf x}.\end{eqnarray*}$$

Let

$$\begin{eqnarray*}&&r=\mathop{\sum }\limits_{{\bf i}}{{r}_{{\bf i}}}{{\phi }_{{\bf i}}}.\end{eqnarray*}$$

Then the changing rate of $\mathcal{E}(f)$ in the direction r is

$$\begin{eqnarray*}\begin{array}{rcl} {{\left. \frac{{\rm d}}{{\rm d}t}\mathcal{E}(f+tr) \right|}_{t=0}} & = & {{\left. \frac{{\rm d}}{{\rm d}t}{{\int }_{{{\mathbb{R}}^{3}}}}E(f+tr)\ {\rm d}{\bf x} \right|}_{t=0}} \\ {} & = & {{\int }_{{{\mathbb{R}}^{3}}}}E^{\prime} (f)r\ {\rm d}{\bf x} \\ {} & = & \mathop{\sum }\limits_{{\bf i}}{{r}_{{\bf i}}}{{\int }_{{{\mathbb{R}}^{3}}}}E^{\prime} (f){{\phi }_{{\bf i}}}\ {\rm d}{\bf x} \\ {} & = & -\mathop{\sum }\limits_{{\bf i}}r_{{\bf i}}^{2}. \\ \end{array}\end{eqnarray*}$$

Hence, $\mathcal{E}(f)$ is non-increasing in the direction r. Let $h={{\sum }_{{\bf i}}}{{h}_{{\bf i}}}{{\phi }_{{\bf i}}}$ be another function, if

$$\begin{eqnarray}&&\mathop{\sum }\limits_{{\bf i}}{{h}_{{\bf i}}}{{r}_{{\bf i}}}\geqslant 0,\end{eqnarray} \tag{ 33 }$$

then it is easy to derive that

$$\begin{eqnarray*}&&{{\left. \frac{{\rm d}}{{\rm d}t}\mathcal{E}(f+th) \right|}_{t=0}}=-\mathop{\sum }\limits_{{\bf i}}{{h}_{{\bf i}}}{{r}_{{\bf i}}}\leqslant 0.\end{eqnarray*}$$

Hence, $\mathcal{E}(f)$ is non-increasing in the direction h.

From the discussion above, we know that J(f) is not increasing in the direction ${{{\bf h}}_{k}}$ at ${{f}^{(k)}}$, as (20) is satisfied. Furthermore, the step-size ${{\tau }_{k}}$ given by (21) makes $J({{f}^{(k+1)}})\leqslant J({{f}^{(k)}})$. Therefore, the function sequence $\{{{f}^{(k)}}\}$ produced by algorithm 4.2 causes the sequence $\{J({{f}^{(k)}})\}$ to decrease monotonically.

Given a volume data $F=\{{{f}_{ijk}}\}$ with size $131\times 131\times 131$, we first generate 5000 projection images $\{{{I}_{i}}\}_{i=1}^{5000}$ with size 131 × 131 at randomly chosen projection directions $\{{{{\bf d}}_{i}}\}_{i=1}^{5000}$ on the unit sphere (figure 2 shows the first five projection images). Then we reconstruct F using different reconstruction algorithms. Let ${{R}^{(l)}}=\{r_{ijk}^{(l)}\}_{ijk=1}^{131}$, $l=1,\ldots ,5$, be the reconstructed volumes using WBP, and our bi-gradient method with $\beta =0$, $\beta =\frac{1}{3}$, $\beta =\frac{1}{2}$ and $\beta =1$, respectively. Table 1 lists the L²-errors ${{E}^{(l)}}$, where ${{E}^{(l)}}$ is defined as

$$\begin{eqnarray*}&&{{E}^{(l)}}=\sqrt{\frac{1}{{{m}^{3}}}\mathop{\sum }\limits_{i=1}^{m}\mathop{\sum }\limits_{j=1}^{m}\mathop{\sum }\limits_{k=1}^{m}{{\left[ {{f}_{ijk}}-r_{ijk}^{(l)} \right]}^{2}}},\ l=1,\ldots ,5,\ m=131.\end{eqnarray*}$$

Since the data is clean, we do not use the regularization term in the reconstruction equation. The L²-errors are monotonically decreasing as the iteration number increases. It can be clearly observed that as $\beta$ increases from 0 to 1, the errors decrease. Choosing $\beta =1$ gives the most accurate result. Following $\beta =1$, the next best result is given by WBP with L²-error 0.08 131. Following WBP are $\beta =\frac{1}{2}$, $\beta =\frac{1}{3}$, and $\beta =0$. More iterations will make the case $\beta =\frac{1}{2}$ better than WBP. For instance, the L²-error is 0.07 855 after 40 iterations. In the next subsection, we will illustrate that the most accurate method may not be the best method for data with high noise.

Zoom In Zoom Out Reset image size

Table 1.  L²-errors ${{E}^{(2)}},\ldots ,{{E}^{(5)}}$ for different iteration numbers for the clean data. ${{E}^{(1)}}=0.08\;131$.

Iteration number	$\beta =0$	$\beta =\frac{1}{3}$	$\beta =\frac{1}{2}$	$\beta =1$
5	0.15 519	0.15 368	0.14 793	0.07 803
10	0.13 726	0.13 419	0.12 313	0.05 790
15	0.12 838	0.12 416	0.10 977	0.05 686
20	0.12 167	0.11 635	0.09 931	0.05 674
25	0.11 733	0.11 115	0.09 246	0.05 672
30	0.11 354	0.10 652	0.08 653	0.05 671

Figure 3 shows the central slices of the reconstructed volume data, where figure (a) is from the initial data. Figures (b–f) are produced from the reconstruction results of the WBP, bi-gradient method with $\beta =0$, $\beta =\frac{1}{3}$, $\beta =\frac{1}{2}$, and $\beta =1$, respectively. These figures show the same conclusions as the numerical results.

Zoom In Zoom Out Reset image size

Remark 5.3. The increase in the accuracy of the bi-gradient method as β increases from 0 to 1 is due to the fact that when $\beta =0$, we search the minimal point in the 1D space ${\rm span}\;[{{{\bf r}}_{k}}]$, but when $\beta =1$, we search the minimal point in the 2D space ${\rm span}\;[{{{\bf r}}_{k}},{{{\bf h}}_{k}}]$. When $\beta \in (0,1)$, the method is an average of the cases $\beta =0$ and $\beta =1$.

At this point, we generate a new set of data by adding additive white Gaussian noise to each of the images I_i, produced in the previous subsection with SNR = 1.0. Figure 4 shows five noisy images of the ones shown in figure 2. We then reconstruct f using different reconstruction algorithms from the noisy images. Because the data is noisy, a regularization term is used. As before, we use ${{R}^{(r)}}$ to represent the reconstructed volume. For our bi-gradient method, we iterate 30 steps. In table 2, we list the energies, defined by (2), for the reconstruction results after 30 iterations. The energies show that the bi-gradient method with $\beta =1$ yields minimal energy. Hence, it is indeed the most accurate method. However, the most accurate method may not lead to the most meaningful results.

Zoom In Zoom Out Reset image size

Table 2.  Energies of different methods after 30 iterations.

SNR	$\beta =0$	$\beta =\frac{1}{3}$	$\beta =\frac{1}{2}$	$\beta =1$
1.0	101 747 81.537 409	101 552 48.877 989	100 916 05.869 403	990 8452.694 144

Figure 5 shows the central slices of the reconstructed volume data. The figures, from the first to the fifth, are the central slices of the reconstruction volumes of WBP, our bi-gradient method with $\beta =0$, $\beta =\frac{1}{3}$, $\beta =\frac{1}{2}$, and $\beta =1$, respectively. It can be clearly observed that the bi-gradient method with $\beta =1$ gives the noisiest result with detailed structures. After that, the bi-gradient method with $\beta =\frac{1}{2}$, $\beta =\frac{1}{3}$ and $\beta =0$ follow successively.

Zoom In Zoom Out Reset image size

In table 3, we list the L²-errors between the reconstructed volumes and the exact initial volume data for iterations 5, 10, 15, 20, 25, and 30. The errors are not monotonically decreasing as the iteration number increases. It is easy to see that the most accurate method ($\beta =1$) for the clean data leads to the maximal L² error, while the other three cases lead to similar L² error bounds. All these three cases are better than WBP for iterations 10, 15, 20, 25, and 30. The iso-surfaces of the reconstructed volume data, as shown in figure 6, demonstrate that:

Zoom In Zoom Out Reset image size

Table 3.  L²-errors ${{E}^{(2)}},\ldots ,{{E}^{(5)}}$ for different iteration numbers. The data are noised with SNR = 1.0. ${{E}^{(1)}}=0.15\;005$.

Iteration number	${{E}^{(2)}}$	${{E}^{(3)}}$	${{E}^{(4)}}$	${{E}^{(5)}}$
5	0.15 521	0.15 379	0.15 071	0.36 863
10	0.13 770	0.13 567	0.14 074	0.45 722
15	0.13 374	0.12 946	0.13 406	0.45 476
20	0.13 394	0.12 671	0.12 680	0.45 221
25	0.13 419	0.12 618	0.12378	0.45 170
30	0.13 428	0.12 571	0.12 203	0.45 145

• (i)  
  The bi-gradient method with $\beta =0$ gives the most coarse level structure, then followed by $\beta =\frac{1}{3}$, $\beta =\frac{1}{2}$ and $\beta =1.0$.
• (ii)  
  The bi-gradient method with $\beta =1$ does not yield valuable iso-surfaces.

Remark 5.4. The noise is added to the 2D images I_i, not to the volume data F. If the 2D image is highly polluted, such that the noise level is much higher than the signal, the reconstructed volume data from these polluted images differs greatly from the initial volume data F. Hence, it is not reasonable to require the reconstruction volumes to be close to the initial volume. Therefore, the L²-error between the initial data and the reconstructed data may not be a reasonable measure to evaluate the reconstruction method. In the next subsection, we explain how we used Fourier shell correlation (FSC) to evaluate the performance of the reconstruction methods.

Remark 5.5. In the previous section, we showed that, for the clean data set, as β increases from 0 to 1, the accuracy of the bi-gradient method also increases. In this subsection, we observe that $\beta =0$ leads to the smoothest result. Now we explain the reason. From our previous work [5], we know that if we choose ${{{\bf f}}^{(0)}}=0$ and search the minimal point in the direction ${{{\bf r}}_{k}}$, ${{{\bf f}}^{(k)}}$ converges to the minimal solution ${{{\bf f}}^{*}}$ of the equation $R{\bf f}=G$ in the space $N{{(R)}^{\bot }}=\{{\bf x}:{{{\bf x}}^{T}}{\bf y}=0,\forall {\bf y}\in N(R)\}$ in the sense of the Euclidean norm $\parallel {\bf f}\parallel =\sqrt{{{{\bf f}}^{T}}{\bf f}}$, where N(R) is the null space of the matrix R. While searching the minimal point in the direction ${{{\bf h}}_{k}}$, ${{{\bf f}}^{(k)}}$ converges to the minimal solution ${{{\bf f}}^{*}}$ of the equation $R{\bf f}=G$ in the space $N{{(R)}^{{{\bot }^{M}}}}=\{{\bf x}:{{{\bf x}}^{T}}M{\bf y}=0,\forall {\bf y}\in N(R)\}$ in the sense of the norm $\parallel {\bf f}{{\parallel }_{M}}=\sqrt{{{{\bf f}}^{T}}M{\bf f}}$. The minimal property of the Euclidean norm makes the case that $\beta =0$ yields the most smooth result. Table 4 lists some Euclidean norms for the cases $\beta =0,\frac{1}{3},\frac{1}{2},1$. It can be seen that as β increases, the Euclidean norm also increases.

Given a volume data $F=\{{{f}_{ijk}}\}$, of size $151\times 151\times 151$, we first generate 200 00 projection images $\{{{I}_{i}}\}_{i=1}^{200\,00}$ with size 151 × 151 at randomly chosen projection directions $\{{{{\bf d}}_{i}}\}_{i=1}^{200\,00}$ on the unit sphere. These images are polluted by adding additive white Gaussian noise with SNR=0.1 (figure 7 shows five of these 'noised' images). These noised images are split randomly into datasets A and B, with each dataset containing 100 00 images. Then, we reconstruct F using our bi-gradient method with $\beta =0,\frac{1}{3},\frac{1}{2},\frac{2}{3}$ for each of the datasets (in this section, we do not take $\beta =1$, as the examples in the previous subsection have shown that taking $\beta =1$ does not lead to desirable results). Because the data are extremely noisy, a regularization term is used. The aim of splitting the dataset into two parts is to allow us to examine the correlation of the reconstructed results using FSC.

Zoom In Zoom Out Reset image size

Table 4.  Euclidean norms of the reconstruction results. For WBP, it is 4.067 898e+02.

Iteration number	$\beta =0$	$\beta =\frac{1}{3}$	$\beta =\frac{1}{2}$	$\beta =1$
5	3.046 299e+02	3.048 806e+02	3.061 394e+02	3.759 363e+02
10	3.225 103e+02	3.231 880e+02	3.270 741e+02	4.281 541e+02
15	3.321 209e+02	3.335 966e+02	3.378 084e+02	4.239 439e+02
20	3.340 637e+02	3.356 739e+02	3.410 555e+02	4.227 851e+02
25	3.355 638e+02	3.377 481e+02	3.444 411e+02	4.227 860e+02
30	3.360 086e+02	3.384 513e+02	3.451 244e+02	4.227 148e+02

As before, we use ${{R}^{(r)}}$ to represent the reconstructed volume. For each value of $\beta$ in our method, the algorithm iterates 30 steps. In table 5, we list the energies for the reconstruction results after 30 iterations. These energies are decreasing as the value of $\beta$ increases and $\beta =\frac{2}{3}$ yields minimal energies for the two sets of data.

Table 5.  Energies of different methods after 30 iterations.

Data	$\beta =0$	$\beta =\frac{1}{3}$	$\beta =\frac{1}{2}$	$\beta =\frac{2}{3}$
Data set A	135 429 907.779	135 138 283.715	134 896 618.185	134 531 987.307
Data set B	135 409 723.391	135 117 692.553	134 875 769.350	134 510 933.838

Figure 8 shows the central slices of the volume data. Figures on the first to the fifth columns are from the reconstruction results of WBP, our bi-gradient method with $\beta =0,\frac{1}{3},\frac{1}{2},\frac{2}{3}$. It can be clearly observed that the case $\beta =\frac{2}{3}$ gives the noisiest result with detailed structures. After that, the cases $\beta =\frac{1}{2},\frac{1}{3},0$ follow successively.

Zoom In Zoom Out Reset image size

In tables 6 and 7, we list the L²-errors between the reconstructed volumes and the exact initial volume data for iterations 5, 10, 15, 20, 25, and 30. It is easy to see that the bi-gradient method with $\beta =0$ leads to the minimal L² error, followed by the cases $\beta =\frac{1}{3}$, $\beta =\frac{1}{2}$, and $\beta =\frac{2}{3}$, and all of these are better than WBP after 30 iterations. The iso-surfaces of the reconstructed volume data, as shown in figure 9, demonstrate that the bi-gradient method with $\beta =0$ gives the coarsest level structure, followed by $\beta =\frac{1}{3}$, $\beta =\frac{1}{2}$, and $\beta =\frac{2}{3}$.

Zoom In Zoom Out Reset image size

Table 6.  L²-errors ${{E}^{(2)}},\ldots ,{{E}^{(5)}}$ for different iteration numbers and for data set A. ${{E}^{(1)}}=0.29\;193$.

Iteration number	${{E}^{(2)}}$	${{E}^{(3)}}$	${{E}^{(4)}}$	${{E}^{(5)}}$
5	0.17 066	0.16 985	0.17 464	0.19 615
10	0.14 953	0.16 200	0.19 536	0.28 074
15	0.14 286	0.15 583	0.19 542	0.29 516
20	0.14 033	0.14 994	0.18 813	0.28 671
25	0.13 990	0.14 878	0.18 605	0.28 404
30	0.13 981	0.15 085	0.18 804	0.28 388

Table 7.  L²-errors ${{E}^{(2)}},\ldots ,{{E}^{(5)}}$ for different iteration numbers and for data set B. ${{E}^{(1)}}=0.29\;208$.

Iteration number	${{E}^{(2)}}$	${{E}^{(3)}}$	${{E}^{(4)}}$	${{E}^{(5)}}$
5	0.17 064,	0.16 984	0.17 466	0.19 621
10	0.14 956	0.16 210	0.19 553	0.28 105
15	0.14 290	0.15 600	0.19 576	0.29 572
20	0.14 035	0.15 018	0.18 841	0.28 768
25	0.13 992	0.14 895	0.18 635	0.28 452
30	0.13 984	0.15 097	0.18 836	0.28 435

Table 8 lists the resolutions for each of the reconstructed volume pairs after 10, 20 and 30 iterations. The resolution is computed at the place where the value of FSC is 0.5. It can be seen that for most of the cases, the resolution of each reconstruction result is higher than that produced by WBP. Figure 10 shows the FSC curves for the five cases and different numbers of iterations.

Zoom In Zoom Out Reset image size

Table 8.  Resolution of different methods. The resolution of WBP is 12.362 563 Å.

Iteration number	$\beta =0$	$\beta =\frac{1}{3}$	$\beta =\frac{1}{2}$	$\beta =\frac{2}{3}$
10	12.055 628 Å	12.107 415 Å	12.152 935 Å	12.227 004 Å
20	12.079 924 Å	12.199 175 Å	12.280 208 Å	12.391 474 Å
30	12.108 918 Å	12.334 633 Å	12.476 160 Å	12.527 178 Å

Remark 5.6. Both the numerical and illustrative results in this subsection reveal that the reconstruction results from datasets A and B are very similar, giving further evidence that the bi-gradient method is robust.

The reconstruction in the previous subsection is conducted for each β. The computation conducted in this way allows examination of the performance of our method. In a real construction, we do not need to compute the reconstruction result for each β. We have noticed that as β increases from 0 to 1, structures from coarse to fine are captured. We now combine these computations into one loop, aiming to obtain a sequence of multi-scale reconstruction results. For a given $K\gt 1$, we assume we are going to reconstruct a sequence of volume data ${{F}_{0}},{{F}_{1}},\ldots ,{{F}_{K}}$ from coarse to fine. Then we carry out the following:

• (i)  
  Construct an initial volume data F₀ using the bi-gradient method with $\beta =0$. The iteration numbers M and N in the the bi-gradient method are taken as 12 for both the first and second stages.
• (ii)  
  For $k=1,2,\ldots ,K$ we do the following: set $\beta =\frac{\alpha }{1+\alpha }$ with $\alpha ={{2}^{k+1-K}}$ and reconstruct volume data F_k using the bi-gradient method with ${{F}_{k-1}}$ as the initial value. The iteration numbers M and N in the bi-gradient method are taken as 5 for both the first and second stages.

For instance, if we take K = 3, and reconstruct the volume data ${{F}_{0}},{{F}_{1}},{{F}_{2}},{{F}_{3}}$ using the data in the previous section, we obtain the results corresponding to $\beta =0,\frac{1}{3},\frac{1}{2},\frac{2}{3}$. The obtained results were very close to the results obtained in the previous subsection; hence, they are not presented again.