Distortion analyses in NISP simulations

Author: Yannick Copin y.copin@ipnl.in2p3.fr

%matplotlib inline

# try:
#     import mpld3  # Interactive plots
#     mpld3.enable_notebook()
# except ImportError:
#     pass

try:
    import seaborn
    seaborn.set_style('ticks', {"axes.grid": True})
except ImportError:
    pass

import warnings
warnings.filterwarnings("ignore")

import numpy as N
from matplotlib import pyplot as P
from spectrogrism import spectrogrism as S
from spectrogrism import nisp
from spectrogrism import distortion as D

datapath = "../spectrogrism/data/"
simulations = S.Configuration([
    ("name", "Zemax"),
    (1, datapath + "run_190315.dat"),            # 1st-order dispersed simulation
    (0, datapath + "run_011115_conf2_o0.dat"),   # 0th-order dispersed simulation
    (2, datapath + "run_161115_conf2_o2.dat"),   # 2nd-order dispersed simulation
    ('J', datapath + "run_071215_conf6_J.dat"),  # J-band undispersed simulation
])
zmx = nisp.ZemaxPositions(simulations)
print(zmx)

Simulations 'Zemax': 4 mode(s)
  Order #1: ../spectrogrism/data/run_190315.dat
  Order #0: ../spectrogrism/data/run_011115_conf2_o0.dat
  Order #2: ../spectrogrism/data/run_161115_conf2_o2.dat
  Band   J: ../spectrogrism/data/run_071215_conf6_J.dat
  Wavelengths: 13 steps from 1.20 to 1.80 µm
  Coords: 289 sources

Photometric modes

Band J

xda = zmx.xda.sel(mode='J')
xda.to_dataset(name='zcoord').info()

xarray.Dataset {
dimensions:
        coordinate = 289 ;
        wavelength = 13 ;

variables:
        <U1 mode() ;
        float64 wavelength(wavelength) ;
        complex128 coordinate(coordinate) ;
        complex128 zcoord(wavelength, coordinate) ;

// global attributes:
}

Restructure mean position (averaged over wavelengths) into an ad-hoc structured grid:

xy = xda.mean(axis=0).values  # Mean (complex) positions (averaged over wavelengths) (289,)
xy = xy.reshape(-1, 17)      # (17, 17)
grid = D.StructuredGrid(xy)
print(grid)

Structured grid: 17 × 17 = 289 positions, y-x-

Estimate effective undistorted grid properties (step, rotation, minimal distortion position) from triangulation analysis (no optical modeling):

step, angle, offset, center = grid.estimate_parameters(fig=True)
print(f"Offset: ({offset.real:.3f}, {offset.imag:.3f}), step: {step:.3f}, angle: {N.rad2deg(angle):.1f}°")
print(f"Center of distortion: ({center.real:.3f}, {center.imag:.3f})")

Offset: (0.001, 0.000), step: 0.011, angle: 0.0°
Center of distortion: (-0.001, -0.012)

Construct effective undistorted grid used as reference:

refgrid = D.StructuredGrid.create(grid.nx, grid.ny, step=step, rotation=angle, offset=offset)
refgrid.reorder(grid.signature)
print(refgrid)

Structured grid: 17 × 17 = 289 positions, y-x-

Compute offset RMS between observed distorted (grid) and reference undistorted (refgrid) grids:

rms = grid.rms(refgrid.xy)
print("RMS = {:.3g} mm = {:.3g} px".format(rms / 1e-3, rms / 18e-6))

ax = grid.plot_offsets(refgrid, units=(18e-6, 'px'))
ax.figure.suptitle("Distortion wrt. ref. grid");

RMS = 0.313 mm = 17.4 px

It is already apparent that the distortion will not be well approximated by a radial expansion. This is an hint of off-axis optics.

dxy = grid.xy - refgrid.xy
r = N.abs(grid.xy - center).ravel()
dr = N.hypot(dxy.real, dxy.imag).ravel()

fig, ax = P.subplots(subplot_kw={'title': "Radial distortion wrt. reference grid",
                                 'xlabel': "r [mm]", 'ylabel': "dr [mm]"})
ax.plot(r / 1e-3, dr / 1e-3, '.');

Adjust the distortion with a simple r² radial expansion:

dist = D.GeometricDistortion(center, Kcoeffs=[0.], Pcoeffs=[])
print(dist)

Geometric distortion: center=(-0.000581568, -0.0117205), K-coeffs=[0.], P-coeffs=[]

minuit = refgrid.adjust_distortion(grid, dist, scale=True, offset=True)

#minuit.print_matrix()
adjrms = minuit.fval ** 0.5  # objfun is RMS**2
print("RMS = {:.4f} mm = {:.2f} px".format(adjrms / 1e-3, adjrms / 18e-6))

RMS = 0.1720 mm = 9.56 px

Reconstruct adjusted distorted grid from fit parameters:

adjdist = D.GeometricDistortion.from_kwargs(**minuit.values.to_dict())
adjgrid = D.StructuredGrid(adjdist.forward(refgrid.xy))

ax = grid.plot_offsets(adjgrid, units=(18e-6, 'px'))
ax.figure.suptitle("Distortion wrt. adjusted grid");

Given the fact that the off-axis distortion pattern will not be well approximated by a radial expansion, one can test a generic polynomial expansion:

from spectrogrism import polyfit

poldeg = (3, 3)
print("Total number of coeff.: {}".format((poldeg[0] + 1) * (poldeg[1] + 1)))
px = polyfit.fit_legendre2D(refgrid.x, refgrid.y, dxy.real.ravel(), deg=poldeg)
py = polyfit.fit_legendre2D(refgrid.x, refgrid.y, dxy.imag.ravel(), deg=poldeg)

Total number of coeff.: 16

polygrid = D.StructuredGrid(refgrid.xy +
                            px(refgrid.xy.real, refgrid.xy.imag) +
                            py(refgrid.xy.real, refgrid.xy.imag) * 1j)

ax = grid.plot_offsets(polygrid, units=(18e-6, 'px'))
ax.figure.suptitle("Distortion wrt. polynomial grid {}".format(poldeg));

Spectroscopic modes

Order 1

df = zmx.xda.sel(mode=1)
grid = D.StructuredGrid(df.mean(axis=0).values.reshape(-1, 17))

step, angle, offset, center = grid.estimate_parameters(fig=True)
print(u"Offset: ({0.real:.3f}, {0.imag:.3f}), scale: {1:.3f}, angle: {2:.1f}°"
     .format(offset, step, N.rad2deg(angle)))
print(u"Center of distortion: ({0.real:.3f}, {0.imag:.3f})".format(center))

Offset: (0.001, -0.003), scale: 0.011, angle: 0.0°
Center of distortion: (-0.002, -0.014)

refgrid = D.StructuredGrid.create(grid.nx, grid.ny, step=step, rotation=angle, offset=offset)
refgrid.reorder(grid.signature)  # Reference (undistorted) grid

rms = grid.rms(refgrid.xy)
print("RMS = {:.3g} mm = {:.3g} px".format(rms / 1e-3, rms / 18e-6))

ax = grid.plot_offsets(refgrid, units=(18e-6, 'px'))
ax.figure.suptitle("1st order mean position - Distortion wrt. ref. grid");

RMS = 0.289 mm = 16.1 px

dist = D.GeometricDistortion(center, Kcoeffs=[0.])
minuit = refgrid.adjust_distortion(grid, dist, scale=True, offset=True) #, print_level=2)

adjdist = D.GeometricDistortion.from_kwargs(**minuit.values.to_dict())
adjgrid = D.StructuredGrid(adjdist.forward(refgrid.xy))

ax = grid.plot_offsets(adjgrid, units=(18e-6, 'px'))
ax.figure.suptitle("1st order mean position - Distortion wrt. adj. grid");

Order 0

df0 = zmx.xda.sel(mode=0)
grid0 = D.StructuredGrid(df0.mean(axis=0).values.reshape(-1, 17))  # Observed grid

step, angle, offset, center = grid0.estimate_parameters()  # Initial parameters
refgrid0 = D.StructuredGrid.create(grid0.nx, grid0.ny, step=step, rotation=angle, offset=offset)
refgrid0.reorder(grid0.signature)  # Reference (undistorted) grid

# Adjust minimal mapping: scale + offset
dist0 = D.GeometricDistortion()
minuit0 = refgrid0.adjust_distortion(grid0, dist0, scale=True, offset=True) #, print_level=2)

adjdist0 = D.GeometricDistortion.from_kwargs(**minuit.values.to_dict())
adjgrid0 = D.StructuredGrid(adjdist0.forward(refgrid0.xy))  # Adjusted (distorted) grid

rms = grid0.rms(adjgrid0.xy)
print("RMS = {:.3g} mm = {:.3g} px".format(rms / 1e-3, rms / 18e-6))

ax = grid0.plot_offsets(adjgrid0, units=(18e-6, 'px'))
ax.figure.suptitle("0th order mean position - Distortion wrt. ref. grid");

RMS = 0.326 mm = 18.1 px

This page was generated from NISP_distortions.ipynb.