add SurfaceViewer to vizualize surfaces

[drpaprika]

Following the implementation of Zernike surfaces I was curious about their shape. In freeform manufacturing, deviation to the best sphere is crutial.

Now, using lens.draw_surface(surface_index=1) the user can view the surface shape in 2D or 3D :

Note that this computed using .sag(): for aspheric surfaces it includes the contributions from the base conic and the polynomials.
I'm not sure that the deviations are actually correct, they seem to small but maybe I'm mistaken?

[codecov]

Codecov Report

Attention: Patch coverage is 71.25000% with 23 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
optiland/visualization/visualization.py	70.00%	21 Missing ⚠️
optiland/optic.py	77.77%	2 Missing ⚠️

@@            Coverage Diff             @@
##           master      #59      +/-   ##
==========================================
- Coverage   96.93%   96.24%   -0.70%     
==========================================
  Files         105      105              
  Lines        6208     6421     +213     
==========================================
+ Hits         6018     6180     +162     
- Misses        190      241      +51     

Files with missing lines	Coverage Δ
optiland/visualization/__init__.py	100.00% <100.00%> (ø)
optiland/optic.py	97.19% <77.77%> (-1.34%)	⬇️
optiland/visualization/visualization.py	84.52% <70.00%> (-10.38%)	⬇️

... and 13 files with indirect coverage changes

[HarrisonKramer]

Solid addition!

We can merge after you update the limits in the call to sag in the SurfaceViewer class.

However, you mention that deviation to the sphere is crucial, but you're plotting deviation to the plane. Would it make sense to add an argument to SurfaceViewer.view like subtract_sphere? This gives the user the option to subtract a sphere or a plane (i.e., subtract nothing). You'd have to add a condition for the title too (deviation from sphere vs. plane).

Here are some (untested) methods that you could incorporate into SurfaceViewer to compute the deviation from the best-fit sphere. Perhaps there are better approaches.

import numpy as np
from scipy.optimize import minimize

class SurfaceViewer:

    # insert other methods here...

    @staticmethod
    def _sphere_sag(x, y, R):
        """Compute the sag of a sphere with radius R.
        
        Args:
            x, y: 2D arrays of coordinates.
            R: Sphere radius.

        Returns:
            2D array of sag values.
        """
        return R - np.sqrt(R**2 - x**2 - y**2)

    def _best_fit_sphere(x, y, z):
        """Find the best-fit sphere radius
        
        Args:
            x, y: 2D arrays of coordinates.
            z: 2D array of sags.

        Returns:
            Optimal sphere radius.
        """
        def error(R):
            z_s = self._sphere_sag(x, y, R)
            return np.sum((z - z_s) ** 2)  # RMS error

        R_guess = np.max(np.sqrt(x**2 + y**2))  # Initial guess
        res = minimize(error, R_guess)
        return res.x[0]

    def _compute_deviation(x, y, z):
        """Compute deviation from the best-fit sphere.
        
        Args:
            x, y: 2D arrays of coordinates.
            z: 2D array of sags.

        Returns:
            2D array of deviation values.
        """
        R = self._best_fit_sphere(x, y, z)
        z_s = self._sphere_sag(x, y, R)
        return z - z_s

[HarrisonKramer]

• You'll need to run self.optic.update_paraxial, then retrieve surface.semi_aperture to get the full extent to use for x & y.
• sag accepts (non-normalized) lens units, so your linear space should span (-surface.semi_aperture, surface.semi_aperture)

[drpaprika]

Following your comment the view() method uses surface.semi_aperture to compute the sag. The figure is now in non-normalized coordinates

I'll take a look at plotting the deviation to the BFS later

[HarrisonKramer]

Sounds good. I leave it up to you whether you want to add the deviation to the BFS in this PR or another. I'll quickly review the code now, then should be good to go.

[HarrisonKramer]

It looks like semi_aperture is always passed as an argument, so the "else" condition will never occur. If that's the case, then you can remove the if statement entirely.

[HarrisonKramer]

Feel free to merge when you're ready to do so. Thanks again.

[drpaprika]

I update the code with an option to display the deviation to the BFS. I runs without error but the result is surprising, as the deviation to a plane is lower than a deviation to a BFS:

The radius of the surface is $R =-100mm$, however the radius of the BFS as computed by _compute_deviation_to_best_fit_sphere() is $R_{BFS} = 5mm$.

There is probably a mistake here, so this PR is not ready to merge before it has been solved

[HarrisonKramer]

I quickly looked at this. First, I didn't have the sag equation quite right. It was missing a sign factor to account for negative radii of curvature. This seems to work as I expect:

return R - np.sign(R) * np.sqrt(R**2 - x**2 - y**2)

Then, it seems to help when changing the algorithm to Nelder-Mead and improving the initial guess:

initial_guess = np.mean(z + (x**2 + y**2) / (2 * z))  # first order approx
res = optimize.minimize(rms_error, initial_guess, method='Nelder-Mead')

I do not know if this will be generally robust, but it retrieves the correct radius of curvature for every surface of the Cooke triplet. By the way, I see it either fails or plots empty axes if you use the object or image surface. Maybe good to catch and prevent that, or update the code to handle it.

[drpaprika]

Implemented your suggestions, the result looks better. However the radius of the bfs is -1500mm while the radius of the surface is -100mm

I'll take a closer look when I have a moment :)

[HarrisonKramer]

I wonder if this has something to do with the Zernike geometry's sag definition. I see you used normalized radial coordinates (rho) for the base conic, so your BFS might be closer to the nominal radius * norm_radius.

Now that I think about this - the base class of the Zernike geometry, NewtonRaphsonGeometry, actually assumes the base conic is defined with a non-normalized radial component. It should work fine still, but maybe it converges more slowly. Do you think we should quickly change that? I don't know if there are other implications, or if it doesn't follow typical freeform standards. We can put that into a separate issue if needed.

[drpaprika]

I tried to see if the results are improved by fitting a conic instead of a sphere : not really. I don't see an issue if we delete this to keep a best fit sphere, it is exploratory.

Regarding your point I think that it's better to have the same convention everywhere to avoid unecessary complexity. Please implement the changes you have in mind, maybe I'm missing something.

I'm still not convinced that the plotted deviations have any physical meaning at the moment, I suggest to not close this PR until we have a better idea of what's happening.