Fast and robust reconstruction algorithm for fluorescence diffuse optical tomography assuming a cuboid target

TL;DR

This paper introduces a rapid and reliable algorithm for fluorescence diffuse optical tomography that accurately locates cuboid-shaped targets, validated through numerical and in vivo experiments.

Contribution

The paper presents a novel fast and robust reconstruction algorithm specifically designed for fluorescence diffuse optical tomography assuming cuboid targets.

Findings

Algorithm successfully locates cuboid targets in simulations.

Validated with in vivo experimental results.

Demonstrates improved speed and robustness over existing methods.

Abstract

A fast and robust algorithm for fluorescence diffuse optical tomography is proposed. We identify the location of a fluorescence target by assuming a cuboid. The proposed numerical method is verified by a numerical experiment and an in vivo experiment.

Equations

$$\left(\frac{1}{c}\frac{\partial }{\partial t}-D\Delta +{\mu }_a\right)u_e=0,\left(\frac{1}{c}\frac{\partial }{\partial t}-D\Delta +{\mu }_a\right)u_m=F,$$

$$F(\mathrm{r},t)=\frac{n(\mathrm{r})}{\tau }{\int }_0^t{e}^{-(t-s)/\tau }{u}_e(\mathrm{r},s)ds.$$

$$-\frac{\partial }{\partial z}{u}_e+\beta u_e=f(t)\delta (\rho -{\rho }_s),-\frac{\partial }{\partial z}{u}_m+\beta u_m=0.$$

$$U_e({\mathrm{r}}_d,t;{\mathrm{r}}_s)={\int }_0^{t}R(t-s)u_e({\mathrm{r}}_d,s)ds,U_m({\mathrm{r}}_d,t;{\mathrm{r}}_s)={\int }_0^{t}R(t-s)u_m({\mathrm{r}}_d,s)ds,$$

$$q(t)={\int }_0^{t}R(t-s)f(s)ds.$$

$$n(\mathrm{r})=\{\begin{array}{rl}M,& \mathrm{r}\in {\Omega }_c,\\ 0,& \mathrm{r}\notin {\Omega }_c,\end{array}$$

$$U_m({\mathrm{r}}_d,t;{\mathrm{r}}_s)={\int }_0^{t}Q(s){\int }_0^{t-s}{\int }_{{\Omega }_0}n({\mathrm{r}}^{\prime })G({\mathrm{r}}_d,{\mathrm{r}}^{\prime };t-s-s^{\prime })G({\mathrm{r}}^{\prime },{\mathrm{r}}_s;s^{\prime })d\mathrm{r}d{s}^{\prime }ds.$$

$$U_m({\mathrm{r}}_d,t;{\mathrm{r}}_s)=M{\int }_0^{t}Q(s){\int }_0^{t-s}{f}_1({\rho }_d,{\rho }_s,t-s,t^{\prime };x_1,x_2,y_1,y_2)f_2(t-s,t^{\prime };z_1,z_2)d{t}^{\prime }ds,$$

$$f_1({\rho }_d,{\rho }_s,t,s;x_1,x_2,y_1,y_2)=\frac{e^{-{\mu }_{a}ct}}{4^3{\pi }^2{D}^{2}t\sqrt{(t-s)s}}{e}^{-\frac{(x_d-x_s{)}^2+(y_d-y_s{)}^2}{4Dct}}$$

$$\times \left[h(x_d,x_s,t,s;x_2)-h(x_d,x_s,t,s;x_1)\right]\left[h(y_d,y_s,t,x;y_2)-h(y_d,y_s,t,s;y_1)\right],$$

$$F(\mathrm{a})=\sqrt{\sum _{i=1}^{N_{\mathrm{SD}}}\sum _{k=1}^{N_t}{\left(U_m({\mathrm{r}}_d^i,t_k^i;{\mathrm{r}}_s^i)-U_m^{\exp }({\mathrm{r}}_d^i,t_k^i;{\mathrm{r}}_s^i)\right)}^2}.$$

$${\mathrm{a}}_{\ast }^{\mathrm{cubic}}=(x_0^{\prime },y_0^{\prime },z_0^{\prime },l^{\prime },M^{\prime })=(0.0,0.0,11.24,4.089,0.0086).$$

$${\mathrm{a}}_{\ast }^{\mathrm{cuboid}}=(x_1^{\ast },x_2^{\ast },y_1^{\ast },y_2^{\ast },z_1^{\ast },z_2^{\ast },M^{\ast })=(-1.074,1.036,-2.146,2.166,9.908,12.02,0.029).$$

$$(x_0^{\prime },y_0^{\prime },z_0^{\prime })=(-4.12,7.72,17,25),l^{\prime }=\mathrm{3.92.}$$

$$x_1^{\ast }=-5.16,x_2^{\ast }=-3.11,y_1^{\ast }=3.83,y_2^{\ast }=12.03,z_1^{\ast }=16.05,z_2^{\ast }=\mathrm{16.34.}$$

Full text

††journal: boe††articletype: Research Article

1 Introduction

Fluorescence diffuse optical tomography (FDOT) is one type of optical tomography which makes use of fluorescence light from fluorophore. In FDOT, diffuse light from fluorophore such as Indocyanine Green (ICG) is detected on the boundary of biological tissue to obtain tomographic images [8, 9]. FDOT has been verified in vivo [4, 7] and also in clinical research for breast cancer [1].

To obtain reconstructed images with good quality, a large number of source-detector pairs is necessary. The number of source-detector pairs can be increased by time-dependent experiments [3]. The superiority of the time-resolved approach over the continuous-wave approach in FDOT was concluded [2] although the quality of images can be improved even for a time-independent experiment if a large data set is used [10].

Even if many sources and detectors are used, the resolution of tomographic images from optical tomography cannot be compared to that from X-ray CT. In this paper, we assume a cuboid or rectangular parallelepiped in the medium. Then we try to find the cuboid instead of giving up reconstructing the shape of the target. Thus we can develop a fast and robust numerical algorithm using a cuboid.

The paper is organized as follows. In Sec. 2, we develop the formulation of our FDOT and give an analytical formula for the emission light. The proposed numerical scheme is described in Sec. 3 with a numerical example of an ellipsoidal target. In Sec. 4, we validate our numerical method with a beef experiment. Finally, Sec. 5 is devoted to conclusion.

2 Formulation

Let us suppose that a fluorescence target is embedded in biological tissue occupying the half space ($-\infty<x<\infty$, $-\infty<y<\infty$, $0<z<\infty$). Let $c$ be the speed of light in the medium. Let $u_{e}(\boldsymbol{\mathrm{r}},t)$, $u_{m}(\boldsymbol{\mathrm{r}},t)$ be the energy densities of the excitation light and emission light, respectively. Here, $\boldsymbol{\mathrm{r}}=(\boldsymbol{\mathrm{\rho}},z)$ with $\boldsymbol{\mathrm{\rho}}=(x,y)$ is position and $t$ is time. We assume that the reduced scattering coefficient $\mu_{s}^{\prime}$ and the absorption coefficient $\mu_{a}$ are constants everywhere in the medium. In the medium ($z>0$), $u_{e}(\boldsymbol{\mathrm{r}},t)$ and $u_{m}(\boldsymbol{\mathrm{r}},t)$ obey the following diffusion equations.

[TABLE]

where $D=1/(3\mu_{s}^{\prime})$ and

[TABLE]

Here, $n(\boldsymbol{\mathrm{r}})$ is proportional to the fluorophore concentration and $\tau$ is the fluorescence lifetime. We assumed that $u_{e}=0$ and $u_{m}=0$ at $t=0$ but the sample is illuminated at position $\boldsymbol{\mathrm{r}}_{s}=(\boldsymbol{\mathrm{\rho}}_{s},0)$ by a pencil beam of the temporal profile $f(t)$ ($t>0$) in the $x$-$y$ plane. Thus, $u_{e},u_{m}$ satisfy the Robin boundary conditions at $z=0$ as

[TABLE]

The parameter $\beta$ is given by $\beta=\frac{1}{2D}(1-2\int_{0}^{1}R(\mu)\mu\,d\mu)/(1+\int_{0}^{1}R(\mu)\mu^{2}\,d\mu)$ with the Fresnel reflectance $R(\mu)$, which depends on the refractive index of the medium. Suppose the out-going light is detected at $\boldsymbol{\mathrm{r}}_{d}$ on the boundary. The excitation light and emission light are detected through a response function $R$ as

[TABLE]

where the function $R$ is determined by the detector. These $U_{e},U_{m}$ correspond to experimentally measured light $U_{e}^{\rm exp}(\boldsymbol{\mathrm{r}}_{d},t;\boldsymbol{\mathrm{r}}_{s})$, $U_{m}^{\rm exp}(\boldsymbol{\mathrm{r}}_{d},t;\boldsymbol{\mathrm{r}}_{s})$.

Let $G(\boldsymbol{\mathrm{r}},\boldsymbol{\mathrm{r}}^{\prime};t)$ be the Green’s function which satisfies $(\frac{1}{c}\frac{\partial}{\partial t}-D\Delta+\mu_{a})G(\boldsymbol{\mathrm{r}},\boldsymbol{\mathrm{r}}^{\prime};t)=\delta(\boldsymbol{\mathrm{r}}-\boldsymbol{\mathrm{r}}^{\prime})\delta(t)$, with the Robin boundary condition $-\frac{\partial}{\partial z}G+\beta G=0$. The Green’s function is obtained as $G(\boldsymbol{\mathrm{r}},\boldsymbol{\mathrm{r}}^{\prime};t)=c(4\pi Dct)^{-3/2}e^{-\mu_{a}ct}\exp(-\frac{(x-x^{\prime})^{2}+(y-y^{\prime})^{2}}{4Dct})g(z,z^{\prime};t)$, where $g(z,z^{\prime};t)=\exp(-\frac{(z+z^{\prime})^{2}}{4Dct})+\exp(-\frac{(z-z^{\prime})^{2}}{4Dct})-2\beta\sqrt{\pi Dct}\exp\left(\beta(z+z^{\prime})+\beta^{2}Dct\right)\mathop{\mathrm{erfc}}(\frac{z+z^{\prime}+2\beta Dct}{\sqrt{4Dct}})$. We introduce

[TABLE]

This $q(t)$ is called the instrument response function.

Let $\Omega_{c}$ be a cuboid specified by $x_{1},x_{2},y_{1},y_{2},z_{1},z_{2}$ as $x_{1}<x<x_{2}$, $y_{1}<y<y_{2}$, $z_{1}<z<z_{2}$. We assume that

[TABLE]

where $M>0$ is a constant. Then we obtain

[TABLE]

where we introduced $Q(t)=\frac{D}{\tau}\int_{0}^{t}e^{-t^{\prime}/\tau}q(t-t^{\prime})\,dt^{\prime}$. Equation (7) can be rewritten as

[TABLE]

where

[TABLE]

with $h(x_{d},x_{s},t,s;x)=\mathop{\mathrm{erf}}\left(\sqrt{\frac{t}{4Dc(t-s)s}}\left(x-\frac{sx_{d}+(t-s)x_{s}}{t}\right)\right)$ and $f_{2}(t,s;z_{1},z_{2})=\int_{z_{1}}^{z_{2}}g(0,z^{\prime};t-s)g(z^{\prime},0;s)\,dz^{\prime}$. We can compute $U_{m}(\boldsymbol{\mathrm{r}}_{d},t;\boldsymbol{\mathrm{r}}_{s})$ using (8).

3 Identification of a cuboid

To illustrate our proposed inversion scheme, we will reconstruct $n(\boldsymbol{\mathrm{r}})$ using solely $U_{m}^{\rm exp}$. The reconstruction will be done with the Levenberg-Marquardt method [5, 6]. The choice of initial guesses is important.

Let $N_{\rm SD}$ be the number of source-detector pairs. The positions of each source and detector are denoted by $\boldsymbol{\mathrm{r}}_{s}^{i},\boldsymbol{\mathrm{r}}_{d}^{i}$ ($i=1,\dots,N_{\rm SD}$). At each detector, light is measured at $N_{t}$ temporal points.

We first look for a rough location of the target in the $x$-$y$ plane or a region of interest $\Gamma$ on the boundary, under which the target is likely to be embedded, by observing $I_{i}=\int_{0}^{T}U_{m}^{\rm exp}(\boldsymbol{\mathrm{r}}_{d}^{i},t;\boldsymbol{\mathrm{r}}_{s}^{i})\,dt$ ($i=1,\dots,N_{\rm SD}$). We call this step of setting $\Gamma$ the topography process. Then we assume a cubic target whose location and size are determined by four parameters $x_{0},y_{0},z_{0},l$ such that $x_{1}=x_{0}-l/2$, $x_{2}=x_{0}+l/2$, $y_{1}=y_{0}-l/2$, $y_{2}=y_{0}+l/2$, $z_{1}=z_{0}-l/2$, $z_{2}=z_{0}+l/2$. Choosing the initial guess for $x_{0},y_{0}$ inside $\Gamma$, we solve the inverse problem of determining $x_{0},y_{0},z_{0},l,M$ by the Levenberg-Marquardt method and obtain $\boldsymbol{\mathrm{a}}_{*}^{\rm cubic}=(x_{0}^{\prime},y_{0}^{\prime},z_{0}^{\prime},l^{\prime},M^{\prime})$. We refer to this step as the cubic tomography. Finally with the obtained values $\boldsymbol{\mathrm{a}}_{*}^{\rm cubic}$ as the initial guess $\boldsymbol{\mathrm{a}}_{0}^{\rm cuboid}$, we solve the inverse problem of determining a cuboid target. In this way, we can identify $\Omega_{c}$ for $n(\boldsymbol{\mathrm{r}})$ by obtaining reconstructed values $\boldsymbol{\mathrm{a}}_{*}^{\rm cuboid}=(x_{1}^{*},x_{2}^{*},y_{1}^{*},y_{2}^{*},z_{1}^{*},z_{2}^{*},M^{*})$. This last step is called the cuboid tomography. The algorithm is summarized as follows.

Step 1.

(Topography process) Find $\Gamma$ on the boundary.

Step 2.

(Cubic tomography) By searching underneath $\Gamma$, obtain reconstructed values $\boldsymbol{\mathrm{a}}_{*}^{\rm cubic}$.

Step 3.

(Cuboid tomography) Find $\boldsymbol{\mathrm{a}}_{*}^{\rm cuboid}$ starting with $\boldsymbol{\mathrm{a}}_{*}^{\rm cubic}$.

Below, we will demonstrate our algorithm with a numerical experiment. We set $\tau=0$, $f(t)=R(t)=\delta(t)$. Let us assume an ellipsoidal target $\{x^{2}/1.5^{2}+y^{2}/3^{2}+(z-11)^{2}/1.5^{2}\leq 1\}$, where the unit of length is ${\rm mm}$, with $n(\boldsymbol{\mathrm{r}})=0.02\,{\rm mm}^{-1}$. In the numerical experiment, $5\%$-noise is added and $U_{m}^{\rm exp}=U_{m}(1+0.05\varepsilon)$, where $\varepsilon$ is drawn from the standard Gaussian distribution. We use the following $N_{\rm SD}=32$ source-detector pairs $p_{i}=(\boldsymbol{\mathrm{r}}_{s}^{i},\boldsymbol{\mathrm{r}}_{d}^{i})$ ($i=1,\dots,N_{\rm SD}$) on the boundary: $p_{4j-3}=(\xi_{j},\zeta_{j}+10\sqrt{3};\xi_{j}-10,\zeta_{j})$, $p_{4j-2}=(\xi_{j},\zeta_{j}+10\sqrt{3};\xi_{j}+10,\zeta_{j})$, $p_{4j-1}=(\xi_{j},\zeta_{j}-10\sqrt{3};\xi_{j}-10,\zeta_{j})$, $p_{4j}=(\xi_{j},\zeta_{j}-10\sqrt{3};\xi_{j}+10,\zeta_{j})$, $j=1,\dots,8$, where $(\xi_{1},\zeta_{1})=(-10,10)$, $(\xi_{2},\zeta_{2})=(-10,0)$, $(\xi_{3},\zeta_{3})=(-10,-10)$, $(\xi_{4},\zeta_{4})=(0,-10)$, $(\xi_{5},\zeta_{5})=(10,-10)$, $(\xi_{6},\zeta_{6})=(10,0)$, $(\xi_{7},\zeta_{7})=(10,10)$, $(\xi_{8},\zeta_{8})=(0,10)$. We set $N_{t}=20$ and $t_{k}^{i}=t_{0}^{i}+(k-1)\Delta t$ ($k=1,\dots,N_{t}$), where $\Delta t=6.67\,{\rm ps}$. Here, $t_{0}^{i}=(k_{0}^{i}-10)\Delta t$, where $U_{m}^{\rm exp}(\boldsymbol{\mathrm{r}}_{d}^{i},t;\boldsymbol{\mathrm{r}}_{s}^{i})$ has the peak at $(k_{0}^{i}-1)\Delta t$ ($i=1,\dots,N_{\rm SD}$). Thus we have $N_{\rm SD}N_{t}$ ($=640$) measured values $U_{m}^{\rm exp}(\boldsymbol{\mathrm{r}}_{d}^{i},t_{k}^{i};\boldsymbol{\mathrm{r}}_{s}^{i})$ ($i=1,\dots,N_{\rm SD}$, $k=1,\dots,N_{t}$). Let $\boldsymbol{\mathrm{a}}$ be a vector which contains unknown parameters to be reconstructed. We define

[TABLE]

By using the Levenberg-Marquardt method, we try to find $\boldsymbol{\mathrm{a}}=\boldsymbol{\mathrm{a}}_{*}$ that minimizes $F(\boldsymbol{\mathrm{a}})$.

(Step 1) We find $I_{4}=I_{10}=I_{17}=I_{27}=5.9\times 10^{-7}$, whereas $I_{1}=I_{11}=I_{20}=I_{26}=1.0\times 10^{-9}$. Hence we can set $\Gamma=\{-10<x_{0}<10,-10<y_{0}<10\}$. At this moment, we have no knowledge about the depth at which the target is embedded. We suppose $0<z_{0}<30$ together with $0<l<\min(20,2z_{0})$ and $0<M<10$.

(Step 2) We set the initial values of $x_{0},y_{0}$ in $\Gamma$. For example, we can start from $\boldsymbol{\mathrm{a}}_{0}^{\rm cubic}=(x_{0},y_{0},z_{0},l,M)=(2,2,5,4,0.1)$. The reconstructed values are

[TABLE]

We note that choosing $x_{0},y_{0}$ in $\Gamma$ is important for fast convergence. If we set $(x_{0},y_{0})=(-15,-15)$, more than $100$ iterations are needed whereas about $10$ iterations are sufficient for suitable $x_{0},y_{0}$.

(Step 3) Now we give $\boldsymbol{\mathrm{a}}_{0}^{\rm cuboid}$ from $\boldsymbol{\mathrm{a}}_{*}^{\rm cubic}$ as $\boldsymbol{\mathrm{a}}_{0}^{\rm cuboid}=(x_{1},x_{2},y_{1},y_{2},z_{1},z_{2},M)=(x_{0}^{\prime}-l^{\prime}/2,x_{0}^{\prime}+l^{\prime}/2,y_{0}^{\prime}-l^{\prime}/2,y_{0}^{\prime}+l^{\prime}/2,z_{0}^{\prime}-l^{\prime}/2,z_{0}^{\prime}+l^{\prime}/2,M^{\prime})$. Then assuming a cuboid, we obtain reconstructed values $\boldsymbol{\mathrm{a}}_{*}^{\rm cuboid}$ as

[TABLE]

The reconstructed cuboid is shown in Fig. 1 with the ellipsoidal target. The cuboid (blue) and ellipsoid (red) are shown in two-dimensional planes in Fig. 2. We emphasize that the narrowing process $\Gamma\to\boldsymbol{\mathrm{a}}_{0}^{\rm cubic}\to\boldsymbol{\mathrm{a}}_{0}^{\rm cuboid}$ is essential. If $\boldsymbol{\mathrm{a}}_{0}^{\rm cubic}$ is used for the initial guess for the cuboid tomography, the iteration of the Levenberg-Marquardt method does not converge except for some special cases.

4 Beef experiment

Let us reconstruct a tube which contains ICG in beef with our reconstruction scheme. Figure 3 shows how the tube was placed in the beef. The tube has the shape of a cylinder of length $8\,{\rm mm}$ and diameter $2\,{\rm mm}$. We performed time-dependent measurements using a holder placed on the top of the beef. Four optical fibers (two are for sources and the other two are for detectors) are attached to the holder. The table on which the beef is placed changes positions while the holder is fixed.

Optical parameters for the beef are $\mu_{s}^{\prime}=0.92\,{\rm mm}^{-1}$, $\mu_{a}=0.023\,{\rm mm}^{-1}$. The refractive index is set to $1.37$. Moreover, $\tau=0.6\,{\rm ns}$. After subtracting the background fluorescence, we could use $N_{\rm SD}=16$ source-detector pairs $p_{i}=(\boldsymbol{\mathrm{r}}_{s}^{i},\boldsymbol{\mathrm{r}}_{d}^{i})$ ($i=1,\dots,N_{\rm SD}$): $p_{1}=(-5-10\sqrt{3},0;-5,10)$, $p_{2}=(-5+10\sqrt{3},0;-5,10)$, $p_{3}=(-5-10\sqrt{3},5;-5,15)$, $p_{4}=(-5+10\sqrt{3},5;-5,15)$, $p_{5}=(-10\sqrt{3},0;0,10)$, $p_{6}=(-10\sqrt{3},5;0,15)$, $p_{7}=(5-10\sqrt{3},5;5,-5)$, $p_{8}=(5-10\sqrt{3},5;5,15)$, $p_{9}=(5-10\sqrt{3},0;5,10)$, $p_{10}=(5-10\sqrt{3},-5;5,5)$, $p_{11}=(-10\sqrt{3},-5;0,5)$, $p_{12}=(-10+10\sqrt{3},0;-10,10)$, $p_{13}=(-15+10\sqrt{3},0;-15,10)$, $p_{14}=(-15+10\sqrt{3},5;-15,-5)$, $p_{15}=(-15+10\sqrt{3},5;-15,15)$, $p_{16}=(-10+10\sqrt{3},5;-10,15)$, where the unit of coordinates is ${\rm mm}$. We set $N_{t}=20$ and $\Delta t=6.1\,{\rm ps}$. The peak of $U_{m}^{\rm exp}(\boldsymbol{\mathrm{r}}_{d}^{i},t;\boldsymbol{\mathrm{r}}_{s}^{i})$ comes at $t=(k_{0}^{i}-1)\Delta t$ ($i=1,\dots,N_{\rm SD}$). Measured times used for reconstruction are $t_{k}^{i}=t_{0}^{i}+(k-1)\Delta t$ ($k=1,\dots,N_{t}$), where $t_{0}^{i}=(k_{0}^{i}-11)\Delta t$.

First we set $\Gamma=\{-10<x<0,5<y<20\}$ since measured values were large in this region such as $I_{5}=3.3\times 10^{4}$, $I_{6}=4.0\times 10^{4}$, $I_{12}=2.9\times 10^{4}$, $I_{16}=3.4\times 10^{4}$. Although the value of $n(\boldsymbol{\mathrm{r}})$, i.e., $M$ for the cube and cuboid, is a parameter to be reconstructed, $M$ is determined only up to a constant which comes from the property of the detector. We set $(x_{0},y_{0},z_{0})=(-5,10,7)$, $l=2$, and obtain

[TABLE]

With the above values as the initial guess we perform the iterative method once again. The reconstructed values are obtained as

[TABLE]

The position of the cuboid is what we expected. For example, we see that the orientation of the cylinder should be almost parallel to the $y$-axis. Thus the position of the fluorescence tube is identified with the proposed numerical scheme.

5 Conclusion

By giving up shape reconstruction, we can identify the location of the target by reconstructing only several unknown parameters. Even for these several parameters, the Levenberg-Marquardt algorithm is not stable unless good initial guesses are used. Thus the proposed procedure of narrowing target domains as $\Gamma\to\boldsymbol{\mathrm{a}}_{0}^{\rm cubic}\to\boldsymbol{\mathrm{a}}_{0}^{\rm cuboid}$ is important for the iterative method to work.

Although we assumed the half space in this paper and made use of an analytical solution to the diffusion equation, the proposed algorithm works also in more general cases where diffusion equations must be solved numerically by finite difference method or finite element method [11]. The proposed algorithm can be applied not only to the Levenberg-Marquardt method but also to other iterative schemes such as the conjugate gradient method and the Gauss-Newton method.

In this paper, the algorithm was explained using a single target. The generalization of the method for multiple targets is straightforward at least if the number of targets is known.

Funding

NSFC (No.11421110002, No.11531005, No.91730304) (to J. J. Liu); Grant-in-Aid for Scientific Research (17K05572 and 17H02081) of the Japan Society for the Promotion of Science (JSPS) (to M. Machida); Grant-in-Aid for Scientific Research (15K21766 and 15H05740) of JSPS (to G. Nakamura). Also, the JSPS A3 foresight program: Modeling and Computation of Applied Inverse Problems.

Disclosures

The authors declare that there are no conflicts of interest related to this article.