Imaging through multimode fiber (MMF) provides high-resolution imaging through a fiber with cross section down to tens of micrometers. It requires interferometry to measure the full transmission matrix (TM), leading to the drawbacks of complicated experimental setup and phase instability. Reference-less TM retrieval is a promising robust solution that avoids interferometry, since it recovers the TM from intensity-only measurements. However, the long computational time and failure of 3D focusing still limit its application in MMF imaging. We propose an efficient reference-less TM retrieval method by developing a nonlinear optimization algorithm based on fast Fourier transform (FFT). Furthermore, we develop an algorithm to correct the phase offset error of retrieved TM using defocused intensity images and hence achieve 3D focusing. The proposed method is validated by both simulations and experiments. The FFT-based TM retrieval algorithm achieves orders of magnitude of speedup in computational time and recovers 2286 × 8192 TM of a 0.22 NA and 50 μm diameter MMF with 112.9 s by a computer of 32 CPU cores. With the advantages of efficiency and correction of phase offset, our method paves the way for the application of reference-less TM retrieval in not only MMF imaging but also broader applications requiring TM calibration. |
1.IntroductionImaging through multimode fibers (MMFs) of tens to hundreds of micrometers enables high-resolution imaging by a hair-thin instrument. It provides minimally invasive high-resolution imaging for locations deep inside living organisms1 without traumatic tissue slices.2 Its broad applications include in vivo endoscopes,3–5 optical tweezers over cellular area,6,7 and remote time-of-flight 3D depth sensing.8 MMF imaging is achieved by exploiting the property of transmission matrix (TM). With the TM, one can collect the feedback signal after rapidly scanning foci on the sample5,9 or directly inverse the scattering process.10,11 However, the calibration of TM requires measuring the transmitted complex fields after sending probing incident complex fields. With both amplitude and phase, the complex fields cannot be measured directly by a camera. Conventionally, external reference methods12–14 build a complicated experimental setup to interferometrically measure the transmitted complex field with an external reference beam. These methods suffer from phase instability of the reference beam, easily caused by mechanical variation and thermal drift. The internal reference methods9,15,16 set parts of the modulation modes as an internal reference. It reduces the number of effective modulation modes and uses speckle reference, which contains dark reference points.1,16 Its retrieved TM has the phase offset error, excluding applications that require 3D focusing.6 Therefore, it is desirable to develop a simple and stable method to measure the full TM in MMF imaging. The reference-less TM retrieval methods computationally recover the TM from the intensity images of the transmitted complex fields measured by a reference-less experimental setup17 [Fig. 1(a)]. This provides a promising solution to robustly measure the TM in MMF imaging. However, the application of reference-less TM retrieval in MMF imaging still faces issues of long computational time and phase offset error. Algorithms based on the Bayesian approach,17,18 Gerchberg–Saxton,19 semi-definite programming,20 Kalman filter,21 prVAMP,22 and reweighted amplitude flow23 have been proposed for reference-less TM retrieval. When the TM is large, the computational time of these algorithms could be hours.19,22 Another issue is that the TM acquired by the reference-less TM retrieval has the error of phase offset.6 The reference-less TM retrieval only maps the incident fields with the intensity images at one fixed camera plane, causing the loss of relative phase information between different pixels. The transmitted complex field predicted using the acquired TM with the phase offset error has correct amplitudes but wrong phases. It causes failure of generating 3D light patterns, such as 3D foci or light sheet, which is essential in volumetric imaging24 and light-sheet imaging by MMF.25 Here we propose a fast and phase-offset-error-free reference-less TM retrieval method and demonstrate with experimental setup of MMF imaging. First, we design the probing incident complex fields with Fourier transform matrix and develop a nonlinear optimization algorithm to solve the TM from the intensity-only measurements. Compared to the scheme of random probing incident fields,19–21 our scheme allows the inverse algorithm to be implemented with FFT [Fig. 1(b)], which greatly reduces the computational complexity. Second, our method measures a set of intensity images at a defocus plane from the distal end of the MMF and develops an algorithm to recover the phase offset from the defocus intensity images [Fig. 1(c)]. The simulation shows the TM retrieval algorithm with FFT has a speedup in computational time compared to that of the TM retrieval algorithm without FFT. The proposed TM retrieval method recovers TM for an MMF of 0.22 NA and diameter with 112.9 s. We build the experimental setup with MMF and verify the proposed methods by evaluating the foci in both 2D and 3D. In addition to MMF imaging, TM measurement is important to applications that require characterizing the scattering property. The applications range from imaging through scattering media or tissue,26 high-capacity communication through optical fiber,27 optical computing systems,28,29 quantum networks,30,31 and thin-lens imaging by nano-optics,32 to holographic display with scattering media.33 As an alternative to the interferometry-based external reference methods, reference-less TM retrieval provides a simple and robust way to acquire the TM in these applications. However, the issues of long computational time and phase offset error commonly exist, especially when the TM is large. Our method may be adapted to solve the issues of reference-less TM retrieval for these applications. 2.Methods2.1.Experimental SetupThe reference-less experimental setup is shown in Fig. 2. An MMF of 0.22 NA and diameter (ChunHui CCS50/125H-F-F-1) is used. It aims to generate incident complex fields impinging on the proximal end of the MMF and measure the intensity images of the transmitted complex fields of the MMF. A collimated laser of 488 nm (Precilasers SF-488-0.5-CW) is directed on a digital micromirror device (DMD, Vialux v9501). It displays the binary hologram obtained by the Lee hologram method.34 A set of half-wave plates and quarter-wave plates interposed between the DMD and the polarization beam splitter PBS1 turns the light reflected from the DMD into circular polarized light. The s and p lights from PBS1 pass through the mirrors M2 and M3 separately, combine by the polarization beam splitter PBS2, and impinge on the iris after being Fourier-transformed by lens L3. The s and p lights are programmed with different carrier frequencies on the binary hologram, which determines the locations the diffraction order on the iris. Tuning M2 and M3 shifts the diffraction order of the s and p lights through the pinhole on the iris. The telescope system formed by lens L4 and the objective lens OBJ1 (Olympus NA 0.25) focuses the light emitted from the iris on the proximal end of the MMF. Thus the incident complex fields with desired phases are generated for both polarizations. Dual polarization modulation allows better control of the transmitted complex fields, resulting in scanning foci of improved quality for MMF imaging.13 The 4f system formed by the objective lens OBJ2 (Olympus NA 0.25) and the lens L5 magnifies the transmitted complex field. A CMOS camera (Basler acA720-520um) captures the intensity image after the light passes through a linear polarizer. The synchronization of the DMD (refresh rate of 16.7 kHz) and the CMOS camera enables measuring the intensity images at a high frame rate up to . By moving the OBJ2 with a stage (Thorlabs CT1P), the CMOS camera captures the intensity images at a defocus plane. 2.2.TM RetrievalOur method recovers the TM from the intensity images measured from the reference-less experimental setup in Fig. 2. The data-acquisition measures the intensity images of the transmitted complex fields for TM retrieval after sending in the phase-only incident complex fields generated by modulating the DMD. Next, we develop an efficient nonlinear optimization algorithm to recover the TM from the intensity measurements. 2.2.1.Pixel-wise inverse problemThe intensity images are denoted as , where , is the total number of the measured intensity images, and are the spatial coordinates. Each image contains pixels. The incident complex field is raster-scanned into a vector , which has modulation modes for both polarizations. The forward model of the intensity measurement is written as where is a vector raster-scanned from , is the transmission matrix (), and takes absolute square of the complex numbers inside.The optimization problem of retrieving the TM from the intensity images is formulated as where is the squared L2 norm of the vector inside. The cost function is the sum of the squared error between the intensity measurements and the intensity predicted by the forward model. The optimization problem solves by minimizing the cost function.The large number of unknown in makes the optimization problem in Eq. (2) difficult to solve directly. However, it can be broken down into smaller optimization problems. Each problem is formulated based on the intensity measurement at one single pixel, whereHere the column vector is the transpose of the row of , is the element of , denotes transpose, and . The vector contains all of the measurements at the same pixel indexed by on the intensity images. The matrix is the so-called probing matrix; each row of is one of the incident fields. Each optimization problem recovers one row of from the measurements at the corresponding pixel. Thus the whole TM can be recovered by solving these small pixel-wise optimization problems independently. 2.2.2.Inverse algorithm with a designed probing matrixConventional reference-less TM retrieval methods19–21 set the phases of the probing matrix as random numbers. It means that the incident fields are modulated with random phases in these methods. By contrast, our method designs the matrix with Fourier transform matrix, where is the Fourier transform matrix, is a by 1 column vector with its phase set as random numbers, is a diagonal matrix with its diagonal entries from , and . The incident complex field has modulation modes for each polarization, so we have . The total number of measured intensity images is , which is -fold of the size of the unknown complex variable, . We set , since the complex variable in the inverse problem contains both real and imaginary parts. The matrix consists of square matrix in the form of , where is the 2D Fourier transform matrix for matrix. Since is a pure phase matrix, the probing matrix remains as pure phase. So it can be loaded into the phase modulator to generate desired incident fields. When the phases of the probing matrix are random,19–21 the multiplication of with a vector has to be computed with matrix–vector multiplication. Our method designs the probing matrix with the Fourier transform matrix, so the matrix–vector multiplication related to can be computed with FFT or inverse FFT (for its complex transpose ). This advantage can be exploited to accelerate the algorithm for TM retrieval.The optimization problem in Eq. (3) recovers the complex unknown from the intensity measurement and can be viewed as a phase-retrieval problem. We solve the optimization problem using the nonlinear optimization phase-retrieval method,35 which has been developed for phase retrieval from defocus intensity images. It derives the derivative for the complex variable and adopts the limited memory Broyden–Fletcher–Goldfarb—Shanno (L-BFGS) method36,37 to solve the phase-retrieval problem. The optimization is initialized by backpropagation, where takes the element-wise square root of the vector , and denotes the complex transpose. Since is the inverse Fourier transform matrix, the matrix–vector multiplication in can be computed with FFT. We derive the first derivative of with respect to as (more details in Appendix A) where and are the real and imaginary parts of . The matrix–vector multiplication related to in Eq. (7) can be computed with FFT.The procedure of the algorithm to solve the optimization problem in Eq. (3) is summarized in Algorithm 1. The algorithm has inputs of the intensity measurements at th pixel and random phase vectors for the probing matrix . It recovers one row of the TM, . The estimation is initialized by Eq. (6). After obtaining the error [Eq. (3)] and gradient [Eq. (7)], the algorithm updates the estimate of iteratively according to the L-BFGS method. The update iteration stops when a preset maximum iteration number is reached. The matrix–vector multiplication related to in Eqs. (3), (6), and (7) can be efficiently computed with FFT, reducing the computational complexity from to . It also has the benefit of memory efficiency, since there is no need to explicitly store the big matrix when solving the inverse problem. 2.2.3.Recovery of the whole TMOne may apply the method in Algorithm 1 on all the optimization problems in the form of Eq. (3) and recover the entire TM. However, this could be unnecessary due to the physical properties of the MMF. There is negligible transmitted light on the pixels outside the distal end of the fiber. The complex field at the distal end of the MMF has highest frequency limited by , where NA is the numerical aperture of the MMF, and is the wavelength. Therefore, we design a preprocessing procedure to reduce the number of effective pixels, which in turn brings down the number of optimization problems. First, we half-sample the measured intensity images by only keeping the pixels of the odd indices in the images. Without loss of generality, we assume that the pixel size of the measured intensity images, , meets the Nyquist sampling theory, Algorithm 1Optimization of recovering a row of TM tk.
So, the pixel size of the half-sampled images has . It meets the sampling requirement of the transmitted complex field, The half-sample step reduces the number of the optimization problems by a factor of 4. Second, from the half-sampled intensity images, we obtain a binary fiber mask, which masks out the MMF region. The pixels within the distal end of the MMF have a value of 1, while the pixels outside of the MMF have a value of 0. For the pixels in the zero-value region, the vectors are directly set to zeros, without solving the optimization problem of Eq. (3). At the end, our TM retrieval method only solves the optimization problems from the half-sampled intensity images for the pixels inside the distal end of the fiber. Thus the number of optimization problems is greatly reduced.The full procedure of the TM retrieval method is summarized in Fig. 3. The optimization problems of Eq. (3) are independent. So our method solves these optimization problems in parallel with a computer of multiple CPU cores, which can be easily implemented by parfor on MATLAB. 2.3.Phase CorrectionThe optimization problem in Eq. (3) has the issue of phase ambiguity. It means multiplying an optimal solution of the optimization problem with an arbitrary phase term, still gives an optimal solution. The phase ambiguity leads to the error of phase offset between the estimated TM and the true TM. It has where is the true TM, is the estimated TM, and the vector is the phase offset (size ). Note that and are the size of the half-sampled intensity image, due to the preprocessing step in Fig. 3.In order to generate foci at the measurement plane of the intensity images, the incident complex fields are modulated with phases of the conjugate of the estimated TM row by row. The phase offset error of still allows pixel-wise constructive interference, resulting in the generation of 2D scanning foci at the measurement plane. However, to generate foci at other focal planes, the incident complex fields are modulated with a free-space-propagated version of . The phase offset error causes random error in the propagated TM, leading to failure of generating 3D distributed foci at other focal planes. To solve the issue of phase ambiguity, we propose a phase correction method after has been obtained by the TM retrieval method in Fig. 3. Our method corrects the phase offset using multiple defocus intensity images. After applying random phases to modulate the incident fields, our method measures intensity images at a defocus plane that is away from the distal end of the MMF (). The defocus intensity images results from free-space propagation of the transmitted complex field at the distal end of the fiber. These defocus intensity images capture the phase information of the transmitted complex fields. Therefore, it is possible to invert the phase offset of the estimated TM from the defocus intensity images. We build the forward model for the inverse problem of the phase offset recovery from the defocus intensity images. The defocus intensity images are denoted with , where . Each intensity image has a size of with pixel size of . From Eq. (11), the transmitted complex field at the distal end of the MMF can be expressed as where is raster-scanned from the corresponding incident complex field of the defocus intensity image, . Note that the transmitted complex field predicted by the estimated TM has a size of with pixel size of . The vector is the raster-scanned form of the transmitted complex field.The complex field at the defocus plane and the transmitted complex field at the distal end of the MMF are related by defocus propagation. The forward model of the defocus intensity can be expressed as where the vector is raster-scanned from , is the Fourier transform matrix for the matrix, is the defocus propagation kernel (more in Appendix B), is for zero padding in the frequency domain, and is the inverse Fourier transform matrix for the matrix. The measured intensity has a size of , whereas the transmitted complex field has a size of . The zero padding here adds zeros in the frequency domain, which results in doubling the number of pixels in both dimensions. It has the effect of reversing the half-sample step in Fig. 3.The optimization problem of solving the phase offset from the defocus intensity images is formulated as where . The cost function is defined as the squared error between the measured defocus intensity and the intensity predicted with the phase offset.The first derivative of the optimization problem in Eq. (14) is derived as More details can be found in Appendix C. The matrix–vector multiplication related to in Eqs. (14) and (15) can be computed with FFT, without explicitly forming the big matrices. With the cost function in Eq. (14) and the first derivative in Eq. (15), our method uses the L-BFGS method36,37 to recover the phase offset from the defocus intensity images. 3.Results for TM Retrieval Using Intensity Images at One Measurement PlaneIn this section, we verify the TM retrieval algorithm using intensity images measured at one fixed plane (Fig. 3) by both simulations and experiments. In the simulation, a TM of size was used to generate simulated data. The TM had been measured experimentally by the off-axis holography method13 with an external reference beam, for an MMF of 0.22 NA and diameter. The incident complex fields had phase modulation modes for each polarization. The matrix in the probing matrix was set as the Fourier transform matrix for a matrix. The transmitted complex fields at the distal end were sampled with . We simulated seven datasets with to 9, where is the total number of in Eq. (5). We ran the TM retrieval algorithm on each of the simulated datasets. The entire TM was recovered by applying the method in Algorithm 1 on all of the 9216 pixels, without the preprocessing step in Fig. 3. For each dataset, the optimization problems were solved in a parallel manner on a computer with 32 CPU cores (Intel Xeon Gold 5218 2.3 GHz). For the dataset of , it took 376.4 s to retrieve the entire TM of size . Figure 4 compares the error of the recovered TM using the datasets of different measurement sizes. It shows the root-mean-square error (RMSE) of both amplitude [Fig. 4(a)] and phase [Fig. 4(b)] of the recovered TM. Each row of the recovered TM is compared with its true value, and the errors of all of the 9216 rows are organized in grids, which are shown in Fig. 4. The phase error is obtained by subtracting the phase of each row of the recovered TM with the true values after removing the constant phase offset [Eq. (11)]. For , most rows of the recovered TM have large errors. For to 6, a few of the rows of the recovered TM have large errors; there are random bright spots (meaning large errors) in the images at the top row of Figs. 4(a) and 4(b). However, these speckles disappear as increases. The error of the recovered TM becomes negligibly small for to 9. The plots on the bottom left on Figs. 4(a) and 4(b) show the RMSE of both amplitude and phase of the recovered TM converge to zero for to 9. The simulation demonstrates the proposed TM retrieval method is able to efficiently recover the TM from the intensity images measured at one imaging plane with negligible errors. Table 1 shows the speedup of computational time by the proposed TM retrieval with FFT. The central of the of the dataset of were used to access the computational time of the TM retrieval algorithms. The TM retrieval algorithm without FFT replaces the FFT in Algorithm 1 with matrix–vector multiplication. The TM retrieval algorithm without FFT recovers the TM with 43,664.1 s (12.1 h). However, the proposed TM retrieval algorithm implemented with FFT recovers the same-sized TM with 35.4 s. For the proposed algorithm, each row of the TM takes 0.035 s on average. Using FFT, the proposed TM retrieval algorithm achieves 1200× speedup. Table 1TM retrieval algorithm with FFT achieves 1200× speedup.
In the experiment, we used an MMF of 0.22 NA and diameter. The illumination was laser of 488 nm. The DMD achieved phase modulation for each polarization, resulting in 8192 modes in total. We tested the TM retrieval algorithm for the cases of to 8. Each case followed the procedure in Fig. 3 to recover the TM. For each case, we generated the probing matrix with random phase vectors and Fourier transform matrix by Eq. (5). The matrix was set as the Fourier transform matrix for a matrix. The phase of the probing matrix was loaded into the DMD, and a series of images [Fig. 5(a)] were measured. Each image has , with a pixel size of . The preprocessing step half-sampled the measured images and obtained images of [Fig. 5(b)]. From the sum image of the preprocessed images, we calculated the fiber mask [Fig. 5(c)], which covers 99.9% of the total energy. The white region of the mask indicates the distal end of the MMF fiber. Only for the pixels inside the fiber mask, the TM was retrieved by Algorithm 1 from the preprocessed images of each dataset. Here we give an example of the case of . The probing matrix has a size of . The preprocessed intensity images contain 57,344 images of size . The number of pixels inside the white region of the fiber mask is 2286, so the retrieved TM has size of . For each selected pixel, an optimization problem in the form of Eq. (3) is formulated; it has inputs of the intensity measurement at the corresponding pixel (a vector of ) and the random phase vectors used to generate the probing matrix. All these 2286 optimization problems were solved in parallel on the computer with 32 CPU cores. For the case of , the computer takes 112.9 s to solve the optimization problems in TM retrieval. The accuracy of the recovered TM was tested by the ability to generate foci at the measurement plane. After the TM was retrieved for each case, we uploaded the phase of the conjugate complex of the recovered TM into the phase modulator, sequentially modulated the incident field with its phase row by row, and measured the generated intensity images. When the displayed phase of a row of the retrieved TM matches with the true TM of the imaging system, a focus is generated at the camera. In order to evaluate the quality of the focus, we measured two images () for each focus with an exposure time of 70 and (over exposure at the focus) and calculated the power ratio (PR) of the focus by combining these two images. The PR is the ratio of the signal to the total energy.1 The signal is the sum of the near the peak of the focus using the image, while the total energy is the sum of the signal and the background (outside the ), which is calculated using the image. The PR reflects the quality of the focus, and hence experimentally shows the correctness of the retrieved TM. We measured a TM by the off-axis holography method with an external reference beam and acquired the corresponding focal images. The result by the holography method acts as a reference for our method. Figure 5(d) shows the PR of the foci of cases of different and the holography method. For the cases to 6, there are several foci that have a low PR. However, for the cases of , the overall quality of the foci is near to that of the holography method. The average PR of the case of is 0.64, which is slightly lower than that of the case of holography (0.651). Figures 5(e) and 5(f) further compare the cases of and holography by showing the distribution of the PR and a sum projection of several selected foci. The TM retrieval method by has more foci of PR above 0.60 than that of the holography method. And hence, the accuracy of the TM recovered by our proposed reference-less method is validated by comparing with the holography method. 4.Results for TM Retrieval with Phase CorrectionIn this section, we validate the TM retrieval algorithm with phase correction. For simulation, we used a simulated TM of size for an MMF with 0.22 NA and diameter, generated by solving Maxwell’s equations.38 The transmitted complex fields of the MMF are sampled by grids with a pixel size of , and the wavelength of illumination is 532 nm. We designed a probing matrix with 2D Fourier transform matrix for a matrix and . A series of 82,944 images of size were generated at the distal end of the fiber (). Then we simulated 50 defocus images [Fig. 6(a)] at away from the distal end of the MMF. Each image has with a pixel size of . The incident complex fields were obtained by 50 random phases, and the defocus intensity images were generated by Eq. (13). We first followed the preprocessing step and the optimization step of the procedure in Fig. 3 to recover the TM from the simulated data. In the preprocessing step, the half-sample step was not performed, since the pixel size already meets the sampling requirement of the complex field. A fiber mask was generated, resulting 6668 selected pixels inside the white region. For the selected pixels, the optimization problems in the form of Eq. (3) were solved, and the rows of the recovered TM corresponding to the black region in the mask were set to zeros. Thus a recovered TM was obtained but has the error of phase offset, since the measured intensity images were at one fixed plane. Next, the phase offset was solved from the defocus images and the recovered TM by the phase-correction algorithm. The computational time for the TM retrieval and the phase correction was 332.9 and 85.2 s, respectively. Figure 6(b) shows the recovered phase offset by the algorithm. Finally, we compensated the phase offset error of the recovered TM using the recovered phase. The amplitude and phase RMSE of the recovered TM with phase correction is and , respectively. The error between the recovered TM with phase correction and the true TM is small, as shown in Fig. 6(c). We further validated the TM algorithm with phase correction by experimentally displaying 3D foci. In the experiment, we used an MMF of diameter and 0.22 NA, and illumination wavelength of 488 nm. The phase modulation on DMD had modes for each polarization. We designed a probing matrix using Fourier transform matrix for a matrix, and . After modulating the DMD with the phase of the probing matrix, we sequentially measured 65,536 intensity images at the distal end of the MMF (). Each image has pixels with a pixel size of . In order to correct the phase offset, we measured 50 images of at away from the distal end [Fig. 7(a)]. The defocus images were measured after applying 50 random phases on the phase modulator. We first recovered a TM from the intensity images measured at by the proposed method in Fig. 3. In the preprocessing step, we half-sampled the images to a size of and generated a fiber mask, which has 3015 pixels inside the white region of the mask. By solving the optimization problems, the TM retrieval algorithm obtained a TM with the error of the phase offset. Next, the algorithm of phase correction recovered the phase offset [Fig. 7(b)] from the defocus intensity images. The recovered phase offset was used to correct the error of phase offset in the recovered TM. The computational time for the algorithm of TM retrieval and the algorithm of phase correction were 199.3 and 20.6 s, respectively. We tested the recovered TM by generating 3D foci on the imaging system. The propagated TM at a defocus distance could be obtained by adding the recovered TM at with a free-space defocus propagation. We generated the two sets of propagated TMs at , and , using the recovered TM with the error of phase offset and the recovered TM with phase correction. We sequentially applied the phases of complex conjugate of the propagated TM to the DMD and measured intensity images at the corresponding defocus distances. Figure 7(c) compares the intensity images measured at different defocus distances for the foci at the center of the images. For the case of the TM with phase error, the focus could be seen at the image center for , but it quickly scattered into random patterns at other defocus distances [top row of Fig. 7(c)]. The phase offset error causes the failure in generating 3D foci. In contrast, the propagated TM generated using the recovered TM with phase correction successfully generated the foci at defocus distances [bottom row of Fig. 7(c)]. As the defocus distances increase from 0 to , the PRs of the foci reduce from 0.60 to 0.51. The decrease of focus brightness could be caused by the defocus propagation. It adds more correlation for the rows of propagated TM corresponding to the neighborhood pixels. Figure 8 shows the sum projection of selected foci at different defocus distances generated by the recovered TM with phase correction. This validates the algorithm of the TM retrieval with phase correction. 5.DiscussionThe optimization problem in Eq. (3) is a phase-retrieval problem. The cost function of the phase-retrieval problem is formulated based on intensity difference, which is suitable for the assumption that the intensity measurements are polluted by Gaussian noise. With the assumption of Poisson noise, the cost function can be formulated with amplitude difference.39 Many algorithms have been proposed for the phase-retrieval problem, including gradient descent,40 Gerchberg–Saxton,41 Kalman filtering,42 L-BFGS,35,43 modified Gauss Newton,35 Wirtinger flow,44 prVBEM,45 PhaseLift,46 reweighted amplitude flow,47 and PhaseMax.48 The L-BFGS method is a second-order optimization method, which was shown to converge faster than the first-order methods, such as gradient descent or Gerchberg–Saxton in phase retrieval from defocus images35 and Fourier ptychography.39 In this work, we used the intensity-based cost function and the L-BFGS method. A fair assessment of the formulation of the cost function and the optimal choice of the algorithm for the phase-retrieval problem in the TM retrieval is outside the scope of this work. This work proposed to design the probing matrix () with the Fourier transform matrix. Using FFT, the computational complexity of the matrix–vector multiplication related to and reduces from to . Here we give an example of the number of modulation modes and the number of measurements . The matrix has a size of . The computational complexity reduces from to , and it is memory-efficient without storing . The computation related to the probing matrix is mostly inevitable in the algorithms of the phase-retrieval problem in TM retrieval. For example, gradient descent-based algorithms have to compute the cost function and gradient descent. The computational complexity of these algorithms is lower-bounded by , due the matrix–vector multiplication related to . It is higher than that of our proposed method using FFT . However, applying the similar FFT-based scheme in these algorithms could further reduce the computational complexity. 6.ConclusionWe have proposed a fast and phase-offset-error-free method for reference-less TM retrieval. By designing the probing incident complex fields with a Fourier transform matrix, the FFT-based TM retrieval algorithm achieves orders of magnitude of improvement in computational complexity [] versus []. Further, the algorithm corrects the error of phase offset in the TM retrieval using the defocus intensity images. We have tested the proposed TM retrieval method by both simulations and experiments with MMF. It has been experimentally demonstrated speedup in recovering the TM of size with 112.9 s for the MMF of 0.22 NA and diameter by the computer of 32 CPU cores. This work validated the method by evaluating 2D and 3D foci with the experimental setup using MMF. With the advantage of computational efficiency and the correction of phase offset, the method may be adapted to efficiently calibrate the TM of scattering media or imaging systems without reference for applications, such as diffusers,22 optical communication,49 3D holography displays,33 quantum networks,30,31 optical computing system,28 and thin lens imaging system.32 7.Appendix A: Derivation of the First Derivative in the TM RetrievalThe optimization problem to recover one row of the TM is expressed as Next, we define where is a vector. According to the chain rule, the first derivative of the cost function can be written asThus we have the Hermitian of the first derivative as 8.Appendix B: Defocus PropagationAccording to the theory of angular spectrum propagation,50 the defocus propagation kernel in frequency domain is expressed as where is the wavelength of the illumination, and are the spatial frequency coordinates, and is the pupil of the imaging system. The pupil is written asThe vector is raster-scanned from . 9.Appendix C: Derivation of the First Derivative in the Algorithm of Phase CorrectionThe optimization of solving the phase offset from the defocus intensity images is rewritten as where . Next, we defineUsing the chain rule, we have We can have By combining Eqs. (26) and (27), we can get It is easy to obtain So, we have Data AvailabilityAdditional data related to this study are available from the corresponding authors upon reasonable request. AcknowledgmentsThis work was supported by the National Natural Science Foundation of China (Grant Nos. T2293751, T2293752, 61735017, 62020106002, and 62005250), the National Key Basic Research Program of China (Grant No. 2021YFC2401403), and the Major Scientific Research Project of Zhejiang Lab (Grant No. 2019MC0AD02). ReferencesS. Turtaev et al.,
“High-fidelity multimode fibre-based endoscopy for deep brain in vivo imaging,”
Light Sci. Appl., 7
(1), 92 https://doi.org/10.1038/s41377-018-0094-x
(2018).
Google Scholar
Z. Wen et al.,
“Single multimode fibre for in vivo light-field-encoded endoscopic imaging,”
Nat. Photonics, 17
(8), 679
–687 https://doi.org/10.1038/s41566-023-01240-x NPAHBY 1749-4885
(2023).
Google Scholar
T. Čižmár and K. Dholakia,
“Exploiting multimode waveguides for pure fibre-based imaging,”
Nat. Commun., 3
(1), 1027 https://doi.org/10.1038/ncomms2024 NCAOBW 2041-1723
(2012).
Google Scholar
I. N. Papadopoulos et al.,
“High-resolution, lensless endoscope based on digital scanning through a multimode optical fiber,”
Biomed. Opt. Express, 4
(2), 260
–270 https://doi.org/10.1364/BOE.4.000260 BOEICL 2156-7085
(2013).
Google Scholar
A. M. Caravaca-Aguirre and R. Piestun,
“Single multimode fiber endoscope,”
Opt. Express, 25
(3), 1656
–1665 https://doi.org/10.1364/OE.25.001656 OPEXFF 1094-4087
(2017).
Google Scholar
S. Bianchi and R. Di Leonardo,
“A multi-mode fiber probe for holographic micromanipulation and microscopy,”
Lab Chip, 12
(3), 635
–639 https://doi.org/10.1039/C1LC20719A LCAHAM 1473-0197
(2012).
Google Scholar
I. T. Leite et al.,
“Three-dimensional holographic optical manipulation through a high-numerical-aperture soft-glass multimode fibre,”
Nat. Photonics, 12
(1), 33
–39 https://doi.org/10.1038/s41566-017-0053-8 NPAHBY 1749-4885
(2018).
Google Scholar
D. Stellinga et al.,
“Time-of-flight 3D imaging through multimode optical fibers,”
Science, 374
(6573), 1395
–1399 https://doi.org/10.1126/science.abl3771 SCIEAS 0036-8075
(2021).
Google Scholar
S. M. Popoff et al.,
“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 https://doi.org/10.1103/PhysRevLett.104.100601 PRLTAO 0031-9007
(2010).
Google Scholar
S. Popoff et al.,
“Image transmission through an opaque material,”
Nat. Commun., 2
(6), 81 https://doi.org/10.1038/ncomms1078 NCAOBW 2041-1723
(2010).
Google Scholar
Y. Choi et al.,
“Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber,”
Phys. Rev. Lett., 109
(20), 203901 https://doi.org/10.1103/PhysRevLett.109.203901 PRLTAO 0031-9007
(2012).
Google Scholar
Y. Choi et al.,
“Overcoming the diffraction limit using multiple light scattering in a highly disordered medium,”
Phys. Rev. Lett., 107
(2), 023902 https://doi.org/10.1103/PhysRevLett.107.023902 PRLTAO 0031-9007
(2011).
Google Scholar
T. Čižmár and K. Dholakia,
“Shaping the light transmission through a multimode optical fibre: complex transformation analysis and applications in biophotonics,”
Opt. Express, 19
(20), 18871
–18884 https://doi.org/10.1364/OE.19.018871 OPEXFF 1094-4087
(2011).
Google Scholar
S. Li et al.,
“Compressively sampling the optical transmission matrix of a multimode fibre,”
Light Sci. Appl., 10
(1), 88 https://doi.org/10.1038/s41377-021-00514-9
(2021).
Google Scholar
D. B. Conkey, A. M. Caravaca-Aguirre and R. Piestun,
“High-speed scattering medium characterization with application to focusing light through turbid media,”
Opt. Express, 20
(2), 1733
–1740 https://doi.org/10.1364/OE.20.001733 OPEXFF 1094-4087
(2012).
Google Scholar
J. Yoon et al.,
“Measuring optical transmission matrices by wavefront shaping,”
Opt. Express, 23
(8), 10158
–10167 https://doi.org/10.1364/OE.23.010158 OPEXFF 1094-4087
(2015).
Google Scholar
A. Drémeau et al.,
“Reference-less measurement of the transmission matrix of a highly scattering material using a DMD and phase retrieval techniques,”
Opt. Express, 23
(9), 11898
–11911 https://doi.org/10.1364/OE.23.011898 OPEXFF 1094-4087
(2015).
Google Scholar
L. Deng et al.,
“Characterization of an imaging multimode optical fiber using a digital micro-mirror device based single-beam system,”
Opt. Express, 26
(14), 18436
–18447 https://doi.org/10.1364/OE.26.018436 OPEXFF 1094-4087
(2018).
Google Scholar
G. Huang et al.,
“Generalizing the Gerchberg–Saxton algorithm for retrieving complex optical transmission matrices,”
Photonics Res., 9
(1), 34
–42 https://doi.org/10.1364/PRJ.406010
(2021).
Google Scholar
M. N’Gom et al.,
“Controlling light transmission through highly scattering media using semi-definite programming as a phase retrieval computation method,”
Sci. Rep., 7
(1), 2518 https://doi.org/10.1038/s41598-017-02716-x
(2017).
Google Scholar
G. Huang et al.,
“Retrieving the optical transmission matrix of a multimode fiber using the extended Kalman filter,”
Opt. Express, 28
(7), 9487
–9500 https://doi.org/10.1364/OE.389133 OPEXFF 1094-4087
(2020).
Google Scholar
C. A. Metzler et al.,
“Coherent inverse scattering via transmission matrices: efficient phase retrieval algorithms and a public dataset,”
in IEEE Int. Conf. Comput. Photogr. (ICCP),
1
–16
(2017). https://doi.org/10.1109/ICCPHOT.2017.7951483 Google Scholar
S. Cheng, T. Zhong and P. Lai,
“Non-convex optimization for retrieving the complex transmission matrix of a multimode fiber,”
in TENCON 2022-2022 IEEE Region 10 Conf. (TENCON),
1
–5
(2022). https://doi.org/10.1109/TENCON55691.2022.9977923 Google Scholar
Z. Wen et al.,
“Fast volumetric fluorescence imaging with multimode fibers,”
Opt. Lett., 45
(17), 4931
–4934 https://doi.org/10.1364/OL.398177 OPLEDP 0146-9592
(2020).
Google Scholar
M. Plöschner et al.,
“Multimode fibre: light-sheet microscopy at the tip of a needle,”
Sci. Rep., 5
(1), 18050 https://doi.org/10.1038/srep18050
(2015).
Google Scholar
J. Bertolotti and O. Katz,
“Imaging in complex media,”
Nat. Phys., 18
(9), 1008
–1017 https://doi.org/10.1038/s41567-022-01723-8 NPAHAX 1745-2473
(2022).
Google Scholar
J. Carpenter, B. J. Eggleton and J. Schröder,
“Observation of Eisenbud–Wigner–Smith states as principal modes in multimode fibre,”
Nat. Photonics, 9
(11), 751
–757 https://doi.org/10.1038/nphoton.2015.188 NPAHBY 1749-4885
(2015).
Google Scholar
M. W. Matthès et al.,
“Optical complex media as universal reconfigurable linear operators,”
Optica, 6
(4), 465
–472 https://doi.org/10.1364/OPTICA.6.000465
(2019).
Google Scholar
G. Wetzstein et al.,
“Inference in artificial intelligence with deep optics and photonics,”
Nature, 588
(7836), 39
–47 https://doi.org/10.1038/s41586-020-2973-6
(2020).
Google Scholar
S. Leedumrongwatthanakun et al.,
“Programmable linear quantum networks with a multimode fibre,”
Nat. Photonics, 14
(3), 139
–142 https://doi.org/10.1038/s41566-019-0553-9 NPAHBY 1749-4885
(2020).
Google Scholar
N. H. Valencia et al.,
“Unscrambling entanglement through a complex medium,”
Nat. Phys., 16
(11), 1112
–1116 https://doi.org/10.1038/s41567-020-0970-1 NPAHAX 1745-2473
(2020).
Google Scholar
E. Tseng et al.,
“Neural nano-optics for high-quality thin lens imaging,”
Nat. Commun., 12
(1), 6493 https://doi.org/10.1038/s41467-021-26443-0 NCAOBW 2041-1723
(2021).
Google Scholar
H. Yu et al.,
“Ultrahigh-definition dynamic 3D holographic display by active control of volume speckle fields,”
Nat. Photonics, 11
(3), 186
–192 https://doi.org/10.1038/nphoton.2016.272 NPAHBY 1749-4885
(2017).
Google Scholar
W. H. Lee,
“Computer-generated holograms: techniques and applications,”
Prog. Opt., 16 119
–232 https://doi.org/10.1016/S0079-6638(08)70072-6 POPTAN 0079-6638
(1978).
Google Scholar
J. Zhong et al.,
“Nonlinear optimization algorithm for partially coherent phase retrieval and source recovery,”
IEEE Trans. Comput. Imaging, 2
(3), 310
–322 https://doi.org/10.1109/TCI.2016.2571669
(2016).
Google Scholar
J. Nocedal and S. Wright,
“Numerical optimization,”
Springer Sci., 35
(67–68), 7
(1999).
Google Scholar
D. C. Liu and J. Nocedal,
“On the limited memory BFGS method for large scale optimization,”
Math. Program., 45
(1), 503
–528 https://doi.org/10.1007/BF01589116 MHPGA4 1436-4646
(1989).
Google Scholar
M. B. Shemirani et al.,
“Principal modes in graded-index multimode fiber in presence of spatial-and polarization-mode coupling,”
J. Lightwave Technol., 27
(10), 1248
–1261 https://doi.org/10.1109/JLT.2008.2005066 JLTEDG 0733-8724
(2009).
Google Scholar
L.-H. Yeh et al.,
“Experimental robustness of Fourier ptychography phase retrieval algorithms,”
Opt. Express, 23
(26), 33214
–33240 https://doi.org/10.1364/OE.23.033214 OPEXFF 1094-4087
(2015).
Google Scholar
J. R. Fienup,
“Phase retrieval algorithms: a comparison,”
Appl. Opt., 21
(15), 2758
–2769 https://doi.org/10.1364/AO.21.002758 APOPAI 0003-6935
(1982).
Google Scholar
R. W. Gerchberg,
“A practical algorithm for the determination of plane from image and diffraction pictures,”
Optik, 35
(2), 237
–246
(1972).
Google Scholar
L. Waller et al.,
“Phase and amplitude imaging from noisy images by Kalman filtering,”
Opt. Express, 19
(3), 2805
–2815 https://doi.org/10.1364/OE.19.002805 OPEXFF 1094-4087
(2011).
Google Scholar
G. R. Brady and J. R. Fienup,
“Nonlinear optimization algorithm for retrieving the full complex pupil function,”
Opt. Express, 14
(2), 474
–486 https://doi.org/10.1364/OPEX.14.000474 OPEXFF 1094-4087
(2006).
Google Scholar
E. J. Candes, X. Li and M. Soltanolkotabi,
“Phase retrieval via Wirtinger flow: theory and algorithms,”
IEEE Trans. Inf. Theory, 61
(4), 1985
–2007 https://doi.org/10.1109/TIT.2015.2399924 IETTAW 0018-9448
(2015).
Google Scholar
A. Drémeau and F. Krzakala,
“Phase recovery from a Bayesian point of view: the variational approach,”
in IEEE Int. Conf. Acoust. Speech Signal Process. (ICASSP),
3661
–3665
(2015). https://doi.org/10.1109/ICASSP.2015.7178654 Google Scholar
E. J. Candes, T. Strohmer and V. Voroninski,
“Phaselift: exact and stable signal recovery from magnitude measurements via convex programming,”
Commun. Pure Appl. Math., 66
(8), 1241
–1274 CPMAMV 0010-3640
(2013).
Google Scholar
G. Wang et al.,
“Phase retrieval via reweighted amplitude flow,”
IEEE Trans. Signal Process., 66
(11), 2818
–2833 https://doi.org/10.1109/TSP.2018.2818077 ITPRED 1053-587X
(2018).
Google Scholar
T. Goldstein and C. Studer,
“Phasemax: convex phase retrieval via basis pursuit,”
IEEE Trans. Inf. Theory, 64
(4), 2675
–2689 https://doi.org/10.1109/TIT.2018.2800768 IETTAW 0018-9448
(2018).
Google Scholar
L. Gong et al.,
“Optical orbital-angular-momentum-multiplexed data transmission under high scattering,”
Light Sci. Appl., 8
(1), 27 https://doi.org/10.1038/s41377-019-0140-3
(2019).
Google Scholar
J. W. Goodman, Introduction to Fourier Optics, Roberts and Company(
(2005). Google Scholar
BiographyJingshan Zhong is an associate research scientist at Zhejiang Lab, China. He received his BS degree in electronic science and technology from the University of Science and Technology of China, Hefei, China and his PhD in electrical and electronic engineering from Nanyang Technological University, Singapore, in 2010 and 2015, respectively. He got his postdoc training from the Computational Imaging Lab at UC Berkeley, advised by Prof. Laura Waller. His research interests include phase imaging, light-field imaging, imaging through scattering, signal processing, numerical optimization, and machine learning. Zhong Wen obtained his BS degree from Dalian University of Technology in 2016. He is currently pursuing his PhD in optical engineering at Zhejiang University. His research interests include endoscopy imaging and computational imaging. Quanzhi Li received his BS degree in electronic science and technology from Harbin Institute of Technology in 2021. He is currently pursuing his PhD in optical engineering at Zhejiang University. His research interests include multimode fiber imaging and computational imaging. Qilin Deng is a graduate student at Zhejiang University. He completed his undergraduate studies at Jiangnan University, Wuxi, China, earning his BS degree with honors in 2020. His research interests focus on computational imaging, optics, and fiber imaging. Qing Yang received her PhD from the College of Materials Science and Engineering, Zhejiang University, China, in 2006. She worked as a visiting scientist in the Department of Materials Science and Engineering, Georgia Institute of Technology, United States, from 2009 to 2012. She worked as a visiting scientist at the University of Cambridge, United Kingdom, in 2018. Currently, she is a professor at the College of Optical Science and Engineering of Zhejiang University. She is also the director of the Research Center for Humanoid Sensing, Zhejiang Lab, Hangzhou, China. Her research focuses on micro/nanophotonics, nanomaterials, and endoscopy imaging. |