## Abstract

Large-area multiphoton laser scanning microscopy (LMLSM) can be applied in biology and medicine for high sensitivity and resolution tissue imaging. However, factors such as refractive index mismatch induced spherical aberration, emission/excitation absorption and scattering can result in axial intensity attenuation and lateral image heterogeneity, affecting both qualitative and quantitative image analysis. In this work, we describe an image correction algorithm to improve three-dimensional images in LMLSM. The method consists of multiplying the measured nonlinear signal by a three-dimensional correction factor, determined by the use of two-photon images of the appropriate specimens and specimen absorption and scattering properties at the excitation and emission wavelengths. The proposed methodology is demonstrated in correcting multiphoton images of objects imbedded in uniform fluorescent background, lung tissue, and *Drosophila* larva.

© 2008 Optical Society of America

## 1. Introduction

In recent years, multiphoton laser scanning microscopy (MLSM) has evolved into an effective imaging modality applicable to many areas of biology and medicine. The high resolution, superior axial depth discrimination capability, and minimally invasive nature of this technique renders MLSM to be the preferred method for three-dimension imaging of thick tissue and in live animals [1]. However, images acquired by MLSM can demonstrate intensity heterogeneity in both the lateral and axial directions, which can hinder accurate image visualization and analysis. As a result, both qualitative presentation and quantitative description of multiphoton results can be adversely affected. To be specific, the nonlinear nature of the multiphoton excitation process makes the signal particularly sensitive to optical effects such as the refractive index mismatch induced spherical aberration and the absorption and scattering of the excitation and emission photons [2–6]. These effects can contribute to image quality degradation and intensity decay, hindering the qualitative presentation of the images and extraction of quantitative results difficult. In confocal microscopy, various software methods have been developed to correct the intensity attenuation in reconstructed images [7–11] and hardware approaches were suggested for image correction in both confocal and multiphoton images [12–15].

While the image inhomogeneity artifacts are known, their effects on image quality and quantification are not apparent in standard imaging applications where the image of the area under examination is composed of a single scan with an area of approximately 100 by 100 µm2. With recent advances in optical microscopy, the global 2-D or 3-D architecture of a tissue specimen can only be revealed by high resolution, large-area examination across the tissue sample. In large-area MLSM (LMLSM), a series of adjacent, high resolution, small-area multiphoton scans are acquired and assembled for tissue visualization on a global scale. This approach has been successfully applied to address a wide array of biomedical problems including transdermal drug delivery, whole-organ cardiac imaging, tissue engineering, and the global understanding of cornea structures. In addition, it has been demonstrated that pathological conditions such as keratoconus, basal cell carcinoma, and corneal infection can be better diagnosed by the large-area imaging methodology [16–23]. However, in many of these studies, imaging field inhomogeneity contributes to periodic intensity variations that lead to artifacts affecting both the qualitative presentation and quantitative analysis of the images.

In this work, we address the problem of intensity inhomogeneity in LMLSM by the use a computational algorithm. Our methodology is effective for the correction of axial intensity profiles and the homogenization of the image array and that this method may be extended to other optical imaging applications where image stitching is necessary.

## 2. Materials and methods

Experiments on TPF intensity degradation along the lateral and axial coordinates were conducted on a LSM 510 META system coupled with a femtosecond titanium:sapphire (Tsunami, Spectra-Physics, Mountain View, CA) laser operating at 780 nm. Multiphoton images were acquired using three high numerical aperture objectives: the water immersion C-Apochromat 40×/NA 1.2W Corr objective, the Fluar 40×/NA 1.3 Oil objective, and the Plan-Apochromat 63×/NA 1.4 Oil DIC objective. The uniform fluorescent samples include low (10 µM) and high (0.28 mM) concentration of sulforhodamine B (SRB) in water. In addition, a 20 µM SRB solution containing 5% Intrafat (Nihon Pharmaceutical Co., Osaka, Japan) was used as a scattering fluorescent sample. The absorption and scattering properties of the fluorescent samples 0.49 mm in thickness were measured by a spectrophotometer (DU-800, Beckman Coulter). To demonstrate the effectiveness of our methodology in correcting the image field inhomogeneity artifact, multiphoton autofluorescence and/or second harmonic generation (SHG) images of uniformly fluorescent solution, *ex-vivo* human lung tissue and stage 3 of a *Drosophila* larva were acquired and analyzed. The lung tissue specimen was obtained from the tissue bank of National Taiwan University Hospital. For imaging purposes, a thin slice of the lung tissue block (approximately 0.4 mm in thickness) was excised and mounted on a glass slide and covered with a coverglass. The approximate imaging depth for the lung tissue specimen was 30 µm.

## 3. Results

Shown in Fig. 1 is the 3×3 image array illustrating TPF lateral intensity inhomogeneity and axial intensity profiles of the low concentration SRB sample. (Curves 1-5 taken with the Plan-Apochromat 63×/NA 1.4 Oil objective). For comparison, the axial intensity profile acquired using the water immersion objective (C-Apochromat 40×/NA 1.2W Corr) was shown (Curve 6). As Fig. 1 shows, the axial intensity profile of the water objective is mostly depth-independent and attenuation of the detected signal fluorescence photons is less than 3% at a depth of 80 µm. On the other hand, the attenuation of TPF signal (~50%) for the oil-immersion objective was caused by refractive index mismatch induced spherical aberration. The lateral profiles along the x (horizontal) and y (vertical) axes through lines passing through the point of maximal TPF intensity (Region 1) at different imaging depths are shown in Fig. 2. Analysis shows that, in general, the image lateral profiles have a convex and asymmetric shape. We found that the lack of intensity uniformity and symmetry within the scanning field is commonplace, and is due to imaging the specimen away from the objective’s optical axis.

#### 3.1 Three-dimensional image correction model

The majority of intensity correction computational methods are developed for confocal microscopy, where absorption (abs) and scattering (scat) are the most essential signal degradation factors. Exponential models describing signal attenuation have been used for image correction away from the focal point [24]:

In Eq. (1), *I*(*z*,λ) is the collected signal, which is a function of the focal position *z* and the wavelength λ, the emission attenuation coefficient (α_{em}=α_{em,abs}+α_{em,scat}), and the excitation attenuation coefficient (α_{exc}=α_{exc,abs}+α_{exc,scat}). It was shown that, in general, the loss of signal intensity was better described by a bi-exponential model [7].

The original image intensity *Î*(*x̄*) can be recovered by multiplying the observed intensity *I*(*x̄*) by a spatial correction constant*C*(*x̄*),

where

is the inverse of overall attenuation coefficient taking into account the effects due to excitation (γ_{exc}) and emission (γ_{em}) attenuation, and *x̄* is spatial coordinate of the focal point [7,12].

Previously, it was demonstrated that for imaging neuron specimens, the attenuation due to spherical aberration is the most significant factor in signal degradation. [12] It was proposed that a spherical aberration factor γ_{sph} (*z*) be included in Eq. (3). For convenience, an exponential model of signal attenuation as a function of imaging depth was used

For obtaining γ_{exc} and γ_{em}, the following equations have been used:

where ω is objective semi-aperture angle and θ, φ are the angular coordinates in the spherical coordinate system.

However, this model and the associated algorithm developed for image reconstruction are less applicable for multiphoton image correction. First of all, in the case of two-photon signal acquisition, Eq. (1) is converted to

[24]. The corresponding change of α_{exc}→2α_{exc} is due to the quadratic dependence of the signal on excitation intensity. In addition, multiphoton induced signal depends on spatial and temporal distribution of intensity at the focal point and the spherical aberration axial factor cannot be described by a simple exponential function. In general, it is complex in form and depends strongly on the microscope objective used [4–7].

In this work, we propose a method of image field correction in nonlinear microscopy by independent determination of the degradation factors. Specifically, for axial intensity profile correction, the experimental TPF axial intensity profile acquired from a homogenous, weakly absorbing and scattering medium with the refractive index close to that of the sample can be used as the axial spherical aberration correction factor. In addition, the spectrophotometrically measured mean absorption and scattering coefficients of the sample at the emission and excitation wavelengths are needed to determine the γ_{em} and γ_{exc} attenuation factors. For lateral image correction, a lateral intensity inhomogeneity correction factor, *L _{lat}*(

*x̄*), needs to be determined in the correction function (Eq. (3)).

The final correction factor has following form:

The lateral inhomogeneity correction factor *L _{lat}*(

*x̄*) can be estimated using a model of homogenous sample or can be acquired by averaging lateral profiles along both the x and y coordinates.

#### 3.2 TPF intensity axial attenuation correction

To demonstrate our approach for the axial attenuation correction, we used samples with low absorption and scattering coefficients (Sample I: 10 µM SRB aqueous solution), high absorption and low scattering coefficients (Sample II: 0.28 mM SRB aqueous solution), and high scattering coefficient (Sample III: 20 µM SRB in Intrafat-water emulsion). The spectral dependences of specimen optical density (OD) for Samples II and III are shown in Fig. 3. Using these curves and the fluorescence spectrum of SRB, we estimated the attenuation coefficients for emission and excitation to be α_{em,abs}=17 cm^{-1} and α_{exc}=0 for Sample II, and α_{em}=α_{em,abs}+α_{em,scat}=42 cm^{-1} and α_{exc,scat}=20 cm^{-1} for Sample III. Estimates of the attenuation coefficients were obtained by measuring the specimens’ absorption spectra and normalizing to sulforhodamine B’s fluorescence emission spectrum from a previously published work. [25] This was done because our laboratory is not equipped with a standard fluorometer. In this way, the spectrally integrated attenuation coefficients can be obtained.

These data were then used to reconstruct the TPF intensity axial profiles acquired using the two oil-immersion objectives. As Figs. 1 and 4 shows, refractive index mismatch induced spherical aberration causes significant TPF intensity attenuation for these objectives (more than 50 % at the depth of 80 µm depth) and that the axial intensity decay cannot be fitted by exponential functions. In Fig. 4, the TPF intensity axial profiles of SRB samples I and II acquired using the Fluar 40×/NA 1.3 objective are shown. Specifically, Curve 3 represents the emission attenuation factor obtained using the attenuation coefficient α_{em}=17 cm^{-1}, and Curve 4 is the corrected TPF intensity profile. For TPF intensity reconstruction by Eq. (8), experimental data series (Curve 1, Fig. 4) was used as the spherical aberration factor-*P*
_{sph}(*z*). In addition, since Sample II is an aqueous SRB specimen with no scattering components and negligible linear absorption at the excitation wavelength of 780 nm (Curve 1, Fig. 3), the excitation attenuation factor γ_{exc}(*z*) is 1 (from α_{exc}=0) for this specimen. Between the depths of 0 and 80 µm, the corrected data agree within 11% to the uniform profile of the 0.28 mM SRB specimen (Curve 4 Fig. 4).

In Fig. 5, measured and calculated TPF axial intensity profiles for the highly absorbing and highly scattering samples (Plan-Apochromat 63×/NA 1.4 Oil objective) are shown. In conjunction with Eqs. (5) and (6), the values of α_{em}=α_{em,abs}=17 cm^{-1} (Sample II, Curve 2), and α_{em}=42 cm^{-1} and α_{exc}=2α_{exc,scat}=40 cm^{-1} (Sample III, Curve 5) were used to calculate TPF axial intensity attenuation factors due to absorption and scattering. These results, combined with the spherical aberration factor *P*
_{sph}(*z*) (Curve 1, Fig. 1), were used to correct the TPF axial intensity profiles and the reconstructed data for the both samples are shown in Fig. 5 as Curves 3 and 6. For the Plan-Apochromat 63×/NA 1.4 Oil objective, the corrected data agree within 13% between the imaging depths of 0-80 µm.

Therefore, our results demonstrate that axial intensity degradation can be reasonably well corrected.

#### 3.3 Lateral image correction

To obtain the lateral correction factor, TPF intensity lateral profiles of the SRB uniform samples are used (Figs. 2 and 6(A)). It can be seen that the lateral profiles are convex in form. In general, these profiles are asymmetric with respect to the optical axis. We found that one of the simplest analytic curves for fitting of these data series along each of x and y coordinates can be formed by two parabolas, with peaks coinciding with the point of TPF intensity maximum *I*
_{max}(*x*
_{max}, *y*
_{max}) (Fig. 6(A)). In addition, one of the parabolas passes through the point *I*(0)_{i} and second parabola passes through the point *I*(*l*)_{i} where *i* is *x* or *y* coordinate and *I*(0)_{i} and *I*(l)_{i} are the TPF intensity at coordinates *i*=0 and *i*=*l* respectively. Note that *l* is size of the image.

In the first approximation, it is assumed that

where, according to our two-parabola model,

where
${B}_{i}=\sqrt{{I}_{max}-{I}_{i}\left(0\right)}$
when *i*≤*i*
_{max} and
${B}_{i}=\frac{\sqrt{{I}_{max}-I{\left(l\right)}_{i}}}{\left(l-{i}_{max}\right)}$
when *i*
_{max}<*i*≤. Furthermore, *C*
_{i}=*I*
_{max} and *A*
_{i}=*B*
_{i}/*i*
_{max}. In our approach, the two-parabola model described by Eq. (9) can be iteratively applied. The iterative procedure is stopped when the intensity variation between the region of maximum intensity and that of the image edges is less than 10%.

The x-y image of correction factor for the profiles in Fig. 6(A) calculated using Eq. (9) is shown in Fig. 6(B). If needed, the lateral correction factor *L _{lat}*(

*x̄*) can be constructed at every axial coordinate.

The procedures outlined above were demonstrated by correction of TPF image of homogenous SRB sample with two bubbles (Fig. 7.), TPF/SHG images of human lung tissue (Fig. 8), and TPF/SHG image of a Stage 3 *Drosophila larva* (Stage 3) (Fig. 9). These images were acquired using the Fluar 40×/NA 1.3 Oil objective.

To evaluate the efficiency of the lateral correction method, mean value (*I*
_{m}) and standard deviation (*StD*) for the 5×5 TPF image array from Fig. 7 and small uniform area of 40×40 µm2 in the first tile center before and after lateral correction were calculated and the results are listed in the Table 1.

These data show that the standard deviation of the corrected image is close to that of uniform area.

For lateral image correction of the lung specimen in Fig. 8 and the Stage 3 Drosophila larva in Fig. 9, averaged lateral profiles of all images along horizontal and vertical directions were used to determine parabolas’ parameters in Eq. (9). The lateral intensity correction profiles for the TPF and SHG channels were independently determined and applied for image correction.

In addition to the two-parabola approach, other mathematical models such as polynomial functions may be used. Since our method is, in general, iterative, implementing an additional iteration is equivalent to increase the previous polynomial by two orders. In that sense, our approach is polynomial in nature. However, our approach is advantageous in that we do not need to predetermine the orders of polynomial to use. Furthermore, we found that the two-parabola model works well in correcting for the intensity variation across the imaging field. Finally, since our approach is computational in nature, minor mismatches at the image boundaries are expected to occur.

## 4. Discussion and Conclusion

Image distorting artifacts along the lateral and axial axes is a well-recognized problem in both confocal and multiphoton microscopy. The axial intensity decay and lateral image inhomogeneity can be due to factors such as absorption and scattering of excitation and emission photons, spherical aberration. Development of image correction tools in three-dimensional multiphoton imaging is needed for accurate image presentation.

In this work, a methododology is developed to compensate for intensity inhomogeneity in TPF and SHG microscopy. For correction of nonlinear signal intensity axial attenuation, the aberration, absorption and scattering factors were determined. For that the mean absorption and scattering attenuation coefficients of the specimen were measured at excitation and emission wavelengths. The corresponding attenuation factors were calculated in view of the nonlinear nature of signal generation. For estimation of spherical aberration contribution, TPF intensity axial profile was recorded using the same objective and a low absorbing and low scattering sample with refractive index close to that of the examined specimens. Finally, a lateral correction factor was determined using the lateral TPF intensity profiles within the image tile borders of a homogeneous sample. In the case of heterogeneous specimens (lung tissue and *Drosophila* larva), the lateral correction factor can be estimated by averaging of the lateral image profiles.

The proposed method has obvious shortcomings in the highly scattering samples or in-depth imaging applications where multiple scattering events can take place. In such cases, the measurement of scattering and absorption coefficients may be achieved by the use of thinner samples or integrating spheres can be used. Furthermore, the calculation of the attenuation factors when α_{scat}
*z*>1 can be performed using Monte Carlo methods. However, in correcting the image field inhomogeneity, since we obtain the correction function by averaging over many sample images acquired at the same depth, our approach is immune to scattering and refractive mismatch induced spherical aberration.

Our methodology allows the separation and analysis of the contribution of individual degradation factors. The suggested approach may be used for automated image analysis and helps to improve both the qualitative and quantitative presentations of multiphoton imaging results.

## Acknowledgment

This work is supported by the National Research Program for Genomic Medicine (NRPGM) of the National Science Council in Taiwan and was completed in the Optical Molecular Imaging Microscopy Core Facility (A5) of NRPGM. We also like to thank the tissue bank of National Taiwan University Hospital for providing the lung tissue specimen.

## References and links

**1. **W. R. Zipfel, R. M. Williams, and W. W. Webb, “Nonlinear magic: multiphoton microscopy in the biosciences,” Nat. Biotechnol. **21**, 1369–1377 (2003). [CrossRef] [PubMed]

**2. **G. J. de Grauw, J. M. Vroom, H. T. M. van der Voort, and H. C. Gerritsen, “Imaging properties in two-photon excitati-on microscopy and effects of refractive-index mismatch in thick specimens,” Appl. Opt. **38**, 5995–6003 (1999). [CrossRef]

**3. **H. C. Gerritsen and C. J. de Grauw, “Imaging of optically thick specimen using two-photon excitation microscopy,” Microsc. Res. Tech. **47**, 206–209 (1999). [CrossRef] [PubMed]

**4. **A. K. Dunn, V. P. Wallace, M. Coleno, M. W. Berns, and B. J. Tromberg, “Influence of optical properties on two-photon fluorescence imaging in turbid samples,” Appl. Opt. **39**, 1194–1201 (2000). [CrossRef]

**5. **C. Y. Dong, K. Koenig, and P. So, “Characterizing point spread functions of two-photon fluorescence microscopy in turbid medium,” J. Biomed. Opt. **8**, 450–459 (2003). [CrossRef] [PubMed]

**6. **C. Y. Dong, B. Yu, P. D. Kaplan, and P. T. C. So, “Performances of high numerical aperture water and oil immersion objective in deep-tissue, multi-photon microscopic imaging of excised human skin,” Microsc. Res. Tech. **63**, 81–86 (2004). [CrossRef]

**7. **T. Visser, F. Groen, and G. Brakenhoff, “Absorption and scattering correction in fluorescence confocal microscopy,” J. Microsc. **163**, 189–200 (1991). [CrossRef]

**8. **J. Roerdink and M. Bakker, “An FFT-based method for attenuation correction in fluorescence confocal microscopy,” J. Microsc. **169**, 3–14 (1993). [CrossRef]

**9. **J. Markham and J. Conchello, “Artifacts in restored images due to intensity loss in 3D fluorescence microscopy,” J. Microsc. **204**, 93–98 (2001). [CrossRef] [PubMed]

**10. **C. Kervrann, D. Legland, and L. Pardini, “Robust incremental compensation of the light attenuation with depth in 3D fluorescence microscopy,” J. Microsc. **214**, 297–314 (2004). [CrossRef] [PubMed]

**11. **S. C. Lee and P. Bajcsy, “Intensity correction of fluorescent confocal laser scanning microscope images by mean-weight filtering,” J. Microsc. **221**, 122–136 (2006). [CrossRef] [PubMed]

**12. **A. Can, O. Al-Kofahi, S. Lasek, D. H. Szaworski, J. N. Turner, and B. Roysam, “Attenuation correction in confocal laser microscopes:a novel two-view approach,” J. Microsc. **211**, 67–79 (2003). [CrossRef] [PubMed]

**13. **M. Schwertner, M. J. Booth, and T. Wilson, “Characterizing specimen induced aberrations for high NA adaptive optical microscopy,” Opt. Express **12**, 6540–6552 (2004). [CrossRef] [PubMed]

**14. **W. Lo, Y. Sun, S. J. Lin, S. H. Jee, and C. Y. Dong, “Spherical aberration correction in multiphoton fluorescence imaging using objective correction collar,” J. Biomed. Opt. **10**, 034006 (2005). [CrossRef] [PubMed]

**15. **M. Rueckel, J. A. Mack-Bucher, and W. Denk “Adaptive wavefront correction in two-photon microscopy using coherence-gated wavefront sensing,” Proc. Nat. Acad. Sci. U.S.A. **103**, 17137–17142 (2006). [CrossRef]

**16. **B. Yu, K. H. Kim, P. T. C. So, D. Blankschtein, and R. Langer, “Topographic heterogeneity in transdermal transport revealed by high-speed two-photon microscopy: Determination of representative skin sample sizes,” J. Invest. Dermatol. **118**, 1085–1088 (2002). [CrossRef] [PubMed]

**17. **T. Ragan, J. D. Sylvan, K. H. Kim, H. Huang, K. Bahlmann, R. T. Lee, and P. T. C. So, “High-resolution whole organ imaging using two-photon tissue cytometry,” J. Biomed. Opt. **12**, 014015 (2007). [CrossRef] [PubMed]

**18. **S. W. Teng, H. Y. Tan, J. L. Peng, H. H. Lin, K. H. Kim, W. Lo, Y. Sun, W. C. Lin, S. J. Lin, S. H. Jee, P. T. C. So, and C. Y. Dong, “Multiphoton autofluorescence and second-harmonic generation imaging of the ex vivo porcine eye,” Invest. Ophth. Vis. Sci. **47**, 1216–1224 (2006). [CrossRef]

**19. **W. Lo, S. W. Teng, H. Y. Tan, K. H. Kim, H. C. Chen, H. S. Lee, Y. F. Chen, P. T. C. So, and C. Y. Dong, “Intact corneal stroma visualization of GFP mouse revealed by multiphoton imaging,” Microsc. Res. Tech. **69**, 973–975 (2006). [CrossRef] [PubMed]

**20. **H. Y. Tan, Y. Sun, W. Lo, S. J. Lin, C. H. Hsiao, Y. F. Chen, S. C. M. Huang, W. C. Lin, S. H. Jee, H. S. Yu, and C. Y. Dong, “Multiphoton fluorescence and second harmonic generation imaging of the structural alterations in keratoconus ex vivo,” Invest. Ophth. Vis. Sci. **47**, 5251–5259 (2006). [CrossRef]

**21. **H. Y. Tan, Y. Sun, W. Lo, S. W. Teng, R. J. Wu, S. H. Jee, W. C. Lin, C. H. Hsiao, H. C. Lin, Y. F. Chen, D. H. K. Ma, S. C. M. Huang, S. J. Lin, and C. Y. Dong, “Multiphoton fluorescence and second harmonic generation microscopy for imaging infectious keratitis,” J. Biomed. Opt. **12**, 024013 (2007). [CrossRef] [PubMed]

**22. **S. J. Lin, S. H. Jee, C. J. Kuo, R. J. Wu, W. C. Lin, J. S. Chen, Y. H. Liao, C. J. Hsu, T. F. Tsai, Y. F. Chen, and C. Y. Dong, “Discrimination of basal cell carcinoma from normal dermal stroma by quantitative multiphoton imaging,” Opt. Lett. **31**, 2756–2758 (2006). [CrossRef] [PubMed]

**23. **H. S. Lee, S. W. Teng, H. C. Chen, W. Lo, Y. Sun, T. Y. Lin, L. L. Chiou, C. C. Jiang, and C. Y. Dong, “Imaging human bone marrow stem cell morphogenesis in polyglycolic acid scaffold by multiphoton microscopy,” Tissue Eng. **12**, 2835–2841 (2006). [CrossRef]

**24. **A. Diaspro, *Confocal and two-photon microscopy. Foundations, applications, and advances* (Wiley-Liss, Inc., New York, 2002).

**25. **D. S. dos Santos, Jr. and R. F. Aroca, “Selective surface-enhanced fluorescence and dye aggregation with layer-by-layer film substrates,” Analyst , **132**, 450–454 (2007). [CrossRef]