Image heterogeneity correction in large-area , three-dimensional multiphoton microscopy

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 twophoton 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 OCIS codes: (180.6900) Three-dimensional microscopy; (170.5810) Scanning microscopy; (190.4180) Multiphoton processes. 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). 2. G. J. de Grauw, J. M. Vroom, H. T. M. van der Voort and H. C. Gerritsen, “Imaging properties in twophoton excitati-on microscopy and effects of refractive-index mismatch in thick specimens,” Appl. Opt. 38, 5995–6003 (1999). 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). 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). 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). 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). 7. T. Visser, F. Groen and G. Brakenhoff, “Absorption and scattering correction in fluorescence confocal microscopy,” J. Microsc. 163, 189–200 (1991). 8. J. Roerdink and M. Bakker, “An FFT-based method for attenuation correction in fluorescence confocal microscopy,” J. Microsc. 169, 3-14 (1993). 9. J. Markham and J. Conchello, “Artifacts in restored images due to intensity loss in 3D fluorescence microscopy,” J. Microsc. 204, 93–98 (2001). 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). 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). 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). 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). 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). #89103 $15.00 USD Received 26 Oct 2007; revised 27 Dec 2007; accepted 3 Jan 2008; published 28 Mar 2008 (C) 2008 OSA 31 March 2008 / Vol. 16, No. 7 / OPTICS EXPRESS 5107 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). 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). 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). 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). 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). 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). 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). 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). 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). 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).


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][3][4][5][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][8][9][10][11] and hardware approaches were suggested for image correction in both confocal and multiphoton images [12][13][14][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 μm 2 .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, smallarea 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][17][18][19][20][21][22][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.

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 40x/NA 1.2W Corr objective, the Fluar 40x/NA 1.3 Oil objective, and the Plan-Apochromat 63x/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.

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 63x/NA 1.4 Oil objective).For comparison, the axial intensity profile acquired using the water immersion objective (C-Apochromat 40x/NA 1.2W Corr) was shown (Curve 6).As Fig. 1 shows, the axial intensity profile of the water objective is mostly depthindependent 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 oilimmersion 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.

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 ˆ( ) I x can be recovered by multiplying the observed intensity 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 ( ) exp( ) 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][5][6][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, ( ) lat L x , needs to be determined in the correction function (Eq.( 3)).
The final correction factor has following form: The lateral inhomogeneity correction factor ( ) lat L 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.

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 40x/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 63x/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 63x/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.

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, when i≤i max and max max ( ) ( ) 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 ( ) lat L 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 40x/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 μm 2 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.

A B
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 twoparabola 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.

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 threedimensional 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 indepth 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.

Fig. 1 .
Fig. 1.A 3×3 TPF image array of the uniform 10 μM SRB sample and axial intensity profiles at the point of intensity maximum (Curve 1) and at image edges (Curves 2-5) acquired using the Plan-Apochromat 63x/NA 1.4 Oil DIC objective.For comparison, the axial intensity profile of the C-Apochromat 40x/NA 1.2W Corr water lens acquired at the image field intensity maximum is shown (Curve 6).

Fig. 6 .
Fig. 6. (A) Uniform SRB sample TPF intensity X and Y profiles and fitting results using the two-parabola model.(B) x-y image of lateral correction factor, obtained from (A). (Fluar 40x/NA 1.3 Oil objective).

Table 1 .
Mean value (I m ) and standard deviation (StD) in different regions of SRB sample (Fig.7).