New Coalignment API

[Deus1704]

Fixes https://github.com/sunpy/sunkit-image/issues/83 Fixes https://github.com/sunpy/sunkit-image/issues/103

Following our discussions, in this pr I've laid the skeleton for the new coalignment interface. A basic use case is demonstrated, where the match_template method is registered by default. For details on the registration process, please refer to the code.

aia_map1 = sunpy.map.Map(sunpy.data.sample.AIA_193_CUTOUT01_IMAGE)
aia_map2 = sunpy.map.Map(sunpy.data.sample.AIA_193_CUTOUT03_IMAGE)
### Creating a template from aia_map1
bottom_left = SkyCoord(600 * u.arcsec, -500 * u.arcsec, frame=aia_map1.coordinate_frame)
top_right = SkyCoord(800 * u.arcsec, -200 * u.arcsec, frame=aia_map1.coordinate_frame)
submap = aia_map1.submap(bottom_left, top_right=top_right)

coaligned_map = coalignment(aia_map2, submap,"match_template")

Tasks:

• [x] Handle clipping
• [x] Deprecate existing API
• [x] Add tests
• [x] Add tests for checking the newly created functions
• [x] Add gallery examples for the new api
• [x] Write an exhaustive changelog
• [x] Decide on a better name for this method than "match_template"

[nabobalis]

I have a broader comment.

Can we not touch the original code but create a new coalignment folder and have separate files for the interface and each method?

This would make it easier to prototype and review for others.

[Deus1704]

I have a broader comment.

Can we not touch the original code but create a new coalignment folder and have separate files for the interface and each method?

This would make it easier to prototype and review for others.

That works. Let me get that done quickly.

[Deus1704]

Should we consider returning only the shifts(pixel/arcsec) required for aligning? Like an additional argument apply_shifts = False that could be passed to the coalignment_interface and return only the shifts. Would that be useful?

[nabobalis]

On the surface, I would rather return the shifted maps.

It might be useful information to have tho but maybe I would want that to be a seperate function.

[kreardon]

I think optionally outputting the calculated shifts would be very useful, since it is already being calculated in this process.

[kreardon]

I also think it would make more sense to have the inputs be ordered as template_map, input_map, method And maybe template_map isn't a great name, since there may be other, non template_matching coalignment methods that could be used? I would call the two input maps something like reference_map and target_map

[Deus1704]

I also think it would make more sense to have the inputs be ordered as template_map, input_map, method And maybe template_map isn't a great name, since there may be other, non template_matching coalignment methods that could be used? I would call the two input maps something like reference_map and target_map

Yeah, agreed. This naming and order would be better.

[nabobalis]

I think optionally outputting the calculated shifts would be very useful, since it is already being calculated in this process.

Would you want the map/cube and the shifts or would you just want the shifts alone?

[kreardon]

I am a bit confused by what sort of "coalignment" this is doing. In the above example, the submap reference is [500, 335] pixels, while the target aia_map2 is [1345, 1217] pixels. The coaligned_map is [956, 1047] pixels and, as can be seen in example images, isn't really coaligned with the reference. There are some x,y ordering and positive/negative sign errors in the code, but even fixing those, the registration process doesn't quite make sense.

The method employed is to use match_template to get the shift and then apply that by "clipping" the array to generate the coaligned map. But the clipping, as applied, always only works on one side of each axis (hence the weird size). Basically it seems like it is only matching the lower left or the upper right corner of the reference array, which isn't quite the expected result.

This all raises some fundamental questions:

1. if the input arrays are different sizes (with the reference smaller), should the output map be the size of the reference or the target array? (my answer? for most use cases, I think output should be size of the target array suitably shifted and padded)
2. the most common usage will likely be two arrays/maps of the same size, and in this case match_template does not work as a method (reference must be smaller). Shouldn't it work for equivalent sized arrays?
3. the algorithm calculates the subpixel offsets, but only applies integer shifts - should interpolating to apply the sub-pixel offsets be an option (or default)?
4. the function takes two maps as inputs, but doesn't use any of the pixel scale or alignment information to generate two arrays that should be directly compared for alignment. Instead of assuming that the plate scales are the same for both maps, shouldn't the arrays be remapped onto a common grid before the cross-correlation calculation?
5. perhaps related to question (1), depending on how it is implemented, should it be possible to perform the shift calculation on one portion of a map/array and then apply those shift to the full map/array? (my answer? yes, this is a common usage scenario)

[nabobalis]

I am a bit confused by what sort of "coalignment" this is doing.

Ideally when we have the interface in place, we can either re-write this method to be useful or we remove it and implement something new.

I am not too worried about the fact that this really old code from sunpy isn't actually good.

Are you aware of any papers about co-alignment methods we could look at implementing here? We want to try and have only methods in this library that are backed by publications if possible.

[kreardon]

Thanks, Nabil. I realize the actual alignment code is just one implementation. But I think my questions were mostly non-implementation-specific, about what inputs and outputs the users would like to employ with such an interface (e.g. array sizes, non-scale-matched) and the desired alignment behaviour (e.g. sub-pixel interpolation, sub-array alignment target).

[nabobalis]

Thanks, Nabil. I realize the actual alignment code is just one implementation. But I think my questions were mostly non-implementation-specific, about what inputs and outputs the users would like to employ with such an interface (e.g. array sizes, non-scale-matched) and the desired alignment behaviour (e.g. sub-pixel interpolation, sub-array alignment target).

Ah I see, I apologize, I thought it was angled at the really old/bad implementation.

about what inputs and outputs the users would like to employ with such an interface

I think this is a topic we want to get as many user feedback as possible. If you happen to know others who have thoughts, I would love to hear it as well.

(e.g. array sizes, non-scale-matched) and the desired alignment behaviour (e.g. sub-pixel interpolation, sub-array alignment target)

Naively, I would assume that is method dependent? Would users tweak these kind of behaviours as keywords typically or would that be method specific keywords?

I will say this as a person who hasn't had to coalign data before and I am not familiar with current methods.

[nabobalis]

So some basic tests for the match template would be good, I would hope the original ones would work.

On a different note. Do we want the ability for the user to pass keywords arguments to the method? Depending on the algorithm, I can imagine the need for a range of method specific keywords but nothing we could know ahead of time for the main level.

[kreardon]

I think this is a topic we want to get as many user feedback as possible. If you happen to know others who have thoughts, I would love to hear it as well.

Yes, I think a map-aware co-alignment process could be useful for many users. The added layer of complexity is automatically dealing with two maps which don't share the same image scaling (as well as higher-order differences like rotation?). One could think of aligning an IRIS slitjaw image and an AIA 1600 image as a useful example. Do we want this API to apply the appropriate scaling before calculating the coalignment offsets? Or should the API only deal with data arrays that have already been mapped to a common grid through some other procedure?

Things get more complicated when thinking about the rarer (corner) case of aligning SolO/EUI and SDO/AIA images taken from different viewing angles. Sometimes these can be appropriately aligned if the two maps are on a common physical scale (km at solar surface/pixel), other times it might make more sense to reproject the images to a common viewing geometry.

So trying to sketch out the inputs and outputs of such an API, here is what I came up with (including some question marks for options that may depend on implementation):

Inputs:

1. Reference object - data array (including coordinate information?)
2. Unaligned object - data array (including coordinate information?)
   • Unaligned object effective coverage could be bigger, smaller, or same size as Reference object?
3. Alignment Region - (portion of Reference object to use for offset/alignment calculation
4. Alignment Method - predefined list
5. Aligned Object Definition - values defining output map size and coordinates

Alignment Process:

1. Put data arrays on common coordinate grid
2. Calculate pixel coordinate correspondence between data arrays using defined Alignment Region and Method
3. Apply integer pixel shifts to Unaligned object
4. Pad Unaligned object as requested in Aligned Object Definition
5. Apply sub-pixel shifts as requested in Aligned Object Definition
6. Map data to coordinate grid requested in Aligned Object Definition

Outputs:

1. Aligned object (including coordinate information)
2. Calculated image offsets (as pixel or coordinate values?)

[nabobalis]

Yes, I think a map-aware co-alignment process could be useful for many users. The added layer of complexity is automatically dealing with two maps which don't share the same image scaling (as well as higher-order differences like rotation?). One could think of aligning an IRIS slitjaw image and an AIA 1600 image as a useful example. Do we want this API to apply the appropriate scaling before calculating the coalignment offsets? Or should the API only deal with data arrays that have already been mapped to a common grid through some other procedure?

That is a good question, my worry is that if we talk about these kind of topics we get closer to the idea of reinventing reproject.

If the API (or part of the API) instead returned (@Cadair's idea) the shifts required in the WCS only. Then a user can decide how to apply it themselves and take account of the different scaling or any other projection effects first?

Things get more complicated when thinking about the rarer (corner) case of aligning SolO/EUI and SDO/AIA images taken from different viewing angles. Sometimes these can be appropriately aligned if the two maps are on a common physical scale (km at solar surface/pixel), other times it might make more sense to reproject the images to a common viewing geometry.

In this case would we make the user do the reprojection first then do the coalignment? That to me seems like a option but maybe I am barking up the wrong tree.

[wtbarnes]

I'm going to try to summarize the discussion of the direction of this API that has occurred at the last community meeting on 7/10 as well as the last GSoC meeting between myself, @nabobalis and @Deus1704.

We have come to the consensus that this new coalignment API should modify the metadata rather than the actual data array. In the current implementation (i.e. prior to this PR), the shifts were being computed based on the cross-correlation between the reference and target array and then the target array was being shifted. Instead, we should be adjusting the metadata based on these computed shifts (e.g. by adjusting the reference pixel and/or PCij matrix).

The logic here is that coalignment, in this context, only serves to correct pointing metadata in cases where this metadata is thought to be incorrect. If you also want the reference and target images to lie on the same pixel grid, you should then use reproject to perform this transformation.

In particular, any given coalignment method should return an affine transform (scale, rotation, shift) and then our general coalignment framework should handle updating the map metadata appropriately given these returned parameters. This means that any coalignment method need only take in the reference and target arrays and then return the affine parameters such that they need no knowledge of the metadata/WCS.

When it comes to updating the metadata, we should keep In mind that at some point in the future, we may want to update non-FITS-WCS metadata. To this end, we should organize this part of the coalignment function to make it easy to swap out a different kind of metadata update scheme in the future. For now, we need only provide this functionality for FITS-WCS metadata.

[Cadair]

To be extra clear:

• the functions should return any of a shift, scale or affine transform matrix required to align the two images. These parameters would need to be combined with the original metadata (I.e. pc @ pc)
• the pluggable bit of the API returning these parameters should do that in a agnostic way to the WCS framework. So some kind of data class or namedtuple or similar to hold the parameters.

[Cadair]

I can also see some of the use for implementing custom regridding routines based on the more limited transforms we are thinking about here (I.e. a simple affine transform) could be implemented in a much faster manner than pure reproject. But that's definitely a future long-term thing if someone has enough motivation.

[wtbarnes]

Two suggestions for figure tests to add:

• Create a cutout from a full-disk AIA image, make the metadata slightly wrong, coalign with the original full-disk image
• Create a fake map with some "interesting" feature (e.g. a Gaussian), create a smaller map with a similar source but at a slightly different location, coalign the two

[Deus1704]

@nabobalis We'll need to increase the py-310 minimum dependency of Astropy to 6.1.0.

[nabobalis]

@nabobalis We'll need to increase the py-310 minimum dependency of Astropy to 6.1.0.

How come?

[Deus1704]

How come?

The origin_mismatch keyword in the separation method was introduced in the 6.1.0 version. So if we are using it by default in our interface then we should have astropy>=6.1 or it will always error stating,

TypeError: SkyCoord.separation() got an unexpected keyword argument 'origin_mismatch'

[nabobalis]

Ah I see, I don't think we should increase the min astropy version to that until sunpy does as well.

We will have to check what version of astropy is installed, if its 6.1.0, add the keyword otherwise, I think we should catch that warning and ignore it via the warnings module.

[wtbarnes]

I just pushed a half-baked, in-progress attempt at enforcing consistency in ordering and naming of the arrays. There are still a few fixes to be made so hold off on any reviews, tweaks, etc.