Three modifications are shown to improve resolution and reduce noise amplification in endoscopic imaging through multi-mode fiber using optimization-based reconstruction (OBR). First, random sampling patterns are replaced by sampling patterns designed to have more nearly equal singular values. Second, the OBR algorithm uses a point-spread function based on the estimated spatial frequency spectrum of the object. Third, the OBR algorithm gives less weight to modes having smaller singular values. In simulations for a step-index fiber supporting 522 spatial modes, the modifications yield a 20% reduction in image error (l2 norm) in the noiseless case, and a further 33% reduction in image error at a 22-dB shot noise-limited SNR as compared to the original method using random sampling patterns and OBR.
© 2014 Optical Society of America
The use of multi-mode fibers (MMFs) for imaging has long been of interest [1,2], and the past few years have seen particularly productive developments in the field of MMF endoscopic imaging [3–12]. An endoscope employing one MMF provides resolution superior to, and is much more compact than, current commercial fiber endoscopes, which may use either a bundle of single-mode fibers or one single-mode fiber with a scanning probe head . Light signals propagating through a MMF are subject to random scattering by index imperfections, bends and other perturbations. All endoscopic MMF imaging methods, prior to image recording, employ some calibration method that requires access to both ends of the fiber. Since bending changes a fiber’s propagation characteristics, all MMF endoscopes to date have been rigid, although compensation of bending has been investigated [8,9].
One endoscopic MMF imaging method  uses random intensity patterns to sample an object and employs an optimization-based image reconstruction (OBR) algorithm. This method at once enhances the resolution to four times the number of spatial modes propagating in the MMF and requires the least complicated calibration setup among all imaging methods. The OBR algorithm decomposes the set of random intensity sampling patterns into a set of intensity modes using a singular-value decomposition (SVD), and reconstructs the image as a linear combination of the intensity modes based on the powers reflected from the object when sampled by the random intensity patterns. This technique is similar to compressed sensing [14,15].
Despite the advantages of this imaging method, reconstructed images were found to be noisy , suggesting that noise amplification occurs in OBR. In this paper, we quantify the noise amplification occurring in OBR, and describe three methods to reduce it. First, we sample the object using a set of intensity patterns designed to have singular values more equalized than a set of random intensity patterns. Second, we modify the OBR algorithm to use a point-spread function (PSF) based on the estimated spatial frequency spectrum of the object. Third, we diminish the influence of intensity modes with smaller singular values on the reconstructed image. Using simulation, we compare images obtained using noise-reduced OBR to those obtained by localized spot sampling with localized reconstruction in the presence of noise, and subject to constraints on the number of sampling patterns and the power illuminating the object through the MMF.
The remainder of this paper is as follows. In Section 2, we describe the imaging system and a model for the electric fields illuminating the object and captured by the MMF, and briefly review the OBR algorithm described in . In Section 3, we analyze the image noises obtained using OBR or localized sampling/reconstruction and quantify the noise amplification of OBR. In Section 4, we introduce three methods to improve the noise performance of OBR. In Section 5, we present simulated images obtained using OBR or localized sampling/reconstruction under the constraints mentioned above. We discuss the results in Section 6 and provide conclusions in Section 7.
2. Model for object sampling and image reconstruction
2.1 Imaging system and model for electric fields
The imaging system simulated here is identical to the experimental system in , except it uses a step-index MMF instead of the graded-index MMF used in . As shown in Fig. 1, a 1550-nm laser beam illuminates a phase-only spatial-light modulator (SLM). The beam reflected from the SLM is coupled into a step-index MMF having 50-μm core diameter. The MMF is 1 m long, and its output facet is anti-reflection-coated with reflectivity <1%. Calibration is performed by imaging intensity patterns appearing at the MMF output onto a camera with 256 × 256 pixels, as shown in Fig. 1(a). After calibration, an object to be imaged is placed at the MMF output, as shown in Fig. 1(b). A sequence of M different intensity patterns is displayed at the MMF output and is used to sample the object. When each intensity pattern is displayed, the total power reflected from the object and coupled back into the MMF is measured using a power meter. Then, using the sequence of M reflected power values, an image is reconstructed by OBR.
A step-index MMF with core diameter d and numerical aperture NA supports electric field modes per polarization at wavelength for large V . For small NA, we can use the weakly guiding approximation and use analytically computed linearly polarized modes (LP) as a basis. For MMFs with higher NAs, we instead use the numerically computed vector modes of the MMF, which deviate slightly from the analytically computed LP modes. The index ranges from 0 to , while the index m ranges from 1 to , where both and depend on V . Each is two-fold degenerate for , for the and polarizations; and fourfold degenerate for , twice for the and polarizations and twice for a rotation of the mode about the axis by radians (where is the direction of propagation). We assume that both the MMF and sample to be imaged have negligible birefringence; thus a single SLM provides full control of the intensity at the output of the MMF, and we can ignore polarization effects during imaging. For , the complex electric field inside the MMF in one polarization is described in cylindrical coordinates by:
Let P1 denote the MMF output plane and P2 denote the object plane. At P1, the total forward-propagating electric field is a linear combination of the LP modes, and is described by:Eq. (1), and the index has been added to differentiate between the degeneracies of .
Free-space propagation over a distance is represented by , and is modeled here by Rayleigh-Sommerfeld diffraction . At P2 the total forward-propagating electric field is thus:18], meaning that the reflected electric field is the product of the incident electric field and the complex reflectance of the object in each polarization on plane P2:
2.2 Optimization-based reconstruction algorithm
The relation between p and in Eqs. (4) and (5) is not linear. Solving for from p using Eq. (5) requires coherent detection of and is in general computationally difficult for a MMF with the number of modes required for imaging. However, if we instead assume that the fraction of the reflected light recaptured by the MMF is independent of the object, and that the object has equal reflectivity in both polarizations, Eq. (5) can be replaced by:Eqs. (4) and (5) and assuming is a perfect mirror, or measured experimentally. In Eq. (8), is the intensity pattern at P2 measured by the camera during calibration and then multiplied by . We note that in actuality is not independent of , so the approximation used in Eq. (6) effectively introduces an error in p. The impact of this error will be seen in subsequent images, but the effect is small if the used to sample the object are statistically similar (e.g. random patterns).
If we employ a sequence of M intensity patterns during imaging, we obtain a corresponding vector of powers p. If we discretize and to a grid of L points and reshape them into row vectors, we obtain the simple expression:3]:
3. Noise of image reconstruction methods
In the previous section we assumed that the received power specified in Eq. (9) is noiseless. We now introduce noise. Given an apparatus, as the optical power is increased, the limiting noise source progresses from thermal to shot to intensity noise . In the following section, we analyze all three regimes.
3.1 Noise of optimization-based image reconstruction
We assume that the calibration images are noiseless. This assumption is valid if sufficient time is available for calibration using the camera to make noise negligible and there is no movement of the MMF during imaging. The noisy reconstructed image is then:Eq. (21), is the component of the noise in the basis of intensity modes, which is useful in the following derivation. We model as a vector of i.i.d. uniformly distributed random variables between 0 and 1. Solving for , we have:Eqs. (22) and (23), we have used the fact that unitary matrices preserve the l2 norm. The approximation in Eq. (22) is accurate as long as all patterns in have total powers within the same order of magnitude. Similarly, we find:Eq. (23). Using Eqs. (22)–(25), we finally have:
3.2 Noise of localized reconstruction
In this imaging method, the sampling patterns at the MMF output plane P2 are localized spots on a regular grid of points. As each sampling pattern illuminates the object, the power reflected from the object yields the intensity in the corresponding region of the reconstructed image, a process termed localized reconstruction . For localized reconstruction using M spots, the corresponding relation between the noise n, spot sampling patterns and image w is:3], and s is a vector of the reciprocal of the power that couples back into the MMF when imaging a perfect mirror using . The noises n are as given in Eqs. (15)–(17). The corresponding RMSDs of image error due to noise are:Eqs. (27)–(29), to the image error RMSDs using localized reconstruction, given by Eqs. (32)–(34).
4. Techniques to improve noise performance
4.1 Sampling using optimal patterns
The dependence of the image error on the singular values of in Eqs. (27)–(29) suggests that noise amplification can be reduced by choosing a specific set of sampling patterns with more equalized singular values. We heuristically expect a grid of spots (identical to that used in localized sampling and reconstruction) to be a good set of sampling patterns: among all possible sampling patterns satisfying the physical constraint of positive intensity (), the overlap integral between different sampling patterns is minimized for a grid of spots, which is a good characteristic for minimizing. We briefly note that forming specific sampling patterns is necessarily more difficult than forming random patterns, since light undergoes random mode coupling during MMF propagation. Conveniently, methods have been described to form the entire set of spot patterns with total time linearly proportional to the number of modes N [4–6].
The impacts of shot and thermal noise are inversely dependent on the integrated total power in the sampling patterns, so it is desirable to maximize the total power in these patterns. However, this power must be limited to prevent damage to tissue [20–22], which occurs at powers well below those that cause any damage to the MMF . Tissue damage is typically caused by thermal effects, which depend on the irradiance incident on the tissue , so we constrain the maximum irradiance of each sampling pattern such that (here we use the terms irradiance and intensity interchangeably in referring to the power per unit area). This constraint is particularly important when MMFs with high NA are used, as they permit the total power to be concentrated in a high-irradiance spot that is much smaller than the MMF core. However, since thermal effects are not directly localized to the irradiated area, a large and uniform sampling pattern can cause more damage than a smaller one with the same irradiance. To account for this, we additionally constrain the total power of each sampling pattern such that.
We plot the noise amplification for OBR when imaging using random and spot sampling patterns in Fig. 2, for both constraints. The intensity patterns have been generated using the model of Section 2, and the noise amplification is calculated using Eqs. (27)–(29) and Eqs. (32)–(34) of Section 3.
Under the total power constraint shown in Fig. 2(b), when imaging with random patterns the noise amplification increases in proportion to the square root of the number of LP modes, while for imaging with spot patterns it remains constant. This applies to all types of noise. Thus imaging using spot patterns offers a clear improvement over imaging with random patterns under the total power constraint, with the drawback of a more difficult calibration of the imaging apparatus. Conversely when operating under the maximum irradiance constraint of Fig. 2(a), spot patterns have less noise amplification than random patterns only for intensity noise. Thus imaging using random patterns is superior unless intensity noise dominates.
The results of Fig. 2 are for a specific MMF core diameter of 50 μm, but we expect the results to be similar for a constant number of modes independent of core diameter (with the distance to the object scaled according to the Rayleigh range of the MMF ). We note also that MMFs are commonly available with well over 500 modes, but it may not be possible to use all these modes for imaging due to time constraints; for example, for 500 modes and a grid of spots with an oversampling factor of 2 × in the x- and y-dimensions, 8000 spots need to be sequentially formed during imaging to exploit all the modes.
4.2 Exploiting known properties of the object
In , the image reconstruction algorithm was chosen under the assumption that nothing is known about r, the object to be imaged. However, at a minimum we do know that r must be positive and have a limited spatial frequency bandwidth. If instead r were i.i.d., it would have a uniform spatial frequency spectrum up to the spatial sampling rate of the camera used to record the intensity patterns. For realistic objects, the spatial frequency spectrum drops off well before this. With this motivation, we modify the OBR algorithm to exploit any information known about the spatial frequencies of the object. This requires a more general objective function than the one described in , which was:
Instead of Eq. (35), we minimize the l2 norm between the object intensity reflectance and the image:Eq. (35) is equivalent to Eq. (36) with the constraint that when r is i.i.d. with mean 0. For Eq. (35):Eq. (36):Equation (38) has the same result as Eq. (37) above.
We note that the continuous form of Eq. (36) allows a simplification if we assume that the point-spread function of the image is spatially invariant. In that case:
In minimizing Eq. (36), we restrict ourselves to solutions of the form:Eq. (10), and remains to be determined. We chose Eq. (42) since it can be simplified to , which allows us to interpret as the point-spread function of the image at pixel i, where is the ith row of . We know the objective function for the optimal PSF from Eq. (41), and so we use the following equation to solve for every :Equations (42) and (43) are all the equations required to reconstruct the image, and replace Eq. (5).
Using the electric field model described in Section 2, we simulate imaging of an object in MATLAB. We compare the images obtained with the original and new algorithms in Fig. 3.
We see that the image quality is substantially improved when using the new algorithm, which reduces the image error RMSD from 9.3% to 7.5% of the maximum image brightness – or equivalently, an increase in the peak signal-to-noise ratio from 20.6 to 22.5 dB. We do note that to use Eq. (43), it is necessary to know r. This seems to present a paradox: if we know r, we already have a perfect image of the object and imaging is unnecessary. In practice, however, the improvement shown in Fig. 3 can be obtained even if we have only an estimate of the spatial frequency spectrum of r. The r used in Eq. (43) is shown in Fig. 3(d), and differs substantially from the actual r shown in Fig. 3(c).
We also briefly mention the issue of computation time. Equation (43) is a quadratic optimization problem with N variables, which can be solved efficiently by convex optimization methods. In practice, the amount of time to solve Eq. (43) on a personal computer is under 1 s. However, the number of separate that must be solved is equal to L, the number of pixels in the final image (e.g., L = 2562 for a 256 × 256 pixel image). Hence, real-time computation of Eq. (43) for requires solving for the different in parallel. We finally note that is not unitary, and thus the noise amplification no longer depends on the singular values of but instead on the singular values of the matrix , which are not the diagonal elements of D. In practice, the changes to the noise amplification are minor.
4.3 Noise-reduced image reconstruction
Noise amplification can be further mitigated during image reconstruction by adding a constraint to the elements of each :Eq. (43), this becomes:Eq. (43), the singular values are more equalized and thus the contribution to noise from low singular values is limited. However, by limiting the contribution to the image by intensity modes with low singular values, the resulting image has reduced resolution. The choice of determines the tradeoff between noise amplification and resolution; when α is sufficiently high, Eq. (45) reduces to Eq. (43). We denote this value of as , given by , where is the solution to Eq. (43).
We observe that the image reconstructed using Eqs. (44) and (45) is less noisy, with a reduction in image error RMSD from 14.0% to 9.4%. For every noise SNR there is an optimal value of that balances the image errors due to reduced resolution and noise, and thus minimizes the total image error. When r is unknown, the optimal value of is also unknown, but can be chosen such that the image error RMSD is subjectively minimized, which is how has been chosen here. Using Eqs. (44) and (45) reduces thermal, shot and intensity noise equally, and also reduces the error introduced due to the approximation of Eq. (6). Since this is always present even in the noiseless case, conservative (large) values of lead to reductions in image error RMSD even when there is no noise.
5. Comparison of image quality
To directly compare the image error between images reconstructed using OBR and localized reconstruction, we simulate imaging an object using both methods in the presence of noise and plot the image error versus the power noise in Fig. 5. The apparatus and noise model are as described in Sections 2-3. The imaged object is shown previously in Fig. 3(d). In the plots of Fig. 5, we have normalized the noise values: at any SNR, , where denotes the average measured power for imaging using OBR. Equal values of , , and equal numbers of sampling patterns were used for both OBR and localized reconstruction. Results for both the maximum irradiance and total power constraints are shown in Fig. 5. For the maximum irradiance constraint the sampling patterns for OBR are random, while for the total power constraint the sampling patterns for OBR are spots, for reasons described in Section 4.1.
The noise amplification curves shown in Fig. 5 behave as predicted based on Fig. 2. For low levels of noise, images reconstructed using OBR have a lower image error RMSD than those using localized reconstruction, because the former have higher resolution. For high levels of noise, images reconstructed using OBR have higher image error RMSD than localized reconstruction, due to noise amplification. Under the maximum irradiance constraint, the image error using OBR is less than that of localized reconstruction for all thermal noise levels, and for noises up to −18 dB and −19 dB for shot and intensity noises, respectively. When imaging using spot patterns under the total power constraint, the image error using OBR is less than that of localized reconstruction for noises up to −15 dB, −15 dB and −16 dB for thermal, shot and intensity noises, respectively.
In Fig. 6, we show images corresponding to three regimes on the plot of Fig. 5(b) as an example. These images are for the total power constraint and various levels of shot noise. OBR achieves superior resolution compared to localized reconstruction for low noise, as seen in Figs. 6(a) and 6(d), but noise amplification occurs for higher amounts of noise as seen in Figs. 6(c) and 6(f). In Fig. 6(a), for 25 dB SNR the image is still slightly noisy; this is due to the error introduced by the assumption made in Eq. (6).
The advantages of OBR using random sampling patterns are faster calibration time and higher resolution compared to localized reconstruction, which come at the expense of noise amplification. By using the improved methods described in Section 4, an image can be reconstructed even with a moderate level of noise. If the noise level is high, averaging of the measured power values and reconstructed images may need to be done (as in ).
In , it was assumed that the apparatus can receive all the power reflected from the object. In Fig. 6, we demonstrate that this assumption introduces an error to the image. Future work may investigate exciting only the lower-order modes of high-NA MMFs for imaging, which is possible if the coupling between higher and lower order modes is limited . By doing so, we would expect an increased fraction of the reflected power to couple back into the MMF due to the availability of higher-order modes, minimizing this error.
In this work, we have imaged objects that were placed 5 μm from the end of the MMF. As the separation between the MMF and the object increases, less power couples back into the MMF, and resolution decreases for both localized reconstruction and OBR, even in the absence of noise. For a MMF with larger core radius (but smaller NA such that the total number of modes is preserved), this distance increases.
We have quantified the noise amplification for the optimization-based multi-mode fiber endoscopic imaging method originally described in , and described three improvements to reduce noise amplification. A comparison of images obtained using noise-reduced OBR and the original OBR yielded a 20% reduction in image error for the noiseless case, and a further 33% reduction in image error at 22 dB SNR for the shot noise-limited case. As a comparison, we simulated imaging an object using both OBR and localized reconstruction under two constraints: on the maximum irradiance or on the total power of the sampling patterns. When imaging using random sampling patterns under the maximum irradiance constraint, we found that the image error using OBR was less than that of localized reconstruction for all SNRs when thermal noise-limited, and down to 18 dB and 19 dB SNR for shot and intensity noise-limited cases, respectively. When using spot sampling patterns under the total power constraint, we found that the image error using OBR was less than that of localized reconstruction for SNRs down to 15 dB, 15 dB and 16 dB for thermal, shot and intensity noise-limited cases, respectively.
Support from the Samsung GRO Program is gratefully acknowledged.
References and links
1. A. Yariv, “Three-dimensional pictorial transmission in optical fibers,” Appl. Phys. Lett. 28(2), 88 (1976). [CrossRef]
2. B. Fischer and S. Sternklar, “Image transmission and interferometry with multimode fibers using self-pumped phase conjugation,” Appl. Phys. Lett. 46(2), 113 (1985). [CrossRef]
4. I. N. Papadopoulos, S. Farahi, and C. Moser, “High resolution, lenseless endoscope based on digital scanning through a multimode optical fiber,” Opt. Express 4(2), 17598–17603 (2013).
7. Y. Choi, C. Yoon, M. Kim, T. D. Yang, C. Fang-Yen, R. R. Dasari, K. J. Lee, and W. Choi, “Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber,” Phys. Rev. Lett. 109(20), 203901 (2012). [CrossRef] [PubMed]
8. S. Farahi, D. Ziegler, I. N. Papadopoulos, D. Psaltis, and C. Moser, “Dynamic bending compensation while focusing through a multimode fiber,” Opt. Express 21(19), 22504–22514 (2013). [CrossRef] [PubMed]
10. I. N. Papadopoulos, S. Farahi, C. Moser, and D. Psaltis, “Increasing the imaging capabilities of multimode fibers by exploiting the properties of highly scattering media,” Opt. Lett. 38(15), 2776–2778 (2013). [CrossRef] [PubMed]
11. S. Bianchi, V. P. Rajamanickam, L. Ferrara, E. Di Fabrizio, C. Liberale, and R. Di Leonardo, “Focusing and imaging with increased numerical apertures through multimode fibers with micro-fabricated optics,” Opt. Lett. 38(23), 4935–4938 (2013). [CrossRef] [PubMed]
14. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). [CrossRef]
15. E. J. Candes and T. Tao, “Near-optimal signal recovery from random projections: universal encoding strategies?” IEEE Trans. Inf. Theory 52(12), 5406–5425 (2006). [CrossRef]
16. J. A. Buck, Fundamentals of Optical Fibers, 2nd ed. (John Wiley, 2004).
17. J. Mertz, Introduction to Optical Microscopy (Roberts and Company, 2010).
18. J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts and Company, 2004).
19. G. P. Agrawal, Fiber-Optic Communication Systems, 4th ed. (John Wiley, 2010), Chap. 3.
20. A. J. Welch, M. J. C. Van Gemert, W. M. Star, and T. Optics, Optical-Thermal Response of Laser-Irradiated Tissue, 2nd ed. (Springer, 2011), Chap. 3.
21. A. J. Welch, “The thermal response of laser irradiated tissue,” IEEE J. Quantum Electron. 20(12), 1471–1481 (1984). [CrossRef]
24. M. B. Shemirani, W. Mao, R. A. Panicker, and J. M. Kahn, “Principal modes in graded-index multimode fiber in presence of spatial- and polarization-mode coupling,” J. Lightwave Technol. 27(10), 1248–1261 (2009). [CrossRef]