{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Fitting the ESIS-I distortion\n", "\n", "This notebook fits the alignment/distortion parameters of the ESIS-I optical\n", "model to a Level-1 flight image, using a synthetic scene built from AIA\n", "imagery captured during the flight.\n", "\n", "Inverting an ESIS image requires a model whose **distortion mapping** — where\n", "each point on the Sun lands on each of the four detectors — matches the real\n", "instrument. Metrology alone does not pin the alignment down precisely enough,\n", "so we recover it by adjusting the model's distortion parameters until the\n", "modeled image reproduces a real Level-1 frame.\n", "\n", "It starts from the per-channel best fit recorded in\n", "`esis.flights.f1.optics.distortion_fit()` and performs a **coupled polish** of\n", "all four channels at once: parameters that correspond to a single physical\n", "part (the field stop, the primary mirror, and the pointing of the payload)\n", "are shared between the channels, while the per-grating parameters remain\n", "independent.\n", "\n", "Requirements:\n", "\n", "* The `JSOC_EMAIL` environment variable must be set to an email address\n", " [registered with JSOC](http://jsoc.stanford.edu/ajax/register_email.html)\n", " (only needed the first time; the AIA scene is cached afterwards).\n", "* Set `quick = False` below for a production fit (hours); the default\n", " settings run a minutes-long demonstration fit." ] }, { "cell_type": "code", "execution_count": null, "id": "1", "metadata": {}, "outputs": [], "source": [ "import dataclasses\n", "import pathlib\n", "from datetime import datetime\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import astropy.units as u\n", "import named_arrays as na\n", "import optika\n", "import esis" ] }, { "cell_type": "code", "execution_count": null, "id": "2", "metadata": {}, "outputs": [], "source": [ "# fast demonstration settings; set False for a production fit\n", "quick = True\n", "\n", "num_scene = 201 if quick else 801" ] }, { "cell_type": "markdown", "id": "3", "metadata": {}, "source": [ "## Data\n", "\n", "Load the Level-1 frame that the original distortion fit was optimized\n", "against, and the AIA frame captured closest to it." ] }, { "cell_type": "code", "execution_count": null, "id": "4", "metadata": {}, "outputs": [], "source": [ "time_index = 15\n", "\n", "l1 = esis.flights.f1.data.level_1()[dict(time=time_index)]\n", "time_esis = l1.inputs.time_start[dict(channel=0)].ndarray\n", "time_esis" ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [], "source": [ "scene = esis.flights.f1.data.synth.scene_aia()\n", "scene = scene[np.argmin(np.abs((scene.inputs.time - time_esis).mean(\"wavelength\")))]\n", "scene.outputs.shape" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "Downsample the scene onto a regular grid to make the raytrace cheaper.\n", "The scene radiance is regridded conservatively so that the total flux is\n", "preserved." ] }, { "cell_type": "code", "execution_count": null, "id": "7", "metadata": {}, "outputs": [], "source": [ "inputs_scene = na.TemporalSpectralPositionalVectorArray(\n", " time=scene.inputs.time,\n", " wavelength=scene.inputs.wavelength,\n", " position=na.Cartesian2dVectorArray(\n", " x=na.ScalarLinearSpace(\n", " scene.inputs.position.x.min(),\n", " scene.inputs.position.x.max(),\n", " axis=\"detector_x\",\n", " num=num_scene,\n", " ),\n", " y=na.ScalarLinearSpace(\n", " scene.inputs.position.y.min(),\n", " scene.inputs.position.y.max(),\n", " axis=\"detector_y\",\n", " num=num_scene,\n", " ),\n", " ),\n", ")\n", "scene_fit = scene(inputs_scene, axis=(\"detector_x\", \"detector_y\"), method=\"conservative\")\n", "scene_fit.outputs.shape" ] }, { "cell_type": "code", "execution_count": null, "id": "8", "metadata": {}, "outputs": [], "source": [ "spectral_lines = [\"O V\", \"Mg X\", \"He I\"]\n", "\n", "fig, ax = plt.subplots(ncols=3, figsize=(12, 4), constrained_layout=True)\n", "for i in range(3):\n", " index = dict(wavelength=i, velocity=0)\n", " ax[i].imshow(scene_fit.outputs[index].ndarray.value.T, origin=\"lower\")\n", " ax[i].set_title(spectral_lines[i]);" ] }, { "cell_type": "markdown", "id": "9", "metadata": {}, "source": [ "## Model\n", "\n", "Start from the hard-coded per-channel best fit and idealize the materials,\n", "so that the modeled counts are not modulated by the wavelength-dependent\n", "coating and filter efficiencies (the correlation merit is insensitive to\n", "per-channel scale factors anyway)." ] }, { "cell_type": "code", "execution_count": null, "id": "10", "metadata": {}, "outputs": [], "source": [ "model = esis.flights.f1.optics.distortion_fit(num_distribution=0)\n", "\n", "model.grating.material = optika.materials.Mirror()\n", "model.primary_mirror.material = optika.materials.Mirror()\n", "model.filter.material = None\n", "model.camera.sensor.material = optika.sensors.materials.IdealSensorMaterial()\n", "\n", "# a single pupil cell is sufficient for distortion mapping\n", "pupil = na.Cartesian2dVectorLinearSpace(\n", " start=-0.25,\n", " stop=0.25,\n", " axis=na.Cartesian2dVectorArray(\"pupil_x\", \"pupil_y\"),\n", " num=2,\n", ")" ] }, { "cell_type": "markdown", "id": "11", "metadata": {}, "source": [ "## Starting comparison\n", "\n", "Render the model image and compare it against the Level-1 data." ] }, { "cell_type": "code", "execution_count": null, "id": "12", "metadata": {}, "outputs": [], "source": [ "kwargs_image = dict(\n", " pupil=pupil,\n", " axis_wavelength=\"velocity\",\n", " axis_field=(\"detector_x\", \"detector_y\"),\n", " noise=False,\n", ")\n", "\n", "\n", "def render(instrument):\n", " image = instrument.system.image(scene_fit, **kwargs_image)\n", " axis_extra = tuple(set(image.outputs.shape) - set(na.shape(l1.outputs)))\n", " return image.outputs.sum(axis_extra)\n", "\n", "\n", "image_start = render(model)\n", "image_start.shape" ] }, { "cell_type": "code", "execution_count": null, "id": "13", "metadata": {}, "outputs": [], "source": [ "def compare(image_model, title):\n", " fig, ax = na.plt.subplots(\n", " figsize=(12, 10),\n", " constrained_layout=True,\n", " axis_rows=\"channel\",\n", " nrows=l1.shape[\"channel\"],\n", " axis_cols=\"column\",\n", " ncols=2,\n", " sharex=True,\n", " sharey=True,\n", " )\n", " fig.suptitle(title)\n", " data = (l1.outputs - l1.outputs.mean()) / l1.outputs.std()\n", " na.plt.pcolormesh(\n", " l1.inputs.pixel.x,\n", " l1.inputs.pixel.y,\n", " C=data.value,\n", " ax=ax[dict(column=0)],\n", " vmax=np.percentile(data.value, 99),\n", " )\n", " image_model = (image_model - image_model.mean()) / image_model.std()\n", " na.plt.pcolormesh(\n", " l1.inputs.pixel.x,\n", " l1.inputs.pixel.y,\n", " C=image_model.value,\n", " ax=ax[dict(column=1)],\n", " vmax=np.percentile(image_model.value, 99),\n", " )\n", " na.plt.text(\n", " x=0.01,\n", " y=0.97,\n", " s=l1.channel,\n", " transform=na.plt.transAxes(ax[dict(column=0)]),\n", " ax=ax[dict(column=0)],\n", " ha=\"left\",\n", " va=\"top\",\n", " color=\"white\",\n", " )\n", " ax.ndarray[0, 0].set_title(\"Level-1 data\")\n", " ax.ndarray[0, 1].set_title(\"model\")\n", "\n", "\n", "compare(image_start, \"starting model (hard-coded per-channel fit)\")" ] }, { "cell_type": "markdown", "id": "14", "metadata": {}, "source": [ "## Coupled parameters\n", "\n", "The hard-coded fit lets every parameter vary independently per channel, but\n", "several of them describe a single physical part, so their per-channel scatter\n", "is unphysical:\n", "\n", "* `roll_field_stop` — there is only one field stop,\n", "* `displacement_primary` — there is only one primary mirror,\n", "* `pitch`, `yaw`, `roll` — the payload has a single pointing.\n", "\n", "Replace those with their channel means to form the coupled starting point.\n", "The per-grating parameters (`yaw_grating`, `pitch_grating`, `roll_grating`,\n", "`spacing_rulings`) remain independent, since each channel has its own\n", "grating. This reduces the fit from 36 to 21 degrees of freedom." ] }, { "cell_type": "code", "execution_count": null, "id": "15", "metadata": {}, "outputs": [], "source": [ "parameters = esis.optics.DistortionParameters.from_instrument(model)\n", "\n", "parameters.roll_field_stop = parameters.roll_field_stop.mean()\n", "parameters.displacement_primary = parameters.displacement_primary.mean()\n", "parameters.pitch = parameters.pitch.mean()\n", "parameters.yaw = parameters.yaw.mean()\n", "parameters.roll = parameters.roll.mean()\n", "\n", "print(\"degrees of freedom:\", na.pack(parameters).shape[\"pack\"])\n", "parameters" ] }, { "cell_type": "markdown", "id": "16", "metadata": {}, "source": [ "## Polish bounds\n", "\n", "Since this is a polish around a known-good starting point, the bounds are\n", "tight absolute windows around the start, sized to comfortably contain the\n", "per-channel scatter of the original fit." ] }, { "cell_type": "code", "execution_count": null, "id": "17", "metadata": {}, "outputs": [], "source": [ "delta = esis.optics.DistortionParameters(\n", " yaw_grating=2 * u.arcmin,\n", " pitch_grating=2 * u.arcmin,\n", " roll_grating=0.25 * u.deg,\n", " roll_field_stop=1 * u.deg,\n", " spacing_rulings=5e-4 * u.um,\n", " displacement_primary=3 * u.mm,\n", " pitch=15 * u.arcsec,\n", " yaw=15 * u.arcsec,\n", " roll=0.5 * u.deg,\n", ")\n", "\n", "\n", "def shift(parameters, delta, sign):\n", " return esis.optics.DistortionParameters(**{\n", " field.name: getattr(parameters, field.name) + sign * getattr(delta, field.name)\n", " for field in dataclasses.fields(parameters)\n", " })\n", "\n", "\n", "bounds = (shift(parameters, delta, -1), shift(parameters, delta, +1))" ] }, { "cell_type": "markdown", "id": "18", "metadata": {}, "source": [ "## How the fit is scored\n", "\n", "Each trial parameter set is graded by `esis.optics.DistortionObjective`, which\n", "reproduces the scene through the perturbed instrument and scores the result as\n", "\n", "$$\\mathrm{merit} = -1000 \\cdot \\mathrm{correlation}(\\mathrm{model}, \\mathrm{data}) + \\mathrm{penalty}_\\mathrm{off\\text{-}target}.$$\n", "\n", "* **correlation** is the Pearson correlation between the standardized modeled\n", " and observed images. It rewards getting structure in the *right place* and is\n", " blind to overall brightness — which is why the idealized materials in the\n", " Model section are harmless, and why per-channel brightness differences are\n", " handled by averaging the correlation across channels (`axis_channel`).\n", "* **off-target penalty** is the squared mean ray position, which keeps the\n", " field from walking off the detector.\n", "\n", "Lower is better. The ray trace samples photons randomly, so the merit is\n", "**stochastic even with `noise=False`**; gradient-based optimizers would chase\n", "that noise, so the fit uses a derivative-free global optimizer\n", "(`scipy.optimize.differential_evolution`)." ] }, { "cell_type": "markdown", "id": "19", "metadata": {}, "source": [ "## Coupled polish\n", "\n", "Fit all four channels at once. With `axis_channel=\"channel\"` the correlation\n", "merit is computed independently per channel and averaged, so per-channel\n", "brightness differences in the data do not suppress it.\n", "\n", "Convergence is logged to a timestamped directory\n", "(`convergence_data.csv`, `full_output.log`, `convergence_plot.png`)." ] }, { "cell_type": "code", "execution_count": null, "id": "20", "metadata": {}, "outputs": [], "source": [ "directory = pathlib.Path(\n", " f\"distortion_polish_{datetime.now().strftime('%Y%m%d_%H%M%S')}\"\n", ")\n", "\n", "x0 = na.pack(parameters).ndarray\n", "lower = na.pack(bounds[0]).ndarray\n", "upper = na.pack(bounds[1]).ndarray\n", "\n", "if quick:\n", " rng = np.random.default_rng(42)\n", " init = lower + (upper - lower) * rng.uniform(size=(5, x0.size))\n", " init[0] = x0\n", " kwargs_optimizer = dict(init=init, maxiter=2, tol=1e6, polish=False)\n", "else:\n", " kwargs_optimizer = dict(\n", " popsize=10,\n", " maxiter=200,\n", " workers=-1,\n", " updating=\"deferred\",\n", " polish=True,\n", " )\n", "\n", "fitted = esis.optics.fit_distortion(\n", " instrument=model,\n", " scene=scene_fit,\n", " observation=l1.outputs.value,\n", " bounds=bounds,\n", " parameters=parameters,\n", " pupil=pupil,\n", " axis_wavelength=\"velocity\",\n", " axis_field=(\"detector_x\", \"detector_y\"),\n", " axis_channel=\"channel\",\n", " directory=directory,\n", " kwargs_optimizer=kwargs_optimizer,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "21", "metadata": {}, "outputs": [], "source": [ "fitted" ] }, { "cell_type": "markdown", "id": "22", "metadata": {}, "source": [ "## Result\n", "\n", "The `compare()` plots below put the Level-1 data next to the modeled image for\n", "each channel, and the final cell prints the merit at the coupled starting point\n", "versus after the polish (more negative is a tighter match).\n", "\n", "A useful sanity check is that collapsing the shared parameters\n", "(`roll_field_stop`, `displacement_primary`, and the `pitch`/`yaw`/`roll`\n", "pointing) from their per-channel values to a single channel mean barely changes\n", "the merit: the per-channel scatter in those parameters was not really earning\n", "its keep, which confirms that they are genuinely shared and that the 21-degree-\n", "of-freedom coupled model is the honest one to polish from. The demonstration\n", "settings (`quick=True`) only exercise the mechanics in a couple of iterations;\n", "the real improvement comes from the `quick=False` production run." ] }, { "cell_type": "code", "execution_count": null, "id": "23", "metadata": {}, "outputs": [], "source": [ "model_fitted = fitted.to_instrument(model)\n", "image_fitted = render(model_fitted)\n", "\n", "compare(image_fitted, \"coupled polish\")" ] }, { "cell_type": "code", "execution_count": null, "id": "24", "metadata": {}, "outputs": [], "source": [ "objective = esis.optics.DistortionObjective(\n", " instrument=model,\n", " parameters=parameters,\n", " scene=scene_fit,\n", " observation=l1.outputs.value,\n", " pupil=pupil,\n", " axis_wavelength=\"velocity\",\n", " axis_field=(\"detector_x\", \"detector_y\"),\n", " axis_channel=\"channel\",\n", ")\n", "\n", "print(\"merit at coupled start:\", objective(na.pack(parameters).ndarray))\n", "print(\"merit after polish: \", objective(na.pack(fitted).ndarray))" ] }, { "cell_type": "markdown", "id": "25", "metadata": {}, "source": [ "## Faster polish (merit scans)\n", "\n", "The differential-evolution polish above is robust but slow, because each\n", "evaluation of the merit is stochastic: the imaging model jitters every ray\n", "randomly within its scene cell. Freezing that jitter with a fixed `seed`\n", "makes the merit a deterministic function of the parameters, and convolving\n", "the modeled image with a one-pixel Gaussian point-spread function\n", "(`sigma_psf`) both models the blur of the real optics and makes the sparse\n", "raytraced image respond smoothly to sub-pixel parameter changes.\n", "\n", "Even so, the smoothed correlation surface is locally too flat for\n", "derivative-based optimizers: `scipy.optimize.least_squares` was tried and\n", "terminates after moving a tiny fraction of the distance to the peak.\n", "`esis.optics.fit_distortion_scan` instead reads the shape of the merit basin\n", "directly: each parameter is scanned along a coarse-to-fine schedule of\n", "offsets and moved to the peak of its correlation curve, refined with a\n", "parabola. With `axis_channel` given, every channel's peak is read from its\n", "own curve during a single pass of coherent scans, so one scan fits all four\n", "channels at once (here `coherent=True` instead applies one shared offset,\n", "matching the coupled single-payload pointing above). Strongly-coupled pairs\n", "can be scanned jointly by keying a round with a tuple of field names.\n", "\n", "This is the engine behind the committed reference and pointing fits; see\n", "`esis.flights.f1.optics.fit_distortion_reference` and\n", "`esis.flights.f1.optics.fit_distortion_pointing` to reproduce them." ] }, { "cell_type": "code", "execution_count": null, "id": "26", "metadata": {}, "outputs": [], "source": [ "grids = [\n", " dict(\n", " pitch=np.linspace(-10, 10, 21) * u.arcsec,\n", " yaw=np.linspace(-10, 10, 21) * u.arcsec,\n", " roll=np.linspace(-0.4, 0.4, 11) * u.deg,\n", " ),\n", " dict(\n", " pitch=np.linspace(-1, 1, 11) * u.arcsec,\n", " yaw=np.linspace(-1, 1, 11) * u.arcsec,\n", " roll=np.linspace(-0.05, 0.05, 11) * u.deg,\n", " ),\n", "]\n", "if quick:\n", " grids = [\n", " dict(\n", " pitch=np.linspace(-10, 10, 5) * u.arcsec,\n", " yaw=np.linspace(-10, 10, 5) * u.arcsec,\n", " ),\n", " ]\n", "\n", "fitted_scan = esis.optics.fit_distortion_scan(\n", " instrument=model,\n", " scene=scene_fit,\n", " observation=l1.outputs.value,\n", " grids=grids,\n", " parameters=parameters,\n", " pupil=pupil,\n", " axis_wavelength=\"velocity\",\n", " axis_field=(\"detector_x\", \"detector_y\"),\n", " axis_channel=\"channel\",\n", " coherent=True,\n", " sigma_psf=1.0,\n", " directory=pathlib.Path(\n", " f\"distortion_polish_scan_{datetime.now().strftime('%Y%m%d_%H%M%S')}\"\n", " ),\n", ")\n", "fitted_scan" ] }, { "cell_type": "code", "execution_count": null, "id": "27", "metadata": {}, "outputs": [], "source": [ "model_scan = fitted_scan.to_instrument(model)\n", "image_scan = render(model_scan)\n", "\n", "compare(image_scan, \"pointing polish (merit scans)\")" ] }, { "cell_type": "code", "execution_count": null, "id": "28", "metadata": {}, "outputs": [], "source": [ "print(\"merit, differential evolution:\", objective(na.pack(fitted).ndarray))\n", "print(\"merit, scans: \", objective(na.pack(fitted_scan).ndarray))" ] }, { "cell_type": "markdown", "id": "29", "metadata": {}, "source": [ "## The committed result\n", "\n", "The production fits made with this machinery are committed to the repository\n", "as plain-text data files, loaded by `esis.flights.f1.optics.distortion_fit`:\n", "\n", "* `_data/distortion_reference.ecsv` — the per-channel reference parameters,\n", "* `_data/distortion_pointing.ecsv` — the per-frame payload pointing offsets\n", " (available via `distortion_fit(axis_time=\"time\")`).\n", "\n", "Blinking the fitted model against every Level-1 frame of the flight shows the\n", "quality of the registration as the payload pointing drifts:\n", "\n", "![model and data blink comparison](../_static/blink.gif)" ] } ], "metadata": { "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.10" } }, "nbformat": 4, "nbformat_minor": 5 }