WIP: Added initial image reprojection functionality

[astrofrog]

This was discussed during SciPy as a possible feature to add to astropy. This routine would act as a high-level wrapper for interpolating reprojection, but also a wrapper to the drizzle algorithm for flux-conserving reprojection, and even a Montage-like algorithm. This is of course work in progress, but feedback about this idea is already welcome (will post to the list once it's a little more functional).

[astrofrog]

This includes the changes from #2730 which should be removed if this is ever merged.

[astrofrog]

@mdboom - at some point, could you send me information about how to get the core drizzle c code? As you can see here I've already started working on the steps to compute the pixel values to pass to drizzle, but I'm not sure exactly what inputs are required?

[Cadair]

I am trying to give this a go to reproject from HPLN/HPLT-TAN to HPLN/HPLT-CAR using this code:

https://gist.github.com/Cadair/f92481222021d5d258b1

It doesn't seem to actually be reprojecting anything, as the output is still circular :p

[mdboom]

The core drizzle C code is exposed to the outside world in a way that is very specific to how drizzlepac wants to use it at the moment. It will take some time to pick things apart. But if you want to see it:

https://svn.stsci.edu/svn/ssb/stsci_python/drizzlepac/trunk

[astrofrog]

@mdboom - no problem, this was just an idea following a conversation with @perrygreenfield, but it is not urgent. Also, the function here could be modified to better deal with what the drizzle code expects, what is here currently was just based on my understanding of what was needed.

[astrofrog]

@keflavich @cdeil - want to try this out? :)

[cdeil]

@ellisowen and I have cube data ... I think this handles only 2D, but @keflavich's FITS_tools.regrid_cube seems to work nicely for us.

In general I'm a huge fan of starting astropy.image in the core package as you do here. But what's the scope? Would functionality that is now e.g. in python-reprojection, FITS_tools, gammapy.image or ... be welcome? Is using scikit-image like we decided for photutils OK?

In my experience pull requests for Astropy take very much discussion, because we want to "get it right" before users try them out. Maybe using one of the existing repos (FITS_tools or python-reprojection) for ~ one months and then letting users try it out for ~ one months and then proposing the astropy.image sub-package would be more effective for developers and users?

[astrofrog]

This only handles image data for now, but indeed we can add 3D. Regarding the scope, there isn't one defined yet, this was just a bit of fun at the scipy sprint. What I really wanted here was reprojection tools, and since astropy.image precludes 3D for example, and astropy.image is quite broad, I'd actually like to change my suggestion to be to add astropy.reprojection. The aim here was really just to start having a place to centralize different reprojection algorithms, and we have a separate package for e.g. convolution. I think astropy.image has too broad a scope (and at the same time too narrow, since it excludes cubes).

The plan was that once something works here, I would post a message to the list to propose the idea, and we can then try and work from there, but I thought it would be nice to have some functionality examples to start discussing since previous discussions on the topic of image analysis did not lead far in the past.

[astrofrog]

Actually the more I think about it the more I think astropy.reprojection is the best way to go, otherwise it's waaaay too broad!

[astrofrog]

So if we do that, the scope is the following: 2- and 3-d data reprojection with:

• interpolation
• montage-like algorithm (adapted from python-reprojection) - 2d only
• drizzle - 2d only

@cdeil - would you be more comfortable with that scope?

[cdeil]

It would be nice if reprojection worked for n-d data (i.e. it operates on the two spatial axes, but there can be any number of other axes that are treated as independent planes that are reprojected one by one).

And I think the downsample and upsample functions from photutils and the block_reduce_hdu function from gammapy could fit there, too.

Personally I'd much prefer to throw this together in the separate python-reprojection package, make a 0.1 release in a few weeks, let users try it out and give feedback and then propose to move it to astropy.reprojection (or some similar name).

[keflavich]

I also like pushing this into python-reprojection, and separately getting a suite of image manipulation tools (downsample, upsample, etc.) into some astropy repository. Is astropy.image_tools too generic? Maybe image_utils? or array_utils?

[astrofrog]

@cdeil @keflavich - I guess my thinking was that if you just consider the problem of reprojection, it's a pretty straightforward problem API wise at the core level. You have an array, and two WCSs, and maybe a couple of options. So there's no so much the issue of trying to try out different APIs, etc.

But point taken, we can develop faster in python-reprojection and then propose to merge it in to the core. However, I'd like to make this a much faster process than for other packages because I think reprojection would be a killer feature for astropy 1.0. I'll revive development of python-reprojection and I think we could aim for an implementation that includes interpolation and the Montage algorithm in a timescale of a couple of weeks. Then that gives a few months of testing and we can still merge it in to the core package in time for 1.0 (if there are no issues of course).

Does that sound reasonable?

[keflavich]

Yes.

[cdeil]

👍 to doing this in the python-reprojection repo in the coming weeks.

We could decide on a case-by-case basis what is in scope and whether this would become astropy.reprojection or astropy.image or astropy.array_utils later.
I kind of like the idea of an astropy.image that will grow useful image manipulation functions over time, and as scikit-image it doesn't have to be limited to 2-d data ...it would be OK to operate on 3-d or n-d data where appropriate, i.e. the term "image" in the name doesn't have to prevent us from handling "cubes".

[astrofrog]

Ok, sounds good. I have a student looking for fun projects, so I'm going to talk to him about working on python-reprojection later this week. I'll update that package to the latest package-template and will move the stuff here into there.

[mdboom]

I've started looking at this in the context of fitting drizzle into it. I think there's one important API change we need. The resampling routine should return both the resampled image and a weighted mask showing the amount of contribution to each pixel (it's basically an alpha mask) -- it will end up being 1 most places in the image, and 0 outside, but some value in between along the edges. We probably don't want to do mosaicing in the first pass, but if we ever want to, you need that information in order to subsequently combine multiple images together later, and we'll kick ourselves by not outputting that information from all of the projections at this point.

@cdeil: Where is the python-reprojection repo? Can't find it.

[keflavich]

@mdboom - https://github.com/astrofrog/reproject/
The mask return is a good idea, but it belongs in a low-level method, since not all reprojection routines have weights (i.e., the interpolative ones do not)

[mdboom]

@keflavich: I'm not sure I follow. Any reprojection routine, regardless of method, is going to have a mask. Unless the reprojection fills the output image exactly, in which case the mask could be a scalar value of 1 -- but I would argue it should still return that so that images can be combined in a uniform way.

[keflavich]

@mdboom - point taken, I think you're right.

[astrofrog]

@mdboom - I agree, that sounds like a good idea!

[astrofrog]

Closing this since it is being developed in https://github.com/astrofrog/reproject/ - we will open a new PR once it is ready to put in the core package (if everyone agrees it belongs there).