Q2d_nm_c_to_a_b needs a non-zero positive and negative m to work

import numpy as np
from prysm.polynomials.qpoly import compute_z_zprime_Q2d, Q2d_nm_c_to_a_b, Q2d
from prysm.polynomials import fringe_to_nm, nm_to_name
from prysm import geometry, coordinates
import matplotlib.pyplot as plt
from matplotlib import colors

def plot_q2d(fringes, coeffs):

    x, y = coordinates.make_xy_grid(2048, diameter=2)
    r, t = coordinates.cart_to_polar(x, y)

    z = np.zeros((2048, 2048))
    dr = np.zeros((2048, 2048))
    dt = np.zeros((2048, 2048))

    nms = []
    for J, C in zip(fringes, coeffs):
        (n, m) = fringe_to_nm(J)
        print(f"(n: {n}, m: {m}) - {C}")
        nms.append((n, m))

    (cms, ams, bms )= Q2d_nm_c_to_a_b(nms, coeffs)
    (_z, _dr, _dt) = compute_z_zprime_Q2d(cms, ams, bms, r, t)

    z  = (_z)
    dr = (_dr)
    dt = (_dt)

    # Mask out a circle for plotting
    mask = geometry.circle(1, r)

    z[mask!=1]=np.nan
    dr[mask!=1]=np.nan
    dt[mask!=1]=np.nan

    #  Plotting
    width = 20
    fig = plt.figure(figsize=(width,width/3))

    color_norm=colors.TwoSlopeNorm(vcenter=0.0)
    z_ax = fig.add_subplot(1, 3, 1)
    z_ax.imshow(z, norm=color_norm, cmap="RdBu", extent=[-1, 1, -1, 1])
    z_ax.set_title("Phase")
    fig.colorbar(z_ax.get_images()[0])

    color_norm=colors.TwoSlopeNorm(vcenter=0.0)
    dr_ax = fig.add_subplot(1, 3, 2)
    dr_ax.imshow(dr, norm=color_norm, cmap="RdBu", extent=[-1, 1, -1, 1])
    dr_ax.set_title("Radial derivative of phase")
    fig.colorbar(dr_ax.get_images()[0])

    color_norm=colors.TwoSlopeNorm(vcenter=0.0)
    dt_ax = fig.add_subplot(1, 3, 3)

    dt_ax.imshow(dt, norm=color_norm, cmap="RdBu", extent=[-1, 1, -1, 1])
    dt_ax.set_title("Tangential derivative of phase")
    fig.colorbar(dt_ax.get_images()[0])

    fig.suptitle("Q2D stuff")
    plt.show()
    plt.close(fig)

zernikes = np.arange(1, 6)
coeffs = np.zeros(zernikes.shape)
#coeffs[5] = 1
#coeffs[4] = 0.5
coeffs[1] = 10

plot_q2d(zernikes, coeffs)

So I've been trying to wrap my head around all things Q2d, and so far compute_z_z_prime seems to the be function I need. Right now the Q2d function will take any n and m values in terms of zernike-esque indices. However trying to do that with a combination of Q2d_nm_c_to_a_b and compute_z_zprime_Q2d you run into errors if negative and positive m terms don't match.

(Ignore the variable name and fact this is based off of zernike code for a sec)

In this example if I dozernikes = np.arange(1, 3), which corresponds to the following nms in zernike indexing: (n: 0, m: 0) (n: 1, m: 1)

It will error in Q2d_nm_c_to_a_b with:

Traceback (most recent call last):
  File "c:\Users\alist\Documents\GitHub\SympleyWavey\simulations\q2d_polynomials_sanity.py", line 84, in <module>
    plot_q2d(zernikes, coeffs)
  File "c:\Users\alist\Documents\GitHub\SympleyWavey\simulations\q2d_polynomials_sanity.py", line 24, in plot_q2d
    (cms, ams, bms )= Q2d_nm_c_to_a_b(nms, coeffs)
  File "C:\Users\alist\miniconda3\envs\prysm\lib\site-packages\prysm\polynomials\qpoly.py", line 1257, in Q2d_nm_c_to_a_b
    max_m_b = max(list(bc.keys()))
ValueError: max() arg is an empty sequence

Now this obviously isn't a very imposing issue to work around, but it is one I had to peak at the source of qpoly.py to figure out, and the docstrings don't list the restriction.

So it looks like either a case the function should handle or have listed as a restriction?

Thanks for the bug report. A mild tangent to the bug -- the way prysm computes all of its polynomials is derived from the teachings of Greg Forbes, by using recurrence relations. That is in contrast to how, say, POPPY computes Zernikes using Rodrigues' form. The recurrence relations are generally straightforward to differentiate, and I wanted to make a raytracer at the end of 2021, so there be derivatives.

Forbes is truly a genius, and buried in the appendix of his papers is a method for computing not only the surface sag, but also its derivatives, in a single pass without ever computing either the polynomials or their derivatives directly. I actually asked him about computing the derivatives directly over e-mail in 2018 or 2019, and he said he never bothered because for raytracing or surface fitting, you don't actually ever have to do that.

Since I am very bad at math and the creator never wrote down how to do do it, there is not a Q2D_der function.

If you pass a single n/m to compute_z_zprime, it will fail because it needs all of the coefficients from n=0..N, even if the lower several are zero. This is a contrivance of the clever math Forbes used, and I typed into the computer.

The reason your code is crashing is perhaps hidden in the docstring of compute_z_zprime_Q2d,

    cm0 : iterable
        surface coefficients when m=0 (inside curly brace, top line, Eq. B.1)
        span n=0 .. len(cms)-1 and mus tbe fully dense
    ams : iterable of iterables
        ams[0] are the coefficients for the m=1 cosine terms,
        ams[1] for the m=2 cosines, and so on.  Same order n rules as cm0
    bms : iterable of iterables
        same as ams, but for the sine terms
        ams and bms must be the same length - that is, if an azimuthal order m
        is presnet in ams, it must be present in bms.  The azimuthal orders
        need not have equal radial expansions.
        For example, if ams extends to m=3, then bms must reach m=3
        but, if the ams for m=3 span n=0..5, it is OK for the bms to span n=0..3,
        or any other value, even just [0].

When you use range(1,3) you are doing Z1 (piston) and Z2 (tip). This will lead to cm being a length 1 list, ams would be a length 1 list, and bms a length 0 list.

The Q2D routines could be modified to remove this constraint, but it seems to only manifest in teaching contexts. Most real usages will expand as high as, say, 3 terms in each c, a, and b.

Since I am very bad at math and the creator never wrote down how to do do it, there is not a Q2D_der function.

Right there with you with you only struggling through these maths heavy optics papers, really thinking I should have paid more attention in my vector calc classes... Interesting to hear that he never bothered to write the derivatives on their own, based especially based on other papers for orthogonal zernike derivatives alone being published (eg Vector polynomials orthogonal to the gradient of Zernike polynomial, A. Gavrielides, 1982 or Determination of phase mode components in terms of local wave-front slopes: an analytical approach, E. Acosta, 1995, and many more I've found on my shackhartmann adventures).

But I guess I'll give him some slack given the combined wavefront and orthogonal derivatives computation, given reconstructing a wavefront is generally what derivatives are intended to be used for(?) And of course for coming up with such a useful set :)

Agreed on the compute_z_zprime_Q2d restriction - there's no top level prysm wrapper that puts the output of Q2d_nm_c_to_a_b straight into compute_z_zprime_Q2d, so the first returning invalid data for the latter is the user's problem not prysm's.

This is probably my fault for starting to write out something for that case however before changing my mind. But the part which actually seems like a bug is Q2d_nm_c_to_a_b in isolation, with the same caveat applies of only likely to appear in trivial learning examples - but the only restriction listed in that docstring is that coefs and nms must be the same length. However if you don't put a negative 'm' in you get the following error:

(prysm) C:\Users\alist>python
Python 3.10.4 | packaged by conda-forge | (main, Mar 30 2022, 08:38:02) [MSC v.1916 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> from prysm.polynomials.qpoly import Q2d_nm_c_to_a_b
>>> Q2d_nm_c_to_a_b(nms=[(0,0), (1,1)], coefs=[1,1])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\alist\miniconda3\envs\prysm\lib\site-packages\prysm\polynomials\qpoly.py", line 1257, in Q2d_nm_c_to_a_b
    max_m_b = max(list(bc.keys()))
ValueError: max() arg is an empty sequence
>>>

(and likewise if you only have negative m it will error in the same way on max_m_a = max(list(ac.keys())) instead

[...] But I guess I'll give him some slack given the combined wavefront and orthogonal derivatives computation, given reconstructing a wavefront is generally what derivatives are intended to be used for(?) And of course for coming up with such a useful set :)

I think the predominant use for derivatives in optics is for raytracing, where only the surface sag and its derivatives matter, and the individual basis functions have no intrinsic meaning on their own. It is true that there are a number of wavefront sensors (ok, maybe half or so of them) that fundamentally sense gradients -- shack hartmanns, modified hartmann masks ("quadriwave lateral shearing interferometer"), pyramid WFS, etc. But there exist zonal reconstruction algorithms for those sensors, which AFAIK are the predominant methods used.

Anyway.... the "integration" of derivatives in prysm is fairly poor. I wrote all of that code as part of the experimental/raytracing module circa, well, the last time I touched the raytracer. And, the Q2Ds are very dissimilar to Zernikes in how they are structured, with no monoindexing scheme defined that I am aware of (Noll, or ANSI j, or Fringe, etc, being monoindexing schemes).

compute_z_zprime_Q2D expects one vector and two matricies (though the matricies could be ragged lists of lists), which count coefficients from n=0 at a given m, where the row index in the matrix is m.

For example,

import numpy as np

from prysm import (
    coordinates,
    geometry,
    polynomials,
)

from matplotlib import pyplot as plt

x, y = coordinates.make_xy_grid(256, diameter=2)
r, t = coordinates.cart_to_polar(x, y)

mask = geometry.circle(1, r)

# Q spherical n=0, n=1, n=2 all with coef of 1
cms = [1,1,1]

# Q coma of orders 0, 1, 2 with coef 1
# and Q ast ...
# and Q trefoil
ams = np.ones((3,3))
bms = ams.copy()

z, dr, dt = polynomials.qpoly.compute_z_zprime_Q2d(cms, ams, bms, r, t)
z[~mask]=np.nan
dr[~mask]=np.nan
dt[~mask]=np.nan
fig, axs = plt.subplots(ncols=3)

axs[0].imshow(z)
axs[1].imshow(dr)
axs[2].imshow(dt)

If you want to extract the derivative w.r.t. a specific term, just set all of the elements of cms/ams/bms to zero except that one element.

For example, to do this for the second-order astigmatism,

cms = []
ams = np.zeros((3,3))
ams[1,1] = 1
bms = np.zeros((3,3))

Technically bms could be [[], [], []] and the code would run slightly faster, but it's pretty zippy as-is, and the syntax there is obtuse.

Realistically, the Q2d_nm_c_to_a_b is probably a mistake, since for example n=4,m=0 is spherical in Zernikes, but it's n=0,m=0 that is most similar in Qs. And for "primary astigmatism" with Zernike it is n=2, m=2, while for Qs it is n=0, m=2. Trefoil n=m=3 vs n=0,m=3, and so on.

I guess it could still be removed, since prysm is <= v1.0...

This is certainly an area where pull requests (especially to enhance docs!) would be appreciated :)

Realistically, the Q2d_nm_c_to_a_b is probably a mistake, since for example n=4,m=0 is spherical in Zernikes, but it's n=0,m=0 that is most similar in Qs. And for "primary astigmatism" with Zernike it is n=2, m=2, while for Qs it is n=0, m=2. Trefoil n=m=3 vs n=0,m=3, and so on.

I guess it could still be removed, since prysm is <= v1.0...

I've gone and reread the original forbes paper and see what you mean now my brain has engaged a bit more. I guess given the matrix form there's no sensible way to get to a monomial indexing at all? Each could just extend off forever, besides the "m + 2n ≤ T" relations described just above eq 2.3.

It might take me a while, but I'll see if I can come up with a bit of example code that uses a set of modes from that basis to say fit a random wavefront? Specifically using the nomenclature in the paper where its just the A term's matrix that is used. Could be a nice little learning exercise (for both myself and anyone coming across the lib) then maybe I'll be able to contribue something more useful to the code/api docs themself.

The matrix is just a dense way of representing coefficients[m,n]. Strictly speaking, one could also give coefficients as a 3-tuple (n,m,coef) -- that is precisely what Q2d_nm_c_to_a_b is for!

The reason the code uses this matrix to represent the coefficients is firstly just because that's how Forbes presented it in the paper, which is pretty trite. But the "core" reason is that the mathematical formalism Forbes teaches to compute orthogonal polynomials in general, but especially for Qs, starts from, say, m=0,n=0, and uses that along the way to computing m=0,n=1, then uses that to compute m=0,n=2, and so on.

It would be completely equivalent to make a Zernike interface in the same formalism, in fact zernike_nm_sequence does exactly that here

Forbes does have a paper on fitting Q polynomials. It uses yet more exotic mathematics and not least squares, and I have not implemented it (since you can just use least squares).

I don't have any knowledge of fitting based on gradient measurements -- you might send an email to Forbes about that; he probably would be interested in writing a paper on the topic.

brandondube / prysm

Q2d_nm_c_to_a_b needs a non-zero positive and negative m to work #87