Improving image quality in diffuse optical tomography

Diffuse optical tomography (DOT) modality uses diffusion equation solution of low power laser distribution in the imaging tissue. DOT method has three main running mode and three geometric subbranches. Run modes are: Continuous Wave (CW), Time Resolved (TR) and Frequency Domain (FD) modes. These run modes might have transmission-through, back-reflected or ring geometry depend on the source and detector placements on tissue surface. In this work, we tested our novel imaging method on the frequency-domain back-reflected DOT simulation model. Our novelty is: selecting a great number of multi-frequency data in the wide frequency spectrum which was not tested before in the literature. Most of the literature works consist of narrow frequency band with the limited number of frequencies such as maximum 12 multi-frequency-data which cover from 100 up to 700 MHz. In our work, we tested and compared reconstructed inclusion images generated by our 500 wide spectrum multi-frequency data versus reconstructed inclusion images generated by 20 narrow band multi-frequency data which was generally used in FD DOT studies. We observed that our novel imaging method which used wide spectrum multi-frequency data is superior to the traditionally used narrow-band multi-frequency method based on the reconstructed inclusion image localization errors. Since most of the homogenous geometry consists of fat tissue, we selected background absorption coefficient µa = 0.2 cm−1 and tissue scattering coefficient μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\upmu }$$\end{document}s = 80 cm−1 for 800 nm laser source wavelength. Three-dimensional cubic imaging simulation media has 75-mm grid sizes in each direction. Two different imaging scenarios were tested. In the first scenario, one inclusion with absorption coefficient µa = 0.7 cm−1 was put inside the three-dimensional imaging geometry at 27 mm depth. In the second scenario one inclusion with absorption coefficient µa = 0.7 cm−1 was embedded in 48 mm depth location. Inclusions are cubic and they have 6 mm xyz length at each direction.


Introduction
We demonstrated in this work that selection of multi-frequency data for Frequency-Domain Diffuse Optical Tomography (FDDOT) is the important step for image reconstruction level. In literature, narrow band multi-frequency method with limited number of data was used such as between 100 and 700 MHz frequencies with maximum 12 data points up to now. In our work, we selected wide frequency spectrum which has 500 multi-frequency data. We tested our wide spectrum data model versus narrow-band 20 multi-frequency data model. Our imaging method covers selecting a great number of multi-frequency data in the wide frequency spectrum which was not tested before in the literature. Most of the literature works consist of narrow frequency band with the limited number of frequencies such as maximum 12 multi-frequency-data which cover from 100 up to 700 MHz. In our work, we tested and compared reconstructed inclusion images generated by our 500 wide spectrum multi-frequency data versus reconstructed inclusion images generated by 20 multi-frequency data which was generally used in FDDOT studies. We observed that our novel imaging method which used wide frequency spectrum data is superior to the traditionally used narrow-band multi-frequency method based on the reconstructed inclusion image localization errors. FDDOT imaging method was studied in literature (Doulgerakis et al. 2019;Brian et al. 1997;Pogue et al. 1997;Ko et al. 2007;Ren et al. 2006;Kazanci 2021). Multi-frequency technique was first mentioned and used early in optical diffusion tomography (Milstein et al. 2004). In that work, basically three different modulation frequencies were used, modeled in simulations and tested with experimental method. They used 78.4 MHz, 314 MHz and mix of these frequencies. In another work, 7 different modulation frequencies were used and compared with each other (Unlu et al. 2006). They used 110 MHz up to 280 MHz modulation frequencies with 30 MHz increasing steps. In that work, they used different frequency selection configurations, such as 2 configurations with 110 MHz and 260 MHz, 110 MHz and 280 MHz or 3 multi-frequency configurations. They also used 6 different modulation frequencies. According to the authors, by selecting the most appropriate modulation frequencies combination, they improved the reconstructed inclusion image location and contrast accuracy. Especially when they used 6 frequencies data configuration, their reconstructed inclusion image location and absorption coefficient errors were significantly reduced. But on the other hand, according to the authors, background scattering and absorption coefficients were degraded and eventually contrast to noise ratio (CNR) was decreased. In another study, different coupling of 5 modulation frequencies were compared with each other (Shifa et al. 2021). They used 200 MHz up to 700 MHz modulation frequencies. They tested with different number of groups such as 500 MHz alone, 500 MHz + 600 MHz, 500 MHz + 600 MHz + 700 MHz, 400 MHz + 5 00 MHz + 600 MHz + 700 MHz, and 200 MHz + 400 MHz + 500 MHz + 600 MHz + 700 MHz. They compared all configurations based on the error and noise analysis. According to their comparison, most successful method is: Using all frequencies data. In another simulation study, researchers used 100 MHz, 150 MHz, 200 MHz, and 250 MHz modulation frequencies and their different number of configurations (Chen et al. 2015). They tested their novel inverse problem solution image reconstruction algorithm with multi-frequency data. In one of the earlier work, different number of frequencies were conducted which are starting from 20 upto 500 MHz (Intes and Chance 2005). Different frequency configurations were used and compared in that study. According to the authors, if 7 or more high frequencies were combined and used, better image quality and enhancements were obtained. Multi frequency works and studies in literature is covering maximum 12 different modulation frequencies and their different number of combinations. They usually used the 100 MHz modulation core frequency and around of it up to 700 MHz. We have seen that by increasing the number of modulation frequencies in forward model is improving the image quality, reducing reconstructed inclusion image location and contrast errors. On the other hand, by increasing the frequency numbers, tissue background scattering and absorption coefficients noise and artefacts are also increasing, thus effecting the overall image reconstruction process. But we should know that all the works in literature were limited to the exact number of frequencies such as 5 or 6 definite frequencies, these include 100 MHz upto 700 MHz, and they are not covering GHz frequency range. In our work, we compared two different cases which have 20 and 500 multi-frequency data. In the first case, 110 MHz upto 300 MHz with 10 MHz increasing modulation frequency steps were added to the forward model weight functions. In the second scenario, 500 data which starts from 100 MHz up to 5 GHz frequency with 10 MHz increasing step were added to the forward model. Instead of using 5-7 frequencies which are around 100 MHz, we used 500 wide spectrum multi-frequency data. We finally realized that, by increasing the frequency spectrum, reconstructed inclusion image quality was improved, inclusion location error was reduced, and background absorption coefficient noises were eliminated. Up to now in literature, researchers could not eliminate background artefacts and noises by increasing the frequencies to the limited numbers such as 12 different frequencies. On the contrary, we showed that by using wide frequency spectrum, we eliminated the tissue background absorption artefacts and noises. By comparison of two different multi-frequency cases, we demonstrated that 500 multi-frequency combination is superior to the 20 multi-frequency combination, and it is worthy to use wide multi-frequency spectrum.

Method
Three-dimensional tissue simulation imaging phantom was prepared. We placed 8 lasers and 8 photo-detectors on the tissue surface with the pencil beam setup. Three-dimensional imaging geometry has 25 × 25 × 25 voxels and its dimensions are 75 × 75 × 75 mm. 8 lasers and 8 photo-detectors are generating 64 laser-photodetector matches on the tissue surface. Two separate multi-frequency data cases for two different scenarios were compared with each other. 800 nm laser source wavelength was selected. Background tissue absorption coefficient µ a is 0.2 cm −1 and tissue scattering coefficient μ s is 80 cm −1 . In the first scenario, we put 1 inclusion inside the homogenous tissue background at the 27 mm depth. Center coordinate of inclusion: (z, y, x) = (27, 0, 0) mm. Inclusion has 6 mm x × y × z grid dimensions and it is cubic. In the second scenario, we put one inclusion at the 48 mm depth which has 6 mm x × y × z grid dimensions. Center coordinates of inclusion: (z, y, x) = (48, 0, 0) mm. Absorption coefficient of inclusion was selected as µ a = 0.7 cm −1 . We used 20 multi-frequency data (upto 300 MHz) for case 1 and 500 multi-frequency data (upto 5 GHz) for case 2. Threedimensional tissue simulation geometry has 25 × 25 × 25 grid elements. Back-reflection imaging geometry has total of 15,625 voxel elements. We have 1280 × 15,625 (64 source-detector match × 20 multi-frequency data = 1280 (64 × 20) × 15,625 (25 × 25 × 25) voxel elements. Each source-detector match in forward model has 15,625 voxel components, so 64 × 20 = 1280 source-detector match with 20 frequencies has 1280 × 15,625 elements.) matrix elements for case 1, and 32,000 × 15,625 (64 source-detector match × 500 multi-frequency data = 32,000 (64 × 500) × 15,625 (25 × 25 × 25) voxel elements for case 2. Each source-detector match in the forward model has 15,625 voxel components, so 64 × 500 = 32,000 source-detector match with 500 frequencies has 32,000 × 15,625 elements.) matrix elements for case 2. Tissue background absorption coefficient µ a = 0.2 cm −1 , and scattering coefficient µ s = 80 cm −1 . Imaging coordinate system has 75 mm length at each direction. 8 lasers and 8 photo-detectors were grouped as illustrated in Fig. 1. Sources are blue diamonds and detectors are yellow squares. Lasers and photodetectors were placed as on chessboard. Forward model weight function was given in Eq. 1 (Wang and Wu 2007) and Eq. 2 as A matrix. A matrix consists of m modulation frequency index, distance SV each source to each voxel distance, distance VD each voxel to each detector distance, c light speed inside the imaging tissue media, D diffusion coefficient. Forward model is the matrix multiplication of A weight matrix and Δ a perturbation differences over homogeneous tissue background. Forward model was represented in Eq. 2. Δ a perturbation vector has 25 × 25 × 25 = 15,625 vector elements which its dimension is (15,625,1). Perturbation vector b has 1280 group elements for case 1 which has 20 multi-frequency elements and 32,000 group elements for case 2 which has 500 multi-frequency elements. As a result, b perturbation vector has 1280 × 15,625 × 15,625 × 1 = 1280 vector elements for case 1, b perturbation vector has 32,000 × 15,625 × 15,625 × 1 = 32,000 vector elements for case 2. To solve the unknown Δ a absorption coefficients, the forward model should be solved by applying inverse problem solution algorithm. Since we are only testing our novel multifrequency method in the homogeneous simulation phantom, we do not need to apply ill-posed regularization or iterative algebraic reconstruction methods while solving the inverse problem. Ill-posed problems are only existing in real world experiments which has heterogenous tissue conditions and boundary problems. It can also be accepted that breast tissue model is mostly homogenous problem since it has more fat tissue percentage. For these reasons we only applied the pseudoinverse problem solution algorithm which was shown in Eq. 3.

Results
Two different scenarios for two different cases were compared with each other. In the first scenario and first case, 20 multi-frequency data was used for one inclusion scenario. One inclusion was put inside the three-dimensional geometry as illustrated in Fig. 4A for XZ and in Fig. 4B for XY view. Two reconstructed inclusion images were drawn for XZ view in Fig. 4C and XY view in Fig. 4D for the first case. Center coordinates of one inclusion: (z, y, x) = (27, 0, 0) mm. It can roughly be seen that reconstructed inclusion images are in correct depth location. Original inclusion position is at 27 mm depth, and it has 6 mm grid sizes on each direction. Reconstructed inclusion image is also at 27 mm depth, only xy sizes are wider than original inclusion xy sizes. Reconstructed inclusion image is bigger than original inclusion. Detailed localization errors were given in Table 1 based on the Eq. 4.
In the first scenario and second case, 500 multi-frequency data was used. One inclusion was put inside the three-dimensional geometry as illustrated in Fig. 5A for XZ and in Fig. 5B for XY view. One embedded perturber was reconstructed and shown for XZ view in Fig. 5C and XY view in Fig. 5D for the second case. Reconstructed inclusion image is in the correct depth location. Original inclusion is at the 27 mm depth, and reconstructed inclusion image is also at the same 27 mm depth position. The only difference is: Reconstructed inclusion is much likely oval. Detailed error analysis was done by using Eq. 4. Detailed position and absorption coefficient errors were given by using Eqs. 4, and 5 in Table 1.
In the second scenario and first case, 20 multi-frequency data was used for one inclusion scenario. One absorber was embedded inside the three-dimensional geometry as shown in   Fig. 6B for XY view. Center coordinates of one inclusion: (z, y, x) = (48, 0, 0) mm. One inclusion image was reconstructed and illustrated for XZ view in Fig. 6C and XY view in Fig. 6D for the first case. It can roughly be seen that reconstructed inclusion image is in wrong depth position (38 mm). Original inclusion position is at the 48 mm depth, and it has 6 mm grid sizes in each direction. On the contrary, the reconstructed inclusion image is at 38 mm depth. Detailed localization errors were given based on the Eq. 4 in Table 1. In the second scenario and second case, 500 multi-frequency data was used. One absorber was put inside the three-dimensional geometry as shown in Fig. 7A for XZ and in Fig. 7B for XY view. One inclusion was reconstructed and illustrated for XZ view in Fig. 7C and XY view in Fig. 7D for the second case. Reconstructed inclusion is at the correct depth location. Original inclusion is at the 48 mm depth, and reconstructed inclusion is at the same 48 mm depth position. Detailed position error analysis was done by using Eq. 4. Detailed position error and absorption coefficient concentration errors were given by using Eq. 4, and Eq. 5 in Table 1. Position errors were calculated by selecting region of interest areas. 75 mm × 75 mm × 75 mm cartesian coordinate system was given in Fig. 8. Since we were trying to reconstruct the inclusions in deep tissue, we put nearest neighbor source and detector to long-distance. Nearest neighborhood source-detector distance is 20 mm. In literature, according to most of the back-reflection DOT works, nearest neighbor source-detector distance was generally selected as 3 mm. But in their works, they could not reconstruct the inclusion images in deep tissue. Reconstructed inclusion images were becoming superficial where depth information was lost.

Discussion
Two different number of multi-frequency data was added to the forward model for two different imaging scenarios. In the first scenario, we put one inclusion inside the imaging geometry at 27 mm depth. For case 1, only 20 multi-frequency data between 110 and 500 MHz was used to construct the forward model weight matrix coefficients. For case 2, 500 multi-frequency data between 110 MHz and 5 GHz (wide frequency spectrum) was used to construct the forward model weight matrix coefficients. When we applied the pseudo-inverse problem solution algorithm, we observed that reconstructed inclusion images for case 1 (which has 20 multi-frequency data) has % 12 position error (real inclusion objects are at 27 mm depth with 6 mm grid dimensions on each direction, reconstructed inclusion objects are at 27 mm depth with bigger dimensions, we used Eq. 4 for localization error). Reconstructed inclusion image for case 2 (which has 500 multi-frequency data) have same % 12 position error (real inclusion objects are at 27 mm depth with 6 mm grid dimensions on each direction, reconstructed inclusion objects are at the same 27 mm depth but with bigger dimensions. We used Eq. 4 for localization error). Similarly, we tested our method with one inclusion model with deeper location in the second scenario. Reconstructed inclusion image for case 1 (which has 20 multi-frequency data) has % 100 position error (real inclusion object is at 48 mm depth with 6 mm grid dimensions on each direction, reconstructed inclusion object is at 38 mm depth with bigger dimensions, we used Eq. 4 for localization error). On the other hand, reconstructed inclusion image for case 2 (which has 500 multi-frequency data) has only % 5 position error (real inclusion object is at 48 mm depth with 6 mm grid dimensions on each direction, reconstructed inclusion object is at the same 48 mm depth but with bigger dimensions. We used Eq. 4 for localization error). In literature works, researchers have tried to improve their image quality by only using maximum 12 multi-frequency data in the narrow multi-frequency band, such as between 100 and 700 MHz. In our novel model we increased the frequency upto 5 GHz, and we also increased the multi-frequency data numbers to the 500. We tested our model versus narrow frequency band data which its frequencies cover between 110 and 300 MHz, and it has 20 multi-frequency data. We observed that by adding more multi-frequency data, it improved the reconstructed image quality and decreased the position error, and background noises and artefacts. This work demonstrated that adding more multi-frequency data in the wide frequency spectrum is important step to improve the reconstructed image quality.