High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering

Thick biological tissues give rise to not only the multiple scattering of incoming light waves, but also the aberrations of remaining signal waves. The challenge for existing optical microscopy methods to overcome both problems simultaneously has limited sub-micron spatial resolution imaging to shallow depths. Here we present an optical coherence imaging method that can identify aberrations of waves incident to and reflected from the samples separately, and eliminate such aberrations even in the presence of multiple light scattering. The proposed method records the time-gated complex-field maps of backscattered waves over various illumination channels, and performs a closed-loop optimization of signal waves for both forward and phase-conjugation processes. We demonstrated the enhancement of the Strehl ratio by more than 500 times, an order of magnitude or more improvement over conventional adaptive optics, and achieved a spatial resolution of 600 nm up to an imaging depth of seven scattering mean free paths.


Supplementary Note 1. Theoretical framework for the CLASS microscopy
As discussed in the main text, the maximization of the total intensity of CLASS image leads to the finding of angle-dependent aberrations. In this section, we provide detailed mathematical description of its working principle and conditions for the convergence.

Mathematical formulation of the aberration correction for the input
Let us begin with the first round of iteration for the aberration correction in the illumination beam path. After applying for initial arbitrary angle-dependent phase corrections ( ⃗ ), the spectrum of the CLASS image is written as For each ⃗ , we identified ( ⃗ ) that maximizes the total intensity of the corrected CLASS image: The maximum of the total intensity occurs when the sum of the cross-terms between singlescattered waves in Eq. (2) are the largest. This is the case when the following cross-term originating from the two arbitrary incident wavevectors, ⃗ and ⃗ ( ⃗ ), are real valued. The second condition that determines the convergence of iteration is the relative intensity of multiple light scattering with respect to that of single scattering. In assessing the cross-term of the single-scattered waves shown in Eq. (3), the cross-term of the multiple-scattered waves described in the second term in Eq. (1) serve as noise. If is too large, then the phase of this crossterm will not be finite in width, but uniformly distributed between -π and π. See the detailed numerical test of the convergence at Supplementary Note 2.

Closed-loop iteration by the phase-conjugation operation
After identifying ( ⃗ ) that maximizes the total intensity of the CLASS image (Eq. (4) Therefore, finding ( ⃗ ) that maximizes the total intensity of CLASS image is identical to replacing input aberration ( ⃗ ) with Φ ( ⃗ ). Under the convergence conditions discussed in 1.2, the uniform random phase distribution of initial aberrations is reduced to the finite width of the distribution of Φ ( ⃗ ). However, there will not be a further reduction without additional independent operation to correct the sample-induced aberrations in the collection process. In order to form a closed-loop correction, we considered a phase conjugation process in which the wave is incident from ⃗ and reflected to ⃗ ( ⃗ Δ ⃗ ). With this reverse process, the CLASS spectrum after the initial correction in Eq. (5) can be written as And this output correction replaces the output pupil function with This is equivalent to replacing output aberration ( ⃗ ) with Φ ( ⃗ ) . Note that the Φ ( ⃗ ) is even more narrowly peaked around zero than Φ ( ⃗ ) because the input correction ( ⃗ ) has converted a to a more slowly varying function . As a consequence, the error of the first round of output aberration correction is smaller than that of the illumination correction.
With the first round of input and output corrections in place, the second round of input and output corrections, ( ⃗ ) and ( ⃗ ) respectively, are applied to progressively reduce the residual error of the aberration corrections. And the iteration keeps going until the residual error is made below certain level, which is typically set as 0. 1 radians. We can find the recursion relation for both input and output pupil functions after n iterations.
with initial condition 0 a , and 0 a .
As the iteration number is increased, the phase angles Φ ( ⃗ ) and Φ ( ⃗ ) converge to zero on the conditions described in 1.2. Consequently, the pupil functions for the input and output in Eqs. (11) and (12) converge to the ideal pupil function ( ⃗ ).
After all, we can find the initial aberrations by accumulating phase corrections identified in all the iterations: Note that the complex conjugation operation is not critical from the mathematical point of view in the second step of CLASS algorithm. It doesn't make any difference in the result except for the sign of the output phase correction map. What is important is the transpose operation, which enables us to swap input and output, and to correct the specimen-induced aberration from the output side. The reason we take phase conjugation (transpose + complex conjugation) rather than just the transpose is to impose a physical meaning to the operation. In this way, the second step is equivalent to physically sending the waves from the output side and measuring the backscattered waves at the input side.

Supplementary Note 2. Numerical simulation study
For the complete understanding of the effect of CLASS algorithm on correcting the aberrations of single-scattered waves and attenuating the contribution of multiple light scattering, we performed a numerical simulation study by assuming a representative example of sample-induced aberrations.

Applying CLASS algorithm to the numerically prepared aberrations
We set up a numerical simulation at the condition of 0 and m 1, 4 , which corresponds to the number of free modes for a 20×20 m 2 field of view. Also, we introduced arbitrary aberrations ( ⃗ ) and ( ⃗ ) as depicted in Supplementary Figures 1a and 1b, respectively, which led to 1 400 and 1 600. The amplitude of the cross-correlation map of these two complex pupil functions (Supplementary Figure 1c ) was well below unity, suggesting that the accumulation of single scattering would be compromised. Supplementary   Figures 1d and 1e show CASS images without and with aberrations, respectively, in the absence of multiple scattering. The target objects were a pair of point particles separated by 600 nm, which corresponds to the diffraction-limit resolution for 0.8 NA at the source wavelength  = 800 nm. As expected, the aberrations made the two particles completely indistinguishable.  Figure 2b shows the CLASS image reconstructed after applying this phase correction. The existence of particles becomes better visualized than before. However, the resolving power has not been sufficiently recovered to distinguish the two particles because the output aberration has not yet been addressed.

Supplementary
This first round of maximization operation is incomplete because only the aberration arising from the incident wave can be dealt with. With the correction in place, we apply the phase correction to the output and identified the that would maximize the total intensity of the phase-conjugated CLASS image. Similar to the correction for the illumination path, this iteration leads to the convergence of ( ⃗ ) to ( ⃗ ). In fact, this correction for the reflection process converges faster than that of to ( ⃗ ). Because of the prior correction , the width of the phase histogram of Φ ( ⃗ ) is narrower than that of Φ ( ⃗ ) (see

Reduction of phase correction error with the increase of iteration number
In this section, we present the numerical simulation data that visualizes the way the iterative maximization operations (Eqs. (2) and (8) In order to quantitatively assess the effectiveness of convergence, we computed the correlation between the original aberrations and the aberration correction maps identified from CLASS algorithm. To this end, the correlation between the obtained aberration maps and the initial aberration maps in Supplementary Figures 1a and 1b was computed for each iteration number. As shown in Supplementary Figure 4, only 3 iterations were good enough to recover the initial aberration map with the accuracy of more than 97%. Note that the output correlation for the n th iteration step is always larger than n th input correlation because the n th output correction was performed after n th input correction.

The effect of aberration corrections to the transfer function of CLASS image
By Eqs. (1), (6), and (10), the spectrum of single-scattered waves in CLASS image after times of input-output phase corrections can be written as The amplitude transfer function (ATF) of CLASS image is modified to ⋆ after n th round of iteration process. The width and height of ATF are related to the spatial resolution and signal strength of single-scattered waves, respectively. In Supplementary Figure 5, we present for various n for the numerical simulation data. As expected, converged to the ideal ATF without aberration as n was increased. The width of ATF was broadened to that of the ideal ATF, indicating that the diffraction-limited spatial resolution was almost recovered.
More importantly, the average height of ATF was increased by about 20 times after the correction, which corresponds to the increase of the intensity of single-scattered waves by about 400 times. This is critical in recovering single-scattering signals when there exists strong multiple-scattering background. Indeed, our correction method not only recovers spatial resolution, but also enhances the signal strength of single-scattered waves.

The intensity of single-and multiple-scattered waves with the increase of iteration number
Since our aberration correction method is based on enhancing the coherent summation of singlescattered waves, the total intensity of CLASS image should increase by the iteration process. In the simulation data, we can separate out single-and multiple-scattered waves in the CLASS image. Therefore, we could analyze the way the correction process affects to the intensities of single-and multiple-scattered waves. Supplementary Figure 6 shows the intensities of the singleand multiple-scattered waves as well as the total intensity of the CLASS image as a function of the iteration number. Initially, multiple-scattered waves was about 20 times larger than the single scattered waves because of the strong aberration. Total intensity of CLASS image was almost equal to that of multiple-scattered waves, and this is the reason the target was initially invisible.
When the iteration number was increased, the intensity of single-scattered waves increased significantly by up to about 400 times while that of the multiple-scattered waves didn't change much. As a result, the intensity ratio between single-and multiple-scattered waves was reversed, respectively. Plots were normalized by the initial contribution of single-scattered waves.

Relation between and Strehl ratio S
In adaptive optics, the degree of aberration is measured by Strehl ratio, which is the peak intensity of the point spread function with aberration relative to that without aberration. For the case when both the input and output aberration matter, this parameter can be described by the amplitude transfer function, Here a ⋆ a , and 0 ( ⃗ ) is the ideal amplitude transfer function. For a relatively mild aberration and weak multiple scattering noise, faithfully describes the effect of aberrations to the signal strength and image contrast. But it's use is limited in case of severe aberration or strong multiple scattering noise. In our study, we introduced the parameter in Eq. (3), which is the ratio of the total intensity of single-scattered waves in the presence of aberration with respect to that in the absence of aberration. In terms of the amplitude transfer function, can be simplified to 〈| ( ⃗ )| 〉 〈| 0 ( ⃗ )| 〉 16 In other words, is the ratio of the total intensity of the point spread function with aberration to that without aberration. As described in the main text, is responsible for the degradation of the signal to noise ratio of CASS microscopy from ⁄ m to × ⁄ m . And the main aim of CLASS microscopy is to increase enough to raise the signal to noise ratio above unity. And the initial will set the limit that the CLASS microscopy can work.
The and are related to a certain extent. The increase of aberration accompanies the simultaneous reduction in the peak height and the total intensity of the point-spread-function. In general, the peak height is more susceptible than the total intensity to the aberration. Therefore, is usually smaller than . But their relation depends on the type of aberration, making it difficult to find a general relation between them. For a heuristic purpose, here we gave an example for the aberration maps assumed in Supplementary  On the other hand, describes the reduction of single-scattering intensity in the CLASS microscopy where images are acquired over wide area. Unlike confocal microscopy where signal is collected only at the confocal pinhole, the entire PSF contributes to the signal intensity in the wide-field detection. Therefore, the degradation of SNR in the CLASS microscope is determined by , not by S. Likewise, the main role of CLASS algorithm is to increase single scattering intensity by 1 , not by 1/S.

The effect of multiple light scattering to the convergence of CLASS algorithm
As we discussed in section 1, the second condition that determines the convergence of iteration is the relative intensity of multiple light scattering with respect to that of single scattering. In assessing the cross-term of the single-scattered waves shown in Eq.
(3), the cross-term of the multiple-scattered waves described in the second term in Eq. (1) serve as noise. If is too large, then the phase of this cross-term will not be finite in width, but uniformly distributed between -π and π. This can be seen in Supplementary Figures 8a-d in which the width of distribution was gradually increased as was increased. But the convergence of the iteration was guaranteed even when ~16,000.

Supplementary Figure 8. Effect of multiple scattering on the phase correction error, Φ ( ⃗ ). The histogram of Φ ( ⃗ ) when the relative intensity of multiple-to single-scattered waves increases. a-d,
Intensity ratio between multiple-to single-scattered waves, are 0, 20, 40, and 45, respectively and 1 400 ⁄ . Therefore, the initial intensity ratio of the multiple-to single-scattered waves, , are 280, 8000, 16,000, and 17,550, respectively. Up to 16,000, the histogram has a finite width such that we could successfully characterize the input and output aberration. Detailed experimental setup of CLASS microscopy is shown in Supplementary Figure 9. The basic layout of the setup is identical to the CASS microscopy that we previously reported 1 .

Supplementary Note 3. Detailed experimental setup and additional analysis of experimental data
However, light source, objective lens and camera were upgraded, and random phase patterns were used for the illumination. The backbone of the setup was a Mach-Zehnder interferometer, and low coherence light source and off-axis geometry were used for the time-resolved and wide-field complex-field imaging.
A femtosecond laser pulse from mode-locked Ti-Sapphire laser (center wavelength 800 nm, bandwidth 30 nm, and pulse width 100 fs) was divided into incidence and reference waves at BS1.
The transmitted light from BS1 was delivered to the spatial light modulator (SLM, Hamamatsu X10468-02). By writing random phase patterns on the SLM, we generated random speckle fields and use them as incidence beams to the sample (inset on the left of the SLM). By two 4-f imaging systems formed by L3, L4, L5, and OL1 (objective lens, Nikon CFI-Apo-40X-W-NIR, 0.8 NA water dipping), we delivered the image of the speckled incidence beam into the sample plane. As shown in the upper right of the figure, we used water as an index matching medium to minimize the refractive index mismatch between the medium and the biological tissue. The reflection from the sample was collected by the same objective lens, and then delivered to the camera (pco.edge 4.2 scientific CMOS camera) placed on the conjugate image plane of the sample plane.
The reference wave reflected at BS1 traveled through similar optics as the sample wave. A scanning mirror was placed in the reference beam path to tune the path length of reference wave.
To generate an off-axis interferogram, we selected the 1 st order diffracted beam from a diffraction grating (Ronchi ruling, 72 lp/mm, Edmund optics) placed on the conjugate image plane of the camera. At BS6, the reference wave was combined with the reflected wave from the sample. Inset on the right of the camera in Supplementary Figure 9 shows the typical image of the interference pattern. Stripe patterns appearing on the speckled pattern of the sample wave were due to this offaxis interference. By applying Hilbert transform on the measured interferogram, we could obtain the time-gated complex field map of the reflected wave from the sample.
To construct a time-gated reflection matrix, typically 2,800 complex images were taken for the pre-determined set of random patterns of illumination written on the SLM. The number of images was determined in such a way to cover the maximum number of orthogonal channels for 0 × 0 view field with the collection angle of 0.8 NA.

Measurement of CLASS image using random basis.
As discussed in the online Methods, we recorded a time-resolved complex field map ⃗ , 0 for the j th random pattern written on the SLM by placing an ideal mirror at the sample plane for The use of random basis enabled us to distinguish the uncontrolled phase shifts of interferometric detection induced by the fluctuation of the relative beam path between sample and reference waves from the angle-dependent aberrations by the sample itself. In this case, we applied for an additional phase correction step to deal with the uncontrolled phase shifts among measurements for different speckled illuminations. To this end, we constructed the time-gated reflection matrix on the momentum difference domain, ( ⃗ ⃗ ⃗ ) as discussed in the method section. By multiplying the input basis matrix, ( ⃗ ), the input plane wave basis was converted into the random speckled illumination basis, The individual columns in ( ⃗ ⃗ ) correspond to the object spectrums in the momentum difference for the illumination of j th speckled pattern. Therefore, the summation of the matrix along row direction (i.e. summation of elements with the same Δ ⃗⃗ ⃗⃗ ⃗⃗ ) will lead to the coherent accumulation of single-scattered waves on the basis of speckled illumination. In order to deal with uncontrolled phase shifts, we added additional phase correction factor, , for each illumination pattern. And the set of 's that maximize the total intensity of the coherent summation will correspond to the uncontrolled phase shifts.
For example, we present the obtained the uncontrolled phase shifts for the measurement in Supplementary Figure 11. Once this process is done, we multiply to ( ⃗ ⃗ ) to which ( ⃗ ) was multiplied to the time-gated reflection matrix ( ⃗ ⃗ ) free of uncontrolled phase shifts. The aberration correction method discussed in the main manuscript was performed afterwards.
Supplementary Figure 11. The uncontrolled phase shifts for the experimental measurement in Fig. 2. Fig. 2 For the experimental data in Fig. 2, CLASS images are shown in Supplementary Figure 12a for different iteration number. The image quality was greatly enhanced only after the 1 st round of the iteration. As the number of iterations was further increased, the pattern boundary was sharpened.

Detailed analysis of iteration process for the experimental data in
In addition, the intensity of CLASS image was increased up to 25 times compared with the initial image. The plot in Supplementary Figure 12b shows the quantified increase of total intensity of CLASS image. Due to the proper accumulation of single-scattered waves by means of aberration correction, the signal to background intensity ratio was also enhanced (Supplementary Figure 12c) by about 25 times. The line profiles of images before and after the correction shown in Supplementary Figure 12d confirmed that line spacing of 600 nm, which corresponds to the diffraction limit of the system, was clearly resolved and the contrast of line structures was greatly enhanced. Fig.3. a, CLASS image after n th iteration (n=0, 1, 2, and 10). n=0 corresponds to the initial CLASS image. Color scales indicate normalized intensity by the maximum intensity for n=0. Scale bar, 4 μm. b, Total CLASS intensity integrated over the entire view field as a function of the number of iterations. Data was normalized by the initial intensity (n=0). c, Signal to background intensity ratio. Signal intensity was the average intensity at the patterned area, and background intensity was the average intensity for the rest. d, Line profiles of images before and after aberration correction (dashed lines in a for n=0, and n=10). For better comparison, the line profile for n=0 case (blue line) is multiplied by a factor of 10.

Transfer function analysis of experimental data
For the experimental data in Fig. 3, we could perform similar transfer function analysis performed for the simulation data in section 2. If we assume that the initial aberration maps of the target were identical to the obtained aberration maps in Fig. 3i-j, we can estimate the initial amplitude transfer function of single-scattered waves in CASS image by Eq. (14). As shown in Supplementary Figure 13a, the initial ATF is not uniform, and has less than 10% of amplitude compared to the aberration-free case. Also, the effective bandwidth at k-space is reduced about half of ideal case, which indicates the reduction in resolving power. By the application of CLASS algorithm, the ATF converged to the aberration-free case (Supplementary Figure 13b). Indeed, both the amplitude and effective bandwidth of ATF have increased by the iteration in accordance with the theoretical expectation. Fig. 3. a, Initial amplitude transfer function obtained by the cross-correlation of input and output aberration maps in Figs. 3i-j. Color scale is normalized by the maximum value of ideal transfer function. Scale bar, 0 . b, Line profiles of the amplitude transfer function at 0 for n=0, 1, 2, 3, and 10.

Comparison with Zernike-based adaptive optics
The major advancement that CLASS microscopy has made is the ability to characterize and eliminate steep variation of angle-dependent phase shifts induced by a thick scattering layer. The conventional adaptive optics experiments usually concern a few low-order Zernike modes, thereby confining their applicability to slowly varying angle-dependent phase shifts caused by mild aberrations of a thin scattering layer. To distinctively visualize the difference between CLASS microscopy and adaptive optics, we applied only the first 15 Zernike modes identified from the input and output aberration maps identified in Figs. 3i and 3j. As shown in Supplementary Figures 14a and 14b, fine details of aberration maps in Figs. 3i and 3j were missing for Zernike modes-based reconstruction. In addition, the coherent summation of singlescattered waves is not as effective as the result shown in Fig. 3l as indicated by the color scale.
This suggests that even 15 orders of Zernike modes were not enough to deal with the aberration of the scattering layer that we used. Indeed, CLASS microscopy can address steep variation of specimen-induced aberrations that conventional adaptive optics cannot handle.

The estimation of in Fig. 3
There were about 20 times increase in signal intensity at the target in Fig. 3l in comparison with the initial intensity in Fig. 3d. But this doesn't mean that η is 1/20. In fact, the real η that we estimated by applying Eq. (16) to the measured aberration maps shown in Figs. 3i and 3j was 1 04. This discrepancy arises mainly because multiple-scattered waves as well as singlescattered waves contributed to the signal intensity at the target in Fig. 3d. Before the application of CLASS algorithm, the intensity of CLASS image at the target is the sum of single-scattered waves and the multiple-scattered waves that survived the time-gating and spatial coherence gating, By inserting the estimated from the measured aberration maps, we can obtain the initial singleto-multiple scattering ratio in the position basis, m 19 1 4 ⁄ 0 1 0 Therefore, multiple scattering was about 10 times stronger than the single scattering before the aberration correction, and this was responsible for the discrepancy mentioned above. Since the initial SNR of m 0 1 was far less than unity, multiple scattering has strongly degraded the resolving power of imaging. By the aberration correction, this SNR was increased to m 1 1 by selectively enhancing the intensity of single-scattered waves.

The difference between input and output aberration maps
In the ideal case, the aberrations for the input and output should be identical in the reflection geometry since waves travel back and forth through physically the same sample. In the real practice, however, the optical system in the illumination beam path from the light source to the beam splitter and that in the collection beam path from the beam splitter to the camera are not perfectly the same. Therefore, input and output paths have different system aberrations. In addition, the slight misalignment of the optical axis or image focus can make the difference even larger. In general, the specimen-induced aberrations are much stronger than these system aberrations such that aberration maps for the input and output looked largely the same. For example, the normalized cross-correlation of input and output pupil functions given by the aberration maps in Figs. 3i and 3j, respectively, was 0.64. But when we looked at the phase difference between the two (Supplementary Figure 15), their difference was clearly visible. We could observe slowly varying phase difference across the pupil due to the difference of the system aberrations between the input and the output.
Supplementary Figure 15. Phase difference between input and output aberrations obtained in Fig. 3.
Note that we also tried to correct aberration with the assumption that the input and output aberrations are identical. Supplementary Figure 16  The measurements of the scattering mean free paths of scattering media In order to determine the scattering mean free path of the scattering medium used in our study, we measured the intensity of ballistic photons as a function of the thickness of the scattering layer.
We placed a stack of phantom scattering layers of known thickness on a flat mirror, and then measured the decay of the intensity of ballistic components depending on the number of phantom layers. Typical thickness of individual layer was 160 µm. Specifically, we illuminated a normally incident plane wave and recorded the complex field map of the backscattered waves. We then extracted the intensity of the normally reflecting plane wave component from the recorded map.
By using the exponential curve fitting, we could obtain the decay constant of 51.2 µm for the roundtrip (Supplementary Figure 17). Therefore, the scattering mean free path of the scattering layer was 102.4 µm.
Supplementary Figure 17. Measurement of the scattering mean free path of the scattering layers used in Fig. 1.
Here one can notice that the first four points up to the thickness of 500 µm (1 mm in travel distance including the returning path, which corresponds to about 10 l s ) fit well to the exponential curve while the intensity taken at thicker depth was slightly larger than the fitted curve. In our measurement of the ballistic photons, we used both angular-and time-gating by the time-resolved complex field imaging. Like any other ballistic photon measurement methods, our method is effective up to a certain thickness of the sample. If the scattering medium becomes too thick, some of the multiple-scattered waves can pass through these gating operations, and this is source of discrepancy that we observed at the thicknesses larger than 500 µm. In fact, the conventional method mostly uses the angular gating in the transmission geometry (for instance, see Ref. µm in our case, will ensure the accurate estimation of the scattering mean free path of individual 160 µm-thick scattering layers. And this estimation will reasonably be applied to the multiple stacked layers.
In the case of rat brain tissues, scattering properties of the tissue slices vary with imaging position. Therefore, we estimated the scattering mean free path in situ from the total intensity of the CLASS image, which is close to the intensity of the ballistic photons, and the thickness of the tissue layer measured from the time delay. We took logarithm of this total intensity of the ballistic photons normalized by that measured with the mirror sample for the calculation of the scattering mean free path.

Supplementary Note 4. Additional experiments
In the main text, we used a resolution target having spatial variation of reflectance for demonstrating the CLASS microscopy. In this section, we provided experimental data for other types of targets such as a resolution target with phase contrast and point particles.

CLASS imaging of phase objects
In this section, we present experimental study of applying CLASS microscopy for phase-contrast targets. As a test target, we fabricated a gold layer coated on a slide glass in such a way that the thickness of the coating varies in space (Supplementary Figure 18a). In particular, the part of gold layer exhibiting the pattern of the resolution target was thicker than the rest by about 100 nm.
Since the gold layer were coated over entire area, reflectance was almost uniform across the field of view. On the other hand, the reflected waves experience less phase retardation at the pattern than that at the rest. The height difference of about 100 nm corresponds to the phase retardation of 2 radians. Therefore, this target acts as a phase-contrast object.
Similar with Fig. 1b, we put both the aberrating layer and 5l s -thick scattering layer on the target as shown in Supplementary Figure 18a. Supplementary Figures 18b and 18c show the amplitude and phase maps, respectively, of CLASS image before aberration correction in which targets were invis ible due to the aberration and multiple scattering. After applying for the CLASS algorithm, we could clearly identify target structures both in amplitude and phase maps ( Supplementary   Figures 18d and 18e). In the amplitude image ( Supplementary Figure 18d)

CLASS imaging of point particles
In this section, we demonstrated the CLASS microscopy for point-like objects. Gold nanoparticles with the diameter of 400 nm were placed underneath the same set of aberrating and scattering layers used in section 4. In this particular experiment, the cylindrical groove in the aberrating layer was rotated by about 60 degrees with respect to horizontal axis. As shown in Supplementary Figure 19a, the image of particles were elongated along the orthogonal direction