1. Introduction

Imaging with second harmonic generated (SHG) light enables label free imaging of non-linear structures. This intrinsic contrast mechanism, which relies on the lack of inversion symmetry, allows selective imaging of particular features, while eliminating background. Leveraging this advantage, SHG microscopy is continuously growing as a valuable resource for the study of biomedical and material systems [1–3]. In biological tissues, light undergoes second harmonic scattering when interacting with non-centrosymmetric molecules that are ordered spatially so that coherent nonlinear second harmonic scattering from the tissues add constructively to produce a measurable SHG signal [4–8]. SHG has proven to be a valuable method for identifying a wide range of diseases [9,10], including to quantify the alignment of collagen surrounding tumors to grade metastatic potential [11]. SHG microscopy has even been used for mapping cell lineage in embryos by tracking cell division using SHG generated by the mitotic spindle during mitosis [12]. SHG microscopy has found significant use in materials science [13] and investigating two-dimensional materials [14].

Standard SHG imaging is based on laser scanning microscopy, in which an incident laser beam at the fundamental wavelength is focused tightly into a sample. A portion of the SHG power is collected in either the forward- or backward-scattered direction at each focal point [15]. An SHG image is built from assigning the measured power to a location in a matrix corresponding the spatial location of the focused fundamental beam. Unfortunately, this leads to slow image formation, since each point in the image must be collected sequentially. The signal to noise ratio (SNR) also suffers because the signal is collected from each spatial point in an image only for the time that the laser beam dwells on each focal point. The SHG signal power is proportional to $|\chi ^{(2)}|^2$, where $\chi ^{(2)}$ is the nonlinear susceptibility responsible for SHG signal generation. Conventional SHG microscopy does not directly reveal the desired spatial map of $\chi ^{(2)}$, with only the magnitude of the susceptibility that depends both on the spatial distribution and sign of the susceptibility distribution within the focal volume of the focused fundamental beam. Complex image information, notably the sign of $\chi ^{(2)}$, which indicates the orientation of the SHG-active molecules, can be obtained by interferometric single-pixel detection SHG imaging [16,17]. However, the lack of a stable reference phase from a repeated set of measurements prevents an improvement in the image SNR that would be possible with averaging the image fields, rather than the image intensity [18].

While such conventional nonlinear laser scanning microscopy benefits from the non-linear spatial filtering that helps with forming three-dimensional images and imaging within scattering media, optical aberrations degrade this imaging method. The SNR, image quality, and spatial resolution of SHG imaging are affected by these optical aberrations introduced by the imaging system itself and from specimen variations in the refractive index [19–21]. In SHG microscopy, the distortions introduced by the optics, particularly the objective lens, and the specimen broaden the size of the focused beam, worsening the ability of the microscope to image fine spatial features and reducing the signal level. Adaptive optics methods [20] have been applied to improve imaging with point-scanning nonlinear microscopy [22,23], including wavefront shaping for polarization-resolved SHG imaging within tissues [24].

Imaging speed and SNR are significantly improved with widefield SHG holographic imaging [25–31]. Speed is increased with SHG holography for two reasons. The first is that widefield images are recorded on a camera, so that each pixel benefits from signal being recorded for the entire imaging time. Thus, even for faster imaging, the SNR of the image is improved. Secondly, the hologram is formed from the interference of a signal and a reference beam, producing a heterodyne signal amplification that allows for optimization of the SHG imaging speed [30]. This amplification allows even very weak SHG signal fields to be detected at the shot noise limit. Additionally, holography allows for extraction of the complex field, so that amplitude and phase information is available, and the nonlinear susceptibility can be extracted by solving the inverse scattering problem [31]. Widefield SHG imaging speed can be further increased using a laser with very high pulse energy (lower repetition rate) enabling the capture of an image per pulse [32].

Widefield SHG imaging has been restricted to a trans-illumination geometry because generally SHG fields that are scattered in the forward direction are much stronger than the backward direction in biological tissues. Point scanning images that are collected in the backscattered direction consist of a combination of directly backscattered SHG radiation [1,33] and forward-scattered SHG light that is re-directed in the backward direction so that it can be collected in a epi configured microscope [34]. The ratio of forward and backward scattered SHG power of ex-vivo tissues has proven useful as a biomarker for distinguishing healthy and cancerous tissues [9,35]. While conventional laser scanning SHG microscopy can be deployed favorably in biological tissues that highly scatter fundamental and SHG light, widefield SHG imaging has been degraded by optical scattering, which is dominated by randomization of the phase of the SHG field [30].

Measuring widefield SHG holographic images in a transmission and epi configuration would be extremely valuable for imaging collagen and muscle in tissues in a minimally invasive manner. While point scanning SHG imaging can be performed in an epi direction, such a conventional approach suffers from very weak signals [33,34], limiting practical use. Holographic widefield SHG in a epi configuration will enable improved detection of weak backscattered signals as a result of heterodyne amplification. Furthermore, imaging in a backscattered configuration would allow for direct optically-sectioned imaging because the low-coherence interferometry will gate only backscattered SHG light over an axial depth of the coherence length of the SHG light – exactly analogous to depth sectioning achieved with optical coherence tomography.

In this Article, we demonstrate the first epi collected widefield SHG images leveraging the heterodyne signal enhancement provided by holographic measurements to mitigate the weak backscattered signal strength. Additionally, we exploit phase information to coherently superimpose measured fields obtained from a set of illumination angles to implement synthetic aperture coherent nonlinear holographic imaging for second harmonic generation (SHG) scattering from samples. In synthetic aperture holography, [36–38] complex spatial frequency information from multiple field measurements is combined to produce a net complex field image with spatial frequency support that is expanded up to a factor of two, improving imaging resolution. Aberrations, represented as a phase variation across the pupil, can severely distort the synthetic aperture image [37].

We introduce a robust and computationally efficient algorithm to estimate and correct the pupil phase distortions responsible for aberrations in SHG imaging. The acquired data contain sufficient redundancy to allow estimation of the imaging system aberrations directly from the recorded data. Redundancy in the field was used to identify conserved coherent field amplitudes to selectively suppress noise in the estimated image. When phase corrections are applied, we observe drastic improvements in SNR and image quality of the SHG images. Utilizing the linear properties of wave propagation and synthetic time reversal, the pupil phase distortions of both the input and imaging pupil planes can be compensated, thereby correcting system as well as sample induced aberrations. The result is a diffraction limited SHG image with a spatial frequency support twice that present in a single holographic SHG image, or four times the spatial frequency support of the fundamental field. Finally, we demonstrate synthetic aperture SHG holography on transmitted SHG fields in addition to the first back scattered SHG fields collected in the epi direction of the SHG holographic microscope.

2. Theory

The experiments described here involve imaging a thin SHG-active sample when illuminated with a fundamental plane wave. The SHG scattered fields are captured in both the transmitted and epi directions as the input plane wave propagation direction is varied across the aperture of the condenser lens. Referring to Fig. 1, we see that microscope consists of a pair of matched objective lenses. The illumination for both the epi and transmission configurations thus pass through the same condenser objective lens with a pupil phase $\phi _{\rm 1}(\mathbf {x}_i)$ at the fundamental beam wavelength $\lambda _1$. Plane wave illumination means that the fundamental beam passes through a small point in the input pupil plane located at $\mathbf {x}_i$, which maps to an input spatial frequency of $\mathbf {u}_i = (\lambda \, f_c)^{-1} \, \mathbf {x}_i$ with wavenumber $\lVert {\mathbf {u}_i}\rVert=1/\lambda$ and $f_c$ denoting the condenser lens focal length. As SHG scattering is driven by the square of the fundamental illumination beam, the effective input pupil phase is $\phi _i=2\, \phi _1$. These input aberrations are transmitted to the scattered field and distort the image. In the case of synthetic aperture holography, these distortions are replicated across the image field spatial frequency distribution, as illustrated in Fig. 1(d).

The propagation of light being linear, we can describe the relationship of a given light field from one plane to another with a simple matrix operation (reflection or transmission matrix depending on the configuration). The choice of input and output planes, and thus the basis of this matrix, is chosen to be the input and output pupil planes, $P_i$ and $P_o$ shown in Fig. 1. A given input angle conveniently corresponds to a point, $\mathbf {u}_i$, in the input pupil plane. At the output pupil plane, the input plane wave is scattered by the sample into many angles, each given by a point in the output pupil, $\mathbf {u}_o$. This scattered field is proportional to the spatial frequency map of the second order susceptibility $\hat {\chi }^{(2)}({{\mathbf{q}}})$ of the sample, where the object spatial frequency, ${{\mathbf{q}}}$, will also be used to denote the scattering vector.

The imaged SHG field can be described in the output spatial frequency plane with coordinates $\mathbf {u}_o$. By invoking the assumption of a thin specimen and assuming that the fundamental field is not depleted appreciably in the nonlinear scattering process, we may write the scattered field in the output pupil plane as

The thin specimen is described by a two-dimensional second order nonlinear susceptibility distribution, $\chi ^{(2)}(\mathbf{r})$, that lies in the sample plane with coordinates $\mathbf {r}$. Light is scattered at the second harmonic frequency of the incident fundamental beam at frequency $\omega _1$, with a Green’s function, $G(\mathbf{r},\mathbf{u}_i)$, describing the square of the fundamental field incident on the sample. This function maps input spatial frequencies for each point $\mathbf {u}_i$ to the SHG driving term at the sample plane. The scattered field is collected by the objective and mapped from the sample plane $\mathbf {r}$ to the output imaging pupil $\mathbf {u}_o$ with the Green’s function, $H(\mathbf{u}_o,\mathbf{r})$, for the SHG field at optical frequency $\omega _2 = 2 \omega _1$. This Green’s function can be used to describe imaging of the forward-scattered field in a trans-SHG holographic microscope or to image the back-scattered field in an epi-SHG holographic microscope.

Within an isoplanatic spatial imaging region, the imaging point spread function is spatially invariant, which allows the transfer function to be modelled with the pupil function, $P({{\mathbf{u}}}) = |P({{\mathbf{u}}})| \, \exp \left [i \, \phi ({{\mathbf{u}}}) \right ]$, where the spatial frequency support is $|P({{\mathbf{u}}})|$ and $\phi ({{\mathbf{u}}})$ accounts for aberrations. In addition to aberrations, there are random phase shifts due to air currents and mechanical vibrations inherent in the measurement process which must also be accounted for in the synthetic image reconstruction. This perturbation adds another phase term for the input pupil function $P_i(\mathbf{u}_i)=|P_i(\mathbf {u}_i)| \, \exp \left [i \, \phi _i(\mathbf {u}_i) \right ]\exp \left [i \, \phi _d(\mathbf {u}_i) \right ]$, where $\phi _d(\mathbf {u}_i)$ is the experimental phase drift, with total phase $\phi _t(\mathbf {u}_i) = \phi _i(\mathbf {u}_i) + \phi _d(\mathbf {u}_i)$.

As shown in Appendix A, for a thin specimen, the illumination and SHG fields propagate through free space, so that input and output Green functions read $G(\mathbf{r},\mathbf{u}_i)=P_i(\mathbf{u}_i) \, e^{-i \, 2 \pi \, \mathbf{u}_i \cdot \, \mathbf{r}}$ and $H(\mathbf{u}_o,\mathbf{r})=P_o(\mathbf{u}_o) \, e^{- i \, 2 \pi \, \mathbf{u}_o \cdot \, \mathbf{r}}$, respectively. Under the conditions outlined here, the SHG field for a given input frequency $\mathbf {u}_i$, measured in the output pupil plane is given by

A reflection or transmission matrix, for backscattered or transmitted SHG fields, respectively, is defined by sampling the continuous scattering operator in Eq. (2) over the discrete input and output spatial frequency coordinates. The reflection matrix defined for the epi imaging condition can be written as the product of three matrices, $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i} = \mathbf {P}_o \, \hat {\mathbf {\chi }}^{(2)}_{{{\mathbf{q}}}} \, \mathbf {P}_i$. The input and output pupil matrices are defined by the discrete form of the pupil functions, $\mathbf {P}_i = \mathrm {diag}\{ P_i(\mathbf{u}_i) \}$ and $\mathbf {P}_o = \mathrm {diag}\{ P_o(\mathbf{u}_o) \}$, respectively. The object susceptibility spectrum is a Toeplitz structure that is given by $\hat {\mathbf {\chi }}^{(2)}_{{{\mathbf{q}}}} = \hat {\chi }^{(2)}({{\mathbf{q}}})$ taking the form $\hat {\mathbf {\chi }}^{(2)}_{{{\mathbf{q}}}} = \hat {\chi }^{(2)}(\mathbf {u}_o+\mathbf {u}_i)$. This matrix can alternatively be constructed with $\hat {\mathbf {\chi }}^{(2)}_{{{\mathbf{q}}}} = \mathbf {F} \, \mathrm {diag}\{ \chi ^{(2)}(\mathbf{r}) \} \, \mathbf {F}^{-1}$, where the susceptibility matrix has been flattened into a one-dimensional vector before being placed on the matrix diagonal. Here, $\mathbf {F}$ and $\mathbf {F}^{-1}$ are the discrete Fourier and inverse Fourier transforms operators, respectively.

These reflection and transmission matrices map the input spatial frequency coordinate, $\mathbf {u}_i$, to the output spatial frequency coordinate, $\mathbf {u}_o$. Scattering from the object probes the object spatial frequency so that in the output pupil plane the scattered field is proportional to the complex spatial frequency distribution of the second order susceptibility of the sample, but is shifted according to the tilt of the input plane wave. Once the transmission or reflection matrix is constructed, we can obtain the synthetic SHG image field from a shifted form of the matrix.

The synthetic SHG image field can be constructed by shifting the columns of the reflection matrix to line up the scattered fields, $E_{\rm SHG} (\mathbf{u}_o,\mathbf{u}_i)$, with respect to $\hat {\chi }^{(2)}({{\mathbf{q}}})$. This shifted operator reads

Similarly, as illustrated in Fig. 2, a reversal synthetic aperture object spectrum can be formed in the input pupil coordinates by first taking the transpose of the reflection matrix $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}$ (i.e., by swapping the input and output spaces) and then shifting the columns again in the same manner as discussed above. This operation will will produce the operator

Optical aberrations appearing in the form of phase aberrations in the input pupil, $\phi _i(\mathbf{u}_i)$, and the output pupil, $\phi _o(\mathbf{u}_o)$, lead to distortions in the synthesized image. These phase distortions can be estimated and corrected using redundancy in the reflection and transmission matrices. Previous work in linear scattering has demonstrated that correlations of the output spatial frequency spectrum between closely spaced input spatial frequency measurements provides a good estimate of input pupil phase difference at the mean of the two input spatial frequency points [39–41].

Here, we present a straightforward and effective algorithm that estimates and corrects aberrations in the synthetic aperture holographic images by determining the input and output pupil phase. The estimation of pupil phase involves utilizing the singular value decomposition (SVD) of the matrices $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$ and $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$. The formation of these matrices introduces strong correlations among the columns over a wide range. Consequently, the SVD is well-suited for this scenario as it identifies the eigenvectors of the correlation matrices [41]. The SVD is given by $\mathbf {D} = \mathbf {U} \, \Sigma \, \mathbf {V}^\dagger = \sum _j \sigma _{j} \, \mathbf {u}_{j} \, \mathbf {v}_{j}^\dagger$. The left singular vectors, $\mathbf {v}_{j}$, are columns in $\mathbf {V}$, and are eigenvectors of the correlation matrix $\mathbf {D}^{\dagger }\mathbf {D}$. Similarly, $\mathbf {u}_{j}$, are columns in $\mathbf {U}$, and are the right eigenvectors of the other correlation matrix $\mathbf {D}\mathbf {D}^{\dagger }$. These singular vectors are paired with the singular values, $\sigma _{j}$, which are listed in decreasing order along the diagonal of $\Sigma = \mathrm {diag}\{ \sigma _{j} \}$, with eigenvalues given by $(\sigma _{j})^2$.

The matrices $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$ and $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$ are arranged such that the synthetic aperture spectrum is reconstructed (in either the forward or reversed direction) by simply summing the columns. Unfortunately, each field (column) of the raw experimental data is out of phase with one another according to both ${{\mathbf{P}}}_o$ and ${{\mathbf{P}}}_i$, as illustrated in Fig. 1. Choosing two neighboring columns of $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$: $\mathbf {d}_{q,\mathbf {u}_i^{1}}={{\mathbf{P}}}_o(\mathbf {q}+\mathbf {u}_i^1)\hat {\chi }^{(2)}({{\mathbf{q}}}){{\mathbf{P}}}_i({{\mathbf{u}}}_i^1)$ and $\mathbf {d}_{q,\mathbf {u}_i^{2}}={{\mathbf{P}}}_o(\mathbf {q}+\mathbf {u}_i^2)\hat {\chi }^{(2)}({{\mathbf{q}}}){{\mathbf{P}}}_i({{\mathbf{u}}}_i^2)$, if the difference in input angle between the two columns is sufficiently small such that $P_o(\mathbf {q}+\mathbf {u}_i^1) \approx P_o(\mathbf {q}+\mathbf {u}_i^2)$, then the phase difference between the two columns is approximately just a piston phase shift set by the phase difference between $P_i({{\mathbf{u}}}_i^1)$ and $P_i({{\mathbf{u}}}_i^2)$. This then nearly isolates the input and output pupils and allows for the problem to be written as a simple matrix operation: $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}\mathbf {a}=E_{\rm SHG}^{\rm s} ({{\mathbf{q}}})$, which gives an explicit expression of the discrete summation form of synthetic SHG field spectrum with phase correction imparted by $\mathbf {a}$ so that the phase of each column is shifted to eliminate aberrations: $\mathbf {a}=e^{i \, \mathbf{\phi _c}(\mathbf {u}_i)}$, with $\phi _c$ being the phase correction. We would like to find $\mathbf {a}$ such that it maximizes the total intensity of $E_{\rm SHG}^{\rm s} ({{\mathbf{q}}})$. When the total intensity is maximum, all the columns (fields) are in phase. This occurs when $\mathbf {a} = P_i(\mathbf{u}_i)^*$, implying that $\phi _c= - \phi _t(\mathbf {u}_i)$, thereby correcting the aberrations imparted by the input pupil. A comparison of the performance of the cross correlation approach and the SVD phase estimate is provided in Appendix D.

 

Fig. 1. Conceptual diagram of SHG synthetic aperture holography represented in a transmission geometry, but which equally applies to a reflection geometry. a) The input fundamental field is focused to a point in the input pupil plane at the input spatial frequency coordinate $\mathbf {u}_i$. When this illumination is at the origin of the input pupil plane coordinates, the fundamental illumination beam is a normally incident plane wave. The scattered field, analogous to a linear transillumination field, is collected by the output pupil and the complex signal field is recorded. b) A second input field example shows an SHG darkfield configuration in which the fundamental beam is incident on the sample at an angle determined by $\mathbf {u}_i$. c) The input illumination angle is scanned across the input pupil to collect SHG scattered fields from a range of object spatial frequency distribution with each scattered spectra aligned to $\hat {\chi }^{(2)}({{\mathbf{q}}})$ showing an enhanced frequency support. d) The full spatial frequency spectrum of the object, $\hat {\chi }^{(2)}({{\mathbf{q}}})$, is estimated from the coherent sum of the recorded spectral field. Aberrations from the input, $P_i(\mathbf{u}_i)$, and output, $P_o(\mathbf{u}_o)$, pupils distort the estimated object spectrum and must be corrected to produce aberration-free images.

Download Full Size | PDF

To motivate this algorithm, we consider an infinitesimal scattering point on axis. Such a scatterer produces a uniform spatial frequency distribution $\hat {\chi }^{(2)}({{\mathbf{q}}}) = 1$. Consequentially, the reflection matrix is rank one and formed by the outer product $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i} = \mathbf {\mathbf {P}_o} \, \mathbf {\mathbf {P}_i}^T$, where the pupils are represented as vectors after flattening with suitable lexicographic ordering. For this simple case, the right singular vector is associated with the input pupil function and the left singular vector is associated with the conjugate of the output pupil function. It is evident in this case that the input and output pupils can be obtained directly from the SVD. A more complex scattering spectrum prevents such clean separation of the input and output pupils from the reflection matrix so that we have to use the shifted reflection matrices to highlight the correlations.

Using the method of Lagrange multipliers it can be shown that the vector $\mathbf {a}$, that maximizes the total intensity of $E_{\rm SHG}^{\rm s} ({{\mathbf{q}}})$ (subject to the constraint that $\mathbf {a}$ is a unit vector), is the left singular vector of $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$ corresponding to the largest singular value. As shown in Appendix C, if ${{\mathbf{v}}}_1$ is the dominant left singular vector of $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$ then optimal phase conjugate occurs for $\mathbf {a}={{\mathbf{v}}}_1$. This proof clarifies the role of the SVD algorithm in providing an excellent phase correction for the input pupil. We find that this algorithm performs extremely well, even under very low SNR conditions, which are difficult to avoid when measuring backward generated SHG; see Appendix D. To find the phase correction for the output pupil, the same process is carried out, except we transform $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$ to $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$ and then take the dominant left singular vector of $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$ as the estimate of the output pupil correction.

Because the input and output pupils are only approximately separable using the shifted representations of the reflection matrix, the algorithm proceeds iteratively, with iteration index denoted by $k$. At each iteration, the reflection matrix, $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}^{(k)}$, is corrected with the input and output pupil phase estimated obtained in this iterative step. For the first iteration, the reflection matrix is initialized with the reflection matrix obtained from the data, $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}^{(0)} = \mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}$. Then, the input pupil phase is estimated by taking the phase argument of the dominant left singular vector given by the SVD of the shifted reflection matrix $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$. The estimated input pupil phase is taken from the phase of the dominant left singular vector of $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$, $\delta \tilde {\phi }^{(k)}_{c,i} = \angle \mathbf {v}_{1}$. The estimated phase correction is then applied and then the matrix is transformed from $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$ to $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$. The output pupil phase is then estimated similarly by the dominant left singular vector of $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$, $\delta \tilde {\phi }^{(k)}_{c,o} = -\angle \mathbf {v}_{1}$. Transforming $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$ back to $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}$ after applying the output pupil phase correction provides the corrected reflection matrix $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}^{(k+1)}$. This pair of steps is counted as one iteration. These operations are shown graphically in Fig. 3. The total phase correction is estimated from $\tilde {\phi }_{c,i} = \sum _k \delta \tilde {\phi }^{(k)}_{c,i}$ and similarly, $\tilde {\phi }_{c,o} = \sum _k \delta \tilde {\phi }^{(k)}_{c,o}$ for the input and output pupils, respectively. A stopping criterion for the number of iterations is established once a maximum total intensity for the reconstructed image is reached. Other metrics have been investigated such as sharpness, but it was found this these reach a peak at the same point.

3. Methods

The strategy outlined above was implemented experimentally after validation and testing of the algorithm through simulations. The experimental system allows for epi and transmission synthetic SHG holograms to be recorded. The SHG field is extracted from the set of holographic intensity patterns and used to build the reflection (and transmission) matrices. These data are then processed according to the algorithm discussed in the previous section to synthesize an enhanced SHG spectrum that is free from optical aberrations, producing aberration-free SHG field images in forward and backscattered configurations with resolution higher than the diffraction limit.

3.1 Experimental setup

The experimental system, as shown in Fig. 4, is driven by a home built Yb:fiber-amplifier system that produces ultrafast laser pulses at a 62 MHz repetition rate that are centered at a $1050$-nm wavelength, with a bandwidth that supports $< 35$-fs transform-limited pulses and an average power of $3.5 - 4$ W ($56-64$ nJ pulse energy). Power in the beam is split into signal and reference arms with a combination of a half waveplate and a polarizing beam splitter. The signal arm is sent through two galvanometric mirrors (Nutfield QS-7) that are relay imaged to one another with a 4-f telescope ($1.33 \times$ magnification). Finally the beam is brought to a focus in the back focal plane of a condenser lens with another lens ($f=200$ mm) to generate a plane wave incident on the sample. To avoid damage in the back focal plane of the lens, as we observed with a microsope objective, an aspheric lens (New Focus 5726, $NA=.16$) serves as the condenser. In the epi direction, the same lens is used to collect the backscattered SHG light. An identical aspheric lens is used in transmission. In both the epi and transmission arms, the SHG signal is isolated with a dichroic optical filter (Semrock FF875-Di01) – rejecting the pump pulse. Typically $3 - 3.4$ W of average power is used for illumination in the fundamental arm. The beam is focused to a collimated beam, i.e., approximately a plane wave, with a transverse width of $500 \, \mu$m. The polarizer in front of the camera is a Meadowlark GPM-150-UNC broadband polarizer. Specimens used for the experiments include a sparse field of $\sim 80-100$-nm diameter Bismuth ferrite, or BiFeO$_3$, (BFO) nanoparticles and a 10$\mu m$ thick section of sheep tendon. The BFO nanoparticles have an average particle size of 80-100 nm and were purchased from Nanoshel and used without any further purification. BFO nanoparticles are suspended in DI water to make 0.5 mg/mL colloidal solutions and sonicated for 15 minutes. A small portion ($<10 \mu L$) of the colloidal solution is then drop cast onto a glass slide and left to dry in ambient conditions before imaging. The result is that the BFO particles are well dispersed on the slide, but still in a high enough density to see constellations of BFO nanoparticles in the images, see Fig. 5. To eliminate the possibility of reflections coming back into the camera on the epi side, the samples are mounted on a glass slide without a coverslip and oriented so that the sample lies on the distal side of the glass relative to the condenser lens.

 

Fig. 2. Conceptual diagram of SHG synthetic aperture holography. a) The input fundamental field is focused to a single point in the input pupil plane. This field produces plane wave illumination of the sample. The optical imaging system filters the SHG field spatial frequency support by the output pupil, $P_o(\mathbf{u}_o)$, which is applied to a portion of the SHG object spectrum centered on $\mathbf {u}_i$. b) Each measured SHG spectrum is flattened into a vector and stacked into a matrix $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}$. c) The conjugate transpose of $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}$ behaves conceptually as a time-reversal experiment. This time-reversal matrix describes a scenario interpreted as an SHG field from the output pupil that backpropagates through the system to the input pupil plane. d) The time-reversal of the data can be realized by taking the conjugate transpose of $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}$.

Download Full Size | PDF

 

Fig. 3. Conceptual diagram of operations needed to find pupil phase corrections. a) Taking the constructed reflection matrix and aligning the output spectra (aligning the columns of $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}$) constructs the matrix $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$. The aberrated input pupil phase is estimated from the phase of the dominant left singular vector of SVD of $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$, giving $\phi _i(\mathbf{u}_i) \approx \angle \mathbf {v}_{1}$. b) Similarly for the output pupil phase correction first the conjugate transpose of $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}$ is taken, then the spectra are again aligned (aligning the columns of $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}^\dagger$) to form $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$. The output pupil phase is estimated from the phase of the dominant left singular vector of SVD of $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$, which reads $\phi _o(\mathbf{u}_o) \approx - \angle \mathbf {v}_{1}$.

Download Full Size | PDF

 

Fig. 4. Diagram of the experimental layout for measuring the SHG reflection and transmission matrix for a set of plane wave illuminations at different angles set by the galvanometric scan mirrors.

Download Full Size | PDF

 

Fig. 5. Experimental results for SHG synthetic aperture with aberration correction in the epi and transmission configurations for a field of BFO nanoparticles. a) uncorrected synthetic aperture SHG image intensity in transmission. b) corrected synthetic aperture SHG image intensity in transmission. c) estimated input and output pupil phase. Because each input angle is a separate measurement, there is a uniformly distributed random phase on top of the optical pupil phase of the illumination condenser optic. d) uncorrected synthetic aperture SHG image intensity in reflection. e) corrected synthetic aperture SHG image intensity in reflection. f) estimated output and input pupil phase in reflection. The estimated SNR from uncorrected to corrected for transmission is 5dB to 27dB respectively. Similarly for reflection the estimated SNR values go from 4dB to 37dB. Scale bar is 50$\mu m$

Download Full Size | PDF

Meanwhile, the reference beam is generated by frequency doubling approximately $35-70$ mW of the laser power by placing the BBO crystal in the center of a telescope with magnification equal to one. The power of the reference beam can be further attenuated with a half-wave plate before being combined with the signal field at the beam splitter. The telescope also collimates the SHG reference beam and a dichroic filter beamsplitter as well as several filters (Semrock FF875-Di01 Dichroic beam splitter, Semrock 720 short pass filter and a BG39 filter) that are used to isolate the reference beam from the residual un-converted fundamental beam. The same models of dichroic beam splitter and filters are used to isolate the fundamental and SHG field generated at the sample plane before imaging. The reference beam is sent through a mechanical delay line that allows for adjustment of the arrive time of the reference SHG pulse to gate the holographic measurement to a particular depth. The signal SHG field is combined with the reference beam with a non-polarizing beam splitter. An image of the SHG field is formed with a tube lens ($f=200$ mm) in the epi arm and ($f=250$ mm) in the transmission arm. Off-axis holographic images are captured with a Hammamatsu ORCA Quest C15550 in the epi arm and a Teledyne prime 95B in the transmission arm.

3.2 Constructing the reflection matrix

For synthetic aperture SHG holography, we record a sequence of $M$ SHG holograms, from which we extract the SHG field in the output spatial image plane for a sequence of input illumination angles, denoted by $\mathbf {u}_i$. Each of the $M$ holograms are captured for a distinct point in the input pupil plane, $\mathbf {u}_i$, which corresponds to a particular incident angle on the specimen. The pair of galvanometer scan mirrors are used to control the incident angle of the fundamental beam by relay imaging the surface of the second galvanometer mirror to the sample plane. The incident angle is controlled by setting the voltage on each galvanometer. A calibration of voltage to resulting plane wave tilt in spatial frequency units is implemented by finding the voltage applied to each galvanometric-mirror required to reach each edge of the pupil. This voltage then corresponds to a spatial frequency of $\mathbf {u}_i = \mathrm {NA}/\lambda _1$, that also outlines the condenser lens pupil boundary. The offset voltage needed to center the illumination path on the objective pupil plane is determined. The pupil boundary voltages are then used to map the control voltages to the input spatial frequency values for each measured hologram that is acquired for a specific set of galvanometer control voltages. The knowledge of the input spatial frequency decoded from the control voltage applied to the galvanometer is used to compute the required shift to align each of the output SHG spectra.

Holograms of the scattered SHG field are recorded for each incident fundamental illumination angle. The holograms are captured with a planar off-axis reference beam, so that the SHG field is extracted with a conventional holographic processing algorithm [25]. The extracted field is spatially cropped to limit the image field of view to a total number of pixels $N$. Each cropped output SHG field is transformed to the output pupil plane, with coordinates $\mathbf {u}_o$. These measured fields are flattened according to lexicographic order into a linear vector truncated to length $N$ in order to manage computer memory usage. Each SHG field (now represented as a vector) is stored in the columns of a transmission or reflection matrix for the transmission and epi SHG fields, respectively. The columns of the matrix are filled with the input spatial frequencies ordered in the same lexicographic order as the flattened output spatial frequency vectors. In this way, a column of the reflection matrix corresponds to a measured output field due to a certain input spatial frequency. A row corresponds to the detected scattered complex-valued SHG field of a certain output spatial frequency (pixel) according to the input spatial frequency. Remapping a row of ${{\mathbf{R}}}$ to a 2D array in the correct ordering yields the spatial frequency spectrum corresponding to the reversal of a plane wave through the system. In other words, mimicking the process of sending a second harmonic plane wave from the output pupil plane to the input pupil plane. The resulting matrix is of size $N \times M$, with columns mapping the output spatial frequency coordinates and the rows indicating the input spatial frequency for illumination.

4. Results

Both samples were imaged in epi and transmission and configurations. In all cases, the SHG field reflection or transmission matrix is recorded for a set of input spatial frequencies, corresponding to a sweep of incident angle of the fundamental beam. These data are then processed to estimate and correct for the input and output pupil phases after which aberration-free synthetic SHG spatial frequency spectrum and images are obtained.

The scattered SHG field from sub-wavelength particles is of similar power in the forward and backward directions. In contrast, for sheep tendon the scattered SHG power is reduced by at least an order of magnitude in the backward direction as compared to forward scattered SHG. As SHG scattering is already a weak process, the scattered SHG fields are relatively weak, and particularly weak in the backscattered direction from the sheep tendon samples. Thus, measuring the backward scattered SHG field is a challenge. On average, the signal in the epi direction from sheep tendon is about 2000 to 4000 photons per second. Whereas in the transmission direction is is around 30000 to 60000 photons per second. These numbers where arrived by taking widefield intensity images and converting the gray scale value (ADU) to number of photons according to the specs of each camera. Fortunately, SHG holography allows for the measurement of a weak field by leveraging the heterodyne enhancement from a strong reference beam [30]. This enhancement enabled us to record SHG widefield holograms in the epi direction. Coherent addition of the fields also aids in increasing the SHG scattered field strength. The synthetic summation of SHG scattered field spectra (over the illumination angles) enables the SHG signal to grow linearly with the number, $M$, of coherent fields added together.

4.1 Coherent averaging

While these strategies enable measurement of backward emitted SHG from the sample in a widefield imaging configuration, the signal is still quite low. We can take another step to further boost this signal by exploiting phase information in hand. Normally, to increase a signal measured on a camera one could just take the average of many measurements or increase the exposure time on the camera. Unfortunately, in the case of holography, increasing exposure time usually fails because the signal of interest is retrieved by analyzing the fringe pattern produced on the camera due to the interference of the signal and reference beam. This fringe pattern is extraordinarily sensitive to air currents, vibrations, and other perturbations to the accumulated relative phase between the reference and signal fields in the non-common path regions of the reference and signal arms. The fringe visibility, and thus the signal, degrades when averaging several holograms or increasing the exposure time due to the fluctuation in the relative phase over the integration timescale.

We have developed a simple strategy to mitigate these relative phase fluctuations – enabling a significant boost in the signal-to-noise ratio of the SHG holographic field measurement. This strategy again leverages coherent field summation. To implement the coherent sum boost, a sequence of holograms with short camera exposure times is taken for each incident angle, i.e., the same input spatial frequency, for the fundamental illumination beam. The SHG field is extracted using standard off-axis holographic processing [25]. These fields are cropped, flattened, and stacked into the columns of a matrix, ${{\mathbf{A}}}$. Since these measurements are all taken at the same input spatial frequency, their output spatial frequency spectra are already aligned. Thus, each SHG spectrum is nominally identical, except for the changes in overall phase that arises as a result of the shot-to-shot changes in relative signal-reference beam phase. Taking inspiration from the algorithm we developed for aberration correction, we find the phase offset between each measured SHG field them using the SVD. The phase correction vector, ${{\mathbf{c}}} = \mathbf {v}_{1}/|\mathbf {v}_{1}|$, is taken from the dominant singular vector $\mathbf {v}_{1}$, of the SVD of ${{\mathbf{A}}}$. Further improvement in the SNR can be obtained by filtering out noise in ${{\mathbf{A}}}$ with a truncated SVD. The coherent sum of the SHG spectra for each input angle is given simply by ${{\mathbf{A}}} \, {{\mathbf{c}}}^*$, which boosts the SHG signal field by the number of hologram measurements. To obtain an equivalent enhancement in signal by averaging the intensity on the camera would require perfectly stable fringes. With the approach outlined here, we are able to form a widefield image while also correcting aberrations even with an extraordinarily weak signal, highlighting the power of this technique. This coherent averaging process is repeated for each fundamental beam illumination angle, i.e., each input spatial frequency $\mathbf {u}_i$, and the coherently summed SHG field for each angle is stacked into a matrix to build the reflection/transmission matrix $\mathbf {R}_{\mathbf{u}_o,\mathbf{u}_i}$.

4.2 Imaging conditions

For the data collection of SHG synthetic aperture imaging of the sheep tendon the imaging conditions are as follows. For both the epi and transmission data a total of 2601 holograms were recorded which evenly sampled the full support of the illumination aperture. The exposure time on the camera was $.1$ s for both configurations. The signal is much weaker in the epi direction so for each angle 15 holograms were coherently summed. For the transmission data 5 holograms were coherently summed for each angle. Under these conditions the total acquisition time for epi collection was 65 minutes for epi and 22 minutes for transmission. We stress that this system was not optimized for speed. To increase the speed, the exposure time could be reduced if a laser with a very high pulse energy were used as in [32]. The number of acquired holograms can also be reduced using compressed sensing techniques. Similarly, the conditions used for the BFO nanoparticle data were the same in the epi and transmission directions. The exposure time set to $.1 sec$ and 10 holograms were coherently summed for each angle with a total of 2601 holograms evenly sampling the input aperture. This leads to a total acquisition time of 43 minutes for both the epi and transmission data for the BFO particles.

4.3 Reconstructions

Before and after images illustrating the correction of aberrations are shown in Fig. 5 for SHG-active BFO nanoparticles and Fig. 6 for thin sheep tendon slices. The input and output pupil phase estimates are also shown. Due to experimental phase drifts in the system (air currents, mechanical vibrations), each measurement is dephased with each other measurement by a random offset phase. With no correction of the relative random phase fluctuation, the resulting synthetic aperture reconstruction has very low SNR (shown in Figs. 5 and 6). This experimental phase drift is corrected along with the input pupil phase correction simultaneously. The resulting input pupil phase map is the superposition of the optical aberrations and the phase drift.

 

Fig. 6. Experimental results for SHG synthetic aperture with aberration correction in the epi and transmission configurations for a 10$\mu m$ thick section of sheep tendon. Images a-d correspond to the epi configuration and e-h in transmission. Images a) and e) show the reconstruction before any corrections are applied. Images b) and f) show the intensity after correction with the input and output pupil phase corrections inset. Panels c) and g) show the phase of the synthetic aperture reconstructions after the corrections are applied. The reconstructed spectrum in images d) and h) have white dashed outlines showing the original spatial frequency support of the system shifted to different positions stitching together an expanded spatial frequency support containing more information. The dashed circles only show a few example positions, in reality there were 2601 measurements taken. For the synthetic aperture reconstruction in the transmission configuration (images e and f) the estimated SNR goes from 12.3dB to 35.8dB after correcting for experimental phase drift and aberrations. Similarly for the epi configuration (a and b) the SNR goes from 13.5dB to 34.8dB. SNR estimates were made by considering the singular values up to an optimal truncation point as signal and the rest noise. The optimal truncation point is determined using the convention established by Gavish and Donahoe [42]. Scale bar is 50$\mu m$

Download Full Size | PDF

BFO nanoparticles were used to validate the performance of the estimation and correction of the input and output pupil phases. In addition, these BFO nanoparticles are well below the imaging resolution, and thus isolated BFO nanoparticle images serve to report the amplitude spread function of the synthetic aperture SHG holographic imaging, highlighting improvements as we correct aberrations. Fig. 5(a) shows the uncorrected synthetic aperture SHG image intensity when the input and output pupil phases are not corrected for the transmission image of BFO nanoparticles and Fig. 5(d) shows corresponding image obtained in the epi direction. The aberration-free images obtained by estimating and correcting for the input and pupil phases are shown in Figs. 5(b) and 5(e) for the transmitted and backscattered SHG field, respectively.

The uncorrected and corrected aberration corrected synthetic aperture reconstructions of thin sheep tendon samples in the epi and transmission directions are shown in Fig. 6. The results show a dramatic improvement in image quality after processing with the SVD-based aberration correction algorithm. While the estimation of SNR from an image is a difficult problem, the distribution of SVD values have been shown to provide a robust strategy for SNR estimation [42]. The SNR values of the images are estimated based on this strategy, with results shown in the caption of Fig. 6.

The results of complex aberration-free synthetic aperture SHG fields recovered from SHG synthetic aperture holographic imaging imaging of a thin sheep tendon in transmission and reflection geometries of the same region are shown in Fig. 6. We display the estimated input and output pupil phases inset within the corrected intensity images Figs. 6(b) and 6(f). Note that these measurements were taken in independent data runs, so there is no correlation between the shot-to-shot random phase fluctuations of the input pupil phases. These results show extremely robust performance of synthetic aperture SHG holographic imaging.

4.4 Resolution

We compare the resolution shown in Fig. 7 with the expected resolution based on the properties of the optical imaging system and the wavelengths used in the experiment. In our experiment, we have a numerical aperture of a input condenser lens $\mathrm {NA}_i=.16$ and a fundamental illumination wavelength $\lambda _1=1.05 \, \mu$m and an output objective lens with numerical aperture $\mathrm {NA}_o=.16$ and an imaging second harmonic wavelength $\lambda _2=\lambda _1/2=.525 \, \mu$m. The NA for this system is low partly to have a large working distance to add large samples and for future experiments imaging in deep tissue.

 

Fig. 7. Spatial frequency maps comparing the resolution of a single hologram on axis a), the synthesized aperture without phase corrections b), and the aberration corrected synthetic aperture image c) for the sheep tendon sample in transmission. The dashed lines show the approximate location of the cutoff frequency for each case. This demonstrates that the resolution is increased by a factor of two for the synthesized aperture.

Download Full Size | PDF

For a single hologram image, the spatial frequency support of the image, and thus the spatial resolution, are set by $\Delta u_o=2 \, \mathrm {NA}_o/\lambda _2$ and $\delta x_o=1/\Delta u_o=\lambda _2/(2 \, \mathrm {NA}_o)$, respectively. The radial spatial frequency support for a single hologram image is given by $\Delta u_o/2=\mathrm {NA}_o/\lambda _2$. In our experimental setup, these number evaluate to $\Delta u_o=.61 \, \mu$m$^{-1}$, with radial spatial frequency support of $.305 \, \mu$m$^{-1}$ (this is the radius of the image spectral amplitude) and a spatial resolution of $\delta x_o=1.6 \, \mu$m. As can be seen in the data shown in Fig. 7, we obtain the expected spatial frequency support and resolution. A calibration was done for the magnification of the system using a Ronchi ruling of a known line spacing in order to have a calibrated pixel size.

By scanning the input angle, we can collect portions of the object spectrum which are higher than this cutoff to record synthetic aperture holography data. In the case of a synthetic aperture experiment, the input spatial frequency support of the input fundamental beam can be scanned over a range of $\pm \mathrm {NA}_i/\lambda _1$ of the input fundamental wavevector, but this effect is doubled in the SHG scattering process. Thus, the spatial frequency support for the synthetic aperture SHG image spans $\Delta u_o=2 \, \mathrm {NA}_o /\lambda _2 + 4 \mathrm {NA}_2/(4 \lambda _1)$. The experimental spatial frequency support for the synthetic aperture SHG field image evaluates to $\Delta u_o=1.22 \, \mu m^{-1}$, with radial spatial frequency support of $.61 \, \mu m^{-1}$, and a spatial resolution of $\delta x_o=.8 \, \mu$m. We see that the aberration-corrected image in Fig. 7(c) matches the expected resolution criterion, whereas the uncorrected, and thus aberrated, image does not provide full spatial frequency support or achieve the minimum spatial resolution. This is a result of the overlapping portions of the spectrum for each hologram being added are out of phase.

5. Discussion

We have adapted methods that were first developed to improve imaging distorted by optical aberrations and linear coherent scattering [37,39–41] for application to nonlinear holographic imaging. [25,27,29–31] Nonlinear scattering is dominated by forward scattering due to phase-matching considerations [31]. As a result, the backward scattered field is quite weak – presenting a significant challenge for detection of the SHG field imaged in an epi configuration. Our new imaging strategy makes use of two coherent summations of the complex-valued SHG field that is recovered from off-axis SHG holography in the forward-scattered (transmission) and backward-scattered (epi) geometries. These coherent combining methods rely critically on the ability to estimate phase differences between a set of measurements. The coherent summation can overcome the inherently weak SHG field strength to produce aberration-free complex-valued SHG field images with increased spatial frequency support, and thus improved spatial resolution.

We have successfully been able to form an aberration corrected synthetic aperture image for the back-scattered as well as forward-scattered SHG from BFO nanoparticles and sheep tendon. The detected signal is boosted using coherent amplification of the field that occurs from heterodyne mixing between the signal field and a reference field in a holographic measurement [30]. The epi-SHG holographic images shown here provide the first complex-valued nonlinear bacskcattered optical field measurements. To ensure that there was no contamination from forward-scattered SHG radiation that is directed in a backward direction, we eliminated all material from the distal end of the sample to prevent Fresnel reflection of forward-scattered SHG fields into the epi SHG holographic imaging system. In addition, the SHG signal and reference fields are broad-bandwidth, with a cross-coherence length of $\sim 8.84\mu m$. This epi SHG image carries the advantage of gating out any stray reflections and will ultimately enable three-dimensional imaging because the low-coherence interferometry provides optical sectioning of the backscattered SHG field.

Even with the coherent heterodyne amplification provided by holographic imaging with a strong reference field, the SNR of the field extracted from epi-SHG holography was relatively low. One strategy to boost the SNR is to simply increase the integration time of the camera, if there is still some dynamic range left of the camera sensor. Unfortunately, this strategy is infeasible in our configuration because the signal and reference beams are not common path. The lack of common path propagation leads to random relative phase fluctuations over the camera integration time. These random fluctuations degrade the fringe visibility and thus the SNR of the extracted SHG field. To combat this SNR degradation, we implemented a new coherent summation strategy for the SHG field for a set of nominally identical SHG holograms. This algorithm is based on estimating the random relative phase variations of the set of the SHG field extracted from multiple holographic measurements. Implementation of this protocol provides a significant boost in the SHG field SNR. We note that a similar strategy has been adopted for linear holographic imaging through turbulent media [43].

Despite this boost in signal SNR, the SHG images are still quite degraded due to a combination of aberration phases introduced by the input and output optics, as well as due to specimen-induced aberrations from propagation of the fundamental and second harmonic fields in the sample. The imaging distortion from aberrations is exacerbated by the fact that we use an aspheric lens for imaging. While such aspheres are generally avoided due to the presence of strong optical aberrations (off-axis), we used these optics since there are no optical components in the pupil plane of the lens. As this plane is inside of conventional optical objectives, these objectives are damaged when employing such a plane-wave illumination. Our synthetic aperture strategy, in which we accurately extract the input and output pupil phase, allows for the use of low quality, and less expensive, optics for imaging.

To accurately estimate the input and output pupil phases from a set of data, we need to introduce some redundancy in the measurements. Such redundancy is attained by capturing data over a range of input fundamental incident angles so that the set of SHG field measurements have partially overlapping spatial frequency support. This set of data contains sufficient redundancy (i.e., spatial frequency redundancy) to estimate the input and output pupil phases. These phases can be extracted with a cross-correlation algorithm that has been used for phase correction in linear scattering based microscopy [39], however, we found that this cross-correlation algorithm performed poorly in the limit of low SNR. A simulation of the robustness of cross-correlation phase estimation as a function of data SNR is discussed in Appendix D. This same analysis shows that our new algorithm for input and output pupil phase and correction based on the SVD of suitably shifted reflection (or transmission) matrices performs extremely well in the presence of low SNR data. Furthermore, we show in Appendix D. that the SVD-based algorithm discovers the optimal phase correction to produce aberration-free SHG field imaging even in the presence of high noise levels. In addition, the SVD algorithm is more computationally efficient than computing the cross correlations for phase estimation. On average it requires about a factor of two fewer iterations and in some cases finds the best phase corrections in a single iteration according to SNR and sharpness metrics.

The application of our new aberration-free synthetic aperture imaging strategy to experimental measurements shows excellent performance. The combination of the coherent signal enhancement and large spatial frequency support allows for the un-distorted estimation of the spatial frequency spectrum of the second-order nonlinear optical susceptibility that gives rise to the coherent nonlinear SHG scattering. The corrections shown for sub-resolution nanoparticles exhibits diffraction-limited image points. High quality amplitude and phase images of thin sheep tendon slices illustrate the power of this new imaging modality. In contrast, for the conventional SHG microscopy case where a focused beam is scanned through the system and the SHG power is measured and assigned to that focal spot. The image formation in the conventional case suffers from aberrations on the illumination because the point spread function is distorted – thus distorting the image. There is no refuge in such a conventional measurement to correct this image from information in the data as we demonstrate here. Furthermore, in our new synthetic aperture imaging, we directly access the nonlinear SHG susceptibility. In contrast, conventional point scanning SHG imaging probes information related to the intensity of scattering generated by the focal volume of the fundamental focused field. The total power measured at each point is related in a complicated way to the local spatial morphology of the nonlinear susceptibility and the phase matching conditions.

Image formation in widefield coherent nonlinear imaging deviates from the conventional strategy for nonlinear optical microscopy. While a full discussion of the novel facets of widefield coherent nonlinear scattering are beyond the scope of this article, we will make a few remarks to highlight the differences in the two approaches. Coherent nonlinear imaging is fundamentally an inverse scattering problem that can be described with a model that is homologous to the first-Born approximation model for linear scattering. For this homology to be valid, the optical interaction must both satisfy the undepleted pump approximation, which requires a weak scattered field, and linear scattering for the fundamental and harmonic fields must be negligible. Under these conditions, the scattered field is driven by a source term that is proportional to the $n^{\rm th}$ power of the fundamental field, where $n$ indicates the order of the optical nonlinearity [31]. For the case of an incident fundamental plane wave, such as we present in this article, the information probed by the object, i.e., the spatial variation of the nonlinear optical susceptibility $\chi ^{(n)}({{\mathbf{r}}})$, in the holographic measurement is determined by ${{\mathbf{Q}}} = {{\mathbf{k}}}_2 - 2 \, {{\mathbf{k}}}_1$. Here ${{\mathbf{k}}}_1 = 2 \, \pi \, \mathbf {u}_i$ is the incident fundamental plane wave direction, with wavenumber $\lVert {{{\mathbf{k}}}_1}\rVert = 2 \, \pi \, n_1/\lambda _1 \equiv \beta _1$, where $n_1$ is the refractive index (RI) of the sample at the fundamental wavelength $\lambda _1$. There are a set of wavevectors of the SHG field, ${{\mathbf{k}}}_2 = 2 \, \pi \, \mathbf {u}_o$, that are determined by the observation direction of the set of output spatial frequencies $\mathbf {u}_o$. The unit vector describing the direction of observation is $\hat {{{\mathbf{s}}}}_2 = {{\mathbf{k}}}_2 / \beta _2$. The wavenumber of the SHG field is $\lVert {{{\mathbf{k}}}_2}\rVert = 2 \, \pi \, n_2/\lambda _2 \equiv \beta _2$, where $\lambda _2 = \lambda _1/2$ is the wavelength of the SHG field. The unit vector describing the direction of illumination is $\hat {{{\mathbf{s}}}}_1 = {{\mathbf{k}}}_1 / \beta _1$. In a linear scattering problem described in the first-Born approximation, measurement of the image field samples of the object spatial frequency distribution as described by the Fourier Diffraction Theorem (FDT), which shows the measured information lies on a portion of the Ewald sphere with the apex located at the origin where ${{\mathbf{Q}}} = 0$.

In our experiment, we study the SHG scattered field. In this case, which we call the nonlinear FDT, the information sampled for an input spatial frequency ${{\mathbf{k}}}_1$ and an observation spatial frequency ${{\mathbf{k}}}_2$ sits at the transverse object spatial frequency value of ${{\mathbf{Q}}}_{\perp } = {{\mathbf{k}}}_{2\perp } -2 \, {{\mathbf{k}}}_{1\perp }$. This is the transverse component of the spatial frequency vector that satisfies ${{\mathbf{Q}}}_{\perp } \cdot \hat {z} = 0$, and where $\hat {z}$ is a unit vector along $z$. The fundamental and second harmonic transverse wavevectors, ${{\mathbf{k}}}_{2\perp }$ and ${{\mathbf{k}}}_{1\perp }$, are similarly defined. The axial spatial frequency information collected for the forward scattered is defined as ${{\mathbf{Q}}}_z = {{\mathbf{Q}}} - {{\mathbf{Q}}}_{\perp } = Q_z \, \hat {z}$, where $Q_z = \sqrt {\beta _2^2 - \lVert {{{\mathbf{k}}}_2}\rVert^2} - 2 \sqrt {\beta _1^2 - \lVert {{{\mathbf{k}}}\rVert_1}^2}$. The spatial frequency information of the object sampled by a holographic measurement is given by the nonlinear FDT, $\hat {\chi }^{(2)}({{\mathbf{Q}}}_{\perp }, Q_z)$, Here, the spatial frequency spectrum of the nonlinear susceptibility is defined by the Fourier transform, $\hat {\chi }^{(2)}({{\mathbf{Q}}}) = \mathcal {F} \{ \chi ^{(2)}({{\mathbf{r}}}) \}$. Similar arguments can be made for the backward scattered fields collected in the epi direction. A full discussion of this analysis will be featured in a future article.

With this introduction in hand, we are in a position to discuss the impact of phase mismatch on the SHG holographic imaging process. Nonlinear optics is notoriously afflicted by the problem of phase mismatch that is defined by $\Delta \beta = \beta _2 - 2 \, \beta _1$. Often the phase mismatch is expressed in terms of a coherence length, $\ell _c = \pi /\Delta \beta$, which is the length scale at which new nonlinear scattering interferes destructively with the nonlinear scattered light generated earlier in the interaction. In classical nonlinear frequency conversion, a finite coherence length is a factor the limits conversion efficiency of power from the fundamental to the nonlinear field. Normally in nonlinear imaging, phase mismatch drastically impacts the spatial features that may be imaged. In conventional point scatting nonlinear SHG microscopy, the power scattered from a small volume from a small point focus is recorded for each focal point. The recorded power depends sensitively on the phase matching conditions of the illumination volume in the sample [34]. Remarkably, widefield SHG holography behaves differently. Because each observation point on the camera involves a particular scattered field direction, only the complex values amplitude of the nonlinear susceptibility at ${{\mathbf{Q}}}_{\perp }$ and ${{\mathbf{Q}}}_z$ are probed, but these values are faithfully measured in the experiment. The impact of the phase mismatch can be seen as a high-pass filter of spatial structure in the sample. To understand this filtering, consider the case of forward-scattered SHG holography where the direction of observation, $\hat {{{\mathbf{s}}}}_2$ is the same as the direction of illumination, $\hat {{{\mathbf{s}}}}_1$. The magnitude of the observed scattering vector of the sample is $\lVert {{{\mathbf{Q}}}}\rVert = \lVert {\beta _2 \, \hat {{{\mathbf{s}}}}_2 - 2 \, \beta _1 \, \hat {{{\mathbf{s}}}}_1}\rVert = \Delta \beta$. This condition is fulfilled regardless of the direction of $\hat {{{\mathbf{s}}}}_j$. As the point of the object spatial frequency that is probed when $\hat {{{\mathbf{s}}}}_2 = \hat {{{\mathbf{s}}}}_1$ is the apex of a portion of the Ewald sphere, no points are sampled inside of a circle with a radius of $\Delta \beta$ in the spatial frequency distribution of the sample spectrum. Consequentially, we conclude that the phase mismatch serves as a high-pass spatial frequency filter in the coherent nonlinear holographic imaging process. Moreover, above the frequency $\lVert {{{\mathbf{Q}}}}\rVert = \Delta \beta$, any spatial frequency structure can be sampled up to a roughly $2 \, \pi /\lambda _2$ in the forward direction.

In the epi (backscattered direction), the situation is slightly different. Here, we probe the transverse spatial frequency components ${{\mathbf{Q}}}{\perp } = - {{\mathbf{k}}}_{2\perp } - 2 {{\mathbf{k}}}_{1\perp }$ and axial spatial frequency components $Q_z =- \sqrt {\beta _2^2 - \lVert {{{\mathbf{k}}}_2}\rVert^2} - 2 \sqrt {\beta _1^2 - \lVert {{{\mathbf{k}}}_1}\rVert^2}$. In the backscattering direction, we can consider the scenario where $\hat {{{\mathbf{s}}}}_2 = -\hat {{{\mathbf{s}}}}_1$ to study the impact of phase mismatch. In this scenario, the magnitude of the observed scattering vector of the sample is $\lVert {{{\mathbf{Q}}}}\rVert = \beta _2 + 2 \, \beta _1 = 2 \, \beta _2 - \Delta \beta$. As generally $\Delta \beta \ll \beta _2$, we observe that the impact of the phase mistmatch on the backscattered SHG field collected in the epi direction is negligible. The consequence of the scattering picture is that signal is only collected for spatial structure that is on the order of the wavelength and smaller. This is the same as the case for linear scattering (e.g., as is observed in linear holography and OCT).

Widefield imaging holds an advantage over point scanning microscopy in that since data is collected on a camera, the data collection occurs over the entire time, $T$, over which the image is acquired. Thus, each of the $N$ pixels in the widefield camera image is integrated for the entire image capture time. In contrast, with point scanning imaging, as is used in conventional nonlinear microscopy, the data for each pixel is captured sequentially, and thus the dwell time for each measured point in the image is reduced by the total number of image pixels to $T/N$. Assuming the local illumination intensity is the same for each case, then the signal photon flux, $\rho _s$, is identical for each scenario. It follows that in a shot-noise limited measurement, the SNR at each pixel of the widefield image captured on a camera is $\mathrm {SNR} = \sqrt {\rho _s \, T}$, whereas the point scanning image signal-to-noise ratio is reduced to $\mathrm {SNR} = \sqrt {\rho _s \, T/N}$ because of the reduced data integration time at each pixel. The widefield image SNR is higher by a factor of $\sqrt {N}$, demonstrating a clear advantage.

Holographic widefield imaging is even more advantageous from an SNR perspective. In the case of holographic imaging, the signal is formed through interference with a reference field. Considering a rate of photon detection for the signal of $\rho _s$ and the reference of $\rho _r$, then the signal-to-noise ratio reads $\mathrm{SNR} = 2 \, \sqrt{\rho_r \, \rho_s} \, T/\sqrt{ \left( \rho_r + \rho_s + 2 \, \sqrt{\rho_r \, \rho_2} \right) \, T }$. For the case of nonlinear holographic imaging, where the signal is generally much weaker than the reference, i.e., $\rho _s \ll \rho _r$, then we reach the shot-noise limit for signal detection $\mathrm {SNR} \approx 2 \sqrt {\rho _s \, T}$ even for very weak signal fields due to the heterodyne amplification captured by the signal term $\mathrm{SNR} = 2 \, \sqrt{ \rho_s \, T}$. A larger reference field means faster detection is possible, with a limit for this scaling is set by the dynamic range of the camera [30].

6. Conclusions

We have presented the first epi SHG holographic images that were enabled by a combination of heterodyne-enhanced signal amplification that is able to boost the weak backscattered SHG signal field, along with a coherent summation strategy to boost the SNR of individual holographic field measurements. The fundamental illumination beam is configured as a plane wave where the incident angle is scanned. The full set of data from both transmission and epi SHG holograms are collected into a scattering matrix. These data exhibit overlap in the measured SHG spatial frequency distributions. This redundancy enabled the robust estimation and correction of the input and output pupil phase that leads to a distortion of the SHG hologram images. Once the scattering matrix is corrected, an aberration-free SHG image field spectrum is estimated with an expanded, synthetic aperture. Results are shown for synthetic aperture SHG holography in both the epi and transmission configurations. This demonstration of epi-collected widefield SHG holographic imaging opens a new path for minimally invasive imaging in scattering media with aberration-corrected SHG holography.

This new form of holography opens possibilities for imaging a wide range of samples and in difficult experimental environments. The ability to estimate and correct for aberrations from the measured field data opens the possibility of using less expensive, or generalized optics in the design of widefield nonlinear imaging systems. Such design flexibility could significantly expand the reach of nonlinear microscopy. In addition, distortion from optical scattering could also be corrected with the approach presented here, opening the possibility of computational correction for scattering in nonlinear optical imaging. Future work will expand aberration-free SHG holography imaging to capture images of objects deep inside of scattering media, such as biological tissues. In addition, this approach will be applied to other nonlinear interactions, such as third harmonic generation (THG) and coherent anti-Stokes Raman scattering (CARS). Moreover, the phase information obtained from coherent nonlinear interactions is relatively unexplored and could prove valuable for a wide range of measurements in biological and materials studies.

A. Green’s function for SHG excitation and detection

Greens’s functions used in the formulae for the SHG field forward and backward scattering configurations defined in Eq. (1) are derived. As is evident in Figs. 1 and 4, the input illumination field Green’s function is obtained from a map from the input spatial frequency plane coordinates, $\mathbf {u}_i$, to the coordinates in the sample plane, $\mathbf {r}$ and the output Green’s function is a map from the sample plane coordinates to the output pupil plane spatial frequency coordinates, $\mathbf {u}_o$. In both cases, the map from input coordinates to the sample plane coordinates and from the sample plane coordinates to the output plane coordinates are accomplished with a 2-f optical system. The relevant Green’s functions are defined below using the notation in the classic optical textbook by Mertz [44]. Following this notation, we will use the wavenumber defined by $\kappa _j = 1/\lambda _j$, for a field at the optical wavelength $\lambda _j$.

A.1. Input Green’s function

The input fundamental field, with wavelength $\lambda _1$, is focused within the input pupil spatial coordinates $\mathbf {x}_i$. Suppose that this field is incident on the front focal plane of the illumination condenser lens with focal length $f_c$ and is denoted by $E_{i}(\mathbf {x}_i)$. The fundamental field in the sample plane that is incident on the sample placed in the back focal plane is given by

(5)$$E_1({{\mathbf{r}}}) ={-} i \, \frac{\kappa_1}{f_c} \, \int \, P_i(\mathbf{x}_i) \, E_{i}(\mathbf{x}_i) \, \exp\left( - i \, 2 \, \pi \, \frac{\kappa_1}{f_c} \, \mathbf{x}_i \cdot {{\mathbf{r}}} \right) \, d^2 \mathbf{x}_i.$$

The fundamental field at the sample plane excites a second-order dipole oscillation that drives scattering at the second harmonic frequency of $\omega _2 = 2 \, \omega _1$, which appears at a wavelength of $\lambda _2 = \lambda _1/2$. Synthetic aperture holographic imaging uses an illumination with a point focus in the input pupil $P_i(\mathbf {x}_i)$ at a spatial coordinate $\mathbf {x}_s$ with a field that is approximated as a 2-D Dirac delta function, $E_{i}(\mathbf {x}_i) = \delta ^{(2)}(\mathbf {x}_i - \mathbf {x}_s)$. With this input field, we have a fundamental plane wave incident on the sample of

(6)$$E_1({{\mathbf{r}}}) = P_i(\mathbf{x}_i) \, \exp\left( - i \, 2 \, \pi \, \frac{\kappa_1}{f_c} \, \mathbf{x}_s \cdot {{\mathbf{r}}} \right).$$

The scattered SHG field is driven by the square of the incident fundamental field, $E_1^2({{\mathbf{r}}})$, from which we define the input Green’s function

(7)$$G(\mathbf{r},\mathbf{u}_i) = P_i(\mathbf{u}_i) \, e^{- i \, 2 \, \pi \, \mathbf{u}_i \cdot {{\mathbf{r}}}},$$

where we have defined the effective input pupil spatial frequency $\mathbf {u}_i = 2 \, \mathbf {x}_s/(\lambda _1 \, f_c) = \mathbf {x}_s/(\lambda _2 \, f_c)$ and we have assumed that the amplitude support of the input pupil is binary. Here we have suppressed scaling factors in favor of compact notation.

A.2. Output Green’s function

The output SHG field, with wavelength $\lambda _2$, is mapped from the sample plane to the output pupil plane with coordinates $\mathbf {x}_o$ using an objective lens with focal length $f_o$. The form of the Green’s function for this mapping is independent of whether we collect forward or backward scattered light. This output field is given by

(8)$$E_2(\mathbf{x}_o) = P_o(\mathbf{x}_o) \, \int \, \chi^{(2)}(\mathbf{r}) \, G(\mathbf{r},\mathbf{u}_i) \, \exp\left( - i \, 2 \, \pi \, \frac{\kappa_2}{f_o} \, \mathbf{r} \cdot \mathbf{x}_o \right) \, d^2 \mathbf{r},$$

where again we have suppressed constants of proportionality for brevity. Identifying the output pupil spatial frequency at the second harmonic optical frequency as $\mathbf {u}_o = \mathbf {x}_o/(\lambda _2 \, f_o)$, then we define the output Green’s function for the SHG field as

(9)$$H(\mathbf{u}_o,\mathbf{r}) = P_o(\mathbf{u}_o) \, e^{- i \, 2 \, \pi \, \mathbf{u}_o \cdot {{\mathbf{r}}}}.$$

B. Scattered SHG field operators

To establish the scattering field operators, we apply Eq. (1) to the Green’s function derived in the previous section. Using the explicit form of the Green’s functions given in Appendix A, we compute the scattering matrix in transmission and reflection in continuous operator form, for which we find identical expressions.

The set of scattered fields that are mapped to the output pupil, $\mathbf {u}_o$, as a function of the input spatial frequency define the SHG reflection operator $R(\mathbf {u}_o,\mathbf {u}_i)$ that is also identical to the transmission operator $T(\mathbf {u}_o,\mathbf {u}_i)$. Inserting Eq. (7) and Eq. (9) into Eq. (1) leads to the integral definition of the SHG scattering operator in transmission

(10)$$R(\mathbf{u}_o,\mathbf{u}_i) \equiv \int P_i(\mathbf{u}_i) \, e^{- i \, 2 \, \pi \, \mathbf{u}_i \cdot {{\mathbf{r}}}} \, \chi^{(2)}(\mathbf{r}) \, P_o(\mathbf{u}_o) \, e^{- i \, 2 \, \pi \, \mathbf{u}_o \cdot {{\mathbf{r}}}} \, d^2 {{\mathbf{r}}} = T(\mathbf{u}_o,\mathbf{u}_i).$$

Defining the scattering vector in transmission as ${{\mathbf{q}}} = \mathbf {u}_o + \mathbf {u}_i$ allows for a compact representation of the scattering operator as

(11)$$R(\mathbf{u}_o,\mathbf{u}_i) = P_o(\mathbf{u}_o) \, \hat{\chi}^{(2)}({{\mathbf{q}}}) \, P_i(\mathbf{u}_i),$$

where the spatial frequency spectral distribution of the nonlinear susceptibility is $\hat {\chi }^{(2)}({{\mathbf{q}}}) = {{\mathcal {F} \{\chi ^{(2)}(\mathbf{r})\}}}({{\mathbf{q}}})$. Here we have defined the Fourier transform as

$$ {{\mathcal{F} \{f({{\mathbf{x}}})\}}}({{\mathbf{u}}}) = \int_{-\infty}^\infty \, f({{\mathbf{x}}}) \, e^{- i \, 2 \, \pi \, {{\mathbf{u}}} \cdot {{\mathbf{x}}}} \, d^2 {{\mathbf{x}}}. $$

The corresponding inverse transform is given by

$$ {{\mathcal{F}^{{-}1} \{F({{\mathbf{u}}})\}}}({{\mathbf{x}}}) = \int_{-\infty}^\infty \, F({{\mathbf{u}}}) \, e^{ i \, 2 \, \pi \, {{\mathbf{u}}} \cdot {{\mathbf{x}}}} \, d^2 {{\mathbf{u}}}. $$

C. Optimal pupil phase estimation through the singular value decomposition

The estimation and removal of the input and output pupil phases to produce an aberration-free synthetic aperture SHG spectrum can be viewed as a constrained optimization problem to produce an undistorted image. By using the method of Lagrange multipliers to find the optimal correction of the reflection and transmission matrices, we show that the dominant eigenvector of the shifted scattering matrix operators, $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$ and $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$, corresponds to the optimal correction. As shown above, since the structure of the matrices $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}$ and $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_o}$ approximately decouples the input and output pupils, the phase shifting problem problem can be written as a simple matrix operation: $\mathbf {S}_{{{\mathbf{q}}},\mathbf{u}_i}\mathbf {a}=E_{\rm SHG}^{\rm s} ({{\mathbf{q}}})$, where $E_{\rm SHG}^{\rm s} ({{\mathbf{q}}})$ is the reconstructed synthetic aperture spectrum and $\mathbf {a}$ is a unit vector that shifts the phase of each column: $\mathbf {a}=e^{i \, \phi _c(\mathbf {u}_i)}$, with $\phi _c$ being the phase correction. We would like to find $\mathbf {a}$ such that it maximizes the total intensity of $E_{\rm SHG}^{\rm s} ({{\mathbf{q}}})$ with the constraint that $\mathbf {a}^{\dagger } \, \mathbf {a}=1$. When the total intensity is maximum this corresponds to the condition that all the columns (fields) in ${{\mathbf{S}}}$ are in phase. This optimal condition occurs when $\mathbf {a} = P_i(\mathbf{u}_i)^*$, implying that $\phi _c= - \phi _i(\mathbf{u}_i)$, thereby correcting the aberrations imparted by the input pupil. As the total intensity as a function of the vector $\mathbf {a}$ is given by

(12)$$f(\mathbf{a})=[E_{\rm SHG}^{\rm s} ({{\mathbf{q}}})]^{{\dagger}}[E_{\rm SHG}^{\rm s} ({{\mathbf{q}}})]=\mathbf{a}^{{\dagger}} \, {{\mathbf{S}}}^{{\dagger}} {{\mathbf{S}}} \, \mathbf{a},$$

the optimization problem can then be written in the form of a constrained optimization problem,

(13)$$\mathrm{maximize}\, f(\mathbf{a}) \ s.t. \ \mathbf{a}^{{\dagger}} \, \mathbf{a}=1.$$

Using the method of Lagrange multipliers, the maximum or minimum of a function $f$ is the solution to $\nabla f=\lambda \nabla g$ where $g$ is a constraint function, in this case $g(\mathbf {a})=0=\mathbf {a}^{\dagger }\mathbf {a}-1$. Written in a different way, the Lagrangian is $\mathcal {L}=\mathbf {a}^{\dagger } \, {{\mathbf{S}}}^{\dagger } {{\mathbf{S}}} \, \mathbf {a}-\lambda \, (\mathbf {a}^{\dagger }\mathbf {a}-1)$, where $\nabla \mathcal {L}=0$. Since the matrices and vectors here are complex valued, care is needed to properly calculate these derivatives using Wirtinger calculus. Conveniently, the expressions for the derivatives we need are in appendix A of [45]. Taking the Wirtinger derivatives gives the gradient of the Lagrangian as

(14)$$\frac{\partial \mathcal{L}}{\partial \mathbf{a}}=0=({{\mathbf{S}}}^{{\dagger}} {{\mathbf{S}}})^T \, \mathbf{a}^*-\lambda \, \mathbf{a}^*.$$

Rearranging Eq. (14) provides the expression

(15)$$({{\mathbf{S}}}^{{\dagger}} {{\mathbf{S}}})^*\mathbf{a}^*=\lambda \, \mathbf{a}^*.$$

By taking the complex conjugate of Eq. (15), we obtain an the formula

(16)$$({{\mathbf{S}}}^{{\dagger}} {{\mathbf{S}}}) \, \mathbf{a}=\lambda^* \, \mathbf{a},$$

which is an eigenvalue equation where $\mathbf {a}$ is an eigenvector of ${{\mathbf{S}}}^{\dagger } {{\mathbf{S}}}$, with eigenvalue $\lambda ^*$. Since ${{\mathbf{S}}}^{\dagger } {{\mathbf{S}}}$ is a Hermitian matrix, it has real eigenvalues, which implies that $\lambda ^*=\lambda$. It is evident that the eigenvectors of ${{\mathbf{S}}}^{\dagger } {{\mathbf{S}}}$ are the left singular vectors of ${{\mathbf{S}}}$ with $\sqrt {\lambda }$ representing the corresponding singular values. The unit vector $\mathbf {a}$ which maximizes the total intensity of the synthetic aperture image is the left singular vector of ${{\mathbf{S}}}$ corresponding to the largest singular value, which is the dominant singular vector. When the total intensity is maximized, this corresponds to the condition in which each field is added coherently in phase with one another.

D. Comparison of the robustness of phase estimation algorithms

The critical aspect of aberration-free synthetic aperture SHG holographic imaging is to robustly estimate the correct input and output pupil phase and use those to correct the data and estimate an undistorted SHG object spatial frequency spectrum. This process becomes difficult when signal levels are low, which certainly is the case for epi directed SHG from biological samples. Not only is the signal low, but exposure times must be kept short due to the instabilities of the interferometer. Finding and correcting for aberrations amounts to finding phase differences between scattered fields originating from similar input angles. These measurements contain phase information so the phase difference between two fields can be found by taking their cross-correlation. The cross-correlation strategy works well when signal to noise levels are high, but as SNR decreases, noise disrupts this the estimation of the phase from the cross-correlation. Our SVD approach has a distinct advantage as it takes the entire dataset into consideration at once instead of finding phase differences between two neighboring field vectors. As a consequence, we anticipate superior performance for the SVD algorithm under conditions of low SNR. To test the robustness of our new SVD-based algorithm, we have run simulations with varying noise levels to compare the fidelity of estimating the pupil phase with our new SVD algorithm as compared to the cross-correlation algorithm used previously to great effect for linear scattering [39].

In the simulation, a reflection matrix is generated and then a pupil phase distortion is applied. A phase distortion is applied by applying random weights to the first 30 elements of the Zernike basis. Then varying levels of noise were applied to each field so that the noise is uncorrelated from one field to another. The noise was added to the fields in the spatial domain with a uniformly distributed random amplitude and a uniformly distributed random phase from $-\pi$ to $\pi$. The noise level was changed by varying the amplitude. To quantify the error of the phase map reconstruction, the recovered pupil map is first transformed into the spatial domain by an inverse Fourier transform. Then each reconstruction is compared to the actual phase distortion using a normalized mean squared error calculation: $\mathrm {NMSE}=\left [ \Vert x_{\rm ref}-x\Vert \right ]/\left [\Vert x_{\rm ref}-\mathrm {mean}(x_{\rm ref})\Vert \right ]$. Inspection of Fig. 8 indicates more robust performance of theSVD algorithm at low SNR values.

 

Fig. 8. Performance of the SVD algorithm compared to the cross-correlation algorithm for estimating the pupil phase under varying levels of SNR. The actual pupil phase is shown in the top left with a black border. At selected SNR levels the recovered pupil phase maps for each technique are shown. The result using the cross-correlation method is shown with dashed red borders and dashed blue borders for the SVD method.

Download Full Size | PDF

Funding

National Science Foundation (DMS-2006416); Chan Zuckerberg Initiative.

Acknowledgments

We would like to thank Tia Tedford, histo technician at the Orthopaedic and Bioengineering Research Lab at Colorado State University, for her expertise in sample preparation of the sheep tendon tissues that we imaged in this work.

Disclosures

“The authors declare no conflicts of interest.”

Data availability

Data and processing scripts written in Matlab are available at [46,47].

References

1. D. S. James and P. J. Campagnola, “Recent advancements in optical harmonic generation microscopy: Applications and perspectives,” BME Front. 2021, 3973857 (2021). [CrossRef]  

2. F. S. Pavone and P. J. Campagnola, Second Harmonic Generation Imaging (CRC Press, 2014).

3. J. J. Field, K. A. Wernsing, S. R. Domingue, A. M. Allende Motz, K. F. DeLuca, D. H. Levi, J. G. DeLuca, M. D. Young, J. A. Squier, and R. A. Bartels, “Superresolved multiphoton microscopy with spatial frequency-modulated imaging,” Proc. Natl. Acad. Sci. U. S. A. 113(24), 6605–6610 (2016). [CrossRef]  

4. S. Fine and W. Hansen, “Optical second harmonic generation in biological systems,” Appl. Opt. 10(10), 2350–2353 (1971). [CrossRef]  

5. I. Freund and M. Deutsch, “Second-harmonic microscopy of biological tissue,” Opt. Lett. 11(2), 94–96 (1986). [CrossRef]  

6. P. J. Campagnola, A. Lewis, L. M. Loew, et al., “High-resolution nonlinear optical imaging of live cells by second harmonic generation,” Biophys. J. 77(6), 3341–3349 (1999). [CrossRef]  

7. A. Zoumi, A. Yeh, and B. J. Tromberg, “Imaging cells and extracellular matrix in vivo by using second-harmonic generation and two-photon excited fluorescence,” Proc. Natl. Acad. Sci. U. S. A. 99(17), 11014–11019 (2002). [CrossRef]  

8. P. J. Campagnola, A. C. Millard, M. Terasaki, P. E. Hoppe, C. J. Malone, and W. A. Mohler, “Three-dimensional high-resolution second-harmonic generation imaging of endogenous structural proteins in biological tissues,” Biophys. J. 82(1), 493–508 (2002). [CrossRef]  

9. P. J. Campagnola and C.-Y. Dong, “Second harmonic generation microscopy: principles and applications to disease diagnosis,” Laser Photonics Rev. 5, 13–26 (2011). [CrossRef]  

10. R. Ambekar, T.-Y. Lau, M. Walsh, R. Bhargava, and K. C. Toussaint, “Quantifying collagen structure in breast biopsies using second-harmonic generation imaging,” Biomed. Opt. Express 3(9), 2021–2035 (2012). [CrossRef]  

11. P. P. Provenzano, K. W. Eliceiri, J. M. Campbell, D. R. Inman, J. G. White, and P. J. Keely, “Collagen reorganization at the tumor-stromal interface facilitates local invasion,” BMC Med. 4(1), 38 (2006). [CrossRef]  

12. N. Olivier, M. A. Luengo-Oroz, L. Duloquin, E. Faure, T. Savy, I. Veilleux, X. Solinas, D. Débarre, P. Bourgine, A. Santos, N. Peyriéras, and E. Beaurepaire, “Cell lineage reconstruction of early zebrafish embryos using label-free nonlinear microscopy,” Science 329(5994), 967–971 (2010). [CrossRef]  

13. T. Manaka and M. Iwamoto, “Optical second-harmonic generation measurement for probing organic device operation,” Light: Sci. Appl. 5(3), e16040 (2015). [CrossRef]  

14. H. Ma, J. Liang, H. Hong, K. Liu, D. Zou, M. Wu, and K. Liu, “Rich information on 2d materials revealed by optical second harmonic generation,” Nanoscale 12(45), 22891–22903 (2020). [CrossRef]  

15. P. Bianchini and A. Diaspro, “Three-dimensional (3d) backward and forward second harmonic generation (shg) microscopy of biological tissues,” J. Biophotonics 1, 443–450 (2008). [CrossRef]  

16. S. Yazdanfar, L. H. Laiho, and P. T. C. So, “Interferometric second harmonic generation microscopy,” Opt. Express 12(12), 2739–2745 (2004). [CrossRef]  

17. D. Raanan, M. S. Song, W. A. Tisdale, and D. Oron, “Super-resolved second harmonic generation imaging by coherent image scanning microscopy,” Appl. Phys. Lett. 120(7), 071111 (2022). [CrossRef]  

18. S. Liu, M. R. E. Lamont, J. A. Mulligan, and S. G. Adie, “Aberration-diverse optical coherence tomography for suppression of multiple scattering and speckle,” Biomed. Opt. Express 9(10), 4919–4935 (2018). [CrossRef]  

19. M. J. Booth, “Adaptive optics in microscopy,” Phil. Trans. R. Soc. A. 365(1861), 2829–2843 (2007). [CrossRef]  

20. M. J. Booth, “Adaptive optical microscopy: the ongoing quest for a perfect image,” Light: Sci. Appl. 3(4), e165 (2014). [CrossRef]  

21. J. Mertz, H. Paudel, and T. G. Bifano, “Field of view advantage of conjugate adaptive optics in microscopy applications,” Appl. Opt. 54(11), 3498–3506 (2015). [CrossRef]  

22. N. Olivier, D. Débarre, and E. Beaurepaire, “Dynamic aberration correction for multiharmonic microscopy,” Opt. Lett. 34(20), 3145–3147 (2009). [CrossRef]  

23. N. Ji, D. E. Milkie, and E. Betzig, “Adaptive optics via pupil segmentation for high-resolution imaging in biological tissues,” Nat. Methods 7(2), 141–147 (2010). [CrossRef]  

24. D. Nuzhdin, E. G. Pendleton, E. B. Munger, L. J. Mortensen, and S. Brasselet, “In-depth polarisation resolved shg microscopy in biological tissues using iterative wavefront optimisation,” J. Microsc. 291(1), 57–72 (2023). [CrossRef]  

25. O. Masihzadeh, P. Schlup, and R. A. Bartels, “Label-free second harmonic generation holographic microscopy of biological specimens,” Opt. Express 18(10), 9840–9851 (2010). [CrossRef]  

26. E. Shaffer, C. Moratal, P. Magistretti, P. Marquet, and C. Depeursinge, “Label-free second-harmonic phase imaging of biological specimen by digital holographic microscopy,” Opt. Lett. 35(24), 4102–4104 (2010). [CrossRef]  

27. E. Shaffer, P. Marquet, and C. Depeursinge, “Real time, nanometric 3d-tracking of nanoparticles made possible by second harmonic generation digital holographic microscopy,” Opt. Express 18(16), 17392–17403 (2010). [CrossRef]  

28. D. R. Smith, D. G. Winters, P. Schlup, and R. A. Bartels, “Hilbert reconstruction of phase-shifted second-harmonic holographic images,” Opt. Lett. 37(11), 2052–2054 (2012). [CrossRef]  

29. D. G. Winters, D. R. Smith, P. Schlup, and R. A. Bartels, “Measurement of orientation and susceptibility ratios using a polarization-resolved second-harmonic generation holographic microscope,” Biomed. Opt. Express 3(9), 2004–2011 (2012). [CrossRef]  

30. D. R. Smith, D. G. Winters, and R. A. Bartels, “Submillisecond second harmonic holographic imaging of biological specimens in three dimensions,” Proc. Natl. Acad. Sci. U. S. A. 110(46), 18391–18396 (2013). [CrossRef]  

31. C. Hu, J. J. Field, V. Kelkar, B. Chiang, K. Wernsing, K. C. Toussaint, R. A. Bartels, and G. Popescu, “Harmonic optical tomography of nonlinear structures,” Nat. Photonics 14(9), 564–569 (2020). [CrossRef]  

32. L. Vittadello, C. Kijatkin, J. Klenen, D. Dzikonski, K. Kömpe, C. Meyer, A. Paululat, and M. Imlau, “In-vivo tracking of harmonic nanoparticles: a study based on a tiger widefield microscope invited,” Opt. Mater. Express 11(7), 1953–1969 (2021). [CrossRef]  

33. R. M. Williams, W. R. Zipfel, and W. W. Webb, “Interpreting second-harmonic generation images of collagen i fibrils,” Biophys. J. 88(2), 1377–1386 (2005). [CrossRef]  

34. R. LaComb, O. Nadiarnykh, S. S. Townsend, and P. J. Campagnola, “Phase matching considerations in second harmonic generation from tissues: effects on emission directionality, conversion efficiency and observed morphology,” Opt. Commun. 281(7), 1823–1832 (2008). [CrossRef]  

35. O. Nadiarnykh, R. B. LaComb, M. A. Brewer, and P. J. Campagnola, “Alterations of the extracellular matrix in ovarian cancer studied by second harmonic generation imaging microscopy,” BMC Cancer 10(1), 94 (2010). [CrossRef]  

36. S. A. Alexandrov, T. R. Hillman, T. Gutzler, and D. D. Sampson, “Synthetic aperture fourier holographic optical microscopy,” Phys. Rev. Lett. 97(16), 168102 (2006). [CrossRef]  

37. A. E. Tippie, A. Kumar, and J. R. Fienup, “High-resolution synthetic-aperture digital holography with digital phase and pupil correction,” Opt. Express 19(13), 12027–12038 (2011). [CrossRef]  

38. P. Rajaeipour, A. Dorn, K. Banerjee, H. Zappe, and Ç. Ataman, “Extended field-of-view adaptive optics in microscopy via numerical field segmentation,” Appl. Opt. 59(12), 3784–3791 (2020). [CrossRef]  

39. S. Kang, P. Kang, S. Jeong, Y. Kwon, T. D. Yang, J. H. Hong, M. Kim, K.-D. Song, J. H. Park, J. H. Lee, M. J. Kim, K. H. Kim, and W. Choi, “High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering,” Nat. Commun. 8(1), 2157 (2017). [CrossRef]  

40. A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, “Smart optical coherence tomography for ultra-deep imaging through highly scattering media,” Sci. Adv. 2(11), e1600370 (2016). [CrossRef]  

41. A. Badon, V. Barolle, K. Irsch, A. C. Boccara, M. Fink, and A. Aubry, “Distortion matrix concept for deep optical imaging in scattering media,” Sci. Adv. 6(30), eaay7170 (2020). [CrossRef]  

42. M. Gavish and D. L. Donoho, “The optimal hard threshold for singular values is $4/\sqrt {3}$,” IEEE Trans. Inform. Theory 60(8), 5040–5053 (2014). [CrossRef]  

43. A. Kanaev, A. Watnik, D. Gardner, C. Metzler, K. Judd, P. Lebow, K. Novak, and J. Lindle, “Imaging through extreme scattering in extended dynamic media,” Opt. Lett. 43(13), 3088–3091 (2018). [CrossRef]  

44. J. Mertz, Introduction to Optical Microscopy (Cambridge University Press, 2019), 2nd ed.

45. R. F. H. Fischer, Precoding and Signal Shaping for Digital Transmission (Wiley-IEEE Press, 2002).

46. N. Ji, T. R. Sato, and E. Betzig, “Characterization and adaptive optical correction of aberrations during in vivo imaging in the mouse cortex,” Proc. Natl. Acad. Sci. U. S. A. 109(1), 22–27 (2012). [CrossRef]  

47. R. A. Bartels, “SHGSyntheticApertureHologrpahy,” GitHub (2023), https://github.com/RandyBartelsCSU/SHGSyntheticApertureHologrpahy.