LM-SAL / aiapy

Python library for AIA data analysis
https://aiapy.readthedocs.io/en/stable/
BSD 3-Clause "New" or "Revised" License
7 stars 3 forks source link

Use reproject for image registration #317

Open wtbarnes opened 1 year ago

wtbarnes commented 1 year ago

Fixes #112

This is a first pass at using the reproject package to perform image registration by constructing a WCS for the aligned, translated, and scaled image and reprojecting the input map into that new WCS. This also now works for submaps. Note that this does not add a dependency as reproject is an optional dependency of sunpy that we already pick up by requiring sunpy[map].

I've marked this as a draft because it still needs quite a bit of cleaning up and I'd like several pairs of eyes on it before we even think about merging (if we do at all).

There is also the issue of performance. This is markedly slower for full-disk images by an OOM.

In [22]: %%timeit
    ...: m_l15_cutout = aiapy.calibrate.prep.register_reproject(m_l1_cutout)
1.6 s ± 15.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [23]: %%timeit
    ...: m_l15 = aiapy.calibrate.prep.register_reproject(m_l1)
33.5 s ± 749 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [24]: %%timeit
    ...: m_l15 = aiapy.calibrate.prep.register(m_l1)
WARNING: SunpyUserWarning: Integer input data has been cast to float64. [sunpy.image.transform]
3.34 s ± 10.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

However, this version does have the advantage that prep on submaps is possible. Additionally, reproject offers several different algorithms for performing the actual reprojection and [there are proposed performance improvements.](). There is also progress being made on parallelizing the map_coordinates function with Dask which is the underlying scipy function called by reproject_interp. There's also a cupy version of map_coordinates.

All my comments above are assuming that the performance bottleneck is the map_coordinates function, i.e. the actual interpolation of the array onto the new coordinates. While this is almost certainly true, it should be confirmed with some profiling.

Some TODOs: