Imaging 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 $1200\times$ 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 $50\mu \mathrm{m}$ 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.

The reference-less experimental setup is shown in Fig. 2. An MMF of 0.22 NA and $50\mu \mathrm{m}$ 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 $-1\mathrm{st}$ diffraction order on the iris. Tuning M2 and M3 shifts the $-1\mathrm{st}$ 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 $10\times$ 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 $20\times$ 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 $525\text{frame}/\mathrm{s}$. By moving the OBJ2 with a stage (Thorlabs CT1P), the CMOS camera captures the intensity images at a defocus plane.

Our 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.

The intensity images are denoted as $I_n(x,y)$, where $n=1,\dots ,N$, $N$ is the total number of the measured intensity images, and $x,y$ are the spatial coordinates. Each image contains $N_x\times N_y$ pixels. The incident complex field is raster-scanned into a vector $e^{j{\mathit{\theta }}_n}$, which has $N_k$ modulation modes for both polarizations. The forward model of the intensity measurement is written as $${\mathbf{I}}_n={|\mathbf{T}{e}^{j{\mathit{\theta }}_n}|}^2,$$(1)where ${\mathbf{I}}_n$ is a vector raster-scanned from $I_n(x,y)$, $\mathbf{T}$ is the transmission matrix ($N_x*N_y\times N_k$), and $|·|$ takes absolute square of the complex numbers inside.

The optimization problem of retrieving the TM from the intensity images is formulated as $$\underset{\mathbf{T}}{\min }\sum _n{\Vert {\mathbf{I}}_n-{|\mathbf{T}{e}^{j{\mathit{\theta }}_n}|}^2\Vert }_2^2,$$(2)where ${|·|}_2^2$ 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 $\mathbf{T}$ by minimizing the cost function.

The large number of unknown in $\mathbf{T}$ makes the optimization problem in Eq. (2) difficult to solve directly. However, it can be broken down into $N_x*N_y$ smaller optimization problems. Each problem is formulated based on the intensity measurement at one single pixel, $$\underset{{\mathbf{t}}^k}{\min }f({\mathbf{t}}^k)={\Vert {\mathbf{I}}^k-{|{\mathbf{Qt}}^k|}^2\Vert }_2^2,$$(3)where $${\mathbf{I}}^k=\left[\begin{array}{c}{\mathbf{I}}_1^k\\ {\mathbf{I}}_2^k\\ ⋮\\ {\mathbf{I}}_N^k\end{array}\right],\mathbf{Q}=\left[\begin{array}{c}{e}^{j{\mathit{\theta }}_{1}T}\\ e^{j{\mathit{\theta }}_{2}T}\\ ⋮\\ e^{j{\mathit{\theta }}_{N}T}\end{array}\right].$$(4)

Here the column vector ${\mathbf{t}}^k$ is the transpose of the $k^{th}$ row of $\mathbf{T}$, ${\mathbf{I}}_n^k$ is the $k^{th}$ element of ${\mathbf{I}}_n$, $T$ denotes transpose, and $k=1,\dots ,N_x*N_y$. The vector ${\mathbf{I}}^k$ contains all of the measurements at the same pixel indexed by $k$ on the intensity images. The matrix $\mathbf{Q}$ is the so-called probing matrix; each row of $\mathbf{Q}$ is one of the incident fields. Each optimization problem recovers one row of $\mathbf{T}$ from the measurements at the corresponding pixel. Thus the whole TM can be recovered by solving these small pixel-wise optimization problems independently.

Conventional reference-less TM retrieval methods19^–21 set the phases of the probing matrix $\mathbf{Q}$ as random numbers. It means that the incident fields are modulated with random phases in these methods. By contrast, our method designs the matrix $\mathbf{Q}$ with Fourier transform matrix, $$\mathbf{Q}=\left[\begin{array}{c}\mathbf{K}\mathrm{diag}(e^{j{\mathit{\psi }}_1})\\ \mathbf{K}\mathrm{diag}(e^{j{\mathit{\psi }}_2})\\ ⋮\\ \mathbf{K}\mathrm{diag}(e^{j{\mathit{\psi }}_M})\end{array}\right],$$(5)where $\mathbf{K}$ is the Fourier transform matrix, $e^{j{\mathit{\psi }}_m}$ is a $N_k$ by 1 column vector with its phase set as random numbers, $\mathrm{diag}(e^{j{\mathit{\psi }}_m})$ is a diagonal matrix with its diagonal entries from $e^{j{\mathit{\psi }}_m}$, and $m=1,\dots ,M$. The incident complex field has $N_{kx}\times N_{ky}$ modulation modes for each polarization, so we have $N_k=2N_{kx}*N_{ky}$. The total number of measured intensity images is $N=M*N_k$, which is $M$-fold of the size of the unknown complex variable, ${\mathbf{t}}^k$. We set $M\ge 3$, since the complex variable in the inverse problem contains both real and imaginary parts. The matrix $\mathbf{Q}$ consists of $M$ square matrix in the form of $\mathbf{K}\mathrm{diag}(e^{j{\mathit{\psi }}_m})$, where $\mathbf{K}$ is the 2D Fourier transform matrix for $N_{kx}\times 2N_{ky}$ matrix. Since $\mathbf{K}$ is a pure phase matrix, the probing matrix $\mathbf{Q}$ 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 $\mathbf{Q}$ are random,19^–21 the multiplication of $\mathbf{Q}$ with a vector has to be computed with matrix–vector multiplication. Our method designs the probing matrix $\mathbf{Q}$ with the Fourier transform matrix, so the matrix–vector multiplication related to $\mathbf{Q}$ can be computed with FFT or inverse FFT (for its complex transpose ${\mathbf{Q}}^H$). 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, $$\begin{gathered} {\mathbf{t}}_0^k={({\mathbf{Q}}^H\mathbf{Q})}^{-1}{\mathbf{Q}}^H\sqrt{{\mathbf{I}}^k} \\ =\frac{1}{M}{\mathbf{Q}}^H\sqrt{{\mathbf{I}}^k}, \end{gathered}$$(6)where $\sqrt{{\mathbf{I}}^k}$ takes the element-wise square root of the vector ${\mathbf{I}}^k$, and $H$ denotes the complex transpose. Since ${\mathbf{K}}^H$ is the inverse Fourier transform matrix, the matrix–vector multiplication in ${\mathbf{Q}}^H\sqrt{{\mathbf{I}}^k}$ can be computed with FFT. We derive the first derivative of $f({\mathbf{t}}^k)$ with respect to ${\mathbf{t}}^k$ as (more details in Appendix A) $${\frac{\partial f}{\partial {\mathbf{t}}^k}}^H=-4{\mathbf{Q}}^H\mathrm{diag}({\mathbf{Qt}}^k)({\mathbf{I}}^k-{|{\mathbf{Qt}}^k|}^2),$$(7)$$\frac{\partial f^H}{\partial {\mathbf{t}}^k}=\frac{{\partial f}^{\mathrm{T}}}{\partial {\mathbf{t}}_x^k}+j\frac{{\partial f}^{\mathrm{T}}}{\partial {\mathbf{t}}_y^k},$$(8)where ${\mathbf{t}}_x^k$ and ${\mathbf{t}}_y^k$ are the real and imaginary parts of ${\mathbf{t}}^k$. The matrix–vector multiplication related to $\mathbf{Q}$ 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 $k$th pixel ${\mathbf{I}}^k$ and random phase vectors for the probing matrix $e^{j{\mathit{\psi }}_m}$. It recovers one row of the TM, ${\mathbf{t}}^k$. The estimation is initialized by Eq. (6). After obtaining the error [Eq. (3)] and gradient [Eq. (7)], the algorithm updates the estimate of ${\mathbf{t}}^k$ 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 $\mathbf{Q}$ in Eqs. (3), (6), and (7) can be efficiently computed with FFT, reducing the computational complexity from $\mathrm{\Theta }(N{N}_k)$ to $\mathrm{\Theta }(N\log N_k)$. It also has the benefit of memory efficiency, since there is no need to explicitly store the big matrix $\mathbf{Q}$ when solving the inverse problem.

One may apply the method in Algorithm 1 on all the $N_x*N_y$ 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 $\mathrm{NA}/\lambda$, where NA is the numerical aperture of the MMF, and $\lambda$ 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, ${\mathrm{ps}}_{\text{intensity}}$, meets the Nyquist sampling theory, $${\mathrm{ps}}_{\text{intensity}}\le \lambda /4\mathrm{NA}.$$(9)

So, the pixel size of the half-sampled images has ${\mathrm{ps}}_{\text{field}}={\mathrm{ps}}_{\text{intensity}}*2$. It meets the sampling requirement of the transmitted complex field, $${\mathrm{ps}}_{\text{field}}\le \lambda /2\mathrm{NA}.$$(10)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 ${\mathbf{t}}^k$ 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.

The 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, $e^{j{\mathit{\varphi }}_0}$ 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 $${\mathbf{T}}_{\text{true}}=\mathrm{diag}(e^{j\mathit{\varphi }}){\mathbf{T}}_{\mathrm{est}},$$(11)where ${\mathbf{T}}_{\text{true}}$ is the true TM, ${\mathbf{T}}_{\mathrm{est}}$ is the estimated TM, and the vector $e^{j\mathit{\varphi }}$ is the phase offset (size $N_x^{\text{half}}*N_y^{\text{half}}\times 1$). Note that $N_x^{\text{half}}$ and $N_y^{\text{half}}$ 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 ${\mathbf{T}}_{\mathrm{est}}$ 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 ${\mathbf{T}}_{\mathrm{est}}$. 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 ${\mathbf{T}}_{\mathrm{est}}$ 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 $z_d$ away from the distal end of the MMF ($z=0$). 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 $I_{n_d}(x,y,z_d)$, where $n_d=1,\dots ,N_d$. Each intensity image has a size of $N_x\times N_y$ with pixel size of ${\mathrm{ps}}_{\text{intensity}}$. From Eq. (11), the transmitted complex field at the distal end of the MMF can be expressed as $${\mathbf{c}}_{n_d}=\mathrm{diag}({\mathbf{T}}_{\mathrm{est}}{e}^{j{\mathit{\theta }}_{n_d}})e^{j\mathit{\varphi }},$$(12)where $e^{j{\mathit{\theta }}_{n_d}}$ is raster-scanned from the corresponding incident complex field of the defocus intensity image, $I_{n_d}(x,y,z_d)$. Note that the transmitted complex field predicted by the estimated TM has a size of $N_x^{\text{half}}\times N_y^{\text{half}}$ with pixel size of ${\mathrm{ps}}_{\text{field}}$. The vector ${\mathbf{c}}_{n_d}$ 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 $${\mathbf{I}}_{n_d}={|{\mathbf{K}}_2^H\mathbf{P}\mathrm{diag}(\mathbf{h}){\mathbf{K}}_1\mathrm{diag}({\mathbf{T}}_{\mathrm{est}}{e}^{j{\mathit{\theta }}_{n_d}})e^{j\mathit{\varphi }}|}^2,$$(13)where the vector ${\mathbf{I}}_{n_d}$ is raster-scanned from ${\mathbf{I}}_{n_d}(x,y,z_d)$, ${\mathbf{K}}_1$ is the Fourier transform matrix for the $N_x^{\text{half}}\times N_y^{\text{half}}$ matrix, $\mathbf{h}$ is the defocus propagation kernel (more in Appendix B), $\mathbf{P}$ is for zero padding in the frequency domain, and ${\mathbf{K}}_2^H$ is the inverse Fourier transform matrix for the $N_x\times N_y$ matrix. The measured intensity has a size of $N_x\times N_y$, whereas the transmitted complex field has a size of $N_x^{\text{half}}\times N_y^{\text{half}}$. 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 $$\underset{\mathit{\varphi }}{\min }g(\mathit{\varphi })=\sum _{n_d}{\Vert {\mathbf{I}}_{n_d}-{|{\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }}|}^2\Vert }_2^2,$$(14)where ${\mathbf{A}}_{n_d}={\mathbf{K}}_2^H\mathbf{P}\mathrm{diag}(\mathbf{h}){\mathbf{K}}_1\mathrm{diag}({\mathbf{T}}_{\mathrm{est}}{e}^{j{\mathit{\theta }}_{n_d}})$. 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 $${\frac{\partial g}{\partial \mathit{\varphi }}}^H=\sum _{n_d}\text{real}[-4\mathrm{diag}(-j{e}^{-j\mathit{\varphi }}){\mathbf{A}}_{n_d}^H\mathrm{diag}({\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }})({\mathbf{I}}_{n_d}-{|{\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }}|}^2)].$$(15)

More details can be found in Appendix C. The matrix–vector multiplication related to ${\mathbf{A}}_{n_d}$ 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.

In 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 $9216\times 8192$ 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 $50\mu \mathrm{m}$ diameter. The incident complex fields had $64\times 64$ phase modulation modes for each polarization. The matrix $\mathbf{K}$ in the probing matrix was set as the Fourier transform matrix for a $64\times 128$ matrix. The transmitted complex fields at the distal end were sampled with $96\times 96\text{pixels}$. We simulated seven datasets with $M=3$ to 9, where $M$ is the total number of $\mathbf{K}\mathrm{diag}(e^{j{\mathit{\psi }}_m})$ 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 $M=8$, it took 376.4 s to retrieve the entire TM of size $9216\times 8192$. 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 $96\times 96$ 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 $M=3$, most rows of the recovered TM have large errors. For $M=4$ 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 $M$ increases. The error of the recovered TM becomes negligibly small for $M=7$ 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 $M=7$ 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 $32\times 32\text{pixels}$ of the $96\times 96\text{pixels}$ of the dataset of $M=8$ 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 $1024\times 8192$ 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.

In the experiment, we used an MMF of 0.22 NA and $50\mu \mathrm{m}$ diameter. The illumination was laser of 488 nm. The DMD achieved $64\times 64$ phase modulation for each polarization, resulting in 8192 modes in total. We tested the TM retrieval algorithm for the cases of $M=4$ 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 $\mathbf{K}$ was set as the Fourier transform matrix for a $64\times 128$ matrix. The phase of the probing matrix was loaded into the DMD, and a series of $M*8192$ images [Fig. 5(a)] were measured. Each image has $128\times 128\text{pixels}$, with a pixel size of $0.47\mu \mathrm{m}$. The preprocessing step half-sampled the measured images and obtained images of $64\times 64$ [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 $M=7$. The probing matrix has a size of $\mathrm{57,344}\times 8192$. The preprocessed intensity images contain 57,344 images of size $64\times 64$. The number of pixels inside the white region of the fiber mask is 2286, so the retrieved TM has size of $2286\times 8192$. 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 $\mathrm{57,344}\times 1$) 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 $M=7$, 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 ($128\times 128\text{pixels}$) for each focus with an exposure time of 70 and $1400\mu \mathrm{s}$ (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 $7\times 7\text{pixels}$ near the peak of the focus using the $70\mu \mathrm{s}$ image, while the total energy is the sum of the signal and the background (outside the $7\times 7\text{pixels}$), which is calculated using the $1400\mu \mathrm{s}$ 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 $M$ and the holography method. For the cases $M=4$ to 6, there are several foci that have a low PR. However, for the cases of $M=\mathrm{7,8}$, the overall quality of the foci is near to that of the holography method. The average PR of the case of $M=8$ 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 $M=8$ and holography by showing the distribution of the PR and a sum projection of several selected foci. The TM retrieval method by $M=8$ 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.

In this section, we validate the TM retrieval algorithm with phase correction. For simulation, we used a simulated TM of size $\mathrm{16,384}\times 9216$ for an MMF with 0.22 NA and $100\mu \mathrm{m}$ diameter, generated by solving Maxwell’s equations.38 The transmitted complex fields of the MMF are sampled by $128\times 128$ grids with a pixel size of $1.1667\mu \mathrm{m}$, and the wavelength of illumination is 532 nm. We designed a probing matrix with 2D Fourier transform matrix for a $96\times 96$ matrix and $M=9$. A series of 82,944 images of size $128\times 128$ were generated at the distal end of the fiber ($z=0\mu \mathrm{m}$). Then we simulated 50 defocus images [Fig. 6(a)] at $z_d=50\mu \mathrm{m}$ away from the distal end of the MMF. Each image has $256\times 256\text{pixels}$ with a pixel size of $0.5833\mu \mathrm{m}$. 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 $6.4\times {10}^{-10}$ and $3.9\times {10}^{-5}$, 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 $50\mu \mathrm{m}$ diameter and 0.22 NA, and illumination wavelength of 488 nm. The phase modulation on DMD had $64\times 64$ modes for each polarization. We designed a probing matrix using Fourier transform matrix for a $64\times 128$ matrix, and $M=8$. 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 ($z=0\mu \mathrm{m}$). Each image has $192\times 192$ pixels with a pixel size of $0.4182\mu \mathrm{m}$. In order to correct the phase offset, we measured 50 images of $192\times 192$ at $40\mu \mathrm{m}$ 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 $z=0\mu \mathrm{m}$ by the proposed method in Fig. 3. In the preprocessing step, we half-sampled the images to a size of $96\times 96$ 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 $z=0\mu \mathrm{m}$ with a free-space defocus propagation. We generated the two sets of propagated TMs at $z=0,-20,-40,-60,-80$, and $-100\mu \mathrm{m}$, 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 $z=0\mu \mathrm{m}$, 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 $100\mu \mathrm{m}$, 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.

The 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 $Q$ ($N\times N_k$) with the Fourier transform matrix. Using FFT, the computational complexity of the matrix–vector multiplication related to $\mathbf{Q}$ and ${\mathbf{Q}}^H$ reduces from $\mathrm{\Theta }(N{N}_k)$ to $\mathrm{\Theta }(N\log N_k)$. Here we give an example of the number of modulation modes $N_k=8192$ and the number of measurements $N=\mathrm{65,536}$. The matrix $\mathbf{Q}$ has a size of $\mathrm{65,536}\times 8192$. The computational complexity reduces from $\mathrm{\Theta }(\mathrm{65,536}\times 8192)$ to $\mathrm{\Theta }(\mathrm{65,536}\times 13)$, and it is memory-efficient without storing $\mathbf{Q}$. The computation related to the probing matrix $\mathbf{Q}$ 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 $\mathrm{\Theta }(N{N}_k)$, due the matrix–vector multiplication related to $\mathbf{Q}$. It is higher than that of our proposed method using FFT $\mathrm{\Theta }(N\log N_k)$. However, applying the similar FFT-based scheme in these algorithms could further reduce the computational complexity.

We 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 [$\mathrm{\Theta }(N\log N_k)$] versus [$\mathrm{\Theta }(N{N}_k)$]. 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 $1200\times$ speedup in recovering the TM of size $2286\times 8192$ with 112.9 s for the MMF of 0.22 NA and $50\mu \mathrm{m}$ 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

The optimization problem to recover one row of the TM is expressed as $$\underset{{\mathbf{t}}^k}{\min }f({\mathbf{t}}^k)={\Vert {\mathbf{I}}^k-{|{\mathbf{Qt}}^k|}^2\Vert }_2^2.$$(16)

Next, we define $$\mathbf{F}={\mathbf{I}}^k-{|{\mathbf{Qt}}^k|}^2,$$(17)$$f({\mathbf{t}}^k)={\mathbf{F}}^H\mathbf{F},$$(18)where $\mathbf{F}$ is a vector. According to the chain rule, the first derivative of the cost function can be written as $$\begin{gathered} \frac{\partial f}{\partial {\mathbf{t}}^k}=-\frac{\partial f}{\partial \mathbf{F}}\frac{\partial \mathbf{F}}{\partial {\mathbf{t}}^k} \\ =-\frac{\partial f}{\partial \mathbf{F}}\frac{\partial {|{\mathbf{Qt}}^k|}^2}{\partial {\mathbf{t}}^k} \\ =-4{\mathbf{F}}^H\mathrm{diag}(\mathrm{conj}({\mathbf{Qt}}^k))\mathbf{Q}. \end{gathered}$$(19)

Thus we have the Hermitian of the first derivative as $$\begin{gathered} {\frac{\partial f}{\partial {\mathbf{t}}^k}}^H=-4{\mathbf{Q}}^H\mathrm{diag}({\mathbf{Qt}}^k)\mathbf{F} \\ =-4{\mathbf{Q}}^H\mathrm{diag}({\mathbf{Qt}}^k)({\mathbf{I}}^k-{|{\mathbf{Qt}}^k|}^2). \end{gathered}$$(20)

According to the theory of angular spectrum propagation,50 the defocus propagation kernel in frequency domain is expressed as $$h(u,v,z_d)=\exp \left[j\frac{2\pi }{\lambda }\sqrt{1-{(\lambda u)}^2-{(\lambda v)}^2}{z}_d\right]p(u,v),$$(21)where $\lambda$ is the wavelength of the illumination, $u$ and $v$ are the spatial frequency coordinates, and $p(u,v)$ is the pupil of the imaging system. The pupil is written as $$P(u,v)=\{\begin{array}{cc}1,& \lambda \sqrt{u^2+v^2}\le \mathrm{NA},\\ 0,& \lambda \sqrt{u^2+v^2}>\mathrm{NA}.\end{array}$$(22)

The vector $\mathbf{h}$ is raster-scanned from $h(u,v,z_d)$.

The optimization of solving the phase offset from the defocus intensity images is rewritten as $$\underset{\mathit{\varphi }}{\min }g(\mathit{\varphi })=\sum _{n_d}{\Vert {\mathbf{I}}_{n_d}-{|{\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }}|}^2\Vert }_2^2,$$(23)where ${\mathbf{A}}_{n_d}={\mathbf{K}}_2^H\mathbf{P}\mathrm{diag}(\mathbf{h}){\mathbf{K}}_1\mathrm{diag}({\mathbf{T}}_{\mathrm{est}}{e}^{j{\mathit{\theta }}_{n_d}})$. Next, we define $${\mathbf{G}}_{n_d}={\mathbf{I}}_{n_d}-{|{\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }}|}^2,$$(24)$$g_{n_d}={\mathbf{G}}_{n_d}^H{\mathbf{G}}_{n_d}.$$(25)

Using the chain rule, we have $$\begin{gathered} \frac{\partial g_{n_d}}{\partial e^{j\mathit{\varphi }}}=-\frac{\partial g_{n_d}}{\partial {\mathbf{G}}_{n_d}}\frac{\partial {\mathbf{G}}_{n_d}}{\partial e^{j\mathit{\varphi }}} \\ =-\frac{\partial g_{n_d}}{\partial {\mathbf{G}}_{n_d}}\frac{\partial {|{\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }}|}^2}{\partial e^{j\mathit{\varphi }}} \\ =-4{\mathbf{G}}_{n_d}^H\mathrm{diag}[\mathrm{conj}({\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }})]{\mathbf{A}}_{n_d}. \end{gathered}$$(26)

We can have $$\frac{\partial e^{j\mathit{\varphi }}}{\partial \mathit{\varphi }}=\mathrm{diag}(j{e}^{j\mathit{\varphi }}).$$(27)

By combining Eqs. (26) and (27), we can get $$\frac{\partial g_{n_d}}{\partial \mathit{\varphi }}=\text{real}\{-4{\mathbf{G}}_{n_d}^H\mathrm{diag}[\mathrm{conj}({\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }})]{\mathbf{A}}_{n_d}\mathrm{diag}(j{e}^{j\mathit{\varphi }})\}.$$(28)

It is easy to obtain $$\frac{{\partial g_{n_d}}^H}{\partial \mathit{\varphi }}=\text{real}[-4\mathrm{diag}(-j{e}^{-j\mathit{\varphi }}){\mathbf{A}}_{n_d}^H\mathrm{diag}({\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }})({\mathbf{I}}_{n_d}-{|{\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }}|}^2)].$$(29)

So, we have $$\begin{gathered} \frac{\partial g^H}{\partial \mathit{\varphi }}=\sum _{n_d}\frac{{\partial g_{n_d}}^H}{\partial \mathit{\varphi }} \\ =\sum _{n_d}\text{real}[-4\mathrm{diag}(-j{e}^{-j\mathit{\varphi }}){\mathbf{A}}_{n_d}^H\mathrm{diag}({\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }})({\mathbf{I}}_{n_d}-{|{\mathbf{A}}_{n_d}{e}^{j\mathit{\varphi }}|}^2)]. \end{gathered}$$(30)

Jingshan 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.