\documentclass{article}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amsbsy}
\title{A 3D finite element based diffuse optical image reconstruction algorithm \-- theory and implementation}

\author{Qianqian Fang}

\date{\today}


\begin{document}
\maketitle              % you need to define \title{..}

Redbird-m is a MATLAB toolbox designed for solving the forward and inverse
problems for diffuse optical tomography (DOT). Redbird-m was ported from 
Redbird - a FORTRAN90 based software written by the same author.
In the forward modeling, a 3D finite element method (FEM)
facilitated by an iterative multi-right-hand-side quasi-minimal residual (QMR) solver was used to model
the forward solution for RF and CW light diffusion. The 3D images for tissue chromorphore 
concentrations and tissue scattering coefficients were
reconstructed with an iterative Gauss-Newton method where the
Jacobian matrix was build by the nodal adjoint method. A simultaneous 
source-detector (SD) coupling coefficient estimation process was 
implemented in conjunction with the optical parameter reconstruction. 
This document summarizes the mathematical treatment for the forward 
and inverse problems used in this algorithm.

\section{Mathematical derivations for the forward model}
The diffusion equation in the time-domain can be expressed as
\begin{equation}
\label{eq:diffusionTD} -\nabla\cdot
D(r)\nabla\Phi(r,t)+\mu_a(r)\Phi(r,t)+\frac{1}{c}\frac{\partial
\Phi(r,t)}{\partial t}=S_0(r,t)
\end{equation}
where $D(r)=\frac{1}{3(\mu_a+\mu_s')}$ is the diffusion
coefficient (unit cm); $c=\frac{c_0}{n}$ is the speed of light in the
medium (unit cm/s); $S_0(r,t)$ is the source. With assumed time dependence
$\exp(j\omega t)$, the frequency-domain (FD) diffusion equation can be
written as
\begin{equation}
\label{eq:diffusionFD} -\nabla\cdot
D(r)\nabla\Phi(r)+\left(\mu_a(r)+\frac{j\omega}{c}\right)\Phi(r)=S_0(r)
\end{equation}
where $\omega$ is the angular frequency and $\Phi(r)$ is the
phasor of $\Phi(r,t)$.

Integrating both sides of (\ref{eq:diffusionFD}) with a set of weight functions $\varphi_j(r) (j=1,...,M)$
where $M$ is the total number of basis functions, over the forward space $\Omega$, subsequently, applying the following vector
identity
\begin{equation}
\label{eq:veciden}\nabla\cdot f(r)\vec{g}(r)=f(r)\nabla\cdot
\vec{g}(r)+\vec{g}(r)\cdot\nabla f(r)
\end{equation}
we get
\begin{eqnarray}
\label{eq:femstep1}
-\int_{\Omega}{(D(r)\nabla^2\Phi(r))\varphi_j(r)dr}&-&\int_{\Omega}{\nabla D(r)\cdot\nabla\Phi(r)
\varphi_j(r)dr}\\\nonumber &+&\int_{\Omega}
{(\mu_a(r)+\frac{j\omega}{c})\varphi_j(r)\Phi(r)dr}=\int_{\Omega}
{S_0(r)\varphi_j(r)dr}
\end{eqnarray}
Assume $D(r)$ is constant for each forward element, i.e. element-based properties, 
we have $\nabla D(r)=0$ and the second term in (\ref{eq:femstep1}) becomes zero. 
However, if we define the optical properties on the nodes using a reconstruction mesh,
the second term will be non-zero. In this case, we expand $D(r)$ and $\mu_{a}(r)$ by
piece-wise linear basis functions
\begin{equation}
\label{eq:parammesh}D(r)=\sum_{i}D_i*\phi_i(r) \\
\mu_{a}(r)=\sum_{i}\mu_{a_i}*\phi_i(r)
\end{equation}

as well as the Green's first identity
\begin{equation}
\label{eq:Green1} \int_{\Omega}{u(r)\nabla^2
v(r)dr}=-\int_{\Omega}{\nabla u(r)\cdot\nabla
v(r)dr}+\int_{\partial\Omega} {u(r)\nabla v(r)dr}
\end{equation}
we can rewrite (\ref{eq:diffusionFD}) as
\begin{eqnarray}
\label{eq:fem}
\int_{\Omega}{D(r)\nabla\varphi_j(r)\cdot\nabla\Phi(r)dr}&-&\int_{\partial\Omega}{D(r)\varphi_j(r)\nabla\Phi(r)\cdot
d\hat{r}}\\\nonumber &+&\int_{\Omega}
{(\mu_a(r)+\frac{j\omega}{c})\varphi_j(r)\Phi(r)dr}=\int_{\Omega}
{S_0(r)\varphi_j(r)dr}
\end{eqnarray}
With Galerkin's method, i.e. the basis function is identical as
the weight function, $\Phi(r)$ is expanded as
$\Phi(r)=\sum_{i=1}^{4}\Phi_i\varphi_i(r)$ over each linear
Largarange forward element, and Equ. (\ref{eq:fem}) becomes
\begin{eqnarray}
\label{eq:galerkin}\nonumber&&\sum_{i}\Phi_i\left(\left\langle
D(r)\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle -
\left\langle
D(r)\nabla\varphi_i(r)\varphi_j(r)\right\rangle_{\partial\Omega}\right.
\\&&+ \left.\left\langle
\left(\mu_a(r)+
\frac{j\omega}{c}\right)\varphi_i(r)\varphi_j(r)\right\rangle\right)=\left\langle
S_0(r)\varphi_j(r)\right\rangle
\end{eqnarray}
where $\langle u(r)\rangle$ denotes $\int_{\Omega}u(r)dr$ and
$\langle u(r)\rangle_{\partial\Omega}$ denotes
$\int_{\partial\Omega}u(r) dr$.

Expanding $\mu_a(r)$ and $D(r)$ on the parameter mesh basis as
$\mu_a(r)=\sum_{k=1}^{4}\mu_a^k\phi_k(r)$ and
$D(r)=\sum_{k=1}^{4}D_k\phi_k(r)$, we get the Galerkin weak form
equation as
\begin{eqnarray}
\label{eq:weakform}\nonumber&&\sum_{i}\Phi_i\left(\sum_{k}D_k\left\langle
\phi_k(r)\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle -
\sum_{k}D_k\left\langle
\phi_k(r)\nabla\varphi_i(r)\varphi_j(r)\right\rangle_{\partial\Omega}\right.
\\&&\left.+\sum_{k}\left(\mu_{a}^k+\frac{j\omega}{c}\right)\left\langle\phi_k(r)\varphi_i(r)\varphi_j(r)\right\rangle\right)=\left\langle
S_0(r)\varphi_j(r)\right\rangle
\end{eqnarray}
where $\phi_k$ is the basis function of the parameter mesh.

\section{Boundary condition}

If the extrapolation boundary condition is used, the boundary
integration term in (\ref{eq:weakform}) is then dropped out,
leaving
\begin{eqnarray}
\label{eq:weakformzerobd}\sum_{i}&&\Phi_i\left(\sum_{k}D_k\left\langle
\phi_k(r)\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle\right.\\\nonumber
+&&\left.\sum_{k}\left(\mu_{a}^k+\frac{j\omega}{c}\right)\left\langle\phi_k(r)\varphi_i(r)\varphi_j(r)\right\rangle\right)=\left\langle
S_0(r)\varphi_j(r)\right\rangle
\end{eqnarray}

However, it has been shown in Haskell 1994, the partial current boundary
condition is more accurate. In this boundary condition, the fluence
satisfies the following condition on the boundary $\partial\Omega$

\begin{equation}
\label{eq:partialcurrent} \left(\boldsymbol\Phi=l_s\nabla{\boldsymbol\Phi}\right)|_{\partial\Omega}
\end{equation}
where $l_s$ is defined in Eq. 2.4.1 in Haskell, as
\begin{equation}
l_s=\frac{1+R_{eff}}{1-R_{eff}}2D
\end{equation}
and $R_{eff}=\frac{R_\phi+R_j}{2-R_\phi+R_j}$ is the effective reflection coefficient, and 
\begin{eqnarray}
R_\phi &=& \int_0^{\pi/2}2\sin \theta \cos\theta R_{Fresnel(\theta)}d\theta\\
R_j &=& \int_0^{\pi/2}3\sin \theta \cos^2\theta R_{Fresnel(\theta)}d\theta \\
R_{Fresnel}(\theta) &=& \frac{1}{2}\left(\frac{n \cos\theta' - n_{out} \cos\theta}{n \cos\theta' + n_{out} \cos\theta}\right)^2 \\
   &+&\frac{1}{2}\left(\frac{n \cos\theta - n_{out} \cos\theta'}{n \cos\theta + n_{out} \cos\theta'}\right)^2
\end{eqnarray}

Replacing $\nabla{\boldsymbol\Phi}=\boldsymbol\Phi/l_s$ to the boundary term in (\ref{eq:weakform}), we
have

\begin{eqnarray}
\label{eq:weakformpartial}\nonumber&&\sum_{i}\Phi_i\left(\sum_{k}D_k\left\langle
\phi_k(r)\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle -
\sum_{k}\frac{1-R_{eff}}{2(1+R_{eff})}\left\langle
\phi_k(r)\varphi_j(r)\right\rangle_{\partial\Omega}\right.
\\&&\left.+\sum_{k}\left(\mu_{a}^k+\frac{j\omega}{c}\right)\left\langle\phi_k(r)\varphi_i(r)\varphi_j(r)\right\rangle\right)=\left\langle
S_0(r)\varphi_j(r)\right\rangle
\end{eqnarray}
and the boundary integration $\langle\phi_k(r)\varphi_j(r)\rangle_{\partial\Omega}=\frac{A}{6}$ if $k=j$ 
or $\frac{A}{12}$ if $k\ne j$. The general expression for the surface element integration is
\begin{equation}
\int_{\Omega_e}\varphi_1^l\varphi_2^m\varphi_3^ndr=\frac{l!m!n!}{(l+m+n+2)!}2A_e
\end{equation}


Equation (\ref{eq:weakformpartial}) is the final formula used in
the reconstruction code.

\section{Widefield illumination}
When a light source is located outside of the domain, i.e. a wide-field source, at the boundary, the 
in-bound flux, $J^-$ is no longer zero, thus, the boundary condition should be modified as
\begin{equation}
\label{eq:widefieldsrc} \left(\boldsymbol\Phi-l_s\nabla{\boldsymbol\Phi}\right)|_{\partial\Omega}=J^-
\end{equation}
Thus, we should replace $\nabla{\boldsymbol\Phi}$ by $(\boldsymbol\Phi - J^-)/l_s$ in (\ref{eq:weakform}), and we get
\begin{eqnarray}
\label{eq:weakformpartial}\nonumber&&\sum_{i}\Phi_i\left(\sum_{k}D_k\left\langle
\phi_k(r)\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle -
\sum_{k}\frac{1-R_{eff}}{2(1+R_{eff})}\left\langle
\phi_k(r)\varphi_j(r)\right\rangle_{\partial\Omega}\right.
\\&&\left.+\sum_{k}\left(\mu_{a}^k+\frac{j\omega}{c}\right)\left\langle\phi_k(r)\varphi_i(r)\varphi_j(r)\right\rangle\right)=\\
&& \left\langle S_0(r)\varphi_j(r)\right\rangle + \sum_kJ^-_kD_k/l_s\langle\phi_k(r)\varphi_j(r)\rangle|_{\partial\Omega}
\end{eqnarray}

\section{Nodal adjoint method for computing the Jacobian matrix}
Equation (\ref{eq:weakformpartial}) for each element written in the
matrix form is
\begin{equation}
\label{eq:matrix} \mathbf{A}{\boldsymbol\Phi}=\mathbf{b}
\end{equation}
where the element of $\mathbf{A}$ can be written as
\begin{equation}
a_{i,j}=\sum_{i,j\in e}\left(\sum_{k=1}^4D_k\left\langle
\phi_k(r)\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle+\sum_{k=1}^4\left(\mu_{a}^k+\frac{j\omega}{c}\right)\left\langle\phi_k(r)\varphi_i(r)\varphi_j(r)\right\rangle\right)
\end{equation}
for node-based FEM matrix, and 
\begin{equation}
a_{i,j}=\sum_{i,j\in e}\left(D_e\left\langle
\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle+\left(\mu_{a}^{e}+\frac{j\omega}{c}\right)\left\langle\varphi_i(r)\varphi_j(r)\right\rangle\right)
\end{equation}
for element-based FEM matrix, where $D_e$ and $\mu_a^{e}$ are the diffusion and absorption coefficients 
of the $e$-th element.

Taking derivative of $\mu_a$ and $D$ on both sides of
(\ref{eq:matrix}), we get
\begin{eqnarray}
\frac{\partial \mathbf{A}}{\partial
\mu_a^k}{\boldsymbol\Phi}&=&-\mathbf{A}\frac{\partial {\boldsymbol\Phi}}{\partial \mu_a^k}\\
\frac{\partial \mathbf{A}}{\partial
D_k}{\boldsymbol\Phi}&=&-\mathbf{A}\frac{\boldsymbol\Phi}{\partial
D^k}
\end{eqnarray}
where the elements of matrix $\frac{\partial \mathbf{A}}{\partial
\mu_a^k}$ and $\frac{\partial \mathbf{A}}{\partial D_k}$ are
written as
\begin{eqnarray}
k_{i,j}=\frac{\partial
a_{i,j}}{\mu_a^k}&=&\left\langle\phi_k(r)\varphi_i(r)\varphi_j(r)\right\rangle\\
h_{i,j}=\frac{\partial a_{i,j}}{D_k}&=&\left\langle
\phi_k(r)\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle
\end{eqnarray}
respectively.

Similarly, taking derivatives with respect to $D_e$ and $\mu_a^{e}$ yields
\begin{eqnarray}
k_{i,j}=\frac{\partial
a_{i,j}}{\mu_a^{e}}&=&\left\langle\varphi_i(r)\varphi_j(r)\right\rangle\\
h_{i,j}=\frac{\partial a_{i,j}}{D_e}&=&\left\langle
\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle
\end{eqnarray}
respectively. Again, the angular bracket $\langle\cdot\rangle$ represents volume integration inside
an element, i.e.
\begin{equation}
\langle f(r) \rangle=\int_{\Omega_e}f(r)dr
\end{equation}
Therefore, the deriatives to $\mu_a^e$ can be calculated using this relationship:
\begin{equation}
\int_{\Omega_e}\varphi_1^l\varphi_2^m\varphi_3^n\varphi_4^sdr=\frac{l!m!n!s!}{(l+m+n+s+3)!}6V_e
\end{equation}
As a result, for element-based $\mu_a$ Jacobian, we have
\begin{equation}
\mathbf{K}_e=\left[\frac{\partial a_{i,j}}{\mu_a^{e}}\right]=\sum_e\frac{V_e}{20}\left(
\begin{array}{llll}
2 & 1 & 1 & 1 \\
1 & 2 & 1 & 1 \\
1 & 1 & 2 & 1 \\
1 & 1 & 1 & 2 
\end{array}
\right)
\end{equation}
For element-based $D$ Jacobian, we have
\begin{eqnarray}\nonumber
\mathbf{H}_e&=&\left[\frac{\partial
a_{i,j}}{D^{e}}\right]\\&=&\sum_e\frac{1}{(6V_e)^2}\left(
\begin{array}{llll}
a_1a_1+b_1b_1+c_1c_1 & a_1a_2+b_1b_2+c_1c_2 & a_1a_3+b_1b_3+c_1c_3 & a_1a_4+b_1b_4+c_1c_4 \\
a_2a_1+b_2b_1+c_2c_1 & a_2a_2+b_2b_2+c_2c_2 & a_2a_3+b_2b_3+c_2c_3 & a_2a_4+b_2b_4+c_2c_4 \\
a_3a_1+b_3b_1+c_3c_1 & a_3a_2+b_3b_2+c_3c_2 & a_3a_3+b_3b_3+c_3c_3 & a_3a_4+b_3b_4+c_3c_4 \\
a_4a_1+b_4b_1+c_4c_1 & a_4a_2+b_4b_2+c_4c_2 & a_4a_3+b_4b_3+c_4c_3 & a_4a_4+b_4b_4+c_4c_4
\end{array}
\right)
\end{eqnarray}
where $a_i,b_i,c_i (i=1,4)$ are the linear coefficients in the expressions of the basis functions $\varphi_i$ as
\begin{equation}
\left(
\begin{array}{l}
\varphi_1\\
\varphi_2\\
\varphi_3\\
\varphi_4
\end{array}\right)=\frac{1}{6V_e}\left(
\begin{array}{llll}
6V_{01} & a_1 & b_1 & c_1 \\
6V_{02} & a_2 & b_2 & c_2 \\
6V_{03} & a_3 & b_3 & c_3 \\
6V_{04} & a_4 & b_4 & c_4 
\end{array}
\right)\left(
\begin{array}{l}
1\\
x\\
y\\
z
\end{array}\right)
\end{equation}
Here, the matrix is the inversion of the Jacobian ($\mathbf{J}_e$) of the tetrahedron, i.e.
\begin{equation}\frac{1}{6V_e}
\left(
\begin{array}{llll}
6V_{01} & a_1 & b_1 & c_1 \\
6V_{02} & a_2 & b_2 & c_2 \\
6V_{03} & a_3 & b_3 & c_3 \\
6V_{04} & a_4 & b_4 & c_4 
\end{array}
\right)=\textrm{inv}\left(
\begin{array}{llll}
1 & 1 & 1 & 1 \\
x_1 & x_2 & x_3 & x_4 \\
y_1 & y_2 & y_3 & y_4 \\
z_1 & z_2 & z_3 & z_4  
\end{array}
\right)=\textrm{inv}(\mathbf{J}_e)
\end{equation}
and the element volume $V_e=\textrm{det}(\mathbf{J}_e)/6$

For the Jacobian corresponding to the $\mu_a^e$ in the $e$-th element, the 
element-based Jacobian can be expressed as
\begin{eqnarray}\nonumber
J_{\mu_a}((s,r),e)&=&\frac{\partial
\Phi_{s,r}}{d\mu_a^{e}}\\&=&\left(\mathbf{K}_e{\boldsymbol\Phi}_s^e\right)^T{\boldsymbol\Phi}_r^e
\end{eqnarray}
where vectors $\boldsymbol{\Phi}_s^e$ and $\boldsymbol{\Phi}_r^e$
are defined by $\boldsymbol{\Phi}_s^e=\{\Phi_s(\vec{p}_k)\}$,
$\boldsymbol{\Phi}_r^e=\{\Phi_r(\vec{p}_k)\}, (k=1,2,3,4)$,
respectively.

For the Jacobian corresponding to the $\tau$-th $\mu_a$, the nodal
adjoint method was applied and the Jacobian for the measurement at
the $r$-th detector with illumination of the $s$-th source is
expressed as
\begin{eqnarray}\nonumber
J_{\mu_a}((s,r),\tau)&=&\frac{\partial
\Phi_{s,r}}{d\mu_a^{\tau}}\\&=&\sum_{n\in\Omega_\tau}\left(\frac{\sum_{e\in\Omega_n}
V_e}{4}\right)\phi(\vec{p}_n)\Phi_s(\vec{p}_n)\Phi_r(\vec{p}_n)
\end{eqnarray}
where $\sum_{n\in\Omega_\tau}$ refers to the summation over the
forward nodes which fall inside $\Omega_\tau$ and
$\sum_{e\in\Omega_n}$ refers to the summation over the forward
elements that share the n-th forward node; $V_e$ is the volume of
the element, $\Phi_s(\vec{p}_n)$ and $\Phi_r(\vec{p}_n)$ are the
forward field and the adjoint field at the n-th forward node,
respectively.

The Jacobian with respect to $D$ is computed by the traditional
element-based adjoint method, which is written as

\begin{eqnarray}\nonumber
J_{D}((s,r),\tau)&=&\frac{\partial
\Phi_{s,r}}{dD_{\tau}}\\&=&\sum_{e\in\Omega_\tau}\left(\mathbf{H}_\tau^e{\boldsymbol\Phi}_s^e\right)^T{\boldsymbol\Phi}_r^e
\end{eqnarray}
where $\sum_{e\in\Omega_\tau}$ denotes the summation over all
elements that in the immediate vicinity of the parameter node $\tau$;
matrix $\mathbf{H}_\tau^e$ has form of
\begin{equation}
h^\tau_{i,j}=\left\langle
\phi_\tau(r)\nabla\varphi_i(r)\cdot\nabla\varphi_j(r)\right\rangle
\end{equation}

\section{Log-magnitude/unwrapped phase form of Jacobian}
The complex form update equation is
\begin{equation}
\left(\mathbf{J}_{\mu_a},\mathbf{J}_{D}\right)\left(
\begin{array}{l}
\boldsymbol{\Delta}{\mu_a}\\
\boldsymbol{\Delta}{D}
\end{array}
\right)=\boldsymbol{\Phi}^{meas}-\boldsymbol{\Phi}^{calc}
\end{equation}
or
\begin{equation}
\label{eq:cpxupdate} \mathbf{J}\left(
\begin{array}{l}
\boldsymbol{\Delta}{\mu_a}\\
\boldsymbol{\Delta}{D}
\end{array}
\right)=\boldsymbol{\Delta}\boldsymbol{\Phi}
\end{equation}
then the real-form of (\ref{eq:cpxupdate}) is given by
\begin{equation}
\label{eq:realupdate}\left(
\begin{array}{l}
\Re e{\mathbf{J}}_{\mu_a}, \Im m{\mathbf{J}}_{D}
\end{array}\right)
\left(
\begin{array}{l}
\boldsymbol{\Delta}{\mu_a}\\
\boldsymbol{\Delta}{D}
\end{array}
\right)=\left(
\begin{array}{l}
\Re e{\boldsymbol{\Delta\Phi}}\\
\Im m{\boldsymbol{\Delta\Phi}}
\end{array}\right)
\end{equation}

\end{document}
