## A Parallel Product-Convolution approach for representing depth varying Point Spread Functions in 3D widefield microscopy based on principal component analysis

Optics Express, Vol. 18, Issue 7, pp. 6461-6476 (2010)

http://dx.doi.org/10.1364/OE.18.006461

Acrobat PDF (8785 KB)

### Abstract

We address the problem of computational representation of image formation in 3D widefield fluorescence microscopy with depth varying spherical aberrations. We first represent 3D depth-dependent point spread functions (PSFs) as a weighted sum of basis functions that are obtained by principal component analysis (PCA) of experimental data. This representation is then used to derive an approximating structure that compactly expresses the depth *variant* response as a sum of few depth *invariant* convolutions pre-multiplied by a set of 1D depth functions, where the convolving functions are the PCA-derived basis functions. The model offers an efficient and convenient trade-off between complexity and accuracy. For a given number of approximating PSFs, the proposed method results in a much better accuracy than the strata based approximation scheme that is currently used in the literature. In addition to yielding better accuracy, the proposed methods automatically eliminate the noise in the measured PSFs.

© 2010 Optical Society of America

## 1. Introduction

1. D. Agard, Y. Hiraoka, P. Shaw, and J. Sedat, “Fluorescence microscopy in three dimensions,” Methods Cell Biol. **30**, 353–377 (1989). [CrossRef] [PubMed]

2. P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Sig. Proc. Mag. **23**, 32–45 (2006). [CrossRef]

3. S. Gibson and F. Lanni, “Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy, ” J. Opt. Soc. Am. A **8**, 1601–1613 (1991). [CrossRef]

*N*×

_{x}*N*×

_{y}*N*is the image size with

_{z}*N*being the axial dimension, then a rigorous representation amounts to using an

_{z}*N*number of 3D depth-varying PSFs (DV-PSFs). For a particular imaging set-up, if one has measured the PSF at zero depth, all DV-PSFs are predictable with a reasonable approximation by modifying the phase-retrieved pupil function [4

_{z}4. B. Hanser, M. Gustafsson, D. Agard, and J. Sedat, “Phase-retrieved pupil functions in wide-field fluorescent microscopy,” J. Microsc. **216**, 32–48 (2004). [CrossRef] [PubMed]

5. J. Shaevitz and D. Fletcher, “Enhanced three-dimensional deconvolution microscopy using a measured depth-varying point-spread function,” J. Opt. Soc. Am. A **24**, 2622–2627 (2007). [CrossRef]

8. E. Hom, F. Marchis, T. Lee, S. Haase, D. Agard, and J. Sedat, “AIDA: an adaptive image deconvolution algorithm with application to multi-frame and three-dimensional data,” J. Opt. Soc. Am. A **24**, 1580–1600 (2007). [CrossRef]

## 2. Developing the model

### 2.1. Tensor product PCA of 3D images

**X**

_{1},

**X**

_{2},…,

**X**

_{N}} be the given set of images. Let

**X**̄ be the mean image given by

**v**

_{i}be the 1D vector obtained by scanning the mean removed image

**X**

^{′}_{i}=

**X**

_{i}-

**X**̄. Let

**R**be the matrix defined by

**R**is an

*M*×

*N*matrix where

*M*is the total number of pixels in each of the images {

**X**

_{1},

**X**

_{2},… ,

**X**

_{N}}. Let the following be the singular value decomposition of

**R**.

**U**and

**V**are orthogonal matrices and

**D**is a diagonal matrix. Without the loss of generality, let us assume that the diagonal elements of

**D**are ordered with decreasing magnitude. Let

**u**

_{i}be the

*i*th column of

**U**, and

**P**

_{i}be the corresponding image obtained by reverse scanning

**u**

_{i}. Then in PCA, one choses the first

*B*column vectors of

**U**and expresses an input image

**X**

_{i}as

*c*

_{i,j}= 〈

**P**

_{j},

**X**

^{′}_{i}〉, with 〈,〉 representing the inner product, i.e., point-wise multiplication and summation. The user-defined parameter,

*B*, in the above expression represents the number of basis functions, which is chosen according to the required level of approximation. The diagonal elements of

**D**provide an estimate of approximation error. Specifically, the ratio,

*d*s represent the diagonal elements of

_{i}**D**. Typically, for an adequate level of approximation,

*B*is much smaller than

*N*.

**S**

_{1}(

*x*,

*y*,

*z*),…,

**S**

_{N}(

*x*,

*y*,

*z*)} be the set of 3D images. Let

*N*be the number of

_{z}*z*sections. For each value of

*z*, we form a set of 2D images given by 𝓢

^{(z)}= {

**S**

^{(z)}

_{1}(

*x*,

*y*),…,

**S**

^{(z)}

_{N}(

*x*,

*y*)}, where

**S**

^{(z)}

_{i}(

*x*,

*y*) =

**S**

_{i}(

*x*,

*y*,

*z*). We then apply PCA on each of the subsets to get basis functions

**P**

^{(z)}

_{1}(

*x*,

*y*),…,

**P**

^{(z)}

_{1}(

*x*,

*y*), where

*B*

_{1}is a user parameter. The images are expressed as

**S**̄

^{(z)}(

*x*,

*y*) is the mean image of the set 𝓢

^{(z)}and

^{′}to get basis functions

**P**

^{′}

_{1}(

*z*,

*j*),…,

**P**

^{′}

_{B2}(

*z*,

*j*), where

*B*

_{2}is the second user parameter. The images are approximated as

*c*

^{′}

_{i,k}= 〈

*C*(

_{i}*z*,

*j*),

**P**

^{′}

_{k}(

*z*,

*j*)〉. Substituting Eqs. (3) and (2) in Eq. (1) gives

*B*

_{1}is chosen such that the approximation error is negligible or zero in the first stage, and the second user parameter is chosen according to the desired final approximation error. By this way, the overall approximation error is governed only by the z-dependent second stage.

*N*+ 1 2D PCA decompositions. Even though the resulting basis function are not as optimal as the ones that might be obtained by direct 3D PCA (the impractical method), we will show experimentally that this approach works well for approximating 3D PSFs.

_{z}### 2.2. The PCA based imaging model

*z*or equivalently different positions of the focal plane [Fig. 1(a)]. Note that

*z*represents the position of the focal plane with respect to the interface plane. If there is no refractive index mismatch, the 3D image is expressed in terms of the object by a simple convolution with a 3D PSF. However, to express the image in terms of the object under the coverslip when there is a mismatch, we need the 4D function,

*h*

^{′}(

*x*,

*y*,

*z*,

*z*

^{′}), which represents the series of 2D image resulting from positioning the focal plane at

*z*and a point source at

*z*

^{′}[Fig. 1(b)].

*f*(

*x*,

*y*,

*z*) be the object (fluorescent dye structure) and

*g*(

*x*,

*y*,

*z*) be the image acquired by the microscope. The forward model is given by

*z*is the axial image coordinate and

*z*

^{′}is the axial object coordinate and the system response is a function of both.

*z*-

*z*

^{′}.In other words, the depth invariance corresponds to following relation:

*h*

_{0}(

*x*,

*y*,

*z*) is the standard depth independent PSF given by

*h*

_{0}(

*x*,

*y*,

*z*) =

*h*

^{′}(

*x*,

*y*,

*z*,0).

12. S. Wiersma, P. Torok, T. Visser, and P. Varga, “Comparison of different theories for focusing through a plane interface,” J. Opt. Soc. Am. B **14**, 1482–1490 (1997). [CrossRef]

*h*

^{′}. The term focus shift signifies the fact that, for a given value of object depth variable

*z*

^{′}, the maximum intensity is observed when the image depth variable

*z*is equal to

*az*

^{′}, where

*a*is a constant that depends on the refractive indices of the mounting and immersion mediums. To account for the focus shift, we consider the followed transformation [Fig. 1(c)]:

*z*is the plane containing the bead with focus shift. This representation is advantageous in the sense that, for a fixed value of

*x*,

*y*, and

*z*, the variation with respect to axial object coordinate

*z*

^{′}is lower compared to the former representation.

*a*in the above expression disappears if it is translated into discrete summations assuming that discretization step size of

*z*is

*a*times larger than that of

*z*

^{′}. Assuming now that the variables represent discrete indices, the above model becomes

*P*

_{1}(

*x*,

*y*,

*z*),…,

*P*(

_{B}*x*,

*y*,

*z*) be the selected principal components [

*Q*’s in Eq. (5)]. Then, the function

_{k}*h*(·,·,·,·) can be approximated as

*h*̄(

*x*,

*y*,

*z*) is the function obtained by averaging

*h*(

*x*,

*y*,

**z**,

*z*

^{′}) along

*z*

^{′}, and

*h*replaced by

*h*

^{(B)}

_{a}, where the latter is given by

*N*×

_{x}*N*×

_{y}*N*be the image size. A quick examination of Eq. (10) might indicate that its implementation will require

_{z}*B*+1 number of 3D FFTs and a 3D IFFT. However, observing the fact that input to

*B*+ 1 convolutions differ only axially, it can be verified that

*B*+1 number of 3D FFTs can be replaced by

*N*times 2D FFTs and then (

_{z}*B*+1)

*N*times 1D FFTs, providing additional computational efficiency.

_{x}N_{y}## 3. Results

7. C. Preza and J. Conchello, “Depth-variant maximum likelihood restoration for three-dimensional fluorescence microscopy,” J. Opt. Soc. Am. A **21**, 1593–1601 (2004). [CrossRef]

*z*

^{′}-axis is divided into some number of intervals (strata) of thickness

*T*and the DV-PSFs at the edge of the intervals are used to approximate all the DV-PSFs using bilinear interpolation.

5. J. Shaevitz and D. Fletcher, “Enhanced three-dimensional deconvolution microscopy using a measured depth-varying point-spread function,” J. Opt. Soc. Am. A **24**, 2622–2627 (2007). [CrossRef]

*z*

^{′}) using an optical trap (refer to [5

5. J. Shaevitz and D. Fletcher, “Enhanced three-dimensional deconvolution microscopy using a measured depth-varying point-spread function,” J. Opt. Soc. Am. A **24**, 2622–2627 (2007). [CrossRef]

*x*and

*y*was 42 nm and that of

*z*

^{′}was 75 nm.

*z*step size was chosen appropriately to compensate for the focus shift. The objective numerical aperture was 1.49 NA. We also used the same specifications to generate theoretical DV-PSFs using the method described in [4

4. B. Hanser, M. Gustafsson, D. Agard, and J. Sedat, “Phase-retrieved pupil functions in wide-field fluorescent microscopy,” J. Microsc. **216**, 32–48 (2004). [CrossRef] [PubMed]

*xz*sections for various values

*z*

^{′}with

*y*chosen to be at the center of the PSFs, and

*xz*

^{′}sections for various values

*z*with

*y*chosen similarly. In all the visualizations, the top left corner corresponds to the origin of the axes. To make low intensities visible, we display the fourth root of the intensity whenever the images are positive. The exceptions are the approximation error images, and the PCA basis functions.

*B*, we define the normalized approximation error as

### 3.1. Application to theoretical PSFs

*z*

^{′}=

*i*× 75nm,

*i*= 0,…, 40}. Figure 2 displays the

*xz*and

*xz*

^{′}sections and Figure 3 shows few first PCA basis functions.

7. C. Preza and J. Conchello, “Depth-variant maximum likelihood restoration for three-dimensional fluorescence microscopy,” J. Opt. Soc. Am. A **21**, 1593–1601 (2004). [CrossRef]

*h*̄ and the first principal component; whereas, for the strata method, this means using the PSF corresponding to values of

*z*

^{′}equal to zero and bottom of the total range. Figure 5 compares the approximated PSFs with the originals for both methods. In the strata method, for the PSFs at the top and the bottom, the approximated PSFs are identical to the originals since these PSFs are themselves the basis functions. However, for the total range of

*z*

^{′}, the mean approximation error is much higher than that of PCA method (about 17 times). This is evident from the approximated images in the middle (for

*z*

^{′}= 750 nm, 1500 nm, 2250 nm). Figure 6 shows the

*xz*

^{′}sections of the same results. While visually less obvious, comparable behavior extends to any number of basis functions.

### 3.2. Application to measured PSFs

*xz*and

*xz*

^{′}sections and Fig. 8 shows first two PCA basis functions. We observed that, subsequent basis functions contain background noise as well as high frequency structures that lie well beyond the theoretical frequency support of the microscope, which possibly originate from the Poisson noise. This indicates that two bases function are sufficient to represent the essential details of DV-PSFs for the range of depth considered in this experiment.

*z*

^{′}. However, for the other values of

*z*

^{′}, the PCA method yields a smaller error, and all of the approximated PSFs are nearly noise-free. This result is also evident from the displayed error images [Fig. 10(b)]. In contrast to the case of theoretical PSFs, here the PCA method provides only a two fold improvement in the error, due to the inclusion noise. Figure 11 displays the

*xz*

^{′}section of the same result for

*z*= 0. We can see significant line artifacts in the strata approximated PSFs which originates from interpolated noise.

8. E. Hom, F. Marchis, T. Lee, S. Haase, D. Agard, and J. Sedat, “AIDA: an adaptive image deconvolution algorithm with application to multi-frame and three-dimensional data,” J. Opt. Soc. Am. A **24**, 1580–1600 (2007). [CrossRef]

## 4. Conclusions

## Acknowledgements

## References and links

1. | D. Agard, Y. Hiraoka, P. Shaw, and J. Sedat, “Fluorescence microscopy in three dimensions,” Methods Cell Biol. |

2. | P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Sig. Proc. Mag. |

3. | S. Gibson and F. Lanni, “Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy, ” J. Opt. Soc. Am. A |

4. | B. Hanser, M. Gustafsson, D. Agard, and J. Sedat, “Phase-retrieved pupil functions in wide-field fluorescent microscopy,” J. Microsc. |

5. | J. Shaevitz and D. Fletcher, “Enhanced three-dimensional deconvolution microscopy using a measured depth-varying point-spread function,” J. Opt. Soc. Am. A |

6. | C. Preza and J. Conchello, “Image estimation account for point-spread function depth variation in three-dimensional fluorescence microscopy,” Proc. SPIE (2003), pp. 1–8. |

7. | C. Preza and J. Conchello, “Depth-variant maximum likelihood restoration for three-dimensional fluorescence microscopy,” J. Opt. Soc. Am. A |

8. | E. Hom, F. Marchis, T. Lee, S. Haase, D. Agard, and J. Sedat, “AIDA: an adaptive image deconvolution algorithm with application to multi-frame and three-dimensional data,” J. Opt. Soc. Am. A |

9. | C. Vonesch and M. Unser, “A Fast Thresholded Landweber Algorithm for Wavelet-Regularized Multidimensional Deconvolution,” IEEE Tran. Imag. Proc. |

10. | C. Vonesch and M. Unser, “A Fast Multilevel Algorithm for Wavelet-Regularized Image Restoration,” IEEE Tran. Imag. Proc. |

11. | I. Jolliffe, |

12. | S. Wiersma, P. Torok, T. Visser, and P. Varga, “Comparison of different theories for focusing through a plane interface,” J. Opt. Soc. Am. B |

**OCIS Codes**

(110.0110) Imaging systems : Imaging systems

(110.0180) Imaging systems : Microscopy

(110.6880) Imaging systems : Three-dimensional image acquisition

**ToC Category:**

Imaging Systems

**History**

Original Manuscript: January 8, 2010

Revised Manuscript: March 5, 2010

Manuscript Accepted: March 5, 2010

Published: March 15, 2010

**Virtual Issues**

Vol. 5, Iss. 7 *Virtual Journal for Biomedical Optics*

**Citation**

Muthuvel Arigovindan, Joshua Shaevitz, John McGowan, John W. Sedat, and David A. Agard, "A Parallel Product-Convolution approach for representing the depth varying Point Spread Functions in 3D widefield microscopy based on principal
component analysis," Opt. Express **18**, 6461-6476 (2010)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-7-6461

Sort: Year | Journal | Reset

### References

- D. Agard, Y. Hiraoka, P. Shaw, and J. Sedat, "Fluorescence microscopy in three dimensions," Methods Cell Biol. 30, 353-377 (1989). [CrossRef] [PubMed]
- P. Sarder and A. Nehorai,"Deconvolution methods for 3-D fluorescence microscopy images," IEEE Sig. Proc. Mag. 23, 32-45 (2006). [CrossRef]
- S. Gibson and F. Lanni, "Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy, " J. Opt. Soc. Am. A 8, 1601-1613 (1991). [CrossRef]
- B. Hanser, M. Gustafsson, D. Agard, and J. Sedat, "Phase-retrieved pupil functions in wide-field fluorescent microscopy," J. Microsc. 216, 32-48 (2004). [CrossRef] [PubMed]
- J. Shaevitz and D. Fletcher, "Enhanced three-dimensional deconvolution microscopy using a measured depthvarying point-spread function," J. Opt. Soc. Am. A 24, 2622-2627 (2007). [CrossRef]
- C. Preza and J. Conchello, "Image estimation account for point-spread function depth variation in threedimensional fluorescence microscopy," Proc. SPIE 1-8 (2003).
- C. Preza and J. Conchello, "Depth-variant maximum likelihood restoration for three-dimensional fluorescence microscopy," J. Opt. Soc. Am. A 21, 1593-1601 (2004). [CrossRef]
- E. Hom, F. Marchis, T. Lee, S. Haase, D. Agard, and J. Sedat, "AIDA: an adaptive image deconvolution algorithm with application to multi-frame and three-dimensional data," J. Opt. Soc. Am. A 24, 1580-1600 (2007). [CrossRef]
- C. Vonesch and M. Unser, "A Fast Thresholded Landweber Algorithm forWavelet-Regularized Multidimensional Deconvolution," IEEE Tran. Imag. Proc. 17, 539-549 (2008). [CrossRef]
- C. Vonesch and M. Unser, "A Fast Multilevel Algorithm for Wavelet-Regularized Image Restoration," IEEE Tran. Imag. Proc. 18, 509-523 (2009). [CrossRef]
- I. Jolliffe, Principal component analysis (Springer, 2002).
- S. Wiersma, P. Torok, T. Visser, and P. Varga, "Comparison of different theories for focusing through a plane interface," J. Opt. Soc. Am. B 14, 1482-1490 (1997). [CrossRef]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.