Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2339

Simulation study on compressive laminar optical tomography for cardiac action potential propagation TAKUMI H ARADA , 1,* N AOKI TOMII , 1 S HOTA M ANAGO, 2 E TSUKO K OBAYASHI , 2 AND I CHIRO S AKUMA 1,3 1 Department

of Bioengineering, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan 2 Department of Precision Engineering, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan 3 Medical Device Development and Regulation Research Center, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan * [email protected]

Abstract: To measure the activity of tissue at the microscopic level, laminar optical tomography (LOT), which is a microscopic form of diffuse optical tomography, has been developed. However, obtaining sufficient recording speed to determine rapidly changing dynamic activity remains major challenges. For a high frame rate of the reconstructed data, we here propose a new LOT method using compressed sensing theory, called compressive laminar optical tomography (CLOT), in which novel digital micromirror device-based illumination and data reduction in a single reconstruction are applied. In the simulation experiments, the reconstructed volumetric images of the action potentials that were acquired from 5 measured images with random pattern featured a wave border at least to a depth of 2.5 mm. Consequently, it was shown that CLOT has potential for over 200 fps required for the cardiac electrophysiological phenomena. c 2017 Optical Society of America

OCIS codes: (170.6960) Tomography; (170.3880) Medical and biological imaging.

References and links 1. E. M. C. Hillman, and S. A. Burgess, "Sub-millimeter resolution 3D optical imaging of living tissue using laminar optical tomography,"Laser & Photon. Rev. 3(1-2), 159–179 (2009). 2. E. M. C. Hillman, O. Bernus, E. Pease, M. B. Bouchard, and A. Pertsov, "Depth-resolved optical imaging of transmural electrical propagation in perfused heart," Opt. Express 15(26), 17827–17841 (2007). 3. Y. Chen, S. Yuan, J. Wierwille, R. Naphas, Q. Li, T. R. Blackwell, P. T. Winnard Jr. V. Raman and K. Glunde, "Integrated optical coherence tomography (OCT) and fluorescence laminar optical tomography (FLOT)," IEEE J. Sel. Topics Quantum Electron. 16(4), 755–766 (2010). 4. B. Yuan, S. A. Burgess,A. Iranmahboob, M. B. Bouchard, N. Lehrer, C. Bordier, and E. M. Hillman, "A system for high-resolution depth-resolved optical imaging of fluorescence and absorption contrast," Rev. Sci. Instrum. 80(4), 043706 (2009). 5. E. M. C. Hillman, A. Devor, M. B. Bouchard, A. K. Dunn, G.W. Krauss, J. Skoch, B. J. Bacskai, A. M. Dale, and D. A. Boas, "Depth-resolved optical imaging and microscopy of vascular compartment dynamics during somatosensory stimulation," Neuroimage 35(1), 89–104 (2007). 6. S. A. Burgess, D. Ratner, B. R. Chen, and E. M. C. Hillman, "Fiber-optic and articulating arm implementations of laminar optical tomography for clinical applications," Biomed. Opt. Express 1(3), 780–790 (2010). 7. M. Azimipour, M. Sheikhzadeh, R. Baumgartner, P. K. Cullen, F. J. Helmstetter, W. Chang, and R. Pashaie, "Fluorescence laminar optical tomography for brain imaging: system implementation and performance evaluation," J. Biomed Opt. 22(1), 016003–016003 (2017). 8. M. S. Ozturk, C. Chen, R. Ji, L. Zhao, B. B. Nguyen, J. P. Fisher, Y. Chen, and X. Intes, "Mesoscopic Fluorescence Molecular Tomography for Evaluating Engineered Tissues," Ann. Biomed. Eng. 44(3), 667–679 (2016). 9. Q. Tang, J. Lin, V. Tsytsarev, R. S. Erzurumlu, Y. Liu, and Y. Chen, "Review of mesoscopic optical tomography for depth-resolved imaging of hemodynamic changes and neural activities," Neurophotonics. 4(1), 011009–011009 (2017). 10. Q. Tang, V. Tsytsarev, J. Lin, Y. Liu, C. Chen, R. S. Erzurumlu, Y. Chen, "3D mesoscopic imaging of neural connections in sensory and motor cortices," in Proc. Int. Photon. Conf. 2016 IEEE. ME 4.5 (2016). 11. M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F. Kelly, and R. G. Baranluk, "Single-pixel imaging via compressive sampling." IEEE Signal Process. Mag. 25(2), 83 (2008).

#285108 Journal © 2017

https://doi.org/10.1364/BOE.8.002339 Received 19 Jan 2017; revised 25 Feb 2017; accepted 7 Mar 2017; published 24 Mar 2017

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2340

12. V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, "Compressive fluorescence microscopy for biological and hyperspectral imaging," Proc Natl Acad Sci USA 109(26), 1679–1687 (2012). 13. Y. Wu, P. Ye, I. O. Mirza, G. R. Arce, and D. W. Prather, "Experimental demonstration of an Optical-Sectioning Compressive Sensing Microscope (CSM)," Opt. Express 18, 24565–24578 (2010). 14. X. Yuan, and S. Pang, "Structured illumination temporal compressive microscopy," Biomed. Opt. Express 7(3), 746–758 (2016). 15. P. Ye, J. L. Paredes, Y. Wu, C. Chen, G. R. Arce and D. W. Prather, "Compressive confocal microscopy: 3D reconstruction algorithms," Proc. SPIE 7210, 72100G (2009). 16. L. Zhu, W. Zhang, D. Elnatan, and B. Huang, "Faster STORM using compressed sensing," Nat Methods 9(7), 721–723 (2012). 17. Y. August, C. Vachman, Y. Rivenson, and A. Stern, "Compressive hyperspectral imaging by random separable projections in both the spatial and the spectral domains," Appl. Opt. 52(10), D46–D54 (2013). 18. G. R. Arce, D. J. Brady, L. Carin, H. Arguello, and D. S. Kittle, "Compressive coded aperture spectral imaging," IEEE Signal Process. Mag. 31(1), 105–115 (2014). 19. A. Veeraraghavan, D. Reddy, and R. Raskar, "Coded strobing photography: Compressive sensing of high speed periodic videos," IEEE Trans. Pattern Anal. Mach. Intell. 33(4), 671–686 (2011). 20. D. Craven, B. McGinley, L. Kilmartin, M. Glavin, and E. Jones, "Compressed sensing for bioelectric signals: a review," IEEE J. Biomed. Health Inform. 19(2), 529–540 (2015). 21. M. Süzen, A. Giannoula, and T. Durduran, "Compressed sensing in diffuse optical tomography," Opt. Express 18, 23676–23690 (2010). 22. F. Yang, M. S. Ozturk, L. Zhao, W. Cong, G. Wang, and X. Intes, "High-resolution mesoscopic fluorescence molecular tomography based on compressive sensing," IEEE Trans. Biomed. Eng. 62(1), 248–255 (2015). 23. C. Antzelevitch, J. Fish, "Electrical heterogeneity within the ventricular wall," Basic Res. Cardiol. 96(6), 517–527 (2001). 24. C. Antzelevitch, "Heterogeneity and cardiac arrhythmias: an overview," Heart Rhythm 4(7), 964–972 (2007). 25. AG. Kléber, Y. Rudy, "Basic mechanisms of cardiac impulse propagation and associated arrhythmias," Physiol. Rev. 84(2), 431–488 (2004). 26. R. H. Clayton, E. A. Zhuchkova, and A. V. Panfilov, "Phase singularities and filaments: simplifying complexity in computational models of ventricular fibrillation," Prog. Biophys. Mol. Biol. 90(1), 378–398 (2006). 27. Z. Qu, J. Kil, F. Xie, A. Garfinkel, and J. N. Weiss, "Scroll wave dynamics in a three-dimensional cardiac tissue model: roles of restitution, thickness, and fiber rotation," Biophys. J. 78(6), 2761–2775 (2000). 28. O. Berenfeld,M. Yamazaki, D. Filgueiras-Rama, J. Kalifa, "Surface and intramural reentrant patterns during atrial fibrillation in the sheep," Methods Inf. Med. 53(4), 314 (2014). 29. J. I. Laughner, F. S. Ng, M. S. Sulkin, R. M. Arthur, and I. R. Efimov, "Processing and analysis of cardiac optical mapping data obtained with potentiometric dyes," Am. J. Physiol. Heart. Circ. Physiol. 303(7), H753–H765 (2012). 30. A. K. Iranmahboob and E. M. C. Hillman, "Diffusion vs. Monte Carlo for Image Reconstruction in Mesoscopic Volumes," OSA Biomedical Topical Meetings, OSA Technical Digest (2008). 31. C. Zhu, and Q. Liu, "Review of Monte Carlo modeling of light transport in tissues," J. Biomed. Opt. Circ. Physiol. 18(5), 050902–050902 (2013). 32. R. G. Baraniuk, "Compressive sensing," IEEE Signal Process. Mag. 24(4), 118–120 (2007). 33. E. J. Candés, "Compressive sampling," in Proc. Int. Congr. Math. 3, 1433–1452 (2006). 34. E. J. Candés and M. B. Wakin, "An introduction to compressive sampling," IEEE Signal Process. Mag. 25(2), 21–30 (2013). 35. Y. C. Eldar and G. Kutyniok, Compressed Sensing Theory and Applications (Cambridge, 2012). 36. H. Niu, Z. J. Lin, F. Tian, S. Dhamne, H. Liu, "Comprehensive investigation of three-dimensional diffuse optical tomography with depth compensation algorithm." J. Biomed Opt. 15, 046005 (2010). 37. A. Mahajan, Y. Shiferaw, D. Sato, A. Baher, R. Olcese, L. Xie, M. Yang, P. Chen, J. G. Restrepo, A. Karma, A. Garfinkel, Z. Qu, and J. N. Weiss, "A Rabbit Ventricular Action Potential Model Replicating Cardiac Dynamics at Rapid Heart Rates," Biophys. J. 94(2), 392–410 (2008). 38. Q. Fang, and D. A. Boas, "Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units," Opt. Express 17(22), 20178–20190 (2009). 39. D. A. Boas, C. Pitris, and N. Ramanujam, Handbook of Biomedical Optics (CRC Press, 2011), pp. 92. 40. J. M. Bioucas-Dias, and M. A. T. Figueiredo, "A New TwIST: Two-Step Iterative Shrinkage/Thresholding Algorithms for Image Restoration." IEEE Trans. Image Process. 16(12), 2992–3004 (2008). 41. J. I. Laughner, F. S. Ng, M. S. Sulkin, R. M. Arthur, and I. R. Efimov, "Processing and analysis of cardiac optical mapping data obtained with potentiometric dyes," Am. J. Physiol. Heart. Circ. Physiol. 303, H753–H765 (2012). 42. R. C. Aster, B. Borchers, and C. H. Thurber, Parameter Estimation and Inverse Problems, Second Edition (Academic Press, 2011). 43. M. N. Do, and M. Vetterli, "The contourlet transform: an efficient directional multiresolution image representation," IEEE Trans. Image Process. 14(12), 2091–2106 (2005). 44. Y. M. Lu, and M. N. Do, "Multidimensional directional filter banks and surfacelets," IEEE Trans. Image Process. 16(4), 918–931 (2007). 45. J. Ma, and G. Plonka, "The curvelet transform," IEEE Signal Process. Mag. 27(2), 118–133 (2010).

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2341

1.

Introduction

Fluorescence tomography is a useful tool for investigating volumetric distributions of tumors and other features in three-dimensional biological tissue. Laminar optical tomography (LOT), focusing on the microscopic activity of tissue, can be applied to the optical mapping of electrophysiological action potential propagation by utilizing voltage-sensitive dyes [1, 2]. Conventional tomography techniques, which use some photomultipliers and two galvano mirrors (GMs), usually focus on a static biological feature, such as a tumor that does not change rapidly over time. The conventional LOT method uses raster scanning, which involves measurement at each point in the tissue and has a low recording speed, mainly due to the mechanical speed limitation of the GMs and the low intensity of the fluorescence to be detected [1, 3–6]. In this approach, there is a trade-off between the recording speed and the spatial resolution of the system. Therefore, the system cannot measure rapid biological activities without a reduction in spatial resolution. However, phenomena such as electrical signal propagation and ion channel activities of neurons and cardiac muscles occur at high speed. To investigate the activity of the whole brain or cardiac muscles, the development of high-frame-rate fluorescence tomography with a wide and deep field of view compared with that in the conventional microscopic technique, such as a confocal microscope or a two-photon microscope, is anticipated. In recent years, LOT has mainly been used to obtain three-dimensional images of optogenetic protein expression in the cortical regions of the brain tissue [7], those of cellular distribution in thick tissue engineering constructs [8], and those of hemodynamic changes and neural activities [9]. In order to achieve fast three-dimensional mesoscopic imaging, a time-resolved angled fluorescence laminar optical tomography (aFLOT) system in combination with voltage-sensitive dye imaging has been used [9, 10]. This system works based on line illumination and has a charge coupled device (CCD) camera, and this can measure the area of the line illumination at 200 fps. However, it requires the repetitive position change of the object by the motor stage to complete the entire field of view (FOV) scanning. Therefore, measuring the aperiodic fast activity of tissue still remains difficult. CS techniques have recently become popular for single-detector systems and rapid video recording techniques and have been used in many imaging applications [11–20]. The fundamental concept of CS involves reconstructing the original information from fewer measurements. This does not involve conventional data compression after data acquisition, but instead resembles the previous method of data compression before data acquisition. The CS enables sub-Nyquist measurements and a drastic reduction in the number of measurements. This advantage enables an increase in the frame rate when using certain types of imaging system. To date, several studies of CS-based tomography have been reported [21, 22]. For example, illumination and detection patterns have been designed to improve the reconstruction quality of the sparse fluorophore distribution using a CS framework. However, they are not intended to increase the frame rate. In other applications, coded masks are usually used for CS-based imaging systems because of the ease with which they can be modulated [11–14, 17, 18]. However, these masks decrease the detected light intensity and this reduction of intensity critically affects the reconstruction accuracy because diffuse optical tomography requires measurement of the distribution of light of low intensity on the surface. Thus, the mask strategy is not suitable for tomography applications. The structure of the heart wall is composed of three layers: the epicardium, myocardium, and endocardium. Many studies have reported that the heterogeneity of the electrophysiological properties of the heart wall structures affects the development of arrhythmia [23, 24]. Cardiac excitation during arrhythmia propagates in an intricate three-dimensional pattern through the cardiac muscle tissue [25–28]. However, a detailed understanding of these propagation mechanisms remains elusive. The elucidation of such propagation mechanisms is of substantial clinical significance to medical cardiac treatments such as cardiac ablation and drug therapy. In the field of cardiac electrophysiology, there is an optical mapping technique known as 2D imaging for measuring membrane potential using voltage-sensitive dye [29]. The measurement system

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2342

of optical mapping comprises an excitation LED or laser, a long-pass filter, and a high-speed camera. The excitation light is irradiated onto the cardiac muscle tissue widely and uniformly, and the dye, which is near the tissue surface, is excited. This fluorescence imaging technique can image only the surface of the tissue, but imaging deeper regions by the optical mapping technique remains difficult. To improve the frame rate of LOT in comparison with the conventional tomography technique, we propose compressive laminar optical tomography (CLOT), which is a form of LOT using patterned illumination based on CS. A high frame rate can be achieved via following two components of the proposed method. (1) two-dimensional (2D) detection of an optical system with a fast change in the illumination patterns for fast data acquisition in the FOV; and (2) a reduction in the number of the recorded images in a single reconstruction of the target by CS framework. In addition, the new aspect of the proposed method is the active modulation of the mathematical model for measurement by using the illumination patterns, which is aimed at efficient application of CS. In this paper, we present simulation experiments investigating the transformation performance for the cardiac action potential distribution, the reconstruction error upon changing the random excitation pattern, the influence of the number of measurements, and the robustness against noise compared with those in the conventional raster scanning system. Finally, we discuss the disadvantages and advantages of CLOT. 2.

System design

To date, the sensitivity and resolution of CCD sensors have been improved, and the image sensors can be used to measure the spatial information of the weak distribution of the fluorescence in the FOV at each pixel simultaneously. This data-recording method can achieve a higher recording speed than the conventional raster scanning method. In addition, a digital micromirror device (DMD) can produce the various illumination patterns at a high speed and is now widely used in projectors. In the proposed method, a camera has been used to detect the fluorescence and the DMD has been used to produce the changes of the illumination patterns. A schematic diagram of this system is shown in Fig. 1. The mathematical model is described in the following section.

Fig. 1. Schematic diagram of optical system

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2343

2.1.

System model

The intensity of the k-th illumination pattern Pk (u,v) generated by the DMD is expressed as follows: Pk (u, v) = |he x (u, v) ∗ |Mo |uk (Mo u, Mo v)| 2 (1) where (u, v) is the coordinate system of the object surface, he x is the excitation point spread function (PSF) of the objective lens, and uk is the geometric image of the k-th illumination mirror pattern. Mo is a magnification factor of the object lens. After excitation, the fluorescence from the specimen at time t, Ak (u, v, t), is represented as follows: Ak (u, v, t) = Pk (u, v) • O(u, v, t) (2) $ =  (x, y, z, t)c(x, y, z)τ(x, y, z, t)H (x, y, z)E(u − x, v − y, z)dxdydz (3) O, H, and E are the intensity response of the specimen, the 3D distribution of excitation light in the tissue generated by the pattern illumination, and the 3D distribution of emission light in the tissue, respectively.  (x, y, z, t), c(x,y,z), and τ(x, y, z, t) are the extinction coefficient of the dye, the concentration of the dye, and the quantum yield of the dye’s fluorescence, respectively. The k-th fluorescent image Ik (n, m, t) on the detector is expressed as follows: Ik (n, m, t) = |hem (n, m)| 2 ∗ | Ak (n/Md , m/Md , t)/|Md || 2

(4)

where (n,m) is the coordinate system of the detector and hem is the emission PSF of the imaging lens. Md is the magnification factor of the imaging lens. The k-th measurement image f k (a, b, t) at detector position (a,b) is expressed as follows: " n m f k (a, b, t) = dndm Ik (n, m, t)rect(a − ,b− ) + ρab (5) ∆det ∆det ρab is the system noise at (a,b). The measurement matrix Y can be written simply as follows: ∆F = J∆V + ρ

(6)

where J is the sensitivity matrix, which represents the spatial sensitivity of the fluorescence intensity measurement at each detector with the illumination pattern and can be obtained by multiplying H and E at every voxel. The Monte Carlo method has been widely used to model photon propagation in scattering samples in many applications of LOT [1, 7–9]. The radiative transport equation (RTE) or the diffusion approximation to the RTE can be used in place of the Monte Carlo method; however, it is known that the diffusion approximation does not provide accurate results in the case of the small source-detector separations used in LOT in general [30]. In the Monte Carlo method, a large number of photons need to be simulated in order to obtain the stochastic accuracy. Various methods that can speed up the Monte Carlo simulations have been developed, and the simulated results can be obtained by using the developed mehods [31]. The membrane potential V is the dominant factor in the intensity response of the specimen; therefore, V is represented such that sV ≈ cτ in Eq. (3) where s is a scaling constant. ρ is the noise matrix. The details of this linear function (6) are described elsewhere [2]. 3.

CS theory and the reconstruction method

Followed by the process of recording the images, the three-dimensional data of the target are reconstructed by using CS framework. The authors would like to emphasize that the real-time calculation of the reconstruction is not the aim of this study, but the reduction in the number of recorded images for a single reconstruction after the rapid image recording is the main objective. This section describes the proposed reconstruction method using CS theory [32-35].

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2344

3.1.

CS theory

The aim of CS is to determine the desired coefficients from fewer measurements with sub-Nyquist sampling. The expression of the linear equation is y = Φx

(7)

where y is the M ×1 measurement vector, Φ is the M ×N basis matrix, and x is the N×1 coefficient vector. This coefficient vector x is sparse if at most K(N) entries are nonzero. To reconstruct the coefficient vector from few measurements, the ` 1 minimization argmink xk1

subject to y = Φx

(8)

has been used using the size of measurements M ≥ O(Klog(N/K )). For instance, the number of measurements can be reduced if the number of nonzero entries is small. Then, the number of nonzero entries affects the compression performance and is an important factor to apply CS to the actual measuring system model. 3.2.

Reconstruction method

In the previous section, the system’s mathematical model is described as a linear equation, which is applicable for CS. However, the membrane potential V is not a sparse matrix, so a transformation to make the matrix sparse is needed. To achieve this, Eq. (6) can be written as follows: F = JV = JΨα (9) where Ψ is the transformation matrix to make the matrix sparse and α is the sparse coefficient matrix. To take JΨ as the new basis matrix, the system’s mathematical model becomes a linear equation with the sparse matrix, which is applicable for CS. First, the sparse signal α is recovered from the measurements Y . Second, the obtained sparse signal is transformed to the desired membrane potential V by inverse transformation of Ψ. In the case of the pattern illumination, the intensity of the excitation light is nonuniform at each voxel and information of each voxel is unevenly obtained. To compensate this ununiformity, the sensitivity matrix J ∈ R M × N is normalized at each voxel in the reconstruction process. Matrix W ∈ R N × N is difined by W = [w w . . . w] (10) w = [c1 , c2 , . . . , c N ]0 cj =

1 max1≤i ≤ M ai j

(11) (12)

where ai j is the element in the i–th row and j–th column of the matrix J. In the inverse reconstruction process, the sparse coefficient matrix is recovered from the following equation: F = JW Ψα

(13)

This is a sensitivity compensation method for the pattern illumination method and a 2–dimensional expansion of the depth compensation algorithm for the conventional DOT described in [36]. The previous research uses a raster scanning and the sensitivity matrix J is fixed [22]. On the other hand, in the optical set-up of the proposed method, J could be modulated by the change of the illumination pattern on every recorded frame. This modulation enables applying CS to the proposed imaging system efficiently. As a result, a reduction of the number of the recorded images in a single reconstruction of the target could be achieved.

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2345

4. 4.1.

Simulation method Measurement objects

Six representative 3D models of the propagation of rabbit action potentials were prepared (Fig. 2). The rabbit action potentials were generated from [37]. The scale of the model was 20 × 20 × 5 mm (32 × 32 × 8 voxels). Each side of the voxel was 625µm each side. These models represent perpendicular propagation, horizontal propagation with different wavefront angles, and the ectopic focus at the center of the tissue. The wave angles of the horizontal propagation models are 90◦ , 45◦ , and -45◦ . The conduction velocity was assumed to be 30 cm/s. 4.2.

Sparsity and energy evaluation

The discrete cosine transform (DCT, DCT-II) and the discrete wavelet transform (DWT) were performed on the ∆V of the models. In the DWT, we prepared Haar, Coiflet, and Daubechies wavelets with a multilevel decomposition structure. The transformations were performed in the x and y directions and the z direction separately and the same level decomposition was applied in the x and y directions. When there are poor sparsity coefficients with a few large coefficients, it can be considered that the recovery signal from the coefficients is relatively good, even though the number of nonzero factors is large in the coefficient matrix. This situation can be referred to as a substantial sparse matrix. Next, we investigated the energy of the coefficients to determine each transformation property. Here, the energy of the coefficients En is defined as the summation of the squares of P each coefficient α i , which is expressed as En = α 2i . After this sparsity and energy evaluation, we chose three transformation methods that had the highest sparsity, the highest energy, or intermediate levels of each, which were used for the following reconstruction evaluation to examine the effects of the sparsity and energy on the reconstruction.

Vertical 1

Vertical 2

Horizontal 1

Horizontal 2

Horizontal 3

Sphere

Fig. 2. 3D object models

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2346

4.3.

Measurement condition

The wavelength of the excitation light is 532 nm and that of the emission light is 610 nm for Di-4-ANEPPS, the commonly used voltage-sensitive dye. The sensitivity matrix J was made using Monte Carlo simulation [38] using the scattering coefficient, absorption coefficient, and scattering anisotropy of the rat muscle (550 nm: µs =8.82 mm −1 , µa =0.166 mm −1 , g=0.909, 600 nm: µs =8.33 mm −1 , µa =0.095 mm −1 , g=0.926 [39]). In the Monte Carlo simulation, 107 photons were simulated at each irradiation and detection point. The simulations were performed on a machine equipped with a Intel Core i7-3960X CPU, 32GB RAM, and a NVIDIA GeForce GTX 780 Ti GPU. It took around 3 minutes for 107 photons. The illumination pattern in FOV was the random matrix pattern (32 × 32 pixels), in which only the number of irradiation points is fixed and each irradiation laser point has a diameter of 1 mm. The random matrix is widely known to be applicable to CS and the reconstruction can be performed precisely. Every image was recorded with the different illumination pattern. The number of recorded images used in a single reconstruction was varied from 5 to 20. The PSF was assumed to be ideal, the magnification factors of all lenses were 1, and the detector had 32 × 32 pixels in this simulation. The FOV was a 20-mm square and the entire images of the illumination patterns were used for reconstruction. In noise addition evaluation, for comparison with the conventional system, the raster scanning method with Tikhonov regularization was also performed under the same conditions. This raster scanning method involves recording using seven detectors at each pixel representing the tissue surface. The seven detectors were arranged at 625-µm intervals from the laser irradiation point to a distance of 3.75 mm away. 4.4.

Evaluation of illumination patterns and the number of measurements

Considering the principle of LOT, the excitation lasers interfere with each other, so information on tissue deep in the heart cannot be obtained if the excitation lasers are located close together. However, data on the surface of the tissue cannot be obtained efficiently from the one measured image if each excitation laser is located sparsely. The optimal number of irradiation points can be determined based on the density of laser irradiation on the tissue. To investigate this issue, we altered the number of irradiation points from 10 to 90 in increments of 10 and from 100 to 1000 points in increments of 100. The number of recorded images used to reconstruct the object influences the extent of error in the reconstruction. It is considered that the reconstruction error decreases if the number of measurements increases because the amount of information obtained becomes larger. To investigate this issue and evaluate the minimum number of measurements, reconstruction was performed with 5, 10, and 20 measurements. 4.5.

Reconstruction and error

There are many algorithms of ` 1 based optimization, for example, Matching Pursuit (MP), Orthogonal Matching Pursuit (OMP), Stagewise Orthogonal Matching Pursuit (StOMP), Iterative Hard Thresholding algorithm (IHT), Compressive Sampling Matching Pursuit (CoSaMP), and so on [35]. The property of the algorithm must be considered when choosing the algorithm for the practical application. To recover sparse signals, the two-step iterative shrinkage/thresholding algorithm for linear inverse problems (TwIST) solver was used here because the algorithm exhibits fast convergence rate for ill conditioned and ill-posed problems [40]. In the reconstruction process, our main interest is the shape of the wave of membrane potential. Then, after reconstruction by both the proposed and the conventional methods, each reconstructed voxel was determined to have high or low potential, representing a depolarized or a repolarized state. This processing comprised two parts: a binarizing part and a morphology operation part. In the binarizing part, the reconstructed voxel is binarized to high or low using a threshold set as the mean of the maximum and minimum potential of each voxel. In the morphological operation

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2347

part, a pixel at a certain depth is set to high if all adjacent pixels are high and it is set to low if all adjacent pixels are low. This part reduces the appreciable error because, electrophysiologically, it is unlikely that only one voxel (625µm cube) would conflicts with its adjacent voxels. 4.6.

Noise addition

To compare the proposed and conventional methods in terms of robustness against noise, random noise based on the intensity as system noise was added to the recorded signal and reconstruction was performed. The noise intensity was set to 10%, 20%, and 30% of the recorded intensity. 5. 5.1.

Results Sparsity and energy

The number of nonzero entries out of all 8192 voxels in each transformation is shown in Table 1. In this table, the first number is the mean number of nonzero entries at all objects and the second number in parentheses is the standard deviation of each object. Here, the first number after the wavelet name represents the decomposition level of DWT for the x and y directions and the second number represents that for the z direction. DCT had 8143 (25) nonzero entries. As shown in Table 1, the Haar 1.3 transformation had the lowest number of nonzero entries. The Daubechies and Coiflet transformations had high values compared with the Haar transformation. The energy of 50 coefficients from the beginning in descending order is shown in Table 2. The energy of DCT was 90% (14). Haar 3.3 had the highest energy among all of the transformations shown in Table 2. Considering these results and the difference in the transformation method between DCT and DWT, we decided to evaluate the reconstruction performance using Haar 1.3, Haar 3.3, and DCT to compare how the sparsity and the energy influence reconstruction among these different transformations. 5.2.

Reconstruction error

The number of error voxels in each model with 10 measurements was determined, as shown in Fig. 3. The performances of each transformation were almost same for each model. The number

Table 1. Number of nonzero entries in DWT. Standard deviations are in parentheses.

Level 1.1 1.2 1.3 2.1 2.2 2.3 3.1 3.2 3.3 4.1 4.2 4.3 5.1 5.2 5.3

Haar 2321 (997) 2428 (1057) 1539 (1320) 2294 (1475) 2181 (1669) 2118 (1683) 2377 (1449) 2315 (1636) 2023 (1746) 2361 (1521) 2278 (1727) 2119 (1777) 2383 (1520) 2313 (1720) 2084 (1804)

Db2 8084 (81) 8092 (70) 8091 (73) 8162 (31) 8163 (28) 8163 (28) 8165 (15) 8166 (14) 8165 (14) 8168 (15) 8168 (14) 8167 (14) 8167 (18) 8168 (17) 8165 (19)

Db3 8188 (7) 8188 (7) 8188 (7) 8192 (0) 8192 (0) 8192 (0) 8192 (1) 8192 (1) 8192 (1) 8192 (0) 8192 (0) 8192 (0) 8192 (0) 8192 (0) 8192 (0)

Db4 8172 (14) 8172 (15) 8172 (15) 8188 (3) 8187 (4) 8188 (4) 8187 (3) 8187 (3) 8187 (3) 8183 (5) 8183 (5) 8183 (5) 8186 (4) 8186 (4) 8186 (4)

Coif 8190 (2) 8190 (2) 8190 (2) 8190 (2) 8190 (2) 8190 (2) 8188 (4) 8185 (9) 8185 (9) 8188 (3) 8188 (3) 8188 (3) 8189 (3) 8189 (3) 8189 (3)

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2348

Table 2. Energy of 50 coefficients in descending order in DWT. Standard deviations are in parentheses.

Level 1.1 1.2 1.3 2.1 2.2 2.3 3.1 3.2 3.3 4.1 4.2 4.3 5.1 5.2 5.3

Haar [%] 5 (0) 10 (0) 20 (1) 19 (1) 38 (2) 75 (4) 74 (4) 94 (7) 95 (9) 92 (8) 93 (11) 93 (14) 88 (13) 89 (16) 90 (19)

Db2 [%] 5 (0) 15 (1) 11 (1) 22 (1) 65 (4) 46 (4) 67 (5) 80 (9) 75 (10) 50 (8) 68 (15) 59 (15) 42 (8) 63 (16) 53 (16)

Db3 [%] 5 (0) 10 (0) 12 (1) 27 (1) 47 (5) 53 (6) 59 (8) 62 (16) 64 (18) 38 (17) 42 (22) 45 (23) 39 (17) 43 (23) 46 (23)

Db4 [%] 5 (0) 11 (0) 10 (0) 34 (2) 54 (3) 56 (4) 38 (21) 42 (23) 43 (24) 40 (22) 44 (24) 45 (25) 40 (22) 44 (24) 45 (25)

Coif [%] 5 (0) 13 (0) 15 (1) 23 (1) 57 (3) 64 (7) 74 (8) 79 (11) 76 (18) 52 (22) 57 (25) 57 (28) 56 (24) 60 (26) 60 (29)

of irradiation points at which the error was minimized differed for each model. In the vertical conduction models (Vertical 1 and 2), there was a tendency for the error to be decreased with a decrease in the number of irradiation points. In the other models, there was an opposite tendency. To investigate each depth layer and the influence of the number of measurements further, the number of error voxles at each depth layer (Z = 1–8 [voxel]) with a change in the number of measurements is shown in Fig. 4. In the vertical conduction models, the errors that occur at the depth layers that are around the boundary layer between the high and low potentials (Z = 3–5) was high. In the horizontal conduction models (Horizontal 1–3 and Sphere), the error up to a depth of 5 voxel decreased with an increase in the number of measurements. Figure 5(a) shows the number of error voxels in the reconstruction by using all the models at each depth layer in conventional method (raster scanning method that uses Tikhonov regularization). From the figure, it is clear that the error increased from a depth of 5 voxels. However, the error was < 5% even at a depth of 8 voxels. Figure 5(b) shows the error by using the Horizontal 1–3, and Sphere models at each depth layer in the proposed method. There was a tendency for the error to be increased with an increase in the depth. The error was < 15% even at a depth of 4 voxles. To clarify the morphological detail, Figs. 6 and 7 show the reconstructed slice images in the Horizontal 3 and Sphere models as examples. Information about the deeper regions of the object was obtained with an increase in the number of recorded images. DCT and DWT produced almost the same results. The images at deep layers had reversed distribution of the action potential. 5.3.

Noise

Reconstruction error in each depth at different noise levels is shown in Fig. 8, A. The conventional method provided a better result at a weak level of noise of 10%; however, DCT, Haar 1.3 and 3.3 provided better results upto a depth of four voxels (Z = 4) at all noise levels. The increase of erroneous voxels with increasing noise level was smaller in the proposed method than in the conventional method. The change of the reconstructed error with changing the number of measurements is shown in Fig. 8, B. The number of errorneous voxels was decreased by 2.57 % (210.7 voxels) at each additional recorded image.

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2349

Vertical 1

Vertical 2

Horizontal 1

Horizontal 2

Horizontal 3

Sphere

Fig. 3. Reconstruction error in each model using 10 measurements. Purple ∗, brown +, and green 4 indicate the results of DCT, Haar 1.3, and Haar 3.3 .

For morphological evaluation, the reconstructed slice images in the Horizontal 3 and Sphere models are shown in Figs. 9 and 10. The conventional tomography method could obtain information about deeper layers compared with the proposed method, as shown in Fig. 9. However, the images from the conventional method were noisy at shallow layers. In the case of 30% noise addition as shown in Fig. 10, the voltage distribution could not be recognized using the conventional method at all depth layers. The images reconstructed with the proposed method were little noisy; however, the images upto a depth of five voxel (Z = 5) were clear, in which the voltage distribution could be recognized. 5.4.

Intensity of the sensitivity matrix

The intensity of the sensitivity matrix is one of the factors that determine the depth that could be reconstructed. In order to investigate this factor further, Fig. 11 shows the normalized intensity of the sensitivity matrix in the case that the detection at which the excitation light was irradiated, which means the relative detected intensity of the fluorophore at each depth. The intensity decreased by ' 10 −1 at each depth layer (e.g., Z = 5 [voxel]: 10 −5 ; Z = 8: 10 −8 ). 5.5.

Reduction of the number of measurements

The number of measurements required in each method cannot simply be compared because all the measurement methods are different from each other. However, the reduction in the number of measurements for reconstruction is an important factor that is required to measure the fast activity. Figure 12 shows the number of measurements required per second for the frame rate of the reconstructed data. While using the conventional method, for the 32×32 pixel ROI, the number of measurements required for a single reconstruction is 1024, and this number is 51.2, 102.4, and 204.8 times larger than those required in the method proposed in this study, that is 20, 10, and 5 measurements, respectively.

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2350

Vertical 1 (Z = 1–4)

Vertical 1 (Z = 5–8)

Vertical 2 (Z = 1–4)

Vertical 2 (Z = 5–8)

Horizontal 1 (Z = 1–4) Horizontal 1 (Z = 5–8) Horizontal 2 (Z = 1–4) Horizontal 2 (Z = 5–8)

Horizontal 3 (Z = 1–4) Horizontal 3 (Z = 5–8)

Sphere (Z = 1–4)

Sphere (Z = 5–8)

Fig. 4. Reconstruction error in each model at each depth layer (DCT using 400 irradiation points).

Conventional method

Proposed method

Fig. 5. Reconstruction error at each depth. (a) in the conventional method with all models, (b) proposed method with 20 measurements (Horizontal 1–3, and Sphere).

6. 6.1.

Discussion Sparsity in the transformations

As shown in Table 1, the Haar 1.3 transformation showed the best results in terms of sparsity. The cause of this was considered to be that Daubechies and Coiflet have complicated waveform wavelet filters and the many coefficients after these transformations are generated. By contrast, the Haar transformation in DWT is the results of mean and subtraction of values using a simple rectangle filter. The cause of the sparsest result occurring at 3-level decomposition in the z direction in the case of Haar transformation was that the low spatial frequency in the z direction was dominant

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2351

Fig. 6. Reconstructed slice images in Horizontal 3 using 400 irradiation points. The layers in the xy plane are arranged from left to right in order of increasing depth. a: 3D object model, b: slice images of the object model, c: 3D reconstructed data of f by DCT, d: slice images reconstructed from 5 measurement by DCT, e: from 10 measurements by DCT, f: from 20 measurements by DCT, g: 3D reconstructed data of j by Haar 1.3, h: slice images reconstructed from 5 measurement by Haar 1.3, i: from 10 measurements by Haar 3.3, j: from 20 measurements by Haar 1.3, k: 3D reconstructed data of n by Haar 3.3, l: slice images reconstructed from 5 measurement by Haar 3.3, m: from 10 measurements by Haar 3.3, n: from 20 measurements by Haar 3.3.

in each model. Compressibility of high level decomposition is high if the object has almost low frequency. The prepared models have low spatial frequency in the z direction; however, some models have high spatial frequency in the x and y directions. Consequently, the 1-level decomposition in the x and y directions and 3-level decomposition in the z direction at Haar showed the best result regarding sparsity. 6.2.

Energy in the transformations

Haar 3.3 had the highest energy in all transformations, as shown in Table 2. The high-pass coefficients are not decomposed in DWT. Hence, the high-pass coefficients with high values

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2352

Fig. 7. Reconstructed slice images in Sphere using 400 irradiation points. The layers in the xy plane are arranged from left to right in order of increasing depth. a: 3D object model, b: slice images of the object model, c: 3D reconstructed data of f by DCT, d: slice images reconstructed from 5 measurement by DCT, e: from 10 measurements by DCT, f: from 20 measurements by DCT, g: 3D reconstructed data of j by Haar 1.3, h: slice images reconstructed from 5 measurement by Haar 1.3, i: from 10 measurements by Haar 3.3, j: from 20 measurements by Haar 1.3, k: 3D reconstructed data of n by Haar 3.3, l: slice images reconstructed from 5 measurement by Haar 3.3, m: from 10 measurements by Haar 3.3, n: from 20 measurements by Haar 3.3.

remain when the low-pass coefficients cannot be decomposed efficiently from a certain level decomposition. Consequently, the 3-level decomposition in the x and y directions and the 3-level decomposition in the z direction at Haar had the highest energy. 6.3.

Effects of the sparsity and energy

As shown in Fig. 3, the reconstruction errors were almost the same between DCT, Haar 1.3 and Haar 3.3. In addition, Figs. 6 and 7 show that the depth that could be reconstructed did not differ in these transformations and the number of measurements was related to the ability to obtain information on deeper layers. Hence, it is considered that the sparsity and energy in

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2353

Fig. 8. Reconstruction error stacked bar graph at each depth. A: from 20 measurements at different noise levels. B: by Haar 1.3 with 10% noise.

these transformations had only little effect on the reconstruction of the target by CS framework in the 2D detection system and on the physical condition, and scattering and attenuation of the light mainly affected the reconstruction in this case. In comparison with the fluorophore that exists locally, the distributions of the action potentials mainly have the direct current (DC) component of the spatial frequency, and more number of coefficients were required to represent the distribution;.therefore, these distributions cannot be represented as the sufficient sparse matrices by using DCT and DWT (maximum data compression ratio = 18.8%; Table 1), and the reconstruction results obtained from DCT, Haar 1.3, and Haar 3.3 were similar in the prepared distributions. In order to make the use of CS framework more effective, the other transformation method is required for the distribution of the action potentials. 6.4.

Illumination pattern

In Vertical 1 and 2, there was a trend for fewer irradiation points to be associated with smaller error, as shown in Fig. 3. Vertical 1 and 2 had the same voltage distribution in each depth layer and the main coefficients representing these models were particularly related to the depth information. Thus, the xy plane sampling with the high-density pattern was not effective in these models and the use of fewer irradiation points to obtain information on deeper layers resulted in smaller error. On the other hand, there was a trend for more irradiation points to be associated with smaller error in the other models. The other models mainly had the xy plane structure and the structure changed smoothly as the depth increased. The main coefficients to representing them were focused on the xy plane information and the high-density pattern for obtaining the xy plane information was effective. For these reasons, the optimal number of irradiation points depends on the object structure and the pattern of around 400 irradiation points (1 points/mm2 ) was useful for these models if the object to be measured had an unknown structure. 6.5.

Reconstruction error and the number of measurements

While using the vertical conduction models, the error that occur at the depth layers that are around boundary layer between the high and low potentials (Z = 3–5, Fig. 4. The reconstructed data were mainly influenced by the signals of the shallow layers because the fluorescent signals at the shallow layers were higher than those at the deep layers. As a result, the signals at the shallow layers were overestimated, and the boundary layer shifted to a deeper layer. Hence, the error at a depth of Z = 3–5 voxels was high. In the other models, the error that occur up to a depth of Z = 4 voxel decreased more than that

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2354

Fig. 9. Reconstructed slice images in Horizontal 3 and Sphere with 10% noise. The layers in the xy plane are arranged from left to right in order of increasing depth. a: Horizontal 3 reconstructed by conventional method, b: Sphere reconstructed by conventional method, c: Horizontal 3 reconstructed from 20 measurements using 400 points by DCT, d: Sphere reconstructed from 20 measurements using 400 points by DCT, e: Horizontal 3 reconstructed from 20 measurements using 400 points by Haar 1.3, f: Sphere reconstructed from 20 measurements using 400 points by Haar 1.3, g: Horizontal 3 reconstructed from 20 measurements using 400 points by Haar 3.3, h: Sphere reconstructed from 20 measurements using 400 points by Haar 3.3

at a depth of Z = 5 voxel with an increase in the number of measurements. Considering the light scattering and attenuation, to get information on deep layers is difficult; however; the information can be obtained by increasing the number of measurements required. In order to reconstruct the information on a depth of Z = 5 voxels or deeper layers, the number of measurements required need to be increased. The depth limitaion for the reconstruction of the target within the reasonable number of measurements was Z = 5 voxels.. In general, the errors that occur up to a depth of Z = 2 voxels in the vertical conduction models and the errors that occur up to a depth of Z = 4 voxels in the horizontal conduction models were 30%, 20% and 10% with 5, 10, and 20 measurements. The images at the deep layers had reversed distribution of the action potential (Figs. 6 and 7). Simply, the transformations, which were performed in the present study, represented the structure of the object as the overlapping waves. The coefficients representing the shallow layers were

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2355

Fig. 10. Reconstructed slice images in Horizontal 3 and Sphere with 30% noise. The layers in the xy plane are arranged from left to right in order of increasing depth. a: Horizontal 3 reconstructed by conventional method, b: Sphere reconstructed by conventional method, c: Horizontal 3 reconstructed from 20 measurements using 400 points by DCT, d: Sphere reconstructed from 20 measurements using 400 points by DCT, e: Horizontal 3 reconstructed from 20 measurements using 400 points by Haar 1.3, f: Sphere reconstructed from 20 measurements using 400 points by Haar 1.3, g: Horizontal 3 reconstructed from 20 measurements using 400 points by Haar 3.3, h: Sphere reconstructed from 20 measurements using 400 points by Haar 3.3

dominantly obtained and the coefficients at the deep layers were hardly obtained due to light scattering and absorption. The overlapping waves have both the ridge and the valley. Therefore, the reconstructed data at the deep layers have reversed values in contrast to the data at the shallow layers when the coefficients representing the deep layers are not almost obtained. As shown in Figs. 6 and 7, the slice images at least up to a depth of 2.5 mm (Z = 4 [voxels]) could be reconstructed from five recorded images. When comparing the reconstructed images from the five images with those from 10 images, the depth that could be reconstructed did not differ. However, one more deeper layer could be reconstructed from 20 images because the obtained information on the deep layers increased. The depth limitation that could be reconstructed within the reasonable number of measurements was Z = 5 voxels even though the sensitivity compensation method was used. As shown in Fig. 5, the error in the conventional method was smaller than that in the proposed

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2356

Fig. 11. Normalized intensity of the sensitivity matrix at each depth.

Fig. 12. Required number of measurements per second for the frame rate of the reconstructed data (32×32 pixel ROI).

method. In Tikhonov regularization, a smoothness constraint is imposed on the reconstructed data during the calculation of the inverse problem [42]. The object models, which mainly had the DC component of the spatial frequency, had smooth distributions and the error was small in the conventional method. The error increased from a depth of Z = 5 voxels. It was considered that the information on the depth of Z = 5 or deeper layers were hardly obtained. However, the error was underestimated because the object models had the smooth distributions to the depth direction from a depth of Z = 5 voxel and the smooth distributions were reconstructed in spite of not obtaining the information on the deep layers. As shown in Fig. 11, the intensity of the detected signal at a depth of Z = 5 voxel was hundred-thousandth of the intensity at a depth of Z = 1 voxel. The small intensity at a depth of Z = 6 or deeper layers would be underestimated and reduced to zero in the calculation process. 6.6.

Noise robustness

The conventional tomography method using Tikhonov regularization is known to be particularly susceptible to noise. The conventional method could reconstruct the data at all depth layers to reveal the rough distribution of the voltage; however, this method resulted in compartments of noise, as shown in Figs. 9 and 10. The data reconstructed by the novel method proposed in this

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2357

paper were associated with less noise. It can be asserted that the newly proposed method measures the rough features of objects and the small coefficients representing noise are suppressed in the process of ` 1 minimization. In the optical mapping at 500 fps, the noise level was around 10% empirically. When the recording speed is increased, the noise level will be increased. Hence, the proposed method is very useful at a high noise level. On the other hand, the conventional method has small error due to the small noise level and can reconstruct the deep layers; however, the distribution of the action potentials in the reconstructed data cannot be recognized at the high noise level. Considering the practical application for the biological tissue, the usable application in the conventional method is limited by the noise level of the detected signal and the number of measurements required in the proposed method should be determined to decrease the noise for recognizing the distribution (in reference to Fig. 8, B). 6.7.

Reconstructable depth limitation

The ` 1 minimization leads to overfitting in the high-intensity region. The information on the shallow layers, therefore, was overfitted and on the deep layers was underfitted in this simulation. If the data compression ratio increases greatly from that in the other transformation methods, the depth, which could be reconstructed, may increase. However, the limitation of depth reconstruction was strictly determined by the light propagation property and the intensity of the fluorescence signal, and the depth limitation was Z = 5 voxel in this simulation condition. The detected signal at the irradiation point of the excitation light is the largest; therefore, the detectable depth limitation can be determined in terms of the intensity of the detected signal. Considering that the bit depth of the commercially available camera was around 16-bits, it could be considered that around 5 voxels (normalized intensity = 10 −5 ; Fig. 11) is the depth limitation that can be detected to use the dynamic range of the camera effectively at the following conditions: excitation = 550 nm, and emission = 600 nm. This hardware aspect must be taken into account for the application of the proposed method. 6.8.

Frame rate

The frame rate of the reconstructed data of conventional LOT was 66.7 fps within the FOV (3.7 × 3.7 mm; spatial resolution = 370 µm) in [2]. The FOV of the present system is a 20-mm square, which was an area 29.2 times larger than that in the previous study. In order to measure the area of a 20-mm square by using the conventional LOT, it takes approximately 438 ms to obtain one frame (approximately 2.28 fps). Assuming that the typical conduction velocity of the electrical propagation of the rabbit is 30 cm/s, the wave moves 131 mm in 438 ms. Then, each movement of 1.5 mm should be measured at 200 fps. Taking into consideration that the recording speed of optical mapping for cardiac action potential propagation is usually 500–1000 fps, the recording speed of the proposed method could be 1000 fps with the change of the illumination pattern on every frame (the frame rate of the reconstructed data is 200 fps when five recorded images are used). Thus, it is asserted that the proposed system has the potential to achieve 200 fps or an even higher frame rate of the reconstructed data. In the viewpoint of the number of measurements, as shown in Fig. 12, the required number of measurements per second in the conventional method is 204.8 times larger than that in the proposed method with five measurements. Assuming that the typical velocity of the electrical propagation is 30 cm/s, the wave moves 1.5 mm during the recording process of five images at 1000 fps. Hence, the limitation of the motion tracking of the conduction is 1.5 mm in the proposed method. However, further reduction in the number of measurements could be achieved if the object models are represented as more sparse matrice by the other effective transformation methods. In addition, the proposed system has an advantage, that is, its spatial resolution, compared with

Vol. 8, No. 4 | 1 Apr 2017 | BIOMEDICAL OPTICS EXPRESS 2358

the conventional method. For 1024 × 1024 camera pixels, the maximum spatial resolution is 19.5 µm when the FOV is a 20-mm square. In order to achieve the same spatial resolution in the conventional method, the measurement takes at least 4161 times longer than that in the proposed method at a recording speed of 500 fps. Hence, in the case of the high-resolution measurement, the proposed method has a major advantage in terms its recording speed. 7.

Conclusion

We investigated the fundamental performance of CLOT using simulation and clarified that the transformation performance, the optimal number of irradiation points, the reconstruction performance with different numbers of measurements, and the robustness against noise. There are certainly some drawbackss in the proposed method compared with the conventional method. The high power excitation laser is required because the expanded excitation light must be used to reflect on DMD in the proposed method. The reconstructable depth is determined by the dynamic range of the camera. It is difficult to precisely reconstruct the boundary between the high and low potentials in the vertical conduction models because the signals at the shallow layers are overestimated. The error is larger than that in the conventional method under the no noise condition. However, there are some major advantages in the proposed method. This method exhibits noise robustness in the high noise level, which might happen in practical biological experiments and the reduction in the number of measurements for a single reconstruction in the proposed method leads to the high frame rate of the reconstructed data, which is very important for fast activity measurements. The findings, which investigated in this paper, demonstrated that the reduction of the number of the recorded images in a single reconstruction of the target could be achieved by the illumination patterns based on CS. As a result, the voltage distribution to a depth of 2.5 mm (Z = 4) can be reconstructed clearly by the proposed method using 5 patterns of random irradiation from 400 points (1 points/mm2 ). This reflects strong potential for the high-frame-rate of the reconstructed data of cardiac action potential propagation at a notable level of over 200 fps. As future work, the effective transfomation method for the distribution of the cardiac action potentials and the pattern optimization should be investigated with the goal of minimizing the number of measurements. In recent years, various transformation techniques have been developed, and some transformation methods may be more useful than the transformation techniques proposed in the present study [43–45]. Also there should be a more effective pattern to obtain coefficients representing the model and decreasing the number of the pattern sequence would lead to more high frame rate. Additionally, near-infrared voltage-sensitive dye is useful for obtaining information on deep layers. The voltage-sensitive dye Di-4-ANEPPS is the gold standard for the recording of action potentials. However, the limited light propagation when using green excitation light is a critical component of this method for depth reconstruction. Further studies are required to precisely measure the action potential propagation in the practical application. Funding Japan Society for the Promotion of Science (JSPS) (15H01801)