diff --git a/optika/apertures/__init__.py b/optika/apertures/__init__.py index 084c24d..efbeaf2 100644 --- a/optika/apertures/__init__.py +++ b/optika/apertures/__init__.py @@ -16,6 +16,7 @@ OctagonalAperture, AbstractIsoscelesTrapezoidalAperture, IsoscelesTrapezoidalAperture, + footprint_aperture, ) __all__ = [ @@ -32,4 +33,5 @@ "OctagonalAperture", "AbstractIsoscelesTrapezoidalAperture", "IsoscelesTrapezoidalAperture", + "footprint_aperture", ] diff --git a/optika/apertures/_apertures.py b/optika/apertures/_apertures.py index 4d45833..8aeae92 100644 --- a/optika/apertures/_apertures.py +++ b/optika/apertures/_apertures.py @@ -2,6 +2,7 @@ import dataclasses import numpy as np import numpy.typing as npt +import scipy.spatial import matplotlib.axes import matplotlib.lines import matplotlib.pyplot as plt @@ -24,9 +25,80 @@ "OctagonalAperture", "AbstractIsoscelesTrapezoidalAperture", "IsoscelesTrapezoidalAperture", + "footprint_aperture", ] +def footprint_aperture( + position: na.AbstractCartesian2dVectorArray | na.AbstractCartesian3dVectorArray, + where: None | bool | na.AbstractScalar = None, +) -> "PolygonalAperture": + """ + Construct a polygonal aperture from the convex hull of a beam footprint. + + This is useful for building the effective pupil of one channel of a + multi-channel instrument: trace rays through the full system (including + any obscurations) with :meth:`optika.systems.SequentialSystem.footprint`, + then convert the illuminated portion of a surface into an explicit + aperture for a physical-optics model of that channel. + + Parameters + ---------- + position + The positions of the rays intersecting the surface, + expressed in local surface coordinates. + where + An optional boolean mask selecting the unvignetted rays. + + Notes + ----- + The footprint is approximated by its convex hull, which is exact only + for convex footprints. + Where an edge of the hull comes from the projection of a downstream + aperture (instead of a physical edge on this surface), light diffracted + past that projection is clipped, which slightly underestimates the + diffraction pattern of the downstream aperture. + """ + position = position.xy + + if where is None: + where = True + where = na.broadcast_to( + na.as_named_array(where), + na.shape_broadcasted(position.x, position.y, where), + ) + + shape = where.shape + x = na.broadcast_to(position.x, shape).ndarray.ravel() + y = na.broadcast_to(position.y, shape).ndarray.ravel() + mask = where.ndarray.ravel() + + x = x[mask] + y = y[mask] + + unit = na.unit(x) + if unit is not None: + x = x.to_value(unit) + y = y.to(unit).value + + points = np.stack([x, y], axis=~0) + hull = scipy.spatial.ConvexHull(points) + + x = x[hull.vertices] + y = y[hull.vertices] + z = np.zeros_like(x) + if unit is not None: + x, y, z = x << unit, y << unit, z << unit + + return PolygonalAperture( + vertices=na.Cartesian3dVectorArray( + x=na.ScalarArray(x, axes=("vertex",)), + y=na.ScalarArray(y, axes=("vertex",)), + z=na.ScalarArray(z, axes=("vertex",)), + ), + ) + + @dataclasses.dataclass(eq=False, repr=False) class AbstractAperture( optika.mixins.DxfWritable, diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index 9da4dbf..8e56ee5 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -746,6 +746,80 @@ def rayfunction_default(self) -> optika.rays.RayFunctionArray: """ return self.rayfunction() + def footprint( + self, + surface: optika.surfaces.AbstractSurface, + wavelength: None | u.Quantity | na.AbstractScalar = None, + field: None | na.AbstractCartesian2dVectorArray = None, + num: int = 101, + normalized_field: bool = True, + ) -> tuple[na.AbstractCartesian3dVectorArray, na.AbstractScalar]: + """ + The beam footprint on a given surface: the positions where a dense + grid of rays traced through the entire system intersects the surface, + expressed in local surface coordinates, along with a mask selecting + the rays that reach the sensor unvignetted. + + This identifies the used portion of each optic, which is useful for + restricting the sampling region of a physical-optics calculation + (see :meth:`wavefield`) or for constructing the effective pupil of + one channel of a multi-channel instrument + (see :func:`optika.apertures.footprint_aperture`). + + Parameters + ---------- + surface + The surface on which to evaluate the footprint. + Must be an element of :attr:`surfaces`, :attr:`object`, + or the :attr:`sensor`. + wavelength + The wavelengths of the rays. + If :obj:`None` (the default), ``self.grid_input.wavelength`` + will be used. + field + The field positions of the rays, in either normalized or + physical units. + If :obj:`None` (the default), ``self.grid_input.field`` + will be used. + num + The number of rays along each axis of the pupil grid. + normalized_field + A boolean flag indicating whether the `field` parameter is given + in normalized or physical units. + """ + surfaces = self.surfaces_all + index = [i for i, s in enumerate(surfaces) if s is surface] + if not index: + raise ValueError( + f"surface {surface.name!r} is not an element of this system." + ) + index = index[0] + + axis = "_footprint" + pupil = na.Cartesian2dVectorLinearSpace( + start=-1, + stop=1, + axis=na.Cartesian2dVectorArray(f"{axis}_x", f"{axis}_y"), + num=num, + centers=True, + ) + + raytrace = self.raytrace( + wavelength=wavelength, + field=field, + pupil=pupil, + axis=axis, + normalized_field=normalized_field, + ) + + position = raytrace.outputs.position[{axis: index}] + where = raytrace.outputs.unvignetted[{axis: ~0}] + + if surface.transformation is not None: + position = surface.transformation.inverse(position) + + return position, where + def wavefield( self, wavelength: None | u.Quantity | na.AbstractScalar = None, @@ -758,7 +832,10 @@ def wavefield( chunk_size: int = 1024, normalized_field: bool = True, seed: None | int = None, - bound: None | Sequence = None, + bound: None | str | Sequence = None, + padding_relative: float = 0.1, + padding_absolute: u.Quantity | na.AbstractScalar = 0 * u.mm, + num_footprint: int = 101, ) -> optika.wavefields.WavefieldFunctionArray: r""" Propagate a wavefield through this system using physical optics and @@ -814,12 +891,28 @@ def wavefield( bound An optional sequence of sampling regions, one for each surface with an aperture, where each element is either :obj:`None` (use - the bounding box of the surface's aperture) or a tuple of the - lower and upper corners of the region in local surface - coordinates. + the bounding box of the surface's aperture), the string + ``"footprint"`` (use the bounding box of the beam footprint from + a geometric raytrace, padded by `padding_relative` and + `padding_absolute`), or a tuple of the lower and upper corners + of the region in local surface coordinates. + The string ``"footprint"`` may also be given for `bound` itself, + which applies it to every surface. This is important for surfaces near an intermediate focus, where the wavefield is concentrated in a region much smaller than the aperture and would otherwise be unresolved by the samples. + padding_relative + The fraction of the footprint extent by which to expand a + ``"footprint"`` sampling region. + padding_absolute + The minimum distance by which to expand a ``"footprint"`` + sampling region. + At an intermediate focus the geometric footprint collapses to a + point, so this should be at least a few times the expected size + of the diffraction pattern there. + num_footprint + The number of rays along each axis of the pupil grid used to + compute the ``"footprint"`` sampling regions. Notes ----- @@ -838,6 +931,11 @@ def wavefield( filters) are often better excluded from the wave propagation (by setting their aperture to :obj:`None` in a copy of the system) than included as an undersampled propagation step. + The effect of an excluded mask on the pupil can be retained by + replacing the aperture of the nearest powered surface with the + convex hull of its illuminated footprint, + using :meth:`footprint` and + :func:`optika.apertures.footprint_aperture`. """ if wavelength is None: wavelength = self.grid_input.wavelength @@ -874,13 +972,46 @@ def wavefield( f"{len(surfaces)} surfaces with an aperture, got {len(num)}." ) - if bound is None: - bound = [None] * len(surfaces) + if (bound is None) or isinstance(bound, str): + bound = [bound] * len(surfaces) if len(bound) != len(surfaces): raise ValueError( f"`bound` must have one element for each of the " f"{len(surfaces)} surfaces with an aperture, got {len(bound)}." ) + bound = list(bound) + for k, surface in enumerate(surfaces): + if not isinstance(bound[k], str): + continue + if bound[k] != "footprint": + raise ValueError( + f"the only string allowed in `bound` is 'footprint', " + f"got {bound[k]!r}." + ) + position_footprint, where_footprint = self.footprint( + surface=surface, + wavelength=grid.wavelength, + field=grid.field, + num=num_footprint, + normalized_field=False, + ) + position_footprint = position_footprint.xy + unit = na.unit(position_footprint.x) + inf = np.inf if unit is None else np.inf * unit + axis_footprint = ("_footprint_x", "_footprint_y") + lower = na.Cartesian2dVectorArray( + x=np.where(where_footprint, position_footprint.x, +inf), + y=np.where(where_footprint, position_footprint.y, +inf), + ).min(axis_footprint) + upper = na.Cartesian2dVectorArray( + x=np.where(where_footprint, position_footprint.x, -inf), + y=np.where(where_footprint, position_footprint.y, -inf), + ).max(axis_footprint) + padding = np.maximum( + padding_relative * (upper - lower), + padding_absolute, + ) + bound[k] = (lower - padding, upper + padding) axes = [ (f"_{axis[0]}_{k % 2}", f"_{axis[1]}_{k % 2}") for k in range(len(surfaces)) @@ -1076,6 +1207,79 @@ def psf( The two logical axes of the grid of sensor positions. kwargs Additional keyword arguments passed to :meth:`wavefield`. + + Examples + -------- + Compute the Huygens PSF of an ideal parabolic telescope, + which should be an Airy pattern with its first minimum at a radius + of :math:`1.22 \lambda f / D`. + + .. jupyter-execute:: + + import matplotlib.pyplot as plt + import astropy.units as u + import astropy.visualization + import named_arrays as na + import optika + + primary = optika.surfaces.Surface( + name="primary", + sag=optika.sags.ParabolicSag(focal_length=-200 * u.mm), + aperture=optika.apertures.CircularAperture(20 * u.mm), + material=optika.materials.Mirror(), + is_pupil_stop=True, + transformation=na.transformations.Cartesian3dTranslation( + z=200 * u.mm, + ), + ) + + sensor = optika.sensors.ImagingSensor( + name="sensor", + width_pixel=10 * u.um, + axis_pixel=na.Cartesian2dVectorArray("detector_x", "detector_y"), + num_pixel=na.Cartesian2dVectorArray(64, 64), + timedelta_exposure=1 * u.s, + is_field_stop=True, + ) + + system = optika.systems.SequentialSystem( + surfaces=[primary], + sensor=sensor, + grid_input=optika.vectors.ObjectVectorArray( + wavelength=500 * u.nm, + field=na.Cartesian2dVectorLinearSpace( + start=-1, + stop=1, + axis=na.Cartesian2dVectorArray("field_x", "field_y"), + num=5, + centers=True, + ), + pupil=na.Cartesian2dVectorLinearSpace( + start=-1, + stop=1, + axis=na.Cartesian2dVectorArray("pupil_x", "pupil_y"), + num=5, + centers=True, + ), + ), + ) + + psf = system.psf( + field=na.Cartesian2dVectorArray(0, 0), + width_sensor=30 * u.um, + num=101**2, + num_sensor=51, + bound="footprint", + seed=42, + ) + + with astropy.visualization.quantity_support(): + plt.figure() + plt.gca().set_aspect("equal") + na.plt.pcolormesh( + psf.inputs.position, + C=psf.outputs, + ) """ wavefield = self.wavefield(axis=axis, **kwargs) diff --git a/optika/systems/_wavefield_test.py b/optika/systems/_wavefield_test.py index 8479507..5d5cfaa 100644 --- a/optika/systems/_wavefield_test.py +++ b/optika/systems/_wavefield_test.py @@ -448,3 +448,85 @@ def test_wavefield_samples_decentered_aperture(): assert x.max() > decenter + 0.9 * half_width assert y.min() < -0.9 * half_width assert y.max() > +0.9 * half_width + + +def test_footprint(): + """ + The beam footprint on the primary of the test telescope must fill its + aperture, and the footprint must be expressed in local coordinates. + """ + system = _telescope() + primary = system.surfaces[0] + + position, where = system.footprint( + surface=primary, + field=na.Cartesian2dVectorArray(0, 0), + num=21, + ) + + shape = na.shape_broadcasted(position.x, position.y, where) + x = na.broadcast_to(position.x, shape).ndarray + y = na.broadcast_to(position.y, shape).ndarray + mask = na.broadcast_to(where, shape).ndarray + + assert np.any(mask) + radius = np.sqrt(np.square(x[mask]) + np.square(y[mask])) + assert np.all(radius <= primary.aperture.radius) + assert radius.max() > 0.9 * primary.aperture.radius + + with pytest.raises(ValueError, match="not an element"): + system.footprint(surface=optika.surfaces.Surface(name="other")) + + +def test_footprint_aperture(): + """ + The convex hull of the footprint on the primary must accept points + inside the beam and reject points outside it. + """ + system = _telescope() + primary = system.surfaces[0] + + position, where = system.footprint( + surface=primary, + field=na.Cartesian2dVectorArray(0, 0), + num=21, + ) + aperture = optika.apertures.footprint_aperture(position, where) + + assert isinstance(aperture, optika.apertures.PolygonalAperture) + point_inside = na.Cartesian3dVectorArray(0, 0, 0) * u.mm + point_outside = na.Cartesian3dVectorArray(100, 0, 0) * u.mm + assert np.all(na.as_named_array(aperture(point_inside))) + assert not np.any(na.as_named_array(aperture(point_outside))) + + +def test_psf_footprint_bound(): + """ + The Huygens PSF computed with automatic footprint sampling regions must + match the PSF computed with the aperture bounding box. + """ + system = _telescope() + kwargs = dict( + field=na.Cartesian2dVectorArray(0, 0), + width_sensor=0.05 * u.mm, + num=51**2, + num_sensor=33, + seed=0, + ) + + psf_default = system.psf(**kwargs) + psf_footprint = system.psf(bound="footprint", **kwargs) + + def _ee50(psf): + return optika.wavefields.encircled_energy_radius( + intensity=psf.outputs, + position=psf.inputs.position, + axis=("wavefield_x", "wavefield_y"), + ) + + ee50_default = _ee50(psf_default) + ee50_footprint = _ee50(psf_footprint) + assert np.abs(ee50_footprint - ee50_default) < 0.2 * ee50_default + + with pytest.raises(ValueError, match="footprint"): + system.psf(bound="banana", **kwargs) diff --git a/optika/wavefields/__init__.py b/optika/wavefields/__init__.py index 843fce9..b174c8f 100644 --- a/optika/wavefields/__init__.py +++ b/optika/wavefields/__init__.py @@ -8,6 +8,7 @@ WavefieldFunctionArray, ) from ._propagation import rayleigh_sommerfeld +from ._metrics import encircled_energy_radius, ensquared_energy, fwhm __all__ = [ "AbstractWavefieldVectorArray", @@ -15,4 +16,7 @@ "AbstractWavefieldFunctionArray", "WavefieldFunctionArray", "rayleigh_sommerfeld", + "encircled_energy_radius", + "ensquared_energy", + "fwhm", ] diff --git a/optika/wavefields/_metrics.py b/optika/wavefields/_metrics.py new file mode 100644 index 0000000..fc021f0 --- /dev/null +++ b/optika/wavefields/_metrics.py @@ -0,0 +1,190 @@ +"""Scalar image-quality metrics of sampled point-spread functions.""" + +import numpy as np +import astropy.units as u +import named_arrays as na + +__all__ = [ + "encircled_energy_radius", + "ensquared_energy", + "fwhm", +] + + +def _centroid( + intensity: na.AbstractScalar, + position: na.AbstractCartesian2dVectorArray, + axis: tuple[str, str], +) -> na.AbstractCartesian2dVectorArray: + return (position * intensity).sum(axis) / intensity.sum(axis) + + +def encircled_energy_radius( + intensity: na.AbstractScalar, + position: na.AbstractCartesian2dVectorArray | na.AbstractCartesian3dVectorArray, + axis: tuple[str, str], + fraction: float = 0.5, +) -> u.Quantity | na.AbstractScalar: + """ + The radius from the intensity-weighted centroid which encloses a given + fraction of the total energy of a sampled point-spread function. + + Parameters + ---------- + intensity + The intensity of the point-spread function at each sample. + position + The position of each sample. + axis + The two logical axes of the grid of samples. + fraction + The fraction of the total energy to enclose. + """ + if isinstance(position, na.AbstractCartesian3dVectorArray): + position = position.xy + centroid = _centroid(intensity, position, axis) + radius = (position - centroid).length + + shape = na.shape_broadcasted(intensity, radius) + intensity = na.broadcast_to(intensity, shape) + radius = na.broadcast_to(radius, shape) + + shape_outer = {ax: n for ax, n in shape.items() if ax not in axis} + unit = na.unit(radius) + result = na.ScalarArray( + ndarray=np.empty(tuple(shape_outer.values())), + axes=tuple(shape_outer), + ) + if unit is not None: + result = result * unit + + for index in na.ndindex(shape_outer): + r = radius[index].ndarray.ravel() + w = intensity[index].ndarray.ravel() + order = np.argsort(r) + energy = np.cumsum(w[order]) + energy = energy / energy[~0] + result[index] = r[order][np.searchsorted(energy, fraction)] + + if not shape_outer: + result = result.ndarray + + return result + + +def ensquared_energy( + intensity: na.AbstractScalar, + position: na.AbstractCartesian2dVectorArray | na.AbstractCartesian3dVectorArray, + axis: tuple[str, str], + width: u.Quantity | na.AbstractScalar, +) -> na.AbstractScalar: + """ + The fraction of the total energy of a sampled point-spread function + inside a square of the given width centered on the intensity-weighted + centroid. + + This is useful for checking if an optical system is pixel-limited: + if the ensquared energy in one pixel is large, the detector + (and not the optics) sets the resolution. + + Parameters + ---------- + intensity + The intensity of the point-spread function at each sample. + position + The position of each sample. + axis + The two logical axes of the grid of samples. + width + The width of the square. + """ + if isinstance(position, na.AbstractCartesian3dVectorArray): + position = position.xy + centroid = _centroid(intensity, position, axis) + offset = position - centroid + + inside = (np.abs(offset.x) < width / 2) & (np.abs(offset.y) < width / 2) + + return (intensity * inside).sum(axis) / intensity.sum(axis) + + +def fwhm( + intensity: na.AbstractScalar, + position: na.AbstractCartesian2dVectorArray | na.AbstractCartesian3dVectorArray, + axis: tuple[str, str], +) -> na.AbstractCartesian2dVectorArray: + """ + The full width at half maximum of a sampled point-spread function along + each axis, measured from the cuts through the brightest sample. + + Parameters + ---------- + intensity + The intensity of the point-spread function at each sample. + position + The position of each sample. + axis + The two logical axes of the grid of samples. + """ + if isinstance(position, na.AbstractCartesian3dVectorArray): + position = position.xy + + shape = na.shape_broadcasted(intensity, position.x, position.y) + intensity = na.broadcast_to(intensity, shape) + x = na.broadcast_to(position.x, shape) + y = na.broadcast_to(position.y, shape) + + shape_outer = {ax: n for ax, n in shape.items() if ax not in axis} + unit = na.unit(x) + + result = na.Cartesian2dVectorArray( + x=na.ScalarArray( + ndarray=np.empty(tuple(shape_outer.values())), + axes=tuple(shape_outer), + ), + y=na.ScalarArray( + ndarray=np.empty(tuple(shape_outer.values())), + axes=tuple(shape_outer), + ), + ) + if unit is not None: + result = result * unit + + def _fwhm_of_cut(cut: np.ndarray, coordinate: np.ndarray) -> float: + index_peak = int(np.argmax(cut)) + half = cut[index_peak] / 2 + left = np.interp( + half, + cut[: index_peak + 1], + coordinate[: index_peak + 1], + ) + right = np.interp( + half, + cut[index_peak:][::-1], + coordinate[index_peak:][::-1], + ) + return right - left + + for index in na.ndindex(shape_outer): + cut = intensity[index] + index_peak = np.unravel_index( + np.argmax(cut.ndarray), + cut.ndarray.shape, + ) + index_peak = dict(zip(cut.axes, index_peak)) + + cut_x = cut[{axis[1]: index_peak[axis[1]]}].ndarray + coordinate_x = x[index][{axis[1]: index_peak[axis[1]]}].ndarray + cut_y = cut[{axis[0]: index_peak[axis[0]]}].ndarray + coordinate_y = y[index][{axis[0]: index_peak[axis[0]]}].ndarray + + result.x[index] = _fwhm_of_cut(cut_x, coordinate_x) + result.y[index] = _fwhm_of_cut(cut_y, coordinate_y) + + if not shape_outer: + result = na.Cartesian2dVectorArray( + x=result.x.ndarray, + y=result.y.ndarray, + ) + + return result diff --git a/optika/wavefields/_tests/test_metrics.py b/optika/wavefields/_tests/test_metrics.py new file mode 100644 index 0000000..f53e9a2 --- /dev/null +++ b/optika/wavefields/_tests/test_metrics.py @@ -0,0 +1,106 @@ +"""Tests of the point-spread-function metrics against an analytic Gaussian.""" + +import numpy as np +import astropy.units as u +import named_arrays as na +import optika + + +def _gaussian_psf( + sigma: u.Quantity = 5 * u.um, + center: u.Quantity = 1 * u.um, + width: u.Quantity = 80 * u.um, + num: int = 401, +) -> tuple[na.AbstractScalar, na.Cartesian2dVectorArray]: + position = na.Cartesian2dVectorLinearSpace( + start=center - width / 2, + stop=center + width / 2, + axis=na.Cartesian2dVectorArray("psf_x", "psf_y"), + num=num, + ).explicit + r2 = np.square(position.x - center) + np.square(position.y - center) + intensity = np.exp(-r2 / (2 * np.square(sigma))) + return intensity, position + + +def test_encircled_energy_radius(): + """ + The 50% encircled-energy radius of a symmetric Gaussian is + :math:`\\sigma \\sqrt{2 \\ln 2}`. + """ + sigma = 5 * u.um + intensity, position = _gaussian_psf(sigma=sigma) + + result = optika.wavefields.encircled_energy_radius( + intensity=intensity, + position=position, + axis=("psf_x", "psf_y"), + fraction=0.5, + ) + + expected = sigma * np.sqrt(2 * np.log(2)) + assert u.isclose(result, expected, rtol=2e-2) + + +def test_encircled_energy_radius_broadcast(): + """ + The metric must support additional broadcast axes, for example a grid + of point-spread functions over the field of view. + """ + sigma = na.ScalarArray(np.array([3, 6]) * u.um, axes=("field",)) + intensity, position = _gaussian_psf() + r2 = np.square(position.x - 1 * u.um) + np.square(position.y - 1 * u.um) + intensity = np.exp(-r2 / (2 * np.square(sigma))) + + result = optika.wavefields.encircled_energy_radius( + intensity=intensity, + position=position, + axis=("psf_x", "psf_y"), + fraction=0.5, + ) + + assert result.shape == {"field": 2} + expected = sigma * np.sqrt(2 * np.log(2)) + assert np.all(np.abs(result - expected) < 0.05 * expected) + + +def test_ensquared_energy(): + """ + The energy of a symmetric Gaussian inside a centered square of width + :math:`w` is :math:`\\mathrm{erf}^2(w / 2 \\sqrt{2} \\sigma)`. + """ + import scipy.special + + sigma = 5 * u.um + width = 15 * u.um + intensity, position = _gaussian_psf(sigma=sigma) + + result = optika.wavefields.ensquared_energy( + intensity=intensity, + position=position, + axis=("psf_x", "psf_y"), + width=width, + ) + + arg = (width / (2 * np.sqrt(2) * sigma)).to_value(u.dimensionless_unscaled) + expected = np.square(scipy.special.erf(arg)) + assert np.abs(float(na.as_named_array(result).ndarray) - expected) < 2e-2 + + +def test_fwhm(): + """ + The full width at half maximum of a Gaussian is + :math:`2 \\sigma \\sqrt{2 \\ln 2}` along both axes. + """ + sigma = 5 * u.um + intensity, position = _gaussian_psf(sigma=sigma) + + result = optika.wavefields.fwhm( + intensity=intensity, + position=position, + axis=("psf_x", "psf_y"), + ) + + expected = 2 * sigma * np.sqrt(2 * np.log(2)) + assert u.isclose(result.x, expected, rtol=2e-2) + assert u.isclose(result.y, expected, rtol=2e-2) diff --git a/pyproject.toml b/pyproject.toml index 510c910..5efc0e6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,6 +19,7 @@ dependencies = [ "pymupdf", "joblib", "ezdxf", + "scipy", ] dynamic = ["version"]