==== Front
J Biomed Opt
J Biomed Opt
JBOPFO
JBO
Journal of Biomedical Optics
1083-3668
1560-2281
Society of Photo-Optical Instrumentation Engineers

10.1117/1.JBO.28.6.066005
JBO-230064GR
230064GR
Imaging
Paper
Neighbor-based adaptive sparsity orthogonal least square for fluorescence molecular tomography
Yi Huangjian ab†yhj2014@nwu.edu.cn

Ma Sihao ab†msh15109271526@163.com

https://orcid.org/0009-0003-0876-6123
Yang Ruigang abyangrg@stumail.nwu.edu.cn

Zhong Sheng abszhong@nwu.edu.cn

https://orcid.org/0000-0003-2035-874X
Guo Hongbo abguohb@nwu.edu.cn

He Xuelei abxueleihe@nwu.edu.cn

He Xiaowei ab*hexw@nwu.edu.cn

Hou Yuqing ab*houyuqin@nwu.edu.cn

a Northwest University, School of Information Sciences and Technology, Xi’an, China
b The Xi’an Key Laboratory of Radiomics and Intelligent Perception, Xi’an, China
* Address all correspondence to Xiaowei He, yhj2014@nwu.edu.cn; Yuqing Hou, houyuqin@nwu.edu.cn
† These authors contributed equally to this study.

29 6 2023
6 2023
29 6 2023
28 6 06600520 3 2023
8 6 2023
13 6 2023
© 2023 The Authors
2023
The Authors
https://creativecommons.org/licenses/by/4.0/ Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Abstract.

Significance

Fluorescence molecular tomography (FMT) is a promising imaging modality, which has played a key role in disease progression and treatment response. However, the quality of FMT reconstruction is limited by the strong scattering and inadequate surface measurements, which makes it a highly ill-posed problem. Improving the quality of FMT reconstruction is crucial to meet the actual clinical application requirements.

Aim

We propose an algorithm, neighbor-based adaptive sparsity orthogonal least square (NASOLS), to improve the quality of FMT reconstruction.

Approach

The proposed NASOLS does not require sparsity prior information and is designed to efficiently establish a support set using a neighbor expansion strategy based on the orthogonal least squares algorithm. The performance of the algorithm was tested through numerical simulations, physical phantom experiments, and small animal experiments.

Results

The results of the experiments demonstrated that the NASOLS significantly improves the reconstruction of images according to indicators, especially for double-target reconstruction.

Conclusion

NASOLS can recover the fluorescence target with a good location error according to simulation experiments, phantom experiments and small mice experiments. This method is suitable for sparsity target reconstruction, and it would be applied to early detection of tumors.

Keywords:

fluorescence molecular tomography
image reconstruction
neighbor strategy
orthogonal least square
support sets
National Natural Science Foundation of China61906154 61971350 61901374 11871321 12271434 66201459 62271394 Key Research and Development Projects of Shaanxi Province2020SF-036 Youth Innovation Team of Shaanxi Provincial Department of Education21JP123 running-headYi et al.: Neighbor-based adaptive sparsity orthogonal least square for fluorescence…
==== Body
pmc1 Introduction

Fluorescence molecular tomography (FMT) is a promising imaging technology that noninvasively and dynamically offers a 3D visualization of the biological process in-vivo at the cellular and molecular levels.1–4 Consequently, it greatly promotes its application in small animal research and preclinical diagnosis.4,5 However, the reconstruction of FMT is severe ill-posed caused by the strong scattering of near-infrared photons propagation in biological tissues. In addition, the number of measurements available is typically much smaller than the number of unknowns, which aggravate the under-determination of the reconstruction.6,7

To alleviate the ill-posedness, great effort has been made on the reconstruction algorithms. Effective regularization methods are developed, such as Tikhonov regularization, Lp-norm (0<p≤1) regularization, total variation regularization, and hybrid regularization.8–12 FMT is widely used in the early tumor detection. The tumor is small and sparse compared to the entire imaging domain. Tikhonov regularization, however, often provides an over-smoothed result and is absence of details in the local features.13,14 As a consequence, sparsity regularization methods are beneficial in this case. L0-norm is an ideal sparsity regularizer, which can provide the sparsest solution. In this work, L0-norm regularization has been utilized to establish the mathematical model of the inverse problem and a greedy algorithm is developed to solve this model. Greedy algorithm is a signal recovery technology for an underdetermined linear system following the principle that the sparsity can be exploited to recover signals from few samples if the system matrix satisfies the restricted isometry property.15,16 For example, the previously proposed orthogonal matching pursuit (OMP) is simple and effective. It selects the atom which is the most related to the current residual at each iteration.17 However, an inaccurate result will be obtained when the columns in the system matrix have strong correlation. Compressive sampling matching pursuit (CoSaMP) is designed for better atom selection and adopts the backtracking strategy to improve the accuracy.18 These algorithms need to know the sparsity in advance. However, the sparsity is usually unknown in the reconstruction process. Thong et al. proposed the sparsity adaptive matching pursuit algorithm to reconstruct the signal without prior information of the sparsity.19,20 This algorithm has provided the sparsity by increasing the size of the support set with a fixed step size. But it will generate an inaccurate sparsity, which further leads to incorrect results.20 There are substantial ongoing researches related to the field of adaptive sparse method. These studies have been applied to FMT reconstruction and achieved great success.21–23 In addition, the greedy algorithms mentioned above more easily get to a local optimal solution instead of the global optimal solution. In general, these methods cannot do very well in multi-target reconstruction.17

In practical applications, the real sparsity cannot be known in advance. Orthogonal least square (OLS) algorithm needs to input the sparsity manually based on experience.24 In each iteration, only one column index with the highest correlation is selected to be incorporated into the support set, which result in a worse reconstruction accuracy. In this work, neighbor-based adaptive sparsity OLS (NASOLS) is developed. This algorithm can provide the sparsity adaptively without knowing sparsity beforehand. A neighbor expansion strategy based on tetrahedron element is presented to provide the support set, which can preserve the local spatial structure information. This effective scheme has improved the resolution of double-target reconstruction. To verify the performance of NASOLS algorithm, numerical simulation experiments, phantom experiment, and small animal experiments were designed in this study. Adaptive sparsity OLS (ASOLS) without neighbor strategy, OLS, CoSaMP, and generalized OMP (gOMP) are selected as the comparative methods. The experimental results show that the NASOLS algorithm has the potential for reconstruction of FMT, especially for multi-target reconstruction.

The outline of this paper is listed as follows. The mathematical model of reconstruction problem and adaptive neighbor OLS algorithm are described in Sec. 2; Sec. 3 presents the results of simulation experiments, phantom experiments, and small animal experiments with the OLS and comparative methods; in Sec. 4, the discussions of the performance of NASOLS and conclusions of this work are shown.

2 Methods

2.1 Mathematical Model of Reconstruction Problem

In general, the light propagation is described by the radiative transfer equation.25 However, it is difficult to solve this complex integro-differential equation. Diffusion equation (DE) is a popular approximation in highly scattering biological tissues.26,27 For continuous wave FMT model, two coupled diffusion equations should be used to describe the behaviors of the excited and the emitted light. And finite element method (FEM) is utilized to solve the coupled DEs.28 Finally, a linear relationship between the surface photon intensity of the emitted light and the distribution of fluorescent yield inside the object can be constructed as follows: Φm=AX,(1)

where Φm is the surface photon intensity, which is measured by detectors. A is the system matrix. X denotes the distribution of fluorescent yield to be reconstructed. More detailed descriptions can be found in Refs. 12, 13, and 29. Solving Eq. (1) suffers from being ill-posed. The L0 regularization is utilized, and Eq. (1) is rewritten as follows based on compressed sensing theory:30 minx ‖X‖0,subject to  ‖AX−Φm‖22<ε,(2)

where ‖X‖0 is the L0-norm, and ε is a given threshold.

2.2 Neighbor-Based Adaptive Orthogonal Least-Squares Method

In this section, we present a neighbor-based adaptive OLS (NAOLS) to solve Eq. (2). The specific process of NAOLS is as follows:

1. Define column index selection formula: the critical point of this method is to select the most appropriate column indexes to add into the support set. This method sequentially projects columns of A onto a residual vector. Specifically, in the i’th iteration, the proposed method chooses new column indexes according to the equation J=arg max ‖qj‖2,(3)

where qj=ajTriajTtjitji,(4)

tji+1≜aj−∑l=1iajTul‖ul‖22ul=tji−tjiTui‖ui‖22ui,(5)

where, ri represents the residual vector of the i’th iteration, tj0=aj, aj is the j’th column of A, ui represents the orthogonal basis vector of the i’th iteration.

2. Calculate sparsity Ki and parameter Li: conventional OLS methods need to preset the sparsity coefficient K empirically, whereas sparsity may not be available in many practical applications, which greatly limits the practicality of FMT reconstruction. Here, based on the initial sparsity K0, we adopt the sparsity adaptive strategy. A nonlinear function was utilized to adjust the adaptive step size. It can be divided into two parts: (a) fast estimation of large step size in the initial stage, and (b) the completion stage is gradually approached in small steps. It is expressed as follows: Ki=Ki−1+⌈K0·(1i+1)2⌉,(6)

Li=Li−1−⌈L0·(1i+1)2⌉,(7)

where I is the iteration number, Li is the number of columns added to the support set in i’th iteration, L0 is the initial parameter. ⌈·⌉ means round up. In initial stage, the step size takes a suitable initial value, and each time Ki and Li are changed by Eqs. (6) and (7). The advantage of step change of nonlinear function is that the step change in the initial stage is large, and the step change in the completion stage is small, as shown in Fig. 1(a). This means the proposed method can approach the real sparsity more quickly. According to Eqs. (6) and (7), the sparsity Ki increases, the parameter Li decreases.

3. Update the support set: considering that tumor always grows in a certain region, which when compared to the whole domain is sparse enough. To utilize the local spatial structure information sufficiently, we proposed a novel neighbor strategy based on the finite element theory. The neighbor set is constructed according to the tetrahedral element structure, as shown in Fig. 1(b). Suppose Si represents the support set generated from the i’th iteration. Let any node Gk∈Si, then the nodes with tetrahedral edge connection relationship with node Gk are the elements of the neighbor set. According to this rule, finding the neighborhood of all nodes in Si and uniting all neighborhood sets to form the final neighbor set Ei. The neighbor operator process is described by N(·), namely Ei=N(Si). Selecting Li+1 nodes from Ei according to Eq. (3) and adding them into the support set Si to form a new support set Si+1. That is, Si+1=Li+1(Ei)∪Si, where Li+1(Ei) indicates selecting Li+1 nodes from Ei.

4. Update the residual: setting the ri as the residual vector represents the i’th iteration, where the residual vector required for the next iteration was formed as ui+1≜qj,ri+1=ri−ui+1.(8)

5. Terminate the iteration: circulating the Steps (2) to (4), in each cycle i=i+1. The algorithm terminates when the halting condition is satisfied. The error reaching an acceptable range ri<ε, the first Ki columns of the final support set are solution sets.

Fig. 1 (a) Schematic diagram of exponential function. (b) Neighbors of red node are black nodes, the green nodes are the irrelevant nodes of the red node.

Algorithm 1 Neighbor-based adaptive sparsity orthogonal least squares method.

Input: System matrix A, Surface photon distribution Φm,	
threshold ε, initial sparsity K0, initial parameter L0	
Output: recovered support Si, estimated signal x^	
Initialize: support set S0=Ø, iteration number i=1,	
residual vector ri=Φm, tji=aj, qj=ajTriajTtjitji, E1=N(S1)	
S1 was generated from choosing L0 columns by Eq. (3)	
While: ‖ri‖2≥ε do	
 1. Calculate sparsity Ki and parameter Li by Eqs. (6) and (7);	
 2. Select Li+1 largest terms from Ei according to Eq. (3) into Si, utilize neighbor operator N(·) to obtain neighbor set Ei;	
 3. i=i+1;	
 4. Si+1=Li+1(Ei)∪Si, Li+1(Ei) indicates to select Li+1 nodes from Ei;	
 5. Perform (8) Li times to update {Uz}z=1Li and ri;	
 6. tji+1=tji−tjiTui‖ui‖22ui;	
end while	
 7. the solution x^=Ki(Si)	

3 Experiments and Results

In this section, experiments were carried out to evaluate the performance of NASOLS, including numerical simulation experiments with single target and double targets, phantom experiments with double targets, and small animal implanted experiments. The sparsity levels of ASOLS, OLS, CoSaMP, and gOMP are set to 4, 10, 8, and 6. The error in halting condition for different methods is set to 1e−8. The experiment codes were written in MATLAB and all these processes were executed on a desktop computer with 3.20 GHz Intel Processor i7-8700 CPU and 16 GB RAM. To quantify the reconstruction performance, location error (LE), normalized root-mean-square error (NRMSE), and contrast-to-noise ratio (CNR) were adopted in this study.31,32 The detailed formula descriptions can be found in the literature.32 In general, a high-quality reconstructed image possesses LE, NRMSE, close to 0 and a high CNR value.

As shown in Fig. 2, the initial value of L0 have affected the recovered results. The initial value of K0 and L0 was determined based on the results of NASOLS. We set K0 increased by 6 every time from 6 to 30, and L0 increased by 5 every time from 5 to 25. Figure 2 shows the quantitative reconstruction results with different L0 and K0. It can be seen that when L0=10, K0=6, NASOLS achieved small LE, small NRMSE and large CNR. So the initial values of K0 and L0 in NASOLS are set to 6 and 10.

3.1 Numerical Experiments

A 3D digital mouse model was employed,31 as shown in Fig. 3(a). The torso section of the mouse, including heart, lung, liver, stomach, kidneys, and muscle, was the investigated region. The optical parameters of the main organs were the same as used in literature.13,33,34 Figure 3(b) shows the reconstruction mesh of the torso section for the inverse problem. It contains 5890 nodes and 29,308 tetrahedral elements. Here, the fluorescent target was excited by eight point sources at different positions in sequence as shown in Fig. 3(c). The black dots represent positions of the excitation point sources, which were modeled as isotropic point sources located one mean free path of photon transport beneath the surface. For each excitation source, the surface data on the opposite side with a 120 deg field of view (FOV) were measurable.

Fig. 2 Quantitative comparison of parameter L0 and K0 test. (a)–(c) LE, NRMSE, and CNR for different K0 and L0. (d) Fixed K0=6, the LE for different L0.

Fig. 3 (a) 3D digital mouse model, (b) reconstruction mesh for the inverse problem, and (c) plane of excitation sources. The black points are the location of the isotropic point sources. For each excitation source, fluorescence is detected at the opposite side with a 120 deg FOV.

3.1.1 Single target reconstruction

Figure 4 shows the reconstructed results with eight projections measurements. In the single target experiment, a sphere with radius of 1 mm was placed in the liver with center at (18.0, 12.0, 16.4 mm). The fluorescent yield of the target was set to be 0.05  mm−1. In the forward problem, the digital mouse model is a mesh with 24,041 nodes and 1,27,248 tetrahedral elements. The first row is the 3D views of the FMT images reconstructed using the NASOLS, ASOLS, OLS, CoSaMP, and gOMP, respectively. The red sphere is the real target, and the green area is the reconstructed target, respectively. The second column is the coronal slices, the third column is the sagittal slices, and the fourth column is the transverse slices. The white circle denotes the actual position of the target. The corresponding quantitative indicators are given in Table 1. It is obvious that the proposed algorithm achieved the smallest LE and NRMSE and largest CNR. It provided better results than the other four approaches, which indicated that our method has the ability for the single target reconstruction.

Fig. 4 (a)–(e) Recovered results by 3D views. The red spheres represent the actual fluorescent target and the green areas denote the reconstructed one. (a1)–(e1) Sagittal slices, (a2)–(e2) coronal slices, and (a3)–(e3) transverse slices. The white circles denote the real target.

Table 1 Quantitative results of single target reconstruction experiment.

Method	Center (mm)	LE (mm)	NRMSE	CNR	
NASOLS	(17.81,12.06,16.02)	0.41	0.0012	2.29	
ASOLS	(17.51,12.16,16.65)	0.57	0.0012	1.61	
OLS	(18.32,12.18,17.03)	0.73	0.0016	1.58	
CoSaMP	(17.59,12.12,16.92)	0.66	0.0014	0.84	
gOMP	(17.09,12.23,16.54)	0.89	0.0014	1.52	
Note: The best results are in bold.

As shown in Fig. 5, different projection numbers (2, 4, 8, 16) simulations have been conducted to show the relationship between the number of projections and the reconstruction accuracies of different algorithms. In the experiments, take projection number of 16 as example, NASOLS provides the smallest LE, NRMSE and the largest CNR compared to the other four methods. With the increasing of projection number, the recovered results became more accurate for most cases. However, NASOLS provided largest CNR with 8 projections, not 16 projections, so did ASOLS, OLS, and gOMP. The reason may be that the redundant information appeared with the increasing of projection numbers, which in turn affected the reconstruction accuracy.

Fig. 5 The quantitative analysis of single target with different projection numbers by NASOLS, ASOLS, OLS, CoSaMP, and gOMP. (a)–(c) LE, NRMSE, and CNR of five methods.

Experiments with different levels of noise (5%, 10%, 15%, 20%, and 25% Gaussian noise) of the measurements were carried out to study the effect of noise on the reconstruction results. It was clear that the recovered results became worse with the increasing of the noise levels, as shown in Fig. 6. However, NASOLS performed best of five methods according to LE, NRMSE, CNR.

Fig. 6 Quantitative results of single target with different levels of noise (5%, 10%, 15%, 20%, and 25% Gaussian noise). (a)–(c) LE, NRMSE, and CNR of five methods.

3.1.2 Double-target reconstruction

In this section, double targets with different edge to edge distance (EED) were designed to further show the performance of the algorithm. Two spheres with radius of 1 mm were placed in the center of (16.0, 14.5, 14.0 mm) and (20.0, 14.5, 14.0 mm), respectively.

The EED of the two spheres is 2 mm, and the forward mesh had 24,275 nodes and 1,28,635 tetrahedral elements. When the EED was set to be 3 mm, the center of two targets were (15.0, 14.5, 14.0 mm) and (20.0, 14.5, 14.0 mm). Its forward mesh of digital mouse concluded 19,306 nodes and 1,01,362 tetrahedral elements. When the EED was set to be 4 mm, the center of two targets were (16.0, 14.5, 14.0 mm) and (22.0, 14.5, 14.0 mm). In total, 19,339 nodes and 1,01,562 tetrahedral elements were included in the forward mesh. The difficulty of reconstruction increases with the decrease of EED, so the double-target experiments with different EEDs could reflect the spatial resolution. Figures 7(a)–7(e) show the reconstructed results when the EED was 2 mm. When the EED was 2 mm, ASOLS, gOMP, and OLS cannot distinguish two targets. Both of NASOLS and CoSaMP had the ability to distinguish two targets. But, there were a lot of artifacts by CoSaMP reconstruction. Figures 7(f)–7(j) show the recovered results when the EED was 3 mm. gOMP and OLS could only reconstruct one target, which was close to the middle of the two targets. ASOLS could reconstruction one target accuracy. NASOLS and CoSaMP could reconstruct two targets. However, the artifact of one target reconstructed by CoSaMP is large, the location error of two targets by CoSaMP is larger than NASOLS. Figures 7(k)–7(o) show the reconstruction results when the EED was 4mm. OLS could only recover one target. NASOLS could recover the two targets very well, which was better than the other four algorithms. From the quantitative results in Table 2, NASOLS had the smallest LE, the smallest NRMSE, and the largest CNR when EED was 2, 3, or 4 mm.

Fig. 7 Double-target recovered results with different EEDs. (a)–(e) 3D view results when the EED is 2 mm by NASOLS, ASOLS, OLS, CoSaMP, and gOMP. (a1)–(e1) Transverse slices at z=14.0  mm. (f)–(j) 3D view results when the EED is 3 mm by the five methods. (f1)–(j1) Transverse slices at z=14.0  mm. (k)–(o) 3D view results when the EED is 4 mm. (k1)–(o1) Transverse slices at z=14.0  mm. The white circles represent the real targets.

Table 2 Quantitative results of different EED experiments.

EED (mm)	Methods	Center (mm)	LE1 (mm)	LE2 (mm)	NRMSE	CNR	
2	NASOLS	(19.10,13.57,14.09)	1.29	1.02	0.0035	1.84	
(16.20,13.50,13.82)	
ASOLS	(18.79,13.55,14.15)	1.54	2.13	0.0041	1.33	
(17.76,13.29,13.91)	
OLS	(18.92,13.87,14.89)	1.57	3.51	0.0042	1.32	
(13.00,15.52,12.50)	
CoSaMP	(18.43,14.02,14.11)	1.64	1.28	0.0064	0.82	
(17.01,13.83,13.56)	
gOMP	(18.39,14.34,15.53)	2.20	2.01	0.0049	0.63	
(17.83,13.67,13.65)	
3	NASOLS	(19.60,14.22,14.85)	1.00	0.88	0.0027	2.28	
(15.58,13.95,14.46)	
ASOLS	(19.37,13.59,13.82)	1.12	1.65	0.0033	1.89	
(16.35,13.59,13.74)	
OLS	(18.60,13.79,13.88)	1.51	2.89	0.0035	1.88	
(17.85,14.14,13.68)	
CoSaMP	(19.22,13.39,13.59)	1.43	1.01	0.0053	1.01	
(15.07,14.55,15.06)	
gOMP	(18.22,14.08,13.77)	1.81	1.82	0.0031	0.95	
(16.78,14.70,13.71)	
4	NASOLS	(22.13,7.10,14.10)	0.44	0.70	0.0011	2.83	
(16.66,7.77,14.02)	
ASOLS	(20.89,8.02,13.95)	1.22	0.97	0.0023	2.39	
(16.93,7.76,14.08)	
OLS	(20.07,7.78,14.21)	1.95	1.46	0.0025	2.29	
(17.07,8.03,14.83)	
CoSaMP	(20.73,7.46,13.92)	1.26	0.89	0.0022	1.62	
(16.70,7.27,13.49)	
gOMP	(20.51,7.52,14.61)	1.45	1.47	0.0027	1.55	
(17.19,6.98,13.32)	
Note: The best results are in bold.

Double targets (EDD at 4 mm) with different projection numbers (2, 4, 8, 16) were also carried out and the recovered results are shown in Fig. 8. Take projection number of 16 as an example, NASOLS has also performed best according to LE, NRMSE, and CNR compared to the other four methods. In addition, we have found that when the projection number is too small, ASOLS and OLS cannot distinguish two targets. For example, OLS could only recover one target with projection numbers 2 and 4, so did ASOLS with 2 projections. Then, it can be seen that there is no result for these cases, as shown in Fig. 8(a).

Fig. 8 The quantitative analysis of double targets with different projection numbers by NASOLS, ASOLS, OLS, CoSaMP, and gOMP. (a)–(d) LE1 (Target 1), LE2 (Target 2), NRMSE, and CNR of five methods.

Double-target (EDD at 4 mm) experiments with different levels of noise (5%, 10%, 15%, 20%, and 25% Gaussian noise) of the measurements have also been conducted. From Fig. 9(a), it is clear that OLS cannot distinguish two targets if there was noise on the data. It indicated that OLS is sensitive to noise. When noise level is 15%, 20%, and 25%, ASOLS cannot distinguish either. Anyway, NASOLS provided best results according to LE, NRMSE, and CNR.

Fig. 9 The quantitative analysis of double targets with different levels of noise by NASOLS, ASOLS, OLS, CoSaMP, and gOMP. (a)–(d) LE1 (Target 1), LE2 (Target 2), NRMSE, and CNR of five methods.

3.2 Real Experiments

In this section, physical phantom and small animal experiments were conducted to further evaluate the algorithm. First, the real fluorescence data was collected by the FMT/micro-computed tomography (CT) dual modality imaging system, as shown in Fig. 10. It is a full-angle, non-contact imaging system. This system could obtain fluorescence signals and structural information simultaneously. The excitation light with 670 nm wavelength was provided by a continuous wave laser source (CrystaLaser, Reno, Nevada, and Model No. CL671-050-O). The rotational stage could rotate by the computer. And a highly sensitive charge-coupled device (CCD) camera (Princeton Instruments PIXIS 2048B, Roper Scientific, Trenton, New Jersey, United States), which was cooled to −80°C to reduce the effects of thermal noise, was used to collect the emitted light. Micro-CT system includes an x-ray detector (1512N-C90-HRCC, Dexela, United States) and an x-ray generator (L9181-02 MT2195, Hamamatsu Photonics, China). The control module has a 360-deg motorized turntable (RAK100, Zolix, China) and a controller (Zolix Instruments Co., Beijing, China). The x-ray source and x-ray detector are placed on the same line. In addition, The CCD was perpendicular with this line. Different projections of fluorescence images could be obtained by rotating the stage. In general, optical images were collected first, and then the CT data were followed by the micro-CT systems. Main organs could be segmented from the CT data, and a heterogeneous mouse model could be obtained. The optical parameter of main organs were the same as used in literature.13 A mesh for the inverse problem was provided by discretizing heterogeneous mouse model. The fluorescence image was mapped to the mesh to obtain the measurement data.35

Fig. 10 FMT/micro-CT dual modality imaging system.

3.2.1 Physical phantom experiment

The phantom, a cube with a side length of 25 mm, was made of polyformaldehyde. Its optical parameters for both excitation and emission wavelengths are the same as illustrated in Ref. 36, which were determined by diffuse optical tomography. Two small holes with 1 mm radius were drilled to emplace the Cy5.5 solution. The Cy5.5 solution was injected into the holes, which became two cylindrical targets with a height of 2.0 mm. Their centers were (5.0, 10.0, 15.0 mm) and (5.0, 15.0, 15.0 mm) with EED of 3 mm, as shown in Figs. 11(a) and 11(b). The mesh for the inverse reconstruction was discretized into 6031 nodes and 32,430 tetrahedral elements. Here, four projections were acquired by rotating the phantom with an angular increment of 90 deg. The fluorescent targets were excited by point sources from four different positions at the z=15.0  mm plane and CCD acquired data at four different views, as shown in Fig. 11(d).

Fig. 11 (a) Physical phantom; (b) geometric structure of the cubic phantom with double targets (2 mm in diameter and 2 mm in height); (c) reconstruction mesh for the inverse problem; (d) x-y view on the z=15.0  mm plane, where the black dots represent the excitation point source positions. Four degrees show the direction of the CCD camera during data acquisition.

Figure 12 shows reconstruction results with 3D and transverse views. It is clear that the proposed method could distinguish two targets while ASOLS, OLS, and gOMP could not do from the 3D views, shown in Figs. 12(a)–12(e). For CoSaMP method, one of the reconstructed targets was in the middle of two actual targets, and the artifacts appeared, which would mislead researchers. The transverse views in Figs. 12(f)–12(j) had also given the similar conclusion. From Table 3, the LE by NASOLS is the smallest, NRMSE is the smallest, and CNR is the largest. This means NASOLS provided the best results of the five methods.

Fig. 12 The phantom experiments results. (a)–(e) 3D view of results using NASOLS, ASOLS, OLS, CoSaMP, and gOMP, respectively. (f)–(j) Transverse views at Z=15  mm slices.

Table 3 Quantitative results of double-target reconstruction in physical phantom experiment.

Methods	Center (mm)	LE1 (mm)	LE2 (mm)	NRMSE	CNR	
NASOLS	(4.48,14.92,15.44)	0.68	0.80	0.0082	2.82	
(5.36,10.65,14.70)	
ASOLS	(4.38,11.62,14.53)	1.81	—	0.0092	2.42	
—	
OLS	(5.44,13.50,14.60)	1.61	—	0.0099	2.42	
—	
CoSaMP	(4.48,14.92,15.44)	0.68	2.37	0.0087	1.56	
(5.62,12.21,14.43)	
gOMP	(4.37,12.32,14.14)	2.88	0.96	0.011	2.31	
(5.09,10.95,15.06)	
Note: The best results are in bold.

The phantom experiments with different projections (1, 2, 3, 4) were carried out too. Figure 13 showed the recovered results for five methods. With the increasing of projection numbers, the reconstruction results were improved according to the evaluation index for five methods. It is clear that the quantitative indicators of the results with four projections were better than those of the results with one projection. And NASOLS provided the best results. In addition, ASOLS and OLS algorithms can only reconstruct one target, then there is no result in Fig. 13(b).

Fig. 13 The quantitative analysis of phantom experiments (two targets) with different projections (1, 2, 3, 4) for five methods: (a) and (b) LEs, (c) NRMSE, and (d) CNR.

3.2.2 In-vivo experiments

To further evaluate the performance of the proposed algorithm, implanted in-vivo experiments were performed on two small mice. All animal studies were performed in accordance with the Fourth Military Medical University Guide for the Care and Use of Laboratory Animals formulated by the National Society for Medical Research. To relieve the pain of the mouse, the experiment was conducted under isoflurane gas anesthesia. In addition, Cy5.5 solution with concentration of 4000 nM was injected into the glass tube with a diameter of 2.1 mm and a height of 2.8 mm. The glass tube was implanted into the adult mouse (excitation spectrum at 671 nm, emission spectrum at 710 nm). First, we collected the fluorescence images and then the CT data. After segmentation, we can obtain the main organs of mouse, including heart, lung, liver, kidney, and muscle, etc.

Figure 14 shows the recovered results on 3D views and 2D views combined with CT for two mice. For two mice, the real center locations of the targets are at (19.8, 27.1, 8.1 mm) and (14.6, 19.4, 7.1 mm), respectively. The inverse mesh of Mouse 1 has 3878 nodes and 18,866 tetrahedral elements while the mesh of Mouse 2 has 9466 nodes and 47,631 tetrahedral elements. Table 4 shows the quantitative results for two mice by NASOLS, ASOLS, OLS, CoSaMP and gOMP, respectively. In the 3D views, the real target and reconstructed target were depicted with red and blue regions respectively. It is clear that OLS and gOMP algorithms has provided a deviated result with great error, whereas NSAOLS can locate the target center more accurately according to the 3D views and transverse views for both two mice. Similar conclusions can also be drawn from LE, NRMSE, CNR of Table 4. It can be concluded that NASOLS shows good advantages compared with the other four methods.

Fig. 14 In-vivo experiment results on small animals. (a)–(e) Recovered results for Mouse 1 by NASOLS, ASOLS, OLS, CoSaMP, and gOMP, respectively. (g)–(j) Reconstruction results for Mouse 2 by NASOLS, ASOLS, OLS, CoSaMP and gOMP, respectively. In the 3D views, the red region represents the real target and the blue region is the recovered results.

Table 4 Quantitative results of in-vivo experiment.

 	Methods	Center (mm)	LE (mm)	NRMSE	CNR	
Mouse 1	NASOLS	(19.6,27.4,6.9)	1.2	0.0013	0.49	
ASOLS	(19.3,26.2,5.8)	2.4	0.0015	0.31	
OLS	(16.9,23.6,6.7)	4.7	0.0024	0.29	
CoSaMP	(19.6,27.1,5.9)	2.1	0.0017	0.19	
gOMP	(18.8,25.4,5.4)	3.3	0.0020	0.20	
Mouse 2	NASOLS	(14.0,20.7,7.3)	1.4	0.0010	1.02	
ASOLS	(13.1,20.9,7.5)	2.1	0.0012	0.66	
OLS	(12.1,17.5,5.7)	3.4	0.0013	0.47	
CoSaMP	(13.5,21.4,7.4)	2.2	0.0010	0.87	
gOMP	(13.5,21.6,7.5)	2.7	0.0015	0.58	
Note: The best results are in bold.

4 Discussion and Conclusion

In this paper, we proposed an adaptive sparsity with neighbor strategy orthogonal least-squares to solve L0-norm regularization problem in FMT. An adaptive strategy was proposed to enhance the practicability of our algorithm. A nonlinear function was presented to adjust the sparsity step, so the estimated sparsity can quickly approach the real sparsity. Then, the nearest neighbor strategy based on a finite-element tetrahedral mesh structure was used to form the support set, which has enhanced the local connection characteristic. Actually, the optimal solution based on FEM was usually distributed around the maximum energy node. The nearest neighbor idea was taking advantage of this property, so the developed method could provide more accurate recovered results, especially for the double-target reconstruction. From the double-target simulation experiments, NASOLS could clearly distinguish two targets with EED at 2, 3, and 4 mm. In fact, FMT is heavily influenced by light scattering and it is difficult to achieve such resolutions in the real biological applications. Here, the double-target reconstructions were carried out in a completely ideal environment. The external noise, the error caused by the real experimental measurement and so on were not considered here. In the real experiments, all these external factors will also interfere with the reconstruction results and it is hard to achieve such resolutions. Different projection numbers in the single-target and double-target simulation. It demonstrated that the experimental results had become better with the increasing of projection numbers. However, the results became worse when the number of projections reached a certain number. For example, the recovered results with 16 projections were not better than those of eight projections in the simulations according to the evaluating indicators. The reason might be that the redundant information appeared with the increasing of projection numbers, which in turn affected the reconstruction accuracy. In addition, different levels of noise were added to the measurements and the results showed that as the levels of noise increased, the recovered results has become worse according to the evaluating indicators.

The accuracy of photon propagation model could be improved by combining anatomical information acquired from CT or magnetic resonance imaging into FMT reconstructions. In this work, the measurement data were detected by a FMT/micro-CT dual imaging system. The anatomical information were provided by micro-CT. Some researchers have introduced diffusion optical tomography (DOT) into FMT to provide optical parameters, which further improved the quality of FMT.37 However, the inverse problem of DOT is also severely ill-conditioned and ill-posed, and its reconstruction results are still susceptible to noise.

Apart from the traditional iteration-based regularization methods, deep learning has become one of the fastest-growing breakthrough technologies in recent years. There are generally two types of deep learning methods for FMT, namely end-to-end deep neural network and post processing methods. Traditional methods can provide the mathematical model from physical theory while deep learning is difficult to give a theoretical explanation for the FMT forward model. In addition, the design of network architecture and training schemes determine the performance of deep learning-based methods. Even so, deep learning based methods still has unprecedented advantages in FMT.

In summary, NASOLS algorithm improves the reconstruction performance of FMT from the experimental results. Therefore, this new method can promote the practical application of FMT in clinic. However, there are still some deficiencies. First, the selection of nearest neighbor nodes might be repeated in the reconstruction process, which resulted in time consuming of reconstruction process. Second, the NASOLS algorithm could improve the target location but its ability to reconstruct the shape of target was poor. As we all know, the shape reconstruction of FMT is very important and it is still a great challenge in practical applications. We will continue to work towards this shape reconstruction of target in the future.

Acknowledgments

The authors would like to thank the anonymous reviewers for their helpful suggestions. This research was supported by the National Natural Science Foundation of China (Grant Nos. 61906154, 61971350, 61901374, 11871321, 12271434, 66201459, and 62271394); Key Research and Development Projects of Shaanxi Province (Grant No. 2020SF-036); the Youth Innovation Team of Shaanxi Provincial Department of Education (Grant No. 21JP123).

Huangjian Yi received her BS degree in applied mathematics and her PhD in intelligent information processing from Xidian University, Xi’an, China, in 2008 and 2013, respectively. She is currently an associate professor in the School of Information Sciences and Technology, Northwest University, Xi’an, China. Her research interests include fluorescence molecular tomography and medical image processing.

Sihao Ma received his BS degree from Xidian University, Xi’an, China, in 2020. He is a master’s student in the School of Information Sciences and Technology, Northwest University, Xi’an, China. His current research interests include tomography and image processing.

Ruigang Yang received his BS degree from Northwest University, Xi’an, China, in 2020. He is a master’s student in the School of Information Sciences and Technology, Northwest University. His current research interest is optical imaging and image processing.

Sheng Zhong received his BS degree from Xi’an Shiyou University in 2020. Currently, he is pursuing a MS degree in Electronic Information at Northwest University, Xi’an, China. His research interest is fluorescence molecular tomography.

Hongbo Guo received his PhD from Northwest University, Xi’an, China, in 2018. Currently, he is an associate professor in the School of Information Sciences and Technology, Northwest University, Xi’an, China. His research interests include optical molecular tomography, medical image processing and visualization.

Xuelei He received his PhD from Northwest University, Xi’an, China, in 2021. He is currently a lecturer in the School of Information Sciences and Technology, Northwest University, Xi’an, China. His research interests include reconstruction of optical molecular image, and biomedical images analysis.

Xiaowei He received his MS degree from Xi’an Jiaotong University in 2005 and his PhD from Xidian University in 2011. He is a professor in the School of Information Sciences and Technology, Northwest University, Xi’an, China. His current research interest is optical molecular imaging.

Yuqing Hou is a professor in the School of Information Sciences and Technology, Northwest University, Xi’an, China. Her current research interest is optical molecular imaging.

Disclosures

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
==== Refs
References

1. Ntziachristos V. et al. , “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. 8 (7 ), 757–761 (2002).1078-8956 10.1038/nm729 12091907
2. Kumar S. Richards-Kortum R. , “Optical molecular imaging agents for cancer diagnostics and therapeutics,” Nanomedicine 1 (1 ), 23–30 (2006).1743-5889 10.2217/17435889.1.1.23 17716206
3. Koch M. Symvoulidis P. Ntziachristos V. , “Tackling standardization in fluorescence molecular imaging,” Nat. Photonics 12 (9 ), 505–515 (2018).NPAHBY 1749-4885 10.1038/s41566-018-0221-5
4. Zhang P. et al. , “A review on imaging methodology advances in fluorescence molecular tomography,” Phys. Med. Biol. 67 , 10TR01 (2022).
5. Corlu A. et al. , “Three-dimensional in vivo fluorescence diffuse optical tomography of breast cancer in humans,” Opt. Express 15 (11 ), 6696–6716 (2007).OPEXFF 1094-4087 10.1364/OE.15.006696 19546980
6. Darne C. Lu Y. Sevick-Muraca E. M. , “Small animal fluorescence and bioluminescence tomography: a review of approaches, algorithms and technology update,” Phys. Med. Biol. 59 (1 ), R1 (2013).PHMBA7 0031-9155 10.1088/0031-9155/59/1/R1 24334634
7. Wang L. V. Wu H. , Biomedical Optics: Principles and Imaging, John Wiley and Sons (2012).
8. Cao X. et al. , “An adaptive Tikhonov regularization method for fluorescence molecular tomography,” Med. Biol. Eng. Comput. 51 (8 ), 849–858 (2013).MBECDY 0140-0118 10.1007/s11517-013-1054-5 23504309
9. Okawa S. et al. , “Reconstruction of localized fluorescent target from multi-view continuous-wave surface images of small animal with LP sparsity regularization,” Biomed. Opt. Express 5 (6 ), 1839–1860 (2014).BOEICL 2156-7085 10.1364/BOE.5.001839 24940544
10. Dutta J. et al. , “Joint L1 and total variation regularization for fluorescence molecular tomography,” Phys. Med. Biol. 57 (6 ), 1459 (2012).PHMBA7 0031-9155 10.1088/0031-9155/57/6/1459 22390906
11. Guo H. et al. , “Total variation constrained graph manifold learning strategy for Cerenkov luminescence tomography,” Opt. Express 30 (2 ), 1422–1441 (2022).OPEXFF 1094-4087 10.1364/OE.448250 35209303
12. Mohajerani P. Ntziachristos V. , “An inversion scheme for hybrid fluorescence molecular tomography using a fuzzy inference system,” IEEE Trans. Med. Imaging 35 (2 ), 381–390 (2015).ITMID4 0278-0062 10.1109/TMI.2015.2475356 26340771
13. Yi H. et al. , “Reconstruction algorithms based on L1-norm and L2-norm for two imaging models of fluorescence molecular tomography: a comparative study,” J. Biomed. Opt. 18 (5 ), 056013 (2013).JBOPFO 1083-3668 10.1117/1.JBO.18.5.056013
14. Fiks I. I. Turchin I. V. , “Reconstruction of fluorophore concentration distribution in diffuse fluorescence tomography based on Tikhonov regularisation and nonnegativity condition,” Quantum Electron. 51 (5 ), 400 (2021).QUELEZ 1063-7818 10.1070/QEL17560
15. Lampe B. et al. , “Restricted entropy and spectrum properties for the compressively sensed domain in hyperspectral imaging,” IEEE Trans. Geosci. Remote Sens. 58 (8 ), 5642–5652 (2020).IGRSD2 0196-2892 10.1109/TGRS.2020.2968077
16. Ranjan S. Vidyasagar M. , “Rank restricted isometry property implies the rank robust null space property in compressed sensing for matrix recovery,” in Sixth Indian Control Conf. (ICC), IEEE, pp. 224–225 (2019).10.1109/ICC47138.2019.9123237
17. Wang D. et al. , “Model recovery for Hammerstein systems using the auxiliary model based orthogonal matching pursuit method,” Appl. Math. Model. 54 , 537–550 (2018).AMMODL 0307-904X 10.1016/j.apm.2017.10.005
18. Zhou H. Narayanan R. M. , “A dual-mesh microwave reconstruction method based on compressive sampling matching pursuit algorithm,” Progr. Electromagn. Res. 166 , 43–57 (2019).PELREX 1043-626X 10.2528/PIER19090203
19. Song S. et al. , “An improved reconstruction method based on auto-adjustable step size sparsity adaptive matching pursuit and adaptive modular dictionary update for acoustic emission signals of rails,” Measurement 189 , 110650 (2022).0263-2241 10.1016/j.measurement.2021.110650
20. Yang F. et al. , “High-resolution mesoscopic fluorescence molecular tomography based on compressive sensing,” IEEE Trans. Biomed. Eng. 62 (1 ), 248–255 (2014).IEBEAX 0018-9294 10.1109/TBME.2014.2347284 25137718
21. Kong L. et al. , “Reconstruction for fluorescence molecular tomography via adaptive group orthogonal matching pursuit,” IEEE Trans. Biomed. Eng. 67 (9 ), 2518–2529 (2020).IEBEAX 0018-9294 10.1109/TBME.2019.2963815 31905129
22. Zhang S. et al. , “Robust reconstruction of fluorescence molecular tomography based on sparsity adaptive correntropy matching pursuit method for stem cell distribution,” IEEE Trans. Med. Imaging 37 (10 ), 2176–2184 (2018).ITMID4 0278-0062 10.1109/TMI.2018.2825102 29993826
23. Ye J. et al. , “Fast and robust reconstruction for fluorescence molecular tomography via a sparsity adaptive subspace pursuit method,” Biomed. Opt. Express 5 (2 ), 387–406 (2014).BOEICL 2156-7085 10.1364/BOE.5.000387 24575335
24. Soussen C. et al. , “Joint k-step analysis of orthogonal matching pursuit and orthogonal least squares,” IEEE Trans. Inf. Theory 59 (5 ), 3158–3174 (2013).IETTAW 0018-9448 10.1109/TIT.2013.2238606
25. Bi B. et al. , “Image reconstruction for diffuse optical tomography based on radiative transfer equation,” Computat. Math. Methods Med. 2015 , 286161 (2015).10.1155/2015/286161
26. Arridge S. R. Schotland J. C. , “Optical tomography: forward and inverse problems,” Inverse Probl. 25 (12 ), 123010 (2009).INPEEY 0266-5611 10.1088/0266-5611/25/12/123010
27. Zhou L. Yazici B. , “Discretization error analysis and adaptive meshing algorithms for fluorescence diffuse optical tomography in the presence of measurement noise,” IEEE Trans. Image Process. 20 (4 ), 1094–1111 (2010).IIPRE4 1057-7149 10.1109/TIP.2010.2083677 20923735
28. Cong A. Wang G. , “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express 13 (24 ), 9847–9857 (2005).OPEXFF 1094-4087 10.1364/OPEX.13.009847 19503194
29. Song X. et al. , “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15 (26 ), 18300–18317 (2007).OPEXFF 1094-4087 10.1364/OE.15.018300 19551128
30. Elad M. , Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing, 1st ed., Springer-Verlag, New York, New York, United States (2010).
31. He X. et al. , “Nonconvex Laplacian manifold joint method for morphological reconstruction of fluorescence molecular tomography,” Mol. Imaging Biol. 23 (3 ), 394–406 (2021).10.1007/s11307-020-01568-8 33415678
32. Yi H. et al. , “Reconstruction for limited-projection fluorescence molecular tomography based on a double-mesh strategy,” BioMed Res. Int. 2016 , 5682851 (2016).10.1155/2016/5682851 27830148
33. Dogdas B. et al. , “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52 (3 ), 577–587 (2007).PHMBA7 0031-9155 10.1088/0031-9155/52/3/003 17228106
34. Alexandrakis G. Rannou F. R. Chatziioannou A. F. , “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50 (17 ), 4225 (2005).PHMBA7 0031-9155 10.1088/0031-9155/50/17/021 16177541
35. Chen X. et al. , “3D reconstruction of light flux distribution on arbitrary surfaces from 2D multi-photographic images,” Opt. Express 18 (19 ), 19876–19893 (2010).OPEXFF 1094-4087 10.1364/OE.18.019876 20940879
36. He X. et al. , “Half thresholding pursuit algorithm for fluorescence molecular tomography,” IEEE Trans. Biomed. Eng. 66 (5 ), 1468–1476 (2018).IEBEAX 0018-9294 10.1109/TBME.2018.2874699 30296209
37. Tan Y. Jiang H. , “DOT guided fluorescence molecular tomography of arbitrarily shaped objects,” Med. Phys. 35 (12 ), 5703–5707 (2008).MPHYA6 0094-2405 10.1118/1.3020594 19175127