Incoherent sum of PSF from configurations

  • 21 September 2023
  • 8 replies

I have a system where the PSF is constructed from several different configurations (>2) simultaneously. I see that the Huygens PSF and MTF can do a coherent sum from different configurations (which is awesome). However in my case I expect them to be mostly incoherent to one another. Is there a way to directly generate incoherent sum PSF and MTF results?


Best answer by Hui Chen 3 October 2023, 08:14

View original

8 replies

Userlevel 5
Badge +3

Do each Huygens PSF separately, each at a different single wavelength or coherent bandwidth. Then add their intensities, weighted by the wavelength-dependent system transmission.  Image delta must also be scaled by wavelength to hold the overall transform array size constant.  Normalize the sum to unity at the peak. That gives the "white" light incoherent intensity PSF.

Userlevel 7
Badge +3


Unfortunately, I think the answer is no -- there is no simple way to have the Huygens PSFs add incoherently across configurations via the Huygens analysis tool inside the GUI.  I’ve tried multiple ways without any success.  

I know you indicated that you are using > 2 configurations, but if you did have just two, then you could use polarization ray tracing and put the two configs in quadrature so that the coherent sum mimics the incoherent sum as shown here:


Otherwise, I think you have to write a script, say using the ZOS-API, in order to extract the PSF for each configuration separately and then add them on an intensity basis. 

It would be nice if there were a simple switch to toggle between coherent and incoherent sums, but that would be a new feature request. 



Thanks Jeff & Mike! It would indeed be nice if this was directly available, maybe in a future version.  I the meantime the script route looks reasonable. Thanks again!

Userlevel 6
Badge +2

Hey @kiteflyer ,

One other thing to add onto Mike J’s answer is you also have to convert the PSF grid from “image pixel” space to “image micron” space.  Basically, for the Huygens PSF calculation, OpticStudio will give you a centroid offset from the Image Vertex but then all the other values for the PSF are in terms of deltas (PSF Spacing).  So, if you have 2 PSFs at different wavelengths with offset centroids, you’ll need to transfer the PSF’s into a common reference w.r.t. the image vertex:

This can easily be done with a bilinear interpolation based on the centroid offset (if the sampling is high enough, bilinear is fine as a bicubic can overshoot the interpolated PSF points and give aliased bilinear will be faster).  Matlab has a built-in function called griddedInterpolant that can help with this, but since I try to avoid Matlab for several reasons (speed, GUI interfaces, licenses, multi-threading) & a bilinear algorithm is really easy, the steps below show how to do this:

The psuedo C# code to do this is:

public double[,] IncoherentPSF(double[,,] psfstack, double[] centerx, double[] centery, double[] weight, double dataspacing)
int ni = psfstack.GetLength(0); // number of rows
int nj = psfstack.GetLength(1); // number of columns
int nw = psfstack.GetLength(2); // number of wavelengths (or configs)

double[,] psf = new double[psfstack.GetLength(0), psfstack.GetLength(1)];

// loop through all wavelengths and interpolate the pixels
for (int w = 0; w < nw; w++)
for (int i = 0; i < ni; i++)
for (int j = 0; j < nj; j++)
psf[i, j] += Bilinear(psfstack, i, j, w, centerx[w], centery[w], weight[w], dataspacing);

return psf;

private double Bilinear(double[,,] psf, int i, int j, int w, double centerx, double centery, double weight, double dataspacing)
// convert from centroid offset in microns to centroid offset in "pixels"
double pixelx = centerx / dataspacing;
double pixely = centery / dataspacing;

// shifted i-th and j-th
int si = (int)Math.Floor(i + pixelx);
int sj = (int)Math.Floor(j + pixely);

// make sure the centroid offset is still inside the PSF grid, otherwise return original PSF
// there should be enough 0-padding that the PSF is not at the edge of the grid
if ((si < 0 && i == 0) || (si > 0 && i == ni - 1))
return res;
if ((sj < 0 && j == 0) || (sj > 0 && j == nj - 1))
return res;

// get the next pixel index (we always go 1 down and 1 right)
int ii = si + 1;
int jj = sj + 1;

// bilinear interpolation is a 3 step linear interpolation process:
// 1. get the horizontal linear interpolation between [i, j] & [ii, j]
// 2. get the horizontal linear interpolation between [i, jj] & [ii, jj]
// 3. get the vertical linear interpolation on this new line

double xratio = (float)(pixelx - Math.Floor(pixelx));
double yratio = (float)(pixely - Math.Floor(pixely));

double t = psf[i, j, w] + (psf[ii, j, w] - psf[i, j, w]) * xratio;
double b = psf[i, jj, w] + (psf[ii, jj, w] - psf[i, jj, w]) * xratio;

return (t + (b - t) * yratio) * weight;

The psfstack input is simply the PSFs to combine incoherently, the centerx/centery/dataspacing are determined from the text output of the Huygens PSF, and the weight will be the respective wavelength weighting to get the PSF intensity correct.

Also, just to emphasize, this can only work with Huygens PSF and not FFT PSF, since the FFT is calculated in pupil space while the Huygens is calculated in image space. 

Userlevel 4
Badge +1

Hi @kiteflyer 
The GUI Huygens PSF with "Configuration: All" does a coherent sum across all configurations. There isn't a direct or easy way to do incoherent sum of all PSFs in the GUI. As Jeff and Michael pointed out, you'll have to do this using a customized API script. 
My colleague and I are developing a Matlab script to do incoherent PSF sum. We are doing this as part of an effort to develop a jitter analysis workflow for an imaging system subject to random vibrations. Due to fast jitter impact, lens position may oscillate within one sensor integration time, leading to a bigger blur on the image sensor and therefore lower MTF and lower resolution. Summing PSFs from lens vibration helps us analyze the jitter impact on MTF and overall imaging performance. 
The tricky part in PSF sum lies in two aspects. First by default PSF grid is centered on chief ray. Since chief ray intercept likely varies across configurations, you need to first align these PSF grids or recenter them before you can sum their intensity point by point. Second, if the point spacing is not the same across all configurations, you will need to resample the grid. In our API workflow, we first shift these PSF grids based on the chief ray intercept in each configuration, then we resample these PSFs using a bigger grid, then we sum them point by point. 

I agree that implementing an Incoherent PSF Sum across configurations may be a useful feature to add to the GUI. We have proposed this to our Dev team for evaluation. We'll keep you posted if there is any update on that front. Let me know if you have any questions. 


Userlevel 7
Badge +3

Hi @Hui Chen,

Can you please provide a copy of the Matlab script you describe above for incoherent PSF sums?  Also, can you provide an update from the Dev team on the possibility of adding this feature as a built-in option?



Userlevel 4
Badge +1

Hi @Jeff.Wilde 

Thank you for your inquiry! Also, it was nice to meet you in person at the Photonics West!

I have not heard anything from Dev regarding this feature. I’ll check if there is a feature request already in place. If not, I’ll create one to log this formally in our customer Feature Request system and put your name down for it. 

As for the script, one of my colleagues has been working on it so she’ll have more info. Could you please create a support ticket for the request, and I’ll tag her to help you?

Thank you and have a great weekend!


Userlevel 7
Badge +3

Hi @Hui Chen ,

Thanks for your response.  Yes, it was good to finally meet face-to-face at Photonics West!

I just submitted a support ticket, so hopefully a copy of the Matlab script can be made available.  I appreciate your help.