-
光在多模光纤中的传播由光纤的传输矩阵描述。光纤的远端作为输入端,光纤的近端作为输出端,如图1所示。只考虑一个偏振方向,经过空间离散化后,输入端和输出端光场可以用矢量
$ {\boldsymbol{E}}^{\mathrm{i}\mathrm{n}},{\boldsymbol{E}}^{\mathrm{o}\mathrm{u}\mathrm{t}} $ 表示。设光纤的传输矩阵为$ \boldsymbol{A} $ ,则输入光场和输出光场的关系如下:$$ \begin{array}{c}{\boldsymbol{E}}^{\mathrm{o}\mathrm{u}\mathrm{t}}={\boldsymbol{A}{{\boldsymbol{E}}}}^{\mathrm{i}\mathrm{n}}\end{array} $$ (1) 位于光纤输出端的相机只能探测到光强,设探测到的光强为
$ {I}_{m} $ ,${{\boldsymbol{A}}}_{mn}$ 为传输矩阵的元素,下标表示相机上不同位置的像素,光强与输入端光场的关系为:$$ \begin{array}{c}{I}_{m}={\left|{\displaystyle\sum }_{n}{{\boldsymbol{A}}}_{mn}{{\boldsymbol{E}}}_{n}^{\mathrm{i}\mathrm{n}}\right|}^{2}\end{array} $$ (2) 在已知输出端光强和传输矩阵的情况下,可以使用相位恢复算法求解以上方程,计算出输入端光场。文中中使用的算法是相位恢复变分贝叶斯期望最大化算法[15](Phase Recovery Variational Bayes Expectation-Maximization, prVBEM),并使用谱初始化方法进行初始化[12-13]。
相机上每一个像素都对应一个
$ {I}_{m} $ ,也对应方程组中的一个方程。方程组的方程达到一定数量时,就足以正确地解出输入光场${{\boldsymbol{E}}}_{n}^{\mathrm{i}\mathrm{n}}$ ,方程数量过多则可能消耗不必要的计算资源。因此可以对光强进行适当的下采样,只使用其中一部分像素的光强。文中只考虑均匀采样的情况,测试采样方案主要由采样间隔决定。采样点的位置如下式:$$ \begin{array}{c}{\boldsymbol{r}}={n}_{1}{\boldsymbol{e}}_{1}+{n}_{2}{\boldsymbol{e}}_{2}\end{array} $$ (3) 式中:
$ {n}_{1},{n}_{2} $ 为整数;$ {\boldsymbol{e}}_{1},{\boldsymbol{e}}_{2} $ 为固定的两个向量。图2中展示了选择不同的$ {\boldsymbol{e}}_{1},{\boldsymbol{e}}_{2} $ 时采样点的分布情况。图 2 采样方案示意图,图中箭头表示向量
$ {\boldsymbol{e}}_{1},{\boldsymbol{e}}_{2} $ 的方向。(a)采样间隔为$ \sqrt{{1}^{2}+{3}^{2}}=3.16 $ pixel,$ {\boldsymbol{e}}_{1} $ =(3,1),$ {\boldsymbol{e}}_{2} $ = (−1,3);(b)采样间隔为3 pixel,$ {\boldsymbol{e}}_{1} $ =(3,0),$ {\boldsymbol{e}}_{2} $ =(0,3)Figure 2. Schematics of different sampling manners, the arrows in the figure represent the directions of the vectors
$ {\boldsymbol{e}}_{1},{\boldsymbol{e}}_{2} $ . (a) Sampling interval is$ \sqrt{{1}^{2}+{3}^{2}}=3.16 $ pixel,$ {\boldsymbol{e}}_{1} $ =(3,1),$ {\boldsymbol{e}}_{2} $ = (−1,3); (b) Sampling interval is 3 pixel,$ {\boldsymbol{e}}_{1} $ =(3,0),$ {\boldsymbol{e}}_{2} $ =(0,3)实验装置如图3所示。光束从532 nm激光器(长春新产业,MSL-FN-532-100 mW)出射经过扩束镜L1(大恒光电,GCO-140113)扩大后被用于照明。照明光束被DMD (Texas Instruments,DLP LightCrafter Display 2000, 640 × 360 pixel,像素尺寸7.56 μm)反射,然后通过200 mm焦距凸透镜TL1和显微物镜OL1(Olympus Plan N 10×)耦合进入多模光纤MMF(大恒光电,DH-FMM105-FC-1 A,长度1 m,纤芯直径105 μm,0.1 NA)中。TL1和OL1将DMD上的图案成像在光纤的端面上。光在光纤中传播一段距离后从另一端射出,显微物镜OL2 (Olympus Plan N 20×)和镜头TL2(Nikon NIKKOR 80~200 mm,对焦距离调为无穷远,焦距调为200 mm)组成的成像系统将光纤端面放大后成像在相机感光面(pco.edge 4.2 Camera Link,像素阵列2 048×2 048,像素尺寸6.5 µm×6.5 µm)上。在TL2与OL2之间插入了偏振片P1,使相机只接收一个偏振方向的光场。根据参考文献[25]中描述的方法对DLP LightCrafter Display 2000模块进行了简单改装,拆除了原有的照明模块,只保留其中的DMD芯片。将3×3个DMD像素组成1个大像素,选择DMD上大小为28×28=784个大像素的区域用于显示输入图像。重建图像使用的CPU为Intel Core i7-6700。算法用matlab实现,prVBEM算法采用了参考文献[15]的作者所公布的代码。
-
成像前需要事先标定光纤传输矩阵,标定传输矩阵也使用相位恢复算法。标定时,首先测得一系列不同输入图像
${{\boldsymbol{E}}}_{nk}^{{\rm{in}}}$ 对应的输出光强${{{I}}}_{mk}^{{\rm{out}}}$ ,设光纤的传输矩阵的元素是${{\boldsymbol{A}}}_{mn}$ ,下标n,m表示不同像素,下标k表示的是不同的输入图像,则有:$$ \begin{array}{c}{{{I}}}_{mk}^{{\rm{out}}}={\left|{\displaystyle\sum }_{n}{{{\boldsymbol{A}}}_{mn}{\boldsymbol{E}}}_{nk}^{{\rm{in}}}\right|}^{2}\end{array} $$ (4) 已知
$ {{\boldsymbol{E}}}_{nk}^{{\rm{in}}} $ 和${{{I}}}_{mk}^{{\rm{out}}}$ ,此时将传输矩阵${{\boldsymbol{A}}}_{mn}$ 的某行(例如${{\boldsymbol{A}}}_{1 n}$ )视为未知数,则可以用相位恢复算法从方程组得到传输矩阵的一行,逐行运行相位恢复算法可以得到完整的传输矩阵。完成传输矩阵的标定后进行图像重建的实验, 不同采样间隔的重建效果如图4所示。不同的采样间隔对应着不同的${\boldsymbol{e}}_{1}\mathrm{、}{\boldsymbol{e}}_{2}$ ,间隔2.24 pixel时的$ {\boldsymbol{e}}_{1}\mathrm{、}{\boldsymbol{e}}_{2} $ 为$ \left(\mathrm{2,1}\right)\mathrm{和}(-\mathrm{1,2}) $ , 3.0、3.61、4.12 pixel时的$ {\boldsymbol{e}}_{1}\mathrm{、}{\boldsymbol{e}}_{2} $ 为$ \left(\mathrm{3,0}\right)\mathrm{和}\left(\mathrm{0,3}\right) $ ,$ \left(\mathrm{3,2}\right)\mathrm{和}(-\mathrm{2,3}) $ ,$ \left(\mathrm{4,1}\right)\mathrm{和}(-\mathrm{1,4}) $ 。可以看到,在采样点数不变的情况下,随着采样间隔增大,恢复结果逐渐变得更接近真实图像。图4中使用了余弦相似度定量描述图像重建的效果,假设真实图像为${{{I}}}_{m}$ ,重建图像为${{{I}}}_{m}^{\text{'}}$ ,下标m表示不同像素,则余弦相似度的定义如下:图 4 不同采样间隔下的结果。 (a)、(b)采样点数 2 400 时的重建图像和采样点分布; (c)、(d)采样点数 2 900 时的重建图像和采样点分布;(a)和(c)下方的数值为重建结果与真实图像的余弦相似度; (b)和(d)中灰色区域为光强包络区域,白点为采样点。右侧插入显示了放大的红色方块内的点分布
Figure 4. Results at different sampling intervals. (a), (b) Reconstruction results and the position distribution of sampling points when the number of sampling points is 2 400; (c), (d) Reconstruction results and the position distribution of sampling points when the number of sampling points is 2 900; Value below (a) and (c) is the cosine similarity between the reconstruction results to the ground truth; Grey area in (b) and (d) is the intensity envelope region, and the white dots are the sampling points. Zoomed-in view of the distribution of points in the red box is shown on the right
$$ \begin{array}{c}similarity=\dfrac{\displaystyle\sum _{m}{{{I}}}_{m}{{{I}}}_{m}^{\text{'}}}{\sqrt{{\displaystyle\sum }_{m}{{{I}}}_{m}^{2}}\sqrt{\displaystyle\sum _{m}{{{I}}}_{m}^{\text{'}2}}}\end{array} $$ (5) 重建过程使用的采样点数越少,则需要的时间和空间通常越少,标定过程也更简单。为了确定不同采样方案需要的采样点数,需要得出图像重建效果与各个条件的关系。从MNIST数据集中选择10张图像,在不同条件下重建图像,并计算10张图像的重建结果与真实图像的余弦相似度的平均值。结果如图5(a)所示,可以看到相似度随着采样点数提高而增加。由于光纤截面的像素数量有限,总采样点数不能无限增加,所以采样间隔较大时采样点数只能在一个小范围内变化。将重建图像与真实图像的平均相似度达到0.95定义为恢复成功,使用光滑样条对图中曲线进行插值,可求出不同采样间隔下成功恢复所需要的采样点数。从图5(b)中可以看到,随着采样间隔的增大,需要的采样点数首先快速下降,之后趋于稳定。在采样间隔为5.39时需要的采样点数取得最小值2 709,约为输入图像像素数量的3.5倍。采样点数量受到光纤截面面积的限制,采样间隔越大则采样点数量上限越低。如图4、图5(c)所示,还可以看到,虽然采样间隔大时需要的采样点数较少,但是采样区域的面积依然逐渐增大。
图 5 (a)不同采样间隔下相似度与采样点数的关系;(b)相同图像恢复质量下,最少采样点数随采样间隔的变化,虚线表示采样点数量上限;(c) 最小采样面积随采样间隔的变化,虚线表示光纤面积
Figure 5. (a) Similarity versus number of sampling points at different sampling intervals; (b) Variation of minimum number of sampling versus sampling interval under the same image quality, the dashed line represents the upper bound of the number of sampling points; (c) Variation of minimum sampling area versus sampling interval, the dashed line represents the area of the fiber
从理论上对所需采样点数的变化进行解释,相位恢复算法通常要求传输矩阵的不同行是独立同分布的高斯型随机变量。采样间隔较小时,独立性比较差。随着采样间隔的提高,传输矩阵的不同行变得更接近独立,因此相同采样数情况下的恢复效果变好。不同行变为统计独立后,继续增加采样间隔不再有明显效果。
光纤输出端光强的散斑大小是不同位置的光强相关性的一个直观衡量,可以借助散斑大小判断传输矩阵不同行间的统计独立程度。任选一张光纤输出端光强的图像,取中心一个30×30 pixel大小的区域的光强,平移这个区域,计算它与周围其他区域的光强的相关性,得到相关性与平均距离的关系。使用函数
$y= A \cdot \mathrm{e}\mathrm{x}\mathrm{p}\left(-\dfrac{{x}^{2}}{{b}^{2}}\right)+c$ 对实验数据进行拟合,得到的参数$ b $ 可以被视为散斑的大小,计算结果如图6所示。拟合使用了matlab中的fit函数,通过Trust-Region算法求解非线性最小二乘问题。从实验拟合数据可以得出参数$ b $ 为3.79 pixel。与图5图4相互印证,参数$ b $ 可用于判断采样间隔的适当大小。当采样间隔大于$ b $ 时,最少采样点数的变化趋于平缓。实际应用中关心的另一个重要因素是算法消耗的时间,而采样点数少并不等同于消耗的时间更少。相位恢复算法一轮迭代的时间随着采样点数量增加而增加,但算法收敛需要的迭代次数可能随采样点数量增加而减少,因此最终消耗的时间并不随采样点数量增加而单调变化。图7(a)中画出了重建图像与真实图像的相似度大于0.9的情况下算法消耗的时间。当采样间隔较小时,可以看到消耗的时间是先减小后增大的。可以预计,如果采样间隔较大时采样点数量也足够多,也应该存在一个最小值,只是由于采样点数量不足,目前无法下降到最小值。在采样间隔固定的情况下,可以找出相位恢复所需要的最短时间。从图7(b)中可以看到,在采样间隔较小时,最少时间随着采样间隔增大而变少。当采样间隔较大时,由于采样点不足,消耗的时间无法下降到应有的最小值,于是消耗的最短时间又会上升。从图7(a)中的趋势可以预测,如果采样点数量足够,则最短计算时间会趋于稳定。以本文的实验结构为例,如果以计算时间短为目标,则应该令采样间隔小于5 pixel,保证采样点数量足够的同时使采样间隔大于散斑大小。
Effect of sampling methods on distal field reconstruction through a multimode fiber (invited)
-
摘要: 利用相位恢复算法可以从光纤近端的光强分布求解光纤远端的场强分布。光纤的响应可以用传输矩阵描述。实验上则是在不同的输入情况下对输出端的光强分布进行足够数量的采样来测量传输矩阵。显然,采样点的位置分布,包括采样点数目和间隔,影响着传输矩阵的测量,而相位恢复算法的精度和效率与传输矩阵有关。文中提出采样间隔应该大于出射散斑大小,以满足传输矩阵不同行的统计独立性,在保证图像重建质量的条件下减少采样点数,提高重建效率。实验结果表明,当采样间隔小于散斑大小时,相同的图像重建质量下,随着采样间隔的增大,光场重建所需的采样点数量明显下降。当采样间隔大于散斑时,所需的采样点数量变化缓慢,约为输入图像像素数量的3.5倍。采样间隔固定时,随着采样点数的增加,相位恢复算法消耗的时间先减小后增大,因此存在一个最佳的采样间隔与采样点数。Abstract: The phase retrieval algorithms can be used to recover the field at the distal end of a fiber from the intensity at the proximal end of the fiber. The response of the fiber can be described by the transmission matrix. In the experiment, a sufficient number of samples are sampled from the output intensity distribution with different input conditions to measure the transmission matrix. Obviously, the position distribution of sampling points, including the sampling number and interval, affects the measurement of the transmission matrix, and the accuracy and efficiency of the phase retrieval algorithm are related to the transmission matrix. We propose that the sampling interval should be greater than the speckle size to ensure the independence of different rows of the transmission matrix; therefore, image quality can be guaranteed with fewer sampling points at higher reconstruction efficiency. The experimental results show that when the sampling interval is less than the speckle size, the number of sampling points required for light field reconstruction decreases significantly with increasing sampling interval under the same image reconstruction quality. When the sampling interval is greater than the speckle size, the number of sampling points required changes slowly and finally remains approximately 3.5 times the number of pixels of the input image. When the sampling interval is fixed, with the increase in sampling points, the time consumed by the phase retrieval algorithm first decreases and then increases, so there is an optimal sampling interval and sampling points.
-
Key words:
- scattering medium imaging /
- multimode fiber /
- phase retrieval
-
图 2 采样方案示意图,图中箭头表示向量
$ {\boldsymbol{e}}_{1},{\boldsymbol{e}}_{2} $ 的方向。(a)采样间隔为$ \sqrt{{1}^{2}+{3}^{2}}=3.16 $ pixel,$ {\boldsymbol{e}}_{1} $ =(3,1),$ {\boldsymbol{e}}_{2} $ = (−1,3);(b)采样间隔为3 pixel,$ {\boldsymbol{e}}_{1} $ =(3,0),$ {\boldsymbol{e}}_{2} $ =(0,3)Figure 2. Schematics of different sampling manners, the arrows in the figure represent the directions of the vectors
$ {\boldsymbol{e}}_{1},{\boldsymbol{e}}_{2} $ . (a) Sampling interval is$ \sqrt{{1}^{2}+{3}^{2}}=3.16 $ pixel,$ {\boldsymbol{e}}_{1} $ =(3,1),$ {\boldsymbol{e}}_{2} $ = (−1,3); (b) Sampling interval is 3 pixel,$ {\boldsymbol{e}}_{1} $ =(3,0),$ {\boldsymbol{e}}_{2} $ =(0,3)图 4 不同采样间隔下的结果。 (a)、(b)采样点数 2 400 时的重建图像和采样点分布; (c)、(d)采样点数 2 900 时的重建图像和采样点分布;(a)和(c)下方的数值为重建结果与真实图像的余弦相似度; (b)和(d)中灰色区域为光强包络区域,白点为采样点。右侧插入显示了放大的红色方块内的点分布
Figure 4. Results at different sampling intervals. (a), (b) Reconstruction results and the position distribution of sampling points when the number of sampling points is 2 400; (c), (d) Reconstruction results and the position distribution of sampling points when the number of sampling points is 2 900; Value below (a) and (c) is the cosine similarity between the reconstruction results to the ground truth; Grey area in (b) and (d) is the intensity envelope region, and the white dots are the sampling points. Zoomed-in view of the distribution of points in the red box is shown on the right
图 5 (a)不同采样间隔下相似度与采样点数的关系;(b)相同图像恢复质量下,最少采样点数随采样间隔的变化,虚线表示采样点数量上限;(c) 最小采样面积随采样间隔的变化,虚线表示光纤面积
Figure 5. (a) Similarity versus number of sampling points at different sampling intervals; (b) Variation of minimum number of sampling versus sampling interval under the same image quality, the dashed line represents the upper bound of the number of sampling points; (c) Variation of minimum sampling area versus sampling interval, the dashed line represents the area of the fiber
-
[1] Bianchi S, Di Leonardo R. A multi-mode fiber probe for holographic micromanipulation and microscopy [J]. Lab on a Chip, 2012, 12(3): 635-639. doi: 10.1039/C1LC20719A [2] Choi Y, Yoon C, Kim M, et al. Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber [J]. Physical Review Letters, 2012, 109(20): 203901. doi: 10.1103/PhysRevLett.109.203901 [3] Caravaca-Aguirre A M, Piestun R. Single multimode fiber endoscope [J]. Optics Express, 2017, 25(3): 1656. doi: 10.1364/OE.25.001656 [4] Oh G, Chung E, Yun S H. Optical fibers for high-resolution in vivo microendoscopic fluorescence imaging [J]. Optical Fiber Technology, 2013, 19(6): 760-771. doi: 10.1016/j.yofte.2013.07.008 [5] Ohayon S, Caravaca-Aguirre A, Piestun R, et al. Minimally invasive multimode optical fiber microendoscope for deep brain fluorescence imaging [J]. Biomedical Optics Express, 2018, 9(4): 1492. doi: 10.1364/BOE.9.001492 [6] Caramazza P, Moran O, Murray-Smith R, et al. Transmission of natural scene images through a multimode fibre [J]. Nature Communications, 2019, 10(1): 2029. doi: 10.1038/s41467-019-10057-8 [7] Zhao T, Ourselin S, Vercauteren T, et al. Seeing through multimode fibers with real-valued intensity transmission matrices [J]. Optics Express, 2020, 28(14): 20978-20991. doi: 10.1364/OE.396734 [8] Chandra R, Goldstein T, Studer C. PhasePack: A phase retrieval library[C]//2019 13th International conference on Sampling Theory and Applications (SampTA), 2019: 1-5. [9] Fannjiang A, Strohmer T. The numerics of phase retrieval [J]. Acta Numerica, 2020, 29: 125-228. doi: 10.1017/S0962492920000069 [10] Candès E J, Strohmer T, Voroninski V. PhaseLift: Exact and stable signal recovery from magnitude measurements via convex programming [J]. Communications on Pure and Applied Mathematics, 2013, 66(8): 1241-1274. doi: 10.1002/cpa.21432 [11] Goldstein T, Studer C. PhaseMax: Convex phase retrieval via basis pursuit [J]. IEEE Transactions on Information Theory, 2018, 64(4): 2675-2689. doi: 10.1109/TIT.2018.2800768 [12] Candes E J, Li X, Soltanolkotabi M. Phase retrieval via Wirtinger flow: Theory and algorithms [J]. IEEE Transactions on Information Theory, 2015, 61(4): 1985-2007. doi: 10.1109/TIT.2015.2399924 [13] Netrapalli P, Jain P, Sanghavi S. Phase retrieval using alternating minimization [J]. IEEE Transactions on Signal Processing, 2015, 63(18): 4814-4826. doi: 10.1109/TSP.2015.2448516 [14] Schniter P, Rangan S. Compressive phase retrieval via generalized approximate message passing [J]. IEEE Transactions on Signal Processing, 2015, 63(4): 1043-1055. doi: 10.1109/TSP.2014.2386294 [15] Dremeau A, Krzakala F. Phase recovery from a Bayesian point of view: The variational approach[C]//2015 IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), 2015: 3661-3665. [16] Shechtman Y, Eldar Y C, Cohen O, et al. Phase retrieval with application to optical imaging: A contemporary overview [J]. IEEE Signal Processing Magazine, 2015, 32(3): 87-109. doi: 10.1109/MSP.2014.2352673 [17] Drémeau A, Liutkus A, Martina D, et al. Reference-less measurement of the transmission matrix of a highly scattering material using a DMD and phase retrieval techniques [J]. Optics Express, 2015, 23(9): 11898-11911. doi: 10.1364/OE.23.011898 [18] Rajaei B, Tramel E W, Gigan S, et al. Intensity-only optical compressive imaging using a multiply scattering material and a double phase retrieval approach[C]//2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2016: 4054-4058. [19] Lee K, Park Y. Exploiting the speckle-correlation scattering matrix for a compact reference-free holographic image sensor [J]. Nature Communications, 2016, 7(1): 13359. doi: 10.1038/ncomms13359 [20] N’Gom M, Lien M B, Estakhri N M, et al. Controlling light transmission through highly scattering media using semi-definite programming as a phase retrieval computation method [J]. Scientific Reports, 2017, 7(1): 2518. doi: 10.1038/s41598-017-02716-x [21] Deng L, Yan J D, Elson D S, et al. Characterization of an imaging multimode optical fiber using a digital micro-mirror device based single-beam system [J]. Optics Express, 2018, 26(14): 18436-18447. doi: 10.1364/OE.26.018436 [22] Huang G, Wu D, Luo J, et al. Retrieving the optical transmission matrix of a multimode fiber using the extended Kalman filter [J]. Optics Express, 2020, 28(7): 9487-9500. doi: 10.1364/OE.389133 [23] Shlezinger N, Dabora R, Eldar Y C. Measurement matrix design for phase retrieval based on mutual information [J]. IEEE Transactions on Signal Processing, 2018, 66(2): 324-339. doi: 10.1109/TSP.2017.2759101 [24] Gao B, Liu H, Wang Y. Phase retrieval for sub-Gaussian measurements [J]. Applied and Computational Harmonic Analysis, 2021, 53: 95-115. doi: 10.1016/j.acha.2021.01.001 [25] Panarin S, Müller J, Prabhakar S, et al. Spatial structuring of light for undergraduate laboratories [J]. American Journal of Physics, 2021, 89(2): 210-219. doi: 10.1119/10.0002365