Zernike Polynomials Fitting irregular and non-rotationally - - PDF document

1

Zernike Polynomials

• Fitting irregular and non-rotationally symmetric surfaces
• ver a circular region.
• Atmospheric Turbulence.
• Corneal Topography
• Interferometer measurements.
• Ocular Aberrometry

Background

• The mathematical functions

were originally described by Frits Zernike in 1934.

• They were developed to

describe the diffracted wavefront in phase contrast imaging.

• Zernike won the 1953 Nobel

Prize in Physics for developing Phase Contrast Microscopy.

2

Phase Contrast Microscopy

Transparent specimens leave the amplitude of the illumination virtually unchanged, but introduces a change in phase.

Applications

• Typically used to fit a wavefront or

surface sag over a circular aperture.

• Astronomy - fitting the wavefront entering

a telescope that has been distorted by atmospheric turbulence.

• Diffraction Theory - fitting the wavefront

in the exit pupil of a system and using Fourier transform properties to determine the Point Spread Function.

Source: http://salzgeber.at/astro/moon/seeing.html

3

Applications

• Ophthalmic Optics - fitting corneal

topography and ocular wavefront data.

• Optical Testing - fitting reflected

and transmitted wavefront data measured interferometically.

Surface Fitting

• Reoccurring Theme: Fitting a complex, non-rotationally

symmetric surfaces (phase fronts) over a circular domain.

• Possible goals of fitting a surface:

– Exact fit to measured data points? – Minimize “Error” between fit and data points? – Extract Features from the data?

4

1D Curve Fitting

5 10 15 20 25

• 0.1

0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5

Low-order Polynomial Fit

y = 9.9146x + 2.3839 R

2 = 0.9383 5 10 15 20 25

• 0.1

0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5

In this case, the error is the vertical distance between the line and the data point. The sum of the squares of the error is minimized.

5

High-order Polynomial Fit

5 10 15 20 25

• 0.1

0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5

• 400
• 200

200 400 600 800 1000 1200 1400 1600 1800

• 0.1

0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5

y = a0 + a1x + a2x2 + … a16x16

Cubic Splines

5 10 15 20 25

• 0.1

0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5

Piecewise definition of the function.

6

Fitting Issues

• Know your data. Too many terms in the fit can be numerically

unstable and/or fit noise in the data. Too few terms may miss real trends in the surface.

• Typically want “nice” properties for the fitting function such as

smooth surfaces with continuous derivatives. For example, cubic splines have continuous first and second derivatives.

• Typically want to represent many data points with just a few terms of a
• fit. This gives compression of the data, but leaves some residual error.

For example, the line fit represents 16 data points with two numbers: a slope and an intercept.

Why Zernikes?

• Zernike polynomials have nice mathematical properties.

– They are orthogonal over the continuous unit circle. – All their derivatives are continuous. – They efficiently represent common errors (e.g. coma, spherical aberration) seen in optics. – They form a complete set, meaning that they can represent arbitrarily complex continuous surfaces given enough terms.

7

Orthogonal Functions

• Orthogonal functions are sets of surfaces which have some

nice mathematical properties for surface fitting.

• These functions satisfy the property

∫∫

⎩ ⎨ ⎧ = =

A j j i

Otherwise j i C dxdy ) y , x ( V ) y , x ( V

where Cj is a constant for a given j

Orthogonality and Expansion Coefficients

∑

=

i i i

) y , x ( V a ) y , x ( W ) y , x ( V ) y , x ( V a ) y , x ( V ) y , x ( W

j i i i j

∑

= dxdy ) y , x ( V ) y , x ( V a dxdy ) y , x ( V ) y , x ( W

i j A i i A j

∑ ∫∫ ∫∫

= dxdy ) y , x ( V ) y , x ( W C 1 a

A j j j

∫∫

=

Linear Expansion

8

Orthogonality - 1D Example

( ) ( )

∫

π 2

integers are m' and m where dx x ' m cos mx sin

( ) ( ) [ ]

∫

π

− + + =

2

dx x ' m m sin x ' m m sin 2 1

Consider the integral

( ) ( )

π = =

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − + + + − =

2 x x

' m m x ' m m cos ' m m x ' m m cos 2 1 =

Orthogonality - 1D Example

( ) ( )

∫

π 2

integers are m' and m where dx x ' m sin mx sin

( ) ( ) [ ]

∫

π

− − + − =

2

dx x ' m m cos x ' m m cos 2 1

Two sine terms

( ) ( )

π = =

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − − + + − =

2 x x

' m m x ' m m sin ' m m x ' m m sin 2 1 =

if m ≠ m’ !!!

9

Orthogonality - 1D Example

( ) ( )

∫

π 2

integers are m' and m where dx x ' m sin mx sin

( )

∫

π

=

2 2

dx mx sin

Two sine terms with m = m’

( )

π = =

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − =

2 x x

m 4 mx 2 sin 2 x π =

Similar arguments for two cosine terms

Orthogonality - 1D Example

• dd

j even j x 2 1 j sin x 2 j cos Vj ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ =

From the previous arguments, we can define

( ) ( )

π = = ∫

∫

π π

2 dx dx x V x V

2 2

Note that when j = 0, V0 = 1 and so the constant Cj = 2π for j = 0, and Cj = π for all other values of j.

10

Orthogonality - 1D Example

The functions satisfy the orthogonality condition over [0, 2π]

( ) ( ) ( )

Otherwise j i 1 dx x V x V

i 2 j i

= ⎩ ⎨ ⎧ π δ + =

∫

π

where δ is the Kronecker delta function defined as

Otherwise j i 1

ij

= ⎩ ⎨ ⎧ = δ

Fit to W(x) = x

10

• 0.4

9 8

• 0.5

7 6

• 0.6666

5 4

• 1

3 2

• 2

1 π aj j

• 1

1 2 3 4 5 6 7 1 2 3 4 5 6 7 W(x) 1 term 3 terms 11 terms

11

Extension to Two Dimensions

In many cases, wavefronts take on a complex shape defined over a circular region and we wish to fit this surface to a series of simpler components.

Wavefront Fitting

=

• 0.003 x

+ 0.002 x + 0.001 x

12

Unit Circle

x y ρ θ 1 Divide the real radial coordinate by the maximum radius to get a normalized coordinate ρ

Orthogonal Functions on the Unit Circle

• Taylor polynomials (i.e. 1, x, y, x2, xy, y2,….) are not
• rthogonal on the unit circle.

⎩ ⎨ ⎧ = = θ ρ ρ θ ρ θ ρ

∫∫

π

Otherwise j i C d d ) , ( V ) , ( V

j j 2 1 i

where Cj is a constant for a given j Many solutions, but let’s try something with the form

( ) ( )

θ Θ ρ = θ ρ

i i i

R ) , ( V

13

Orthogonal Functions on the Unit Circle

• The orthogonality condition separates into the product of

two 1D integrals

⎩ ⎨ ⎧ = = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ θ θ Θ θ Θ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ρ ρ ρ ρ

∫ ∫

π

Otherwise j i C d ) ( ) ( d ) ( R ) ( R

j j 2 i j 1 i

This looks like the 1D Example, so sin(mθ) and cos(mθ) a possibility This has extra ρ in it, so we need different functions

ANSI Standard Zernikes

⎪ ⎩ ⎪ ⎨ ⎧ < θ ρ − ≥ θ ρ = θ ρ m for ; m sin ) ( R N m for ; m cos ) ( R N ) , ( Z

m n m n m n m n m n

Double Index n is radial order m is azimuthal frequency Normalization Radial Component Azimuthal Component

ANSI Z80.28-2004 Methods for Reporting Optical Aberrations of Eyes.

14

ANSI Standard Zernikes

• nly depends
• n |m| (i.e. same

for both sine & cosine terms Powers of ρ

[ ] [ ]

s 2 n 2 / ) m n ( s s m n

! s ) m n ( 5 . ! s ) m n ( 5 . ! s )! s n ( ) 1 ( ) ( R

− − =

ρ − − − + − − = ρ

∑

Constant that depends

• n n and m

2 n 2 d ) ( R ) ( R

' nn m ' n 1 m n

+ δ = ρ ρ ρ ρ

∫

The Radial polynomials satisfy the orthogonality equation

ANSI Standard Zernikes

constant that depends on n & m

m m n

1 2 n 2 N δ + + =

15

Orthogonality

⎩ ⎨ ⎧ = = = θ ρ ρ θ ρ θ ρ

∫∫

π

Otherwise m' m ; n' n for C d d ) , ( Z ) , ( Z

m n, m n 2 1 ' m ' n

( ) ( ) ( ) ( )

⎩ ⎨ ⎧ = = = ρ ρ ρ ρ × θ ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ θ − θ ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ θ − θ

∫ ∫

π

Otherwise m' m ; n' n for C d ) ( R ) ( R d m sin m cos ' m sin ' m cos N N

m n, m n 1 ' m ' n 2 m n ' m ' n

Orthogonality

( ) [ ]

⎩ ⎨ ⎧ = = = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + δ π δ + δ Otherwise m' m ; n' n for C 2 n 2 1 N N

m n, n ' n m m ' m m n ' m ' n

This equals zero unless m=m’ This equals zero unless n=n’ So this portion is satisfied What happens when n=n’ and m=m’?

16

Orthogonality

( ) (

) [ ]

m n, m 2 m n

C 2 n 2 1 1 N = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + π δ +

When n=n’ and m=m’

( ) [ ]

m , n m m0

C 2 n 2 1 1 1 2 2n = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + π δ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ δ + + π =

m , n

C

Orthogonality

m ' m n ' n m n 2 1 ' m ' n

d d ) , ( Z ) , ( Z δ πδ = θ ρ ρ θ ρ θ ρ

∫∫

π

17

The First Few Zernike Polynomials

( ) ( ) ( ) ( )

( )

( )

( )

( )

( )

θ ρ = θ ρ − ρ = θ ρ θ ρ = θ ρ θ ρ = θ ρ θ ρ = θ ρ = θ ρ

− −

2 cos 6 , Z 1 2 3 , Z 2 sin 6 , Z cos , Z sin , Z 1 , Z

2 2 2 2 2 2 2 2 1 1 1 1

Zernike Polynomials

Azimuthal Frequency, θ

Radial Polynomial, ρ

Z0 Z

1 1

Z

11 −

Z2 Z31

−

Z3

1

Z4 Z4

2

Z2

2

Z42

−

Z33

−

Z3

3

Z4

4

Z44

−

Z22

−

18

Caveats to the Definition of Zernike Polynomials

• At least six different schemes exist for the Zernike polynomials.
• Some schemes only use a single index number instead of n and m.

With the single number, there is no unique ordering or definition for the polynomials, so different orderings are used.

• Some schemes set the normalization to unity for all polynomials.
• Some schemes measure the polar angle in the clockwise direction from

the y axis.

• The expansion coefficients depend on pupil size, so the maximum

radius used must be given.

• Some groups fit OPD, other groups fit Wavefront Error.
• Make sure which set is being given for a specific application.

Another Coordinate System

x y ρ φ Normalized Polar Coordinates:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = φ = ρ

−

y x tan r r

1 max

1 ρ ranges from [0, 1] φ ranges from [-180°, 180°]

NON- STANDARD

19

Zernike Polynomials - Single Index

Azimuthal Frequency, θ

Radial Polynomial, ρ Z0 Z1 Z4 Z5 Z3 Z9 Z8 Z7 Z6 Z10 Z11 Z12 Z13 Z14 Z2

ANSI STANDARD

Starts at 0 Left-to-Right Top-to-Bottom

Other Single Index Schemes

Z1 Z3 Z4 Z6 Z5 Z10 Z8 Z7 Z9 Z15 Z13 Z11 Z12 Z14 Z2

NON- STANDARD

Starts at 1 cosines are even terms sines are odd terms

Noll, RJ. Zernike polynomials and atmospheric turbulence. J Opt Soc Am 66; 207-211 (1976). Also Zemax “Standard Zernike Coefficients”

20

Other Single Index Schemes

Z1 Z3 Z4 Z5 Z6 Z10 Z7 Z8 Z11 Z18 Z13 Z9 Z12 Z17 Z2

NON- STANDARD

Starts at 1 increases along diagonal cosine terms first 35 terms plus two extra spherical aberration terms No Normalization!!!

Zemax “Zernike Fringe Coefficients” Also, Air Force or University of Arizona

Other Single Index Schemes

• Born & Wolf
• Malacara
• Others??? Plus mixtures of non-normalized, coordinate

systems.

NON- STANDARD

Use two indices n, m to unambiguously define polynomials. Use a single standard index if needed to avoid confusion.

21

Examples

Example 1: 4m 0.25 D of myopia for a 4 mm pupil (rmax = 2 mm)

( ) ( ) ( )

θ ρ + θ ρ = ρ = ρ = = , Z 3 4000 1 , Z 4000 1 2000 8000 2 8000 r W

2 2 2 2

4mm

Examples

Example 2: 1m 1.00 D of myopia for a 2 mm pupil (rmax = 1 mm)

( ) ( )

θ ρ + θ ρ = ρ = = , Z 3 4000 1 , Z 4000 1 2000 2000 r W

2 2 2

2mm Same Zernike Expansion as Example 1, but different rmax. Always need to give pupil size with Zernike coefficients!!

22

RMS Wavefront Error

• RMS Wavefront Error is defined as

( )

∑ ∫∫ ∫∫

>

= φ ρ ρ φ ρ ρ φ ρ =

m all , 1 n 2 m , n 2 WFE

a d d d d ) ) , ( W ( RMS

Zeroth Order Zernike Polynomials

Z0

This term is called Piston and is usually ignored. The surface is constant over the entire circle, so no error or variance exists.

23

First Order Zernike Polynomials

Z1

1

Z11

−

These terms represent a tilt in the wavefront.

max 11 max 1 1 11 1 1 1 1 11 1 1 1 1

r x a r y a cos a sin a ) , ( Z a ) , ( Z a + = θ ρ + θ ρ = θ ρ + θ ρ

− − − −

Combining these terms results in a general equation for a plane, thus by changing the coefficients, a plane at any orientation can be created. This rotation of the pattern is true for the sine/cosine pairs of Zernikes

Second Order Zernike Polynomials

Z2 Z2

2

Z22

− These wavefronts are what you would expect from Jackson crossed cylinder J0 and J45 and a spherical lens. Thus, combining these terms gives any arbitrary spherocylindrical refractive error.

24

Third Order Zernike Polynomials

Z31

−

Z3

1

Z33

−

Z3

3

The inner two terms are coma and the outer two terms are trefoil. These terms represent asymmetric aberrations that cannot be corrected with convention spectacles or contact lenses.

Fourth Order Zernike Polynomials

Z4 Z4

2

Z42

−

Z4

4

Z44

−

Spherical Aberration 4th order Astigmatism 4th order Astigmatism Quadroil Quadroil These terms represent more complex shapes of the wavefront. Spherical aberration can be corrected by aspheric lenses.

25

Discrete data

• Up to this point, the data has been continuous, so we can

mathematically integrate functions to get expansion coefficients.

• Real-world data is sampled at discrete points.
• The Zernike polynomials are not orthogonal for discrete

points, but for high sampling densities they are almost

• rthogonal.

Speed

• The long part of calculating Zernike polynomials is

calculating factorial functions.

⎪ ⎩ ⎪ ⎨ ⎧ < θ ρ − ≥ θ ρ = θ ρ m for ; m sin ) ( R N m for ; m cos ) ( R N ) , ( Z

m n m n m n m n m n

[ ] [ ]

s 2 n 2 / ) m n ( s s m n

! s ) m n ( 5 . ! s ) m n ( 5 . ! s )! s n ( ) 1 ( ) ( R

− − =

ρ − − − + − − = ρ

∑

26

Speed

• Chong et al.* developed a recurrence relationship that

avoids the need for calculating the factorials.

• The results give a blazing fast algorithm for calculating

Zernike expansion coefficents using orthogonality.

*Pattern Recognition, 36;731-742 (2003).

Least Squares Fit

⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛

− − − −

) y , x ( f ) y , x ( f ) y , x ( f a a a a ) y , x ( Z ) y , x ( Z ) y , x ( Z ) y , x ( Z ) y , x ( Z ) y , x ( Z ) y , x ( Z ) y , x ( Z ) y , x ( Z ) y , x ( Z ) y , x ( Z ) y , x ( Z

N N 2 2 1 1 nm 11 1 1 00 N N m n N N 1 1 N N 1 1 N N 2 2 m n 2 2 1 1 2 2 1 1 2 2 1 1 m n 1 1 1 1 1 1 1 1 1 1

M M L M L M M M L L

( )

F Z Z Z A F Z ZA Z F ZA

T 1 T T T −

= = =

27

Gram-Schmidt Orthogonalization

• Examines set of discrete data and creates a series of

functions which are orthogonal over the data set.

• Orthogonality is used to calculate expansion coefficients.
• These surfaces can then be converted to a standard set of

surfaces such as Zernike polynomials. Advantages

• Numerically stable, especially for low sampling density.

Disadvantages

• Can be slow for high-order fits
• Orthogonal functions depend upon data set, so a new set

needs to be calculated for every fit.

Elevation Fit Comparison

0.125 Seconds Chong Algorithm 10 Seconds Gram-Schmidt

32 Orders or 560 total polynomials

28

Shack-Hartmann Wavefront Sensor

Plane Wavefronts Eye Lenslet Array Aberrated Wavefronts Eye Lenslet Array Perfect wavefronts give a uniform grid of points, whereas aberrated wavefronts distort the grid pattern.

Least Squares Fit

( )

F Z Z Z A F Z ZA Z F ZA

T 1 T T T −

= = =

Again, conceptually easy to understand, although this can be relatively slow for high order fits.

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = dy / ) y , x ( dW dy / ) y , x ( dW dy / ) y , x ( dW dx / ) y , x ( dW dx / ) y , x ( dW dx / ) y , x ( dW a a a dy / ) y , x ( dV dy / ) y , x ( dV dy / ) y , x ( dV dy / ) y , x ( dV dy / ) y , x ( dV dy / ) y , x ( dV dy / ) y , x ( dV dy / ) y , x ( dV dy / ) y , x ( dV dx / ) y , x ( dV dx / ) y , x ( dV dx / ) y , x ( dV dx / ) y , x ( dV dx / ) y , x ( dV dx / ) y , x ( dV dx / ) y , x ( dV dx / ) y , x ( dV dx / ) y , x ( dV Z

N N 2 2 1 1 N N 2 2 1 1 J 2 1 N N J N N 2 N N 1 2 2 J 2 2 2 2 2 1 1 1 J 1 1 2 1 1 1 N N J N N 2 N N 1 2 2 J 2 2 2 2 2 1 1 1 J 1 1 2 1 1 1

M M M L M M M M L L L M M M M L L

29

Example Image Wavefront Reconstruction

PSF