Question about Huygens PSF

[hemkumarsrinivas]

Hi,
I was recently trying to write my own function for the Huygens PSF. I have tried a very simplistic approach so far, where I have a spot diagram and assume each spot in the diagram to be a spherical emitter. A coherent sum over the emitters then gives the PSF.

I have tried the following implementation, with a dummy 'spot diagram', where I just generate a grid of uniformly distributed points.

# function to calculate Huygens contribution
def huygens_contribution(x_spot, y_spot, x_obs, y_obs):
    r = np.sqrt((x_obs - x_spot)**2 + (y_obs - y_spot)**2)
    amplitude = np.exp(1j * (2 * np.pi / wavelength) * r)  # Example phase term
    return amplitude

x_spots = [...]  # List of x coordinates of the spots
y_spots = [...]  # List of y coordinates of the spots
X_obs, Y_obs = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))  # Observation grid
wavelength = 500e-9  # Example wavelength in meters

# Initialize a complex array for coherent summation
amplitude_sum = np.zeros(X_obs.shape, dtype=np.complex128)

# Loop over the spots
for x_spot, y_spot in zip(x_spots, y_spots):
    # Calculate contribution  for each observation point
    for i in range(X_obs.shape[0]):
        for j in range(X_obs.shape[1]):
            amplitude_sum[i, j] += huygens_contribution(x_spot, y_spot, X_obs[i, j], Y_obs[i, j])

# Calculate the intensity from  coherent sum
psf = np.abs(amplitude_sum)**2 

At the end of it, I am able to generate a Huygens PSF. Here's an example for it:

One main challenge I face is its speed. I wanted to ask, if there are any ideas, suggestions to improve this ? Moreover, how can I extract the coordinates of the spots from the spot diagram to use with such an implementation ?

It would be great if you could share some ideas on this. I will be happy to work further on this and contribute to the package.

Edited by Kramer on Feb. 1st, 2025 to format code as Python

[HarrisonKramer]

This is great, thanks for sharing! The Huygens PSF has been on the roadmap for a long time. It's really about time to get it implemented.

Speed-up possibilities

Regarding speed, you can use NumPy broadcasting to avoid for-loops entirely. That turns your code into, for example:

# use broadcasting here
X_spots = x_spots[:, np.newaxis, np.newaxis]  # shape (N, 1, 1)
Y_spots = y_spots[:, np.newaxis, np.newaxis]  # shape (N, 1, 1)
R = np.sqrt((X_obs - X_spots) ** 2 + (Y_obs - Y_spots) ** 2)  # shape (N, M, M)

# phase term
amplitude = np.exp(1j * k * R)

# coherent sum
amplitude_sum = np.sum(amplitude, axis=0)

# compute PSF
psf = np.abs(amplitude_sum) ** 2

Alternatively, you can use numba for big speed-ups and parallelization. I do this in the scatter module.

Or, the operations could be moved to GPU using torch or cupy, among others. I started writing a configurable backend for Optiland, but it's not yet complete.

On the implementation

You can trace rays through the system over a uniform pupil grid as follows. This traces 64 rays per axis, so really it's closer to ≈64^2 rays in total. Many other distributions are possible, see the distribution module or API docs.

lens = CookeTriplet()
rays = lens.trace(Hx=0, Hy=1, wavelength=0.5876, num_rays=64, distribution='uniform')

The returned rays instance has the information you need:

x_points = rays.x
y_points = rays.y
z_points = rays.z
dir_cos_x = rays.L  # x component of the direction cosines
dir_cos_y = rays.M  # y component of the direction cosines
dir_cos_z = rays.N  # z component of the direction cosines
opd = rays.opd  # optical path difference

It would be great if you could get a prototype working. Eventually, the completed implementation can go into its own class in the psf module.

Happy to discuss more, help where needed, or provide feedback as you investigate. I'll keep this issue open to monitor the progress of the implementation.

Thanks,
Kramer

[hemkumarsrinivas]

Hi,
Thanks a lot for the response. I think I will stick to broadcasting arrays on numpy for the time being. I will then try making a prototype using the information from the rays instance. Perhaps, when things work well, I will try implementing it with namba.

[hemkumarsrinivas]

Hi,
So, I managed to work further based on the previous discussion. I have to tried to generate the Huygens PSF for the CookeTriplet example.

from optiland.samples.objectives import CookeTriplet
import numpy as np
from optiland import analysis
import matplotlib.pyplot as plt
from scipy.fftpack import fftshift, fft2
lens = CookeTriplet()

def huygens_psf(lens, field, wavelength):
    rays = lens.trace(Hx=field[0], Hy=field[1], wavelength=wavelength, num_rays=64, distribution='uniform')
    x_spots = rays.x  # List of x coordinates of the spots
    y_spots = rays.y  # List of y coordinates of the spots
    X_obs, Y_obs = np.meshgrid(np.linspace(min(x_spots),max(x_spots),100), np.linspace(min(y_spots),max(y_spots),100))  # Observation grid
  
    X_spots = x_spots[:, np.newaxis, np.newaxis]  # shape (N, 1, 1)
    Y_spots = y_spots[:, np.newaxis, np.newaxis]  # shape (N, 1, 1)
    R = np.sqrt((X_obs - X_spots) ** 2 + (Y_obs - Y_spots) ** 2)  # shape (N, M, M)
    k = 2*np.pi/wavelength
    # phase term
    amplitude = np.exp(1j * k * R)
    # coherent sum
    amplitude_sum = np.sum(amplitude, axis=0)
    # compute PSF
    huygens_psf = np.abs(amplitude_sum) ** 2 
    ax = plt.subplot(111)
    h = ax.pcolormesh(X_obs,Y_obs,huygens_psf/np.max(huygens_psf))
    plt.colorbar(h,label='Normalized Intensity')
    ax.set_xlabel('X (mm)')
    ax.set_ylabel('Y (mm)')
    ax.set_title('Hx='+str(field[0])+', Hy='+str(field[1]))

huygens_psf(lens,(0,0),0.550)
huygens_psf(lens,(0,0.7),0.550)
huygens_psf(lens,(0,1),0.550)

The results are as follows

When compared to the spot diagrams, what I see is, there is not much intensity in the PSF's y-direction for the off-axis field points, inspite of having a few spots in the y-direction due to the coma. Any comments on this ?

Also, I would like to generate the MTF plots based on the Huygens PSF. It would be great if you could share some ideas or point me to some literature with relevant details. I can eventually define all of this in a class as you have suggested and copy the structure that you have for the FFT and Geometric PSF classes.

Regards,
Hem

[HarrisonKramer]

Hi Hem,

Thanks for sharing!

A few remarks:

• I notice that the plots have a very narrow range (≈0.99-1.0). I wonder what they look like over a larger grid.
• You'll also need a phase term that takes into account the optical path length for each ray - something like $\exp{(i \cdot k \cdot OPL)}$.
• I now also realize you'll have to apply a field-dependent correction when tracing off-axis rays. I needed to do this for the Wavefront class. Essentially, this is to assure that if rays are traced from the object point to the reference sphere in object space, that they have the same starting phase. I think this is adding too much complexity for now. We can come back to this when needed.

For now, I would recommend looking at the on-axis field point only to keep it simpler. You should be able to reproduce an Airy disk pattern when the rays have identical OPL and merge to a point. Then, if using real rays from the Cooke triplet, you should be able to generate approximately the same PSF as is computed with the FFTPSF class. Lastly, you can check if it's the same when the image plane is defocused.

I think the MTF calculation will be relatively straightforward once the PSF is working, as they're related by a Fourier transform.

As for references, the classic book (and what I used in school) is Fourier Optics by Goodman. I list many more references in the documentation as well. And here's a useful blog post that summarizes the mathematics a bit.

Let me know if anything isn't clear or if you have any other questions. Happy to help where I can.

Regards,
Kramer

[hemkumarsrinivas]

Hi,
Thanks a lot for the inputs. I now mainly modified the 'amplitude' term as follows : . Firstly, I realized that I didn't factor the 1/R dependence in the calculation, which explains almost the constant intensity throughout the plots I posted earlier. Second, I tried to take the OPD itself as the phase contribution coming from the OPL. My understanding in this case is, we are effectively propagating backwards from the image plane to the exit pupil and the spots in the image plane are displaced from their ideal/reference positions by their OPDs.

So, far, it looks like definitely a step in the right direction, since the Airy disk diameter is more or less the same when compared to the FFT psf. The off-axis points also look similar. I tried generating the plots in my case with 128 rays and and observation grid of 128x128. The FFT psf also has a grid of 128x128 points.

I wonder, if the same number of rays were used in the FFT as well. If it were more than 128 rays, then would that explain lack of fine structures in the Huygens PSF ? I did try to generate the Huygens PSF using 256 rays, but that seems to take forever and clearly needs numba or another implementation. Let me know if the method I've used makes sense and how I could improve this to match the FFT PSF.

PS: The images are from the Cooke Triplet example once again.

[HarrisonKramer]

Hi,

Thanks for the update and good to see you're making progress here.

Thoughts on the Huygens PSF

I've had a bit more time to dig into this one. I agree we need to propagate from the image plane back to the exit pupil. This operation is already done in the Wavefront class, so you can likely use that implemention, or at least part of it. Ultimately, you need the pupil function:

$$ P(x,y) = A(x,y) , e^{i, \phi(x,y)}, $$

where $(x, y)$ is a point in the exit pupil, $A$ is the amplitude, and $\phi(x, y) = k \cdot OPD(x, y)$.

This pupil function is calculated in the class optiland.wavefront.Wavefront, but the $(x, y)$ points are not available by default, so you'd need to return them. Rays are traced through the system from object to image, then traced backwards to the exit pupil. Retrieving the OPD on the exit pupil can be done as follows for the Cooke triplet:

from optiland import wavefront

lens = CookeTriplet()
wavelength = 0.55
wf = wavefront.Wavefront(lens, fields=[(0, 1)], wavelengths=[wavelength], distribution='uniform', num_rays=32)

field_idx = 0  # only one field
wave_idx = 0  # only one wavelength
opd_waves, intensity = wf.data[field_idx][wave_idx]
opd_mm = opd_waves * (wavelength * 1e-3)  # wavelength is in microns, convert to mm

As you've largely already implemented, you can then compute the complex field at each point on the image plane as

$$ U(x_i, y_i) = \iint_{\text{pupil}} P(x,y) \frac{e^{ikr(x,y,z; x_i,y_i,z_i)}}{r(x,y,z; x_i,y_i,z_i)} dA. $$

where $r(x,y,z; x_i,y_i,z_i)$ is the distance from point $(x, y, z)$ in the exit pupil to image point $(x_i, y_i, z_i)$. Of course, you'll be computing the discrete equivalent. We also need to think about normalization, but that can wait for now.

On computational speed

• One trick you might consider for calculating $r$ is that if the wavefront is sufficiently flat on the image plane, you can assume a fixed r for all points. That could speed things up for large arrays. You'd have to check if it's actually flat first though.
• I question whether numba will be significantly faster than NumPy, although I'd be happy to be proven wrong. Ideally we can find a reasonable grid size in the pupil and image plane that don't result in an overly slow calculation.

On the FFT PSF implementation

• By default, this class uses ≈128x128 rays across the pupil and a 1024x1024 image grid, but this is resampled to better center the PSF in the image. This resampling method probably won't work for the Huygens PSF. You'll need some other approach - maybe similar to what you've already done.

I hope this helps you move things along. Please let me know if anything isn't clear, or if you think I missed anything. I still am not sure how best to implement this.

Kramer

[hemkumarsrinivas]

Hi Kramer,
after a hectic month in between and a brief hiatus from exploring Optiland, I am back to pursuing this topic once again. Thank you very much for sharing the information about calculating the pupil function using the wavefront class and thereby extending it to calculate the Huygens PSF.

On the PSF calculation

I did try implementing the method you had suggested and I get the following results :

The result for the on-axis field point seems to be the same as before, but the results for the off-axis field points are drastically different. I am not sure if this is because of a scaling factor involved or something else might be wrong with the implementation. Perhaps you have an idea why this could be ?

In addition, I was also thinking of an alternative method to calculate the PSF. As I understand, the spot diagram can be generated physically by introducing a tiny aperture and scanning it through the entire pupil (effectively a Hartmann test). Now, each point would in principle generates a Bessel function as a PSF. If we sum up the PSFs generated by each spot, we would then be generating the Huygens PSF.

Now, I tried this approach and obtained the following results :

I once again have ignored the computational speed aspect and have done the calculations with 128 rays and a grid of 128 observation points. For the on-axis field point, the Airy disc structure seems to emerge finally and begins looking similar to the FFT PSF. The off-axis points once again seem quite different and made me wonder if the high frequency modulation can be produced in the Huygens PSF at all, since the ringing structure is characteristic of FFTs and more the sampling points, the more prominent they are.

In any case, I found this to be interesting and wanted to ask if this method might make sense after all. Let me know your thoughts on this. Looking forward to your inputs.

Regards,
Hem

[HarrisonKramer]

Hi Hem,

No worries. No real rush. Thanks for continuing to investigate this.

I am still a bit puzzled why the direct integration method does not work as expected. I will also investigate this a bit. Can you share the code you used here? I will also work through the mathematics in more detail to make sure I'm not missing anything.

On your proposed approach - this is interesting, and it good to see that it recovers the expected PSF for the on-axis field. How specifically are you generating the PSF for each spot in the pupil? Perhaps you can share code here as well, or just quickly describe the approach.

I think we're getting closer! I will follow up here as soon as I have any updates. Thanks again.

Regards,
Kramer

[HarrisonKramer]

Hi Hem,

I worked through this formulation a bit. The equation I wrote above was missing the obliquity factor. I believe the full (discrete) equation should be

$E(x, y) = \frac{1}{i\lambda} \displaystyle\sum_j A_j \exp(ikL_j) \cdot \frac{\exp(ikR_j(x, y))}{R_j(x, y)} \cdot \frac{1 + \cos(\theta_j(x, y))}{2} \Delta S_j$,

where

• $k$ is the wavenumber.
• $\lambda$ is the wavelength.
• $E(x, y)$ is the complex field at point $(x, y)$ on the image surface.
• $A_j$ is the amplitude for point $j$ in the pupil.
• $L_j$ is the optical path difference.
• $R_j(x, y)$ is the distance from point $j$ in the pupil to point $(x, y)$ on the image surface.
• $\theta_j$ is the angle between the wavefront surface normal and the propagation direction (connecting pupil point and image point)
• $\Delta S_j$ is the pupil area element.

I put this into a python function and used parallelized numba for a speed-up. This runs in ≈1 second on my machine for a 128x128 pupil and image grid.

@njit(parallel=True, fastmath=True)
def compute_psf(image_x, image_y,
                pupil_x, pupil_y, pupil_z,
                pupil_amp, pupil_ops, dS,
                z_img, wavelength, Rp):
    """
    Compute the point spread function using the Huygens–Fresnel diffraction integral

    Args:
        image_x, image_y (np.ndarray): 2D arrays of image plane coordinates.
        pupil_x, pupil_y, pupil_z (np.ndarray): 1D arrays of pupil plane coordinates.
        pupil_amp (np.ndarray): 1D array of pupil plane amplitudes.
        pupil_opd (np.ndarray): 1D array of optical path difference.
        dS (np.ndarray): 1D array of pupil plane areas.
        z_img (float): Image plane distance from the pupil plane.
        wavelength (float): Wavelength of the light.
        Rp (float): Radius of the pupil plane.

    Returns:
        np.ndarray: 2D array of the point spread function.
    """
    k = 2.0 * np.pi / wavelength  # wavenumber
    Nx, Ny = image_x.shape
    field = np.zeros((Nx, Ny), dtype=np.complex128)

    # Constant factor
    const = 1.0 / (1j * wavelength)

    # Loop over all image points
    for ix in prange(Nx):
        for iy in range(Ny):
            x = image_x[ix, iy]
            y = image_y[ix, iy]
            sum_val = 0.0 + 0.0j  # initialize field sum for image point (x, y)

            # Loop over all pupil points
            for j in range(pupil_x.shape[0]):
                u = pupil_x[j]
                v = pupil_y[j]
                w = pupil_z[j]

                # Compute distance R from the pupil point to the image point
                dx = x - u
                dy = y - v
                dz = z_img - w
                R = np.sqrt(dx*dx + dy*dy + dz*dz)

                # Spherical propagation kernel:
                wave = np.exp(1j * k * R) / R

                # Compute the unit normal at the pupil point
                nux = u / Rp
                nuy = v / Rp
                nuz = w / Rp

                # Compute the cosine of the angle between (P - Q) and the pupil normal.
                # P is the image point, Q is the pupil point
                dot = dx * nux + dy * nuy + dz * nuz
                cos_theta = dot / R
                Q_obliq = 0.5 * (1.0 + cos_theta)  # obliquity factor

                # Pupil function
                pupil_phase = np.exp(1j * k * pupil_opd[j])

                # Add contribution to the field sum
                sum_val += pupil_amp[j] * pupil_phase * wave * Q_obliq * dS[j]

            field[ix, iy] = const * sum_val

    # Compute psf as the squared magnitude of the field
    psf = np.abs(field)**2
    return psf

I used the approach that I outlined above to retrieve the coordinates on the exit pupil, as well as the OPD. Here's a comparison of the FFT vs. Huygens approach for the Cooke triplet and the on-axis field:

and for the max field:

Seems to approximately match! I am not certain how best to verify that the implementation is correct. Also, I am still mostly ignoring the normalization, so that'll need to be worked out.

Here's the very rough workaround code I wrote to generate the required pupil info. This is overly complex and inefficient, so needs to be largely reworked.

def get_pupil_data(lens, field, wavelength=0.55, distribution='uniform', num_rays=128):
    """
    This is an inefficient workaround that largely replicates calculations in the wavefront.Wavefront class.
    This is needed to retrieve the pupil (x, y, z) coords for the PSF calculation.
    """
    wf = wavefront.Wavefront(lens, fields=[field], wavelengths=[wavelength], distribution=distribution, num_rays=num_rays)

    # distance in z from the last surface to the exit pupil
    pupil_z = (lens.paraxial.XPL() + lens.surface_group.positions[-1, 0])

    # trace the chief ray through the system
    wf._trace_chief_ray(field, wavelength)

    # reference sphere center and radius
    xc, yc, zc, R = wf._get_reference_sphere(pupil_z)

    # trace all rays
    rays = lens.trace(*field, wavelength, num_rays, distribution)

    field_idx = 0  # only one field
    wave_idx = 0  # only one wavelength
    opd_waves, intensity = wf.data[field_idx][wave_idx]
    opd_mm = opd_waves * (wavelength * 1e-3)  # wavelength is in microns, convert to mm

    t = wf._opd_image_to_xp(xc, yc, zc, R, wavelength)

    x = rays.x - t * rays.L
    y = rays.y - t * rays.M
    z = rays.z - t * rays.N

    return x, y, z, intensity, opd_mm, R.item()

field = (0, 0)
wavelength = 0.55
distribution = 'uniform'
num_rays = 128

image_x = np.linspace(-16e-3, 16e-3, 128)  # mm, image plane coordinates
image_x, image_y = np.meshgrid(image_x, image_x)

pupil_x, pupil_y, pupil_z, pupil_amp, pupil_opd, Rp = get_pupil_data(lens, field, wavelength, distribution, num_rays)

dS = np.full_like(pupil_x, (2*np.pi*Rp) / pupil_x.size)  # approximate area per pupil point

z_img = lens.surface_group.positions[-1, 0]

# compute Huygens-Fresnel PSF
psf = compute_psf(image_x, image_y, pupil_x, pupil_y, pupil_z, pupil_amp/pupil_amp.size, pupil_opd, dS, z_img, wavelength*1e-3, Rp)

I still wonder if numba is absolutely needed here. Perhaps it's still possible to do this efficiently using numpy only. What do you think of this implementation? Maybe it can be improved a bit or sped up in other ways. Eventually, we can build the final implementation into its own class e.g., HuygensPSF in psf.py.

edit: I also realized that I am assuming a planar, non-tilted image surface. This should also be made more general.

Interested to hear your thoughts.

Regards,
Kramer

[hemkumarsrinivas]

Hi Kramer,
Wow, the results look great and it seems that the obliquity factor is the key ingredient that I'd missed! Thanks a lot for sharing this and the explanation. It is very helpful in identifying where I missed out/what assumptions were wrong.

The PSF calculation

Firstly the role of the obliquity factor : I did read about the obliquity factor in Goodman's book, but was not sure how exactly I could get it. My naive assumption was that the theta could be got from the direction cosines for each point(using the rays class). Upon looking at your code, I see that you take the normals at the pupil points to calculate the theta.

Second, the area element dS in my case was also simply got from the grid dimensions. I didn't know about the possibility of getting the reference sphere radius and using its circumference to scale the area element. The rest of the calculation is more or less the same as I do and hence I'd say that I have "only" missed out the obliquity factor and the scaling of the area elements.

Regarding comparing/cross-checking the accuracy of the method: I was wondering if it could be done with a 3rd party software (like Zemax for eg.) with exactly the same parameters. From what I see in literature, there are barely any examples that we could implement and check and hence, I feel the only way might be to compare it with an existing platform.

Let me know if there might be a way to cross-check and if I could do something for this.

The workaround code

I think the wavefront class could have a function like pupil_data(), which returns all the parameters returned by get_pupil_data(). I feel that the scope of this function could be wider and can also be used in other cases. Perhaps doing this might also be able to reduce the number of variables that need to be defined adhoc for the purpose of using it in the compute_psf function. Only thing that might be a challenge while integrating them in the wavefront class would be the dependency on the rays class. Do you think this might be okay after all to do so ?

With regards to the implementation : I think for now, the numba implementation seems like a good method, since it is able to keep things fast. I see that with numpy alone, we might have to resort to broadcasting, which in my opinion might make the definition of variables a bit challenging. I could give this a shot and compare if an implementation with numpy alone might be helpul.

Next Steps

Now, having got a basic method to calculate the Huygens PSF, I was wondering how this could be taken further. I was thinking perhaps I'd start with a function to be included in the psf class. In addition, you'd also mentioned two other points to work on, namely the tilted image surfaces and normalization.

For the normalization, I was wondering if a cumulative distribution function might be enough. Is there any other method you think might be appropriate ?

For the tilted image surface, I would expect that the points on the image plane will have different reference spheres. Perhaps we should calculate the OPDs with the respective reference spheres for each z value in the image plane?

Let me know your thoughts on this. I will be happy to work more on this and contribute.

Regards,
Hem

[HarrisonKramer]

Hi Hem,

Good to hear this is nearly consistent with your implementation. I had not gone through these mathematics in quite some time and had forgotten much of this, but I'm glad the initial results look reasonable. It took a bit of trial and error. A few thoughts here:

On verifying the implementation

I agree it probably makes sense to check against other codes or commercial programs. I already have some confidence by comparing the FFT and Huygens approach above. The FFT approach is computed on a plane perpendicular to the chief ray, which I think explains why the Huygens PSF looks "stretched" in the y direction - this seems due to the projection onto the image plane, which is not perpendicular to the chief ray.

The validity of the calculation also relies on having sufficient pupil and image grid sampling. I suppose we can simply let the user decide what "sufficient" means and document this.

On the normalization

I suspect that the discrete nature of the integral will make normalization more difficult. The easiest approach might be to first compute the PSF for the system, ignoring any constant or normalization factors, then compute an on-axis "perfect" reference. This could be computed for a single on-axis point and it could use a perfect spherical wavefront with the same radius as derived above. This would define the 100% reference. As long as the calculation methods are otherwise consistent, the normalization should work in this way, I expect. I don't think this will add too much overheard either.

Implementation in Optiland

We start to move away from the physics and into some of the subtleties of the program. I'd like the implementation to follow a similar style to the FFT approach, e.g., you should be able to call it in a similar way, so that you can "swap out" the FFT method for the Huygens method if you wish. I like your proposal to update the Wavefront class to provide the necessary information for the PSF calculation. I see several things we'll probably need to do:

• Update the optiland.wavefront.Wavefront class to save the intersection points in the exit pupil and the reference sphere radius. In this class, all data is generated in the __init__ method. Ideally, this should also save the intersection points and reference sphere radius. In my implementation above, I am calling several methods twice. This is inefficient.
• Write a new class in psf.py called HuygensPSF. This will probably subclass Wavefront, as is done for FFTPSF. This should have similar 2D/3D plotting functionality to that implemented in FFTPSF.
• Add a static method with numba that computes the PSF (function above).
• Update the method for normalization, as discussed above. This could call the same PSF function, but the image grid would be a single on-axis point. If the method above is used, we can ignore any constant factors in the integral.
• Generalize the calculation for tilted or curved image planes. This should be straightforward, as I already have a function in optiland.visualization.utils called transform that takes (x, y, z) points and transforms them into the coordinate system of a given surface. Once we have this, we can pass it to the PSF function, replacing the z_img input with the image surface coordinates in 3D: (image_x, image_y, image_z).
• As is standard for any new functionality, add docs and tests.

These are just my initial thoughts. Please let me know if you see things differently. I am open to other approaches as well.

Feel free to work on this if you like. This could be a considerable amount of work and might require some further investigation. I am happy to work on this together, or provide guidance. Let me know your thoughts.

Regards,
Kramer

[hemkumarsrinivas]

Hi Kramer,
The sequence that you've proposed for implementation seems good. I am happy to start with updating the wavefront class and then moving on to the HuygensPSF class. I will work on this in the coming days and reach out to you if there are any questions.

Regards,
Hem

[HarrisonKramer]

Hi Hem,

I hope you're doing well! I just wanted to check in and see if you’ve had a chance to look into the updates for the Wavefront and HuygensPSF classes. No rush - just wanted to touch base and see how it's going or if there’s anything I can do to help.

Best,
Kramer