## Abstract

Optical imaging through or inside scattering media, such as multimode fiber and biological tissues, has a significant impact in biomedicine yet is considered challenging due to the strong scattering nature of light. In the past decade, promising progress has been made in the field, largely benefiting from the invention of iterative optical wavefront shaping, with which deep-tissue high-resolution optical focusing and hence imaging becomes possible. Most of the reported iterative algorithms can overcome small perturbations on the noise level but fail to effectively adapt beyond the noise level, e.g., sudden strong perturbations. Reoptimizations are usually needed for significant decorrelation to the medium since these algorithms heavily rely on the optimization performance in the previous iterations. Such ineffectiveness is probably due to the absence of a metric that can gauge the deviation of the instant wavefront from the optimum compensation based on the concurrently measured optical focusing. In this study, a square rule of binary-amplitude modulation, directly relating the measured focusing performance with the error in the optimized wavefront, is theoretically proved and experimentally validated. With this simple rule, it is feasible to quantify how many pixels on the spatial light modulator incorrectly modulate the wavefront for the instant status of the medium or the whole system. As an example of application, we propose a novel algorithm, the dynamic mutation algorithm, which has high adaptability against perturbations by probing how far the optimization has gone toward the theoretically optimal performance. The diminished focus of scattered light can be effectively recovered when perturbations to the medium cause a significant drop in the focusing performance, which no existing algorithms can achieve due to their inherent strong dependence on previous optimizations. With further improvement, the square rule and the new algorithm may boost or inspire many applications, such as high-resolution optical imaging and stimulation, in instable or dynamic scattering environments.

© 2021 Chinese Laser Press

## 1. INTRODUCTION

Optical focusing through or within scattering media has been long desired yet has been considered challenging until the invention of optical wavefront shaping [1,2]. As reported so far, wavefront shaping techniques can be categorized into optical phase conjugation (OPC) [3–6] and iterative wavefront shaping (including the transmission matrix measurement method) [7–9]. OPC can achieve rapid optical focusing ($\sim 1\text{\hspace{0.17em}}\mathrm{ms}$), but diffusive light can only be focused back to the original light source, and the system requires accurate alignment to fulfill rigorous phase conjugation. In comparison, iterative wavefront shaping can function more freely, even though more iterations and measurements are needed. With a much simpler setup to control the light in a larger field of view, iterative wavefront shaping allows arbitrary image transmission [10] and raster scanning of the focal point [11] through/inside scattering media. In this technology, photons are modulated by a spatial light modulator (SLM) before they are multiply scattered and become diffusive in the medium. A series of modulation patterns is displayed on the SLM under the control of an iterative optimization algorithm, and an optimal phase pattern is obtained when the feedback signal is maximized, indicating that the scattering-induced phase distortions are compensated as much as possible [12–14]. This technology has broad applications in biomedicine, for instance, in deep tissue imaging [6,15], phototherapy [16], laser surgery, and photoacoustic imaging [8]. Moreover, wavefront shaping has been used to control the light transmission through a multimode fiber (MMF), in which speckles arise due to the intermodal interference and mode dispersion [17,18]. MMF-based biomedical applications, such as endoscopic imaging [19,20] and optogenetics [21,22], are therefore progressively developed. Demonstrations based on MMF for wavefront shaping are therefore of interest since all these applications might be affected by the perturbations to the MMF on different levels, including consistently environmental noises (e.g., temperature change, pressure) and sudden strong perturbations (e.g., mechanical disturbances or biological motions). Those perturbations are due to the regular and/or irregular motion of the scattering medium or the system. Such instability is also a major factor that impedes wavefront-shaping-assisted MMF from wider applications, since the modulated optical field through the MMF will accordingly decorrelate. The decorrelation essentially indicates that the transmission matrix (TM) of the MMF is highly susceptible to perturbations, such as bending, twisting, or temperature change [23,24].

To combat against the instability, a carefully designed setup can be utilized to separate the target photon from the perturbations due to the instable medium [25]. On the other hand, an efficient optimization algorithm for wavefront shaping is also required. Various algorithms have been reported and discussed, such as the continuous sequential algorithm (CSA) [26,27], the stepwise sequential algorithm (SSA) [26], the partitioning algorithm (PA) [26], the genetic algorithm (GA) [9,28–30], particle swarm optimization (PSO) [31–33], the simulated annealing algorithm (SA) [34,35], and some artificial-intelligence-assisted algorithms [36–38]. These methods have realized superior light focusing and even noise-resistance ability. But their optimization mechanisms originate from numerical optimization without specific consideration about the physics behind the strong scattering process. These methods heavily rely on (or learn from) the net performance accumulated from previous iterations, especially the GA, which needs a large pool containing many random phase/amplitude masks to initiate optimization. The “learning experience” is highly specific to the current status of the medium or the system, and hence the optimization can merely adapt and generalize to the subtle instability of the medium, e.g., on the environmental noise level [28]. Once the perturbations to the medium are further strengthened beyond the noise level, e.g., a sudden perturbation, the transmission matrix of the medium may be altered significantly, and the resultant optical focus probably fades or even disappears. In this regard, another optimization is inevitable since the modulation patterns optimized in previous steps (without perturbations) are now weakly correlated with the new optimization condition that matches the state of the perturbed medium. An adaptive optimization algorithm is therefore desired, aiming to avoid strong dependency on the optimization from previous iterations.

Probably the fact that existing algorithms are less adaptive to perturbations is because of the absence of a practical metric that can directly relate the instant focusing status with the accuracy of the optimized wavefront. Such possibility is theoretically explored in this study based on the fully developed speckle patterns, which can be easily observed within or behind strong scattering media, such as an MMF or biological tissue. The intensity of a fully developed speckle pattern is governed by the negative exponential decay. As a result, phasors or elements in the corresponding TM of the medium follow a circular Gaussian distribution as discussed by Goodman [39]. Based on this plain assumption, we derive and define a metric called error rate (denoted as $r$), specific for the binary-amplitude modulation, to estimate how many pixels on the SLM are wrongly set based on the concurrently measured optical focusing: $r$ is physically related to the focusing performance as measured by peak-to-background ratio (PBR) by a simple square rule. This metric can imply how far the optimization has gone towards the theoretically optimal phase compensation, namely, the ideal single-point focusing. Therefore, based on the real-time probed error rate instead of parameters inherited or accumulated from previous iterations, a novel algorithm called the dynamic mutation algorithm (DMA) is developed in this study as an application of the proposed practical metric. The optimization based on the error rate can automatically adapt strong perturbations: compared with other algorithms, the proposed DMA is advantageous, in both simulations and experiments, by its high adaptability and unique recovery ability against dynamic changes. Also, the diminished focus under perturbations (for example, twisting an MMF) can be effectively regained without additional operations, such as repeating the iterative optimization. With such a guiding metric, the new algorithm may inspire further optimization of optical focusing in dynamic media from a practical perspective and boost more applications of wavefront shaping in living biological tissue.

## 2. PRINCIPLE

The DMA is an optimization method based on real-time experimental data instead of results from previous iterations, leading to high adaptability to dynamic changes. The key of the algorithm is to estimate the error rate of the corresponding amplitude mask, which implies the extent that the measured results deviate from the theoretical one, and then guide the optimization towards the optimal solution. In this section, we will elucidate the concept of the error rate and the DMA, followed by how the error rate can be used to modify the amplitude masks to achieve adaptive focusing. The steps involved in the optimization will also be explained.

#### A. Error Rate and Square Rule

Beginning from the ideal focusing, the optimized wavefront, modulated by a binary-amplitude mask, is unique due to the deterministic transmission matrix of a medium ($T$), with elements *t _{mn}*. For simplicity, the input optical field is considered a plane wave (with all phases set to zero) and modulated by a binary-amplitude-only SLM, i.e., digital micromirror devices (DMDs). Denoting the optimal wavefront ${E}_{\mathrm{op},1}={[{e}_{1},\cdots ,{e}_{N}]}^{T}$ as the optimal modulation for optical focusing at the first output channel, the optical field (${U}_{f,1}$) with the first output channel focused is governed by

*p*; the rest of the $(1-r)N$ input channels are reordered with index $q$. Considering the binary-amplitude modulation, the optimal amplitude of the elements in ${E}_{\mathrm{op},1}$ is determined by the first row in $T$ and element (${e}_{i}$) and only turned “on” if the real part of ${t}_{1i}$, denoted as ${R}_{i}$, is greater than zero:

*I*) of $\mathit{T}$ are statistically independent and follow the probability distribution density $f(x)=\mathrm{exp}(-{x}^{2}/2{\sigma}^{2})/\sqrt{2\pi {\sigma}^{2}}$, where $\sigma $ is the standard deviation of the distribution. With Eq. (3), only elements in $T$ with positive real parts ($R>0$) will be selected in the calculation. Since the imaginary part is independent of the real part, the selected elements with $R>0$ have the imaginary part ($I$) fulfilling $-\infty <I<+\infty $. The terms $\u27e8{t}_{1p}{e}_{p}\u27e9$ and $\u27e8{t}_{1q}{e}_{q}\u27e9$ from Eq. (2) can thus be expressed as

*r*-ratio incorrect modulation, can be obtained by

*r*-ratio of pixels oppositely modulates the wavefront in any given modulating masks rather than only in the optimal mask. This generalization will be proved experimentally in Section 2.E.

Mathematically, the relative PBR is essentially related to Pearson’s correlation coefficient (${\rho}_{r}$) between the optimal mask (${E}_{\mathrm{op},1}$) and a mask with $\mathit{r}$-ratio incorrect modulation (${E}_{r,1}$). Their elements, i.e., ${e}_{i}$ and ${e}_{r,i}$ for ${E}_{\mathrm{o}\mathrm{p},1}$ and ${E}_{r,1}$, respectively, obey the symmetric Bernoulli distribution, i.e., $\mathit{e}\sim $Bernoulli ($p=0.5$), so that both $\u27e8{e}_{i}\u27e9$ and $\u27e8{e}_{r,i}\u27e9$ are 0.5 and $\u27e8{e}_{i}{e}_{r,i}\u27e9=0.5(1-r)$. By defining $\delta f=f-\u27e8f\u27e9$, ${\rho}_{r}$ can be formulated as

Experimentally, the estimated error rate ($r$) of any optimized pattern is obtained as follows: an $N$-element DMD is used to modulate the input wavefront, and then a “theoretical PBR” is calculated by ${\eta}_{0}=N/2\pi $; the “experimental PBR” (${\eta}_{\mathrm{ex}}$) is obtained from the instant speckle pattern, i.e., a focus at the target position; finally, the experimental error rate for each focusing optimization can be obtained by substituting Eq. (10) into Eq. (11):

#### B. Mutation Rate

Benefitting from Eq. (9), the percentage of DMD elements with incorrect modulation in the mask can be directly estimated. An ideal solution, or mask, can be obtained if the incorrect elements are corrected. That said, the exact positions/indices of these elements are unknown. Inspired by the mutation process in the GA, having a suitable mutation rate to change the state of modulating elements regarding the error rate may be able to improve the optimization.

In the GA, the mutation rate is usually preset in a decaying manner, and it gradually becomes smaller regardless of the actual optimization performance or status. In the proposed DMA, the mutation rate is adjusted dynamically according to the error rate so that the information about how the instable medium instantly affects the optimized mask can be considered. One more benefit of integrating the error rate is that, in every iteration, the number of elements to be mutated ($N\mu $) can be well controlled below the number of the wrongly modulating elements ($Nr$), i.e., $N\mu \le Nr$. As an example to scale the mutation rate, the mutation rate [$\mu (s)$] in the $s$th iteration can be simply set to be proportional to the instant error rate [$r(s)$] with a mutation constant ($C$) that is greater than unity [Eq. (13)]. To generate new DMD patterns in the $s$th iteration, a total number of $N\mu (s)$ elements in the DMD mask generated in the $(s-1)$th iteration are randomly selected and mutated by reversing the element state of on/off or equivalently following Eq. (14):

The function with respect to the mutation constant, Eq. (13), is to bound the mutation rate between $0.5/C$ and 0, if the error rates at the beginning and at the end of the optimization are assumed to be 0.5 and 0, respectively. The mutation constant is the only parameter needed to be set before optimization; a smaller mutation constant is suggested in instable environments to provide a larger range of mutation rates. The mutation rate is autotuned by the algorithm within the range in response to the actual situations, so as to increase the chance for the incorrect elements to be mutated and lead the optimization towards the theoretical result.#### C. DMA Workflow

The block diagram in Fig. 1 shows the typical workflow of the DMA for the binary-amplitude modulation with DMD. First, all DMD pixels are set to be 1 (“on” state). The error rate is found according to Eq. (12), and the mutation rate is computed through Eq. (13). Then the mask is mutated to generate a new DMD pattern with Eq. (14). Then the error rate is assessed again based on the instantly measured focus performance, i.e., the PBR. If the error rate becomes smaller, which means the PBR is improved, a new mutation rate is calculated according to the error rate. If there is an abrupt rise of the error rate, probably caused by the changes of medium of interest, the mutation rate will be updated as well. Otherwise, if the error rate is just slightly fluctuating and shows no improvement in the PBR, the current mask is mutated again with the same mutation rate. The mutation rate is not updated in this case in order to minimize the error rate fluctuations. This process is repeated until the PBR of optical focusing saturates or plateaus.

#### D. Experimental Setup

The experimental setup used in this study is shown in Fig. 2(a). A continuous-wave 532 nm laser source (EXLSR-532-300-CDRH, Spectra Physics, USA) is used to illuminate the DMD (DLP4100, Texas Instruments Inc., USA). A pair of convex lenses (L1 and L2) is used to expand the light beam, such that it covers all pixels on the DMD. Another pair of convex lenses (L3 and L4) is used to demagnify the beam after it is modulated by the DMD. After that, the shrunk modulated light is focused with a $40\times $ objective lens (NA = 0.65) onto a scattering sample. A CMOS camera (Blackfly S BFS-U3-04S2M-CS, FLIR, Canada) is placed behind the scattering sample, and an image is captured in each measurement for the calculation of the instant PBR and the error rate. The PBR can be calculated by dividing the intensity of the target mode by the average background intensity. A 1 m bare optical MMF (SUH200, Xinrui, China, with diameter $=$ 200 μm, NA $=$ 0.22) is chosen as the scattering sample, with two collimators (PAF2-A4A, Thorlabs, USA) and a fiber rotator (HFR007, Thorlabs, USA). During the optimization for optical focusing, $64\times 64$ input modes are used ($16\times 25$ pixels on the DMD are grouped as a mode, and each pixel is $10.8\text{\hspace{0.17em}}\mathrm{\mu m}\times 10.8\text{\hspace{0.17em}}\mathrm{\mu m}$), and every algorithm is run for 10,000 measurements without stop, even if rotation to the fiber is applied.

In the following sections, the parameters of the investigated algorithms used for the simulation and experiment are set to be the same as follows. For the DMA, the mutation constant in Eq. (13) is set to 200. For the GA, the population size is 20 and the offspring size is 10. The initial mutation rate is 0.1, which decays exponentially to a final value at 0.001. For the CSA, there is no preset parameter, whose input modes are optimized one by one with a linear rastering manner [26,27]. These initial parameters related to specific algorithms are summarized in Table 1. All these algorithms are implemented for 10,000 measurements, and each measurement is set to spend 0.2 s.

As an example, speckle before and after DMA optimization is shown in Fig. 2. Figure 2(b) shows the speckle field before optimization, where the PBR at the target position (central point) is around 3. After wavefront shaping optimization guided by the DMA, the focus has a PBR enhanced to 120 as shown in Fig. 2(c). The full width at half-maximum (FWHM) of the focus is 15.2 and 14.5 μm in the horizontal and vertical directions, respectively.

#### E. Verification of the Square Rule

Numerical and experimental proofs of the square rule have been done to validate its practical implications. The square rule can be simple and straightforward [denoted as “theoretical prediction” in Figs. 3(a) and 3(b)], and it can be computationally re-created in simulation by using a TM whose elements follow the circular Gaussian distribution [denoted as “ideal simulation” in Fig. 3(b)]. Yet, to prove the square rule experimentally, a series of error rates is selected, i.e., 0%, 10%, 20%, 30%, 40%, and 50%, and each data point is the average of five executions. First, we use a TM-based method [42] to generate an optimal DMD mask (set as $r=0$) for an optimized focus through an MMF, and then we mutate the optimal mask with *r*-ratio of pixels to test the corresponding focusing performance, i.e., PBR. As shown in Fig. 3(a), the experimental ${\eta}^{\prime}\text{\u2212}r$ curve is shaped like a parabola with a right shift away from the theoretical curve, ${\eta}^{\prime}={(1-2r)}^{2}$. Such a shift may be attributed to the instability of the MMF, since the instability effect can be accumulated during the process of TM measurement. The measured TM may be subjected to deviation from the ideal assumption that the elements in the TM of the scattering sample follow an ideal circular Gaussian distribution. In view of this, the real part of the measured transmission matrix of the MMF is found to have a right-shifted Gaussian distribution (mean $=$ 0.002 and standard deviation $=$ 0.04) [inset in Fig. 3(a)]. Based on this TM measurement, a TM whose real part follows the Gaussian distribution with mean $=$ 0.002 and standard deviation $=$ 0.04 is generated to form an optimized focus with a series of error rates to degrade the performance. The induced ${\eta}^{\prime}\text{\u2212}r$ curve [denoted as “simulation” in Fig. 3(a)]matches well with the experimental one. Notably, in our experimental setup, the TM measurement needs around 30 min to complete, and the obtained real part of TM is used to generate the optimal mask. Decorrelation of the medium is hardly ignored and inevitably coupled into the measured TM. In addition, all the mutated masks (with $r=10\%\text{\u2212}50\%$) are generated from the same optimal mask after the TM measurement is completed. These masks therefore carry the information of medium instability from the TM measurement. Furthermore, these masks are sequentially displayed on the DMD to modulate the wavefront, which also costs time. When the case of $r=50\%$ is tested, the medium has been altered and decorrelated from the status when the case of $r=0$ is tested. This may imply that the right shift of the experimental ${\eta}^{\prime}\text{\u2212}r$ curve, or the unsatisfactory optimization result, can be attributed to the instable medium represented by a biased or shifted TM.

Nevertheless, the instable medium (on the noise level) can still be governed by the square rule to some extent, but it fails to maintain accuracy. That is because the square rule is based on an unchanged and stable medium. Therefore, to experimentally recreate the square rule, at least, the time span between the generation of the optimal mask and mutated mask is limited within the decorrelation window of the medium. For example, an optimal mask is generated for every error rate investigation during the experiment. By doing so, the instable effect due to the practical medium can be almost eliminated, as shown in Fig. 3(b), where the ${\eta}^{\prime}\text{\u2212}r$ curve attained from experiments matches well with the theoretical curve as well as the ideal simulation curve. Therefore, the square rule functions well in a real-time representation, and in other words, the error rate can provide an effective instant metric to evaluate the distance to the ideal optimal optimization. That will provide a plain yet universal perspective to analyze the imperfect focusing performance for the scattering medium, even if it is heavily perturbed.

## 3. RESULTS

#### A. Simulation

Simulations are done to evaluate the performance of the DMA, which is compared with two representative existing algorithms, i.e., the GA and the CSA. In addition to their popularity, the GA and CSA are selected because they share some similarities with the DMA. Both the DMA and GA have a mutation process and target optimization in instable environments. Meanwhile, the DMA and CSA are straightforward algorithms that do not rely much on previous results. Simulations with the DMA, GA, and CSA have been performed under various conditions of different levels of noise, and the results are compared based on the PBR throughout a fixed number of measurements (or iterations). Whenever the intensity of the target mode is measured, it is counted as one measurement, and the GA usually needs several measurements for one iteration. Each curve in the plots is an averaged result of 50 executions with a new transmission matrix generated to simulate the scattering process for every execution. $N=64\times 64$ input modes (modulating elements for binary-amplitude modulation) are used, and the output mode at the center is chosen for optimization.

### 1. Influence of Noise Level

In this section, the algorithms are compared under different levels of noise. Additive white Gaussian noise is added in every intensity measurement to mimic the instability of the optical system in the actual environment [28]. The Gaussian noises, with standard deviations of 30% and 60% of the initial average intensity $\u27e8{I}_{0}\u27e9$, are set to represent the low-noise and high-noise situations, respectively.

Figures 4(a)–4(c) show the simulation results under the noise-free, low-noise ($0.3\u27e8{I}_{0}\u27e9$), and high-noise ($0.6\u27e8{I}_{0}\u27e9$) conditions, respectively. In the noise-free and low-noise situations, the DMA has the fastest growth of the initial PBR and achieves a high PBR. Although the CSA obtains the highest final PBR in the noise-free case, its performance declines drastically with increased noise level. In the high-noise situation, the GA, benefitting from its large population of optimizing masks, exhibits its superior noise-resisting ability. Meanwhile, the DMA can also reach a comparable level of PBR without such a large population as presented in the GA. Wavefront optimization guided by the error rate is therefore immune to the need for a large population.

### 2. Influence of Transmission Matrix Change

Apart from the noise caused by the instability of the optical system, optimization results can be greatly affected by the instability or slight movement of the scattering medium. To simulate this situation, a 25% right shift of the scattering medium was implemented at the 5000th measurement. The scattering medium was represented by an $\mathit{M}\text{\hspace{0.17em}}\text{output}\text{\hspace{0.17em}}\text{modes}\times \mathit{N}\text{\hspace{0.17em}}\text{input}\text{\hspace{0.17em}}\text{modes}$ transmission matrix. $M\times 0.25N$ new elements, following the same circular Gaussian distribution, are generated and inserted to the left of the matrix. The right 25% of the original matrix elements are removed so that a new $M\times N$ transmission matrix to mimic the shift of the medium is formed. Different from the noise addition process, which only affects the intensity measurement, the shift also leads to changes in the transmission matrix. Simulations are done in noise-free, low-noise, and high-noise situations. The simulation parameters used for the algorithms here are the same as those mentioned in Section 3.A.1. Figures 4(d)–4(f) show how the DMA, GA, and CSA respond to the shift in the transmission matrix.

As seen, right after the transmission matrix is changed, there is a sudden drop of PBR in all three algorithms. The GA fails to adapt to the sudden change of the TM shift for all three investigated noise conditions, which is associated with the mechanism of the GA: the offspring (amplitude masks) with lower cost (intensity) in the population is replaced, and the best offspring is always kept [28]. Also, the whole population is generated based on the medium status before the TM shift, whose dependency and correlation regarding the largely shifted TM are relatively low. Therefore, when it encounters a relatively large enhancement drop, it is hard to produce offspring with cost (i.e., the PBR) larger than the former best one, which makes the optimization trapped in a local maximum. In contrast, the DMA and CSA successfully recover the focus after the TM shift under the noise-free and low-noise conditions. Such achievement is probably due to their absence of a pool with a large population, whose information is strongly related to the status of the medium before the TM shift. Or equivalently, the population size of the DMA and CSA is 1, so the modulating mask can be instantly guided by the information from the sudden shift without constraints from the other masks in the pool. Under the high-noise condition, the DMA is the only algorithm that can adapt to the sudden change and bounce back to the original level after $\sim 5000$ measurements. The CSA, limited by its weak noise-resisting ability, fails to tackle the high-noise conditions.

As a comparison, the square rule, providing the DMD error rate from the instant PBR, shows its advances to deal with different levels of noise conditions and sudden changes. In the next section, experimental performance will be further discussed.

#### B. Experiment

### 1. Focusing Against Strong Noise

With external perturbations on the noise level, the optimizations for single-point focusing via the DMA, GA, and CSA are shown in Fig. 5. The final PBRs that the algorithms achieved can be similarly divided into two groups, i.e., DMA and GA with effective focusing, and CSA with weak effectiveness, compared to the simulation results in the situation with noise of $0.6\u27e8{\mathit{I}}_{0}\u27e9$ [Fig. 4(c)]. Both the DMA and GA demonstrate their robustness in a noisy environment. The DMA shows a higher initial PBR rate and reaches its optimal state after around 5000 measurements. The GA transcends the DMA at around the 7500th measurement. The CSA has a slow initial PBR rate, as the optimization starts from the pixels at the edge of the mask, which contributes less to the optimization due to the Gaussian beam used in the experiment. The contribution from the modulating element increases when it comes to the central part of the mask, and then it slows down again and eventually reaches its maximum when the process approaches another end of the mask. The CSA is sensitive to strong noise [33], so it cannot obtain a PBR as high as the other two algorithms.

### 2. Focusing Against Strong Perturbations

With more apparent perturbations, e.g., a slight movement, a bending, or a small rotation to the MMF, the corresponding TM can be significantly changed and the speckle patterns decorrelated. If that occurs during the experiment, the optimization process is disrupted, and the resultant focal spot may be ruined. In this section, experiments were done to study how perturbations, with rotation to the MMF as an example, affect the optimization, and how different algorithms respond to such heavy instability. The same experimental setup was used with an additional fiber rotator, which can rotate the MMF with various angles.

The relationship between the degree of rotation to the fiber and the corresponding PBR drop is shown in Fig. 6(a). Different degrees of rotation (2.5°, 5°, 7.5°, 10°, 12.5°, and 15°) were introduced when the PBR reached 100. Notably, the optimization algorithms were kept running during the whole optimization and not commanded to stop before and/or after the rotations. As the MMF is altered by the rotation, the PBR drops immediately as seen in Fig. 6(c). Moreover, the more the fiber is rotated, the larger the percentage drop for the PBR. This indicates that the corresponding TM is altered significantly due to the fiber rotation. By applying the DMA to optimize the focus, the optimization can automatically adapt to the degraded focus without re-running the optimization process. As shown in Fig. 6(c), the DMA adapts to the focus degradation with various degrees of fiber twisting, i.e., 2.5°, 5°, and 7.5° rotation corresponding to $\sim 20\%$, $\sim 40\%$, and $\sim 60\%$ PBR drop. The DMA is always able to recover the PBR to the value before perturbation regardless of how much the PBR has dropped. Figure 6(b) shows that the number of measurements required for the PBR to rebound to the level before perturbations increases with respect to the rotation angle. It again validates that strong perturbations can significantly change the status of the medium, which poses challenges to any iterative algorithm. The DMA, on the other hand, shows critical advantages over current popular algorithms.

As an example to study how the DMA, GA, and CSA combat against the heavy instability of the MMF, a 5° rotation for the MMF was implemented at the 5000th measurement for these three algorithms. Figure 7 shows how the algorithms perform throughout the experiments, and Fig. 7 shows the focal spots before optimization (zeroth measurement), right before fiber rotation (5000th measurement), right after fiber rotation (5001st measurement), and after reoptimization (15,000th measurement). The experimental results agree well with the simulation results in Section 3.A: the DMA and CSA show their recovery abilities. The DMA rebounds soon after the rotation and takes around 4000 measurements to regain the PBR it achieved before the rotation, and the resultant optical focus is as bright as, if not brighter than, the focus before perturbation (Fig. 8). Comparably, the CSA does not recover right away after the rotation. It starts to recover after around 6250 measurements.

The recovery efficiency of the CSA may depend on when the perturbation occurs. In the experiment, as the perturbation is induced when the optimizing elements are near the edges of the DMD, the recovery speed is slow. More importantly, merely changing one element for each measurement in the CSA is not efficient to overcome the instability since the positive contribution from one element is probably below the noise level. In contrast, the GA cannot obtain further PBR after the rotation, as the optimization is trapped in the local maximum due to three reasons. First, the optimizing masks in the population library are merely based on the medium status before perturbations. Second, the decorrelation due to 5° rotation cannot be tolerated or generalized from that of the population generated via GA. Third, the mutation process in the GA is not adaptive to the sudden perturbations during optimization since the mutation rate in the GA is exponentially decayed regardless of the focus degrading. Collectively, the DMA can effectively battle those defects inherently embedded in the CSA and GA, which are also shared by most of the popular algorithms.

## 4. DISCUSSION

As seen, the simulation and experiment results have demonstrated the high adaptability of the DMA. It performs comparably with the GA in a noisy environment and overcomes the heavy perturbations with a robust recovery ability. Apart from this distinctive adaptability, the ease of implementation is another advantage of the DMA. For the GA, several key parameters, such as population size, offspring size, and mutation rates, need to be adjusted appropriately at the beginning [28]. However, in DMA, only one parameter is required to be set, which is the mutation constant, a number bounding the mutation rate. Benefitting from the error rate and the square rule, the optimization error in the modulating mask can be easily quantified, causing the whole optimization process to be straightforward. It is simply based on real-time measurements, and it is less likely to be affected by the improper selection of parameters. Also, the error-rate-based DMA does not strongly depend on the modulating mask in previous iterations, and therefore adaptability can be effectively achieved to deal with dynamic media. As an example of the square rule’s applications, the DMA does show its capability to adapt to strong and/or sudden perturbations benefitting from the use of the practical metric, the error rate. Note that the metric can also be incorporated into other optimization methods, such as the GA, to improve adaptability.

Although only one of the numerical optimization methods, i.e., the GA, is chosen for demonstration in this study, others, such as the SA and PSO, are similar. Including in their pools a large population of masks, the mechanisms to generate new modulating patterns originate from the philosophy of numerical optimization: (1) the portion of the modulating elements to be mutated is preset or decays with certain rules; (2) the monitored PBR is used to produce an acceptance probability of the newly generated masks or to update the mask. These two mechanisms ensure the generalization to the noise, even on the scale of $\u27e8{I}_{0}\u27e9$ [33], around a specific equilibrium of the medium. However, a sudden strong perturbation behaves differently since the equilibrium of the medium has been considerably changed, and the experience based on the previous equilibrium is “out-of-date” as discussed in the last section. Therefore, it can be considered that the physics-based square rule probably enhances the adaptability of the existing methods. Notably, if a new parameter is incorporated in those methods, other parameters probably need to be further tuned to match the function of the square rule, which is a nontrivial manipulation. Incorporation of other optimization algorithms with the square rule is therefore beyond the scope of this paper.

## 5. CONCLUSION

In this study, a simple square rule of binary-amplitude-modulation-based wavefront shaping optical focusing based on universal strong scattering media has been theoretically obtained. With this rule, the real-time error in the modulating mask can be simply calculated from the concurrently measured PBR of the optical focus. Based on such a real-time metric, a novel feedback-based wavefront shaping algorithm, the DMA, has been proposed. Both the simulation and experimental results have demonstrated its high adaptability and unique recovery ability that no other existing algorithms can achieve: focusing of diffused light can be regained without re-running the optimization even after a 60% drop of the PBR. This is due to the application of the square rule, which guides the optimization with universal physics knowledge about the strong scattering process instead of a random guess. Notably, the square rule assumes that the transmission matrix of a medium follows a circular Gaussian distribution. It can be easily fulfilled when the transmitted medium is a strong scattering medium: photons are multiply scattered, and most of these scattering events are independent [39]. Therefore, the square rule between the DMD error rate and the degraded focus can be generally applied to the process in a strong scattering regime. The algorithm is therefore particularly suitable to be used in heavily instable or motional scattering environments. Note that the DMA in this study merely serves as an example application to utilize the error rate and square rule to optimize the single-point focusing. On the other hand, MMF is used as the example of scattering media in this study, so that by applying rotation of certain degrees we can induce controllable, repeatable, and quantifiable perturbations to the resultant speckle patterns. This is necessary for the current phase of the proof of principle, although it is not ideal. With further improvement, we believe the study may boost or inspire many applications of wavefront shaping with instable media or even living biological tissue.

## Funding

National Key Research and Development Program of China (2017YFA0700401); National Natural Science Foundation of China (81627805, 81671726, 81827808, 81930048); Research Grants Council, University Grants Committee (25204416); Innovation and Technology Commission (GHP/043/19SZ, GHP/044/19GD, ITS/022/18); Guangdong Science and Technology Department (2019A1515011374, 2019BT02X105); Science, Technology and Innovation Commission of Shenzhen Municipality (JCYJ20170818104421564); Youth Innovation Promotion Association of the Chinese Academy of Sciences (2018167).

## Disclosures

The authors declare no conflicts of interest.

## REFERENCES

**1. **I. M. Vellekoop and A. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett. **32**, 2309–2311 (2007). [CrossRef]

**2. **E. N. Leith and J. Upatnieks, “Holographic imagery through diffusing media,” J. Opt. Soc. Am. **56**, 523 (1966). [CrossRef]

**3. **Z. Yaqoob, D. Psaltis, M. S. Feld, and C. Yang, “Optical phase conjugation for turbidity suppression in biological samples,” Nat. Photonics **2**, 110–115 (2008). [CrossRef]

**4. **Z. Yu, J. Huangfu, F. Zhao, M. Xia, X. Wu, X. Niu, D. Li, P. Lai, and D. Wang, “Time-reversed magnetically controlled perturbation (TRMCP) optical focusing inside scattering media,” Sci. Rep. **8**, 2927 (2018). [CrossRef]

**5. **Z. Yu, M. Xia, H. Li, T. Zhong, F. Zhao, H. Deng, Z. Li, D. Li, D. Wang, and P. Lai, “Implementation of digital optical phase conjugation with embedded calibration and phase rectification,” Sci. Rep. **9**, 1537 (2019). [CrossRef]

**6. **Y. Liu, P. Lai, C. Ma, X. Xu, A. A. Grabar, and L. V. Wang, “Optical focusing deep inside dynamic scattering media with near-infrared time-reversed ultrasonically encoded (TRUE) light,” Nat. Commun. **6**, 5904 (2015). [CrossRef]

**7. **S. Popoff, G. Lerosey, R. Carminati, M. Fink, A. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett. **104**, 100601 (2010). [CrossRef]

**8. **P. Lai, L. Wang, J. W. Tay, and L. V. Wang, “Photoacoustically guided wavefront shaping (PAWS) for enhanced optical focusing in scattering media,” Nat. Photonics **9**, 126–132 (2015). [CrossRef]

**9. **Y. Luo, S. Yan, H. Li, P. Lai, and Y. Zheng, “Focusing light through scattering media by reinforced hybrid algorithms,” APL Photon. **5**, 016109 (2020). [CrossRef]

**10. **S. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Image transmission through an opaque material,” Nat. Commun. **1**, 81 (2010). [CrossRef]

**11. **S. Ohayon, A. Caravaca-Aguirre, R. Piestun, and J. J. DiCarlo, “Minimally invasive multimode optical fiber microendoscope for deep brain fluorescence imaging,” Biomed. Opt. Express **9**, 1492–1509 (2018). [CrossRef]

**12. **I. M. Vellekoop, “Feedback-based wavefront shaping,” Opt. Express **23**, 12189–12206 (2015). [CrossRef]

**13. **A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nat. Photonics **6**, 283–292 (2012). [CrossRef]

**14. **Z. Li, Z. Yu, H. Hui, H. Li, T. Zhong, H. Liu, and P. Lai, “Edge enhancement through scattering media enabled by optical wavefront shaping,” Photon. Res. **8**, 954–962 (2020). [CrossRef]

**15. **R. Horstmeyer, H. Ruan, and C. Yang, “Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue,” Nat. Photonics **9**, 563–571 (2015). [CrossRef]

**16. **H. Ruan, T. Haber, Y. Liu, J. Brake, J. Kim, J. M. Berlin, and C. Yang, “Focusing light inside scattering media with magnetic-particle-guided wavefront shaping,” Optica **4**, 1337–1343 (2017). [CrossRef]

**17. **N. Takai and T. Asakura, “Statistical properties of laser speckles produced under illumination from a multimode optical fiber,” J. Opt. Soc. Am. A **2**, 1282–1290 (1985). [CrossRef]

**18. **M. Plöschner, T. Tyc, and T. Čižmár, “Seeing through chaos in multimode fibres,” Nat. Photonics **9**, 529–535 (2015). [CrossRef]

**19. **Y. Choi, C. Yoon, M. Kim, T. D. Yang, C. Fang-Yen, R. R. Dasari, K. J. Lee, and W. Choi, “Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber,” Phys. Rev. Lett. **109**, 203901 (2012). [CrossRef]

**20. **I. N. Papadopoulos, S. Farahi, C. Moser, and D. Psaltis, “High-resolution, lensless endoscope based on digital scanning through a multimode optical fiber,” Biomed. Opt. Express **4**, 260–270 (2013). [CrossRef]

**21. **J. Yoon, M. Lee, K. Lee, N. Kim, J. M. Kim, J. Park, H. Yu, C. Choi, W. Do Heo, and Y. Park, “Optogenetic control of cell signaling pathway through scattering skull using wavefront shaping,” Sci. Rep. **5**, 13289 (2015). [CrossRef]

**22. **A. M. Aravanis, L.-P. Wang, F. Zhang, L. A. Meltzer, M. Z. Mogri, M. B. Schneider, and K. Deisseroth, “An optical neural interface: *in vivo* control of rodent motor cortex with integrated fiberoptic and optogenetic technology,” J. Neural Eng. **4**, S143–S156 (2007). [CrossRef]

**23. **O. Tzang, A. M. Caravaca-Aguirre, K. Wagner, and R. Piestun, “Adaptive wavefront shaping for controlling nonlinear multimode interactions in optical fibres,” Nat. Photonics **12**, 368–374 (2018). [CrossRef]

**24. **T. Zhong, Z. Yu, H. Li, Z. Li, and P. Lai, “Active wavefront shaping for controlling and improving multimode fiber sensor,” J. Innov. Opt. Health Sci. **12**, 1942007 (2019). [CrossRef]

**25. **J. Yang, L. Li, J. Li, Z. Cheng, Y. Liu, and L. V. Wang, “Fighting against fast speckle decorrelation for light focusing inside live tissue by photon frequency shifting,” ACS Photon. **7**, 837–844 (2020). [CrossRef]

**26. **I. M. Vellekoop and A. Mosk, “Phase control algorithms for focusing light through turbid media,” Opt. Commun. **281**, 3071–3080 (2008). [CrossRef]

**27. **J. Thompson, B. Hokr, and V. Yakovlev, “Optimization of focusing through scattering media using the continuous sequential algorithm,” J. Mod. Opt. **63**, 80–84 (2016). [CrossRef]

**28. **D. B. Conkey, A. N. Brown, A. M. Caravaca-Aguirre, and R. Piestun, “Genetic algorithm optimization for focusing through turbid media in noisy environments,” Opt. Express **20**, 4840–4849 (2012). [CrossRef]

**29. **D. Wu, J. Luo, Z. Li, and Y. Shen, “A thorough study on genetic algorithms in feedback-based wavefront shaping,” J. Innov. Opt. Health Sci. **12**, 1942004 (2019). [CrossRef]

**30. **B. R. Anderson, P. Price, R. Gunawidjaja, and H. Eilers, “Microgenetic optimization algorithm for optimal wavefront shaping,” Appl. Opt. **54**, 1485–1491 (2015). [CrossRef]

**31. **H.-L. Huang, Z.-Y. Chen, C.-Z. Sun, J.-L. Liu, and J.-X. Pu, “Light focusing through scattering media by particle swarm optimization,” Chin. Phys. Lett. **32**, 104202 (2015). [CrossRef]

**32. **B.-Q. Li, B. Zhang, Q. Feng, X.-M. Cheng, Y.-C. Ding, and Q. Liu, “Shaping the wavefront of incident light with a strong robustness particle swarm optimization algorithm,” Chin. Phys. Lett. **35**, 124201 (2018). [CrossRef]

**33. **Z. Fayyaz, N. Mohammadian, M. Reza Rahimi Tabar, R. Manwar, and K. Avanaki, “A comparative study of optimization algorithms for wavefront shaping,” J. Innov. Opt. Health Sci. Sci. **12**, 1942002 (2019). [CrossRef]

**34. **L. Fang, H. Zuo, Z. Yang, X. Zhang, J. Du, and L. Pang, “Binary wavefront optimization using a simulated annealing algorithm,” Appl. Opt. **57**, 1744–1751 (2018). [CrossRef]

**35. **Z. Fayyaz, F. Salimi, N. Mohammadian, A. Fatima, M. R. R. Tabar, and M. R. Avanaki, “Wavefront shaping using simulated annealing algorithm for focusing light through turbid media,” Proc. SPIE **10494**, 104946M (2018). [CrossRef]

**36. **S. Cheng, H. Li, Y. Luo, Y. Zheng, and P. Lai, “Artificial intelligence-assisted light control and computational imaging through scattering media,” J. Innov. Opt. Health Sci. **12**, 1930006 (2019). [CrossRef]

**37. **R. Horisaki, R. Takagi, and J. Tanida, “Learning-based focusing through scattering media,” Appl. Opt. **56**, 4358–4362 (2017). [CrossRef]

**38. **A. Turpin, I. Vishniakou, and J. D. Seelig, “Light scattering control in transmission and reflection with neural networks,” Opt. Express **26**, 30911–30929 (2018). [CrossRef]

**39. **J. W. Goodman, *Speckle Phenomena in Optics: Theory and Applications* (Roberts and Company, 2007).

**40. **D. Akbulut, T. J. Huisman, E. G. van Putten, W. L. Vos, and A. P. Mosk, “Focusing light through random photonic media by binary amplitude modulation,” Opt. Express **19**, 4017–4029 (2011). [CrossRef]

**41. **Y. Shen, Y. Liu, C. Ma, and L. V. Wang, “Sub-Nyquist sampling boosts targeted light transport through opaque scattering media,” Optica **4**, 97–102 (2017). [CrossRef]

**42. **H. Yu, K. Lee, and Y. Park, “Ultrahigh enhancement of light focusing through disordered media controlled by mega-pixel modes,” Opt. Express **25**, 8036–8047 (2017). [CrossRef]