Indrasen Bhattacharya: Difference between revisions

Latest revision as of 00:37, 8 December 2020

Introduction

Aberrations in the wavefront profile degrade the resolving power of imaging systems. It is highly desirable to measure and compensate any aberrations present in an imaging system. Such imaging systems may include microscopes, telescopes, lithography steppers, iPhone cameras, AR/VR goggles as well as the human eye. We shall consider the more traditional case of optical microscopes and lithography steppers in this work, but we do not rule out extensions to other imaging systems. The point spread function (PSF) at the image plane determines the resolving power of microscopes and the feature size in lithography systems. It is important to have a narrow PSF that is preferably diffraction limited by the numerical aperture (NA) of the imaging system. Wavefront aberrations reduce the sharpness of the PSF and lead to undesirable asymmetries and artifacts. The wavefront aberrations are typically expressed in an orthogonal radial basis on the unit circular pupil called the Zernike polynomials. The Zernike polynomials $R_{n}^{m} (ρ)$ , together with azimuthal sinusoids, form a complete basis on the unit circle: any 2-dimensional phase function can be expressed as an appropriate linear combination of these polynomials. The weighting factors $α_{n}^{m}$ of the Zernike polynomials are termed as Zernike coefficients. Since any wavefront aberration profile can be fitted to a series of Zernike coefficients: the pupil function and Zernike coefficients are equivalent representations. Conveniently, the Zernike polynomials can be interpreted as specific physical aberrations such as spherical aberration, coma, astigmatism, trefoil and others, which makes the Zernike coefficients an intuitive framework to work with.

More explicitly, if the phase due to aberrations is expressed in normalized Cartesian basis as: $e^{i Φ (κ_{x}, κ_{y})} = e^{i Φ (ρ, θ)}$ , the Zernike polynomial expansion of the aberration phase $Φ$ is given by:

$Φ (ρ, θ) = Σ_{n, m} α_{n}^{m} R_{n}^{m} (ρ) c o s (m θ)$

where $0 \leq ρ \leq 1, 0 \leq θ \leq 2 π$ . For the purpose of this analysis, we consider only the cosine terms without loss of generality.

It is conceptually straightforward to calculate the point spread function from the wavefront profile: they are related by a Fourier transform where the aberrations occur in the phase of the pupil function. However, in this problem, we are interested in the reciprocal calculation of determining what is wrong with the system based on how the PSF appears. This is the inverse problem: determining the Zernike coefficients from the aberrated point spread function. It is important to explicitly state the inputs and outputs of the computational procedure we are hoping to undertake. If the wavefront aberration is given as a function of the NA-normalized pupil coordinates (in Cartesian basis) as $Φ (κ_{x}, κ_{y})$ , the point spread function can be expressed as:

$U (x, y; f) = \frac{1}{π} \int_{κ_{x}^{2} + κ_{y}^{2} \leq 1} \int e^{i f (κ_{x}^{2} + κ_{y}^{2})} e^{i Φ (κ_{x}, κ_{y})} e^{[2 π i κ_{x} x + 2 π i κ_{y} y]} d κ_{x} d κ_{y}$

Note that this formula neglects the Radiometric correction factor, which becomes relevant for high numerical aperture systems. It also approximates the defocus phase for the small NA case. However, these approximations typically suffice for NA smaller than 0.6, which captures many relevant systems.

We do not expect to be able to retrieve a complex phase from the intensity point spread function at a single imaging plane, unless we make additional prior assumptions regarding the phase. Instead, we use the approach proposed by Van Der Avoort in [2], which relies on using multiple point spread functions. In the approach proposed by these authors, we hope to acquire a focal progression of point spread functions: $U (x, y, - f_{m a x}) ..., U (x, y, - δ f), U (x, y, 0), U (x, y, δ f), ... U (x, y, f_{m a x})$ and perform computations that will let us closely recover the Zernike coefficients $α_{n}^{m}$ for the phase function. The authors have shown that under a linearizing assumption, the aberration phase has a one to one correspondence with a focal progression of intensity point spread functions: it can be exactly inverted to obtain the Zernike coefficients $α_{n}^{m}$ . This approach is feasible for lithography and optical microscopy imaging systems, and indeed has been used in these applications. It will require some thought and experiment design to extend these ideas to the human eye - we leave this for future work.

As a concrete example, a focal progression is displayed in the figure below. The focus increases in steps of 1 normalized unit (see 'Notations and Conventions') starting at -6 on the top left corner, and increasing rightwards. This was generated using a numerical implementation of the fourier transform formula above, using the DFT method. The units for the X and Y axes are normalized according to the conventions specified in the section on 'Notations and Conventions'. The system was simulated to be aberrated by vertical astigmatism and horizontal coma. It can be qualitatively observed that the impacts of the aberration become much more apparent at a higher defocus, particularly the coma.

Fig. 1: Focal progression of intensity for an optical system with two higher order aberrations: vertical astigmatism $(α_{2}^{2} = 0.5)$ and horizontal coma $(α_{3}^{1} = 0.3)$

Notation and Conventions

In this project, we have focused on one particular method: the one proposed in [2]. We have used python for a computational implementation of the technique, following the numerical conventions used in the paper. In particular, all spatial units $(X, Y, Z)$ have been normalized to the numerical aperture of the system as follows:

$x = X \times \frac{2 π N A}{λ}, y = Y \times \frac{2 π N A}{λ}, f = Z \times \frac{π N A^{2}}{λ}$

These are the spatial quantities used in the Fourier transform expression above.

We are also using the convention for Zernike polynomials used in reference [1], appendix VII and section 9.2. Please note that this convention is different from the normalized convention used on the Wikipedia page. As is the case for most conventions, we use m to denote the azimuthal degree of the Zernike polynomial and n the radial degree. n takes values m, m+2, m+4, ... Essentially, a certain azimuthal order can only occur for a sufficiently high radial order. Also, the difference between the radial and azimuthal order of the polynomials is always even. These are certain mathematical aspects of the Zernike polynomial expansion of the pupil phase that are useful to have at reference.

Certain choices regarding discretization and infinite series truncation have also been made in this implementation. Where these choices are different from the implementation in the original paper, I have clarified in the text in the relevant sections below.

Analytical Forward Model

An important aspect of the inverse problem solution is the development of a linearized model for the field point spread function. As a result, it will turn out that the field PSF is linear in the Zernike coefficients and the intensity PSF is quadratic. The error in the field PSF is of third order in the Zernike coefficients. This linearization will later turn out to be greatly beneficial for solving the inverse problem. We omit the derivation here but specify the key assumptions and results.

If the aberration phase is small, we can linearize it as follows: $e^{[i Φ (ρ, θ)]} \sim 1 + i Φ (ρ, θ) = 1 + i Σ_{n, m} α_{n}^{m} R_{n}^{m} (ρ) c o s (m θ)$

Substituting this simplification into the Fourier transform formula and using certain mathematical identities, we obtain a closed form solution for the field point spread function in normalized radial coordinates:

$U (r, ϕ; f) = 2 V_{0}^{0} (r, f) + 2 i Σ_{n, m} i^{m} α_{n}^{m} V_{n}^{m} (r, f) c o s (m ϕ)$

We note that this formula is linear in the Zernike coefficients $α_{n}^{m}$ . The special functions $V_{n}^{m}$ used here are given by:

$V_{n}^{m} (r, f) = \int_{0}^{1} ρ e^{(i f ρ^{2})} R_{n}^{m} (ρ) J_{m} (2 π r ρ) d ρ$ , where $J_{m} (x)$ is the Bessel function corresponding to the azimuthal order.

The authors further express the special functions as an infinite summation using Bessel functions. We omit that expression here, but it can be found in [2], equation (10). This closed form analytical expression has been implemented in python. Please note that the infinite summation in equation (10), using the indexing variable l requires termination for a computational implementation. We will discuss this in the next section.

Note that the intensity PSF is simply given as: $I = | U |^{2}$ and can be written based on the above special functions as well. The intensity PSF will turn out to be quadratic in the Zernike coefficients, under this linearization assumption.

Verification of ENZ Forward Model

The authors of [2] have termed the analytical closed form solution for the field as the extended Nieboer-Zernike approach (ENZ approach). Since the inverse problem solution relies on these special functions and the accuracy of the forward model, we decided to independently test the forward model first. This will also provide some insight on how to terminate the infinite summation used in equation (10) of [2] for expressing the special functions $V_{n}^{m} (r, f)$ .

In order to obtain a ground truth comparison, we have implemented the Fourier transform expression for the PSF using in-built DFT packages in python:

$U (x, y; f) = \frac{1}{π} \int_{κ_{x}^{2} + κ_{y}^{2} \leq 1} \int e^{i f (κ_{x}^{2} + κ_{y}^{2})} e^{i Φ (κ_{x}, κ_{y})} e^{[2 π i κ_{x} x + 2 π i κ_{y} y]} d κ_{x} d κ_{y}$

The phase and resulting PSF are shown in Fig. 2 below. This was for a case with both astigmatism, coma and defocus. The amount of defocus is 3 normalized units; astigmatism and coma are same as in Fig. 1.

Fig. 2: Ground truth PSF calculation using the DFT approach. Top panel: pupil phase and 1D cuts, Bottom panel: resulting intensity PSF in real space and 1D cuts

In order to verify the analytical formula based on the special functions, we have implemented this and compared it with the ground truth in Fig. 2. For the same defocus and aberration parameters, the resulting point spread function image is shown below in Fig. 3. The difference with the ground truth, normalized to the maximum of the ground truth PSF is shown in Fig. 4. As observed, the intensity PSF resulting from the analytical approach is a good approximation for the DFT based calculation. The error is on the order of 5% close to the peak of the PSF. This is a good sanity check to show that the analytical approximation is accurate enough for us to continue with this approach.

Fig. 3: PSF resulting from the closed form analytical solution, using the same aberration and defocus parameters as Fig. 2.

Fig. 4: Difference in PSF, normalized by maximum of the PSF. As can be seen, errors on the order of 5% can be observed.

We note that this small error was able to be achieved by terminating the infinite summation in equation (10) of [2] with 10 terms. Therefore, we continue with this choice for the inverse problem solution. Please note that the aberrations specified in Zernike coefficients are in radians (vertical astigmatism $(α_{2}^{2} = 0.5)$ and horizontal coma $(α_{3}^{1} = 0.3)$ ). These can be converted to wavelength units, resulting in astigmatism of 79.6 milliwaves, and horizontal coma of 47.7 milliwaves. These are fairly realistic numbers and this case (among other sanity checks) will be used for the rest of this study.

Inverse Problem Solution

Here we discuss how to retrieve the pure phase aberrations from a defocus series of PSFs $I (r, ϕ; f)$ . In a real world situation, these PSFs would be obtained empirically. However, for the sake of our computational experiments, I have used the Fourier transform method to generate a ground truth PSF progression. Next, we summarize the key steps in retrieving the Zernike coefficients. For a flavor of the derivation, please see section 3 of reference [2]. The key details involve:

1. Expressing the intensity PSF in the analytical form $I = | U |^{2}$

2. Inner product of the measured PSFs $I_{m e a s} (r, ϕ; f)$ with the azimuthal sinusoid in order to obtain azimuth-free special functions $Ψ_{m e a s}^{m} (r, f)$

3. Calculation of intensity basis functions $Ψ_{n}^{m} (r, f)$ from the original special functions $V_{n}^{m} (r, f)$

4. For each azimuthal order $m$ , the formulation of a system of equations for obtaining the Zernike coefficients for all valid radial orders corresponding to that azimuthal order: $α_{n}^{m}, n = m, m + 2, m + 4, .... m + 2 M$ ; where $M$ is used to express an appropriately chosen termination of the Zernike polynomial radial order

5. The inversion of the system of equations (by inverting the corresponding Gram matrix) in order to obtain the Zernike coefficients

The procedure appears to be rather involved, but it can be summarized as follows. The measured intensity PSF is written in terms of unknown Zernike coefficients. Then, the formula is projected into the basis of the azimuthal sinusoids and the intensity basis functions. Once this is done, the quadratic term is rigorously found to be zeroed out, leaving only terms linear in the Zernike coefficients. The linear expression can be solved for the Zernike coefficients just like any other system of linear equations.

Here we state the formula for the measurement projection (equation 23 in [2]):

$Ψ_{m e a s}^{m} (r, f) = \frac{1}{2 π} \int_{0}^{2 π} I_{m e a s} (r, ϕ; f) c o s (m ϕ) d ϕ$

The intensity basis functions can be pre-calculated and are given by (equation 24 in [2]):

$Ψ_{n}^{m} (r, f) = - 8 ϵ_{m}^{- 1} I m [i^{m} V_{n}^{m} (r, f) V_{0}^{0 *} (r, f)]$

And each linear equation for solving for the Zernike coefficients is given simply by:

$Σ_{n = m, m + 2, ... m + 2 M} α_{n}^{m} (Ψ_{n}^{m}, Ψ_{n^{'}}^{m}) = (Ψ_{m e a s}^{m}, Ψ_{n^{'}}^{m})$

The full system of equations can be expressed by vectorizing over the radial order $n^{'}$ , leading to a matrix formulation. There is one matrix for each azimuthal order, and the size of the matrix is [M, M], given by how we originally chose to truncate the Zernike polynomial radial order. By inverting this matrix equation we obtain a vector of size [M,1] containing all the relevant Zernike coefficients for that azimuthal order. This needs to be done for each azimuthal order that we choose to model. I have explained the discretization and truncation choices in the next section.

Results

We have implemented the ENZ inverse calculation on a small scale problem and compared with ground truth. The following three cases were tested:

Case 1: Sanity check with no aberrations, all Zernike coefficients expected to be zero

Case 2: Only vertical astigmatism ( $α_{2}^{2}$ ) is present

Case 3: Vertical astigmatism and horizontal coma ( $α_{3}^{1}$ ) are both present

The results below (with retrieved aberration in shown as red) demonstrate that there is a reasonable qualitative correspondence in the retrieved aberrations and ground truth. However, the quantitative correspondence was found to be lacking. Some of the reasons for this are explained at the end of this section. Please note that the aberrations are specified by both the descriptive name and the (m, n) order of the Zernike coefficient in the comparison graphs below.

Fig. 5: Results for Case 1 (sanity check with no aberrations)

Fig. 6: Results for Case 2 (only one form of aberration)

Fig. 7: Results for Case 3 (two dominant aberrations: same as the case shown in the PSF progression)

Due to time and resource constraints, the simulation settings were rather coarse. For instance, the choice of the quantity M, the number of radial orders per azimuthal order could only be 2. Ideally, this number would be on the order of 10 in order to consider a sufficiently large number of radial orders. The Gram matrix method attempts a best fit of the data to the basis functions included in the projection, up to M. Hence, the results may be poor if M is rather low. Another constraint was the spatial discretization, which could not be too fine due to memory and run time constraints on my personal computer. It would be very interesting to perform a more detailed study exploring the impact of focus and spatial discretization, as well as M on the accuracy of the results.

Limitations of Method and Future Work

There are some well known limitations of the approach as presented in this write-up (which have been addressed in later publications [3]-[5]):

1. The approach is valid only for pure phase aberrations. There is an extension to phase and amplitude aberrations using complex Zernike coefficients that may be worth investigating.

2. The radiometric correction factor has been neglected in the current formulation.

3. A vectorial formulation is required for higher numerical aperture systems.

4. The effect of finite source extent has been neglected, although it is fairly straightforward to incorporate.

Ideally, a little more investigation into the accuracy, numerical stability and a few more test cases would have been great - but hopefully this is a good starting point.

References

[1] M. Born and E. Wolf, Principles of Optics, (Cambridge University Press, Cambridge, 2002)

[2] C. Van Der Avoort, J. J. M. Braat, P. Dirksen, A. J. E. M. Janssen, "Aberration retrieval from the intensity point-spread function in the focal region using the extended Nijboer-Zernike approach", J. Mod. Opt., 52, 12, 2005, pp. 1695-1728

[3] A. J. E. M. Janssen, "Extended Nijboer-Zernike approach for the computation of optical point-spread functions", J. Opt. Soc. Am. A, 19, 5, 2002

[4] J. J. M. Braat, P. Dirksen, A. J. E. M. Janssen, A. S. van des Nes, "Extended Nijboer-Zernike representation of the vector field in the focal region of an aberrated high-aperture optical system", J. Opt. Soc. Am. A, 20, 12, 2003

[5] A. J. E. M. Janssen, J. J. M. Braat, P. Dirksen, "On the computation of the Nijboer-Zernike aberration integrals at arbitrary defocus", J. Mod. Opt., 51, 5, 2004, pp. 687-703

[6] P. Dirksen, J. Braat, A. J. E. M. Janssen, "Estimating resist parameters in optical lithography using the extended Nijboer-Zernike theory", J. Microlith., Microfab., Microsyst. 5(1), 013005, 2006

[7] Extended Nijboer-Zernike (ENZ) Analysis and Aberration Retrieval research site (Netherlands group): https://www.nijboerzernike.nl/

Appendix

The computational implementation (jupyter notebook) and results for this project can be found in this zipped folder: File:Indrasen final.zip

@@ Line 1: / Line 1: @@
 == Introduction ==
-Contrast sensitivity function (CSF) is a subjective measurement of an ability of the visual system to detect a low contrast pattern stimuli. The stimuli  considered usually are vertical sinusoidal Gabor patches of  decreasing shades of black to grey. This use of sine wave gratings was first introduced in vision by Schade [1] and was subsequently used by early investigators to measure basic visual sensitivity [2]. The resulting measurement is used to validate the representation of the eyes' visual performance as it complements the visual acuity. CSF  is synonyms with an audiogram where a person’s highest detectable pitch is measured and well as the ability to hear all lower pitches [3].
+Aberrations in the wavefront profile degrade the resolving power of imaging systems. It is highly desirable to measure and compensate any aberrations present in an imaging system. Such imaging systems may include microscopes, telescopes, lithography steppers, iPhone cameras, AR/VR goggles as well as the human eye. We shall consider the more traditional case of optical microscopes and lithography steppers in this work, but we do not rule out extensions to other imaging systems. The point spread function (PSF) at the image plane determines the resolving power of microscopes and the feature size in lithography systems. It is important to have a narrow PSF that is preferably diffraction limited by the numerical aperture (NA) of the imaging system. Wavefront aberrations reduce the sharpness of the PSF and lead to undesirable asymmetries and artifacts. The wavefront aberrations are typically expressed in an orthogonal radial basis on the unit circular pupil called the Zernike polynomials. The Zernike polynomials <math>R^m_n(\rho)</math>, together with azimuthal sinusoids, form a complete basis on the unit circle: any 2-dimensional phase function can be expressed as an appropriate linear combination of these polynomials. The weighting factors <math>\alpha^m_n</math> of the Zernike polynomials are termed as Zernike coefficients. Since any wavefront aberration profile can be fitted to a series of Zernike coefficients: the pupil function and Zernike coefficients are equivalent representations. Conveniently, the Zernike polynomials can be interpreted as specific physical aberrations such as spherical aberration, coma, astigmatism, trefoil and others, which makes the Zernike coefficients an intuitive framework to work with.
-The CSF measurements are usually acquired with  small patches of sinusoidal grating designed to fall within   few central degrees of the visual field. It is well known that the CSF decreases as one measures contrast sensitivity at increasingly peripheral locations in the visual field. The reasons for such decreased CSF is attributed to a number of neural factors [4]. The human eye is structured such that the distribution of the cone mosaic falls off rapidly as a function of visual eccentricity, so that there are fewer sensors available to detect and encode the incoming stimuli. Towards the periphery the amount of retinal ganglion cells’ density falls as well. This structure of the cones is also a key factor in deciding the CSF. In particular, for this  project, we make an effort to explore the role of  inference engines in shaping the CSF.
+More explicitly, if the phase due to aberrations is expressed in normalized Cartesian basis as: <math> e^{i\Phi(\kappa_x, \kappa_y)} = e^{i\Phi(\rho, \theta)} </math>, the Zernike polynomial expansion of the aberration phase <math>\Phi</math> is given by:
-== Background ==
+<math>\Phi(\rho, \theta) = \Sigma_{n, m} \alpha^m_n R^m_n(\rho) cos(m\theta)</math>
+where <math>0 \le \rho \le 1, 0 \le \theta \le 2\pi </math>. For the purpose of this analysis, we consider only the cosine terms without loss of generality.
+It is conceptually straightforward to calculate the point spread function from the wavefront profile: they are related by a Fourier transform where the aberrations occur in the phase of the pupil function. However, in this problem, we are interested in the reciprocal calculation of determining what is wrong with the system based on how the PSF appears. This is the inverse problem: determining the Zernike coefficients from the aberrated point spread function. It is important to explicitly state the inputs and outputs of the computational procedure we are hoping to undertake. If the wavefront aberration is given as a function of the NA-normalized pupil coordinates (in Cartesian basis) as <math>\Phi(\kappa_x, \kappa_y)</math>, the point spread function can be expressed as:
+<math> U(x,y; f) = \frac{1}{\pi}\int_{\kappa_x^2 + \kappa_y^2 \le 1} \int e^{if(\kappa_x^2 + \kappa_y^2)} e^{i\Phi(\kappa_x, \kappa_y)} e^{[2\pi i \kappa_x x + 2\pi i \kappa_y y] } d\kappa_x d\kappa_y </math>
+Note that this formula neglects the Radiometric correction factor, which becomes relevant for high numerical aperture systems. It also approximates the defocus phase for the small NA case. However, these approximations typically suffice for NA smaller than 0.6, which captures many relevant systems.
+We do not expect to be able to retrieve a complex phase from the intensity point spread function at a single imaging plane, unless we make additional prior assumptions regarding the phase. Instead, we use the approach proposed by Van Der Avoort in [2], which relies on using multiple point spread functions. In the approach proposed by these authors, we hope to acquire a focal progression of point spread functions: <math> U(x, y, -f_{max}) ... , U(x, y, -\delta f), U(x, y, 0), U(x, y, \delta f), ... U(x, y, f_{max}) </math> and perform computations that will let us closely recover the Zernike coefficients <math>\alpha^m_n</math> for the phase function. The authors have shown that under a linearizing assumption, the aberration phase has a one to one correspondence with a focal progression of intensity point spread functions: it can be exactly inverted to obtain the Zernike coefficients <math>\alpha^m_n</math>. This approach is feasible for lithography and optical microscopy imaging systems, and indeed has been used in these applications. It will require some thought and experiment design to extend these ideas to the human eye - we leave this for future work.
+As a concrete example, a focal progression is displayed in the figure below. The focus increases in steps of 1 normalized unit (see 'Notations and Conventions') starting at -6 on the top left corner, and increasing rightwards. This was generated using a numerical implementation of the fourier transform formula above, using the DFT method. The units for the X and Y axes are normalized according to the conventions specified in the section on 'Notations and Conventions'. The system was simulated to be aberrated by vertical astigmatism and horizontal coma. It can be qualitatively observed that the impacts of the aberration become much more apparent at a higher defocus, particularly the coma.
+[[File:Fig1_indrasen.png]]
+Fig. 1: Focal progression of intensity for an optical system with two higher order aberrations: vertical astigmatism <math>(\alpha_2^2 = 0.5)</math> and horizontal coma <math>(\alpha_3^1 = 0.3)</math>
+== Notation and Conventions ==
+In this project, we have focused on one particular method: the one proposed in [2]. We have used python for a computational implementation of the technique, following the numerical conventions used in the paper. In particular, all spatial units <math>(X, Y, Z) </math> have been normalized to the numerical aperture of the system as follows:
+<math>x = X \times \frac{2\pi NA}{\lambda}, y = Y \times \frac{2\pi NA}{\lambda}, f = Z \times \frac{\pi NA^2}{\lambda}</math>
+These are the spatial quantities used in the Fourier transform expression above.
+We are also using the convention for Zernike polynomials used in reference [1], appendix VII and section 9.2. Please note that this convention is different from the normalized convention used on the Wikipedia page. As is the case for most conventions, we use m to denote the azimuthal degree of the Zernike polynomial and n the radial degree. n takes values m, m+2, m+4, ... Essentially, a certain azimuthal order can only occur for a sufficiently high radial order. Also, the difference between the radial and azimuthal order of the polynomials is always even. These are certain mathematical aspects of the Zernike polynomial expansion of the pupil phase that are useful to have at reference.
+Certain choices regarding discretization and infinite series truncation have also been made in this implementation. Where these choices are different from the implementation in the original paper, I have clarified in the text in the relevant sections below.
+== Analytical Forward Model ==
+An important aspect of the inverse problem solution is the development of a linearized model for the field point spread function. As a result, it will turn out that the field PSF is linear in the Zernike coefficients and the intensity PSF is quadratic. The error in the field PSF is of third order in the Zernike coefficients. This linearization will later turn out to be greatly beneficial for solving the inverse problem. We omit the derivation here but specify the key assumptions and results.
+If the aberration phase is small, we can linearize it as follows:
+<math>e^{[i\Phi(\rho, \theta)]} \sim 1 + i\Phi(\rho, \theta) = 1 + i\Sigma_{n,m} \alpha^m_n R^m_n(\rho) cos(m\theta)</math>
+Substituting this simplification into the Fourier transform formula and using certain mathematical identities, we obtain a closed form solution for the field point spread function in normalized radial coordinates:
+<math>U(r,\phi; f) = 2V^0_0(r,f) + 2i\Sigma_{n,m}i^m \alpha^m_n V^m_n(r,f) cos(m\phi)</math>
+We note that this formula is linear in the Zernike coefficients <math>\alpha^m_n</math>. The special functions <math>V^m_n</math> used here are given by:
+<math>V^m_n(r,f) = \int_0^1 \rho e^{(if\rho^2)} R^m_n(\rho) J_m(2\pi r \rho) d\rho</math>, where <math>J_m(x)</math> is the Bessel function corresponding to the azimuthal order.
-The human visual system’s CSF or the modulation transfer function (MTF) fundamentally characterises eye’s spatial frequency response and thus, one may think of it as a bandpass filter. This bandpass nature, however, is determined by a variety of factors such as the front-end optics, cone distribution geometry and the neural mechanism that is responsible for this interpretation which is typically inferred using the machine learning techniques such as the kNN and SVM [5].
+The authors further express the special functions as an infinite summation using Bessel functions. We omit that expression here, but it can be found in [2], equation (10). This closed form analytical expression has been implemented in python. Please note that the infinite summation in equation (10), using the indexing variable l requires termination for a computational implementation. We will discuss this in the next section.
-The methods of CSF measurement usually employ sine-wave gratings of a fixed frequency. By varying the frequency, a set of stimuli are constructed. The response of the eye to these stimuli are then determined. This procedure is repeated for a large number of grating frequencies and the corresponding responses are recorded.  A mechanism that emulates the brain inference is applied to these responses to find the resulting curve called the contrast sensitivity function and is illustrated in Fig. 1. [6].
+Note that the intensity PSF is simply given as: <math>I = |U|^2</math> and can be written based on the above special functions as well. The intensity PSF will turn out to be quadratic in the Zernike coefficients, under this linearization assumption.
-[[File:Fig1_okkeun.png]]
+== Verification of ENZ Forward Model ==
-Fig 1. The plot of CSF varies with retinal eccentricity.
+The authors of [2] have termed the analytical closed form solution for the field as the extended Nieboer-Zernike approach (ENZ approach). Since the inverse problem solution relies on these special functions and the accuracy of the forward model, we decided to independently test the forward model first. This will also provide some insight on how to terminate the infinite summation used in equation (10) of [2] for expressing the special functions <math>V^m_n(r,f)</math>.
-From the measurements, contrast sensitivity score is determined such that it is equal to (1/threshold) for the given spatial form.  The contrast sensitivity scores obtained for each of the  sine-wave gratings examined are then plotted as a function of target spatial frequency yielding the contrast sensitivity function (CSF) [7]. Some typical CSF's are depicted in Figure 1 below shows the typical inverted-U shape of the CSF on logarithmic axes.
+In order to obtain a ground truth comparison, we have implemented the Fourier transform expression for the PSF using in-built DFT packages in python:
-The Robson JG and Graham N [8] have demonstrated that the CSF attains its maximum while using foveal vision as indicated in Fig 2. They observed that the CSF reduces with increasing retinal eccentricity. Further, the decline of the CSF with eccentricity is greater for higher spatial frequencies  but more gradual for lower spatial frequencies [9].
+<math> U(x,y; f) = \frac{1}{\pi}\int_{\kappa_x^2 + \kappa_y^2 \le 1} \int e^{if(\kappa_x^2 + \kappa_y^2)} e^{i\Phi(\kappa_x, \kappa_y)} e^{[2\pi i \kappa_x x + 2\pi i \kappa_y y] } d\kappa_x d\kappa_y </math>
-[[File:Fig2_okkeun.png]]
+The phase and resulting PSF are shown in Fig. 2 below. This was for a case with both astigmatism, coma and defocus. The amount of defocus is 3 normalized units; astigmatism and coma are same as in Fig. 1.
-Fig 2. The plot of CSF as function of Log Spatial Frequency for various spatial locations in the retina.
+[[File:Fig2_indrasen.png]]
-== Methods ==
+Fig. 2: Ground truth PSF calculation using the DFT approach. Top panel: pupil phase and 1D cuts, Bottom panel: resulting intensity PSF in real space and 1D cuts
-In order to determine the CSF for various spatial frequency, we have used the readily available tool from the MATALB called The Image System Engineering Toolbox for Biology (ISETBio). The ISETBio is available as an open source project [10]. The toolbox is primarily used to design image sensors as well as for modeling image formation in biological systems and thus used widely for commercial and research purposes [11].
+In order to verify the analytical formula based on the special functions, we have implemented this and compared it with the ground truth in Fig. 2. For the same defocus and aberration parameters, the resulting point spread function image is shown below in Fig. 3. The difference with the ground truth, normalized to the maximum of the ground truth PSF is shown in Fig. 4. As observed, the intensity PSF resulting from the analytical approach is a good approximation for the DFT based calculation. The error is on the order of 5% close to the peak of the PSF. This is a good sanity check to show that the analytical approximation is accurate enough for us to continue with this approach.
-We have used ISETBio first to simulate the cone photoreceptor absorptions, which is the first stage processing in human vision and thus the CSF to some extent is determined by the photoreceptors. One of the factors that  still remains unclear is  the relation between size of the photoreceptors and the CSF as they are small in the fovea and much larger in the periphery .   In order to include the cone mosaic distribution factors into account, we generated   two portions of the retina- the center of the fovea, and slightly to the periphery using the “coneMosaic” method in ISETBio. We computed the CSF for six   spatial frequencies from 1 to 1.5 in log space. We also varied the integration time in the retina as a factor in order to know how the integration factor affects the CSF. For each of the frequency, we generated a Gabor patch stimulus and obtained the visual response using the toolbox. The stimulus responses between the null and the Gabor patch stimulus was discriminated by using the inference engines from machine learning. In particular, we were interested to compare the performance of both the support vector machines (SVM) and k-nearest neighbours (kNN) algorithm in terms of their role in inferring the CSF.
+[[File:Fig3_indrasen.png]]
-SVM is a supervised learning method developed at AT&T Bell Laboratories by Vapnik with colleagues [12]. The way we used the SVM in our experiment is as follows:  we marked the responses that belong to stimuli and null stimuli categories as one or zero, respectively.  We then trained the SVM using the training samples and tested the model with the test set. We computed the percentage of correct detection of stimuli by the SVM and compare it against the given threshold to see if for the set of front-end optic parameters and the cone mosaics the contrast change is detected. We repeated this for various spatial frequencies and thus obtained the CSF curves.
+Fig. 3: PSF resulting from the closed form analytical solution, using the same aberration and defocus parameters as Fig. 2.
-[[File:Fig3_okkeun.png]]
+[[File:Fig4_indrasen.png]]
-Figure 3. An example of a separable problem in a 2 dimensional space. The support vectors, marked with grey
+Fig. 4: Difference in PSF, normalized by maximum of the PSF. As can be seen, errors on the order of 5% can be observed.
-squares, define the margin of largest separation between the two classes.
+We note that this small error was able to be achieved by terminating the infinite summation in equation (10) of [2] with 10 terms. Therefore, we continue with this choice for the inverse problem solution. Please note that the aberrations specified in Zernike coefficients are in radians (vertical astigmatism <math>(\alpha_2^2 = 0.5)</math> and horizontal coma <math>(\alpha_3^1 = 0.3)</math>). These can be converted to wavelength units, resulting in astigmatism of 79.6 milliwaves, and horizontal coma of 47.7 milliwaves. These are fairly realistic numbers and this case (among other sanity checks) will be used for the rest of this study.
+== Inverse Problem Solution ==
+Here we discuss how to retrieve the pure phase aberrations from a defocus series of PSFs <math>I(r, \phi; f)</math>. In a real world situation, these PSFs would be obtained empirically. However, for the sake of our computational experiments, I have used the Fourier transform method to generate a ground truth PSF progression. Next, we summarize the key steps in retrieving the Zernike coefficients. For a flavor of the derivation, please see section 3 of reference [2]. The key details involve:
+. Expressing the intensity PSF in the analytical form <math>I = |U|^2</math>
+. Inner product of the measured PSFs <math>I_{meas}(r, \phi; f)</math> with the azimuthal sinusoid in order to obtain azimuth-free special functions <math>\Psi^m_{meas}(r,f)</math>
+. Calculation of intensity basis functions <math>\Psi^m_{n}(r,f)</math> from the original special functions <math>V^m_{n}(r,f)</math>
+. For each azimuthal order <math>m</math>, the formulation of a system of equations for obtaining the Zernike coefficients for all valid radial orders corresponding to that azimuthal order: <math>\alpha^m_n, n = m, m+2, m+4, .... m+2M</math>; where <math>M</math> is used to express an appropriately chosen termination of the Zernike polynomial radial order
+. The inversion of the system of equations (by inverting the corresponding Gram matrix) in order to obtain the Zernike coefficients
+The procedure appears to be rather involved, but it can be summarized as follows. The measured intensity PSF is written in terms of unknown Zernike coefficients. Then, the formula is projected into the basis of the azimuthal sinusoids and the intensity basis functions. Once this is done, the quadratic term is rigorously found to be zeroed out, leaving only terms linear in the Zernike coefficients. The linear expression can be solved for the Zernike coefficients just like any other system of linear equations.
+Here we state the formula for the measurement projection (equation 23 in [2]):
+<math>\Psi^m_{meas}(r, f) = \frac{1}{2\pi} \int_0^{2\pi} I_{meas}(r, \phi; f) cos(m\phi) d\phi</math>
+The intensity basis functions can be pre-calculated and are given by (equation 24 in [2]):
+<math>\Psi^m_n(r,f) = -8 \epsilon_m^{-1} Im[i^m V^m_n(r,f) V^{0*}_0(r,f)]</math>
+And each linear equation for solving for the Zernike coefficients is given simply by:
+<math>\Sigma_{n=m, m+2, ... m+2M} \alpha^m_n (\Psi^m_n, \Psi^m_{n'}) = (\Psi^m_{meas}, \Psi^m_{n'})</math>
+The full system of equations can be expressed by vectorizing over the radial order <math>n'</math>, leading to a matrix formulation. There is one matrix for each azimuthal order, and the size of the matrix is [M, M], given by how we originally chose to truncate the Zernike polynomial radial order. By inverting this matrix equation we obtain a vector of size [M,1] containing all the relevant Zernike coefficients for that azimuthal order. This needs to be done for each azimuthal order that we choose to model. I have explained the discretization and truncation choices in the next section.
 == Results ==
-In our experiments, we vary the integration time of the cones and determined the CSF for the fixed location of the cone mosaic, that is at the fovea. The size of the cones mosaic is set at 149 X 149 with each pixel 2X2 um2 size.  Figures 4 and 5 indicate the effect of integration time on the CSF. We observe from the Fig. 4  that we obtained the desired result for all the integration times except 60 msec for the kNN. That is, the CSF declines towards the higher spatial frequency as expected ( the overshoots are the spline interpolation artefacts and should be ignored). On the other hand, using the SVM inference engine, we were unable to detect such pattern except for the integration times 10 and 20 msec.  We suspect that the potential problem could be some bugs in the use of SVM classifier or the contrast detection. We are in the process of debugging these results. Also, we could not see the exact role of integration times just from the Fig 5. That is, while in kNN it is clear that as the integration time increases the CSF increases and such pattern is not witness in the SVM graph, which we also would like to improve in the future.
+We have implemented the ENZ inverse calculation on a small scale problem and compared with ground truth. The following three cases were tested:
+Case 1: Sanity check with no aberrations, all Zernike coefficients expected to be zero
-[[File:Fig5_okkeun_kNN.png]]
+Case 2: Only vertical astigmatism (<math>\alpha_2^2</math>) is present
-Fig. 4.  Plot of CSF vs. spatial frequency for various integration times using kNN classifier
+Case 3: Vertical astigmatism and horizontal coma (<math>\alpha^1_3</math>) are both present
-[[File:Fig4_okkeun_SVM.png]]
+The results below (with retrieved aberration in shown as red) demonstrate that there is a reasonable qualitative correspondence in the retrieved aberrations and ground truth. However, the quantitative correspondence was found to be lacking. Some of the reasons for this are explained at the end of this section. Please note that the aberrations are specified by both the descriptive name and the (m, n) order of the Zernike coefficient in the comparison graphs below.
-Fig. 5.  Plot of CSF vs. spatial frequency for various integration times using SVM binary classifier
+[[File:Fig5_Indrasen.png|900px]]
-== Conclusions ==
+Fig. 5: Results for Case 1 (sanity check with no aberrations)
-We illustrated that the shape of the CSF is determined by the integration time of the photoreceptors in the retinal ganglion cells.  However, we noted that integration time is one of the key factors that determine the CSF. While obtaining these curves, we assumed the default setting in the ISETBio, such as the fixed front end optics and cone mosaic distribution.
+[[File:Fig6_Indrasen.png|900px]]
-With these fixed settings, we found that the CSF declines strongly as the frequency increases. Interestingly, as the theory suggested the value of the CSF for lower integration time is smaller than that of for the higher integration time. However, using the SVM, currently, we could not support this hypothesis as we have gotten the opposite results.
-== Limitations and Future Work ==
+Fig. 6: Results for Case 2 (only one form of aberration)
-One of the major limitations of  our project  was  that we generated the  cone mosaic pattern  randomly for each  spatial frequency. First, this sets different cone mosaic pattern for each frequency. This is not ideal for a controlled experiment, as one person's eye would only have only one mosaic distribution. Therefore, the results in this study has to be inferred carefully as we did not account the fixed cone mosaic pattern. Second, we implicitly assume that the relative distribution of cone types matters more than the specific location of those cones in the mosaic.
+[[File:Fig7_Indrasen.png|900px]]
-Another limitation of the current appraoch is that we totally ignored the role of rods in shaping the CSF. While the cone distribution matters, we believe that the rod masaic may have their role as well to give a better picutre of the visual system.Thus, we may need to write scripts to include the rod distribution.
+Fig. 7: Results for Case 3 (two dominant aberrations: same as the case shown in the PSF progression)
-Finally, the script used to do the SVM and KNN inference engines for the CSF graphs indicated in this project are  provided in the Appendix  As a future work, we would like to add more functionality metioned in the previous paragraphs.
+Due to time and resource constraints, the simulation settings were rather coarse. For instance, the choice of the quantity M, the number of radial orders per azimuthal order could only be 2. Ideally, this number would be on the order of 10 in order to consider a sufficiently large number of radial orders. The Gram matrix method attempts a best fit of the data to the basis functions included in the projection, up to M. Hence, the results may be poor if M is rather low. Another constraint was the spatial discretization, which could not be too fine due to memory and run time constraints on my personal computer. It would be very interesting to perform a more detailed study exploring the impact of focus and spatial discretization, as well as M on the accuracy of the results.
-== References ==
+== Limitations of Method and Future Work ==
-[1]
+There are some well known limitations of the approach as presented in this write-up (which have been addressed in later publications [3]-[5]):
+. The approach is valid only for pure phase aberrations. There is an extension to phase and amplitude aberrations using complex Zernike coefficients that may be worth investigating.
+. The radiometric correction factor has been neglected in the current formulation.
+. A vectorial formulation is required for higher numerical aperture systems.
+. The effect of finite source extent has been neglected, although it is fairly straightforward to incorporate.
-[2] Westheimer, G. (1960). Modulation thresholds for sinusoidal light distributions on the retina. Journal of Physiology, 152, 67-74
+Ideally, a little more investigation into the accuracy, numerical stability and a few more test cases would have been great - but hopefully this is a good starting point.
-[3] http://www.ssc.education.ed.ac.uk/resources/vi&multi/VIvideo/cntrstsn.html
+== References ==
-[4] Jennings, J. A., & Charman, W. N. (1981). Off-axis image quality in the human eye. Vision Research, 21(4), 445–455.
+[1] M. Born and E. Wolf, Principles of Optics, (Cambridge University Press, Cambridge, 2002)
-[5] Nicolas P. Cottaris, Brian A. Wandell, Fred Rieke, David H. Brainard; A computational observer model of spatial contrast sensitivity: Effects of photocurrent encoding, fixational eye movements, and inference engine. Journal of Vision 2020;20(7):17
+[2] C. Van Der Avoort, J. J. M. Braat, P. Dirksen, A. J. E. M. Janssen, "Aberration retrieval from the intensity point-spread function in the focal region using the extended Nijboer-Zernike approach", J. Mod. Opt., 52, 12, 2005, pp. 1695-1728
-[6] https://foundationsofvision.stanford.edu/chapter-7-pattern-sensitivity/
+[3] A. J. E. M. Janssen, "Extended Nijboer-Zernike approach for the computation of optical point-spread functions", J. Opt. Soc. Am. A, 19, 5, 2002
-[7] http://usd-apps.usd.edu/coglab/CSFIntro.html
+[4] J. J. M. Braat, P. Dirksen, A. J. E. M. Janssen, A. S. van des Nes, "Extended Nijboer-Zernike representation of the vector field in the focal region of an aberrated high-aperture optical system", J. Opt. Soc. Am. A, 20, 12, 2003
-[8] Robson JG, Graham N. Probability summation and regional variation in contrast sensitivity across the visual field. Vision Res. 1981;21:409–18.
+[5] A. J. E. M. Janssen, J. J. M. Braat, P. Dirksen, "On the computation of the Nijboer-Zernike aberration integrals at arbitrary defocus", J. Mod. Opt., 51, 5, 2004, pp. 687-703
-[9] Skalicky S.E. (2016) Contrast Sensitivity. In: Ocular and Visual Physiology. Springer, Singapore. https://doi.org/10.1007/978-981-287-846-5_20.
+[6] P. Dirksen, J. Braat, A. J. E. M. Janssen, "Estimating resist parameters in optical lithography using the extended Nijboer-Zernike theory", J. Microlith., Microfab., Microsyst. 5(1), 013005, 2006
-[10] https://github.com/isetbio/isetbio
+[7] Extended Nijboer-Zernike (ENZ) Analysis and Aberration Retrieval research site (Netherlands group): https://www.nijboerzernike.nl/
-[11] http://isetbio.org/
-[12] Cortes, Corinna; Vapnik, Vladimir N. (1995). "Support-vector networks". Machine Learning. 20 (3): 273–297.
 == Appendix ==
-The codes for this project can be found on this notebook: [[Code_okkeun.zip]]
+The computational implementation (jupyter notebook) and results for this project can be found in this zipped folder: [[File:Indrasen_final.zip]]

Indrasen Bhattacharya: Difference between revisions

Latest revision as of 00:37, 8 December 2020

Contents

Introduction

Notation and Conventions

Analytical Forward Model

Verification of ENZ Forward Model

Inverse Problem Solution

Results

Limitations of Method and Future Work

References

Appendix

Navigation menu

Indrasen Bhattacharya: Difference between revisions

Latest revision as of 00:37, 8 December 2020

Introduction

Notation and Conventions

Analytical Forward Model

Verification of ENZ Forward Model

Inverse Problem Solution

Results

Limitations of Method and Future Work

References

Appendix

Navigation menu

Search