A Kalman-based tomographic scheme for directly reconstructing activation levels of brain function

Bingyuan Wang; Tiantian Pan; Yao Zhang; Dongyuan Liu; Jingying Jiang; Huijuan Zhao; Feng Gao

doi:10.1364/OE.27.003229

1. Introduction

During the past decades, fNIRS has been playing an increasingly important role in researching human brains, because of its non-invasiveness, potential wearability, higher temporal resolution and reasonable spatial resolution [1–3]. It is based on the physiological phenomenon called neuro-vascular-coupling that the vascular will expand and the blood volume will increase when the nearby neuro cells are activated by the brain stimulation [4,5]. More specifically, the oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR) concentrations will increase and decrease respectively at the onset of the stimulation, and return back to the initial values after the stimulation. Since the HbO and HbR are the primary chromophores in the near-infrared window (650 nm~1000 nm), light at multiple near-infrared wavelengths can penetrate into the human brain tissue deep enough to detect the temporal and spatial distribution of the changes in the HbO and HbR concentrations in the cerebral cortex region under investigation, from which the neural activations are inferred by a temporal-space analysis [6]. In the fNIRS process, it is normally assumed that the absorption coefficient ( $μ_{a}$ ) changes by the changes in the HbO and HbR concentrations but the reduced scattering coefficient ( ${μ^{'}}_{s}$ ) almost remains unchanged [2]. Nevertheless, both the optical properties significantly influence the light propagation inside the brain tissues. This finally leads to the changes in the light intensity detected on the scalp surface with the changes in the HbO and HbR concentrations in a profound way, making their accurate assessment a rather nontrivial task [4,7].

To quantitatively relate the fNIRS measurement to the changes in the HbO and HbR concentrations, two image reconstruction schemes, referred to as optical topography (OT) and diffuse optical tomography (DOT), are successively proposed based on the modified-Lambert-Beer-law (MLBL) and the more sophisticated light propagation model, respectively [4,8]. OT scheme is commonly used in the commercial fNIRS systems, where the changes in $μ_{a}$ are separately recovered for each non-overlapping channel (source-detector pair), and finally transformed into a 2-dimensional image covering the cortex surface under probing. The spatial resolution of this scheme is inherently limited by the source-detector separation (SDS), and the quantitativeness depends on the rationality of the MLBL as well as the accuracy of the differential path factor (DPF). DOT schemes are proposed to tomographically reconstruct 3-dimensional spatial distribution of the changes in $μ_{a}$ , and expected to provide improved spatial resolution and quantitativeness [9]. However, the inverse problems involved in DOT are difficult to solve because of the ill-posedness caused by the highly scattering nature of the near-infrared light in the tissues and the much less number of measurements than that of the unknowns [10,11]. To mitigate the adversities, a variety of a priori information has been adopted in the reconstruction process by use of regularization strategies [12].

For fNIRS, the temporal resolution of a measurement is one of the crucial performances. The best way of implementing a high temporal resolution of the measurement is to employ a fully parallel mode where the encoded sources are lit simultaneously and the detectors decode the source signals from the received signals in parallel. Although the fully parallel modes have been widely applied in OT scheme where the non-overlapping channels with the same SDS are deployed, it is normally unsuitable for DOT scheme where an overlapping arrangement of the channels with different SDSs is necessary and the inter-channel crosstalk caused by the channels with shorter SDSs can dominate the signals at the channels with longer SDSs [13,14]. To cope with this problem, a multi-field illumination mode is used where, at each field, part of the sources is lit simultaneously and their associated channels with only the same SDSs are working simultaneously. The data collected at all the fields accomplishes a complete frame of the measurements for all the required channels. The frame is repeated one after another to temporally cover the brain activation process in correspondence to the stimulation.

In the conventional indirect reconstruction approaches, the spatial distribution of the changes in $μ_{a}$ are assumed to be temporally-invariant during a frame and reconstructed by DOT for every frame. For all the sites, the time-courses of the changes in $μ_{a}$ with a frame as the temporal unit are analyzed independently with the help of the temporal a priori information to obtain the spatial distribution of the activation levels [14]. But this temporally-invariant assumption neglects the changes in $μ_{a}$ over fields in the same frame and the temporal resolution is the duration needed by one frame although the measurements are obtained field by field. Therefore, this temporally-invariant assumption is inaccurate and may degrade the calculated spatial distribution of the activation levels.

To overcome the above problems, in this paper we propose a Kalman-based tomographic scheme to directly reconstruct the spatial distribution of the activation levels of the brain function without first reconstructing the spatial distribution of the changes in $μ_{a}$ for all the frames. The estimation of the spatial distribution of the activation levels is updated field by field instead of frame by frame and the temporal resolution is also decreased from the duration needed by a frame to that of a field. In addition, this direct approach avoids the error propagation between the two successive steps in the indirect approaches. Furthermore, it makes better use of the temporal a priori information and avoids solving straightforwardly the ill-posed inverse problems in DOT. To validate the potential performance improvements of this proposed scheme, some dynamic simulative and phantom experiments were conducted, which showed that the proposed scheme provided higher quality images.

2. Methods

The indirect and direct Kalman estimation approaches based on the hemodynamic response function (HRF) model are described.

2.1 HRF model

In the functional magnetic resonance imaging (fMRI), researchers assume the brain to be a linear time-invariant system and introduce the so-called ‘canonical’ HRF (cHRF) to express the time-course of the blood oxygen level dependence (BOLD) signal in response to impulse stimulation [15]. Accordingly, the time-course of the BOLD in response to any stimulation can be calculated as the convolution of the cHRF and the specific stimulations. Likewise, the researchers in the fNIRS also utilize the HRF to model the time-course of the concentration of the HbO and HbR during the brain activation process [16]. In this paper, we utilize the following cHRF whose specific parameters for the HbO are the same as those used by Viola Bonomini et al. [17]

c H R F (k) = α_{1} [Γ_{n} (k, τ_{1}, ρ_{1}) - α_{2} Γ_{n} (k, τ_{2}, ρ_{2})]

where

Γ_{n} (k, τ_{j}, ρ_{j}) = \frac{1}{p! τ_{j}} {(\frac{k - ρ_{j}}{τ_{j}})}^{p} e^{- \frac{k - ρ_{j}}{τ_{j}}} δ (k - ρ_{j})

with

δ (k - ρ_{j}) = {\begin{cases} 1, if k - ρ_{j} \geq 0 \\ 0, o t h e r w i s e \end{cases}

;

α_{1}

regulates the amplitude;

α_{2}

determinates the undershoot;

τ_{1}

and

τ_{2}

regulate the cHRF shape;

ρ_{1}

and

ρ_{2}

tune its scale; these specific parameters are set as follows:

p = 5

;

α_{1} = 1282

;

α_{2} = 0.17

;

τ_{1} = τ_{2} = 1

;

ρ_{1} = - 0.5

;

ρ_{2} = 3.5

. The cHRF of the HbR is inverted and with a maximum set at 1/3 with respect to that of the HbO. Since the cHRF of the HbR and HbO have the same form, the cHRF for

μ_{a}

is defined as the summation of the cHRF of the HbO and HbR multiplied by their own respective extinction coefficients [18]. The HRF of

μ_{a}

is calculated by convolving the cHRF of

μ_{a}

and a box function denoted by

t (k) = {\begin{cases} 1, k \in [10, 20] \\ 0, o t h e r w i s e \end{cases}

. Then the HRF of

μ_{a}

is normalized by its maximum value and expressed as

h (k)

.

2.2 Indirect Kalman estimation

The spatial distribution of the activation levels are separately recovered for every site based on its time-course of the changes in $μ_{a}$ , which is extracted from the images reconstructed independently by DOT for every frame.

2.2.1 DOT reconstruction for each frame

We define the boundary measurements at the s-th (s = 1, 2, …, S) field in the f-th (f = 1, 2, …, F) frame as $S (k)$ with $k = s + S (f - 1)$ . In the f-th frame, $S (k)$ of all the fields are assembled together to form $F (f) = {S {[S (f - 1) + 1]}^{T}, S {[S (f - 1) + 2]}^{T}, \dots, S {[S (f - 1) + S]}^{T}}^{T} \in R^{M \times 1}$ with M being the total number of the measurements in one frame. The boundary measurements of the brain at the rest state are denoted as $F (1)$ . Moreover, the changes in $μ_{a}$ in the f-th (f = 2, 3, …, F) frame are denoted as $x (f) \in R^{N \times 1}$ with $N$ representing the total number of sites. In the first step of the conventional indirect approaches, $x (f)$ (f = 2, 3, …, F) can be reconstructed independently one frame after another by solving the following approximate linear equation

y (f) = J x (f)

where

y (f) = F (f) - F (1)

;

J \in R^{M \times N}

is the sensitivity matrix of

y (f)

to

x (f)

and is computed based on the spatial distribution of the optical properties of the brain at the rest state.

2.2.2 Estimation of the spatial distribution of the activation levels

After implementing Eq. (2) for all the frames, the time-course of the changes in $μ_{a}$ for the n-th (n = 1,2,…,N) site inside the medium $X_{n} \in R^{1 \times (F - 1)}$ can be obtained. In terms of the general linear model (GLM), $X_{n}$ is interpreted as the linear combination of two explanatory variables

{(X_{n})}^{T} = H β_{n}

where

H

is a design matrix composed of

h (k)

and a constant vector;

{(\cdot)}^{T}

denotes the transposition computation;

β_{n} = {[β_{1, n}, β_{2, n}]}^{T}

is a vector with

β_{1, n}

and

β_{2, n}

denoting the coefficients of

h (k)

and the constant vector, respectively. In this paper,

β_{1, n}

is defined as the activation level that indicates how much the site is activated. Generally speaking, bigger

β_{1, n}

denotes stronger activation. Assuming that

β_{n}

is invariant over time, we could obtain the following state space models

{\begin{cases} [^{X_{n} (f)] T} = H (f) β_{n} (f) + w (f) \\ β_{n} (f) = β_{n} (f - 1) + v (f) \end{cases}

where

w (f)

and

v (f)

are the measurement and the process noises, respectively. In the indirect Kalman approach, Eq. (4) is independently solved to obtain

β_{n}

for all the sites to compose the spatial distribution of the activation levels [19].

2.3 Direct Kalman estimation

The state space model and Kalman approach are extended to the direct approaches to simultaneously estimate the spatial distribution of the activation levels of all the sites field by field.

2.3.1 Activation levels-boundary measurements map

To directly reconstruct the spatial distribution of the activation levels, it is essential to establish a map between the spatial distribution of the activation levels and the boundary measurements. The inverse problems in Eq. (2) of all the frames except the first frame are stacked together to compose a matrix equation as follows

Y = J X

where

Y = [y (2), y (3), \dots, y (F)]

and

X = [x (2), x (3), \dots, x (F)]

. According to the GLM, we stack together Eq. (3) of all the sites to obtain

X = B H^{T}

with

B = {[β_{1}, β_{2}, \dots, β_{N}]}^{T}

and then substitute it into Eq. (5) to get the following equation

Y = J B H^{T}

In terms of the property of Kronecker products, the left side of Eq. (6) is vectorized to the vector forms and the right side can be reformulated as follows

v e c (Y) = (H \otimes J) v e c (B)

where

v e c (\cdot)

and

\otimes

denote the vectorization computation and Kronecker products, respectively. So far, the direct map has been explicitly established as a linear equation and the temporal a priori information (HRF) is contained in

H

.

2.3.2 Estimation of the spatial distribution of the activation levels

The whole activation process denoted by Eq. (7) can be divided into many fields and expressed by the following state space models by assuming $B$ to be constant over fields

{\begin{cases} B (k) = I B (k - 1) + W (k) \\ S (k) = [H_{k} \otimes J (s)] B (k) + V (k) \end{cases}

where

B (k)

denotes

B

at the k-th field;

S (k)

is the components in

v e c (Y)

corresponding to the measurements obtained at the k-th field; The matrix

I

is an identity matrix, indicating that the state

B

will follow a random walk with zero drift over time. Thus,

B (k)

and

B (k - 1)

are the same except the process noise term

W (k) \sim Ν (0, Q)

with

0

and

Q

denoting the mean values and the covariance matrix, respectively;

H_{k}

expresses the k-th row of

H

;

J (s)

denotes the part of the sensitivity matrix corresponding to the channels covered by the s-th field in any frames;

V (k) \sim N (0, R)

is the measurement noise term with

R

representing the covariance matrix.

According to the Kalman methods, $B$ can be recursively estimated by the following 5 steps

{\begin{cases} \hat{B} (k) = I B (k - 1) \\ \hat{P} (k) = I P (k - 1) I^{T} + Q \\ K^{g} (k) = \hat{P} (k) H_{k} [^{(^{H_{k}) T} P (k) (^{H_{k}) T} + R] - 1} \\ B (k) = \hat{B} (k) + K^{g} (k) [S (k) - (^{H^{k}) T} \hat{B} (k)] \\ P (k) = [Ι - K^{g} (k) H^{k}] \hat{P} (k) \end{cases}

In the above Eq. (9),

\hat{B} (k)

and

\hat{P} (k)

are respectively the intermediate state and covariance matrix that are predicted by the state space model and will be updated by the boundary measurements at the next field;

K^{g} (k)

is the Kalman gain at the k-th field;

{(\cdot)}^{- 1}

is the inverse operation of the matrix in the bracket. When a new field is finished, the Kalman estimation is put forward one step and the spatial distribution of the activation levels is also updated once. When all the measurements are finished, the Kalman estimation ends and returns back the final spatial distribution of the activation levels.

3. Experiments and results

We first described the multi-field illumination mode and some image metrics. Then the dynamic simulative and phantom experiments were conducted to compare the indirect and direct Kalman approaches.

3.1 Multi-field illumination mode

As shown in Fig. 1, 20 sources and 12 detectors were assigned in an overlapping way to obtain the measurements necessary for DOT. All the sources were located with the same spacing of 19 mm, composing an array of 4 rows and 5 columns. At every midpoint among 4 neighboring sources, a detector was arranged, composing an array of 3 rows and 4 columns. The channels with SDS of 13.44 mm, 30.04 mm and 40.31 mm were called 1-st, 2-nd and 3-rd Nearest-Neighbor (NN) channels, respectively. In the brain tissues, the measurements are generally most sensitive to the changes in $μ_{a}$ nearby the depth equal to half of the SDS under the midpoint of the channels [20]. In this over-lapping source and detector arrangement, the 1-st NN channels mainly carried the optical information of the shallow layer of the tissues. The measured intensities of the 3-rd NN channels were much weaker than those of the 2-nd NN channels. If the 2-nd and 3-rd NN channels sharing the same detector were working simultaneously, the measurements of the 3-rd NN channels would be influenced severely by the simultaneously lit sources in the 2-nd NN channels having much stronger light intensities because of much shorter SDS. Since the 2-nd NN channels already covered the whole area and could provide sufficient overlapping measurements, we only selected them to reconstruct the spatial distribution of the activation levels to achieve the temporal resolution as high as possible.

Fig. 1 The process to illuminate the sources field by field in the multi-field illumination mode. Every channel is denoted by a red arrow whose tail and head denote the source and detector, respectively. The bright and dark red circles are the lit and unlit sources in the field, respectively. The bright and dark blue squares are the detectors working and unworking in the field, respectively.

Download Full Size | PDF

All the sources in the 2-nd NN channels were lit field by field in the multi-field illumination mode, as shown in Fig. 1. Specifically, in the first field, the sources numbered as 1, 5 and 13 (the bright red circles) were lit at the same time and the detectors (the bright blue squares) in the 2-nd NN channels associated with these sources were working simultaneously. Likewise, the remaining sources were selected to work one field after another. The 8 successive fields composed a complete frame (S = 8) and the boundary light intensities were obtained frame by frame till the experiments ended.

3.2 Image quality metrics

To quantitatively compare the performances of the the proposed direct Kalman approach and the conventional indirect Kalman approach, some metrics were adopted [21,22]. Firstly, the relative root mean square error (rRMSE) was defined below to evaluate the difference between the reconstructed and the true images.

r R M S E = \sqrt{| | β^{t r} - β^{r e c} | |_{2}^{2} / | | β^{t r} | |_{2}^{2}}

where

β^{r e c}

and

β^{t r}

were the reconstructed and true activation levels of all the sites. Secondly, the contrast-to-noise ratio (CNR) was used to evaluate the discrimination of the reconstructed targets from the background. To define CNR, we first defined the region of interest (ROI) as the region composed of sites with activation levels bigger than half of the maximum value in the image; the remaining area was the region of background (ROB). The CNR was defined below

C N R = \frac{| Mean (β_{ROI}) - Mean (β_{ROB}) |}{\sqrt{w Var (β_{R O I}) + (1 - w) Var (β_{R O I})}}

where

β_{R O I}

and

β_{R O B}

were the activation levels in the ROI and ROB, respectively;

w = | R O I | / (| R O I | + | R O B |)

with

| \cdot |

representing the area calculation operation;

Mean (\cdot)

and

Var (\cdot)

were the mean and variance operators, respectively. Besides, we drew the X-profile along the horizonal line passing through the centers of the targets and defined the spatial resolution (SR) and gray scale resolution (GSR) based on this X-profile. SR was defined as the ratio of the average of the maximum values of the two targets on the X-profiles to the minimum value on the X-profile between the two targets. GSR was defined as the ratio of the maximum value of the right target on the X-profile to that of the left target. Moreover, to evaluate the area consistency between the true and reconstructed areas of the targets, the area ratio (AR) was defined as the ratio of the area of the reconstructed ROI to that of the true targets. At last, the location error (LE) was defined as the Euclidean distance between the centroid of the true targets and the reconstructed ROI. Generally, smaller rRMSE and LE, bigger CNR and SR, GSR closer to the true values and AR closer to 1 indicated better recovered results. In fNIRS, it is more meaningful to determine the locations of the activated regions than to recover the absolute values of the activation levels [2]. Therefore, in this work we normalized the spatial distribution of the activation levels by their maximum to expediently compare the performances of the two approaches. It is worth noting that we repeated every case of simulative experiment 10 times and used the mean and variance values of the metrics to statistically compare the two approaches.

3.3 Dynamic simulative experiments

We conducted some dynamic simulative experiments where the optical properties varied over the fields to mimic the brain activation process.

3.3.1 Experimental setup

The simulative experiments were conducted using a cuboidal medium with a length (L = 120 mm), width (W = 100 mm) and height (H = 32 mm) as shown in Fig. 2. Inside the medium, two cylindrical targets were embedded to mimic the potentially activated cerebral cortex. Both of the two targets had diameters of 15 mm, heights of 12 mm and depths of 4 mm. Except the two targets, the remaining parts of the medium were called the background with $μ_{a}$ and ${μ^{'}}_{s}$ being 0.0173/mm and 1.0134/mm at the wavelength of 780 nm, respectively [9,23–26]. The initial optical properties of the targets were the same as those of the background. During the simulation process, ${μ^{'}}_{s}$ of the whole medium remained unchanged while $μ_{a}$ inside the targets varied according to $[1 + h (k) β] μ_{a 0}$ with $β$ and $μ_{a 0}$ being the activation levels of the targets and the initial $μ_{a}$ , respectively. The biggest differences in $μ_{a}$ between the background and the targets were $β μ_{a 0}$ . Besides, to simulate the task-related scalp blood flow we globally changed $μ_{a}$ in the scalp (depth $\leq$ 4 mm) according to $[1 + 0.1 h (k) β] μ_{a 0}$ which had a ten times smaller amplitude than that in the activated targets [27–29].

Fig. 2 (a) The 3-Dimensional and (b) the longitudinal section of the cuboidal medium. The small red and green cylinders on the top surface of the medium denote the the sources and detectors, respectively. The top layer denotes the scalp. The two big blue cylinders denote the targets inside the medium.

Download Full Size | PDF

We conducted the simulative experiments for $k \in [1, 40]$ when $h (k)$ finished a complete process including the activation and a subsequent process to return back to the rest state. To avoid the ‘inverse crime’ problem, the simulative boundary measurements at all the channels were generated by using the Monte Carlo (MC) method (10⁸ photon packets). However, the sensitivity matrix in the inverse problem was generated by using the finite element method (FEM) to solve the diffusion equation (DE) based on a mesh with 43073 elements and 8101 nodes [30]. For every field in the multi-field illumination mode, the MC was once used to generate the boundary measurements based on the spatial distribution of the optical properties. For the DE, the Robin boundary conditions were used for all the boundaries. Besides, the signal to noise ratio (SNR) of the measurements was proportional to the squared root of the light intensity at each channel [31]. We added Gaussian white noise to the pure simulated data, achieving the SNR = 10 dB and 20 dB for the channels with the weakest measured intensities.

3.3.2 Spatial resolution

The two targets with the same activation levels of 0.5 were both activated to investigate the SR of the proposed approaches. The X-coordinates of the target centers were 45 mm and 75 mm. The averaged spatial distribution of the activation levels reconstructed by the indirect and direct Kalman approaches over the 10 repetitive simulations were depicted in the left and middle columns in Fig. 3, respectively. Besides, the X-profiles were listed in the right column to help to compare the performances of the approaches. From Fig. 3, we could see that the direct Kalman approach separated better the two targets and suppressed better the artifacts though both two approaches recovered the rough sizes and locations of the two targets. Furthermore, the X-profile of the conventional indirect approach was wider and had less obvious highs for the two targets and less obvious lows between the two highs, compared with those of the direct Kalman approach. The metrics adopted to quantitatively compare the approaches were calculated and listed in Table 1. From this table, it was easy to see that the direct Kalman approach provided smaller rRMSE and LE, bigger CNR, and AR closer to 1. It is worth noting that the noise of SNR = 10 dB was so strong for the indirect Kalman approach that the SR for some experiments were infinite, making that the averaged SR was infinite. Besides, the averaged activation level of the left target was much less than that of the right one. However, the direct Kalman approach provided similar activation levels for both targets for all the cases. For the SNR = 20 dB, the direct Kalman approach provided bigger SR than the indirect Kalman approach. In summary, the direct Kalman approach improved almost all the metrics.

Fig. 3 The indirect Kalman (Left) and direct Kalman (Middle)-calculated images and their X-profiles (Right): (a) SNR = 10 dB and (b) SNR = 20 dB. To investigate the SR, the activation levels of the both targets were set to 0.5. The black dotted circles in the left and middle columns denoted the true locations of the targets

Download Full Size | PDF

Table 1. The metrics for the results recovered by the indirect and direct Kalman approaches in the simulative experiments with two targets having the same activation levels of 0.5.

View Table | View all tables in this article

3.3.3 Gray scale resolution

To investigate the GSR of the proposed approach, some simulative experiments for the two targets with different activation levels were conducted. The activation levels of the left and right targets were 0.5 and 1, respectively. Besides, the X-coordinates of the two target centers were 40 mm and 80 mm, farther than those in Section 3.3.2. They were easier to be separated and better mitigated the mutual influence between the two targets. The remaining experimental setup was the same as those in Section 3.3.2. The averaged spatial distribution of the activation levels computed by the indirect and direct Kalman approaches over the 10 repetitive simulative experiments, and their X-profiles were presented in Fig. 4. From Fig. 4, it was clear that both the two approaches could also calculate the correct sizes and locations of the targets. In addition, the differences between the activation levels of the two targets were recovered. The images reconstructed by the indirect Kalman approach had more artifacts than those of the direct Kalman approach. Besides, the profiles corresponding to the direct Kalman approach were narrower than those of the indirect approach. The quantitative metrics were calculated and presented in Table 2 to further compare the two approaches. Except the LE for the case of SNR = 20 dB, the direct Kalman approach provided smaller rRMSE, LE, bigger CNR, GSR closer to 2 (the true value) and AR closer to 1. In summary, this showed that the direct Kalman approach provided better GSR while improving almost all the other metrics.

Fig. 4 The indirect Kalman (Left) and direct Kalman (Middle)-calculated images and their X-profiles (Right): (a) SNR = 10 dB and (b) SNR = 20 dB. To investigate the GSR of the two approaches, the activation levels of the left and right targets were set to 0.5 and 1, respectively. The black dotted circles in the left and middle columns denoted the true locations of the targets.

Download Full Size | PDF

Table 2. The metrics for the results reconstructed by the indirect and direct Kalman approaches in the simulative experiments with two activated regions having different activation levels.

View Table | View all tables in this article

3.4 Dynamic phantom experiments

The performances of the indirect and direct Kalman approaches were compared realistically based on some dynamic phantom experiments.

3.4.1 Experimental setup

As shown in Fig. 5, a glass tank filled with the mixed solution of the Indian ink and the intralipid solution was used to mimic the brain tissues. The mixed solution was called the background solution and its optical properties were the same as those of the background in the simulative experiments (i.e. $μ_{a} = 0.0173 / m m$ , ${μ^{'}}_{s} = 1.0134 / m m$ ). The length and width of the tank were the same as the medium used in the simulative experiments while the height was 50 mm, higher than that of the medium. The depth of the mixed solution was 32 mm, equal to the height of the medium in the simulative experiments. A 3-D printed plastic cylinder (called the target cylinder henceforth) with thickness of 1 mm and diameter of 15 mm was used to mimic the activated brain regions where the Indian ink concentration of the mixed solution temporally varied while the intralipid concentration remained unchanged. It was worth noting that the bottom of the target cylinder was not next to the bottom of the glass tank and the space between them was 4 mm. The sources and detectors were also aligned in the same overlapping way as that in the simulative experiments. But they were placed below the glass tank, different from the simulative experiments where the sources and detector were placed on the top surface.

Fig. 5 The schematic of the setup of the dynamic phantoms. The constant flow pumps 1 and 2 were used respectively to pump the solution into and out from the target cylinder located in the tank filled with the background solution. The cup A and B were filled with the target and background solution, respectively. The cup C was used to collect the waste solution pumped out form the target cylinder. Switch A and Switch B were used to select a cup from which the solution was pumped into the target cylinder.

Download Full Size | PDF

Before the measurements began, the target cylinder was filled with the background solution to mimic the brain at the rest state. The Indian ink concentration of the target solution in the Cup A was 1.5 times higher than that of the background solution while the intralipid concentration was the same as that of the background solution, i.e., $μ_{a} = 0.02595 / m m$ , ${μ^{'}}_{s} = 1.0134 / m m$ . The Cup B was filled with the background solution i.e., $μ_{a} = 0.0173 / m m$ , ${μ^{'}}_{s} = 1.0134 / m m$ . The flows of the two constant flow pumps were adjusted to be equal so that the volume of the two solutions in the target cylinder remained unchanged during the experimental process. When the measurements got started, we operated the dynamic phantom equipment according to the process listed in Table 3. The mixed solution in the target cylinder would be varied gradually and successively by that contained in the Cup A and B. Correspondingly, $μ_{a}$ of the mixed solution in the target cylinder first increased and then decreased. During the whole experiments, the measurements were obtained field by field and frame by frame like in the simulative experiments.

Table 3. The process of the operations on the two switches, constant flow pumps, changes in μ_a and the simulated states of the brain.

View Table | View all tables in this article

The boundary measurements were obtained by a NIRS-DOT system based on a lock-in photon-counting technique as shown in Fig. 6 [32]. This system not only had the high sensitivity of the photon-counting technique but also had the parallel features of the lock-in technique. The 20 laser diode (LD) sources were modulated by square waves having different frequencies. When multiple sources were irradiating the phantom, the detected light intensities were decoded by the reference signals. The light intensities corresponding to every source could be separated and recovered. A photon counting head (H10682-01, HAMAMATSU, Japan) was used to detect the light intensities diffused from the phantom and had higher sensitivity than the conventional analog signal recording. An individual photon was converted into an electrical pulse. The pulse density of the signal, rather than the signal amplitude, provided the measurements of the light intensities at the input of the photon counting head [31]. This system was working on a field programmable gate array which had parallel nature and could work at enough high frequencies to capture the short impulses. The light intensities were denoted as the total number of the photon detected and decoded in an accumulation period (500 and 1000 ms). It was worth noting that the duration needed by a field was equal to the accumulation period and the duration of a frame was 8 times longer than that of a field.

Fig. 6 The schematic of the NIRS-DOT system based on the lock-in photon counting technique used in the dynamic phantom experiments.

Download Full Size | PDF

3.4.2 Results

The spatial distribution of the activation levels calculated by the indirect and direct Kalman approaches and their X-profiles were presented in Fig. 7. The images reconstructed by the indirect Kalman approach had much more artifacts than those of the direct Kalman approach. Besides, the targets reconstructed by the direct Kalman approach had sizes closer to the true ones compared with those of the indirect Kalman approach, confirmed by the narrower X-profiles of the direct Kalman approach. To investigate quantitatively the spatial resolution of the approaches from the phantom experiments with only one target, the full width half maximum (FWHM) of the X-profiles was adopted. Generally, smaller FWHM indicates better spatial resolution. The metrics including the FWHM were calculated and shown in Table 4, where the direct Kalman approach provided smaller rRMSE, FWHM and LE, bigger CNR and AR closer to 1 than those of the indirect Kalman approach.

Fig. 7 The indirect Kalman (Left) and direct Kalman (Middle)-calculated images and their X-profiles (Right): (a) accumulation period = 500 ms and (b) accumulation period = 1000 ms. The activation level of the target was 0.5. The black dotted circles in the left and middle columns denoted the true locations of the targets.

Download Full Size | PDF

Table 4. The metrics for the results reconstructed by the indirect and direct Kalman approaches in the dynamic phantom experiments with the activation level of 0.5 for accumulation period of 500 ms and 1000 ms.

View Table | View all tables in this article

4. Discussions

The primary differences between the indirect and direct Kalman approaches can be found by comparing Eqs. (4) and (8). Firstly, Eqs. (4) and (8) associate the activation levels with the changes in $μ_{a}$ and the boundary measurements, respectively. Secondly, Eq. (4) separately estimates the activation levels of all the sites while Eq. (8) simultaneously estimates the activation levels of all the sites. Thirdly, Eq. (4) assumes the changes in $μ_{a}$ constant in the same frame and updates the activation levels in the unit of frame. However, Eq. (8) assumes the changes in $μ_{a}$ constant in the same field and updates the activation levels in the unit of field.

In the phantom experiments, to make it easier to mix completely the solution in the target cylinder, we put the outlet of the tube for the solutions above the liquid level as shown in Fig. 5 and the solution from Cup A and B would arrive at the top of the target cylinder. On the other hand, the inlet of the tube for the waste solution was located at the bottom of the target cylinder so that the mixed solution was pumped out from the bottom of the target cylinder. This setup was helpful to replace thoroughly the solution in the target cylinder with that contained in Cup A and B. Nevertheless, because the exchange of the solution needed some time, we could not increase and then decrease $μ_{a}$ in the dynamic phantom experiments so fast as the simulative experiments. The dynamic phantom experiments lasted 152 s, much longer than 40 s given in the simulative experiments. The time unit of $k$ is the duration needed by a field (accumulation period), so the HRF model $h (k) = {\begin{cases} 1, 40 < k < 110 \\ 0, o t h e r w i s e \end{cases}$ and $h (k) = {\begin{cases} 1, 80 < k < 220 \\ 0, o t h e r w i s e \end{cases}$ were used for the phantom experiments with the accumulation period of 1000 ms and 500 ms, respectively. Although $μ_{a}$ in the phantom experiments was not changed accurately according to $h (k)$ in the simulative experiments, it still roughly mimicked the activation and the subsequent process to return back to the rest state. Besides, the accurate optical properties of the target cylinder, the boundary conditions at the glass tank-air, liquid-air, and liquid-target cylinder interfaces were unknown and inevitably affected the light propagation passing through them. However, these effects remained unchanged for the whole dynamic phantom experiments and could be removed in part by the differential imaging method.

For the phantom experiments, the SNR of the measurements increased as the accumulation period became longer. However, the temporal resolution would become too low and even failed to capture the rapidly changing $μ_{a}$ if the accumulation period was too long. In such situations, the direct Kalman approach was very useful to improve the temporal resolution while keeping the SNR of the measurements big enough. Thus, the speed limitation of the instruments could be relaxed by the direct Kalman approach without speeding up instruments and extra hardware cost.

In addition, from Fig. 3 and Fig. 4 of the simulative experiments we could find that the artifacts in the images reconstructed by the indirect Kalman approach of SNR = 20 dB were much fewer than those of SNR = 10 dB. However the differences between the artifacts in the images reconstructed by the direct Kalman approach for different SNRs were not obvious. Likewise, in Table 1 and Table 2 of the simulative experiments, the CNRs of the indirect Kalman approach of SNR = 20 dB were much bigger than those of SNR = 10 dB. However, the CNRs of the direct Kalman approach for different SNRs were almost equal. For the dynamic phantom experiments, the similar conclusions could also be drawn from Fig. 7 and Table 4. This indicated that the direct Kalman approach was more robust to the noise than the indirect one.

It is worth noting that in this paper we utilized a semi-3-dimensional DOT scheme that provided a similar way of data visualization to OT while persisting in the physically rigorous model-based inversion. This scheme could significantly reduce the number of the unknowns by reasonably assuming that the changes in $μ_{a}$ were invariant along the depth direction in the cerebral cortex, and therefore substantially improved the ill-posedness of the inversion process as well as relax the computation burden [33].

In this paper we didn’t consider the physiological noises such as arterial pulsations, breathing cycle, low frequency oscillations of blood pressure et al., which can be easily removed by many methods [2]. Though some researchers assumed that the interferences of the task-related scalp-hemodynamics were global [27,28], it was still necessary to simulate the local ones to investigate the performances of the proposed approaches [34]. In the scalp (depth $\leq$ 4 mm), we simulated a task-related local region where $μ_{a}$ changed as $[1 + C h (k) β] μ_{a 0}$ while ${μ^{'}}_{s}$ remained unchanged; $C$ was set to 0.1 and 1 to simulate different amplitudes. To not overlap the activated regions in the cerebral cortex with centers located at (40 mm, 50 mm) and (80 mm, 50 mm), we set the diameter to 15 mm and the X-Y coordinate of the center to (60 mm, 80 mm) for the local region in the scalp. To better remove the interferences caused by the scalp blood flow, besides to the 2-nd NN channels, we also utilized the 1-st NN channels that mainly accounted for the scalp blood flow. Correspondingly, Eq. (2) was extended as follows

[\begin{array}{l} y (f) \\ y_{1 s t} (f) \end{array}] = [\begin{array}{l} J J_{s c a l p 1} \\ 0 J_{s c a l p 2} \end{array}] [\begin{array}{l} x (f) \\ x_{s c a l p} (f) \end{array}]

where

y_{1 s t} (f)

denoted the boundary measurements of the 1-st NN channels in the f-the frame;

J_{s c a l p 1}

and

J_{s c a l p 1}

denoted the sensitivity matrixes of the boundary measurements to the changes in

μ_{a}

in the scalp;

x_{s c a l p} (f)

is the changes in

μ_{a}

in the scalp in the f-th frame. The Gaussian white noise was added into the pure simulated data obtained by MC, achieving the SNR = 10 dB and 20 dB for the channels with the weakest measured intensities. The rest of the setup was the same as those in Section 3.3.2.

The reconstructed images and metrics were presented in Fig. 8 and Table 5, respectively for $C = 0.1$ , as well as in Fig. 9 and Table 6, respectively, for $C = 1$ . It was evidently observed that both the approaches selectively reconstructed the correct locations and sizes of the activated regions in the cerebral cortex. For all the cases, the direct Kalman approach improved the images and almost all the metrics compared with the indirect Kalman approach. Furthermore, much more artifacts appeared in the results for $C = 1$ as compared to those for $C = 0.1$ . In summary, both approaches could selectively detect absorption changes originating from the cerebral hemoglobin, which were not overlapped by those from the scalp hemoglobin.

Fig. 8 The indirect Kalman (Left) and direct Kalman (Middle)-calculated images and their X-profiles (Right): (a) SNR = 10 dB and (b) SNR = 20 dB. The amplitude of the local activated region $C$ was set to 0.1. The black dotted circles in the left and middle columns denoted the true locations of the targets.

Download Full Size | PDF

Table 5. The metrics for the results reconstructed by the indirect and direct Kalman approaches in the simulative experiments where the amplitude of the local region $C$ was set to 0.1.

View Table | View all tables in this article

Fig. 9 The indirect Kalman (Left) and direct Kalman (Middle)-calculated images and their X-profiles (Right): (a) SNR = 10 dB and (b) SNR = 20 dB. The amplitude of the local activated region $C$ was set to 1. The black dotted circles in the left and middle columns denoted the true locations of the targets.

Download Full Size | PDF

Table 6. The metrics for the results reconstructed by the indirect and direct Kalman approaches in the simulative experiments where the amplitude of the local region $C$ was set to 1.

View Table | View all tables in this article

5. Conclusions

The proposed direct Kalman approach recursively estimated the spatial distribution of the activation levels from the boundary measurements and avoided solving directly the ill-posed inverse problems in DOT. The combination of the temporal a priori information and the Kalman approach was promising to improve the calculated images. More importantly, in the proposed direct Kalman approach, the spatial distribution of the activation levels was updated field by field, while they were updated frame by frame in the conventional indirect Kalman approach. Thus the temporal resolution was improved 8 times because 8 fields composed a complete frame in this paper. Then the performance improvements of the proposed approach were verified by the dynamic simulative experiments. Besides, the dynamic phantom experiments were creatively designed by using the constant flow pumps to mimic the activation and the subsequent process to return back to the rest state. Both the simulative and phantom experiments demonstrated that the proposed approach suppressed better the artifacts so that the targets were easier to be distinguished from the background (bigger CNR), provided images with less differences from the true ones (smaller rRMSE), separated better the adjacent targets (bigger SR), obtained targets with sizes closer to the true ones (AR closer to 1) and recovered more accurately the differences between the activation levels of the two targets (GSR closer to 2). In conclusion, the proposed approach improved the image quality and relaxed the burden to speed up the measurements by using the temporal a priori information.

Funding

National Natural Science Foundation of China (61575140, 61475115, 61475116, 81671728, 81571723, 81771880); Tianjin Municipal Government of China (16JCZDJC31200, 17JCZDJC32700).

Disclosures

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

References

1. T. J. Huppert, S. G. Diamond, M. A. Franceschini, and D. A. Boas, “HomER: a review of time-series analysis methods for near-infrared spectroscopy of the brain,” Appl. Opt. 48(10), D280–D298 (2009). [CrossRef] [PubMed]

2. F. Scholkmann, S. Kleiser, A. J. Metz, R. Zimmermann, J. Mata Pavia, U. Wolf, and M. Wolf, “A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology,” Neuroimage 85(Pt 1), 6–27 (2014). [CrossRef] [PubMed]

3. D. A. Boas, C. E. Elwell, M. Ferrari, and G. Taga, “Twenty years of functional near-infrared spectroscopy: introduction for the special issue,” Neuroimage 85(Pt 1), 1–5 (2014). [CrossRef] [PubMed]

4. M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” Neuroimage 63(2), 921–935 (2012). [CrossRef] [PubMed]

5. S. Lloyd-Fox, A. Blasi, and C. E. Elwell, “Illuminating the developing brain: the past, present and future of functional near infrared spectroscopy,” Neurosci. Biobehav. Rev. 34(3), 269–284 (2010). [CrossRef] [PubMed]

6. S. Tak and J. C. Ye, “Statistical analysis of fNIRS data: a comprehensive review,” Neuroimage 85(Pt 1), 72–91 (2014). [CrossRef] [PubMed]

7. C. Habermehl, J. Steinbrink, K. R. Müller, and S. Haufe, “Optimizing the regularization for image reconstruction of cerebral diffuse optical tomography,” J. Biomed. Opt. 19(9), 96006 (2014). [CrossRef] [PubMed]

8. G. Taga, K. Asakawa, A. Maki, Y. Konishi, and H. Koizumi, “Brain imaging in awake infants by near-infrared optical topography,” Proc. Natl. Acad. Sci. U.S.A. 100(19), 10722–10727 (2003). [CrossRef] [PubMed]

9. X. Wu, A. T. Eggebrecht, S. L. Ferradal, J. P. Culver, and H. Dehghani, “Fast and efficient image reconstruction for high density diffuse optical imaging of the human brain,” Biomed. Opt. Express 6(11), 4567–4584 (2015). [CrossRef] [PubMed]

10. N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15(21), 13695–13708 (2007). [CrossRef] [PubMed]

11. C. Chen, F. Tian, H. Liu, and J. Huang, “Diffuse optical tomography enhanced by clustered sparsity for functional brain imaging,” IEEE Trans. Med. Imaging 33(12), 2323–2331 (2014). [CrossRef] [PubMed]

12. J. Prakash, C. B. Shaw, R. Manjappa, R. Kanhirodan, and P. K. Yalavarthy, “Sparse Recovery Methods Hold Promise for Diffuse Optical Tomographic Image Reconstruction,” IEEE J. Sel. Top. Quant. 20, 74 (2014).

13. W. Chen, X. Wang, B. Wang, Y. Wang, Y. Zhang, H. Zhao, and F. Gao, “Lock-in-photon-counting-based highly-sensitive and large-dynamic imaging system for continuous-wave diffuse optical tomography,” Biomed. Opt. Express 7(2), 499–511 (2016). [CrossRef] [PubMed]

14. A. T. Eggebrecht, S. L. Ferradal, A. Robichaux-Viehoever, M. S. Hassanpour, H. Dehghani, A. Z. Snyder, T. Hershey, and J. P. Culver, “Mapping distributed brain function and networks with diffuse optical tomography,” Nat. Photonics 8(6), 448–454 (2014). [CrossRef] [PubMed]

15. J. Cohen-Adad, S. Chapuisat, J. Doyon, S. Rossignol, J. M. Lina, H. Benali, and F. Lesage, “Activation detection in diffuse optical imaging by means of the general linear model,” Med. Image Anal. 11(6), 616–629 (2007). [CrossRef] [PubMed]

16. Y. Zhang, D. H. Brooks, and D. A. Boas, “A haemodynamic response function model in spatio-temporal diffuse optical tomography,” Phys. Med. Biol. 50(19), 4625–4644 (2005). [CrossRef] [PubMed]

17. V. Bonomini, L. Zucchelli, R. Re, F. Ieva, L. Spinelli, D. Contini, A. Paganoni, and A. Torricelli, “Linear regression models and k-means clustering for statistical analysis of fNIRS data,” Biomed. Opt. Express 6(2), 615–630 (2015). [CrossRef] [PubMed]

18. O. M. L. C. S. Prahl, “Optical Absorption of Hemoglobin” (1999), retrieved https://omlc.org/spectra/hemoglobin/index.html.

19. X. S. Hu, K. S. Hong, S. S. Ge, and M. Y. Jeong, “Kalman estimator- and general linear model-based on-line brain activation mapping by near-infrared spectroscopy,” Biomed. Eng. Online 9(1), 82 (2010). [CrossRef] [PubMed]

20. Y. Zhang, L. Huang, N. Zhang, H. Tian, and J. Zhu, “Experimental study on the sensitive depth of backwards detected light in turbid media,” Opt. Express 26(11), 14700–14709 (2018). [CrossRef] [PubMed]

21. B. Wang, W. Wan, Y. Wang, W. Ma, L. Zhang, J. Li, Z. Zhou, H. Zhao, and F. Gao, “An L_p (0 ≤ p ≤ 1)-norm regularized image reconstruction scheme for breast DOT with non-negative-constraint,” Biomed. Eng. Online 16(1), 32 (2017). [CrossRef] [PubMed]

22. F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluorescence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008). [CrossRef] [PubMed]

23. Y. Zhan, A. T. Eggebrecht, J. P. Culver, and H. Dehghani, “Image quality analysis of high-density diffuse optical tomography incorporating a subject-specific head model,” Front. Neuroenergetics 4, 6 (2012). [CrossRef] [PubMed]

24. T. Shimokawa, T. Kosaka, O. Yamashita, N. Hiroe, T. Amita, Y. Inoue, and M. A. Sato, “Hierarchical Bayesian estimation improves depth accuracy and spatial resolution of diffuse optical tomography,” Opt. Express 20(18), 20427–20446 (2012). [CrossRef] [PubMed]

25. T. Shimokawa, T. Kosaka, O. Yamashita, N. Hiroe, T. Amita, Y. Inoue, and M. A. Sato, “Extended hierarchical Bayesian diffuse optical tomography for removing scalp artifact,” Biomed. Opt. Express 4(11), 2411–2432 (2013). [CrossRef] [PubMed]

26. X. Wu, A. T. Eggebrecht, S. L. Ferradal, J. P. Culver, and H. Dehghani, “Quantitative evaluation of atlas-based high-density diffuse optical tomography for imaging of the human visual cortex,” Biomed. Opt. Express 5(11), 3882–3900 (2014). [CrossRef] [PubMed]

27. F. Scarpa, S. Brigadoi, S. Cutini, P. Scatturin, M. Zorzi, R. Dell’acqua, and G. Sparacino, “A reference-channel based methodology to improve estimation of event-related hemodynamic response from fNIRS measurements,” Neuroimage 72, 106–119 (2013). [CrossRef] [PubMed]

28. T. Sato, I. Nambu, K. Takeda, T. Aihara, O. Yamashita, Y. Isogaya, Y. Inoue, Y. Otaka, Y. Wada, M. Kawato, M. A. Sato, and R. Osu, “Reduction of global interference of scalp-hemodynamics in functional near-infrared spectroscopy using short distance probes,” Neuroimage 141, 120–132 (2016). [CrossRef] [PubMed]

29. T. Funane, H. Atsumori, T. Katura, A. N. Obata, H. Sato, Y. Tanikawa, E. Okada, and M. Kiguchi, “Quantitative evaluation of deep and shallow tissue layers’ contribution to fNIRS signal using multi-distance optodes and independent component analysis,” Neuroimage 85(Pt 1), 150–165 (2014). [CrossRef] [PubMed]

30. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [PubMed]

31. W. Becker, Advanced Time-Correlated Single Photon Counting Techniques (Springer Science & Business Media) (2005).

32. X. Ding, B. Wang, D. Liu, Y. Zhang, J. He, H. Zhao, and F. Gao, “A three-wavelength multi-channel brain functional imager based on digital lock-in photon-counting technique,” SPIE Proc. 10480, 104800S (2018).

33. F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “Optical tomographic mapping of cerebral haemodynamics by means of time-domain detection: methodology and phantom validation,” Phys. Med. Biol. 49(6), 1055–1078 (2004). [CrossRef] [PubMed]

34. N. M. Gregg, B. R. White, B. W. Zeff, A. J. Berger, and J. P. Culver, “Brain specificity of diffuse optical imaging: improvements from superficial signal regression and tomography,” Front. Neuroenergetics 2(14), 14 (2010). [CrossRef] [PubMed]

SNR/dB	10		20
Method	Indirect Kalman	Direct Kalman	Indirect Kalman	Direct Kalman
SR	–	–	1.51 ± 0.295	–
CNR	2.43 ± 1.341	6.03 ± 0.794	4.81 ± 0.29.	6.85 ± 0.185
rRMSE	1.19 ± 0.178	0.71 ± 0.053	0.94 ± 0.064	0.66 ± 0.011
AR	2.45 ± 0.620	1.17 ± 0.158	2.35 ± 0.170	1.31 ± 0.055
LE/mm	8.87 ± 2.855	2.92 ± 1.214	3.99 ± 0.844	2.13 ± 0.341

SNR/dB	10		20
Method	Indirect Kalman	Direct Kalman	Indirect Kalman	Direct Kalman
GSR	4.65 ± 3.524	2.59 ± 0.402	2.82 ± 0.334	2.47 ± 0.057
CNR	3.34 ± 0.931	5.84 ± 0.360	4.41 ± 0.218	5.91 ± 0.049
rRMSE	1.127 ± 0.129	0.65 ± 0.026	0.89 ± 0.019	0.65 ± 0.003
AR	2.29 ± 0.321	1.27 ± 0.086	1.85 ± 0.053	1.33 ± 0.022
LE/mm	7.17 ± 3.964	2.47 ± 0.555	2.07 ± 0.198	2.56 ± 0.098

Time/s	The operations on the two switches and constant flow pumps.	The changes in $μ_{a}$	The simulated state of the brain
1-30	Switch A and Switch B were both closed. Both the constant flow pumps were not working.	Remained unchanged.	The rest state before the onset of the stimulations.
31-60	Switch A and Switch B were open and closed, respectively. Both the constant flow pumps were working.	Gradually increased (smaller than or equal to the value of the target solution).	The activation process of the brain after the onset of the stimulation.
61-80	Switch A and Switch B were both closed. Both the constant flow pumps were not working.	Remained unchanged.	The brain remained the activated state.
81-110	Switch A and Switch B were closed and open, respectively. Both the constant flow pumps were working.	Gradually decreased (bigger than or equal to the value of the background solution).	The brain returned back to the rest state after the stimulation ended.
111-152	Switch A and Switch B were both closed. Both the constant flow pumps were not working.	Remained unchanged.	The brain remained the rest state.

Accumulation period /ms	500		1000
Method	Indirect Kalman	Direct Kalman	Indirect Kalman	Direct Kalman
FWHM/mm	22.0	21.0	25.0	20.0
CNR	5.81	8.02	6.30	7.93
rRMSE	1.35	0.98	1.14	1.03
AR	2.78	2.11	2.33	2.22
LE/mm	6.16	6.14	6.93	6.13

SNR/dB	10		20
Method	Indirect Kalman	Direct Kalman	Indirect Kalman	Direct Kalman
SR	–	–	6.10 ± 5.034	–
CNR	2.72 ± 0.854	3.90 ± 0.730	4.10 ± 0.211	5.09 ± 0.229
rRMSE	1.15 ± 0.197	0.91 ± 0.094	1.05 ± 0.065	0.73 ± 0.011
AR	2.78 ± 0.882	1.76 ± 0.344	2.70 ± 0.313	1.52 ± 0.056
LE/mm	6.03 ± 1.994	5.93 ± 2.074	1.89 ± 0.886	3.15 ± 0.193

A Kalman-based tomographic scheme for directly reconstructing activation levels of brain function

Abstract

1. Introduction

2. Methods

2.1 HRF model

2.2 Indirect Kalman estimation

2.2.1 DOT reconstruction for each frame

2.2.2 Estimation of the spatial distribution of the activation levels

2.3 Direct Kalman estimation

2.3.1 Activation levels-boundary measurements map

2.3.2 Estimation of the spatial distribution of the activation levels

3. Experiments and results

3.1 Multi-field illumination mode

3.2 Image quality metrics

3.3 Dynamic simulative experiments

3.3.1 Experimental setup

3.3.2 Spatial resolution

3.3.3 Gray scale resolution

3.4 Dynamic phantom experiments

3.4.1 Experimental setup

3.4.2 Results

4. Discussions

5. Conclusions

Funding

Disclosures

References

Cited By

Figures (9)

Tables (6)

Equations (12)

Optics Express

SNR/dB	10		20
Method	Indirect Kalman	Direct Kalman	Indirect Kalman	Direct Kalman
SR	–	–	1.40 ± 0.120	121 ± 73.40
CNR	2.35 ± 1.096	3.71 ± 0.445	3.16 ± 0.197	4.76 ± 0.059
rRMSE	1.18 ± 0.142	0.98 ± 0.083	1.25 ± 0.056	0.81 ± 0.006
AR	2.37 ± 0.773	1.86 ± 0.322	3.25 ± 0.160	1.49 ± 0.038
LE/mm	11.0 ± 5.062	8.37 ± 2.110	8.38 ± 0.633	5.88 ± 0.480

SNR/dB	10		20
Method	Indirect Kalman	Direct Kalman	Indirect Kalman	Direct Kalman
SR	–	–	1.40 ± 0.120	121 ± 73.40
CNR	2.35 ± 1.096	3.71 ± 0.445	3.16 ± 0.197	4.76 ± 0.059
rRMSE	1.18 ± 0.142	0.98 ± 0.083	1.25 ± 0.056	0.81 ± 0.006
AR	2.37 ± 0.773	1.86 ± 0.322	3.25 ± 0.160	1.49 ± 0.038
LE/mm	11.0 ± 5.062	8.37 ± 2.110	8.38 ± 0.633	5.88 ± 0.480