13-fold resolution gain through turbid layer via translated unknown speckle illumination

: Fluorescence imaging through a turbid layer holds great promise for various biophotonics applications. Conventional wavefront shaping techniques aim to create and scan a focus spot through the turbid layer. Finding the correct input wavefront without direct access to the target plane remains a critical challenge. In this paper, we explore a new strategy for imaging through turbid layer with a large field of view. In our setup, a fluorescence sample is sandwiched between two turbid layers. Instead of generating one focus spot via wavefront shaping, we use an unshaped beam to illuminate the turbid layer and generate an unknown speckle pattern at the target plane over a wide field of view. By tilting the input wavefront, we raster scan the unknown speckle pattern via the memory effect and capture the corresponding low-resolution fluorescence images through the turbid layer. Different from the wavefront-shaping-based single-spot scanning, the proposed approach employs many spots (i.e., speckles) in parallel for extending the field of view. Based on all captured images, we jointly recover the fluorescence object, the unknown optical transfer function of the turbid layer, the translated step size, and the unknown speckle pattern. Without direct access to the object plane or knowledge of the turbid layer, we demonstrate a 13-fold resolution gain through the turbid layer using the reported strategy. We also demonstrate the use of this technique to improve the resolution of a low numerical aperture objective lens allowing to obtain both large field of view and high resolution at the same time.


Introduction
Imaging through turbid layer holds great promise for many biophotonics applications. Various approaches have been reported in recent years, including wavefront shaping techniques [1][2][3][4][5][6][7][8][9], measurement of the transmission matrix [10], object recovery from its autocorrelation [11,12], accumulation of single-scattered waves [13], among others. In the context of wavefront shaping, a common strategy is to generate a pre-distorted wavefront that creates focus at the target plane. Since there is no direct access to the target plane, different guide stars can be used as reference beacons for wavefront shaping. Once the correct wavefront is found, the focus spot is raster scanned at the target plane using the optical memory effect [14][15][16]. Despite exciting progress on the development of wavefront shaping techniques, finding the correct wavefront without direct access to the target plane remains a challenge. In addition, the field of view is limited by the angular [14] or the translational range [15] of the single focused spot.
In this work, we explore a new strategy for improving imaging resolution through a turbid layer with a large field of view using speckle illumination and iterative recovery. Our setup is shown in Fig. 1(a), where a fluorescence sample is sandwiched between two turbid layers. Instead of generating one focus spot on the sample via wavefront shaping, we use an unshaped beam to illuminate the turbid layer and generate an unknown speckle pattern on the sample. By tilting the input wavefront, we then raster scan the unknown speckle pattern via the memory effect and capture the corresponding low-resolution fluorescence images through the turbid layer (Visualization 1). Different from the wavefront-shaping-based single-spot scanning, the proposed approach employs many spots (i.e., speckles) in parallel. Based on all captured images, we jointly recover the fluorescence object and the unknown speckle pattern. Without direct access to the object plane or knowledge of the illumination pattern, we achieve one order of magnitude resolution enhancement using the reported strategy ( Fig. 1(b)-1(c), will be discussed later). The use of illumination patterns to encode high-resolution information into low-resolution measurements is well-known and has been demonstrated in various types of structured illumination (SI) setups [16][17][18][19][20][21]. Typical linear SI techniques are targeted at 2-fold resolution improvement with known system point spread function (PSF). With certain support constraints, 3-fold resolution gain has been reported [22]. However, getting at least one order of magnitude resolution improvement with neither direct access to the targeted object nor the PSF has not been addressed before. In this work we explore such super resolution by relying on the memory effect which allows to obtain many images from a single unknown speckle pattern. This strategy is different from previous SI demonstrations which uses multiple speckle patterns and allows us to achieve more than one order of magnitude resolution gain without direct access to the object plane.
This paper is structured as follows: in Section 2, we discuss the forward modeling and recovery methods of the reported scheme. Section 3 reports experimental results and demonstrates a 13-fold resolution gain through turbid layers. In Section 4, we discuss the use of the proposed scheme to improve the resolution of a regular fluorescence microscope platform. We demonstrate the use of a 0.1 numerical aperture (NA) objective lens to achieve the resolution of a 0.4 NA. Finally, we summarize the results and discuss future directions in Section 5.

Modeling and simulation
We model the effect of the turbid layer in Fig. 1(a) as an unknown low-pass filter, with the point-spread-function (PSF) denoted as PSF(x, y). In the acquisition process, the captured image can be expressed as where I n (x, y) is the n th fluorescence intensity measurement (n = 1,2,3…), Object(x, y) is the fluorescence object we wish to recover, P unknown (x, u) is the unknown illumination pattern on the sample, (x n , y n ) is the n th positional shift of the illumination pattern, and '*' stands for convolution. In our experiment, we tilt the incident angle of the laser light to shift the unknown speckle pattern to different positions on the sample. For each position of the speckle pattern, the resulting fluorescence signal is detected through the turbid media (i.e., convolved with the PSF). The goal of our imaging scheme is to recover Object(x, y) from many measurements I n (x, y) (n = 1, 2, 3…), a problem referred to as blind deconvolution. In Eq. (1), there are four unknown terms in the right-hand side: the fluorescence object Object(x, y), the illumination pattern generated by the turbid layer P unknown (x, y), the PSF of the turbid layer PSF(x, y), and the step size of the positional shift. To jointly recover the first three unknown terms, we seek to minimize the cost function where 'a', 'b' and 'c' are step sizes. This stochastic gradient descent scheme is able to accelerate the recovery process [23,24] and it is similar to our previous demonstrations of jointly recovering the object and speckle patterns [25,26]. In our implementation, we choose where FT denotes the Fourier transform. These step sizes are chosen based on the ptychographic algorithm [27] and are related to Lipschitz constants [28]. A detailed implementation can be found in the Appendix. We initialize the object by averaging all measurements, i.e., is the total number of acquired images. We initialize the unknown pattern by setting it to an all-ones matrix. To initialize the unknown PSF(x, y), we first perform a Fourier transform on the initialized object. We estimate the cutoff frequency f cutoff to be the point where the Fourier spectrum intensity drops to 5% of its maximum. The PSF is then initialized as an Airy function with a cutoff frequency f cutoff . To obtain the step size of the positional shift, we assume the step size is the same for all measurements so that we only need to recover one parameter. In our implementation, we define the convergence index as the difference between the measurements and the generated low-resolution data from the recovery [29]. We then iterate over different step sizes and pick the one that generates the highest convergence index. Figure 2 summarizes the recovery process.   Fig. 3(b) and 3(c1)-(c3). The resolution improvement factor is shown in Fig. 3(e), where the resolution gain increases as the speckle pattern cutoff frequency increases. In our imaging setting, the detection NA is much smaller than the speckle NA, and thus, the final achievable resolution is determined by the speckle NA. As shown in Fig. 3(e), we can achieve one order of magnitude resolution gain without knowledge of the speckle pattern or the system PSF.
For the unknown incoherent PSF, we assume it is shift-invariant across the entire field of view. As such, we only need to update one PSF in each iteration. In our implementation, the PSF updating process was performed in the Fourier domain, i.e., we updated the optical transfer function (OTF) of the turbid layer in each iteration. Figure 4(a) shows the groundtruth recovery assuming the OTF is known. Figure 4(b) and 4(c) show a comparison between the cases with and without the OTF updating process. In Fig. 4(b), we set the initial OTF cutoff frequency to be 1.5 times larger than that of the ground truth. By updating the OTF in the gradient descent process, we can recover both the object and the correct OTF in Fig. 4(b1) and 4(b2). Figure 4(b3) and 4(b4) show the results without updating the OTF. Similarly, we set the initial OTF cutoff frequency to be half of that of the ground truth in Fig. 4(c). Figure  4(c1) and 4(c2) show the results with the OTF updating process. Figure 4(c3) and 4(c4) show comparisons without updating the OTF.  We set the cutoff frequency of the initial OTF to be 1.5 times larger than that of the ground truth. The results with (b1-b2) and without (b3-b4) the OTF updating process. (c) We set the cutoff frequency of the initial OTF to be half of that of the ground truth. The results with (c1-c2) and w/o (c3-c4) the OTF updating process. Figure 5 shows the process for recovering the step size of the positional shift. Figure  5(a1)-5(a3) shows the recovered object image with different step size errors. In Fig. 5(b1), the convergence index [29] increases as the loop number increases. To recover the step size, we iterate over different step sizes and pick the one that generates the highest convergence index, as shown in Fig. 5(b2).
In Fig. 6, we analyze the performance of the reported scheme with respect to the number of translated patterns and different noise levels. Figures 6(a1)-6(a4) show the recovered object using different numbers of translated patterns. Figure 6(a5) quantifies the result using mean square error (MSE) and structural similarity (SSIM) index. We can see that more translated patterns lead to improved image quality of the reconstruction (lower MSE and higher SSIM). The achievable resolution, on the other hand, is determined by the speckle size and remains the same in Fig. 6(a1)-6(a4). In Fig. 6(b), we analyze the effect of additive noise. Different amounts of speckle noise are added into the simulated raw images. The reconstructed images are shown in Fig. 6(b1)-6(b4). The corresponding MSEs and SSIMs are quantified in Fig. 6(b5). Evidently, the reported scheme is robust against additive noise.

13-fold resolution gain through turbid layers
Next we consider experimental results using the setup in Fig. 1(a), where the fluorescence object is sandwiched between two turbid layers (two scotch tapes as shown in the inset of Fig.  1(a)). A 532-nm laser diode is coupled to a single mode fiber and illuminates the object with turbid layers. The fiber is 1 cm away from the first turbid layer. The distance between the object and the first turbid layer is ~8 cm and the distance between the object and the second turbid layer is about ~2 cm. The 532-nm excitation light forms a random speckle pattern on the sample through the first turbid layer. The resulting fluorescence emission is low-pass filtered by the second turbid layer and detected by a camera with a 550-nm long-pass filter and a photographic lens (Nikon 50mm f/1.2). In our implementation, we mechanically moved the fiber to different positions with a 0.5-µm motion step size and generated slightly tilted wavefronts for illumination. Based on the memory effect, the tilted wavefront laterally translates the unknown speckle pattern on the fluorescence object (Visualization 1), resulting in the forward imaging model of Eq. (1).
In the first experiment, we use a fluorescence USAF target as the object. Figure 7(a) shows the raw fluorescence image through the turbid layer. Clearly no detail can be resolved from the fluorescence USAF resolution target. Figures 7(b) and 7(c) show the recovered image of the resolution target and Fig. 7(d) shows the recovered speckle patterns. Based on Fig. 7(a) and 7(c4), we achieve 13-fold resolution gain through the turbid layer using the reported strategy. In the second experiment, we draw some lines on a fluorescence microsphere slide and use it as the object. Figure 8(a) depicts the raw fluorescence image through the turbid layer. Figure 8

Wide-field, high-resolution fluorescence imaging
For many microscopy imaging applications, it is important to get both wide field of view and high resolution at the same time [30]. The reported scheme provides a potential solution to achieve this goal, as we show next. In Fig. 9(a), we use a low-NA objective lens (4X, 0.1 NA) to acquire raw images of a sample and use a high-NA condenser lens (0.9 NA) to generate a speckle pattern on the object. A diffuser is placed at the back focal plane of the condenser lens in the setup. We also place an aperture stop at the center of the back focal plane of the condenser lens to block the direct-transmitted light to the fluorescence detection system (Fig. 9(a)). In this experiment, we translate the sample to 30 by 30 different positions with 0.3 µm step size and capture the corresponding images through the 0.1 NA objective lens. There is no scattering layer between the sample and the objective lens. The captured images are then used to recover both the high-resolution object and the unknown speckle pattern. Figure 9(b1) shows the captured raw image of a fluorescence resolution target. Figure 9(b2) and (b3) show the recovered object and the speckle pattern. The resolution of the recovered images corresponds to a NA of 0.4. The resolution improvement factor in this experiment is 4, limited by the precision of the employed motorized stage (Newport LTA-HS).
In Fig. 9(c), we demonstrate the use of the reported scheme to recover a wide-field, highresolution fluorescence image of a microsphere sample. Figure 9(c1) depicts the recovered image where the field of view is determined by the employed 4X objective. Figure 9(c2) and 9(c3) show the recovered object image and the speckle pattern and Fig. 9(c4) shows the raw image captured using the 0.1 NA objective lens.

Discussion
In summary, we report a new strategy for improving fluorescence imaging resolution through a turbid layer with a large field of view. Instead of generating one focus spot on the sample via wavefront shaping, we use an unshaped beam to illuminate the turbid layer and generate an unknown speckle pattern on the sample. By tilting the input wavefront, we raster scan the unknown speckle pattern via the memory effect and capture the corresponding low-resolution fluorescence images through the turbid layer. Without direct access to the object plane or knowledge of the illumination pattern, we achieve 13-fold resolution gain using the reported strategy. While the idea of structured illumination is well-known for improving resolution in microscopy platform, previous demonstrations, to the best of our knowledge, are targeted at 2-fold resolution gain with a known system PSF. In this work we explore super resolution by relying on the memory effect which allows to obtain many images from a single unknown speckle pattern. This strategy is different from previous SI demonstrations which uses multiple speckle patterns [18,21] and allows us to achieve more than one order of magnitude resolution gain without direct access to the object plane. Our work shares some roots with Refs [11,12]. which recover the object from its autocorrelation measurements. In our work, we use a joint object-pattern recovery scheme to reconstruct both the object and the pattern. The scattering layer may not be needed in our scheme and we demonstrate the use of a 0.1 NA objective lens to obtain a 0.4-NA resolution with a large field of view.
The limitations of the reported strategy are threefold. First, it is based on the traditional memory effect where the target plane is at a distance from the diffusing layer. For imaging inside thick scattering media, we may need to exploit the translational memory effect [15] for the reported scheme. Second, we need a large number of acquisitions in the reconstruction process. This may be due to the low modulating efficiency for converting the high-frequency information to the low-frequency band, especially when we aim at one order of magnitude resolution gain. Third, we assume the PSF is shift-invariant in the recovery process. If needed, we can divide the captured images into many smaller segments, and jointly recover the object, the illumination pattern and the PSF for each segment. This is similar to the pupil recovery process in Fourier ptychographic microscopy, where the pupil aberrations are recovered at different spatial locations independently [30,31].
There are also a few directions for improving the reported scheme. First, we can use a scanning galvo mirror to better control the tilted beam and achieve better positional accuracy. Second, we can employ motion correction [32] in the reconstruction process to address the positional error of the scanning process. Third, we may incorporate support constraints [22] including signal sparsity [33,34] in the reconstruction to reduce the number of acquisitions. Fourth, multi-layer modeling [35] can be integrated into the reported scheme to handle 3D fluorescence sample. In step 4, we generate the 'mother speckle pattern' corresponding to P unknown (x,y) in the Main Text with which we obtain shifted speckle patterns in the next step. The 'mother speckle pattern' is a random intensity pattern of size 165 by 165. It is 30 pixels ('MaxShiftInPixel') larger than the object image in each direction, because we assume the maximum shifting distance in each direction is 30 (= 0.75 * (41-1)) pixels. In step 5, we generate a sequence of shifted speckle patterns whose size is 135 by 135. Each shifted speckle pattern is generated by laterally shifting 'mother speckle pattern' and cropping the central 135-by-135-pixel area. The 'mother speckle pattern' is laterally shifted along a spiral path (which is represented by a 41-by-41 matrix, as shown in Fig. 12(a)). The spiral path matrix can be easily obtained using the built-in function 'spiral()' by MATLAB. As shown in Fig. 12(a), the spiral path matrix consists of 1681 integers from 1 to 1681. To generate the n th shifted speckle pattern, for instance, we will search for the location of 'n' in the matrix. Assuming the number of 'n' is located at (LocationX, LocationY), the lateral shift will be 'LocationX*ShiftStepSize' along the x direction and 'LocationY*ShiftStepSize' along the y direction. The lateral shifting of the pattern is performed in the Fourier domain by multiplying an equivalent phase factor to the spectrum of the 'mother speckle pattern' (we use the function 'subpixelshift()'). We get the final n th shifted speckle pattern by cropping the central 135-by-135-pixel area from the shifted 'mother speckle pattern'. The number of shifted speckle patterns is 1681 (= 41 by 41) and we save them in 'ShiftedPatterns', which has dimension 135 by 135 by 1681. The relative lateral shifts along both directions are stored In step 6, we generate the low-resolution images captured through the diffuser. We simulate these images by multiplying the true object image with a shifted speckle pattern and then applying the low-pass filter given by the OTF of the diffuser. The resultant lowresolution images are stored in the matrix 'TargetImgs' which is organized as an image stack with dimension 135 by 135 by 1681. Figure 13(a) shows the 9th generated low-resolution target image. %% Step 6: Low-resolution target image generation 82. TargetImgs = zeros(size(ShiftedSpecklePatterns)); 83. nTargetImgs = nPattern;% Total number of target low-resolution images 84. for iTargetImg = 1:nTargetImgs 85.
TargetImgs(:,:,iTargetImg) = abs(ifft2(ifftshift(fftshift(fft2(InputImg .*ShiftedSpecklePatterns(:,:,iTargetImg))).*OTF))); 86. disp(iTargetImg); end; In step 7, we apply the iterative algorithm described in the Main Text to reconstruct both the high-resolution object image, the OTF, and the 'mother speckle pattern'. In our iterations, ImgRecovered, OTF and MotherSpeckleRecovered correspond to Object, PSF and P unknown in the Main Text. The number of iterations are defined by variable 'nIterative'. Assuming we have no any prior about the object image or the speckle pattern, we initialize the object by averaging all measurements and initialize the 'mother speckle pattern' using an all-ones matrix. Figure 13(b) shows the reconstructed high-resolution object image and 'mother speckle pattern' after 30  LRTempTargetImg_AmpUpdated_FT = fftshift(fft2(LRTempTargetImg_AmpUpdated)); 101.
% Update the target image 103.
MotherSpeckleRecovered pause(0.5) 120. end Fig. 13. Simulation result of the resolution improvement. (a) The 9th generated low-resolution target image. (b) The reconstructed high-resolution object image and 'mother speckle pattern'.

Funding
NSF 1555986, and NIH R21EB022378, NIH R03EB022144. Z. Zhang acknowledges the support of the China Scholarship Council.

Disclosures
The authors declare that there are no conflicts of interest related to this article.