# Efficient and robust recurrence relations for the Zernike circle polynomials and their derivatives in Cartesian coordinates Algorithm

Efficient and robust recurrence relations for the Zernike circle polynomials and their derivatives in Cartesian coordinates

Background
There was a recent paper published covering how to efficiently recursively generate zernike polynomial gradients in cartesian coordinates (paper here) [1], written by Torben Andersen. I first want to say thank you to Mr. Andersen for this excellent paper and contribution to the optics community.

I will try to post review of this paper in the near future, as I was pretty impressed with the paper itself, and I strongly encourage interested readers to go and check it out for themselves (it is open source, so anyone can access it). However, the point of this post is the code and algorithm. We utilize the algorithm from this paper to generate our zernike gradient values, which we use to define normal vectors for our zernike surfaces (note, we utilize the jacobi polynomial method, posted here, to generate our Zernike polynomials).

I found the puedo code provided to be more trouble than helpful when trying to utilize the algorithm in the paper, so instead we simply rewrote it from scratch from the excellent derivations provided in the paper.

Testing Method
The zernike values generated have been cross checked against hand calculated zernike value results for the first 4 radial order of Zernike polynomials, as well as against the results from Prysm (if you havent checked out Prysm, please do yourself a favor and give it a look. It is an incredibly powerful and well documented open source optics code base, written by @bdube ).

The gradient values for the first 4 radial orders were similarly checked against hand calculated values and matched to 1E-8 accuracy. I am happy to discuss the test cases in more detail if desired, but I have not posted the tests I used in the code (it’s long enough already!).

Languages Written
The algorithm has been written in python, and in rust, both of which are provided at the end of this post, so feel free to use them!

Feedback
I am interested to hear if you see any errors in the code or how you would write this more efficiently. I don’t view the code I wrote as ‘well written’, it was more just a way to get the algorithm down and working, so I am excited to see what improvements you all suggest, and I hope in general it is helpful to everyone who may need to generate Zernike and Zernike polynomial gradients. @bdube, @Salsbury , @hkang, @henryquach, and @tbrendel in particular curious to hear your thoughts on the code.

References
[1] Torben B. Andersen, “Efficient and robust recurrence relations for the Zernike circle polynomials and their derivatives in Cartesian coordinates,” Opt. Express 26, 18878-18896 (2018)

Code
Note, for both the Rust and Python code, the ‘order’ or ‘index’ refers to the Ansi standard indexing, described here.

Both codes below are written under the MIT license, feel free to use however you would like, you do not need to reference me/this post for where you got the code (although I do hope you include the reference to the original paper!).

Rust Code

[dependencies]
nalgebra = "0.27"


Main Code

use nalgebra as na;
use std::collections::HashMap;

fn main() {
// Define an example of input coordinates
let coordinates = vec![
na::Point2::new(0.0, 1.0),
na::Point2::new(0.0, 0.5),
na::Point2::new(0.0, 0.0),
na::Point2::new(0.0, -0.5),
na::Point2::new(0.0, -1.0),
];

// Define an example maximum radial order, here 4, which corresponds to
// a maximum ansi index for the zernike calculated of 14
let max_order = 4;

// Then, calculate our Zernike function, and gradients of it.

// The output is not ordered, I am sure there is some clever rust way to avoid this initially, but I don't know it, thus we sort
zernike_values.sort_by(|a, b| Ord::cmp(&a.0.j, &b.0.j));

// And lastly, print out our results (if you want, there is no real purpose for printing it out however besides cluttering your terminal)
print!("The zernike val, dx, dy is:");
for index in 0..zernike_values.len() {
for pos in 0..coordinates.len() {
print!(
"z: {}, dz/dx: {}, dz/dy: {} at (x: {}, y: {}) for ansi index: {}, ",
zernike_values[index].1[pos],
coordinates[pos].x,
coordinates[pos].y,
zernike_values[index].0.j,
);
}
}
}

#[derive(Debug, Clone)]
pub struct Ansi {
j: u16,
}

coordinates: &[na::Point2<f64>],
mut max_order: u16,
) -> (
Vec<(Ansi, Vec<f64>)>,
Vec<(Ansi, Vec<f64>)>,
Vec<(Ansi, Vec<f64>)>,
) {
//   Calculates recursively the Zernike and Gradient polynomialt sets
//  The zernike polynomial gradients are ouput as a vector of the ansi order
// and the evaluated gradient value at the provided input coordinates
// This algorithm is taken from the following publication:
// Torben B. Andersen, "Efficient and robust recurrence relations for the Zernike circle polynomials and their derivatives in Cartesian coordinates,"
// Opt. Express 26, 18878-18896 (2018)
// Author : Logan Rodriguez Graves. Date: 1/27/022

//
// Inputs:
// coordinates: (x,y) point coordinates of where you would like to evaluate your Zernike, and Gradients functions
// max_order: maximum radial order to recursively generate the Zernike, and Zernike Gradient Functions to.
// Outputs a tuple of vectors of tuples: (Vec<(Ansi Order, Vec<zernike values>)>, Vec<(Ansi Order, Vec<x gradient values>)>, Vec<(Ansi Order, Vec<y gradient values>)>)
// where each function has been evaluated at the input coordinate points for each ansi order

// Determine the max radial order, if it is less than 2 then set it to at least 2.
if max_order < 2 {
max_order = 2;
}

// Hashmaps are your friend here! They will allow us to directly reference into
// the sets we need, as opposed to an unnecessarily extra complex step of translating nm
// indices to a single order.
let mut zernike_map: HashMap<(u16, u16), (Ansi, Vec<f64>)> = HashMap::new();
let mut gradient_x_map: HashMap<(u16, u16), (Ansi, Vec<f64>)> = HashMap::new();
let mut gradient_y_map: HashMap<(u16, u16), (Ansi, Vec<f64>)> = HashMap::new();

// Define the seeding Ansi order 1 and 2 zernike polynomials
let mut seed_zernike_1 = vec![1.0; coordinates.len()];
let mut seed_zernike_2 = vec![1.0; coordinates.len()];

for pos in 0..coordinates.len() {
seed_zernike_1[pos] *= coordinates[pos].y;
seed_zernike_2[pos] *= coordinates[pos].x;
}

// Define a variable for the ansi order going forward
let mut ansi_order = 2;

// Add in the seed values to our hashmaps
zernike_map.insert((0, 0), (Ansi { j: 0 }, vec![1.0; coordinates.len()]));
zernike_map.insert((1, 0), (Ansi { j: 1 }, seed_zernike_1));
zernike_map.insert((1, 1), (Ansi { j: 2 }, seed_zernike_2));

gradient_x_map.insert((0, 0), (Ansi { j: 0 }, vec![0.0; coordinates.len()]));
gradient_x_map.insert((1, 0), (Ansi { j: 1 }, vec![0.0; coordinates.len()]));
gradient_x_map.insert((1, 1), (Ansi { j: 2 }, vec![1.0; coordinates.len()]));

gradient_y_map.insert((0, 0), (Ansi { j: 0 }, vec![0.0; coordinates.len()]));
gradient_y_map.insert((1, 0), (Ansi { j: 1 }, vec![1.0; coordinates.len()]));
gradient_y_map.insert((1, 1), (Ansi { j: 2 }, vec![0.0; coordinates.len()]));

// Outer loop for radial order index
for n in 2..max_order + 1 {
// Inner loop for azimuthal index
for m in 0..n + 1 {
// Step the ansi order for our next iterations
ansi_order += 1;
// Define the temporary vectors to hold the new calculated zernike and gradient values
// They get consumed in each m loop into a hashmap, hence instantiating them here
let mut new_zern_values = Vec::with_capacity(coordinates.len());

// Handle the various cases of exception
if m == 0 {
// Collect Needed Prior Zernike and Gradient Polynomials
// Index (n - 1, 0)
let (_, zern_n1_0) = zernike_map.get(&(n - 1, 0)).unwrap();
// Index (n - 1, n - 1)
let (_, zern_n1_n1) = zernike_map.get(&(n - 1, n - 1)).unwrap();

// Iterating through every input coordinate point, define the new
// Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in 0..coordinates.len() {
let x: f64 = coordinates[pos].x;
let y: f64 = coordinates[pos].y;

new_zern_values.push(x * zern_n1_0[pos] + y * zern_n1_n1[pos]);
}
} else if m == n {
// Collect Needed Prior Zernike and Gradient Polynomials
// Index (n - 1, 0)
let (_, zern_n1_0) = zernike_map.get(&(n - 1, 0)).unwrap();
// Index (n - 1, n - 1)
let (_, zern_n1_n1) = zernike_map.get(&(n - 1, n - 1)).unwrap();

// Iterating through every input coordinate point, define the new
// Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in 0..coordinates.len() {
let x: f64 = coordinates[pos].x;
let y: f64 = coordinates[pos].y;

new_zern_values.push(x * zern_n1_n1[pos] - y * zern_n1_0[pos]);
new_gradient_y.push(-1.0 * n as f64 * zern_n1_0[pos]);
}
} else if n % 2 != 0 && m == (n - 1) / 2 {
// Collect Needed Prior Zernike and Gradient Polynomials
// Index (n - 1, n - 1 - m)
let (_, zern_n1_n1m) = zernike_map.get(&(n - 1, n - 1 - m)).unwrap();

// Index (n - 1, m - 1)
let (_, zern_n1_m1) = zernike_map.get(&(n - 1, m - 1)).unwrap();

// Index (n - 1, n - m)
let (_, zern_n1_nm) = zernike_map.get(&(n - 1, n - m)).unwrap();

// Index (n - 2, m - 1)
let (_, zern_n2_m1) = zernike_map.get(&(n - 2, m - 1)).unwrap();

// Index (n - 2, m - 1)

// Index (n - 2, m - 1)

// Iterating through every input coordinate point, define the new
// Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in 0..coordinates.len() {
let x: f64 = coordinates[pos].x;
let y: f64 = coordinates[pos].y;

// Calculate new values
new_zern_values.push(
y * zern_n1_n1m[pos] + x * zern_n1_m1[pos]
- y * zern_n1_nm[pos]
- zern_n2_m1[pos],
);

n as f64 * zern_n1_n1m[pos] - n as f64 * zern_n1_nm[pos]
);
}
} else if n % 2 != 0 && m == (n - 1) / 2 + 1 {
// Collect Needed Prior Zernike and Gradient Polynomials

// Index (n - 1, n - 1 - m)
let (_, zern_n1_n1m) = zernike_map.get(&(n - 1, n - 1 - m)).unwrap();

// Index (n - 1, m)
let (_, zern_n1_m) = zernike_map.get(&(n - 1, m)).unwrap();

// Index (n - 1, m - 1)
let (_, zern_n1_m1) = zernike_map.get(&(n - 1, m - 1)).unwrap();

// Index (n - 2, m - 1)
let (_, zern_n2_m1) = zernike_map.get(&(n - 2, m - 1)).unwrap();

// Index (n - 2, m - 1)

// Index (n - 2, m - 1)

// Iterating through every input coordinate point, define the new
// Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in 0..coordinates.len() {
let x: f64 = coordinates[pos].x;
let y: f64 = coordinates[pos].y;

// Calculate new values
new_zern_values.push(
x * zern_n1_m[pos] + y * zern_n1_n1m[pos] + x * zern_n1_m1[pos]
- zern_n2_m1[pos],
);

n as f64 * zern_n1_m[pos]
+ n as f64 * zern_n1_m1[pos]
);
}
} else if n % 2 == 0 && m == n / 2 {
// Collect Needed Prior Zernike and Gradient Polynomials

// Index (n - 1, n - 1 - m)
let (_, zern_n1_n1m) = zernike_map.get(&(n - 1, n - 1 - m)).unwrap();

// Index (n - 1, m)
let (_, zern_n1_m) = zernike_map.get(&(n - 1, m)).unwrap();

// Index (n - 1, m - 1)
let (_, zern_n1_m1) = zernike_map.get(&(n - 1, m - 1)).unwrap();

// Index (n - 2, m - 1)
let (_, zern_n2_m1) = zernike_map.get(&(n - 2, m - 1)).unwrap();

// Index (n - 2, m - 1)

// Index (n - 2, m - 1)

// Iterating through every input coordinate point, define the new
// Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in 0..coordinates.len() {
let x: f64 = coordinates[pos].x;
let y: f64 = coordinates[pos].y;

// Calculate new values
new_zern_values.push(
2.0 * x * zern_n1_m[pos] + 2.0 * y * zern_n1_m1[pos] - zern_n2_m1[pos],
);

}
} else {
// Collect Needed Prior Zernike and Gradient Polynomials

// Index (n - 1, n - m)
let (_, zern_n1_nm) = zernike_map.get(&(n - 1, n - m)).unwrap();

// Index (n - 1, n - 1 - m)
let (_, zern_n1_n1m) = zernike_map.get(&(n - 1, n - 1 - m)).unwrap();

// Index (n - 1, m)
let (_, zern_n1_m) = zernike_map.get(&(n - 1, m)).unwrap();

// Index (n - 1, m - 1)
let (_, zern_n1_m1) = zernike_map.get(&(n - 1, m - 1)).unwrap();

// Index (n - 2, m - 1)
let (_, zern_n2_m1) = zernike_map.get(&(n - 2, m - 1)).unwrap();

// Index (n - 2, m - 1)

// Index (n - 2, m - 1)

// Iterating through every input coordinate point, define the new
// Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in 0..coordinates.len() {
let x: f64 = coordinates[pos].x;
let y: f64 = coordinates[pos].y;

// Calculate new values
new_zern_values.push(
x * zern_n1_m[pos] + y * zern_n1_n1m[pos] + x * zern_n1_m1[pos]
- y * zern_n1_nm[pos]
- zern_n2_m1[pos],
);

n as f64 * zern_n1_m[pos]
+ n as f64 * zern_n1_m1[pos]
);
n as f64 * zern_n1_n1m[pos] - n as f64 * zern_n1_nm[pos]
);
}
}
// Add the new values to our nifty hashmap

zernike_map.insert((n, m), (Ansi { j: ansi_order }, new_zern_values));
} // End of inner azimuthal loop
} // End of outer radial order loop

// There is no real reason to convert out of a hashmap here, other than we use a vec outside of this.
// However, I am pretty fond of the hashmap and may change it to output that instead
let output_zern = zernike_map
.values()
.cloned()
.collect::<Vec<(Ansi, Vec<f64>)>>();

.values()
.cloned()
.collect::<Vec<(Ansi, Vec<f64>)>>();

.values()
.cloned()
.collect::<Vec<(Ansi, Vec<f64>)>>();

}


Python Code

from collections import defaultdict
import numpy as np

"""Recursively generate Zernike and Zernike Gradient polynomials and evaluate at input x and y coordinates
Algorithm and background from this paper:
Torben B. Andersen, "Efficient and robust recurrence relations for the Zernike circle polynomials and their derivatives in Cartesian coordinates,"
Opt. Express 26, 18878-18896 (2018)

// Author : Logan Rodriguez Graves. Date: 1/27/022

Parameters
----------
x : numpy.ndarray
x coordinates to evaluate at
y : numpy.ndarray
y coordinates to evaluate at
max_order : int
polynomial order

Returns
-------
Note: n, m order is the radial and azimuthal order respectively, with m = 0:n
dict
zernike polynomial evaluated at the given points, with key values of the n,m order

dict
zernike x gradient  polynomial evaluated at the given points, with key values of the n, m order

dict
zernike y gradient polynomial evaluated at the given points, with key values of the n, m order
"""

# Determine the max radial order, if it is less than 2 then set it to at least 2.
if max_order < 2:
max_order = 2

num_coords = len(x)

# Dict is your friend here! They will allow us to directly reference into
# the sets we need, as opposed to an unnecessarily extra complex step of translating nm
# indices to a single order.

# declare a default dict
zernike_map = defaultdict(list)

# Assign seed values
zernike_map[(0, 0)].append((0, [1 for i in range(len(x))]))
zernike_map[(1, 0)].append((1, y))
zernike_map[(1, 1)].append((2, x))

gradient_x_map[(0, 0)].append((0, [0 for i in range(len(x))]))
gradient_x_map[(1, 0)].append((1, [0 for i in range(len(x))]))
gradient_x_map[(1, 1)].append((2, [1 for i in range(len(x))]))

gradient_y_map[(0, 0)].append((0, [0 for i in range(len(x))]))
gradient_y_map[(1, 0)].append((1, [1 for i in range(len(x))]))
gradient_y_map[(1, 1)].append((2, [0 for i in range(len(x))]))

# Define a variable for the ansi order going forward
ansi_order = 2

# Outer loop for radial order index
for n in range(2, max_order + 1):
# Inner loop for azimuthal index
for m in range(0, n + 1):
# Step the ansi order for our next iterations
ansi_order += 1

new_zern_values = []

# Handle the various cases of exception
if m == 0:
# Collect Needed Prior Zernike and Gradient Polynomials
# Index(n - 1, 0)
zern_n1_0 = zernike_map.get((n - 1, 0))[0][1]
# Index(n - 1, n - 1)
zern_n1_n1 = zernike_map.get((n - 1, n - 1))[0][1]

# Iterating through every input coordinate point, define the new
# Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in range(0, num_coords):

new_zern_values.append(
x[pos] * zern_n1_0[pos] + y[pos] * zern_n1_n1[pos])

elif m == n:
# Collect Needed Prior Zernike and Gradient Polynomials
# Index(n - 1, 0)
zern_n1_0 = zernike_map.get((n - 1, 0))[0][1]
# Index(n - 1, n - 1)
zern_n1_n1 = zernike_map.get((n - 1, n - 1))[0][1]

# Iterating through every input coordinate point, define the new
# Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in range(0, num_coords):

new_zern_values.append(
x[pos] * zern_n1_n1[pos] - y[pos] * zern_n1_0[pos])

elif n % 2 != 0 and m == (n - 1) / 2:
# Collect Needed Prior Zernike and Gradient Polynomials
# Index(n - 1, n - 1 - m)
zern_n1_n1m = zernike_map.get((n - 1, n - 1 - m))[0][1]

# Index(n - 1, m - 1)
zern_n1_m1 = zernike_map.get((n - 1, m - 1))[0][1]

# Index(n - 1, n - m)
zern_n1_nm = zernike_map.get((n - 1, n - m))[0][1]

# Index(n - 2, m - 1)
zern_n2_m1 = zernike_map.get((n - 2, m - 1))[0][1]

# Index(n - 2, m - 1)

# Index(n - 2, m - 1)

# Iterating through every input coordinate point, define the new
# Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in range(0,  num_coords):

# Calculate new values
new_zern_values.append(y[pos] * zern_n1_n1m[pos] + x[pos] *
zern_n1_m1[pos] - y[pos] * zern_n1_nm[pos] - zern_n2_m1[pos])

n * zern_n1_n1m[pos] - n * zern_n1_nm[pos] + gradient_y_n2_m1[pos])

elif n % 2 != 0 and m == (n - 1) / 2 + 1:
# Collect Needed Prior Zernike and Gradient Polynomials

# Index(n - 1, n - 1 - m)
zern_n1_n1m = zernike_map.get((n - 1, n - 1 - m))[0][1]

# Index(n - 1, m)
zern_n1_m = zernike_map.get((n - 1, m))[0][1]

# Index(n - 1, m - 1)
zern_n1_m1 = zernike_map.get((n - 1, m - 1))[0][1]

# Index(n - 2, m - 1)
zern_n2_m1 = zernike_map.get((n - 2, m - 1))[0][1]

# Index(n - 2, m - 1)

# Index(n - 2, m - 1)

# Iterating through every input coordinate point, define the new
# Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in range(0, num_coords):

# Calculate new values
new_zern_values.append(x[pos] * zern_n1_m[pos] + y[pos] *
zern_n1_n1m[pos] + x[pos] * zern_n1_m1[pos] - zern_n2_m1[pos])

n * zern_n1_m[pos] + n * zern_n1_m1[pos] + gradient_x_n2_m1[pos])

elif n % 2 == 0 and m == n / 2:
# Collect Needed Prior Zernike and Gradient Polynomials

# Index(n - 1, n - 1 - m)
zern_n1_n1m = zernike_map.get((n - 1, n - 1 - m))[0][1]

# Index(n - 1, m)
zern_n1_m = zernike_map.get((n - 1, m))[0][1]

# Index(n - 1, m - 1)
zern_n1_m1 = zernike_map.get((n - 1, m - 1))[0][1]

# Index(n - 2, m - 1)
zern_n2_m1 = zernike_map.get((n - 2, m - 1))[0][1]

# Index(n - 2, m - 1)

# Index(n - 2, m - 1)

# Iterating through every input coordinate point, define the new
# Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in range(0, num_coords):

# Calculate new values
new_zern_values.append(
2.0 * x[pos] * zern_n1_m[pos] + 2.0 * y[pos] * zern_n1_m1[pos] - zern_n2_m1[pos])

2.0 * n * zern_n1_m[pos] + gradient_x_n2_m1[pos])
2.0 * n * zern_n1_n1m[pos] + gradient_y_n2_m1[pos])

else:
# Collect Needed Prior Zernike and Gradient Polynomials

# Index(n - 1, n - m)
zern_n1_nm = zernike_map.get((n - 1, n - m))[0][1]

# Index(n - 1, n - 1 - m)
zern_n1_n1m = zernike_map.get((n - 1, n - 1 - m))[0][1]

# Index(n - 1, m)
zern_n1_m = zernike_map.get((n - 1, m))[0][1]

# Index(n - 1, m - 1)
zern_n1_m1 = zernike_map.get((n - 1, m - 1))[0][1]

# Index(n - 2, m - 1)
zern_n2_m1 = zernike_map.get((n - 2, m - 1))[0][1]

# Index(n - 2, m - 1)

# Index(n - 2, m - 1)

# Iterating through every input coordinate point, define the new
# Zernike, and Zernike Gradient values evaluated at the coordinates input.
for pos in range(0, num_coords):

# Calculate new values
new_zern_values.append(x[pos] * zern_n1_m[pos] + y[pos] * zern_n1_n1m[pos] +
x[pos] * zern_n1_m1[pos] - y[pos] * zern_n1_nm[pos] - zern_n2_m1[pos])

n * zern_n1_m[pos] + n * zern_n1_m1[pos] + gradient_x_n2_m1[pos])
n * zern_n1_n1m[pos] - n * zern_n1_nm[pos] + gradient_y_n2_m1[pos])

# Add the new values to our nifty dictionaries

zernike_map[(n, m)].append((ansi_order, new_zern_values))
#  End of inner azimuthal loop
# End of outer radial order loop

# Define our input coordinates
x = np.array([0, 0, 0, 0, 0])
y = np.array([1, 0.5, 0, -0.5, -1])

# Define our max radial order to calculate to
max_order = 4

# Get our values out

1 Like

FYI, prysm does radial derivatives of Zernike polynomials on the dev branch (and whatever the ‘natural’ derivatives are of all the other polynomials it can work with…). Under experimental/raytracing/surface there are functions which convert cartesian and polar derivatives to each other, too. And Q2D derivatives for even the off-axis conic case, which was brutal and miserable to implement.

I did see Andersen’s paper at the time (Greg Forbes directed me to it), but I also found the pseudo-code too much of a mess to use. And I ended up not wanting to have a whole second set of Zernike code for derivatives, when the cost of converting polar to Cartesian is not all that high. You can still trace a billion rays per second with pure python with those conversions happening, so

A test you can do is to compare the closed-form approach to finite differences of the polynomials themselves – it is straightforward to find cross-comparison test cases for the polyomials. As N\rightarrow \infty the difference between finite differences and closed form should vanish. And, if you have an error in closed form it is almost always extremely obvious.

In your code, you have tuples of integers for keys, which is hard to beat for performance (python’s hash function is optimized for that case). It takes about 50 nanoseconds to do the dict lookups or assignments. Even with 64x64, likely the smallest useful array size, the hashing gets amortized to 12 picoseconds per element, which is a rounding error on a rounding error. The list comp initialization in the dict are not as nice, though – np.zeros_like(x) and np.ones_like(x) would be a lot better (likely more clear, and a ton faster). You’ll also see those operations become literally free if you run your function multiple times within the same python garbage collection generation, for finessed reasons having to do with numpy’s internal memory pool.

The same “write it vectorized instead of in just python” applies there - -the rest of the python reads like it was converted from the Rust code, but shouldn’t it be the other way around?

In the Rust code (I mostly cannot read rust, please bear with me):

• don’t bother with u16 for the ANSI order. On 32 and 64 bit machines the cost of doing things with the native integer type is the same as smaller integer types. While 65k Zernike orders is a lot, it doesn’t really do anything to the code and most systems programmers will see the u16 as important to something, when it isn’t really important to this algorithm.

There is the potential for something interesting in the hash map, since you could use u32s for the indices and bit concatenate them to a u64 for the keys. The hash of integers are themselves and any decent hash function is a no-op on integer types, so you could make your lookups literally free there.

Densely packed memory is important to making numerical algorithms go zoom zoom. In python (well, numpy) that’s always true and out of your control. I don’t know if your Rust code has those guarantees (I presume vec satisfies it, but the entries in the hashmap may be scattered in memory). You should be able to get the recurrence relation down to 21 CPU clocks per element of the array for all orders 0…5 inclusive. If it’s much slower than that, I would presume the problem is the data is too far apart in memory and causes the CPU to spend all of its time waiting to get data out of RAM (which takes about 100-200 cycles…)

I believe it. Thanks for the heads up, definitely more of a fan of the gradients as you derived them in Prysm so I may end up going that route as well; like you said 2 methods for doing the zernikes feels less than ideal.

Ha, you got me there. I feel a lot more comfortable writing in Rust than Python, so I converted from the Rust code to Python.

Thats a great idea, I will look into this as it could lead to a more efficient code.