Comparison between the OpticStudio and Optiland backends

from __future__ import annotations

from dataclasses import asdict
from typing import TYPE_CHECKING, Literal

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.colors import CenteredNorm
from mpl_toolkits.axes_grid1 import make_axes_locatable

import visisipy
from visisipy.backend import BackendSettings
from visisipy.opticstudio import OpticStudioBackend
from visisipy.optiland import OptilandBackend

if TYPE_CHECKING:
    from visisipy.analysis.cardinal_points import CardinalPointsResult
    from visisipy.analysis.refraction import FourierPowerVectorRefraction

Initialize both backends

BACKEND_SETTINGS = BackendSettings(
    field_type="angle",
    fields=[(0, 0)],
    wavelengths=[0.543],
    aperture_type="float_by_stop_size",
    aperture_value=3.0,
)

OpticStudioBackend.initialize(**BACKEND_SETTINGS, mode="standalone", ray_aiming="off")
OptilandBackend.initialize(**BACKEND_SETTINGS)

Define an eye model. Instead of the default material model (which is a fitted model for all wavelengths), we will use the refractive indices of the Navarro model for light of 543 nm.

model = visisipy.EyeModel(geometry=visisipy.models.NavarroGeometry(), materials=visisipy.models.NavarroMaterials543())
model_myopic = visisipy.EyeModel(
    visisipy.models.create_geometry(axial_length=model.geometry.axial_length + 2), visisipy.models.NavarroMaterials543()
)
model_hyperopic = visisipy.EyeModel(
    visisipy.models.create_geometry(axial_length=model.geometry.axial_length - 2), visisipy.models.NavarroMaterials543()
)

Cardinal points analysis

Calculate the cardinal point locations and the difference between the two backends.

def cardinal_points_to_dataframe(cardinal_points: CardinalPointsResult):
    """
    Convert cardinal points to a pandas DataFrame.
    """
    cardinal_points_dict = {
        k: (float("nan"), float("nan")) if v is NotImplemented else v for k, v in asdict(cardinal_points).items()
    }

    return pd.DataFrame.from_dict(cardinal_points_dict, orient="index", columns=["object", "image"])

opticstudio_cardinal_points = visisipy.analysis.cardinal_points(model=model, backend=OpticStudioBackend)
optiland_cardinal_points = visisipy.analysis.cardinal_points(model=model, backend=OptilandBackend)

cardinal_points_comparison = cardinal_points_to_dataframe(opticstudio_cardinal_points).join(
    cardinal_points_to_dataframe(optiland_cardinal_points), lsuffix="_opticstudio", rsuffix="_optiland"
)

cardinal_points_comparison["object_difference"] = (
    cardinal_points_comparison["object_optiland"] - cardinal_points_comparison["object_opticstudio"]
)
cardinal_points_comparison["image_difference"] = (
    cardinal_points_comparison["image_optiland"] - cardinal_points_comparison["image_opticstudio"]
)

cardinal_points_comparison

	object_opticstudio	image_opticstudio	object_optiland	image_optiland	object_difference	image_difference
focal_lengths	-16.467904	22.029115	-16.467899	22.029115	0.000005	3.574475e-07
focal_points	-14.885414	0.000014	-7.253626	0.000014	7.631788	-2.605261e-07
principal_points	1.582490	-22.029102	9.214274	-22.029102	7.631784	3.820264e-07
anti_principal_points	-31.353319	22.029129	-23.721525	22.029129	7.631794	9.692146e-08
nodal_points	7.143701	-16.467890	14.775490	-16.467886	7.631789	4.367338e-06
anti_nodal_points	-36.914530	16.467918	-29.282741	16.467913	7.631789	-4.888390e-06

FFT PSF analysis

Calculate FFT PSFs for both backends and compare the results. The difference plots show the relative difference, but PSF values smaller than 0.001 are not shown to avoid zero division errors.

def plot_dataframe(ax: plt.Axes, df: pd.DataFrame, title: str, cbar_label: str = "Relative intensity", **kwargs):
    im = ax.imshow(
        df.values,
        extent=(df.columns[0], df.columns[-1], df.index[-1], df.index[0]),
        **kwargs,
    )

    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)

    ax.set_title(title)
    ax.set_xlabel("X [μm]")
    ax.set_ylabel("Y [μm]")

    plt.colorbar(im, label=cbar_label, cax=cax)


def compare_fft_psfs(
    model: visisipy.EyeModel,
    field: tuple[float, float],
    sampling: int = 128,
):
    """
    Compare FFT PSFs for OpticStudio and Optiland backends.
    """
    shift_max_peak_distance = 10  # Maximum distance to shift the PSF peak for alignment

    fft_psf_opticstudio = visisipy.analysis.fft_psf(
        model=model,
        sampling=sampling,
        field_coordinate=field,
        backend=OpticStudioBackend,
    )

    fft_psf_optiland = visisipy.analysis.fft_psf(
        model=model,
        sampling=sampling,
        field_coordinate=field,
        backend=OptilandBackend,
    )

    fig, ax = plt.subplots(1, 3, figsize=(14, 4), constrained_layout=True)
    plot_dataframe(ax[0], fft_psf_opticstudio, "OpticStudio FFT PSF")
    plot_dataframe(ax[1], fft_psf_optiland, "Optiland FFT PSF")

    peak_opticstudio = np.unravel_index(np.argmax(fft_psf_opticstudio), fft_psf_opticstudio.shape)
    peak_optiland = np.unravel_index(np.argmax(fft_psf_optiland), fft_psf_optiland.shape)

    # Shift the PSFs to align their peaks if the peaks are close
    peak_difference = (peak_optiland[0] - peak_opticstudio[0], peak_optiland[1] - peak_opticstudio[1])
    print(peak_difference)
    if abs(peak_difference[0]) < shift_max_peak_distance and abs(peak_difference[1]) < shift_max_peak_distance:
        # Align the peaks of the two PSFs
        psf_opticstudio_shifted = np.roll(
            fft_psf_opticstudio.values,
            (peak_optiland[0] - peak_opticstudio[0], peak_optiland[1] - peak_opticstudio[1]),
            axis=(0, 1),
        )
    else:
        # If peaks are not close, use a default value based on sampling differences
        psf_opticstudio_shifted = np.roll(
            fft_psf_opticstudio.values,
            (0, 1),
            axis=(0, 1),
        )

    comparison = pd.DataFrame(
        psf_opticstudio_shifted - fft_psf_optiland.values,
        index=fft_psf_opticstudio.index,
        columns=fft_psf_opticstudio.columns,
    )

    plot_dataframe(ax[2], comparison, "Difference", "OpticStudio - Optiland", cmap="coolwarm", norm=CenteredNorm())

    fig.suptitle(f"FFT PSF Comparison for x={field[0]}°, y={field[1]}°")

compare_fft_psfs(model, field=(0, 0))

(np.int64(0), np.int64(1))

compare_fft_psfs(model, field=(0, 30))

(np.int64(0), np.int64(67))

compare_fft_psfs(model, field=(30, 30), sampling=256)

(np.int64(173), np.int64(-54))

Huygens PSF analysis

Calculate Huygens PSFs for both backends and plot them side by side. Because both backends use different sampling strategies and PSF extents, the PSFs are not directly comparable, but they should be qualitatively similar.

def compare_huygens_psfs(
    model: visisipy.EyeModel, field: tuple[float, float], pupil_sampling: int = 128, image_sampling: int = 128
):
    """
    Compare Huygens PSFs for OpticStudio and Optiland backends.
    """
    huygens_psf_opticstudio = visisipy.analysis.huygens_psf(
        model=model,
        field_coordinate=field,
        pupil_sampling=pupil_sampling,
        image_sampling=image_sampling,
        backend=OpticStudioBackend,
    )

    huygens_psf_optiland = visisipy.analysis.huygens_psf(
        model=model,
        field_coordinate=field,
        pupil_sampling=pupil_sampling,
        image_sampling=image_sampling,
        backend=OptilandBackend,
    )

    fig, ax = plt.subplots(1, 2, figsize=(10, 4), constrained_layout=True)
    plot_dataframe(ax[0], huygens_psf_opticstudio, "OpticStudio Huygens PSF")
    plot_dataframe(ax[1], huygens_psf_optiland, "Optiland Huygens PSF")

    fig.suptitle(f"Huygens PSF Comparison for x={field[0]}°, y={field[1]}°")

compare_huygens_psfs(model, field=(0, 0))

compare_huygens_psfs(model, field=(0, 30), pupil_sampling=128, image_sampling=256)

compare_huygens_psfs(model_myopic, field=(0, 0))

Raytracing analysis

Perform a number of ray traces in both backends and plot the results on top of each other.

RAYTRACE_FIELDS = [(0, y) for y in np.arange(0, 90, step=5).astype(float)]

raytrace_opticstudio = visisipy.analysis.raytrace(
    model=model,
    coordinates=RAYTRACE_FIELDS,
    wavelengths=[0.543],
    field_type="angle",
    pupil=(0, 0),
    backend=OpticStudioBackend,
)

raytrace_optiland = visisipy.analysis.raytrace(
    model=model,
    coordinates=RAYTRACE_FIELDS,
    wavelengths=[0.543],
    field_type="angle",
    pupil=(0, 0),
    backend=OptilandBackend,
)

# Optiland uses different coordinate conventions
raytrace_optiland["z"] -= model.geometry.cornea_thickness + model.geometry.anterior_chamber_depth

Plot the ray trace results on top of each other.

def plot_raytrace_result(raytrace_result: pd.DataFrame, ax: plt.Axes, color="blue", linestyle="-"):
    for _, df in raytrace_result.groupby("field"):
        ax.plot(df.z, df.y, color=color, linestyle=linestyle)


fig, ax = plt.subplots()

visisipy.plots.plot_eye(ax=ax, geometry=model.geometry, lens_edge_thickness=0.5)

plot_raytrace_result(raytrace_opticstudio, ax=ax)
plot_raytrace_result(raytrace_optiland, ax=ax, color="red", linestyle=":")

ax.set_xlim(-5, 25)
ax.set_ylim(-15, 15)
ax.set_aspect("equal")

Calculate the distances between the ray trace results point by point and plot them as a function of the Z-coordinate and the surface index.

_select_columns = ["field", "index", "x", "y", "z"]

raytrace_comparison = pd.merge(
    raytrace_opticstudio[_select_columns],
    raytrace_optiland[_select_columns],
    on=["field", "index"],
    suffixes=("_opticstudio", "_optiland"),
).query("index != 0")

raytrace_comparison.eval(
    "distance = sqrt((x_optiland - x_opticstudio)**2 + (y_optiland - y_opticstudio)**2)", inplace=True
)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 5), sharey=True)

_field_coordinates = [f[1] for f in raytrace_comparison.field]
m = ax[0].scatter(raytrace_comparison["index"], raytrace_comparison.distance, c=_field_coordinates)

ax[0].set_xlabel("surface #")
ax[0].set_ylabel("difference [mm]")
ax[0].set_title("Difference by surface number")

ax[1].scatter(raytrace_comparison.z_opticstudio, raytrace_comparison.distance, c=_field_coordinates)

ax[1].set_xlabel("z [mm]")
ax[1].set_title("Difference by z coordinate")

plt.colorbar(m, ax=ax[1], label="field angle [deg]")

<matplotlib.colorbar.Colorbar at 0x2a8a5d61370>

Strehl ratio analysis

Calculate the Strehl ratio for both backends and all supported PSF types, and compare the results.

psf_types: list[Literal["fft", "huygens"]] = ["fft", "huygens"]

strehl_ratios = {}

for backend, name in zip([OpticStudioBackend, OptilandBackend], ["OpticStudio", "Optiland"]):
    strehl_ratios[name] = {}

    for psf_type in psf_types:
        try:
            strehl_ratios[name][psf_type] = visisipy.analysis.strehl_ratio(
                model=model,
                field_coordinate=(0, 0),
                psf_type=psf_type,
                backend=backend,
            )
        except NotImplementedError:
            strehl_ratios[name][psf_type] = np.nan

pd.DataFrame(strehl_ratios).style.format(precision=3)

	OpticStudio	Optiland
fft	nan	0.677
huygens	0.671	0.276

Zernike standard coefficients

Calculate the Zernike standard coefficients for the Navarro eye model at different field angles. Note that the y-axes in the plots are logarithmic.

zernike_central_opticstudio = visisipy.analysis.zernike_standard_coefficients(
    model=model,
    field_coordinate=(0, 0),
    wavelength=0.543,
    sampling=128,
    maximum_term=45,
    field_type="angle",
    backend=OpticStudioBackend,
)

zernike_central_optiland = visisipy.analysis.zernike_standard_coefficients(
    model=model,
    field_coordinate=(0, 0),
    wavelength=0.543,
    sampling=128,
    maximum_term=45,
    field_type="angle",
    backend=OptilandBackend,
)

def compare_zernike_coefficients(
    zernike_coefficients_opticstudio: dict,
    zernike_coefficients_optiland: dict,
):
    """
    Compare Zernike coefficients from OpticStudio and Optiland.
    """
    differences = {k: (v - zernike_coefficients_optiland[k]) for k, v in zernike_coefficients_opticstudio.items()}

    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 5))
    ax[0].bar(
        zernike_coefficients_opticstudio.keys(),
        zernike_coefficients_opticstudio.values(),
        width=0.5,
        label="OpticStudio",
    )
    ax[0].bar(
        np.fromiter(zernike_coefficients_optiland.keys(), dtype=int) + 0.5,
        zernike_coefficients_optiland.values(),
        width=0.5,
        label="Optiland",
    )

    ax[0].set_yscale("log")

    ax[0].set_title("Zernike coefficients")
    ax[0].set_xlabel("Coefficient")
    ax[0].set_ylabel("Value")
    ax[0].legend()

    ax[1].bar(differences.keys(), differences.values(), width=0.5, label="Difference")
    ax[1].axhline(0, color="black", linestyle="--", linewidth=0.5)
    ax[1].set_title("Difference between OpticStudio and Optiland")

    ax[1].set_yscale("log")

    ax[1].set_xlabel("Coefficient")
    ax[1].set_ylabel("Difference")


compare_zernike_coefficients(zernike_central_opticstudio, zernike_central_optiland)

This is an emmetropic eye, so the aberrations are very small. They will be slightly larger at a nonzero eccentricity:

zernike_10deg_opticstudio = visisipy.analysis.zernike_standard_coefficients(
    model=model,
    field_coordinate=(0, 10),
    wavelength=0.543,
    sampling=128,
    maximum_term=45,
    field_type="angle",
    backend=OpticStudioBackend,
)

zernike_10deg_optiland = visisipy.analysis.zernike_standard_coefficients(
    model=model,
    field_coordinate=(0, 10),
    wavelength=0.543,
    sampling=128,
    maximum_term=45,
    field_type="angle",
    backend=OptilandBackend,
)

compare_zernike_coefficients(zernike_10deg_opticstudio, zernike_10deg_optiland)

zernike_myopic_opticstudio = visisipy.analysis.zernike_standard_coefficients(
    model=model_myopic,
    field_coordinate=(0, 0),
    wavelength=0.543,
    sampling=128,
    maximum_term=45,
    field_type="angle",
    backend=OpticStudioBackend,
)

zernike_myopic_optiland = visisipy.analysis.zernike_standard_coefficients(
    model=model_myopic,
    field_coordinate=(0, 0),
    wavelength=0.543,
    sampling=128,
    maximum_term=45,
    field_type="angle",
    backend=OptilandBackend,
)

compare_zernike_coefficients(zernike_myopic_opticstudio, zernike_myopic_optiland)

zernike_hyperopic_opticstudio = visisipy.analysis.zernike_standard_coefficients(
    model=model_hyperopic,
    field_coordinate=(0, 0),
    wavelength=0.543,
    sampling=128,
    maximum_term=45,
    field_type="angle",
    backend=OpticStudioBackend,
)

zernike_hyperopic_optiland = visisipy.analysis.zernike_standard_coefficients(
    model=model_hyperopic,
    field_coordinate=(0, 0),
    wavelength=0.543,
    sampling=128,
    maximum_term=45,
    field_type="angle",
    backend=OptilandBackend,
)

compare_zernike_coefficients(zernike_hyperopic_opticstudio, zernike_hyperopic_optiland)

RMS-HOA

Calculate the root-mean-square of the higher-order aberrations (RMS-HOA) for the Navarro eye model on the previously defined eye models.

rms_hoa_parameters = {
    "emmetropic": {
        "model": model,
        "field": (0, 0),
    },
    "myopic": {
        "model": model_myopic,
        "field": (0, 0),
    },
    "myopic_10deg": {
        "model": model_myopic,
        "field": (0, 10),
    },
    "hyperopic": {
        "model": model_hyperopic,
        "field": (0, 0),
    },
}

rms_hoa_comparison = []

for name, parameters in rms_hoa_parameters.items():
    rms_hoa_opticstudio = visisipy.analysis.rms_hoa(
        model=parameters["model"],
        field_coordinate=parameters["field"],
        wavelength=0.543,
        sampling=128,
        maximum_term=45,
        field_type="angle",
        backend=OpticStudioBackend,
    )

    rms_hoa_optiland = visisipy.analysis.rms_hoa(
        model=parameters["model"],
        field_coordinate=parameters["field"],
        wavelength=0.543,
        sampling=128,
        maximum_term=45,
        field_type="angle",
        backend=OptilandBackend,
    )

    rms_hoa_comparison.append(
        {
            "name": name,
            "opticstudio": rms_hoa_opticstudio,
            "optiland": rms_hoa_optiland,
            "difference": rms_hoa_opticstudio - rms_hoa_optiland,
        }
    )

pd.DataFrame(rms_hoa_comparison).set_index("name")

	opticstudio	optiland	difference
name
emmetropic	0.025453	0.025453	1.659608e-09
myopic	0.029329	0.029329	-5.773978e-09
myopic_10deg	0.113635	0.113635	2.679752e-07
hyperopic	0.020970	0.020970	-4.476286e-10

Refraction

Calculate the on-axis and off-axis refraction using both backends.

refraction_onaxis_opticstudio = visisipy.analysis.refraction(
    model=model,
    field_coordinate=(0, 0),
    sampling=128,
    wavelength=0.543,
    field_type="angle",
    backend=OpticStudioBackend,
)

refraction_onaxis_optiland = visisipy.analysis.refraction(
    model=model,
    field_coordinate=(0, 0),
    sampling=128,
    wavelength=0.543,
    field_type="angle",
    backend=OptilandBackend,
)

def compare_refractions(
    refraction_opticstudio: FourierPowerVectorRefraction,
    refraction_optiland: FourierPowerVectorRefraction,
):
    """
    Compare refractions from OpticStudio and Optiland.
    """
    df = pd.DataFrame(
        {
            "OpticStudio": asdict(refraction_opticstudio),
            "Optiland": asdict(refraction_optiland),
        }
    )
    df["Difference"] = df["OpticStudio"] - df["Optiland"]

    return df


compare_refractions(refraction_onaxis_opticstudio, refraction_onaxis_optiland)

	OpticStudio	Optiland	Difference
M	0.000046	4.546456e-05	1.316288e-07
J0	0.000000	-7.418248e-12	7.418248e-12
J45	0.000000	-1.580311e-12	1.580311e-12

refraction_offaxis_opticstudio = visisipy.analysis.refraction(
    model=model,
    field_coordinate=(0, 10),
    sampling=128,
    wavelength=0.543,
    field_type="angle",
    backend=OpticStudioBackend,
)

refraction_offaxis_optiland = visisipy.analysis.refraction(
    model=model,
    field_coordinate=(0, 10),
    sampling=128,
    wavelength=0.543,
    field_type="angle",
    backend=OptilandBackend,
)

compare_refractions(refraction_offaxis_opticstudio, refraction_offaxis_optiland)

	OpticStudio	Optiland	Difference
M	0.09206	9.205975e-02	-9.527374e-09
J0	0.11967	1.196702e-01	8.005518e-08
J45	0.00000	2.542934e-13	-2.542934e-13

refraction_myopic_opticstudio = visisipy.analysis.refraction(
    model=model_myopic,
    field_coordinate=(0, 0),
    sampling=128,
    wavelength=0.543,
    field_type="angle",
    backend=OpticStudioBackend,
)

refraction_myopic_optiland = visisipy.analysis.refraction(
    model=model_myopic,
    field_coordinate=(0, 0),
    sampling=128,
    wavelength=0.543,
    field_type="angle",
    backend=OptilandBackend,
)

compare_refractions(refraction_myopic_opticstudio, refraction_myopic_optiland)

	OpticStudio	Optiland	Difference
M	-5.944698	-5.944698e+00	6.935850e-08
J0	0.000000	-5.859356e-12	5.859356e-12
J45	0.000000	-1.092574e-12	1.092574e-12