Greedy algorithms for diffuse optical tomography reconstruction

Optics Communications 410 (2018) 164–173

Contents lists available at ScienceDirect

Optics Communications journal homepage: www.elsevier.com/locate/optcom

Greedy algorithms for diffuse optical tomography reconstruction B.P.V. Dileep *, Tapan Das, Pranab K. Dutta Department of Electrical Eng., IIT Kharagpur, Kharagpur, 721302, India

a r t i c l e

i n f o

Keywords: Diffuse optical tomography (DOT) Single measurement vector (SMV) Multiple measurement vectors (MMV) Greedy algorithms

a b s t r a c t Diffuse optical tomography (DOT) is a noninvasive imaging modality that reconstructs the optical parameters of a highly scattering medium. However, the inverse problem of DOT is ill-posed and highly nonlinear due to the zig-zag propagation of photons that diffuses through the cross section of tissue. The conventional DOT imaging methods iteratively compute the solution of forward diffusion equation solver which makes the problem computationally expensive. Also, these methods fail when the geometry is complex. Recently, the theory of compressive sensing (CS) has received considerable attention because of its efficient use in biomedical imaging applications. The objective of this paper is to solve a given DOT inverse problem by using compressive sensing framework and various Greedy algorithms such as orthogonal matching pursuit (OMP), compressive sampling matching pursuit (CoSaMP), and stagewise orthogonal matching pursuit (StOMP), regularized orthogonal matching pursuit (ROMP) and simultaneous orthogonal matching pursuit (S-OMP) have been studied to reconstruct the change in the absorption parameter i.e, 𝛥𝛼 from the boundary data. Also, the Greedy algorithms have been validated experimentally on a paraffin wax rectangular phantom through a well designed experimental set up. We also have studied the conventional DOT methods like least square method and truncated singular value decomposition (TSVD) for comparison. One of the main features of this work is the usage of less number of source–detector pairs, which can facilitate the use of DOT in routine applications of screening. The performance metrics such as mean square error (MSE), normalized mean square error (NMSE), structural similarity index (SSIM), and peak signal to noise ratio (PSNR) have been used to evaluate the performance of the algorithms mentioned in this paper. Extensive simulation results confirm that CS based DOT reconstruction outperforms the conventional DOT imaging methods in terms of computational efficiency. The main advantage of this study is that the forward diffusion equation solver need not be repeatedly solved. © 2017 Elsevier B.V. All rights reserved.

1. Introduction The diffuse optical tomography (DOT) has emerged as an imaging modality that attracted various research groups in the recent years. In DOT, near infrared (NIR) light is passed through the tissue under study by using one or several sources and the measurements are obtained at the boundary of the tissue with the help of detectors. The main aim of DOT [1] is to reconstruct the optical parameters of a cross section of an illuminated tissue from the boundary data. The optical parameters that are considered in DOT are absorption and scattering coefficients. The primary challenge in DOT is the inverse problem [2] which is ill-posed, nonlinear, and under-determined due to the propagation of photons that diffuses through the scattering medium such as illuminated tissue. In general, the optical parameters [3] varies nonlinearly with the scattered optical flux at the boundary of the tissue.

The photons within the illuminated tissue undergoes multiple scattering as the photons propagates through the medium, and hence, the medium such as tissue is scattering dominant than absorption in nature. The diffusion equation is used as the model for the propagation of photons through the illuminated tissue. Here, the interest is to image the change in the absorption parameter [4] by keeping the scattering parameter as constant throughout the paper. The DOT is used as an important tool [5,6] for many biomedical imaging areas owing to its inherent features like noninvasive, portability and inexpensiveness. DOT has been deployed in potential applications such as breast and brain imaging [7,8]. The recent advances [9,10] in the field of DOT by different research communities was on proposing new methods to deal with the significant difficulties present in the reconstruction problem. The conventional DOT imaging methods to solve this reconstruction problem include iterative [11] and linearization techniques. The

* Corresponding author.

E-mail addresses: [email protected] (B.P.V. Dileep), [email protected] (T. Das), [email protected] (P.K. Dutta). https://doi.org/10.1016/j.optcom.2017.09.056 Received 5 June 2017; Received in revised form 13 September 2017; Accepted 16 September 2017 0030-4018/© 2017 Elsevier B.V. All rights reserved.

B.P.V. Dileep et al.

Optics Communications 410 (2018) 164–173

iterative methods distorted born [12] and Bayesian approaches (based on the born method) [13] are used to deal the nonlinearity that are inherent in the inverse problem. However, these methods fail as the greens function needs to be calculated during each iteration which adds the computational complexity to the problem. On the other hand, several groups have published work on linearization techniques [14,15]. These methods overcome the nonlinearity present in the inverse problem. However, they fail when the contrast between the unknown optical properties and the homogeneous background is more than the born estimated limit [16]. Also, there are other methods proposed in the literature to solve the DOT inverse problem [17]. The recent breakthrough in signal processing has led to new developments and extensive investigations in DOT. The DOT with CS is one among them. The concept of compressive sensing [18] uses significantly less number of measurements for a signal to be recovered from the data. In this scenario, the DOT problem has been converted to single measurement vector (SMV) and multiple measurement vector (MMV) problems [19] using the CS framework. The existing methods to solve the SMV problems in CS are summarized in review article [20], and in this paper, Greedy algorithms like OMP [21], CoSaMP [22], and StOMP [23], ROMP [24], have been used to solve the DOT inverse problem under SMV using CS framework. There are a number of algorithms proposed by different research groups in the literature to solve the MMV problems in CS and few among them are S-OMP [25], P-thresholding [26], and MUSIC [27]. We also have studied S-OMP algorithm [25] in this paper to solve this reconstruction problem under MMV using CS framework. In addition, the conventional DOT methods like least square method and truncated singular value decomposition (TSVD) [1] have been studied. The capability of the Greedy algorithms in this reconstruction has been studied by considering the tissue like paraffin wax rectangular phantom with single inclusion. To validate these, the experimental verification has been performed on a paraffin wax rectangular phantom through a well designed experimental setup. The performance metrics such as mean square error (MSE), normalized mean square error (NMSE), structural similarity index (SSIM) [28], and peak signal to noise ratio (PSNR) [29] are used to assess the performance of the Greedy algorithms in the DOT reconstruction. The advantage of using CS based DOT reconstruction is that the greens function need not be repeatedly solved during each iteration as compared to that of conventional DOT iterative methods. The simulation studies show that CS based DOT reconstruction works well and reliably reconstructs the change in the 2D absorption coefficient map as compared to conventional DOT imaging methods. The remaining sections of this paper are organized as follows. Mathematical preliminaries are given in Section 2. The compressive sensing framework for DOT problem is presented in Section 2.2. Section 3 describes the Greedy algorithms. Section 4 illustrates the experimental studies for DOT. The results and discussion are provided in Section 5 followed by the conclusion in Section 6.

2.2. CS framework for DOT problem We use the concept of imaging physics for imaging geometry in DOT since it provides convenient detection system using the array of massive detectors such as photodetectors. This theory can be extended to frequency and time domain DOT as well. In this section, the formulation of forward and inverse problems of DOT under SMV and MMV problems will be discussed in detail using CS framework [27]. The diffusion equation plays a significant role in these formulations. 2.2.1. Forward problem formulation using diffusion equation Considering the diffused medium as highly scattered with low absorption, the following diffusion equation is used in continuous domain to model the propagation of photons from a continuous wave source 𝑞 (𝐫; 𝑝) for a specific illumination pattern denoted by 𝑝 [14] ∇.𝐷 (𝐫) ∇𝑣 (𝐫, 𝑡 ∶ 𝑝) − 𝛼 (𝐫) 𝑣 (𝐫, 𝑡 ∶ 𝑝) = −𝑞 (𝐫; 𝑝) ,

(1)

Here, 𝑝 = 1, 2, … , 𝑙, the normalized absorption and diffusion coefficients 𝑐 , where 𝜇𝑎 (𝐫) denotes are given by 𝛼 (𝐫) = 𝑐𝜇𝑎 (𝐫) and 𝐷 (𝐫) = 3 𝜇 (𝐫)+𝜇 ′ ( 𝑎 𝑠 (𝐫)) the absorption coefficient, 𝜇𝑠′ (𝐫) is the reduced scattering coefficient, 𝑐 is the speed of light, and 𝐫 is the location inside the medium respectively. In addition, the photon flux from a continuous wave source satisfies the boundary conditions. The extrapolated boundary condition has been used in this formulation and is given by (2)

𝑣 + 𝓁 𝑛.∇𝑣 ̂ =0

where 𝓁 denotes the extrapolation length and is computed by using 𝓁 = 2𝑊̂ 𝐷0 , 𝑊̂ is the coefficient that depends upon the relative refractive index [14], 𝐷0 represents the background diffusion coefficient, and 𝑛̂ indicates the direction of normal to the medium. Imaging the change in the absorption parameter 𝛥𝛼 is considered in this paper because the angiogenesis in the cancer disease changes mainly the absorption parameter in the hemoglobin content and the contrast of this parameter is very important in optical imaging. In general, the scattering coefficient is much greater than the absorption coefficient in the highly scattered medium accordingly 𝜇𝑠′ >> 𝜇𝑎 . Therefore, the diffusion constant can be assumed to be constant such that 𝐷 (𝐫) = 𝐷0 . Similarly, we can express the absorption coefficient as 𝛼 (𝐫) = 𝛼0 + 𝛥𝛼 (𝐫), where 𝛼0 denotes the background absorption coefficient, 𝛥𝛼 is the change in the absorption coefficient. Unlike the iterative DOT imaging methods, we will not assume that the change is small in comparison to the background absorption coefficient 𝛼0 . The optical intensity of scattered diffused optical flux at position 𝐫 for the 𝑝th illumination pattern is given by the following equation [14] ( )2 [ ] 𝓁∗ 𝛾 (𝐫; 𝑝) = 1 + 𝑣0 (𝐫; 𝑝) − 𝑣 (𝐫; 𝑝) , 𝑝 = 1, 2, … , 𝑙 𝓁 (3) ( )2 ( ) ( ) ( ) 𝓁∗ = 1+ 𝑔0 𝐫, 𝐫 ′ 𝑣 𝐫 ′ ; 𝑝 𝛥𝛼 𝐫 ′ ; 𝑝 𝑑𝐫 ′ ∫ 𝓁 𝐷

where 𝓁 ∗ = 3 − 𝑐0 , 𝑘0 is the diffuse wave number and can be computed √𝛼 by 𝑘0 = 𝐷0 . 0 In this problem, we assume that the optical parameters are sparsely distributed in the scattered medium and can be expressed as

2. Mathematical preliminaries and CS framework for DOT problem 2.1. Mathematical preliminaries

𝛥𝛼 (𝐫) =

The mathematical preliminaries for CS framework for DOT problem are considered in this section. The basic definitions of mutual coherence [31] and restricted isometry property (RIP) [32] plays a significant role in understanding the concept of CS to the DOT problem. In order to satisfy the definitions of mutual coherence and RIP property, the measurement matrix 𝑊 must have the linearly independent columns [30]. Throughout this paper, the 𝑚th row and 𝑛th column of the matrix 𝑋 is denoted by 𝐱𝑚 and 𝐱𝑛 , respectively. When the index set is denoted by 𝐻, the corresponding rows of 𝑋 is represented by 𝑋 𝐻 , and the columns of 𝑊 is given by 𝑊 𝐻 . supp𝑋 represents the index set which are nonzero rows in 𝑋. Note that 𝑋 is a vector for SMV and a matrix for MMV problems.

𝑘 ∑

( ) ( ) 𝛥𝛼 𝐫(𝑖) 𝛿 𝐫 − 𝐫(𝑖) , 𝐫 ∈ 𝛺

(4)

𝑖=1

( ) where 𝛿 𝐫 − 𝐫(𝑖) is the delta function, the 2D locations are represented { }𝑘 by 𝐫(𝑖) 𝑖=1 and k is the number of optical parameter targets which is { }𝑛 unknown. Our assumption is that 𝐫(𝑖) 𝑖=1 possible optical parameter { }𝑘 target locations are available on a fixed grid from which the 𝐫(𝑖) 𝑖=1 locations have been selected. The following equation has been used to measure the scattered diffused optical flux from the 𝑝th illumination pattern 𝛽(𝐫; 𝑝) =

𝑘 ∑ 𝑖=1

165

( ) ( ) ( ) 𝑔0 𝐫, 𝐫(𝑖) 𝑣 𝐫(𝑖) ; 𝑝 𝛥𝛼 𝐫(𝑖) , 𝑝 = 1, 2, … , 𝑙

(5)

B.P.V. Dileep et al.

Optics Communications 410 (2018) 164–173

where 𝑣 (𝐫; 𝑝) is the unknown total optical flux that satisfies the following Foldy–Lax equation (details can be found in Appendix).

2.3. Accurate inversion for 𝛥𝛼 Considering that we estimated the active set for support of 𝐶, then the least square solution 𝑋̂ 𝐶 based on 𝐶 can be obtained by using

2.2.2. SMV and MMV formulations for the inverse problem This section describes the DOT imaging problem that has been converted to SMV and MMV problems [19] using the imaging geometry. This geometry is commonly used in most of the DOT imaging problems. We use different illumination patterns for imaging the change in the absorption map. By taking { the }𝑚 measurement of scattered optical flux (5) at detector locations 𝐫𝑑𝑖 , the DOT imaging problem can be expressed in the 𝑖=1 following matrix form

𝑅𝑚×𝑟 ,

and the noise matrix

3. Greedy algorithms (7)

The following Greedy algorithms like OMP [21], CoSaMP [22], StOMP [23], ROMP [24] and S-OMP [25] have been used to study the DOT problem in this paper. These algorithms solve the SMV and MMV based formulated DOT problems using CS framework.

And the measurement matrix 𝑊 ∈ 𝑅𝑚×𝑛 is given by ( ) ( ) ⎛ 𝑔0 𝐫𝑑 ; 𝐫1 … 𝑔0 𝐫𝑑 ; 𝐫𝑛 ⎞ 1 1 ⎜ ⎟ 𝑊 =⎜ ( ⋮ (8) ) ⋱ ( ⋮ )⎟ ⎜𝑔 𝐫 ; 𝐫 ⎟ ⋯ 𝑔 𝐫 ; 𝐫 ⎝ 0 𝑑𝑚 1 0 𝑑𝑚 𝑛 ⎠ ( ) where the general expression for 𝑔0 𝐫, 𝐫 ′ for a homogeneous medium is written as ( ) ( ) exp −𝑘0 ||𝐫 − 𝐫 ′ || (9) 𝑔0 𝐫, 𝐫 ′ = 4𝜋𝐷0 |𝐫 − 𝐫 ′ |

3.1. Orthogonal matching pursuit (OMP) algorithm The main idea behind the OMP algorithm [21] is to pick columns of measurement matrix 𝑊 in the Greedy fashion. At each step, we choose the column of 𝑊 so that this is strongly correlated with the remaining columns of the data vector 𝑀. Then the contribution of 𝑊 is subtracted from 𝑀 and subsequently the residual is calculated. The least square step in OMP adds the computational complexity to the algorithm. The main advantage of this algorithm is that it works better for the low sparse signals. Also, it can be implemented very easily. However, the disadvantage of OMP is that it does not perform well for the high sparsity signals. The implementation steps of this algorithm can be found in Appendix.

where ||𝐫 − 𝐫 ′ || denotes the distance between 𝐫 and 𝐫 ′ , 𝐫 and 𝐫 ′ are locations of the targets at different points on the surface of the medium. Here, 𝐫 = 𝐫𝑑𝑖 and 𝐫 ′ = 𝐫𝑖 . In general, however, greens function depends on the particular boundary condition. Proper choice of the homogeneous greens function may produce a scaled version of the original image with some amount of artifacts. The induced optical current matrix 𝑋 as in (6) is given by ( ) ( ) ⎛ 𝑡 𝐫1 ; 1 … 𝑡 𝐫1 ; 𝑙 ⎞ ⎜ 𝑋= (10) ⋮ ⋱ ⋮ ⎟ ) ( )⎟ ⎜ ( ⎝ 𝑡 𝐫𝑛 ; 1 ⋯ 𝑡 𝐫𝑛 ; 𝑙 ⎠ { }𝑘 The rows of 𝑋 are nonzero when 𝐫𝑖 ∈ 𝐫(𝑖) 𝑖=1 . Therefore, we define the sparsity as 𝑘 = ‖𝑋‖0 = supp𝑋 (support of 𝑋), where ‖𝑋‖0 denotes the number of nonzero elements (targets) for SMV and the number of nonzero rows for MMV in 𝑋. Furthermore, we define the active index set of optical targets which is denoted by 𝐶

3.2. Compressive sampling matching pursuit (CoSaMP) algorithm The CoSaMP algorithm [22] is similar to OMP. The sparsity level 2𝑘 is used in CoSaMP whereas the sparsity level 𝑘 is used in OMP. The algorithm identifies the largest components of the proxy vector 𝑠 and these components are added to the components of the current approximation. Finally, it solves the least square step and forms the new approximation vector and update the residual. The main advantage of this algorithm is that it is stable when the signal is contaminated with noise. Moreover, it also offers the performance guarantee to recover the sparse signal. The implementation steps of this algorithm can be found in Appendix.

{ ( ) } 𝐶 = 𝑗 ∈ {1, 2, … ..𝑛} ∶ 𝛥𝛼 𝐫(𝑗) ≠ 0 (11) ( ) ( ) ( ) ( ) Recall that 𝑡 𝐫𝑗 ; 𝑝 = 𝑣 𝐫𝑗 ; 𝑝 𝛥𝛼 𝐫𝑗 , where 𝑣 𝐫𝑗 ; 𝑝 is the photon flux at position 𝐫𝑗 within the field of view i.e FOV. Since our main focus is to image the optical parameters of the targets inside the FOV of the medium reachable by the induced optical current, the optical photon distribution can be (assumed value throughout the FOV. ) to be a(positive ) More specifically, 𝑡 𝐫𝑗 ; 𝑝 = 0 if 𝛥𝛼 𝐫𝑗 = 0. Therefore, the sparsity can be written as 𝑘 = ‖𝑋‖0 = |𝐶|. From the above equations, the DOT imaging problem can be formulated as (𝑍0) ∶ min ‖𝑋‖0 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑀 = 𝑊 𝑋 + 𝐵

(13)

where the superscript † indicates the pseudo-inverse. However, the main goal here is to obtain the change in the absorption parameter { ( )}𝑘 distribution measure 𝛥𝛼 𝐫(𝑖) 𝑖=1 instead of the induced optical current { ( )}𝑘,𝑙 distribution measure 𝑡 𝐫(𝑖) ; 𝑝 𝑖,𝑝=1 (details can be found in Appendix). ( ) Note that for an obtained estimate 𝑡 𝐫𝑖 ; 𝑝 , we can easily calculate the ( ) total optical flux 𝑣 𝐫𝑖 ; 𝑝 by using the Foldy–Lax equation. Hence, 𝛥𝛼 can be calculated using least squares [33] ))∗ ( ) ∑𝑙 ( ( 𝑡 𝐫(𝑗) ; 𝑝 ( ) 𝑝=1 𝑣 𝐫(𝑗) ; 𝑝 , 𝑗 = 1, … , 𝑘, 𝛥𝛼 𝐫(𝑗) = (14) )|2 ∑𝑙 | ( 𝑝=1 || 𝑣 𝐫(𝑗) ; 𝑝 ||

(6)

𝑀 =𝑊𝑋+𝐵 where the scattered optical flux 𝑀 ∈ 𝐵 ∈ 𝑅𝑚×𝑟 are calculated by ( ) ( ) ⎛ 𝛽 𝐫𝑑 ; 1 … 𝛽 𝐫𝑑 ; 𝑙 ⎞ 1 1 ⎜ ⎟ 𝑀 = ⎜ ( ⋮ ) ⋱ ( ⋮ )⎟ ⎜𝛽 𝐫 ; 1 ⋯ 𝛽 𝐫 ; 𝑙 ⎟ ⎝ ⎠ ( 𝑑𝑚 ) ( 𝑑𝑚 ) ⎛ 𝑏 𝐫𝑑 ; 1 … 𝑏 𝐫𝑑 ; 𝑙 ⎞ 1 1 ⎟ ⎜ 𝐵 = ⎜ ( ⋮ ) ⋱ ( ⋮ )⎟ ⎜𝑏 𝐫 ; 1 ⋯ 𝑏 𝐫 ; 𝑙 ⎟ ⎝ ⎠ 𝑑𝑚 𝑑𝑚

𝑋̂ 𝐶 = 𝑊𝐶† 𝑀

3.3. Stagewise orthogonal matching pursuit (StOMP) algorithm The StOMP algorithm [23] is based on Greedy OMP. The main difference between the two algorithms is that the StOMP uses many vectors to construct a solution vector while OMP uses one vector at a time. The main idea behind the StOMP is that it compares the resultant values of the dot product of 𝑀 with the columns of 𝑊 and then selects all resultant vectors above a preset threshold value. The least square step is used to find the solution vector and the residue vector is updated. The advantage of this algorithm is that it finds a good solution in less number of iterations. The drawback of this algorithm is that the extensive computation of an appropriate threshold value is required. The details of this algorithm can be found in Appendix.

(12)

166

B.P.V. Dileep et al.

Optics Communications 410 (2018) 164–173

Fig. 1. The paraffin wax sample.

3.4. Regularized orthogonal matching pursuit (ROMP) algorithm Fig. 3. The schematic diagram of the DOT imaging setup.

The ROMP algorithm [24] is also based on Greedy OMP. Like StOMP, this algorithm uses several vectors at each step to construct the solution vector. The difference between StOMP and ROMP is that ROMP does not require a preset threshold value rather it uses the resultant vectors which have a dot product similar to the required vector. This algorithm uses all resultant vectors which have a dot product greater than half the size of the biggest dot product. The advantage of this step is that the resultant vectors can make a similar contribution to the resultant vector which is required. Then the residue vector is updated. The preset threshold value is not used in ROMP, so it is a great advantage over StOMP. The details of this algorithm can be found in Appendix.

4.1. Preparation of the rectangular phantom A rectangular phantom (as shown in Fig. 1) with dimensions 0.5 cm× 0.5 cm (cm indicates centimeter) having center at (0.25, 0.25) cm is considered in the experimental study. The paraffin wax is used as the background material in the phantom preparation. We first heat the paraffin wax, and the heated solution is poured into the rectangular box to form the rectangular shape of the phantom. After solidification of this heated solution in the rectangular box, we made a small hole in it which corresponds to the inclusion. Blue ink is added to the hole made on the rectangular phantom for imitating the inclusion with a high absorption parameter. The center of the inclusion is located at (0.45, 0.45) cm with a diameter of 0.08 cm. Here, an optical wavelength of 680 nm has been used as the laser source. The experimentally estimated values of the absorption and scattering coefficients of the paraffin wax are 0.25 cm−1 and 20 cm−1 respectively.

3.5. Simultaneous orthogonal matching pursuit (S-OMP) algorithm The S-OMP [25] is a Greedy algorithm. This algorithm picks all the columns of the measurement matrix 𝑊 at a time to form the solution vector whereas the OMP picks one column at a time. The S-OMP starts with a vector 𝑀 and identifies the largest components in a vector 𝑤∗𝑓 𝑀𝐽 and update the support set 𝑄. Then it uses the orthogonal projection rule to subtract its contribution to 𝑀, and the least square step has been used to compute the solution vector. The main advantage of S-OMP is that it is faster and easier to implement. The implementation steps of this algorithm can be found in Appendix.

4.2. Experimental setup The experimental setup for DOT is illustrated in Fig. 2 and the schematic diagram of the same is shown in Fig. 3. The main components involved in this setup are the laser source, rectangular box to hold the sample under study and the photodetector. The semiconductor laser having 2 mW power with the operating wavelength of 680 nm has been used as the light source. The plastic fibers are attached to the laser source for illuminating the sample at different locations. The length and the core diameter of each plastic fiber used are about 30 cm and 1 mm respectively. The rectangular box contains holes with equal spacing

4. Experimental studies for DOT The imaging setup and rectangular phantom preparation are being described as a part of the experimental studies for DOT. The DOT imaging setup has been used to study the different applications of the DOT like brain and breast imaging.

Fig. 2. The DOT imaging setup.

167

B.P.V. Dileep et al.

Optics Communications 410 (2018) 164–173

Table 1 Imaging geometry for both simulated and experimental study with single inclusion. –

𝐃

𝐒

𝐃

𝐒

𝐃

–

𝐃

𝐓

𝐓

𝐓

𝐓

𝐓

𝐃

𝐒

𝐓

𝐓

𝐓

𝐓

𝐓

𝐒

𝐃

𝐓

𝐓

𝐓

𝐓

𝐓

𝐃

𝐒

𝐓

𝐓

𝐓

𝐓∗

𝐓∗

𝐒

𝐃

𝐓

𝐓

𝐓

𝐓∗

𝐓∗

𝐃

–

𝐃

𝐒

𝐃

𝐒

𝐃

–

between them at different locations of XY-plane. The sample is placed inside this rectangular box and the fibers of source and detector are attached to the holes of this box to form the optical geometry for source–detector. The collection fiber are connected to a circular ring according to their position on the measurement surface of the sample. The photodetector is placed near the collection fiber for obtaining the measurements. The optical intensity of the source is varied so that the detector intensities are within the range of the photodetector values. The laser source is rotated to scan the whole boundary and the measurements are taken by the photodetector. Each time the sample is translated and rotated by hand to change the source position. For each source position, the photodetector is moved to each of the fibers mounted on a circular disk to take one reading at a time. The sample is translated and rotated 8 times which corresponds to 8 source positions. In this way the coordinate space of the sample is fixed. In total, 3 sets of measurement readings are taken and the average of these readings were taken as the final boundary measurements. The whole experiment is performed in a dark room environment so the background count can be assumed as zero. However, the final measurement readings at the collection fiber for the imaging setup used are affected by several other factors like the coupling from a light source to source fiber, the light coupling from source fiber to detector fiber via the sample, and the coupling efficiency of the fibers. These factors may vary a little bit even if the necessary precautions are being taken to keep the angle and position invariant in the measurement surface of the medium. Hence, calibration of the data needs to be performed by considering the whole imaging setup. First the experiment is performed with a homogeneous phantom having known absorption values and the detector readings are measured using the photodetector. The scale factor is defined for calibration, and it has to be changed iteratively until the simulated data is matched with the experimental data with some tolerance. Once the simulated and experimental data are equal, then the corresponding scale factor is the scale factor for the calibration, which take care of the uncertainties that are involved in the source and detector couplings. This value for the scale factor has been verified by measuring the input and output intensities using the photodetector. After the calibration is done, the fibers are kept inside the holes and the entire experiment is carried out ensuring that the imaging setup is not disturbed. The uniformity of the signal is achieved at the collection fiber by adjusting the fibers manually, while the light source is moved from one fiber to another fiber so that the total optical intensity would be the same. The continuous mode measurement is considered in this study to make the DOT system simple, portable, and low-cost.

( ) Fig. 4. Ground truth change in the absorption profile for 𝛥𝜇𝑎 cm−1 for both simulated and experimental case.

denoted by S and 𝐃. The area to be imaged, i.e., target is represented by 𝐓 and the location of the inclusion is denoted by 𝑇 ∗ respectively. The sources and detectors are alternatively placed along the boundary of the geometry. We have used 8 sources and 12 detectors in this geometry. The absorption and scattering parameters of the homogeneous background are 0.25 cm−1 and 20 cm−1 , whereas the maximum change in the absorption coefficient 𝛥𝜇𝑎 is 0.25 cm−1 . The distribution of the values of 𝛥𝜇𝑎 at the targets is called ground truth change in the absorption profile as depicted in Fig. 4. Since 𝛥𝜇𝑎 is sparse, so its minimum value is zero at the targets and its maximum value at the location of the inclusion is same as the background absorption coefficient about 0.25 cm−1 . This ground truth change in the absorption profile is considered as reference image for 𝛥𝜇𝑎 and will be used for the comparison with the reconstructed images of 𝛥𝜇𝑎 . Only single inclusion is considered in this geometry. 5.2. Reconstructed images for 𝛥𝜇𝑎 By employing the above imaging geometry as shown in Table 1, we perform the 2D simulation of the DOT problem using the Greedy algorithms based on compressive sensing framework and compared with conventional methods. The reconstructed images for 𝛥𝜇𝑎 for both simulated and experimental cases are presented in Figs. 5–8. It has been observed from Figs. 5–8 that the reconstructed images for the simulated and experimental cases of the rectangular sample with single inclusion is found to be satisfactory. The red portion in Figs. 5–8 indicates the inclusion. By comparing Fig. 4 with Figs. 5–8, for the OMP case, the location of the inclusion is identified perfectly in simulated and experimental images and matched with the ground truth image, but the reconstructed value for 𝛥𝜇𝑎 does not match with the ground truth. For the CoSaMP case, the location of the inclusion is identified perfectly in simulated and experimental images and matched with the ground truth image and the reconstructed value for 𝛥𝜇𝑎 is near to the ground truth. For the StOMP case, the deviation of the location of the inclusion is observed in simulated image and location of the inclusion is identified perfectly in the experimental image but reconstructed value for 𝛥𝜇𝑎 does not match with the ground truth. For the ROMP case, the location of the inclusion is identified perfectly in simulated and experimental images, but the reconstructed value for 𝛥𝜇𝑎 is very high compared to the ground truth. For the S-OMP case, the location of the inclusion is identified correctly in a simulated image and the deviation of the location of the inclusion is observed in an experimental image but the reconstructed value for 𝛥𝜇𝑎 does not match with the ground truth. For the conventional methods like least square method and TSVD, the location of the inclusion and the reconstructed value for 𝛥𝜇𝑎 in simulated and experimental images does not match with the ground truth image. Overall the results obtained with CoSaMP is much better compared to the other Greedy algorithms as well as conventional methods. These reconstructions show the potential of the Greedy algorithms in DOT.

5. Results and discussion The DOT imaging geometry is elaborated in this section for 2D reconstruction of the DOT problem. 5.1. DOT imaging geometry for 2D simulation The imaging geometry for DOT is illustrated in Table 1. The total geometry is divided into several grids. The sources and detectors are 168

B.P.V. Dileep et al.

Optics Communications 410 (2018) 164–173

(a) OMP.

(a) S-OMP.

(b) CoSaMP.

(b) Least square method.

(c) StOMP.

(c) TSVD. ( ) Fig. 6. Simulated images for 𝛥𝜇𝑎 cm−1 using (a) S-OMP, (b) least square method, and (c) TSVD.

square error (NMSE), structural similarity index (SSIM) and peak signal to noise ratio (PSNR) which are calculated as follows. The MSE is given by 𝑁 )2 1 ∑( 𝑟 𝜇 − 𝜇𝑎𝑎 𝑁 𝑖=1 𝑎

The NMSE is computed as follows 𝑁 ( 𝑟 𝑎 )2 1 ∑ 𝜇𝑎 − 𝜇𝑎 𝑁 𝑖=1 𝜇𝑎𝑎

(d) ROMP. ( ) Fig. 5. Simulated images for 𝛥𝜇𝑎 𝑐𝑚−1 using (a) OMP, (b) CoSaMP, (c) StOMP, and (d) ROMP.

(15)

(16)

where 𝜇𝑎𝑟 and 𝜇𝑎𝑎 are the reconstructed and actual absorption coefficients, 𝑁 is the number of absorption coefficients. The SSIM is used to measure the similarity between the original image and the reconstructed image. For comparing two images 𝑎 (ground truth image is shown in Fig. 4) and 𝑏 (any reconstructed image is shown in Figs. 5–8), SSIM is computed as

5.3. The performance evaluation metrics The performance evaluation metrics play a significant role in the DOT reconstruction. There are various performance metrics that have been used to evaluate the performance of the Greedy algorithms and conventional methods in DOT reconstruction. The performance metrics considered in this paper are mean square error (MSE), normalized mean

𝑆(𝑎, 𝑏) =

(2𝜇𝑎 𝜇𝑏 + 𝑐1 )(2𝜎𝑎𝑏 + 𝑐2 ) (𝜇𝑎2 + 𝜇𝑏2 + 𝑐2 )(𝜎𝑎2 + 𝜎𝑏2 + 𝑐2 )

(17)

where 𝜎𝑎 , 𝜎𝑏 , 𝜇𝑎 , 𝜇𝑏 , 𝜎𝑎𝑏 are the standard deviation, mean, covariance of image a and b, respectively. The calculated SSIM is a scalar value 169

B.P.V. Dileep et al.

Optics Communications 410 (2018) 164–173

(a) OMP.

(a) S-OMP.

(b) CoSaMP.

(b) Least square method.

(c) TSVD. (c) StOMP. ( ) Fig. 8. Experimental images for 𝛥𝜇𝑎 cm−1 using (a) S-OMP, (b) least square method, and (c) TSVD.

The PSNR is calculated as ( ) (𝑝𝑒𝑎𝑘𝑣𝑎𝑙)2 10log10 𝜂

where 𝑝𝑒𝑎𝑘𝑣𝑎𝑙 is the peak value taken from the range of the image, and 𝜂 is the MSE. The PSNR is measured in decibels (dB). The comparison of performance metrics for the simulated and experimental cases are provided in Tables 2 and 3 respectively. It has been evident from Tables 2 and 3 that the MSE and NMSE are significantly less for CoSaMP in comparison to other Greedy algorithms and conventional methods. It is also seen that the SSIM and PSNR are significant for CoSaMP as compared to other Greedy algorithms and conventional methods. Overall the CoSaMP performs better over other Greedy algorithms and conventional methods. The overall performance of the Greedy algorithms concerning the error calculations in this reconstruction is satisfactory. Also, the errors produced by these algorithms are comparatively small. These results again show the applicability of the above mentioned algorithms based on CS as compared to the conventional DOT methods in the DOT imaging. Note that the overall reconstruction for the absorption coefficient corresponding to 𝜇𝑎𝑟 can be achieved by adding the change in the absorption coefficient 𝛥𝜇𝑎 with

(d) ROMP. ( ) Fig. 7. Experimental images for 𝛥𝜇𝑎 cm−1 using (a) OMP, (b) CoSaMP, (c) StOMP, and (d) ROMP.

which lies between −1 and +1. A high value of SSIM indicates that the reconstructed image is much closer to the original image i.e, the structural distortion in the reconstructed absorption image is less. Similarly, the structural dissimilarity is computed by using (1 − 𝑆(𝑎, 𝑏)) 2

(19)

(18) 170

B.P.V. Dileep et al.

Optics Communications 410 (2018) 164–173 Table 2 Comparison of performance metrics for the simulated case. Performance metrics

OMP

CoSaMP

StOMP

ROMP

S-OMP

Least square method

TSVD

MSE NMSE SSIM PSNR (dB)

0.0003 0.0392 0.4447 24.1502

0.0002 0.0349 0.5033 24.8029

0.0007 0.0796 0.3489 23.1802

0.0004 0.0463 0.4361 24.1966

0.0003 0.004 0.4868 24.7667

0.1276 0.4779 0.0394 6.3579

0.1048 0.3977 0.0089 12.1634

Table 3 Comparison of performance metrics for the experimental case. Performance metrics

OMP

CoSaMP

StOMP

ROMP

S-OMP

Least square method

TSVD

MSE NMSE SSIM PSNR (dB)

0.0013 0.0981 0.5378 25.0882

0.0011 0.0964 0.5443 25.14

0.004 0.2296 0.1820 23.0647

0.0054 0.2752 0.3562 19.5168

0.002 0.2271 0.1084 20.3657

0.1413 0.5016 0.0024 7.5270

0.1732 0.6162 0.0308 6.3035

{ ( )}𝑘,𝑙 A.2. Expanded form of 𝑡 𝐫(𝑖) ; 𝑝 𝑖,𝑝=1

the background absorption coefficient 𝜇0 . This also avoids the division by zero problem in computing NMSE. The future work will be in this direction.

( ) ( ) ( ) Since 𝑡 𝐫𝑗 ; 𝑝 = 𝑣 𝐫𝑗 ; 𝑝 𝛥𝛼 𝐫𝑗 , so we have the expanded form of the vector equation ( ) ( ) 𝑣 𝐫(𝑖) ; 1 ⎤ ⎡ 𝑡 𝐫(𝑖) ; 1 ⎤ ( )⎡ ⎢⋮ ⎥ = 𝛥𝛼 𝐫 ⎢⋮ ⎥ , 𝑖 = 1, … , 𝑘, (23) (𝑖) )⎥ )⎥ ⎢ ( ⎢ ( ⎣ 𝑡 𝐫(𝑖) ; 𝑙 ⎦ ⎣ 𝑣 𝐫(𝑖) ; 𝑙 ⎦

6. Conclusion This paper describes diffuse optical tomographic reconstruction using compressive sensing framework. As the optical parameters of inclusion are sparse, the DOT reconstruction problem can be formulated as a SMV and MMV problem using CS framework. In this paper we studied the feasibility of Greedy algorithms for DOT reconstruction with rectangular sample having single inclusion. We have compared the Greedy algorithms with conventional DOT reconstruction methods like least square method and TSVD which shows marked improvement in DOT. The validation of Greedy algorithms has been performed experimentally on a paraffin wax rectangular phantom. The overall reconstruction is adequate for indicating the irregularities in the case of simulated and experimental data. The performance metrics such as MSE, NMSE, SSIM, and PSNR are used to study the performance of the Greedy algorithms and conventional methods in this reconstruction. Extensive numerical studies shows the potential of Greedy algorithms based on CS in DOT. These Greedy algorithms outperform over the conventional DOT imaging methods in terms of computational efficiency. Specially the CoSaMP algorithm has outperformed over all the methods as expected. This reconstruction may not produce accurate images when the complex geometries are considered for the study. The future work will be the application of CS based diffuse optical tomographic reconstruction of absorption profile for different geometries with more number of inclusions.

A.3. Implementation steps of the greedy algorithms The implementation steps of the Greedy algorithms are given below. OMP The steps of OMP are as follows. Input: 1. An input measurement matrix 𝑊 2. The data vector 𝑀 3. The sparsity level 𝑘 Output: 1. 2. 3. 4.

An obtained signal 𝑧̂ A set 𝜆𝑘 contains 𝑘 components from {1, … , 𝑒} The approximation vector 𝑑𝑘 of the data 𝑌 The residual 𝑎𝑘

Appendix Steps: A.1. Foldy–Lax equation 1. Initialize the residual 𝑎0 = 𝑀, the index set 𝛤0 = 0 and the iteration counter 𝑖 = 1. 2. Find the index 𝛾𝑖 that solves the given optimization problem | | easily 𝛾𝑖 = arg max𝑗=1,…,𝑒 |≺ 𝑎𝑖−1, 𝛹𝑗 ≻|. where 𝛹𝑗 is the columns | | of 𝑊 . If the maximum value occurs at the multiple indices, halt the tie immediately. 3. Add the index [set and the ] matrix with chosen atoms 𝜆𝑖 = 𝜆𝑖−1 ∪ { } 𝛾𝑖 and 𝑊𝑖 = 𝑊𝑖−1 𝛹𝛾𝑖 . We assume 𝑊0 to be an empty matrix. 4. Solve the following least square problem for an obtained new ‖ signal 𝑋𝑖 = arg min𝑋 ‖ ‖𝑀 − 𝑊𝑖 𝑋 ‖2 . 5. Find the new approximation vector and the new residual 𝑑𝑖 = 𝑊𝑖 𝑋𝑖 and 𝑎𝑖 = 𝑀 − 𝑑𝑖 . 6. Perform 𝑖 = 𝑖 + 1 and go to step 2 if 𝑖 < 𝑘. 7. The obtained signal 𝑧̂ has the nonzero indices that correspond to the elements in 𝜆𝑘 . The value of 𝑧̂ in the element 𝛾𝑗 equals the 𝑗th element of 𝑋𝑖 .

The Foldy–Lax equation [33–36] is given by ( ) ( ) ∑ ( ) ( ) ( ) 𝑣 𝐫(𝑗) ; 𝑝 = 𝑣0 𝐫(𝑗) ; 𝑝 + 𝑔0 𝐫(𝑗) , 𝐫(𝑖) 𝑣 𝐫(𝑖) ; 𝑝 × 𝛥𝛼 𝐫(𝑖) ,

(20)

𝑖≠𝑗

where 𝑖, 𝑗 = 1, 2, … , 𝑘 and 𝑝 = 1, 2, … , .𝑙. Here, we consider the multiple scattering of photons from all the possible scattering points except the self field for avoiding the singularity that occurred in this problem. ( ) The homogeneous green’s function 𝑔0 𝐫, 𝐫 ′ can be described as ) ( 2 ) ( ) 1 ( 𝛿 𝐫 − 𝐫′ ∇ − 𝑘20 𝑔0 𝐫, 𝐫 ′ = − 𝐷0

(21)

And the corresponding homogeneous optical flux 𝑣0 (𝐫; 𝑝) can be calculated by using 𝑣0 (𝐫; 𝑝) =

∫

( ) ( ) 𝑔0 𝐫, 𝐫 ′ 𝑞 𝐫 ′ ; 𝑝 𝑑𝐫 ′ , 𝑝 = 1, 2, … , 𝑙

(22)

where 𝐫 = 𝐫𝑗 and 𝐫 ′ = 𝐫𝑖 . 171

B.P.V. Dileep et al.

Optics Communications 410 (2018) 164–173

CoSaMP

S-OMP

The steps of the CoSaMP are as follows. Input:

The S-OMP can be implemented by the following steps. Input:

1. An input sensing matrix 𝑊 , the data vector 𝑀 2. The sparsity level 𝑘

1. An input sensing matrix 𝑊 , the data vector 𝑀 2. The sparsity level 𝑘 is the number of indices in the support set 𝑄

Output:

Output:

1. An approximation vector 𝑑 2. The residue vector 𝑎

1. Approximation vector 𝑋 Steps:

Steps: 1. Set 𝑀0 = 𝑀 and 𝑄 = 0, where 𝑄 is the support ( set. ) { }𝐽 2. After 𝐽 iterations, 𝑄𝑗 = 𝑓𝑗 𝑗=1 , and 𝑀𝐽 = 𝐼 − 𝑃𝑄𝐽 𝑀, where { }𝐽 𝑃𝑄𝐽 is the orthogonal projection onto span 𝑤𝑓𝑗 . Here, 𝑤𝑓𝑗

1. Initialize 𝑎 = 𝑀, the iteration counter 𝑖 = 0 and initial approximation 𝑑 0 = 0. 2. Set 𝑖 = 𝑖 + 1 and(form ) a proxy vector s = (W*a) and identify 𝑠𝑢𝑝𝑝 is the support. support 𝛺 = 𝑠𝑢𝑝𝑝 𝑠2𝑘 ( , where ) 3. Set 𝑇 = 𝛺 ∪ 𝑠𝑢𝑝𝑝 𝑑 𝑖−1 , form the approximation vectors as 𝐸𝑇 = 𝑊𝑇† 𝑎 and 𝐸𝑇 𝑐 = 0. where 𝑊)𝑇† is the pseudo inverse. ( 𝑖 4. Set 𝑑 = 𝐸𝑘 and 𝑎 = 𝑀 − 𝑊 ∗ 𝑑 𝑖 .

𝑗=1

is the columns of 𝑊 . ‖ ‖ ‖ ‖ 3. Select 𝑓𝐽 +1 such that ‖𝑤∗𝑓 𝑀𝐽 ‖ = max1≤𝑓 ≤𝑛 ‖𝑤∗𝑓 𝑀𝐽 ‖ , and ‖2 ‖ ‖2 { ‖ 𝐽}+1 also set 𝑄𝐽 +1 = 𝑄𝐽 ∪ 𝑓𝐽 +1 . 4. Obtain 𝑋 = 𝑊𝑄† ∗ 𝑀, where 𝑊𝑄† is the pseudo inverse.

StOMP

References

The StOMP can be implemented by the following steps. Input:

[1] D.A. Boas, D.H. Brooks, E.L. Miller, C.A. DiMarzio, M. Kilmer, R.J. Gaudette, Q. Zhang, Imaging the body with diffuse optical tomography, IEEE Signal Process. Mag. 18 (6) (2001) 57–75. [2] S. Arridge, J. Schotland, Optical tomography: Forward and inverse problems, Inverse Problems 25 (12) (2009). http://dx.doi.org/10.1088/0266-5611/25/12/123010. [3] S.R. Arridge, Optical tomography in medical imaging, Inverse Problems 15 (2) (1999) R41. [4] R.J. Gaudette, D.H. Brooks, C.A. DiMarzio, M.E. Kilmer, E.L. Miller, T. Gaudette, D.A. Boas, A comparison study of linear reconstruction techniques for diffuse optical tomographic imaging of absorption coefficient, Phys. Med. Biol. 45 (4) (2000) 1051. URL http://stacks.iop.org/0031-9155/45/i=4/a=318. [5] R. Weissleder, Scaling down imaging: molecular mapping of cancer in mice, Nature Rev. Cancer 2 (1) (2002) 11–18. [6] R. Weissleder, V. Ntziachristos, Shedding light onto live molecular targets, Nature Med. 9 (1) (2003) 123–128. [7] S.D. Konecky, G.Y. Panasyuk, K. Lee, V. Markel, A.G. Yodh, J.C. Schotland, Imaging complex structures with diffuse light, Opt. Express 16 (7) (2008) 5048–5060. [8] M.A. Franceschini, D.A. Boas, Noninvasive measurement of neuronal activity with near-infrared optical imaging, Neuroimage 21 (1) (2004) 372–386. [9] A. Gibson, J. Hebden, S.R. Arridge, Recent advances in diffuse optical imaging, Phys. Med. Biol. 50 (4) (2005) R1. [10] J. Beuthan, Optical diagnostics—state of the art, Med. Laser Appl. 22 (1) (2007) 43–47. [11] S.R. Arridge, P. van der Zee, M. Cope, D.T. Delpy, Reconstruction methods for infrared absorption imaging, in: Optics, Electro-Optics, and Laser Applications in Science and Engineering, International Society for Optics and Photonics, 1991, pp. 204–215. [12] Y. Yao, Y. Wang, Y. Pei, W. Zhu, R.L. Barbour, Frequency-domain optical imaging of absorption and scattering distributions by a born iterative method, J. Opt. Soc. Am. A 14 (1) (1997) 325–342. [13] J.C. Ye, K.J. Webb, C.A. Bouman, R.P. Millane, Optical diffusion tomography by iterative-coordinate-descent optimization in a bayesian framework, J. Opt. Soc. Am. A 16 (10) (1999) 2400–2412. [14] V.A. Markel, J.C. Schotland, Inverse problem in optical diffusion tomography. ii. role of boundary conditions, J. Opt. Soc. Am. A 19 (3) (2002) 558–566. [15] V.A. Markel, V. Mital, J.C. Schotland, Inverse problem in optical diffusion tomography. iii. inversion formulas and singular-value decomposition, J. Opt. Soc. Am. A 20 (5) (2003) 890–902. [16] V. Ntziachristos, R. Weissleder, Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized born approximation, Opt. Lett. 26 (12) (2001) 893–895. [17] N. Cao, M. Ortner, A. Nehorai, Solutions for diffuse optical tomography using the feynman-kac formula and interacting particle method, in: Proc. SPIE BIOS, vol. 7, 2007. [18] D.L. Donoho, Compressed sensing, IEEE Trans. Inf. Theory 52 (4) (2006) 1289–1306. [19] J. Chen, X. Huo, Theoretical results on sparse representations of multiplemeasurement vectors, IEEE Trans. Signal Process. 54 (12) (2006) 4634–4643. [20] R.G. Baraniuk, Compressive sensing, IEEE Signal Process. Mag. 24 (4) (2007). [21] J.A. Tropp, A.C. Gilbert, Signal recovery from random measurements via orthogonal matching pursuit, IEEE Trans. Inform. Theory 53 (12) (2007) 4655–4666. http: //dx.doi.org/10.1109/TIT.2007.909108.

1. An input sensing matrix 𝑊 , the data vector 𝑀 2. The threshold value 𝑛 Output: 1. An approximation vector 𝑋 or index vector 𝑉 2. The residue vector 𝑎 Steps: 1. Set the residual 𝑎0 = 𝑀, the iteration counter 𝑖 = 0 and index set 𝑉 = 0. 2. Create 𝐽 which consists of the location of all entries in the vector 𝑈 above a preset value 𝑛 and form that vector as 𝑈 = (𝑊 𝑇 ∗ 𝑎𝑖 ), where 𝑊 𝑇 is the transpose of matrix 𝑊 . 3. Set the index set by 𝑉 = 𝑉 ∪ 𝐽 and residual by 𝑋𝑖 = ‖ min𝑋 ‖ ‖𝑀 − 𝑊𝑉 𝑋 ‖2 and 𝑎𝑖+1 = 𝑀 − 𝑊 𝑋𝑖 . ROMP The steps of ROMP are as follows. Input: 1. An input sensing matrix 𝑊 , the data vector 𝑀 2. The sparsity level 𝑘 Output: 1. An approximation vector 𝑋 or index set 𝑉 2. The residue vector 𝑎 Steps: 1. Set the residual 𝑎0 = 𝑀, the iteration counter 𝑖 = 0, and index set 𝑉 = 0. 2. Choose a set 𝐽 consisting of 𝑛 biggest coordinates in the vector 𝑈 = (𝑊 𝑇 ∗ 𝑎𝑖 ), where 𝑊 𝑇 is the transpose of matrix 𝑊 . 3. Regularize and find among all subsets 𝐽0 ⊂ 𝐽 , the maximum ‖ ‖ ‖ ‖ ‖ 𝑈 ‖ and ‖𝑈𝐽0 ‖ , where 𝐽0 is defined as 𝑖, 𝑗 ∈ 𝐽0 if ‖ ‖ 𝑈𝑖 ‖ ≤ 2 ‖ ‖ ‖2 ‖ 𝑗‖ 𝑈𝑖 ∈ 𝑈𝐽𝑜 if 𝑝 ∈ 𝐽0 . 4. Update the index set by 𝑉 = 𝑉 ∪ 𝐽0 and the residual by 𝑋𝑖 = ‖ min𝑋 ‖ ‖𝑀 − 𝑊𝑉 𝑋 ‖2 and 𝑎𝑖+1 = 𝑀 − 𝑊 𝑋𝑖 .

172

B.P.V. Dileep et al.

Optics Communications 410 (2018) 164–173 [29] A. Hore, D. Ziou, Image quality metrics: psnr vs. ssim, in: Pattern Recognition (Icpr), 2010 20th International Conference on, IEEE, 2010, pp. 2366–2369. [30] D.L. Donoho, Neighborly polytopes and sparse solutions of underdetermined linear equations. [31] D.L. Donoho, M. Elad, Optimally sparse representation in general (non-orthogonal) dictionaries via 1 minimization, Proc. Natl Acad. Sci. USA 100 (2002) 2197–2202. [32] E.J. Candès, J. Romberg, T. Tao, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Trans. Inf. Theory 52 (2) (2006) 489–509. [33] J.C. Ye, S.Y. Lee, Y. Bresler, Exact reconstruction formula for diffuse optical tomography using simultaneous sparse representation, in: 2008 5th IEEE International Symposium on Biomedical Imaging: From Nano To Macro, IEEE, 2008, pp. 1621–1624. [34] L. Tsang, J.A. Kong, K.-H. Ding, C.O. Ao, Scattering of Electromagnetic Waves, Numerical Simulations, vol. 25, John Wiley & Sons, 2004. [35] M.I. Mishchenko, L.D. Travis, A.A. Lacis, Multiple Scattering of Light By Particles: Radiative Transfer and Coherent Backscattering, Cambridge University Press, 2006. [36] A.C. Fannjiang, Compressive inverse scattering: I. high-frequency simo/miso and mimo measurements, Inverse Problems 26 (3) (2010) 035008.

[22] D. Needell, J.A. Tropp, Cosamp: iterative signal recovery from incomplete and inaccurate samples, Appl. Comput. Harmon. Anal. 26 (3) (2009) 301–321. [23] D.L. Donoho, Y. Tsaig, I. Drori, J.-L. Starck, Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit, IEEE Trans. Inform. Theory 58 (2) (2012) 1094–1121. [24] D. Needell, R. Vershynin, Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit, Found. Comput. Math. 9 (3) (2009) 317–334. [25] J.A. Tropp, A.C. Gilbert, M.J. Strauss, Algorithms for simultaneous sparse approximation. part i: greedy pursuit, Signal Process. 86 (3) (2006) 572–588. [26] R. Gribonval, H. Rauhut, K. Schnass, P. Vandergheynst, Atoms of all channels, unite! average case analysis of multi-channel sparse recovery using greedy algorithms, J. Fourier Anal. Appl. 14 (5–6) (2008) 655–687. [27] O. Lee, J.M. Kim, Y. Bresler, J.C. Ye, Compressive diffuse optical tomography: noniterative exact reconstruction using joint sparsity, IEEE Trans. Med. Imaging 30 (5) (2011) 1129–1142. [28] Z. Wang, A.C. Bovik, H.R. Sheikh, E.P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE Trans. Med. Imaging 13 (4) (2004) 600–612.

173