## Abstract

We introduce genetic algorithms (GA) for wavefront control to focus light through highly scattering media. We theoretically and experimentally compare GAs to existing phase control algorithms and show that GAs are particularly advantageous in low signal-to-noise environments.

©2012 Optical Society of America

## 1. Introduction

Recently, focusing light through turbid media has been achieved by optimizing wavefronts via a feedback system and spatial light modulators for wavefront control [1]. The resulting waves overcome the effects of multiple scattering and are capable of producing a focal spot after propagation through the scattering media. Such systems optimize the incident wavefront by dividing it into N input modes and setting the phase of each mode, using a phase mask, to maximize the intensity of a spot in the output plane. Several phase control algorithms have been introduced to optimize the phase of each input mode [2–4]. The stepwise sequential algorithm and the continuous sequential algorithm (CSA) optimize each input mode independently and are consequently slow [2]. Furthermore, by modulating a single mode the interference signal detected at the output plane is small, resulting in initial errors in measurement and phase mask selection that remain until the SNR increases. The partitioning algorithm (PA), on the other hand, modulates and maximizes the interference signal of a randomly selected half of the input modes during each iteration [2], which increases the SNR and causes a faster initial enhancement of the focal spot intensity than the sequential algorithms. However, as the PA progresses the magnitude of enhancement improvement slows and the enhancement approaches the theoretical maximum more slowly than the CSA. This algorithm achieves best performance when the number of input modes is much larger than the number of iterations performed during the sample persistence time [2,5].

Other techniques have been developed to compensate for random multiple scattering. For example, an observed transmission matrix (TM) can be measured and used to create a focus through the medium in the output plane [4]. The observed transmission matrix is measured by calculating the complex field response for a set of given input basis modes. The complex field response is measured by interference between each basis element and known phase reference beams (phase shifted by 0, π/2, π, and 3π/2) after both propagate through the scattering medium [4]. Optical phase conjugation (OPC) has also been demonstrated for reverse propagation through highly scattering materials [6–8]. More recently, a method to optimize the input modes in parallel has been introduced [3]. In principle this method has a five-fold reduction in the number of required measurements with respect to the CSA, although more samples were required in the reported experiment to overcome the noisy environments. Hence, this signals to a tradeoff by which longer integration times are required to reduce the number of measurements, which is at odds with time scarce experiments.

Focusing light through scattering materials has potential in various biomedical applications, such as photodynamic therapy [9], fluorescence imaging [10], optical trapping [11], and neuron excitation and imaging [12,13]. Focusing through turbid biological materials, which typically change on a millisecond timescale, requires faster algorithms. Moreover, the intensity measurements required in these techniques are noisier, due to sample vibrations and lower photon counts as a result of higher frame rates. Therefore, in this paper we introduce a method to speed up focusing through scattering materials in low signal-to-noise environments. The method uses a genetic algorithm (GA) approach for adaptively optimizing the phase mask. The GA technique increases the focal spot intensity faster than the CSA or PA while being less susceptible to noise errors.

## 2. Genetic algorithms applied to phase mask optimization

A GA is an optimization algorithm which uses principles inspired in nature to “evolve” toward a best solution. GAs are well-suited for large-scale optimization problems [14] and thus attractive to optimize the phase of the N input modes in the focusing task. GAs have been previously used, among many other applications, in optical pattern recognition for optimal filter design [15–17].

The GA developed for focusing through turbid media begins by generating an initial population of phase masks. Each phase mask is created by selecting each input mode value from a uniform pseudorandom distribution of phase values. Once the population of phase masks is created a cost function for each mask is measured. For this task the cost function is the intensity of a predefined spot in the output plane. The population of masks is then ranked according to the cost function, with a higher intensity receiving a higher ranking.

The algorithm then iteratively optimizes the phase masks through breeding and mutation operations. A random binary breeding template array, *T*, is created before breeding. Subsequently, two parent masks (ma and pa) are randomly selected from the population for breeding. A higher probability of selection is given to higher ranked phase masks. The input modes of the two parent masks are combined using *T* to create a new offspring according to: *Offspring = ma∙T + pa∙(1-T)*. This new mask is then mutated by randomly changing the phase of a set number of input modes. In order to prevent the algorithm from mutating too many optimized phase modes, the percentage of modes mutated, *R,* decreases as the algorithm runs and nears the optimal phase mask. In this implementation *R* was determined as *R = (R _{0} - R_{end})∙exp(-n/λ) + R_{end}*. Where

*R*is the initial mutation rate,

_{0}*R*is the final mutation rate,

_{end}*n*is measurement number, and

*λ*is the decay factor. The cost function of the offspring is measured and used to rank and add the offspring within the existing population through a generational method. For each generation, a specified number of masks, G (in this case half of the population size), is replaced by the new offspring regardless of their cost function, at the expense of the lower ranked masks of the previous generation. (See Fig. 1 .)

## 3. Modeling of GA-enabled focusing through scattering

Modeling and simulation of the focusing algorithms allows for testing and optimizing the GA, as well as comparison with previously developed algorithms. The scattering process was modeled as a linear relation between the complex field at the input modes and the complex field at the output modes [2]:

where*E*is the complex field at output mode m.

_{m}*A*is the amplitude contribution from input mode n, which is assumed to be uniform across the input plane and is defined as:

_{n}*A*.

_{n}= 1/√N*ø*is the phase at input mode n. The optical path through the scattering material is represented by an MxN transmission matrix (

_{n}*t*) generated from a circular Gaussian distribution. The transmission matrix element

*t*relates the field at input mode n to the output mode m. The transmission matrix relates the complex field at N input modes to the complex field at M output modes.

_{mn}The simulation uses the GA phase optimization algorithm to maximize the intensity of a specified output mode. The cost function measurement for the GA is defined as the intensity of the target output mode, *m*, with a given input mode phase distribution:

To simulate experimental conditions the simulation included adding Gaussian noise to the cost function measurement. Varying levels of additive Gaussian noise were simulated to test the robustness of the GA in noisy environments. To further simulate experimental conditions a perturbation to the transmission matrix was included after each measurement to model either a sample whose structure or position changed with time. The perturbation simulations were implemented following the model described in Ref [2].

## 4. Comparison of simulated algorithms

The overall possible enhancement which can be achieved by the GA is dependent on the number of input modes optimized. The enhancement, η, is defined as the ratio between the output mode intensity, I_{m}, and the initial average intensity, 〈I_{0}〉. The enhancement dependence on input modes is the same as the CSA and TM (excluding reference beams) focusing methods: η = π/4(N-1) [2,4], where N is the number of optimized input modes and the persistence time, T_{p}, is much larger than the iteration time, T_{i}. The difference between the GA and CSA method is that the GA initially converges quickly, but approaches the maximum enhancement much more slowly. The PA under its ideal conditions (N>>T_{p}/T_{i}), produces a theoretical enhancement that does not depend on the number of modes, but on the relationship between T_{p} and T_{i}: η = 0.5T_{p}/T_{i} [2]. These optimum enhancements are established for low noise environments, but the achievable maximum enhancement degrades with additional noise.

The algorithms are compared based on the number of measurements. A measurement being defined as the process of measuring the intensity of a specific output mode (*I _{m}*) with a phase mask, as in Eq. (2). The CSA with 10 phase samples performs 10 measurements per input mode. Thus, when optimizing a phase mask with N input modes, there are 10∙N measurements to optimize a full phase mask. In this section the CSA, TM, GA, and PA algorithm simulations were run through 10∙N measurements for comparison. Each algorithm used N=256 input modes with a large persistence time (T

_{p}=50,000 measurements). Section A2 includes simulations comparing all algorithms in the regime where T

_{p}is relatively small and the number of input modes for the GA and PA have been increased for best performance. The GA parameters selected for the simulation are indicated in section A3. The effect of population size and number of input modes in GAs is explored in section A1.

Figure 2
shows the enhancement of a focal spot with the various phase control algorithms for different levels of noise. Each algorithm was simulated 500 times, each time with a new transmission matrix, and then averaged for comparison. As expected, the CSA shows a quadratic improvement in enhancement. The PA initially enhances the intensity of the output mode more quickly than all others, but the rate of improvement rapidly decreases. In contrast, the simulation results indicate that with noise at 30% (Fig. 2(a)) of the initial average intensity, 〈I_{0}〉, the GA has an initially higher enhancement increase and can attain a higher overall enhancement through 10∙N measurements than the CSA, PA, and TM optimization methods. After 4∙N measurements the TM method has measured the transmission matrix and creates a focus, which is initially higher than all other iterative methods.

As the measurement noise level increases, the advantage of the GA becomes more dramatic. The CSA requires more measurements to raise the signal above the noise floor to start consistently making correct phase decisions. The inaccuracies in the measured transmission matrix increase with more noise, resulting in errors to the calculated phase mask. The PA, which modulates half of the input modes for optimization performs well in these conditions. Modulating half of the input modes increases the signal modulation strength and as such raises the modulated signal measurement above the noise level. The GA however, outperforms all algorithms at high noise. The GA is less susceptible to noise, because it measures the merit of an entire mask and as a result, its cost function measurement quickly rises above the noise floor.

In appendix sections A1 and A2 the GA performance under relatively short T_{p} conditions is explored. In section A1 the implications of the population size and number of input modes in these conditions suggests a tradeoff between population size and intensity enhancement consistency and a limitation in the enhancement imposed by the short T_{p}. Section A2 compares the CSA, TM, PA, and GA algorithms with the short T_{p} condition under increasingly noisy conditions, with an increased number of input modes for the PA and GA. As with the long T_{p} case, the GA performance is not severely affected by high noise levels, while the other algorithms performance deteriorates significantly.

## 5. Experiment

The algorithms were tested using a 425 micron thick chicken eggshell as the scattering material. The transport mean free path of the eggshell was 23.6 microns, calculated with the measured total intensity transmission. A HeNe laser and a phase-only spatial light modulator (Holoeye HEO 1080P) generated N=256 propagating input modes, that were focused with an objective lens (100X, NA = 0.75) onto the surface of the scattering sample. A second objective (20X, NA = 0.4) imaged a plane 1 mm behind the eggshell. The target was to create a focus on the opposing side of the eggshell. The cost function intensity measurements were provided by a 12-bit CCD camera (Point Grey Research Grasshopper). For comparison, each algorithm was run through 10∙N measurements, which is equivalent to the number of phase mask measurements the CSA runs through in optimizing all input modes a single time. Each measurement was obtained at a rate of approximately 4.5 Hz. We varied the measurement noise by decreasing the shutter speed of the camera. Each algorithm was run ten times for each noise level. The GA parameters selected for the experiment are shown in section A3. The results are shown in Fig. 3(b) -3(f). A movie comparing experimental results of the measurement phase masks, focal spot, and enhancement for the four algorithms at low noise can be found in the supplemental material (Media 1).

The low noise experiment (Fig. 3d) shows the evolution of the enhancement of the spot intensity with each measurement. The standard deviation of the noise was measured to be at 7 percent of the initial average intensity, 〈I_{0}〉. The experiment agreed well with simulation and shows that the GA achieves a higher enhancement after 10∙N measurements than the CSA, TM, and PA. Similar to the simulations the TM achieved a higher enhancement after 4∙N measurements than the CSA and GA after 4∙N measurements.

The higher noise experiments (Fig. 3e and f) show how the CSA and TM methods suffer with high noise. Noise in Fig. 3e was measured to be 18 percent of the average background intensity; while the noise in Fig. 3f was measured to be 25 percent of the initial average intensity. With increased measurement noise the errors in the measured transmission matrix increased and lowered the enhancement of the focus spot. The CSA struggled to increase the enhancement with higher noise. The PA and GA perform just as well in any of the noise environments, with the PA more severely affected at the highest noise levels. However the GA initially increases faster and achieves a much higher enhancement.

Our experiments varied from the simulations in two areas. First, the experimental noise was signal dependent due to vibrations in the system. In other words, a displacement that reduces the signal to half its maximum, at a given enhancement, will also reduce the signal to half of the maximum for a higher enhancement. Hence, as the intensity of the output mode increases with the wavefront optimization, the measured noise in the system also increases due to vibrations creating shifts in the egg shell position relative to the optimized wavefront [6]. In our case the vibrations were mostly caused by the LC-SLM fan. While these could be reduced, we found it instructive to analyze the behavior under this perturbation. Unlike the other algorithms, the CSA was highly affected by this noise because it measures small variations on a high peak intensity. The TM measured the field before setting the optimized wavefront, typically with low intensities, so it was not significantly affected by this vibration noise. The PA intensity was never high enough to be negatively affected. The GA showed little, if any, degradation in performance from the increased noise, which further demonstrates its advantage under high noise conditions. In fact, the overall enhancement from the GA increased with increasing noise, although this was likely a result of the vibrations lowering the average peak intensity with a longer integration time in the lower noise experiments.

A second deviation from the simulation was a sample (Chicken eggshell) with a varying persistence time most likely caused by thermal drift. The persistence time is a measure of the temporal stability of the scattering sample [2]. We have measured persistence times as long as 5000 seconds (22,500 measurements) and as short as 140 seconds (630 measurements) with the Chicken eggshell sample, with an average persistence time of 2200 seconds (10,000 measurements). This high variability in the persistence time increased the unpredictability of the overall enhancement achieved. (See Fig. 4 .)

Therefore, in the experiments the overall enhancement reached by each algorithm was lower than in the respective simulations by at least 50 percent due to these two non-ideal factors. Nevertheless, the experiments demonstrate the same trend expected from the simulations, the same relative merit of each algorithm, and further illustrate the advantage of the genetic algorithm in non-ideal conditions.

## 6. Conclusion

The GA optimizes modes in parallel which results in a fast initial increase in the intensity enhancement. In addition, the GA does not rely on detecting small changes in the interference signal as the CSA, and TM methods do, but instead measures the response of full phase masks created through breeding and mutation operations. This process reduces the negative effects of measurement error as seen in the linear solutions, like the CSA and TM methods.

The simulations and experimental results show that the genetic algorithm phase mask optimization algorithm works well in high noise environments. Under our experimental conditions with the eggshell, the GA showed an initial faster enhancement of the focus intensity and achieved a higher overall enhancement than the CSA, PA, and TM methods. The GA achieved similar results with noisy environments that severely inhibited the ability of the CSA to work and decreased enhancement using the TM method. The high noise robustness of the GA could be useful as the need for higher speed algorithms increases. Higher speed algorithms would allow for focusing through more temporally dynamic samples such as tissue.

## A1. Performance of the GA as a function of number of input modes and population size

To explore the effects of the various GA parameters with dynamic samples, we analyze the performance with a short persistence time equal to 2560 measurements, or the same number of measurements required to optimize the CSA mask once. To better understand the GA performance dependence on input modes in this environment, optimizations with 256, 1024, and 4096 input modes were simulated. The GA parameters are listed in section A3. 500 simulations were performed and averaged.

The results are shown in Fig. 5a
, with an inset showing the standard deviation of the enhancement data. As expected, the simulations indicate that more input modes results in a higher enhancement. They also indicate that the enhancement diminishes at higher N as a result of the short T_{p}, which changes the transmission matrix such that the seeding and subsequent breeding of the initial mask population degrades in quality as measured by the enhancement factor. Furthermore, by providing more degrees of freedom, using a higher number of input modes, leads to longer improvements before breaking down. The initial rate of increase is about the same for all input mode values for the first 500 measurements. Although not implemented in this report, the GA could be modified to avoid the enhancement degradation via periodic updates of the population measurements.

Further simulations were run to explore the effect of population size on the total enhancement and peak enhancement time (Fig. 5b) with the GA parameters set as listed in section A3. They reveal that a smaller population will optimize more quickly initially, but will saturate and break down sooner than larger population sets as a result of having less diversity in the population. The higher population also allows the algorithm to achieve a more consistent enhancement (Fig. 5b inset). The standard deviation of the enhancement for the 500 runs is smaller for the higher population set, thus indicating a tradeoff between enhancement speed and intensity enhancement consistency.

## A2. Comparison of algorithms with low sample persistence time

A full comparison of all algorithms was simulated under a relatively low persistence time (T_{p} = 2560 measurements) condition. For this comparison the GA used N = 1024 input modes. The PA used N = 2048 input modes, which placed it in the regime of N>>T_{p}/T_{i} where it should perform optimally [2]. The CSA and TM methods still optimized N = 256 modes. The GA parameters are shown in section A3.

As with the results shown in the long persistence time simulations, the increasing noise levels severely degrade the CSA, PA, and TM algorithm performance. Under the lowest noise case, 0.3〈I_{0}〉, the TM method, with its relatively low number of measurements can achieve a high enhancement after 1024 measurements, although this high enhancement quickly decays as the transmission matrix changes. In this implementation the PA fails to achieve the theoretical enhancement of η = 0.5T_{p}/T_{i}, despite N>> T_{p}/T_{i}. The CSA quadratically improves with the number of measurements, but at a slower rate than in the long T_{p} simulations. The GA optimizes more quickly through the first measurements than the other algorithms, before reaching a break down point just before 2000 measurements. As with the long T_{p} case the GA continues to achieve a high enhancement despite increasing noise conditions. (See Fig. 6
.)

## A3. Modeling and experimental GA parameters

The reported simulations and experiments used different GA parameters. These parameters are reported in Table 1 .

## Acknowledgments

We thankfully acknowledge support from NSF awards DGE-0801680 and ECCS-1028714.

## References and links

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

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

**3. **M. Cui, “Parallel wavefront optimization method for focusing light through random scattering media,” Opt. Lett. **36**(6), 870–872 (2011). [CrossRef] [PubMed]

**4. **S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. 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**(10), 100601 (2010). [CrossRef] [PubMed]

**5. **I. M. Vellekoop and C. M. Aegerter, “Focusing light through living tissue,” Proc. SPIE **7554**, 755430, 755430-10 (2010). [CrossRef]

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

**7. **M. Cui and C. Yang, “Implementation of a digital optical phase conjugation system and its application to study the robustness of turbidity suppression by phase conjugation,” Opt. Express **18**(4), 3444–3455 (2010). [CrossRef] [PubMed]

**8. **C. L. Hsieh, Y. Pu, R. Grange, and D. Psaltis, “Digital phase conjugation of second harmonic radiation emitted by nanoparticles in turbid media,” Opt. Express **18**(12), 12283–12290 (2010). [CrossRef] [PubMed]

**9. **T. J. Dougherty, C. J. Gomer, B. W. Henderson, G. Jori, D. Kessel, M. Korbelik, J. Moan, and Q. Peng, “Photodynamic therapy,” J. Natl. Cancer Inst. **90**(12), 889–905 (1998). [CrossRef] [PubMed]

**10. **I. M. Vellekoop and C. M. Aegerter, “Scattered light fluorescence microscopy: imaging through turbid layers,” Opt. Lett. **35**(8), 1245–1247 (2010). [CrossRef] [PubMed]

**11. **T. Čižmár, M. Mazilu, and K. Dholakia, “In situ wavefront correction and its application to micromanipulation,” Nat. Photonics **4**(6), 388–394 (2010). [CrossRef]

**12. **B. A. Wilt, L. D. Burns, E. T. Wei Ho, K. K. Ghosh, E. A. Mukamel, and M. J. Schnitzer, “Advances in light microscopy for neuroscience,” Annu. Rev. Neurosci. **32**(1), 435–506 (2009). [CrossRef] [PubMed]

**13. **V. Nikolenko, K. E. Poskanzer, and R. Yuste, “Two-photon photostimulation and imaging of neural circuits,” Nat. Methods **4**(11), 943–950 (2007). [CrossRef] [PubMed]

**14. **R. L. Haupt and S. E. Haupt, *Practical Genetic Algorithms,* 2nd ed. (John Wiley & Sons, 2004).

**15. **U. Mahlab, J. Shamir, and H. J. Caulfield, “Genetic algorithm for optical pattern recognition,” Opt. Lett. **16**(9), 648–650 (1991). [CrossRef] [PubMed]

**16. **D. B. Conkey, A. Brown, A. Caravaca, and R. Piestun, "Genetic algorithm optimization of phase masks for focusing light through turbid media," in *Novel Techniques in Microscopy*, OSA Technical Digest (CD) (Optical Society of America, 2011), paper NTuA5.

**17. **O. Katz, E. Small, Y. Bromberg, and Y. Silberberg, “Focusing and compression of ultrashort pulses through scattering media,” Nat. Photonics **5**(6), 372–377 (2011). [CrossRef]