Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

Intensity-based iterative reconstruction for helical grating interferometry breast CT with static grating configuration

Jinqiu Xu, Zhentian Wang, Stefano van Gogh, Michał Rawlik, Simon Spindler, and Marco Stampanoni

Open Access Open Access

Abstract

Grating interferometry breast computed tomography (GI-BCT) has the potential to provide enhanced soft tissue contrast and to improve visualization of cancerous lesions for breast imaging. However, with a conventional scanning protocol, a GI-BCT scan requires longer scanning time and higher operation complexity compared to conventional attenuation-based CT. This is mainly due to multiple grating movements at every projection angle, so-called phase stepping, which is used to retrieve attenuation, phase, and scattering (dark-field) signals. To reduce the measurement time and complexity and extend the field of view, we have adopted a helical GI-CT setup and present here the corresponding tomographic reconstruction algorithm. This method allows simultaneous reconstruction of attenuation, phase contrast, and scattering images while avoiding grating movements. Experiments on simulated phantom and real initial intensity, visibility and phase maps are provided to validate our method.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Breast cancer is the most diagnosed malignancy and leading cancer-related cause of death in women [1]. Early detection of this disease has been known to help save women's lives [2]. Currently, 2D digital mammography (DM), digital breast tomosynthesis (DBT), absorption breast computed tomography (CT), ultrasound imaging, and breast magnetic resonance imaging (MRI) all provide tools of breast diagnosis. Nevertheless, these methods have some limitations, such as patient discomfort associated with breast compression and diagnostic performance (sensitivity and specificity) [1]. Grating interferometry CT (GI-CT) based phase contrast and scattering imaging are emerging techniques that may overcome some challenges of current modalities [3–5]. X-ray phase imaging usually provides higher sensitivity to variations within soft tissue, since for a given resolution, the difference between different tissues induced by X-ray refraction can be larger than its absorption counterpart. In addition, scattering imaging can provide complementary information on highly scattering structures in the breast [6,7], such as micro-calcifications. GI-CT has been shown to cope well with the relatively broad spectra provided by conventional X-ray sources and allows to retrieve 3D attenuation, phase, and scattering signals, therefore baring a strong potential for medical applications [8,9].

A typical GI-CT system consists of three gratings: G0, G1, and G2 (Fig. 1(a)). During data acquisition, one of the gratings is usually translated at least three times per angular projection to record a series of interferometric patterns [10]. This procedure is called phase stepping (PS) and generates an intensity curve for each detector pixel which can be used later to retrieve the aforementioned three signals [11,12]. Obviously, the PS approach increases scanning time and system complexity, something we want to address with our present work.

 figure: Fig. 1.

Fig. 1. Schematics of helical GI-BCT. (a) The setup of GI-BCT consists of X-ray source, a source grating G0, a beam-splitter grating G1 and an analyzer grating G2. (b) The geometry of a helical scan. The imaging part moves along the $z - $ direction from bottom to top while it rotates around the sample. S is the X-ray source, R denotes the distance between S and the rotation axis ($z - $ axis), D represents the distance between the rotation axis and the detector, and h is the height of the upward movement of S after one rotation.

Download Full Size | PDF

Several methods have been developed to mitigate this problem. Zhu et al. proposed a so-called reverse projection (RP) method which sets the beam-splitter (or the analyzer) grating at a specific position and exploits the relationship between the projected image at the rotation angle $\varphi $ and its corresponding reverse image at $\varphi + {180^ \circ }$, and allows to extract both attenuation and phase signals [13] directly. This method would eliminate the need for PS at each projection angle, but the phase signal can only be retrieved correctly for small refraction angles while the scattering signal cannot be obtained simultaneously [14]. Both Ritter et al. and Teuffenbach et al. presented solutions based directly on the measured intensity on the detector to simultaneously reconstruct attenuation, phase contrast and scattering without PS and with no intermediate signal retrieval [15–17]. Ritter et al. performed a maximum likelihood reconstruction with gradient descent as the solver [15,16] while Teuffenbach et al. presented a single-shot interferogram-based method relying on a statistical iterative reconstruction approach (IB-SIR) [17,18]. The authors tested two single-shot data acquisition methods: a sliding window acquisition pattern (SW-IBSIR) and a single-shot version of the electromagnetic phase-stepping (SSEPS-IBSIR) and demonstrated that with a suitable acquisition pattern it is possible to simultaneously reconstruct the three object’s properties [17].

The aforementioned approaches were developed for conventional circular GI-CT reconstruction. Currently, however, the field of view (FOV) of GI-CT is mainly limited by the grating size and the size of the high-efficiency photon-counting detector. While it is feasible to stitch a few gratings in a row to expand the FOV horizontally and the same applies for the modular sensors, large area gratings and photon-counting detectors are very challenging and extremely expensive. Thus, a helical CT (rotation + translation) approach would allow continuous reconstruction with a large axial coverage, have moderate requirements to gratings and detector size, reduce examination time and provide isotropic 3D resolution [19]. To date, there are only a few studies about helical GI-CT reconstructions. Fu et al. applied PS to helical GI-CT [20] and Marschner et al. extended a method originally proposed by Kottler et al. [10]. The latter avoids PS during the helical scan but needs intermediate signals retrieval steps and requires each part of the sample to be measured completely at least three times to obtain enough information to retrieve attenuation, phase, and scattering signals [21].

Our group is building a GI breast CT (GI-BCT) prototype relying on helical geometry for the volumetric (3D) measurement of the whole breast in-vivo [22]. Specifically, in this work, we present a reconstruction approach where the analyzer grating has been tilted to form Moiré fringes [23]. Further, inspired by the work of Teuffenbach et al., we present an intensity-based iterative reconstruction method that simultaneously reconstructs attenuation, phase contrast and scattering with static gratings and without intermediate signal extraction. It was worth mentioning that scattering in this paper is considered isotropic as in previous works [24,25]. We evaluated the performance of our method for different pitches (a crucial parameter in helical scans [26]) and different noise levels. The simulated phantom and real initial intensity, visibility and phase maps were used in numerical experiments for validation purposes.

2. Method

2.1 GI-BCT setup

A schematic of a helical GI-BCT is provided in Fig. 1. It consists of a source grating G0, a beam-splitter grating G1, and an analyzer grating G2. To cope with the cone beam geometry of the system, the linear G0, G1 and G2 gratings are bent into a cylindrical shape and this avoids visibility losses away from the optical axis [27]. G0 is an absorbing grating used to transform the large focal spot into a series of line sources satisfying the spatial coherence requirement for an efficient operation of the interferometer with conventional X-ray laboratory sources [28]. G1 is a phase grating that imprints a periodic modulation to the wavefront, which in turn produces an interference pattern at fractional Talbot distances downstream [25]. G2 is a tilted absorbing grating used to generate the Moiré fringes and induce displacements of fringe locations for each detector line, which are crucial for signal recovery during the helical scan [23]. G2 is positioned at one particular fractional Talbot distance and converts the interference pattern created by G1 into intensity modulations recorded by the detector [1,29]. The amplitude, phase and visibility of the intensity modulation provide insights into the attenuation, phase contrast, and scattering properties of the sample.

During a helical scan, the X-ray source and the detector form a cone beam geometry and the imaging system moves along the z-direction while it rotates around the sample [19,30]. From the perspective of the sample, the X-ray source traces a helical trajectory, as shown by the blue line in Fig. 1(b).

2.2 Parametric formulation of the X-ray path and the model for helical scan imaging

Referring to the helical CT geometry shown in Fig. 1(b) let us assume that ${\theta _s}$ is the angle of starting position of X-ray source S, ${\theta _r}$ is the rotation angle of the imaging system (which rotates clockwise). Then the coordinate of S can be defined as $({R\cos \theta ,R\sin \theta ,{{{\theta_r}h} / {({2\pi } )}}} )$, in which $\theta = {\theta _s} + {\theta _r}$. Q is a detector element and its local coordinate in the detector plane is $({u,v} )$. Thus the coordinate of Q in the $oxyz$ coordinate system is $({ - D\cos \theta - u\sin \theta , - D\sin \theta + u\cos \theta ,v + {{{\theta_r}h} / {({2\pi } )}}} )$[31]. Finally, the ray $SQ({x(t ),y(t ),z(t )} )$ can be denoted as follows:

$$\left\{ \begin{array}{l} x(t )= R\cos \theta + t({ - D\cos \theta - u\sin \theta - R\cos \theta } )\\ y(t )= R\sin \theta + t({ - D\sin \theta + u\cos \theta - R\sin \theta } )\\ z(t )= {{{\theta_r}h} / {({2\pi } )+ tv}} \end{array} \right., $$
where $t \in ({0,1} )$.

Let i denote the index of the X-ray path on the detector pixel $({u,v} )$ at angle $\theta $, $i = 1,2,\ldots ,M$, where M is the total number of X-rays. ${I_i}$ denotes the detected intensity on the detector corresponding to the $i - \textrm{th}$ X-ray path. ${T_i}$, ${D_i}$ and ${\phi _i}$ model attenuation, scattering, and phase, respectively. $I_i^0$, $V_i^0$ and $\phi _i^0$ denote the expected mean, expected visibility and phase position of a blank scan without sample present. They are viewed as known values and set to be same at different projection angles in the following experiments. Then ${I_i}$ can be described as [17]

$${I_i} = I_i^0{T_i}[{1 + V_i^0{D_i}\cos ({\phi_i^0 - {\phi_i}} )} ]. $$
${T_i}$, ${D_i}$ and ${\phi _i}$ are independently related to the linear attenuation $\mu (\boldsymbol{s} )$, scattering property $\varepsilon (\boldsymbol{s} )$ and phase contrast $\varphi (\boldsymbol{s} )\;({\boldsymbol{s} = ({x,y,z} )} )$ of the object. Based on the formula of the ray of the helical cone beam geometry in Eq. (1), we thus have
$${T_i} = \textrm{exp} \left( { - \int_0^1 {\mu (\boldsymbol{s} )dt} } \right),\boldsymbol{s} \in {l_i}, $$
$${D_i} = \textrm{exp} \left( { - \int_0^1 {\varepsilon (\boldsymbol{s} )dt} } \right),\boldsymbol{s} \in {l_i}, $$
and
$${\phi _i} = {c_\phi }{\partial _x}\int_0^1 {\varphi (\boldsymbol{s} )dt} ,\boldsymbol{s} \in {l_i}. $$
Where ${l_i}$ denotes the integration path from the X-ray source S to the detector pixel $({u,v} )$ in the helical scan. $\varepsilon (\boldsymbol{s} )= {c_D}\boldsymbol{\zeta }(\boldsymbol{s} )$, $\boldsymbol{\zeta }(\boldsymbol{s} )$ represents the linear diffusion coefficient [29], ${c_D}$ and ${c_\phi }$ are constants accounting for geometrical factors that depend on the imaging system, and ${\partial _x}$ denotes the derivative perpendicular to the grating linear structures, i.e. along the $x - $ axis (see Fig. 1).

For further numerical processing, the above equations need to be discretized. With respect to the field of view, $\mu (\boldsymbol{s} )$, $\varepsilon (\boldsymbol{s} )$ and $\varphi (\boldsymbol{s} )$ are firstly discretized as three-dimensional digital images ${[{{\mu_{r,s,t}}} ]_{R \times S \times T}}$, ${[{{\varepsilon_{r,s,t}}} ]_{R \times S \times T}}$, ${[{{\varphi_{r,s,t}}} ]_{R \times S \times T}}$, while R, S, and T denote the number of the slice, height, and width of the 3D dataset, $r = 1,\ldots ,R$, $s = 1,\ldots ,S$, $t = 1,\ldots ,T$. Let ${\mu _j} = {\mu _{r,s,t}}$, ${\varepsilon _j} = {\varepsilon _{r,s,t}}$ and ${\varphi _j} = {\varphi _{r,s,t}}$, in which j denotes the index of image pixel and $j = (r-1) \times ({S \times T} )+ (s-1) \times T + t$. then Eqs. (3) and (4) can be approximated respectively as

$${T_i} = \boldsymbol{exp} \left( { - \sum\limits_{j = 1}^N {{A_{i,j}}} {\mu_j}} \right)$$
and
$${D_i} = \boldsymbol{exp} \left( { - \sum\limits_{j = 1}^N {{A_{i,j}}} {\varepsilon_j}} \right), $$
where N denotes the total number of image pixels and $N = R \times S \times T$. The coefficient ${A_{i,j}}$ is defined as the intersection length between the X-ray path ${l_i}$ and the $j - \textrm{th}$ pixel. Equation (5) can be approximated by using the finite difference between the integrations along two adjacent X-ray paths along the $u - $ direction. If we take the X-ray paths along one row ($u - $ direction) of the detector as an example, then Eq. (5) can be discretized as
$${\phi _i} = \left\{ \begin{array}{l} {{{c_\phi }\left( {\sum\limits_{j = 1}^N {{A_{i + 1,j}}{\varphi_j} - } \sum\limits_{j = 1}^N {{A_{i,j}}{\varphi_j}} } \right)} / d}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} kU \le i < ({k + 1} )\times U\\ 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i = ({k + 1} )\times U \end{array} \right., $$
where d is the detector pixel size, $k = 0,1,\ldots ,V\;({V = {M / {U - 1}}} )$, and U is the number of detector pixels along $u - $ direction of the detector. To simplify the following expressions, in vector and matrix form, we define $\boldsymbol{\mu } = ({{\mu_1},{\mu_2},\ldots ,{\mu_N}} )$, $\boldsymbol{\varepsilon } = ({{\varepsilon_1},{\varepsilon_2},\ldots ,{\varepsilon_N}} )$ and $\boldsymbol{\varphi } = ({{\varphi_1},{\varphi_2},\ldots ,{\varphi_N}} )$.

2.3 Optimization model

According to the imaging model (2), the reconstruction problem can be formulated as the following optimization problem

$$\arg {\min _{\boldsymbol{x} = ({\boldsymbol{\mu ,\varepsilon ,\varphi }} )}}\left\{ {\sum\limits_i {{{({I_i^{measure} - {I_i}} )}^2}} } \right\}, $$
where $I_i^{measure}$ denotes the measured intensity. Furthermore, we can add regularization terms on Eq. (9) to obtain a smooth solution. In this paper, total variation (TV) regularization is chosen because of its effective noise removal [32], yielding to the following optimization model:
$$\arg {\min _{\boldsymbol{x} = ({\boldsymbol{\mu ,\varepsilon ,\varphi }} )}}\left\{ {\sum\limits_i {{{({I_i^{measure} - {I_i}} )}^2}} + {\beta_{\boldsymbol{\mu}}}{{||\boldsymbol{\mu } ||}_{TV}} + {\beta_{\boldsymbol{\varepsilon}}}{{||\boldsymbol{\varepsilon } ||}_{TV}} + {\beta_{\boldsymbol{\varphi}}}{{||\boldsymbol{\varphi } ||}_{TV}}} \right\}, $$
with ${||\cdot ||_{TV}}$ we denote the TV of the image [33] with parameters ${\beta _{\boldsymbol{\mu}}}$, ${\beta _{\boldsymbol{\varepsilon}}}$ and ${\beta _{\boldsymbol{\varphi}}}$ controlling the regularization strengths.

We use an approach based on the alternating direction method of multipliers (ADMM) to solve the optimization model (10) [34]. Let ${\boldsymbol{x}^{(k )}} = ({{\boldsymbol{\mu }^{(k )}},{\boldsymbol{\varepsilon }^{(k )}},{\boldsymbol{\varphi }^{(k )}}} )$ denote the $k - \textrm{th}$ iterations and ${\boldsymbol{x}^{(0 )}} = ({\boldsymbol{0},\boldsymbol{0},\boldsymbol{0}} )$, then the update for ${\boldsymbol{x}^{({k + 1} )}}$ can be split into two steps. Firstly, we evaluate the intensity model

$${\boldsymbol{x}^{({k + {1 / 2}} )}} = \arg {\min\nolimits _{\boldsymbol{x} = ({\boldsymbol{\mu ,\varepsilon ,\varphi }} )}}\left\{ {\sum\limits_i {{{({I_i^{measure} - {I_i}} )}^2}} } \right\}, $$
and minimize it with the classical L-BFGS-B algorithm [17,35]. Then we evaluate the following three quantities:
$${\boldsymbol{\mu }^{({k + 1} )}} = \arg {\min\nolimits _{\boldsymbol{\mu}}}\{{||{\boldsymbol{\mu } - {\boldsymbol{\mu }^{({k + {1 / 2}} )}}} ||_2^2 + {\beta_{\boldsymbol{\mu}}}{{||\boldsymbol{\mu } ||}_{TV}}} \}. $$
$${\boldsymbol{\varepsilon }^{({k + 1} )}} = \arg {\min\nolimits _{\boldsymbol{\varepsilon}}}\{{||{\boldsymbol{\varepsilon } - {\boldsymbol{\varepsilon }^{({k + {1 / 2}} )}}} ||_2^2 + {\beta_{\boldsymbol{\varepsilon}}}{{||\boldsymbol{\varepsilon } ||}_{TV}}} \}. $$
$${\boldsymbol{\varphi }^{({k + 1} )}} = \arg {\min\nolimits _{\boldsymbol{\varphi}}}\{{||{\boldsymbol{\varphi } - {\boldsymbol{\varphi }^{({k + {1 / 2}} )}}} ||_2^2 + {\beta_{\boldsymbol{\varphi}}}{{||\boldsymbol{\varphi } ||}_{TV}}} \}. $$
$${\boldsymbol{x}^{({k + 1} )}} = ({{\boldsymbol{\mu }^{({k + 1} )}},{\boldsymbol{\varepsilon }^{({k + 1} )}},{\boldsymbol{\varphi }^{({k + 1} )}}} ). $$

The optimization problems (12), (13), and (14) were finally solved by Chambolle’s algorithm [36]. For the optimization problem (11), we mainly present the derivative expression of the object function $f = \sum\limits_i {{{({I_i^{measure} - {I_i}} )}^2}}$ on $\boldsymbol{x}$, then it is easy to solve ${\boldsymbol{x}^{({k + {1 / 2}} )}}$ according to the solving process of L-BFGS-B algorithm [35]. The derivative can be expressed as

$$\frac{{\partial f}}{{\partial \boldsymbol{x}}} = \left( {\frac{{\partial f}}{{\partial \boldsymbol{\mu }}},\frac{{\partial f}}{{\partial \boldsymbol{\varepsilon }}},\frac{{\partial f}}{{\partial \boldsymbol{\varphi }}}} \right). $$

Take $\frac{{\partial f}}{{\partial \boldsymbol{\varphi }}}$ as the example to illustrate the derivative and it is easy to extend to $\frac{{\partial f}}{{\partial \boldsymbol{\mu }}}$ and $\frac{{\partial f}}{{\partial \boldsymbol{\varepsilon }}}$. Let ${f_i} = {({I_i^{measure} - {I_i}} )^2}$, we have

$$\begin{array}{l} \frac{{\partial f}}{{\partial \boldsymbol{\varphi }}} = \sum\limits_i {\frac{{\partial {f_i}}}{{\partial \boldsymbol{\varphi }}}} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_i {\frac{{\partial {f_i}}}{{\partial {\phi _i}}} \cdot \frac{{\partial {\phi _i}}}{{\partial \boldsymbol{\varphi }}}} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = {\kern 1pt} {\kern 1pt} {\left( {\left[ {{{\left( {\frac{{\partial {\phi_1}}}{{\partial \boldsymbol{\varphi }}}} \right)}^T},{{\left( {\frac{{\partial {\phi_2}}}{{\partial \boldsymbol{\varphi }}}} \right)}^T},\ldots ,{{\left( {\frac{{\partial {\phi_{M - 1}}}}{{\partial \boldsymbol{\varphi }}}} \right)}^T}} \right] \cdot {{\left[ {\frac{{\partial {f_1}}}{{\partial {\phi_1}}},\frac{{\partial {f_2}}}{{\partial {\phi_2}}},\ldots ,\frac{{\partial {f_{M - 1}}}}{{\partial {\phi_{M - 1}}}}} \right]}^T}} \right)^T} \end{array}$$
where $\frac{{\partial {f_i}}}{{\partial {\phi _i}}} = 2({{I_i} - I_i^{measure}} )I_i^0{T_i}V_i^0{D_i}\sin ({\phi_i^0 - {\phi_i}} )$. Let $\boldsymbol{P}$ denote ${\left[ {\frac{{\partial {f_1}}}{{\partial {\phi_1}}},\frac{{\partial {f_2}}}{{\partial {\phi_2}}},\ldots ,\frac{{\partial {f_{M - 1}}}}{{\partial {\phi_{M - 1}}}}} \right]^T}$. Based on the definition of ${\phi _i}$ in Eq. (8), then Eq. (17) can be transformed to
$${\left( {\frac{{\partial f}}{{\partial \boldsymbol{\varphi }}}} \right)^T} = \frac{{{c_\phi }}}{d}{\boldsymbol{W}^T}\boldsymbol{P}, $$
where ${\boldsymbol{W}^T} = \left[ \begin{array}{l} {A_{2,1}} - {A_{1,1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{3,1}} - {A_{2,1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {A_{M,1}} - {A_{M - 1,1}}\\ {A_{2,2}} - {A_{1,2}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{3,2}} - {A_{2,2}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {A_{M,2}} - {A_{M - 1,2}}{\kern 1pt} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} \\ {A_{2,N}} - {A_{1,N}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{3,N}} - {A_{2,N}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {A_{M,N}} - {A_{M - 1,N}}{\kern 1pt} \end{array} \right]$.

Then the matrix ${\boldsymbol{W}^T}$ is split into the minus of two matrices $\boldsymbol{W}_1^T$ and $\boldsymbol{W}_2^T$:

$$\boldsymbol{W}_1^T = \left[ \begin{array}{l} {A_{2,1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{3,1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {A_{M,1}}\\ {A_{2,2}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{3,2}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {A_{M,2}}{\kern 1pt} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} \\ {A_{2,N}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{3,N}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {A_{M,N}}{\kern 1pt} \end{array} \right]\;\textrm{and}\; \boldsymbol{W}_2^T = \left[ \begin{array}{l} {A_{1,1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{2,1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {A_{M - 1,1}}\\ {A_{1,2}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{2,2}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {A_{M - 1,2}}{\kern 1pt} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} \\ {A_{1,N}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A_{2,N}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {A_{M - 1,N}}{\kern 1pt} \end{array} \right],$$
thus
$${\left( {\frac{{\partial f}}{{\partial \boldsymbol{\varphi }}}} \right)^T} = \frac{{{c_\phi }}}{d}({\boldsymbol{W}_1^T\boldsymbol{P} - \boldsymbol{W}_2^T\boldsymbol{P}} ). $$

Please note that the elements of ${\boldsymbol{W}_1}$ and ${\boldsymbol{W}_2}$ in the $({k + 1} )\times U$ rows are 0, where $k = 0,1,\ldots ,V\;({V = {M / {U - 2}}} )$. Let ${\overline {\boldsymbol{W}} _1}$, ${\overline {\boldsymbol{W}} _2}$ and $\overline {\boldsymbol{P}}$ represent ${\boldsymbol{W}_1}$, ${\boldsymbol{W}_2}$ and $\boldsymbol{P}$ with the $({k + 1} )\times U$ rows removed, respectively. ${\overline {\boldsymbol{W}} _1}$, ${\overline {\boldsymbol{W}} _2}$ and $\overline {\boldsymbol{P}}$ satisfy

$${\left( {\frac{{\partial f}}{{\partial \boldsymbol{\varphi }}}} \right)^T} = \frac{{{c_\phi }}}{d}({\overline {\boldsymbol{W}}_1^T\overline {\boldsymbol{P}} - \overline {\boldsymbol{W}}_2^T\overline {\boldsymbol{P}} } ). $$

Both ${\overline {\boldsymbol{W}} _1}$ and ${\overline {\boldsymbol{W}} _2}$ can be viewed as the system matrixes with $U - 1$ detector pixels along $u - $ direction of the detector and differ from each other by one pixel shift along $u - $ direction. Thus $\overline {\boldsymbol{W}} _1^T\overline {\boldsymbol{P}}$ and $\overline {\boldsymbol{W}} _2^T\overline {\boldsymbol{P}}$ can be regarded as the back projection of data $\overline {\boldsymbol{P}}$. Finally, $\frac{{\partial f}}{{\partial \boldsymbol{\varphi }}}$ can be calculated directly by twice back projection.

The proposed algorithm is summarized below. Note that the L-BFGS-B and Chambolle’s algorithms use Python’s Scipy and Scikit-image libraries.

Algorithm Process of the proposed iterative algorithm
Input: initial guess ${x^{(0 )}} = ({\boldsymbol{0},\boldsymbol{0},\boldsymbol{0}} )$; measured intensity ${I^{measure}} = [{I_i^{measure}} ]$; maximum number of iteration $K = 3000$; TV regularization parameters ${\beta _{\boldsymbol{\mu}}}$, ${\beta _{\boldsymbol{\varepsilon}}}$ and ${\beta _{\boldsymbol{\varphi}}}$; parameters of L-BFGS-B algorithm: bounds (the value range of iteration $\boldsymbol{x}$), ${maxcor = }{10}$ (the maximum number of variable metric corrections used to define the limited memory matrix), ${maxls} = 10$ (maximum number of line search steps).
Output: ${\boldsymbol{x}^{(K )}} = ({{\boldsymbol{\mu }^{(K )}},{\boldsymbol{\varepsilon }^{(K )}},{\boldsymbol{\varphi }^{(K )}}} )$
Main iteration steps:
For $k = 0:K$ do
      Step 1:
      Use ${\boldsymbol{x}^{(k )}}$ as the initial solver to solve the optimization problem by using L-BFGS-B algorithm:
      ${\boldsymbol{x}^{({k + {1 / 2}} )}} = \arg {\min\nolimits _{\boldsymbol{x} = ({\boldsymbol{\mu ,\varepsilon ,\varphi }} )}}\left\{ {\sum\limits_i {{{({I_i^{measure} - {I_i}} )}^2}} } \right\}$
      Step 2:
      Solving optimization problems by using Chambolle’s algorithm:
      ${\boldsymbol{\mu }^{({k + 1} )}} = \arg {\min\nolimits _{\boldsymbol{\mu}}}\{{||{\boldsymbol{\mu } - {\boldsymbol{\mu }^{({k + {1 / 2}} )}}} ||_2^2 + {\beta_{\boldsymbol{\mu}}}{{||\boldsymbol{\mu } ||}_{TV}}} \}$
      ${\boldsymbol{\varepsilon }^{({k + 1} )}} = \arg {\min\nolimits _{\boldsymbol{\varepsilon}}}\{{||{\boldsymbol{\varepsilon } - {\boldsymbol{\varepsilon }^{({k + {1 / 2}} )}}} ||_2^2 + {\beta_{\boldsymbol{\varepsilon}}}{{||\boldsymbol{\varepsilon } ||}_{TV}}} \}$
      ${\boldsymbol{\varphi }^{({k + 1} )}} = \arg {\min\nolimits _{\boldsymbol{\varphi}}}\{{||{\boldsymbol{\varphi } - {\boldsymbol{\varphi }^{({k + {1 / 2}} )}}} ||_2^2 + {\beta_{\boldsymbol{\varphi}}}{{||\boldsymbol{\varphi } ||}_{TV}}} \}$
      ${\boldsymbol{x}^{({k + 1} )}} = ({{\boldsymbol{\mu }^{({k + 1} )}},{\boldsymbol{\varepsilon }^{({k + 1} )}},{\boldsymbol{\varphi }^{({k + 1} )}}} )$
end for

3. Results

Experiments under ideal conditions (noise-free intensity data, ideal ${I^0}$, ${V^0}$ and ${\phi ^0}$ maps) were first carried out to demonstrate the effectiveness of the proposed algorithm. Then our method was tested with real maps to evaluate the performance with different pitches and noise levels.

Figure 2 shows a horizontal slice of the simulated breast phantom for the three contrast modalities, and the sagittal and coronal slice for attenuation and phase contrast along two blue dotted lines in horizontal slice. The attenuation and phase contrast both contain adipose, glandular, and skin tissues. The modeled tissue shapes had been derived from a 3D tomogram of a real patient, acquired with a breast CT scanner at the University Hospital Zürich under ethically approved patient informed consent [37,38]. The attenuation and scattering contain randomly generated dots with diameters of 3, 5, and 7 voxels (the voxel size is $40 \times 40 \times 40\mu {m^3}$) which represent the breast calcifications. The material of the calcifications was substituted by Aluminum because of the similar attenuation properties [39]. The phase shift coefficients $\varphi $ and attenuation coefficients $\mu $ of all tissues were calculated by the known relation $\varphi = {{2\pi \delta } / \lambda }$ and $\mu = {{4\pi \beta } / \lambda }$, where $\lambda $ is the X-ray wavelength at the energy of 45 keV, $\delta $ is the decrement of the refractive index, $\beta $ is the absorption index and the values of $\delta $ and $\beta $ were calculated using NIST XCOM and f1f2 Kissel.dat of the DABAX library based on the tissue definition in ICRU 46 [37,40]. The value of $\varepsilon $ were independently set as 6, 3.8 and 3.2 $\; $ for the calcifications with diameters of 3, 5, and 7 voxels, and the value unit is $c{m^{ - 1}}$.

 figure: Fig. 2.

Fig. 2. One horizontal slice of simulated breast phantom for three contrast modalities, and the coronal slice and sagittal slice for attenuation and phase contrast along two blue dotted lines in horizontal slice. From left to right: the slice of attenuation, phase contrast, and scattering phantom. Attenuation and phase contrast both contain adipose, glandular, and skin tissue equivalent. Attenuation and scattering contain calcifications as well. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

Download Full Size | PDF

Four regions of interests (ROI 1, ROI 2, ROI 3 and ROI 4) in both the attenuation and phase contrast images were used to calculate the quantitative metrics (contrast, contrast-to-noise ratio (CNR) and signal-to-noise (SNR) [41,42]) of the reconstructed results. Further, the Mean Square Error (MSE) and Structural Similarity (SSIM) metrics were used to assess the quality of the reconstructed images, with the phantom as a reference image. Taking ROI 1 and ROI 2 as an example, the equations of the five metrics are as follows:

$$\textrm{Contrast} = 20{\log _{10}}\left( {\frac{{{m_{\textrm{ROI}\;{2}}}}}{{{m_{\textrm{ROI}\;{1}}}}}} \right), $$
$$\textrm{CNR} = \frac{{|{{m_{\textrm{ROI}\;{2}}} - {m_{\textrm{ROI}\;{1}}}} |}}{{\sqrt {\sigma _{\textrm{ROI}\;{2}}^2 + \sigma _{\textrm{ROI}\;{1}}^2} }}, $$
$$\textrm{SNR} = \frac{{{m_{\textrm{ROI}\;{2}}}}}{{{\sigma _{\textrm{ROI}\;{1}}}}}, $$
$$\textrm{MSE} = \frac{1}{{R \times S \times T}}\sum\limits_{r = 1}^R {\sum\limits_{s = 1}^S {\sum\limits_{t = 1}^T {{{({{g_{r,s,t}} - {{\overline g }_{r,s,t}}} )}^2}} } }, $$
$$\textrm{SSIM}({g,\overline g } )= \frac{{({2{m_g}{m_{\overline g }} + {C_1}} )({2{\sigma_{g\overline g }} + {C_2}} )}}{{({m_g^2 + m_{\overline g }^2 + {C_1}} )({\sigma_g^2 + \sigma_{\overline g }^2 + {C_2}} )}}, $$
where the unit of Contrast is the decibel (dB); symbol ${m_{({\cdot} )}}$ denotes the mean value of corresponding signals; ${\sigma _{({\cdot} )}}$ means the standard deviations of signal amplitudes; g and $\overline g $ are the reference image and the test image, respectively; ${C_1}$ and ${C_2}$ are the constants that are used to prevent the denominator from being 0. The contrast metrics (in Fig. 2) between ROI 1 and ROI 2 are 2.52 (attenuation) and 1.37 (phase contrast), and the contrast metrics between ROI 3 and ROI 4 are 1.90 (attenuation) and 0.91 (phase contrast), respectively, which will be used as the quantitative reference value for the following experimental results.

Table 1 lists the experimental parameters used in the following simulation studies. In the simulated experiments, ${I^0}$, ${V^0}$ and ${\phi ^0}$ were taken as known values while the measured intensity ${{I}^{measure}}({{{I}^{measure}} = [{I_i^{measure}} ]} )$ was simulated by Eq. (2) considering the experimental parameters listed in Table 1. The linear integrals in Eqs. (3), (4) and (5) were generated by the Astra Toolbox [43]. The parameters bounds (mentioned in section 2.3) of attenuation, phase contrast and scattering were set to $({0,3} )$, $({0,2000} )$ and $({0,30} )$ respectively. The bounds were not set closed to the value range of the object properties and only provided rough constraints.

Tables Icon

Table 1. Experimental parameters.

View Table | View all tables in this article

3.1 Simulated maps

For a first benchmark, we implemented an ideal (noise-free) scenario with pitch 1 and compared with the standard PS method. The pitch, is defined as the ratio of h and the detector collimation at the rotation axis, is a crucial parameter in helical scans and affects both reconstructed image quality and scanning time [26]. ${I^0}$ and ${V^0}$ have been set as constant matrices. Further, each row of ${\phi ^0}$ is constant and, because G2 is tilted, a horizontal Moiré fringe is generated on the detector. It follows that each row of the detector (along z-axis), “sees” a different step of the traditional phase stepping curve [23]. The TV regularization parameters ${\beta _{\boldsymbol{\mu}}}$, ${\beta _{\boldsymbol{\varphi}}}$ and ${\beta _{\boldsymbol{\varepsilon}}}$ were determined manually by trial-and-error method. In this experiment, they were set to 0.005, 0.006 and 0.003, respectively. TV regularization term helps the reconstructed image to be smooth. For the PS method, three phase steps were adopted and the initial phase maps were set to ${\pi / 2}$, ${\pi / 2} + {{2\pi } / 3}$ and ${\pi / 2} + {{4\pi } / 3}$, respectively. For each phase step, the scanning geometry parameters were same as those of the proposed method. With the constant matrices ${I^0}$, ${V^0}$ and the three initial phase maps, three intensity data were generated to retrieve intermediate signals. Then Simultaneous Iterative Reconstruction Technique (SIRT) was used to reconstruct attenuation and scattering, Hilbert transform and back projection was adopted to recover phase contrast.

Figure 3 shows the reconstructed horizontal slice of the proposed method and PS method, from left to right, they are the attenuation $\mu $, two regions of interests (ROIs) of the attenuation and phase contrast, phase contrast $\varphi $ and scattering $\varepsilon $. PS method suffers from slight wrapping artifacts around the irregular outer edge in phase contrast but performs better for detail recovery as shown in both ROIs due to three times the intensity data of the proposed method, however, the overall image structures in attenuation and phase contrast, and all calcifications in attenuation and scattering can be recovered with the proposed method.

 figure: Fig. 3.

Fig. 3. Reconstructed results with simulated maps in noise-free case of the proposed method and PS method. The pitch is 1. From left to right: attenuation, ROIs of attenuation and phase contrast, phase contrast, and scattering. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

Download Full Size | PDF

Table 2 lists the MSEs, SSIMs, contrasts and relative errors of contrasts at the horizontal slice, with the standard deviations of MSEs, SSIMs and contrasts’ relative errors computed across 20 slices of the reconstructed volume. For both methods, the MSEs approach 0 and SSIMs are very close to 1. For the proposed method, the contrast metrics for attenuation and phase contrast are close to those of the phantom shown in Fig. 2.

Tables Icon

Table 2. MSEs, SSIMs, contrasts, contrasts’ relative errors of the images shown in Fig. 3 (simulated initial maps in noise-free case with pitch 1), and corresponding standard deviations computed across 20 slices (in the parentheses).

View Table | View all tables in this article

3.2 Real maps

As a first step towards implementing more realistic simulations, we acquired from our GI-BCT prototype the flat-field (i.e. without sample) PS curves by moving G0 to enable the retrieval of realistic values for intensity, phase and visibility maps $.$ Fig. 4 shows the figure of the initial maps. The maps are severely affected by the grating quality and alignment while the four vertical stripes in the maps correspond to dead areas on the detector. For the following experiments, we focused on one single fringe period, and therefore selected the 10 central rows from the three maps as being the ${I^0}$, ${V^0}$ and ${\phi ^0}$ for each projection angle.

 figure: Fig. 4.

Fig. 4. The real intensity, phase and visibility maps. The flat-field data had been measured with a CdTe, single-photon counting Dectris detector. The grating design energy was 46 keV, the tube was run at 70 kVp and 10 mA with an exposure time of 15 s.

Download Full Size | PDF

We studied the performance of our approach for different pitches and different noise levels. Intensive preparative work has shown that pitches of 1, 1.4 and 2 are very representative and were therefore chosen for the following simulations. The photon counting intrinsic features of the detector were simulated by generating Poisson-distributed noise from the intensity ${\boldsymbol{I}^{measure}}$ with average per-pixel count numbers of 8000, 4000 and 2000. For brevity, they were named as $\boldsymbol{I}_{({8000} )}^{measure}$, $\boldsymbol{I}_{({4000} )}^{measure}$ and $\boldsymbol{I}_{({2000} )}^{measure}$, respectively.

3.2.1 Effect of pitch variation

We first evaluated the performance of the proposed method with different pitches in the noise-free case. The TV regularization parameters ${\beta _{\boldsymbol{\mu}}}$, ${\beta _{\boldsymbol{\varphi}}}$ and ${\beta _{\boldsymbol{\varepsilon}}}$ were set to 0.005, 0.006 and 0.004. One reconstructed horizontal slice is shown in Fig. 5, from top to bottom, the pitches are 1, 1.4, and 2 respectively, from left to right are the attenuation, ROIs of attenuation and phase contrast, phase contrast and scattering, respectively. The image quality metrics MSE, SSIM and contrast at the central slice are listed in Table 3. A comparison with the measured values shown in Table 2 and the contrast values of our reference phantom allow us to state that our algorithm performs well with real maps when pitches are 1 and 1.4. If the pitch increases to 2, the edges in attenuation and phase contrast get severely blurred, as shown in the ROIs. However, all calcifications in attenuation and scattering object properties could be recovered at all pitches.

 figure: Fig. 5.

Fig. 5. Reconstructed results with different pitches in the noise-free case. From top to bottom, the pitches are 1, 1.4, and 2, respectively. From left to right: reconstructed attenuation, ROIs of attenuation and phase contrast, phase contrast, and scattering. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

Download Full Size | PDF

Tables Icon

Table 3. MSEs, SSIMs, contrasts, contrasts’ relative errors of the images shown in Fig. 5 (real initial maps in noise-free case with pitches 1, 1.4 and 2), and corresponding standard deviations computed across 20 slices (in the parentheses).

View Table | View all tables in this article

The reasons for increasing image degradation with increasing pitches are twofold: on the one hand, larger pitches mean fewer projection data over a fixed height of the scanned object leading to degradation of the reconstruction quality. On the other hand, since our method directly recovers $\mu $, $\varphi $ and $\varepsilon $ without any grating movement and intermediate signal retrieval, signal’s ($T$, $\phi $ and $D$) cross-talk during iteration cannot be excluded and likely leads to the blurring observed in the ROIs.

3.2.2 Behavior under different noise levels

Following the results presented in section 3.2.1, we have set the pitch to 1. To cope with the influence of noise on iterative results, the TV regularization parameters ${\beta _{\boldsymbol{\mu}}}$, ${\beta _{\boldsymbol{\varphi}}}$ and ${\beta _{\boldsymbol{\varepsilon}}}$ were increased compared with those in section 3.2.1. They were set to 0.007, 0.009 and 0.006, respectively. The MSEs, SSIMs, contrasts, contrasts’ relative errors, CNRs and SNRs at the reconstructed horizontal slice are listed in Table 4. Figure 6(a) displays horizontal slice of the reconstructed results. The sagittal and coronal slice of the iterative results with intensity $\boldsymbol{I}_{({4000} )}^{measure}$ are presented in Fig. 6(b). Figure 6(a) shows all calcifications in both attenuation and scattering reconstructions are visible. From top to bottom, as the noise level increases, blurring at the interfaces is more and more obvious as shown in ROIs, especially for the highest noise level, the blurring is significant. However, the reconstructed images indicate that the overall image quality of the attenuation and phase contrast channels with $\boldsymbol{I}_{({8000} )}^{measure}$ and $\boldsymbol{I}_{({4000} )}^{measure}$ are - as expected - not degraded; this is also supported by the SSIM metrics which are very close to 1; further, the SNRs and CNRs of the phase contrast channel are larger than those of attenuation.

 figure: Fig. 6.

Fig. 6. Reconstructed results with different noise levels. The pitch is 1. (a) From top to bottom, the intensity data are $\boldsymbol{I}_{({8000} )}^{measure}$, $\boldsymbol{I}_{({4000} )}^{measure}$ and $\boldsymbol{I}_{({2000} )}^{measure}$, respectively. From left to right: attenuation, ROIs of attenuation and phase contrast, phase contrast, and scattering. (b) Sagittal and coronal slice of the reconstructed attenuation and phase contrast in noisy case with intensity $\boldsymbol{I}_{({4000} )}^{measure}$ along two blue dotted lines in horizontal slice shown in Fig. 2. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

Download Full Size | PDF

Tables Icon

Table 4. MSEs, SSIMs, contrasts, contrasts’ relative errors, SNRs, CNRs of the images shown in Fig. 6 (real initial maps in noisy case with pitch 1), and corresponding standard deviations computed across 20 slices (in the parentheses).

View Table | View all tables in this article

4. Discussion

4.1 Convergence and running time

The experiments in section 3.2.1 and section 3.2.2 were used to explore the convergence behavior of the proposed method. The energy function (26) is computed to check how the method evolves with the iterations.

$$f= \sum\limits_i {{{({I_i^{measure} - I_i^0{T_i}[{1 + V_i^0{D_i}\cos ({\phi_i^0 - {\phi_i}} )} ]} )}^2}}. $$

The energy curves are plotted in Fig. 7. All curves in noise-free and noisy cases keep decreasing with the iteration numbers and reach stable states after 200 iterations. The global decreasing property of the energy curves suggests that the iterative image sequence is to minimize the energy function. The running time for each iteration step including step 1 and step 2 is 160 seconds. The computer configuration is Intel Xeon Gold 5222 CPU, 163GB of memory; Geforce RTX 2080 Ti independent graphics card, 24GB of memory.

 figure: Fig. 7.

Fig. 7. Energy curves of the proposed method with real initial maps for the noise-free and noisy cases. (a) Noise-free case with pitches 1, 1.4 and 2. (b) Noisy case with pitch 1.

Download Full Size | PDF

4.2 Necessity of tilting G2

Section 3 has shown that our algorithm can simultaneously reconstruct the attenuation, phase contrast, and scattering with static grating configuration under acceptable pitch and input noise levels in a helical GI-CT configuration. It is important to note that, tilting G2 is crucial for avoiding grating movements and significantly helps in dealing with the intrinsic under-determination nature of our optimization problem. In fact, from compressed sensing theory, it is known that the highest probability to correctly recover a signal in an underdetermined system can be achieved if the signal is sampled as nonuniformly as possible [17]. When G2 is tilted, the phase map contains phase shifts between adjacent rows, which increases the nonuniformity of the intensity data, making the retrieval of the three object properties easier. To demonstrate this statement, a simulated experiment with non-tilted grating was implemented. Based on the experiment in section 3.1, we set ${\phi ^0}$ to be a constant matrix while keeping all other reconstruction parameters the same as those in section 3.1. Figure 8 compares the reconstructed results with tilted G2 and non-tilted G2. It can be seen that none of the calcifications in scattering is recovered with non-tilted G2. Therefore, tilting G2 is crucial as it enables a proper signal retrieval, reduces the data acquisition complexity and ultimately the time of helical scan.

 figure: Fig. 8.

Fig. 8. Reconstructed results with tilted G2 and non-tilted G2 in the noise-free case (pitch is 1). From left to right: attenuation, phase contrast, and scattering. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

Download Full Size | PDF

5. Conclusion

In summary, this paper presents a method to simultaneously reconstruct attenuation, phase contrast and scattering for helical GI-BCT without intermediate signal retrieval and grating movements. The algorithm builds upon an optimization model which combines a forward imaging model for helical geometry with regularization terms. The optimization problem is solved iteratively via an alternating minimization technique yielding to a helical GI-BCT iterative reconstruction method. The experiments combining a simulated phantom and real flat-field data demonstrate the robustness of the method and its ability to simultaneously recover a complex, multichannel sample structure. Future work will investigate different optimization strategies and regularization methods to further increase the robustness of our method for higher pitch and noise levels, which would enable us to acquire volumetric scan faster.

Funding

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (Sinergia Grant Nr. CRSII5_183568); ETH-Research Commission (Grant Nr. ETH-12 20-2); PHRT-PIP (Project “CLARINET”); Tsinghua University (Initiative Scientific Research Program, No. 20213080039).

Acknowledgments

We thank Prof. Dr. Andreas Boss from the University Hospital Zürich for his precious advice on the phantom model.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. S. T. Taba, T. E. Gureyev, M. Alakhras, S. Lewis, D. Lockie, and P. C. Brennan, “X-Ray Phase-Contrast Technology in Breast Imaging: Principles, Options, and Clinical Application,” AJR. Am. J. Roentgenol. 211(1), 133–145 (2018). [CrossRef]  

2. F. Elmi, A. F. Movaghar, M. M. Elmi, H. Alinezhad, and N. Nikbakhsh, “Application of FT-IR spectroscopy on breast cancer serum analysis,” Spectrochim. Acta. A. Mol. Biomol. Spectrosc. 187, 87–91 (2017). [CrossRef]  

3. S. D. Auweter, J. Herzen, M. Willner, S. Grandl, K. Scherer, F. Bamberg, M. F. Reiser, F. Pfeiffer, and K. Hellerhoff, “X-ray phase-contrast imaging of the breast–advances towards clinical implementation,” Br. J. Radiol. 87(1034), 20130606 (2014). [CrossRef]  

4. M. Strobl, “General solution for quantitative dark-field contrast imaging with grating interferometers,” Sci. Rep. 4(1), 7243 (2015). [CrossRef]  

5. F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. F. Eikenberry, Ch Brönnimann, C. Grünzweig, and C. David, “Hard-X-ray dark-field imaging using a grating interferometer,” Nat. Mater. 7(2), 134–137 (2008). [CrossRef]  

6. T. Rauch, J. Rieger, G. Pelzer, F. Horn, R. Erber, M. Wunderle, J. Emons, N. Nabieva, N. Fuhrich, T. Michel, A. Hartmann, P. A. Fasching, and G. Anton, “Discrimination analysis of breast calcifications using x-ray dark-field radiography,” Med. Phys. 47(4), 1813–1826 (2020). [CrossRef]  

7. T. Michel, J. Rieger, G. Anton, F. Bayer, M. W. Beckmann, J. Durst, P. A. Fasching, W. Haas, A. Hartmann, G. Pelzer, M. Radicke, C. Rauh, A. Ritter, P. Sievers, R. Schulz-Wendtland, M. Uder, D. L. Wachter, T. Weber, E. Wenkel, and A. Zang, “On a dark-field signal generated by micrometer-sized calcifications in phase-contrast mammography,” Phys. Med. Biol. 58(8), 2713–2732 (2013). [CrossRef]  

8. S. Grandl, M. Willner, J. Herzen, D. Mayr, S. D. Auweter, A. Hipp, F. Pfeiffer, M. Reiser, and K. Hellerhoff, “Evaluation of phase-contrast CT of breast tissue at conventional X-ray sources - presentation of selected findings,” Z. Med. Phys. 23(3), 212–221 (2013). [CrossRef]  

9. S. Grandl, M. Willner, J. Herzen, A. Sztrókay-Gaul, D. Mayr, S. D. Auweter, A. Hipp, L. Birnbacher, M. Marschner, M. Chabior, M. Reiser, F. Pfeiffer, F. Bamberg, and K. Hellerhoff, “Visualizing typical features of breast fibroadenomas using phase-contrast CT: an ex-vivo study,” PLoS One 9(5), e97101 (2014). [CrossRef]  

10. C. Kottler, F. Pfeiffer, O. Bunk, C. Grünzweig, and C. David, “Grating interferometer based scanning setup for hard X-ray phase contrast imaging,” Rev. Sci. Instrum. 78(4), 043710 (2007). [CrossRef]  

11. Z.-Y. Wang and B. Han, “Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms,” Opt. Lett. 29(14), 1671–1673 (2004). [CrossRef]  

12. K. Hashimoto, H. Takano, and A. Momose, “Improved reconstruction method for phase stepping data with stepping errors and dose fluctuations,” Opt. Express 28(11), 16363–16384 (2020). [CrossRef]  

13. P.-P. Zhu, K. Zhang, Z.-T. Wang, Y.-J. Liu, X.-S. Liu, Z.-Y. Wu, S. A. McDonald, F. Marone, and M. Stampanoni, “Low-dose, simple, and fast grating-based X-ray phase-contrast imaging,” Proc. Natl. Acad. Sci. 107(31), 13576–13581 (2010). [CrossRef]  

14. J. Fu, X.-H. Shi, W. Guo, and P. Peng, “Fast X-ray Differential Phase Contrast Imaging with One Exposure and without Movements,” Sci. Rep. 9(1), 1113 (2019). [CrossRef]  

15. A. Ritter, F. Bayer, J. Durst, K. Gödel, W. Haas, T. Michel, J. Rieger, T. Weber, L. Wucherer, and G. Anton, “Simultaneous maximum-likelihood reconstruction for x-ray grating based phase-contrast tomography avoiding intermediate phase retrieval,” arXiv preprint arXiv:1307.7912 (2013).

16. A. Ritter, G. Anton, and T. Weber, “Simultaneous Maximum-Likelihood Reconstruction of Absorption Coefficient, Refractive Index and Dark-Field Scattering Coefficient in X-Ray Talbot-Lau Tomography,” PLoS One 11(10), e0163016 (2016). [CrossRef]  

17. M. V. Teuffenbach, T. Koehler, A. Fehringer, M. Viermetz, B. Brendel, J. Herzen, R. Proksa, E. J. Rummeny, F. Pfeiffer, and P. B. Noël, “Grating-based phase-contrast and dark-field computed tomography: a single-shot method,” Sci. Rep. 7(1), 7476 (2017). [CrossRef]  

18. D. Hahn, P. Thibault, A. Fehringer, M. Bech, T. Koehler, F. Pfeiffer, and P. B. Noël, “Statistical iterative reconstruction algorithm for X-ray phase-contrast CT,” Sci. Rep. 5(1), 10452 (2015). [CrossRef]  

19. P. M. Silverman, C. J. Cooper, D. I. Weltman, and R. K. Zeman, “Helical CT: practical considerations and potential pitfalls,” Radiographics 15(1), 25–36 (1995). [CrossRef]  

20. J. Fu, M. Willner, L.-Y. Chen, R.-B. Tan, K. Achterhold, M. Bech, J. Herzen, D. Kunka, J. Mohr, and F. Pfeiffer, “Helical differential X-ray phase-contrast computed tomography,” Phys. Med. 30(3), 374–379 (2014). [CrossRef]  

21. M. Marschner, M. Willner, G. Potdevin, A. Fehringer, P. B. Noël, F. Pfeiffer, and J. Herzen, “Helical X-ray phase-contrast computed tomography without phase stepping,” Sci. Rep. 6(1), 23953 (2016). [CrossRef]  

22. M. Büchner, Towards the development of an X-Ray phase contrast breast CT scanner. ETH Zurich. (2019).

23. C. Arboleda, Z.-T. Wang, and M. Stampanoni, “Tilted-grating approach for scanning-mode X-ray phase contrast imaging,” Opt. Express 22(13), 15447–15458 (2014). [CrossRef]  

24. Z.-T. Wang, K.-J. Kang, Z.-F. Huang, and Z.-Q. Chen, “Quantitative grating-based x-ray dark-field computed tomography,” Appl. Phys. Lett. 95(9), 094105 (2009). [CrossRef]  

25. A. W. Lohmann and D. E. Silva, “An interferometer based on the Talbot effect,” Opt. Commun. 2(9), 413–415 (1971). [CrossRef]  

26. S. D. Jean, F. M. Michael, K.-S. Carolyn, and F. Keyvan, “Effect of pitch, collimation, and reconstruction interval on low-contrast detectability in spiral computed tomography,” Proc SPIE. 3695, 603–614 (1999). [CrossRef]  

27. V. Revol, C. Kottler, R. Kaufmann, I. Jerjen, T. Lüthi, F. Cardot, P. Niedermann, U. Straumann, U. Sennhauser, and C. Urban, “X-ray interferometer with bent gratings: Towards larger fields of view,” Nucl. Instrum. Methods Phys. Res. A: Accel. Spectrom. Detect. Assoc. Equip. 648, S302–S305 (2011). [CrossRef]  

28. V. Revol, “X-ray phase contrast imaging by grating interferometry X-ray Phase Contrast Imaging by Grating Interferometry Dissertation,” ETH Zurich (2011).

29. M. Bech, O. Bunk, T. Donath, R. Feidenhans’l, C. David, and F. Pfeiffer, “Quantitative x-ray dark-field computed tomography,” Phys. Med. Biol. 55(18), 5529–5539 (2010). [CrossRef]  

30. W. A. Kalender, “Technical foundations of spiral CT,” Semin. Ultrasound. CT. MR. 15(2), 81–89 (1994). [CrossRef]  

31. Z. Xue, L.-L. Zhang, and J.-X. Pan, “Fast real-time image reconstruction with helical cone-beam ART,” in 13th International Conference on Computer Application and System Modeling (Iccasm) (2010), pp. 423–426.

32. D. Strong and T. Chan, “Edge-preserving and scale-dependent properties of total variation regularization,” Inverse Problems 19(6), S165–S187 (2003). [CrossRef]  

33. S. J. LaRoque, E. Y. Sidky, and X.-C. Pan, “Accurate image reconstruction from few-view and limited-angle data in diffraction tomography,” J. Opt. Soc. Am. A 25(7), 1772–1782 (2008). [CrossRef]  

34. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers,” Found. Trends Mach. Learn. 3(1), 1–122 (2010). [CrossRef]  

35. R. H. Byrd, P.-H. Lu, J. Nocedal, and C.-Y. Zhu, “A Limited Memory Algorithm for Bound Constrained Optimization,” SIAM J. Sci. Comput. 16(5), 1190–1208 (1995). [CrossRef]  

36. A. Chambolle, “An Algorithm for Total Variation Minimization and Applications,” Journal of Mathematical Imaging and Vision 20(1), 89–97 (2004). [CrossRef]  

37. S. Van Gogh, Z. Wang, M. Rawlik, C. Etmann, S. Mukherjee, C. B. Schonlieb, F. Angst, A. Boss, and M. Stampanoni, “INSIDEnet: Interpretable NonexpanSIve Data-Efficient network for denoising in Grating Interferometry Breast CT,“ Med. Phys. (accepted).

38. S. Shim, N. Saltybaeva, N. Berger, M. Marcon, H. Alkadhi, and A. Boss, “Lesion Detectability and Radiation Dose in Spiral Breast CT With Photon-Counting Detector Technology: A Phantom Study,” Invest. Radiol. 55(8), 515–523 (2020). [CrossRef]  

39. L. M. Warren, A. Mackenzie, D. R. Dance, and K. C. Young, “Comparison of the x-ray attenuation properties of breast calcifications, aluminium, hydroxyapatite and calcium oxalate,” Phys. Med. Biol. 58(7), N103–N113 (2013). [CrossRef]  

40. D. R. White, R. V. Griffith, and I. J. Wilson, “ICRU Report 46, Photon, Electron, Proton and Neutron Interaction Data for Body Tissues,” International Commission on Radiation Units and Measurements (1992).

41. K. M. Kempski, M. T. Graham, M. R. Gubbi, T. Palmer, and M. A. Lediju Bell, “Application of the generalized contrast-to-noise ratio to assess photoacoustic image quality,” Biomed. Opt. Express 11(7), 3684–3698 (2020). [CrossRef]  

42. M. Welvaert and Y. Rosseel, “On the Definition of Signal-To-Noise Ratio and Contrast-To-Noise Ratio for fMRI Data,” PLoS ONE 8(11), e77089 (2013). [CrossRef]  

43. W. van Aarle, W. Jan Palenstijn, J. De Beenhouwer, T. Altantzis, S. Bals, K. J. Batenburg, and J. Sijbers, “The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography,” Ultramicroscopy 157, 35–47 (2015). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (8)
Fig. 1.
Fig. 1. Schematics of helical GI-BCT. (a) The setup of GI-BCT consists of X-ray source, a source grating G0, a beam-splitter grating G1 and an analyzer grating G2. (b) The geometry of a helical scan. The imaging part moves along the $z - $ direction from bottom to top while it rotates around the sample. S is the X-ray source, R denotes the distance between S and the rotation axis ($z - $ axis), D represents the distance between the rotation axis and the detector, and h is the height of the upward movement of S after one rotation.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. One horizontal slice of simulated breast phantom for three contrast modalities, and the coronal slice and sagittal slice for attenuation and phase contrast along two blue dotted lines in horizontal slice. From left to right: the slice of attenuation, phase contrast, and scattering phantom. Attenuation and phase contrast both contain adipose, glandular, and skin tissue equivalent. Attenuation and scattering contain calcifications as well. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. Reconstructed results with simulated maps in noise-free case of the proposed method and PS method. The pitch is 1. From left to right: attenuation, ROIs of attenuation and phase contrast, phase contrast, and scattering. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. The real intensity, phase and visibility maps. The flat-field data had been measured with a CdTe, single-photon counting Dectris detector. The grating design energy was 46 keV, the tube was run at 70 kVp and 10 mA with an exposure time of 15 s.

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. Reconstructed results with different pitches in the noise-free case. From top to bottom, the pitches are 1, 1.4, and 2, respectively. From left to right: reconstructed attenuation, ROIs of attenuation and phase contrast, phase contrast, and scattering. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

View in Article | Download Full Size | PDF

Fig. 6.
Fig. 6. Reconstructed results with different noise levels. The pitch is 1. (a) From top to bottom, the intensity data are $\boldsymbol{I}_{({8000} )}^{measure}$, $\boldsymbol{I}_{({4000} )}^{measure}$ and $\boldsymbol{I}_{({2000} )}^{measure}$, respectively. From left to right: attenuation, ROIs of attenuation and phase contrast, phase contrast, and scattering. (b) Sagittal and coronal slice of the reconstructed attenuation and phase contrast in noisy case with intensity $\boldsymbol{I}_{({4000} )}^{measure}$ along two blue dotted lines in horizontal slice shown in Fig. 2. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

View in Article | Download Full Size | PDF

Fig. 7.
Fig. 7. Energy curves of the proposed method with real initial maps for the noise-free and noisy cases. (a) Noise-free case with pitches 1, 1.4 and 2. (b) Noisy case with pitch 1.

View in Article | Download Full Size | PDF

Fig. 8.
Fig. 8. Reconstructed results with tilted G2 and non-tilted G2 in the noise-free case (pitch is 1). From left to right: attenuation, phase contrast, and scattering. The gray value units for attenuation, phase contrast and scattering are attenuation coefficient $\mu [{c{m^{ - 1}}} ]$, phase shift coefficient $\varphi [{c{m^{ - 1}}} ]$ and $\varepsilon [{c{m^{ - 1}}} ]$, respectively.

View in Article | Download Full Size | PDF

Tables (4)

Tables Icon

Table 1. Experimental parameters.

View Table | View all tables in this article
Tables Icon

Table 2. MSEs, SSIMs, contrasts, contrasts’ relative errors of the images shown in Fig. 3 (simulated initial maps in noise-free case with pitch 1), and corresponding standard deviations computed across 20 slices (in the parentheses).

View Table | View all tables in this article
Tables Icon

Table 3. MSEs, SSIMs, contrasts, contrasts’ relative errors of the images shown in Fig. 5 (real initial maps in noise-free case with pitches 1, 1.4 and 2), and corresponding standard deviations computed across 20 slices (in the parentheses).

View Table | View all tables in this article
Tables Icon

Table 4. MSEs, SSIMs, contrasts, contrasts’ relative errors, SNRs, CNRs of the images shown in Fig. 6 (real initial maps in noisy case with pitch 1), and corresponding standard deviations computed across 20 slices (in the parentheses).

View Table | View all tables in this article

Equations (27)

Equations on this page are rendered with MathJax. Learn more.

{ x ( t ) = R cos ⁡ θ + t ( − D cos ⁡ θ − u sin ⁡ θ − R cos ⁡ θ ) y ( t ) = R sin ⁡ θ + t ( − D sin ⁡ θ + u cos ⁡ θ − R sin ⁡ θ ) z ( t ) = θ r h / ( 2 π ) + t v ,
I i = I i 0 T i [ 1 + V i 0 D i cos ⁡ ( ϕ i 0 − ϕ i ) ] .
T i = exp ( − ∫ 0 1 μ ( s ) d t ) , s ∈ l i ,
D i = exp ( − ∫ 0 1 ε ( s ) d t ) , s ∈ l i ,
ϕ i = c ϕ ∂ x ∫ 0 1 φ ( s ) d t , s ∈ l i .
T i = e x p ( − ∑ j = 1 N A i , j μ j )
D i = e x p ( − ∑ j = 1 N A i , j ε j ) ,
ϕ i = { c ϕ ( ∑ j = 1 N A i + 1 , j φ j − ∑ j = 1 N A i , j φ j ) / d k U ≤ i < ( k + 1 ) × U 0 i = ( k + 1 ) × U ,
arg ⁡ min x = ( μ , ε , φ ) { ∑ i ( I i m e a s u r e − I i ) 2 } ,
arg ⁡ min x = ( μ , ε , φ ) { ∑ i ( I i m e a s u r e − I i ) 2 + β μ | | μ | | T V + β ε | | ε | | T V + β φ | | φ | | T V } ,
x ( k + 1 / 2 ) = arg ⁡ min x = ( μ , ε , φ ) { ∑ i ( I i m e a s u r e − I i ) 2 } ,
μ ( k + 1 ) = arg ⁡ min μ { | | μ − μ ( k + 1 / 2 ) | | 2 2 + β μ | | μ | | T V } .
ε ( k + 1 ) = arg ⁡ min ε { | | ε − ε ( k + 1 / 2 ) | | 2 2 + β ε | | ε | | T V } .
φ ( k + 1 ) = arg ⁡ min φ { | | φ − φ ( k + 1 / 2 ) | | 2 2 + β φ | | φ | | T V } .
x ( k + 1 ) = ( μ ( k + 1 ) , ε ( k + 1 ) , φ ( k + 1 ) ) .
∂ f ∂ x = ( ∂ f ∂ μ , ∂ f ∂ ε , ∂ f ∂ φ ) .
∂ f ∂ φ = ∑ i ∂ f i ∂ φ = ∑ i ∂ f i ∂ ϕ i ⋅ ∂ ϕ i ∂ φ = ( [ ( ∂ ϕ 1 ∂ φ ) T , ( ∂ ϕ 2 ∂ φ ) T , … , ( ∂ ϕ M − 1 ∂ φ ) T ] ⋅ [ ∂ f 1 ∂ ϕ 1 , ∂ f 2 ∂ ϕ 2 , … , ∂ f M − 1 ∂ ϕ M − 1 ] T ) T
( ∂ f ∂ φ ) T = c ϕ d W T P ,
W 1 T = [ A 2 , 1 A 3 , 1 ⋅ ⋅ ⋅ A M , 1 A 2 , 2 A 3 , 2 ⋅ ⋅ ⋅ A M , 2 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ A 2 , N A 3 , N ⋅ ⋅ ⋅ A M , N ] and W 2 T = [ A 1 , 1 A 2 , 1 ⋅ ⋅ ⋅ A M − 1 , 1 A 1 , 2 A 2 , 2 ⋅ ⋅ ⋅ A M − 1 , 2 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ A 1 , N A 2 , N ⋅ ⋅ ⋅ A M − 1 , N ] ,
( ∂ f ∂ φ ) T = c ϕ d ( W 1 T P − W 2 T P ) .
( ∂ f ∂ φ ) T = c ϕ d ( W ¯ 1 T P ¯ − W ¯ 2 T P ¯ ) .
Contrast = 20 log 10 ( m ROI 2 m ROI 1 ) ,
CNR = | m ROI 2 − m ROI 1 | σ ROI 2 2 + σ ROI 1 2 ,
SNR = m ROI 2 σ ROI 1 ,
MSE = 1 R × S × T ∑ r = 1 R ∑ s = 1 S ∑ t = 1 T ( g r , s , t − g ¯ r , s , t ) 2 ,
SSIM ( g , g ¯ ) = ( 2 m g m g ¯ + C 1 ) ( 2 σ g g ¯ + C 2 ) ( m g 2 + m g ¯ 2 + C 1 ) ( σ g 2 + σ g ¯ 2 + C 2 ) ,
f = ∑ i ( I i m e a s u r e − I i 0 T i [ 1 + V i 0 D i cos ⁡ ( ϕ i 0 − ϕ i ) ] ) 2 .
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use