⇦Chinese Physics Letters, 2020, Vol. 37, No. 2, Article code 024202 Wavefront Shaping for Fast Focusing Light through Scattering Media Based on Parallel Wavefront Optimization and Superpixel Method * Yingchun Ding (丁迎春)1, Xinjing Lv (吕新晶)1, Youquan Jia (贾佑权)1, Bin Zhang (张彬)2, Zhaoyang Chen (陈朝阳)1**, Qiang Liu (柳强)2** Affiliations 1College of Mathematics and Physics, Beijing University of Chemical Technology, Beijing 100029 2State Key Laboratory of Precision Measurement Technology and Instruments, Department of Precision Instruments, Tsinghua University, Beijing 100084 Received 26 November 2019, online 18 January 2020 *Supported by the National Key Research and Development Program of China (Grant No. 2017YFB1104500), the Beijing Natural Science Foundation (Grant No. 7182091), the National Natural Science Foundation of China (Grant No. 21627813), and the Research Projects on Biomedical Transformation of China-Japan Friendship Hospital (PYBZ1801).
**Corresponding author. Email: chenzy@mail.buct.edu.cn; qiangliu@mail.tsinghua.edu.cn Citation Text: Ding Y C, Lv X J, Jia Y Q, Zhang B and Chen C Y et al 2020 Chin. Phys. Lett. 37 024202    Abstract When light travels in biological tissues, it undergoes multiple scattering and forms speckles, which seriously restricts the penetration depth of optical imaging in biological tissues. With wavefront shaping method, by modulating the wavefront of incident light to compensate for the wavefront aberration, light focusing and scanning imaging through scattering media can be achieved. However, wavefront shaping must be accomplished within the speckle decorrelation time. Considering the short speckle decorrelation time of living tissues, the speed of wavefront shaping is rather essential. We propose a new iterative optimization wavefront shaping method to improve the speed of wavefront shaping in which the existing parallel optimization wavefront shaping method is improved and is combined with the superpixel method. Compared with the traditional multi-frequency parallel optimization method, the modulation rate of our method is doubled. Moreover, we combine the high frame rate amplitude modulator, i.e., the digital micromirror device (DMD), with the superpixel method to replace the traditional phase modulator (i.e., spatial light modulator), which further increases the optimization speed. In our experiment, when the number of the optical modes is 400, light focusing is achieved with only 1000 DMD superpixel masks and the enhancement factor reaches 223. Our approach provides a new path for fast light focusing through wavefront shaping. DOI:10.1088/0256-307X/37/2/024202 PACS:42.25.Dd, 42.25.-p, 42.25.Fx, 42.30.Kq © 2020 Chinese Physics Society Article Text Biomedical imaging has been intensively studied over the past decades because it is a powerful tool for diagnosis and treatment of diseases such as cancer.[1] Biological tissues are usually strong scattering media in which light experiences multiple scattering events and rapidly becomes scattered light.[2] This presents a serious obstacle to the application of traditional high resolution optical imaging methods in biological tissue imaging such as confocal microscopy,[3] optical coherence tomography (OCT)[4] and two- or three-photon microscopy.[5] These methods utilize ballistic light to image and the penetration depth is seriously limited, usually one transport mean free path.[6] Wavefront shaping method[7] has been developed to achieve light focusing and scanning imaging through scattering media by modulating the wavefront phase and amplitude of the incident light to compensate for the distortion caused by scattering. The wavefront shaping methods mainly include the transmission matrix method[8,9] which is based on the random matrix theory, the iterative optimization method,[10] and optical phase conjugation (OPC).[11] Compared with the cumbersome procedures in transmission matrix measurement and the complex experimental system in optical phase conjugation, the iterative optimization wavefront shaping method is simple and not susceptible to disturbance. However, the iterative optimization process is time-consuming which limits its in vivo application because wavefront shaping must be completed within the speckle decorrelation time and the decorrelation time of biological tissues is rather short.[12] The iterative optimization process of wavefront shaping is controlled by optimization algorithm. In 2008, Vellekoop et al. proposed a sequence optimization algorithm for iterative optimization of the incident light in which one optical mode was optimized at one time.[13] To perform iterative optimization of $N$ optical modes with $m$ phase steps through the sequence optimization algorithm, $(m-1)N$ phase masks were required. In 2012, Conkey et al. introduced genetic algorithm (GA) for iterative optimization wavefront shaping.[14] Compared with the sequence optimization algorithm, the GA had the advantages of high robustness and initial faster enhancement of the focus intensity, with the expense of the greatly increased number of masks required for optimization. In 2011, Cui et al. proposed a multi-frequency parallel optimization method in which the wavefront of the incident light was determined in a parallel fashion.[15] Compared with the pixel by pixel optimization, only one-fifth of masks were needed in the parallel optimization method in theory.[14] This method could theoretically increase the optimization speed by five times. The enhancement factor is defined as the ratio of the light intensity of the focus and the average light intensity of the background.[16] In their experiment, 441 optical modes were optimized with only 2048 masks and the enhancement factor reached nearly 270. In addition to the iterative optimization algorithm, the speed of iterative optimization wavefront shaping is mainly affected by the refresh rate of the optical modulation devices. Light focusing within the speckle decorrelation time of biological tissues and scanning imaging can be achieved with high refresh rate optical modulation devices.[17] Liquid crystal spatial light modulators (SLM) are most commonly used in wavefront shaping experiments to adjust the wavefront of the incident light.[18–21] However, the low refresh rate of SLMs, usually about 100 Hz, prevents their further applications in biomedical imaging.[22] The refresh rate of digital micromirror devices (DMDs) can reach 32 kHz, which are more suitable for in vivo imaging.[23] As a binary amplitude modulation device, DMDs can not modulate the phase of incident light directly. In 2018, Jia et al. performed phase modulation of incident light with the DMD-based superpixel method and achieved light focusing through scattering media.[24] In this Letter, we propose a new iterative optimization method for wavefront shaping. We improve the multi-frequency parallel optimization algorithm, and the number of masks needed for optimization is reduced by half. Moreover, we combine the improved parallel optimization algorithm with the DMD-based superpixel method to further increase the optimization speed.
cpl-37-2-024202-fig1.png
Fig. 1. Schematic of the experimental principle. (a) Modulation frequency distribution in the first optimization. (b) Fourier phase spectrum of the collected signal. (c) The optimal phase mask corresponding to the modulated half optical modes. (f) Modulation frequency distribution in the second optimization. (e) The whole optimal phase mask obtained after the second optimization. (d) The optimal superpixel mask.
The experimental principle is shown in Fig. 1. To facilitate description, we divide the incident light into 16 optical modes in this demo and the initial phase of all the optical modes is set to 0. The optimization process is divided into two steps. Figure 1(a) shows the modulation method used in the first half of the optimization process. Half of the optical modes are phase modulated at different frequencies, i.e., the $k$th optical mode is modulated at frequency $f_{k}$, while the rest optical modes remain unchanged and serve as the reference light. Interference occurred between the modulated light and the reference light and the interference fringes shifted with the phase change of the modulated light. A detector was placed at the target point to record the change of the light intensity during the whole process. In this process, the signal collected by the detector is processed by Fourier transform, and the processing result is shown in Fig. 1(b). The initial phase difference between each modulated beam and the reference beam can be extracted from the obtained Fourier phase spectrum. Then, the optimal phase mask corresponding to the modulated optical modes can be obtained by calculation, as shown in Fig. 1(c). The modulation method used in the latter half of the optimization process is shown in Fig. 1(f). Similarly, the latter half of the optical modes are phase modulated at different frequencies and the corresponding phase information is extracted by Fourier transform. Finally, the whole optimal phase mask corresponding to all the optical modes is shown in Fig. 1(e). Figure 1(d) shows the superpixel mask obtained by translating the optimal mask through the look-up table method.[24] When the optimal mask is projected onto the optical modulation device, the wavefront of the incident light is phase modulated and light focusing can be achieved at the target point. Multi-frequency parallel optimization method utilizes Fourier transform to extract phase information since there is no interference between signals of different frequencies in the frequency domain. When the light modulated at a designed frequency interferes with the reference light, the light intensity recorded by the detector is a harmonic signal which can be expressed as $$ I=A^{2}+A_{\rm R}^{2} +2AA_{\rm R} \cos (f\times t+\phi),~~ \tag {1} $$ where $A$ denotes the amplitude of the modulated light, $A_{\rm R}$ denotes the amplitude of the reference light, $f$ denotes the modulation frequency and $\phi$ denotes the phase difference between the reference light and the modulated light. We perform Fourier transform on the signal $I$ and obtain the initial phase corresponding to the modulation frequency $f$ in the Fourier phase spectrum. When the optical modes of the incident light are modulated independently at different frequencies, the signal captured by the detector is the superposition of these harmonic signals, which are independent of each other in the frequency domain. A simple example is given to illustrate the process of signal demodulation. The signal $I_{0}$ contains two different frequency components and can be express as $I_{0} = \cos[2\pi t+(\pi /4)]+\cos[4\pi t+(\pi /2)]$. As shown in Figs. 2(a) and 2(b), the Fourier amplitude spectrum and phase spectrum obtained by Fourier transform contain the left-right symmetrical positive frequency components and negative frequency components. In the figures, the components of positive frequencies are marked in blue and the components of negative frequencies are marked in red. The frequency components represented by the two peaks on the left half of the amplitude spectrum are consistent with the signal $I_{0}$, and the initial phase corresponding to the two modulation frequencies can be obtained in the phase spectrum.
cpl-37-2-024202-fig2.png
Fig. 2. (a) The Fourier amplitude spectrum of signal $I_{0}$. (b) The Fourier phase spectrum of signal $I_{0}$.
To further improve the optimization speed, we combined the DMD with the superpixel method to modulate the wavefront of the incident light. In the superpixel method, the DMD is divided into separate superpixel blocks, each of which is composed of $n\times n$ subpixels, and each micromirror on the DMD serves as a subpixel.[24] The light reflected from the DMD surface will pass through a 4$f$ system consisting of two lenses and a space filter. The position of the space filter is specially selected, and only the $+$1 diffraction light can pass it through. The light arriving at the position of the space filter from different subpixels of one superpixel block has different phases. Phase of the superposed superpixel block can be adjusted by changing the subpixel block combination to realize the phase modulation of the incident light with DMD. Although the frame rate of the DMD can reach 32 kHz, it is impossible to realize continuous modulation of the optical modes at different frequencies relying solely on hardware. We realize modulation of the optical modes at different frequencies in discrete form. A certain number of masks are used to control the phase distribution of the wavefront at different times. The parameters corresponding to different modulation frequencies are as follows: $$\begin{align} &y=\exp \left(i 2\pi f t \right),~~ \tag {2}\\ &t=0,i_{\rm s},2i_{\rm s},\ldots,T-i_{\rm s}, \\ &f=a,a+j,a+2j,\ldots,b, \end{align} $$ where $y$ denotes the modulated optical field, $f$ the modulation frequency, $a$ the minimum frequency, $b$ the maximum frequency, $j$ the frequency step, $t$ the sampling time, $i_{\rm s}$ the sampling time interval, and $T$ the total sampling time. To obtain accurate phase information in the Fourier phase spectrum, signals of each frequency must be full period sampled. Accordingly, $$ T\cdot j=1.~~ \tag {3} $$ The sampling frequency of the signal is defined as $$ f_{\rm s} =\frac{1}{i_{\rm s} }.~~ \tag {4} $$ In the spectrum obtained by Fourier transform, the maximum value of the frequency is $f_{\rm s}$. When the modulation frequency $f$ of the signal is less than half of the sampling frequency $f_{\rm s}$, i.e., $$ f < \frac{f_{\rm s} }{2},~~ \tag {5} $$ the components of this modulation frequency always appear in the left half of the spectrum. According to formulae (4) and (5), the modulation frequency $f$ and the sampling time interval $i_{\rm s}$ should follow the formula $$ f\cdot i_{\rm s} < 0.5.~~ \tag {6} $$ Under this condition, there is no interference between the positive frequency components and the negative frequency components that the phase information can be extracted accurately. Our experimental scheme is shown in Fig. 3. The incident light is divided into 400 optical modes, i.e., $N = 400$, and the initial phase of all the optical modes is 0. Each column corresponds to the phase distribution of the 400 optical modes at a given time, and each row corresponds to the phase of one optical mode at different times. The phase mask loaded on the DMD is marked with a red box. It starts from the far left, and as the optimization proceeds, it moves one square at one time to the right-most end. The blue region is the phase distribution on the top half of the DMD, and the yellow region is the phase distribution on the bottom half of the DMD.
cpl-37-2-024202-fig3.png
Fig. 3. Experimental scheme diagram. The DMD is divided into the dynamic half region and the stationary half region. The two regions are marked in blue and yellow, respectively. The phase mask loaded on the DMD is marked with a red box.
The optimization was implemented in two steps. As shown in Fig. 3, in the first step, the blue region is the dynamic half region and the yellow region is the stationary half region. The optical modes in the dynamic half region are frequency modulated while the optical modes locating in the stationary half region are not modulated and served as the reference light. A charge coupled device (CCD) camera is used to record the light intensity at the target point during the whole process. Considering that the two optimization processes are independent of each other, we make an improvement on the multi-frequency parallel optimization algorithm. We design 200 modulation frequencies used in two steps respectively. The relevant parameters are set as follows: $$ t=0,0.1,0.2,\ldots,49.9;~~~ f=0.5,0.52,0.54,\ldots,4.48.~~ \tag {7} $$ The advantage of this design is that the frequency step $j$ could be increased while the sampling time interval $i_{\rm s}$ remains unchanged. According to formulae (2) and (3), the total sampling time $T$ would be reduced. Finally, only 1000 masks are needed to complete the optimization.
cpl-37-2-024202-fig4.png
Fig. 4. Schematic of the experimental setup. DMD: digital micromirror device, L1, L2: lenses, OBJ1, OBJ2: objectives, S: scattering diffuser, CCD: charge coupled device.
The experimental setup is shown in Fig. 4. After beam expansion and collimation, the light irradiated on the DMD surface. The modulated light reflected from the DMD surface passes through a 4$f$ system consisting of two lenses and a space filter. We load the superpixel mask on the DMD. A superpixel block is composed of $4 \times 4$ subpixels. The position of the spatial filter is carefully selected such that only the $+$1 level diffraction light could pass through it. Then the light is focused on the scattering media through the objective and the scattered light is captured by an 8-bit CCD camera. A PC was used to load the superpixel masks onto the DMD and to process the light field distribution maps captured by CCD. After successfully simulating the proposed experimental scheme with MATLAB, we construct the experiment system and conduct the wavefront shaping light focusing experiment with our proposed method. First, 200 optical modes are modulated at different frequencies and the signals collected by the detector during this process are shown in Fig. 5(a). Figures 5(c) and 5(d) are the amplitude and phase spectra of the signals extracted through Fourier transform, respectively. We modulate the remaining 200 optical modes in the same way and the measured signals are shown in Fig. 5(b). Compared with the previously collected signals, it can be seen that the intensity of the signals have been greatly improved because the light intensity of the target point has been improved after the first optimization.
cpl-37-2-024202-fig5.png
Fig. 5. (a) Signals collected in the first optimization. (b) Signals collected in the second optimization. (c) Amplitude spectrum of the signals in (b). (d) Phase spectrum of the signals in (b).
cpl-37-2-024202-fig6.png
Fig. 6. (a) The focusing effect when the first optimization is accomplished. The inset presents the superpixel mask obtained when the first half of the optimization process is completed. (b) The focusing effect when the whole optimization process is completed. The inset presents the superpixel mask obtained when the whole optimization process is completed.
The experimental focusing effect is shown in Fig. 6. Figure 6(a) shows the focusing effect when the first half of the optimization process is completed, and the enhancement factor is 72. Figure 6(b) shows the focusing effect when the whole optimization process is completed with an enhancement factor of 223. In the simulation, the enhancement factors obtained in the two optimization processes are about 90 and 270, respectively. The experimental results are highly consistent with the simulation results. Affected by ambient noise and restricted by experimental equipment, the final enhancement factor is slightly lower than the theoretical value. In summary, we have proposed a new method of multi-frequency parallel optimization. The optimization process is divided into two parts. The first half and the second half of the optimization process would not affect each other, so a set of modulation parameters could be used together. Through the precise design of the modulation scheme, the number of masks used is greatly reduced. In the experiment, when the number of the optical modes is 400, only 1000 masks are needed to complete the optimization, and the enhancement factor of the focus is close to the theoretical value. Moreover, we have combined the superpixel method with a DMD to modulate the phase of the incident light. Benefiting from the extremely high frame rate of DMDs, the optimization speed could be further improved. Our proposed method has been proved to have faster speed of iterative optimization wavefront shaping and will promote the application of wavefront shaping in living tissue imaging. References Optical properties of biological tissues: a reviewLine-scanning reflectance confocal microscopy of human skin: comparison of full-pupil and divided-pupil configurationsOptical coherence tomographyMulti-photon microscopy with a low-cost and highly efficient Cr:LiCAF laserGoing deeper than microscopy: the optical imaging frontier in biologyFocusing coherent light through opaque strongly scattering mediaMeasuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered MediaLabel-Free Microscopic Imaging Based on the Random Matrix Theory in Wavefront Shaping *Research on intelligent algorithms for amplitude optimization of wavefront shapingOptical phase conjugation for turbidity suppression in biological samplesRelation between speckle decorrelation and optical phase conjugation (OPC)-based turbidity suppression through dynamic scattering media: a study on in vivo mouse skinPhase control algorithms for focusing light through turbid mediaGenetic algorithm optimization for focusing through turbid media in noisy environmentsParallel wavefront optimization method for focusing light through random scattering mediaFeedback-based wavefront shapingBreaking the spatial resolution barrier via iterative sound-light interaction in deep tissue microscopyRecent advances in wavefront shaping techniques for biomedical applicationsHigh-efficiency and flexible generation of vector vortex optical fields by a reflective phase-only spatial light modulatorSpatial Light Modulator Based on Surface Plasmon Polaritons for Chromatic DisplayFast generation of controllable equal-intensity four beams based on isosceles triangle multilevel phase grating realized by liquid crystal spatial light modulatorPhotoacoustic Wavefront Shaping with High Signal to Noise Ratio for Light Focusing Through Scattering MediaDigital micromirror device-based ultrafast pulse shaping for femtosecond laserSuperpixel-Based Complex Field Modulation Using a Digital Micromirror Device for Focusing Light through Scattering Media
[1]Cherry S R, Badawi R D and Qi J 2016 Essentials of in vivo Biomedical Imaging (Boca Raton: CRC Press)
[2] Jacques S L 2013 Phys. Med. & Biol. 58 R37
[3] Gareau D S, Abeytunge S and Rajadhyaksha M 2009 Opt. Lett. 34 3235
[4] Huang D, Swanson E A, Lin C P, Schuman J S, Stinson W G, Chang W, Hee M R, Flotte T, Gregory K, Puliafito C A and Fujimoto J G 1991 Science 254 1178
[5] Sakadžić S, Demirbas U, Mempel T R, Moore A, Ruvinskaya S, Boas D A, Sennaroglu A, Kartner F X and Fujimoto J G 2008 Opt. Express 16 20848
[6] Ntziachristos V 2010 Nat. Methods 7 603
[7] Vellekoop I M and Mosk A P 2007 Opt. Lett. 32 2309
[8] Popoff S M, Lerosey G, Carminati R, Fink M, Boccara A C and Gigan S 2010 Phys. Rev. Lett. 104 100601
[9] Yu L, Xu X, Zhang Z, Feng Q, Zhang B, Ding Y and Liu Q 2019 Chin. Phys. Lett. 36 114203
[10] Feng Q, Zhang B, Liu Z, Lin C and Ding Y 2017 Appl. Opt. 56 3240
[11] Yaqoob Z, Psaltis D, Feld M S and Yang C 2008 Nat. Photon. 2 110
[12] Jang M, Ruan H, Vellekoop I M, Judkewitz B, Chung E and Yang C 2015 Biomed. Opt. Express 6 72
[13] Vellekoop I M and Mosk A P 2008 Opt. Commun. 281 3071
[14] Conkey D B, Brown A N, Caravaca-Aguirre A M and Piestun R 2012 Opt. Express 20 4840
[15] Cui M 2011 Opt. Lett. 36 870
[16] Vellekoop I M 2015 Opt. Express 23 12189
[17] Park J, Sun W and Cui M 2015 Proc. Natl. Acad. Sci. USA 112 9236
[18] Yu H, Park J, Lee K R, Yoon J, Kim K D, Lee S and Park Y K 2015 Curr. Appl. Phys. 15 632
[19] Cai M, Wang Z, Liang J, Wang Y, Gao X, Li Y, Tu C and Wang H 2017 Appl. Opt. 56 6175
[20] Wang J, Huang P, Du J, Guo Y, Luo X and Du C 2008 Chin. Phys. Lett. 25 2908
[21] Liu X, Zhang J, Wu L and Gan Y 2011 Chin. Phys. B 20 024211
[22] Sun J, Zhang B, Feng Q, He H, Ding Y and Liu Q 2019 Sci. Rep. 9 4328
[23] Gu C, Zhang D, Chang Y and Chen S 2015 Opt. Lett. 40 2870
[24] Jia Y, Feng Q, Zhang B, Wang W, Lin C and Ding Y 2018 Chin. Phys. Lett. 35 054203