1.

Introduction

Diffuse optical imaging has received heightened attention in the last decade because of its capability to provide functional images of the tissue under investigation using nonionizing near-infrared light (600–1000 nm).^{1, 2} Specifically, imaging of brain and breast has been the primary applications of diffuse optical tomography.^{1, 2, 3} The critical step in obtaining these images is estimating internal distribution of optical properties of the tissue using the measurements made on the tissue boundary.^{4, 5} Because the scattering is the dominant mechanism for near-infrared light (NIR) interaction with tissue, the estimation problem, also known as the inverse problem, is nonlinear, ill-posed, and some times underdetermined.⁴ Thus, solving an inverse problem necessitates the use of computationally intensive models. The computed data using these models are matched with the experimental data iteratively in the least-squares sense to obtain the optical properties of the tissue.⁴

The major challenge in terms of obtaining these optical images in real time is the computational cost associated with these advanced computational models, because these are used repeatedly.^{2, 4, 6} These iterative techniques used in the inverse problem relies on the calculation of modeled data and Jacobian (or its variant) at each iteration to obtain an update of optical properties.^{4, 6, 7} Depending on the model used, these calculations can span up to several hours, especially in three-dimensional (3-D) cases.^{6, 8} There were attempts earlier to accelerate these calculations by using parallel computers, which have been shown to give a speedup of factor n, with n being the number of parallel processors used.⁶ The limitation with the use of parallel computers is the cost associated with achieving higher speedups and also the complex approach in parallelizing these computational models. Here, we aim to take advantage of general purpose graphics processing units (GPUs) in massively parallelizing these calculations. This work specifically aims to present accelerating the frequency domain diffuse optical image reconstruction using a cost-effective (∼$1200) programmable GPU (NVIDIA Tesla C 1060). This is achieved using an open-source GPU-accelerated linear algebra library (CULA) that utilizes the NVIDIA's compute unified device architecture (CUDA).⁹

Earlier works have shown the parallel computation capability of GPUs in performing high-speed Monte Carlo simulation of photon propagation in tissue and proven to give a speedup of ∼100 for simple cases^{10, 11} and ∼10 in heterogeneous tissues¹² when compared to the implementation on a modern central processing unit (CPU). Also, the standard filtered backprojection (Radon transform-based) algorithm used in computed tomography (CT) has been shown to give a speedups of 100 compared to standard CPUs.¹³ These kinds of massively parallel programmable GPUs have been used in the context of optical projection tomography and shown to give 300-fold acceleration compared to traditional CPUs.¹⁴ For iterative algorithms that are used in CT reconstruction, such as simultaneous algebraic reconstruction, the reported speedups are 12-fold per iteration, comparing GPU to CPU.¹⁵ The image reconstruction algorithms that use Fourier transform or its variant have been shown to give speedups on the order of 100 or more using GPUs, and the reconstruction algorithms that mainly involve matrix computations have been shown to give speedups on the order of 10 or more.^{14, 15} The main aim of this work is to demonstrate the parallel computing power of these GPUs for performing diffuse optical tomographic image reconstruction.

Here, we used a finite element (FE) method (FEM)–based computational model that solves the frequency-domain diffusion equation (DE) requiring complex arithmetic. The number of unknowns in the inverse problem will be equal to 2N _N (N _N absorption coefficients and N _N-diffusion coefficients), with N _N being the number of finite element nodes used in the FE mesh to obtain the modeled data (described later). To obtain the modeled data, Jacobian, and update of optical properties at every iteration, a set of linear system of equations (order of 2N _N) needs to be solved along with the matrix-matrix multiplications (order of 2N _N), which will take typically O[(2N _N)³] operations. We will show that these operations could be done in parallel using a GPU (NVIDIA Tesla C 1060) to accelerate the diffuse optical tomographic image reconstruction procedure.

As diffuse optical image reconstruction relies on numerical calculations, specifically system solve (i.e., solving x in Ax = b) and matrix multiplications. Several packages,^{16, 17, 18, 19, 20, 21, 22, 23, 24} listed in Table 1, exist that provide options for carrying out the calculations on GPUs. The FEM-based forward models based on frequency-domain DE results in sparse symmetric linear system of equations (complex type). This, in turn, limits the use of sparse numerical computations on GPUs in these cases (more discussion to follow). The only package that is capable of dealing with sparse system solve for the real type is Cusp,¹⁹ allowing only matrix-vector computations. It will be shown that even with the use of the full (nonsparse) matrices, the GPUs are capable of given an acceleration of up to 7 (for completing a start-to-end single iteration of the diffuse optical image reconstruction) compared to CPU sparse computations. Note that the Linux-based platform is used to carry out the computations performed in this work.

Table 1

Comparison of features of available CUDA-based numerical linear algebra packages along with their source.

2.

Methods

2.1.

Diffuse Optical Tomographic Image Reconstruction

Diffuse optical tomographic image reconstruction is typically performed using Newton-type algorithms, where the modeled data [G(μ) with μ representing the set of optical properties] is matched to the experimental data (y), iteratively, in the least-square sense.^{4, 7} The most popular technique among the full Newton-type algorithms is Levenberg–Marquardt minimization, described in detail in Ref. 7. The FEM-based frequency-domain diffusion model for the calculation of G(μ) is described in Refs. 25, 26; here, it is only briefly reviewed. The frequency-domain DE is given by²⁵

Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \qquad&& -\nabla D(r)\nabla \Phi (r,\omega)\, + \,\left({\mu _a (r)\, + \,\frac{{i\omega}}{c}} \right)\,\Phi (r,\omega)\nonumber\\ \qquad&&\quad= \,q_o (r,\omega), \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill & & -\nabla D\left(r\right)\nabla \Phi (r,\omega )+\left({\mu }_a\left(r\right)+\frac{i\omega }{c}\right)\Phi (r,\omega )\hfill \\ \hfill & & =q_o(r,\omega ),\hfill \end{array}$$

where Φ(r, ω) photon density (complex values) at position r for the light modulation frequency of ω (= 2πf, with f = 100 MHz). The light source, represented by q ₀(r, ω), is modeled as isotropic and c represents the speed of light in tissue. The absorption coefficient is represented by μ_a(r) and the diffusion coefficient by D(r), defined as

Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} D(r)\, = \,\frac{1}{{3[\mu _{\rm a} (r)\, + \,\mu '_{\rm s}(r)]}}, \end{equation}\end{document} $$D\left(r\right)=\frac{1}{3[{\mu }_{\mathrm{a}}\left(r\right)+{\mu }_{\mathrm{s}}^{\prime }\left(r\right)]},$$

with μ'_s(r) representing the reduced scattering coefficient. The μ in this work represents [D(r); μ_a(r)]. A Robin (type III) boundary condition is typically used to take care of the refractive-index mismatch at the tissue boundary.²⁷

In the FEM framework, the imaging domain is discretized into linear triangular elements (for two dimensions) or linear tetrahedral elements (for three dimensions) connected at N _N vertex nodes. Then the computational (forward) model for solving the diffusion equation reduces to^{25, 26}

Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} K\Phi \, = \,q, \end{equation}\end{document} $$K\Phi =q,$$

where K is known as the mass matrix with a dimension of N _N×N _N (symmetric matrix) and is a function of μ [i.e., Φ = K ⁻¹ q = F(μ)], with K assembled over all elements of the finite element mesh. q represents the forcing, including the source term [q ₀(r,ω)] and the boundary condition.²⁵ The modeled data [G(μ)] is obtained by sampling of Φ at the measurement position [i.e., G(μ) = S{Φ} = S{F(μ)}, where S represents the sampling matrix (containing source/detector positions) and F is the forward model].⁷ Note that K is highly sparse (with a banded structure in case of bandwidth-optimized FE meshes) and typically sparse matrix solvers are used to obtain Φ. This process of solving for Φ [Eq. 3] involves a decomposition method (iterative procedure) preambled by a preconditioning step, because K is a large sparse complex matrix with a high numerical condition number.^{5, 25, 28} In this work, the Jacobi preconditioner (simplest of all preconditioners) has been used, where the preconditioned matrix (P) is a diagonal matrix consisting of diagonal values of K.^{29, 30} The inverse of P (which is one over the diagonal values of K) is left multiplied in Eq. 3 to form a better conditioned linear system of equations. This results in an increased rate of convergence in solving for Φ. The preconditioned linear system of equations is given by

Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \,P^{- 1} K\Phi \, = \,\,P^{- 1} q, \end{equation}\end{document} $$P^{-1}K\Phi =P^{-1}q,$$

with P ^–1 defined as^{29, 30}

Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} P_{i,j}^{- 1} \, = \left\{{\begin{array}{@{}l@{\quad}ll} \displaystyle\frac{1}{{K_{i,i}}}&{\rm if}&i = j \\ 0 & {\rm if}&i \ne j. \\ \end{array}} \right. \end{equation}\end{document} $$P_{i,j}^{-1}=\left\{\begin{array}{ccc}\hfill {\displaystyle \frac{1}{K_{i,i}}}& \hfill \mathrm{if}& \hfill i=j\\ \hfill 0& \hfill \mathrm{if}& \hfill i\ne j.\end{array}\right.$$

Because the complex sparse arithmetic is not yet supported in NVIDIA's CUDA, the GPU implementation of the same is performed by converting K to a full (nonsparse) matrix.

The most important steps pertaining to the image reconstruction procedure are given as a flowchart in Fig. 1. The iterative image reconstruction procedure starts with an initial guess for the optical properties (μ₀) typically obtained using the calibration procedure of experimental data (y).^{31, 32} Using this μ₀, the forward model is solved to obtain G(μ) and, more importantly, the Jacobian (J = ∂G(μ)/∂μ), which gives the rate of change in the modeled data with respect to optical properties. J is typically obtained using the adjoint formulation,²⁵ and the most important computations pertained to the calculation of J are given in Fig. 2. The computation times for each step (in percent) in calculation of J with adjoint formulation is also given in Fig. 2. It could be seen that calculation of J needs both Φ and Φ* (adjoint fluence, obtained by interchanging the source and detector positions) and typically takes ∼35% of time in the preconditioning {i.e., calculation of P ⁻¹ K with the order of calculations as O[(2N _N)³]}, as both P ⁻¹ and K are complex [for the real case, it will be O(N _N ³)].²⁹ Obtaining Φ in Eq. 4 is performed using LU decomposition, and the number of operations required are O(2*(2N _N)³/3),^{29, 30} contributing to ∼25% of total computation time (Fig. 2) taken for calculation of J.

Fig. 1

Flowchart presenting the important steps involved in reconstructing optical properties in diffuse optical tomography along with associated computation time (in percent of total time per iteration). The steps that used GPU-based computations are indicated with a double arrow.

Fig. 2

Flowchart showing the important computations for calculation of the Jacobian (J). The computation time in terms of percentage of total time taken in calculation of J at each step is also indicated. Note that the calculation of J accounts for 75% of the total time taken in completing one iteration of diffuse optical image reconstruction (also indicated in Fig. 1).

As the Rytov approximation is used in the work, the frequency domain data becomes, y = [ln(A); θ], where ln(A) is the natural logarithm of amplitude (A) and θ is the phase of the frequency domain signal, making the J a real valued matrix (dimension of 2N _M×2N _N, where N _M is the number of measurements).^{4, 7} The procedure involved in the calculation of real valued J using complex Φ is given in Ref. 25. Note that the calculation of J [and G(μ)] consumes ∼75% of total computation time taken for completing a single iteration (see Fig. 1). To obtain an update of the optical properties (Δμ), the Levenberg–Marquardt (LM) minimization^{33, 34} is used and the objective function for this minimization scheme is given as^{33, 34}

Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Omega = \,\parallel \!y - G(\mu)\!\parallel^{2}. \end{equation}\end{document} $$\Omega =\parallel y-G\left(\mu \right){\parallel }^2.$$

The objective of the LM minimization scheme is to match y with G(μ) in the least-squares sense, by changing μ. The details of LM minimization scheme are discussed in Ref. 7. The update equation (for getting Δμ) for the LM minimization becomes⁷

Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} [{J^{\rm T} J + \lambda I}]\,\Delta \mu = J^{\rm T} [{y -G(\mu)}], \end{equation}\end{document} $$\left[J^{\mathrm{T}}J+\lambda I\right]\Delta \mu =J^{\mathrm{T}}\left[y-G\left(\mu \right)\right],$$

where J ^T represents the transposed J and I is the identity matrix. λ is the regularization parameter, chosen empirically (starts at 10 multiplied by the maximum of the diagonal values of J ^T J and reduced by a factor of 10^0.25 at every subsequent iteration). Computing Δμ with the use of J and G(μ) typically consumes 20% of total time in any given iteration (see Fig. 1). LU factorization is used in solving Eq. 7 to obtain Δμ, with the number of operations of order O(2*(2N _N)³/3).^{29, 30} The procedure for calculation of J, G(μ), and subsequently Δμ is repeated until the relative difference in the objective function [Ω, Eq. 6] does not improve by >2% (the same is indicated in Fig. 1). Because 95% of total time per iteration is spent on calculation of J, G(μ), and Δμ, the emphasis of this work is to accelerate these computations, which involve either matrix-matrix multiplications [calculation of P ⁻¹ K in Eq. 4 and J ^T J in Eq. 7] and solving the system of linear equations [Eqs. 4, 7], using GPUs.

3.

Results and Discussion

There are many potential CUDA (GPU)-based numerical linear algebra packages that could be used to accelerate the diffuse optical image reconstruction procedure, some of these along with their features are listed in Table 1. With the need of matrix-matrix computations and system solve leads to the choice of CULA and Jacket from the ones listed in Table 1. Because CULA basic¹⁶ has an advantage of being open-source, the same is used in the presented work. Note that CULA is also part of Jacket.²⁰ Also, from the Table 1, it is evident that the complex sparse arithmetic is not supported in CUDA.

The difference (in percent) in the L2-norms of the data-model misfit (δ) and estimate optical properties (μ) with respect to iteration number between single- and double-precision computations using CPU are plotted in Fig. 3 for the test object described in Sec. 2.3. It is evident from Fig. 3 that the maximum difference in terms of percentage in either δ or estimated μ is <0.005%, asserting the fact that single precision is sufficient for carrying out the diffuse optical tomographic image reconstruction computations. Because Rytov approximated data are used in this work [i.e., using ln(A) rather than A] the single-precision calculations were adequate. The main advantages of single- over double-precision computations are twofold. First, given the limitation on the available GPU/CPU memory, diffuse optical tomographic image reconstruction could be performed on meshes twice as large. Second, the basic version of CULA (open source), which has only single-precision capability, could be used to perform the GPU computations.

Fig. 3

Difference (in percent) in the L2 norm between the single- and double-precision computation over the iterations in (a) data-model misfit (δ) and (b) estimated optical properties (μ).

The optical images obtained with integration of CULA- based functions (Table 3 with the instances given in Sec. 2.2) into the diffuse optical image reconstruction procedure using the test case described earlier were identical (visually and numerically) compared to images obtained using single-precision CPU-based implementation. The observed difference (in percent) in the L2-norms of the data-model misfit (δ) and estimate optical properties (μ) between the GPU implementation (with CULA basic) and traditional CPU implementation (with double precision) are same as the observed values plotted in Fig. 3.

A 2-D test case reconstruction results along with the target distribution using the three strategies (CPU-Full, CPU-Sparse, and GPU-Full, discussed in Sec. 2.3) are given in Fig. 4. The number of iterations to converge to a solution for each strategy is equal to 10. Obtained μ_a and μ^′ _s distributions were correlated to find a similarity measure among the three strategies, resulting in a correlation coefficient of 1 (with in the limits of single precision), asserting that the reconstructed results are identical.

Fig. 4

Two-dimensional reconstruction results obtained with 1% noisy data using three strategies—CPU-Full, CPU-Sparse, and GPU-Full— discussed in this work. The target 2-D distribution is given in column 1 for (a) μ_a images and (b) μ^′ _s images. The number of nodes along with computational time taken for each strategy in a given iteration is given in Table 4.

The time taken per iteration [start to end, including data transfer between CPU (host) and GPU (device) and vice versa] on GPU (NVIDIA Tesla C 1060) and CPU (Intel Xeon E 5410) for different mesh sizes both in two and three dimensions were plotted in Fig. 5 for the three strategies (GPU-Full, CPU-Full, and CPU-Sparse, discussed in Section. 2.3). In the 2-D case, it is clear from Fig. 5 that the computation time per iteration increases with increasing N _N for all three strategies, with CPU-Full taking the maximum time. Among the three strategies, GPU-Full took the shortest computation time. As stated earlier, because sparse complex arithmetic is not supported in NVIDIA CUDA, the GPU-Sparse (GPU-based computation using sparse representation of matrices for instances 1–4, described Section 2.2) strategy was not attempted. Comparing the three strategies in this 2-D case, the speedup (acceleration) with the use of GPU was up to 7.5 comparing CPU-Full to GPU-Full. The maximum speedup has been lowered to 2 in comparing CPU-sparse to GPU-Full. For the 3-D case [Fig. 5], a similar trend in the 2-D case was observed with the computation time per iteration being higher than 2-D case for the same N _N. In the 3-D case, the speedup (acceleration) with the use of GPU was up to 40 in comparing CPU-Full to GPU-Full, and was lowered to 7 in comparing CPU-sparse to GPU-Full.

Fig. 5

Plots showing the computation time taken for each strategy [indicated in the legend of (a)] per single iteration versus mesh size (node number, N _N) for (a) the 2-D case and (b) 3-D case. The strategies are explained at the end of Section 2.3. The specifications of CPU and GPU are given in Table 2.

The main emphasis of this work is to prove that the GPU computing offers a considerable speedup compared to CPU computing and the comparison has been only performed by the time taken per iteration. The diffuse optical image reconstruction typically requires about 10 iterations, and for a typical three-dimensional problem (Table 4), the least total time taken by the CPU (CPU-sparse) to converge to a solution is 290 min, and 48 min for GPU. The gain in terms of total computation time is significant in nature (comparing ∼5 h to 0.75 h), making the GPU computing very attractive to be used in the real time, especially in multiwavelength cases (the problem gets scaled by the number of wavelengths), where the aim is to obtain functional images.

Table 4

Comparison of computation time (in seconds) taken at each instance (as listed in Section 2.2) between 2-D and 3-D meshes of similar size using both CPU (Intel Xeon E 5410 with 4 GB memory) and GPU (NVIDIA Tesla C 1060) in a given iteration using the three strategies discussed in Section 2.3. The last row indicates the total time taken for the iteration. Because J is full matrix, the CPU-Sparse calculations are not applicable (NA) in instances 5 and 6.

Even though it is not possible to attempt sparse complex arithmetic on NVIDIA CUDA-enable GPUs, we have attempted to test the capabilities of parallelizing sparse real computations (with the aid of Cusp, refer to Table 1) using a steady-state diffusion equation (where K and Φ are real) to estimate μ_a using intensity-based measurements in both 2-D and 3-D cases. Even though the GPU is able to provide speedup (up to 2) for calculation of Φ or Φ* in solving sparse system of linear equations, the overall computation time taken per iteration from start to end (including the overhead of transferring sparse matrices from the CPU memory to GPU memory) is increased (at least by 70%) in comparison of GPU to CPU in these sparse cases (similar to observed trend in Refs. 40, 41). The sparse cases have large overhead for GPU calculations in comparison to the CPU ones because the GPU compute kernels are not fully optimized. This overhead can become negligable when the problem size (number of nonzero entries) becomes bigger (typically >40,000). For fully bandwidth-optmized FEM meshes, the number of nonzero entries for the GPU-limited mesh (with N _N = 13,377) is ∼30,000, making the GPU-sparse calculations not attractive compared to CPU-saprse calculations. Thus, only CPU-sparse implementation is attempted in this work. Also, note that sparse computations on currently available GPUs were mainly centered around matrix-vector computations (encountered in iterative solvers to deduce x in Ax = b), which are not as straightforward to implement and integrate into the existing domain-specific applications, such as optical image reconstruction.^{42, 43} Hopefully, the next generation of GPUs will be able to perform sparse computations with similar ease as current real computations and support sparse complex computations (including matrix-matrix computations). Although it is not attempted for the fair comparison of the CPU and GPU computations, the optimal computations will utilize both CPU (for sparse matrix computations) and GPU (for full matrix computations).

Table 4 gives the computational time taken in a given iteration for the mesh of similar size (N _N ≈ 9200) both in 2-D and 3-D cases. Although the GPU calculations took similar time between 2-D and 3-D calculations, there is a factor-of-20 difference between the 2-D and 3-D calculations with CPU-Full strategy, especially in solving for Φ in Eq. 4. This is primarily due to the higher numerical condition number of 3-D mass matrix [K inEq. 3], which takes more iterations to converge to a solution [Φ in Eq. 4] compared to their counterparts in 2-D cases for the same size K.^{8, 44} This higher condition number in the 3-D case could result from the nonuniform nature of 3-D mesh. The tetrahedral elements are formed using 3-D delaunay, which makes a uniform mesh more complex to generate,⁴⁵ as compared to the 2-D mesh. Also, the condition number of the mass matrix (K) depends on the number of source fibers and their placement. Within this work, where a ring type of data acquisition is performed, the regions far away from the fibers have negligible contribution as opposed to regions close to source/detector fibers. Therefore, the 3-D problem has a higher condition number, requiring more iterations. Hence, this leads to longer computation time, as compared to the 2-D problem, for decomposing P ^–1 K and then solving for Φ in Eq. 4. Because the decomposition primarily involves matrix-vector computations, with GPU offering massive parallelization, results in similar computation time for the same size matrices both in the 2-D and 3-D cases (Table 4) with GPU giving peak performance in the 3-D case. Note that with J being a full matrix, computations for the last two instances (5 and 6) in Table 4 were not applicable to use sparse computations. These computations were carried out using CPU-Full strategy, and the same computation time taken by the CPU-Full strategy were taken into account in calculating the total time taken per iteration for CPU-sparse strategy.

The imaging domains considered in this work are circular for the 2-D case and cylindrical for the 3-D case, with a tumor-to-background contrast being 2. But the trends observed in terms of computational speedup observed between CPU and GPU computing should be similar to any test object (varying tumor-to-background contrast, size, and position of tumor). It is very well known that the contrast recovery is highly dependent on the spatial location of the tumor;³⁸ the images obtained either using single or double precision will be identical (shown in Fig. 3).

As mentioned earlier, due to the available onboard memory (4 GB), the FEM mesh that could be used in GPU calculations has to be <13,377 nodes on a NVIDIA Tesla C1060 GPU card. With the advent of newer GPU architectures and generations, this limit could soon disappear. The current Fermi-based NVIDIA Tesla C2050 GPU card offers 6 GB of onboard memory (with 448 CUDA cores),⁴⁶ which can enable GPU calculation for mesh with the number of nodes as 16,384. Note that the choice of FEM mesh is dependent on the size of the imaging domain, number of measurements, signal-to-noise ratio of the data, and the required stability of the forward solver.⁴⁷ Discussion of the same is beyond the scope the presented work.

Because modern day CPUs have multiple cores, it is possible to take advantage of the parallel computing power of these multicore CPUs. Recent work by Borsic ⁴⁸ for a similar problem (electrical impedance tomography) like diffuse optical image reconstruction, the achievable speedup using a dual quad-core Intel Xeon processor was 7.6. But this required rewriting/optimizing the routines used in the computing of Jacobian and forward problem, where as the GPU calculations performed in this work used already preexisiting routines. The speedup obtained using a multicore processor with out the optimization is only 4.6.⁴⁸ Also, the typical cost of dual quad-core processor is at least twice expensive compared to the top of the line GPU used in this work (Table 2).

It is important to note that traditional parallel computing uses large number of processors connected in parallel to form a cluster, with an overhead of cost and demanding large footprint in the lab settings. These bulky parallel computer systems are not very attractive in the clinical settings. The GPU boards are comparatively low cost and fit into a traditional desktop computer, resulting in an desktop parallel computer (with a very small footprint), could become very attractive in the clinic to perform image reconstruction/analysis tasks. Usage of these GPU cards to accelerate FE meshing using multimodal imaging data that could be used in the NIR imaging studies is currently being explored. The CULA-based mex programs used in this work along with necessary documentation for installation and compiling are provided as open source.³⁶

4.

Conclusions

Because general purpose computing using GPUs is becoming attractive in many areas of medical image reconstruction/processing, this work demonstrated that with the use of GPUs the speedups for a single iteration (start to end) of diffuse optical tomographic image reconstruction could be up to 40 compared to traditional CPU calculations. Also, the GPU computations are carried out using an open-source package (CULA) and appropriate code wrappers were used to integrate these routines into the existing diffuse optical image reconstruction package (NIRFAST). The current limitation to use GPU on a daily basis for diffuse optical tomographic image reconstruction is the available onboard memory (currently it is limited to 4 GB), which puts a restriction on mesh size that could be used. Moreover, the sparse complex calculations are not supported in GPU computing, making it less attractive for use in FE-based models. Because there is tremendous growth in the capabilities of GPU computing, this limitation might vanish, making GPU computing capabilities similar to those of CPU in applications that require more memory and sparse computations. More importantly, along with the cost effectiveness, these GPU cards are integrated with traditional desktop computers, giving massive parallel computing power at the desktop. Specific to diffuse optical tomography, GPU computing holds the promise of making the image reconstruction procedure match the time line, with the data acquisition, leading to dynamic NIR imaging viable in the clinic.