Open Access
Add to BibSonomy
Abstract
A reflection matrix based optical coherence tomography (OCT) is recently proposed and expected to extend the imaging-depth limit twice. However, the imaging depth and hence the image quality heavily depend on the number of primary singular values considered for image reconstruction. To this regard, we propose a method based on correlation between image pairs reconstructed from different number of singular values and corresponding remainders. The obtained correlation curve and another feature curve fetched from the former are then fed to a long short-term memory (LSTM) network classifier to identify the optimized number of primary singular values for image reconstruction. Simulated targets with different combinations of filling fraction and signal-to-noise ratio (SNR) are reconstructed by the developed method as well as two current adopted methods for comparison. The results demonstrate that the proposed method is robust to recover the image with satisfactory similarity close to the reference one. To our knowledge, this is the first comprehensive study on the optimized number of the primary singular values considered for image reconstruction in reflection matrix based OCT.
© 2022 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
The scattering of light in tissue is not only a source of information, but also a limiting factor for deep imaging [1]. With an increased depth penetration, the exponential attenuation of single-scattered photons and the increment of multiple-scattered photons make traditional focusing and imaging techniques based on the Born approximation invalid, and severely limit the imaging depth in tissue for most optical imaging techniques [2]. In order to minimize the multiple-scattered contribution and enhance the single-scattered contribution to optical imaging, a method based on reflection matrix measurement and successive iterative time reversal processing for optical coherence tomography (OCT) is recently proposed [3–5]. The imaging-depth limit can be extended by at least 2 times in comparison with conventional OCT [3], where wide-field detection followed by digital-pinhole filtering and iterative time reversal processing is the underlying mechanism. The iterative time reversal processing is mathematically equivalent to singular value decomposition of the filtered matrix [6,7]. Singular vectors with singular values above a chosen threshold are picked out and reshaped into two matrices according to the scanning trajectory of the focused beam. The image is then reconstructed by coherence summation of a chosen number of primary components, each of them given by the Hadamard product between two reshaped matrices weighted by corresponding singular value. It is noted that insufficient components only provide partial information, while excessive components introduce a lot of noise corresponding to weak singular values. Therefore, it is crucial to determine an optimized number of primary singular values of the filtered matrix for high fidelity image recovery.
To determine the primary singular values of the filtered matrix and hence the eigenstates associated with the target, criterion based on the image quality or the distribution of singular values has been considered. Based on the image quality assessed by the standard deviation (STD) of the reconstructed image, the STD-based method is feasible to determine the best reconstruction and corresponding primary singular values [3]. Another adopted approach is based on the distribution of singular values. Singular values beyond the superior bound of the empirical distribution corresponding to multiple-scattered contributions (described by the random matrix theory) were considered as single-scattered contributions and picked out for reconstruction [8–10]. However, the ideal assumptions and the computational complexity of the empirical distribution method make it only a theoretical illustration and difficult to apply in practice. Alternatively, the maximum ratio method was adopted to determine the primary singular values by delimiting at the point of the maximum ratio between consecutive singular values to pick out the significantly large ones [11,12].
Although successful demonstrations have been reported for some simple targets, image reconstruction based on current methods is not always feasible, especially for targets under complex situations. Therefore, in this paper, we propose a method based on correlation between image pairs reconstructed from chosen primary singular values and the remainders respectively. A correlation curve is thus obtained by calculation of the correlation coefficient versus the number of primary singular values. The key point corresponding to the optimized number is exactly located near the demarcation point of the two zones of the correlation curve. To favor the finding of the demarcation point on the correlation curve, another feature curve for variation-trend characterization is further fetched from the correlation curve. These two curves are then fed to a long short-term memory (LSTM) network classifier to identify the optimized number of primary singular values to recover the target. Image reconstructions are done on targets with different combinations of filling fraction and signal-to-noise ratio (SNR) by the developed method and two current adopted methods for comparison.
2. Method
The matrix-based OCT approach relies on the measurement of the reflection matrix from the scattering medium in real space (point-to-point basis). One focal spot in the target plane corresponds to one frame of 2D optical field in the sensor plane. Each frame is reshaped into a vector for the matrix rearrangement in the order of the scanning trajectory and is recombined into one reflection matrix for post-analysis. According to the fact that multiple-scattered contribution is randomly distributed, and that single-scattered contribution emerges as the diagonal and closed-diagonal elements of the reflection matrix, a digital pinhole is thus used to filter out the off-diagonal elements (mainly associated with multiple-scattered contribution) and transfers the original reflection matrix into a filtered matrix. The residual multiple-scattered contribution in the filtered matrix can be further suppressed to the greatest extent by iterative time reversal processing. Mathematically, it is realized by singular value decomposition of the filtered matrix R expressed by $R = U\sum {V^\dagger }$. Here U and V are unitary matrices whose columns correspond to the input and output singular vectors ${U_i}$ and ${V_i}$, and the superscript $\dagger $ stands for transpose conjugate. $\sum$ is a diagonal matrix consisting of the real positive singular values ${\sigma _i}$ in a decreasing order of i$\{{i \in [{1,N} ]} \}$, ${\sigma _1} > {\sigma _2} > \cdots > {\sigma _N}$, where N is the size of $\sum$. The most important yet challenging task is to determine the optimized number M (M < N) of the primary singular values. Once M is determined, the image can be reconstructed by
where $\overline {{U_i}}$ and $\overline {{V_i}}$ are wavefronts reshaped from vectors ${U_i}$ and ${V_i}$, and the reshaping order is matched with the scanning trajectory of the focused beam.Figure 1 outlines the proposed method to determine the optimized number M for image reconstruction. A typical reflection matrix, a digital-pinhole filtered matrix and its singular value decomposition by two unitary matrices and one diagonal matrix consisting of singular values are shown in Fig. 1(a). With singular values divided into the primary singular values with a number of n and the remainders with a number of N-n, corresponding image pairs (${I_n}$, ${I_{N - n}}$) can be reconstructed. The correlation coefficient C(n) between image pairs is then calculated using the zero-mean normalized cross-correlation (ZNCC) method as previously reported [13].
Fig. 1. Outline of the proposed method. (a) Reflection matrix, filtered matrix and its singular value decomposition. (b) The correlation curve (the blue curve) of correlation coefficients versus the number n, where correlation coefficients marked in red circular points and corresponding singular value decompositions and image pairs for two typical cases are presented. The feature curve (the yellow curve) is fetched from the correlation curve. (c) Diagram depicting key point identification by the classifier and successive image reconstruction based on the optimized M.
It is noticed from Fig. 1(b) that the correlation curve can be divided into a rapid changed zone and a gradual changed zone. In the first rapid changed zone, with the addition of more primary singular values, the correlation coefficient first increases and then decreases with a narrow peak due to the fact that more primary singular values participate in the reconstruction of ${I_n}$ while less primary singular values participate in the reconstruction of ${I_{N - n}}$. In the second gradual changed zone, with subsequent addition of low singular values (mainly associated with the noisy background) to the reconstruction of ${I_n}$ while almost full of noise in ${I_{N - n}}$, the correlation coefficient changes relatively slowly. An optimized point M should exist between above-mentioned two zones, where the correlation between ${I_n}$ and ${I_{N - n}}$ is very low, and the optimized number M of primary singular values are considered to lead to the best reconstruction of ${I_n}$. In order to manifest the variation-trend characteristics of the correlation curve, a feature curve is introduced based on the correlation curve as follows. The slope at each point on the correlation curve is first calculated by
where ${C_0}$ is the correlation coefficient at the local peak or valley point ${n_0}$ closest to n, and C(n) is the correlation coefficient at position n. To suppress local fluctuations on fetching local characteristics of the curve, a sliding window with a width of $(2w + 1)$ is further adopted to the calculation of the variance ${V_C}$ of the correlation coefficient and the average ${\mu _k}$ of the slope given by where ${\mu _C}$ is the average of the correlation coefficients within the window. Thus, the feature parameter at n can be defined by Based on Eq. (6), the feature curve corresponding to the correlation curve can be obtained. As shown in the yellow curve in Fig. 1(b), above-mentioned two zones in the normalized feature curve exhibit distinguishable variation-trend characteristics, which favor the finding of the optimized M.The optimized point M can be identified by setting a threshold for feature curves to separate the two zones. However, drastic local changes in correlation curves may cause fluctuations in feature curves, resulting in failure to correctly locate the optimized point. LSTM networks are good at sequence classification and key point identification, which is achieved by finding and exploiting long range dependencies in the sequential data [14]. Thanks to the memory cell, which can maintain its state over time to store historical information, and nonlinear gating units, which regulate the information flow into and out of the cell. The changes at each point on the correlation curve are not independent but related to each other, then the correlation curve is exactly what LSTM networks are good at dealing with. As a result, we have constructed and trained a LSTM classifier to solve the problem of key point identification of the correlation coefficient sequence. The diagram of the identification of the key point by the classifier and successive image reconstruction based on the optimized M is depicted in Fig. 1(c). The classifier is composed of two LSTM layers and two dropout layers. Both the correlation curve and the feature curve are fed to the classifier as the input. Simulated targets with different combinations of filling fractions (ratio between the number of pixels of the target and the total number of pixels) and SNRs for training. Once the classifier is trained successfully, the key point can be identified based on the results of the two-class classification. With the identified key point as the estimate of the optimized number M of primary singular values, the target image can be reconstructed using Eq. (1).
3. Results
When multiple scattering dominates, the matrix coefficients become random and uncorrelated [15]. Indeed, provided that there is no ballistic wave in the field at the output of the sample, the multiple-scattered matrix is a random matrix and is applicable to the random matrix theory [9]. To demonstrate the distribution of the multiple-scattered contributions in a strongly scattering medium is random, the multiple-scattered light from a target sample behind a scattering layer of human epidermal tissue is simulated by the Tracepro software. The specific parameters used in the Tracepro software are as follows: the total number of simulated rays is 50 million, the type of bulk scattering is human skin epidermis (refractive index: n = 1.335, scattering anisotropy parameter: g = 0.72, scattering coefficient: ${\mu _s} = 16.5m{m^{ - 1}}$, scattering mean free path ${l_s} = 60.6\mu m$, anisotropic scattering model: Henyey-Greenstein model). It is found that, with a thin scattering layer (D = 0.4 mm), the multiple-scattered light is not random but somehow concentrated due to limited scattering events (Fig. 2(a)). However, with a thick scattering layer (D = 0.7 mm), increased scattering events randomize the distribution of multiple-scattered light (Fig. 2(b)). Accordingly, as the matrix-based OCT aims to deep imaging of biological tissue, it is reasonable to take the multiple-scattered contributions of the reflection matrix as random variables in simulations. The multiple-scattered light and single-scattered light from the target sample behind the scattering layer are simulated in the Tracepro software as well and used to calculate SNR corresponding to the random model of multiple-scattered contributions
Fig. 2. Distribution of multiple scattered light with different thickness of the scattering layer and the horizontal and vertical profile curves at the center. The simulated scattering medium is human epidermal tissue. (a) The thickness of the scattering layer D = 0.4 mm. (b) The thickness of the scattering layer D = 0.7 mm.
To demonstrate the dependence of the number of the primary singular values on the target to be recovered, different samples are simulated as the targets with the results given in Fig. 3. Under a constant SNR, the number of eigenstates associated with the target increases with the enlarging filling fraction. Therefore, the peaks of the first zones of the correlation curves become broadened due to the increased number of the primary singular values. Under the same filling fraction, the target signal (single-scattered contribution) is gradually submerged in noise with decreasing SNR. Hence shrunken peaks of the first zones of the correlation curves are noticeable due to the gradually decreased number of the primary singular values. With fixed filling fraction and SNR, the distribution of the target of different letters (Z, J, U) makes little effect on the peak width of the first zones of the correlation curves. These results demonstrate that the peak width of the first zone of the correlation curve depends on filling fraction and SNR rather than distribution of the target, and is reasonable to be taken as a representation of the optimized number of the primary singular values considered for reconstruction.
Fig. 3. Dependence of the optimized number of the primary singular values on the target. Details of the area indicated by dotted rectangles are shown in the left zoom-in views. Normalized correlation curves corresponding to a) different filling fractions under a constant SNR, b) different SNRs (SNR 1=-36.3 dB, SNR 2=-37.1 dB, SNR 3=-37.9 dB, SNR4=-38.8 dB) under the same filling fraction, and c) different letters (Z, J, U) with fixed filling fraction and SNR.
To confirm that the LSTM classifier can better estimate the number of primary singular values compared to setting a threshold for feature curves directly, and both the correlation curve and the feature curve are required as the input of the LSTM network, similarity evaluations on image reconstructions of targets with different combinations of filling fraction and SNR are carried out. The similarity between the reconstructed image and the reference one is calculated by the ZNCC method [16,17], from which the highest similarity score can be obtained as the best similarity ${S_{best}}$, and the similarity score ${S_{method}}$ of the image reconstructed by the method can be obtained as well. Then the similarity evaluation is represented by the bias from the best similarity: ${S_{best}} - {S_{method}}$. Figure 4 demonstrates the results of bias of similarity on different targets reconstructed by the identified number of primary singular values under four cases, which correspond to setting a threshold for feature curves directly (FC_threshold), both the correlation curve and the feature curve as the input of the LSTM classifier (LSTM_CC + FC), the correlation curve alone as the input (LSTM_CC), and the feature curve alone as the input (LSTM_FC). It is found that the performance of the method used feature curves directly is very unstable and the results are greatly affected by the curve fluctuations. However, the reconstruction under the case of “LSTM_CC + FC” is always advantageous over the other cases, especially for the target with low filling fraction and low SNR. One intuitive explanation could be the fluctuations caused by low SNR makes it difficult to identify the demarcation point of the shrunk peak of the correlation curve due to low filling fraction. The training dataset for the LSTM network consists of four kinds of simulated targets with different filling fractions (samples in Fig. 3(a)), each with different levels of noise added.
Fig. 4. Evaluation of image reconstructions based on the identified number of primary singular values by four cases, where “FC_threshold” stands for the method based on setting a threshold for feature curves directly, “LSTM_CC + FC” stands for both the correlation curve and the feature curve as the input of the LSTM classifier, “LSTM_CC” stands for the correlation curve alone as the input, and “LSTM_FC” stands for the feature curve alone as the input. (a) Bias of similarity of the reconstructed images corresponding to targets with increased filling fractions under high and low SNRs (high SNR: -37.1 dB, low SNR: -40.2 dB), respectively. (b) Bias of similarity of the recovered images corresponding to targets with increased SNRs (SNR 1=-40.2 dB, SNR 2=-38.8 dB, SNR 3=-37.9 dB, SNR 4=-37.1 dB) under low and high filling fractions, respectively.
To demonstrate the feasibility and robustness of the developed method, reconstructions are performed on targets with different combinations of filling fraction and SNR by the developed method and two established methods, including the STD-based method and the maximum ratio method. As shown in Fig. 5, four targets including simple sample with low filling fraction under high SNR and low SNR, and complex sample with high filling fraction under high SNR and low SNR are simulated. Samples without any noise and the simulated confocal images corresponding to four targets using digital-pinhole filtering are also presented for reference in the first and second columns of Fig. 5, respectively. Obviously, even under high SNR, the noise intentionally added to the reference sample is not suppressed and hence blurs the quality of the confocal images. To suppress the noisy background, image reconstruction through the iterative time reversal processing by different methods is further performed. The first method for comparison is the STD-based method and the results are shown in the third and fourth columns of Fig. 5. The noise contribution in the image results in monotonously increased STD with the increase of the number of primary singular values used for reconstruction, and thus all the singular values are included in reconstruction under the criterion of maximum STD. Therefore, the images reconstructed by the STD-based method without preprocessing are almost identical to the corresponding confocal cases. To make the STD-based method useful, preprocessing based on denoising and contrast enhancement is applied before STD calculation. The pre-processing operation consists of applying a Gaussian smoothing filter for denoising and then applying a grey scale transformation for image enhancement. In this case, the noisy background in the reconstructed images is almost suppressed. For different samples or different noise levels, preprocessing operations of different intensities are required to obtain satisfactory results (compared to reference samples). However, without appropriate preprocessing, the optimized M cannot be correctly identified, resulting missed information as noticed in the reconstructed images of complex targets. It should be emphasized that the preprocessing of the image cannot always be satisfactory as the targets to be recovered are unknown. The second method for comparison is the maximum ratio method expressed by the criterion:${\sigma _M}/{\sigma _{M + 1}}\textrm{ = }\max ({\sigma _i}/{\sigma _{i + 1}})$, and the results are shown in the fifth column of Fig. 5. The maximum ratio method is only feasible for the first target of simple sample under high SNR and fails for all the other three targets. The reconstructed results achieved by the developed method are presented in the last column. Satisfactory reconstructions are realized, confirming its robustness for all targets with different combinations of filling fraction and SNR. The similarity of the reconstructed images by different methods (STD-based method with processed images, maximum ratio method, the developed method) and the corresponding identified number n are shown in Table 1. The sensitivity of the method can be evaluated by the similarity versus different number n. However, this sensitivity is sample dependence. For example, a simple sample under low SNR in Fig. 5, the STD-based method provides a similarity of 0.6986 at the number n of 34, while the developed method offers a similarity of 0.6953 at the number n of 37. In contrast, another complex sample under high SNR in Fig. 5, the STD-based method provides a similarity of 0.5577 at the number n of 50, while the developed method offers a similarity of 0.7826 at the number n of 204.
Fig. 5. Reconstructions from the identified number of primary singular values by different methods on targets with different combinations of filling fraction and SNR. The size of the aperture in simulated confocal imaging is equal to the size of an Airy disk. High SNR corresponds to SNR=-37.1 dB, and low SNR corresponds to SNR=-40.2 dB.
Table 1. The identified number and corresponding similarity of three methods.
In fact, the scattering medium of biological tissues induces aberrations, which significantly degrades the imaging quality. To demonstrate that the developed method does not suffer from aberrations, two typical aberrations including axisymmetric (spherical aberration) and non-axisymmetric (coma) aberrations have been introduced in simulations. The results demonstrate that the developed method can successfully compensate the aberrations and achieve satisfactory reconstructions (Fig. 6). In addition, the virtual digital pinhole used in matrix pre-processing affects the quality of image reconstructions as well. A virtual digital pinhole of appropriate size is required to filter most of the multiple-scattered contributions, while retaining as many of the single-scattered contributions as possible. The simulations demonstrate that, for high SNR cases, the pinhole size can be enlarged appropriately to retain more single-scattered signals (the first row of Fig. 7). On the contrary, when SNR is low, the pinhole size needs to be reduced to filter out more multiple-scattered signals (the second row of Fig. 7), otherwise the SVD method cannot separate the residual multiple and single-scattered signals effectively. However, under extremely low SNR (SNR=-44.2 dB, corresponding imaging-thickness=$12.38{l_s}$), it is too scattering to reconstruct the image (the third row of Fig. 7).
Fig. 7. Effect of the size of the virtual digital pinhole on reconstructions under different SNRs. High SNR= -35.1 dB; low SNR= -40.2 dB; extremely low SNR=-44.2 dB.
4. Conclusion
The noise due to multiple scattering available in tissue imaging cannot be ultimately suppressed by digital-pinhole filtering based confocal imaging, and further iterative time reversal processing is definitely required to remove the noisy background. However, the quality of the reconstructed image through iterative time reversal processing heavily depends on the accurate determination of the optimized number of singular values, which is related to the filling fraction and SNR rather than the distribution of the target. To this regard, we have developed a method based on correlation between image pair reconstructed from different numbers of singular values and corresponding remainders. The obtained correlation curve and another feature curve fetched from the former are then fed to the LSTM network to identify the optimized number of primary singular values for best reconstruction. It is confirmed that taking both the correlation curve and the feature curve as the input is of great superiority in the two-class classification in contrast to the correlation curve or the feature curve alone.
Due to the fact that noise contribution in the image results in monotonously increased STD, the images reconstructed by the STD-based method without preprocessing are almost identical to the corresponding confocal cases. If appropriate preprocessing is performed before STD calculation, the noisy background in the reconstructed images can be suppressed. However, preprocessing of the image cannot always be satisfactory as the targets to be recovered are unknown. The maximum ratio method is only feasible for the simple sample under high SNR and fails in complex sample or under low SNR. In comparison, satisfactory reconstructions are realized for all simulated targets by the method developed in this study, confirming its robustness to samples with different combinations of filling fraction and SNR.
Future study will contribute to verifying the effectiveness of this method in tissue imaging experiments and making further refinement of the method.
Funding
National Natural Science Foundation of China (11974310, 31927801, 61905214, 62035011); National Key Research and Development Program of China (2017YFA0700501, 2019YFE0113700); Natural Science Foundation of Zhejiang Province (LR20F050001); Fundamental Research Funds for the Central Universities (2020XZZX005-07) .
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request
References
1. A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nat. Photonics 6(5), 283–292 (2012). [CrossRef]
2. S. Kang, S. Jeong, W. Choi, H. Ko, T. D. Yang, J. H. Joo, J.-S. Lee, Y.-S. Lim, Q. H. Park, and W. Choi, “Imaging deep within a scattering medium using collective accumulation of single-scattered waves,” Nat. Photonics 9(4), 253–258 (2015). [CrossRef]
3. A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, “Smart optical coherence tomography for ultra-deep imaging through highly scattering media,” Sci. Adv. 2(11), e1600370 (2016). [CrossRef]
4. Q. Yang, J. Cao, Y. Miao, J. Zhu, and Z. Chen, “Extended imaging depth of en-face optical coherence tomography based on fast measurement of a reflection matrix by wide-field heterodyne detection,” Opt. Lett. 45(4), 828–831 (2020). [CrossRef]
5. Q. Yang, Y. Miao, T. Huo, Y. Li, E. Heidari, J. Zhu, and Z. Chen, “Deep imaging in highly scattering media by combining reflection matrix measurement with Bessel-like beam based optical coherence tomography,” Appl. Phys. Lett. 113(1), 011106 (2018). [CrossRef]
6. L. Borcea, G. Papanicolaou, C. Tsogka, and J. Berryman, “Imaging and time reversal in random media,” Inverse Problems 18(5), 1247–1279 (2002). [CrossRef]
7. A. Aubry and A. Derode, “Detection and imaging in a random medium: A matrix method to overcome multiple scattering and aberration,” J. Appl. Phys. 106(4), 044903 (2009). [CrossRef]
8. J. Veraart, E. Fieremans, and D. S. Novikov, “Diffusion MRI noise mapping using random matrix theory,” Magn. Reson. Med. 76(5), 1582–1593 (2016). [CrossRef]
9. S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett. 104(10), 100601 (2010). [CrossRef]
10. A. Aubry and A. Derode, “Singular value distribution of the propagation matrix in random scattering media,” Waves in Random and Complex Media 20(3), 333–363 (2010). [CrossRef]
11. M. B. Roumeliotis, R. Z. Stodilka, M. A. Anastasio, E. Ng, and J. J. L. Carson, “Singular value decomposition analysis of a photoacoustic imaging system and 3D imaging at 0.7 FPS,” Opt. Express 19(14), 13405–13417 (2011). [CrossRef]
12. K. Konstantinides and K. Yao, “Statistical analysis of effective singular values in matrix rank determination,” IEEE Transactions on Acoustics, Speech, and Signal Processing 36(5), 757–763 (1988). [CrossRef]
13. B. Pan, “Reliability-guided digital image correlation for image deformation measurement,” Appl. Opt. 48(8), 1535–1542 (2009). [CrossRef]
14. I. Sutskever, O. Vinyals, and Q. V. Le, “Sequence to sequence learning with neural networks,” in Advances in neural information processing systems, 2014, 3104–3112.
15. S. Shahjahan, A. Aubry, F. Rupin, B. Chassignole, and A. Derode, “A random matrix approach to detect defects in a strongly scattering polycrystal: How the memory effect can help overcome multiple scattering,” Appl. Phys. Lett. 104(23), 234105 (2014). [CrossRef]
16. B. Pan, H. Xie, and Z. Wang, “Equivalence of digital image correlation criteria for pattern matching,” Appl. Opt. 49(28), 5501–5509 (2010). [CrossRef]
17. H. Hirschmuller and D. Scharstein, “Evaluation of stereo matching costs on images with radiometric differences,” IEEE Trans. Pattern Anal. Mach. Intell. 31(9), 1582–1599 (2009). [CrossRef]
