Add Gaussian kernel, achieve flux conservation in adaptive algorithm

[svank]

This PR applies on top of PR #275, and wraps up my changes to the adaptive algorithm in reproject_adaptive.

This fixes a few bugs, exposes the flux-conserving mode already present in deforest.pyx, and adds a Gaussian kernel which offers much better photometry and anti-aliasing than the default Hann kernel. The bug fixes and Gaussian kernel together achieve very robust photometry---or flux conservation, when that mode is enabled (easily within 1% in my tests).

Summary of each commit:

• "Remove incorrect normalization in Hann window" and "Reverse sign of Jacobian values" both fix items that were incorrect but didn't cause any problems, as described in the commit messages. But fixing or removing things that are incorrect seems like good insurance against future bugs.
• "Use padded singular values..." fixes a bug, in that the padded singular values were not being used as called for by the DeForest (2004) description of the algorithm. Here is a notebook discussing this in more (too much?) detail, demonstrating that the fix does improve the outputs of the algorithm, and showing that the effect in the test cases (and hopefully most real-world cases) is minimal. (There are a few embedded videos that don't show up in the Gist interface, so download the notebook to view those.)
• "Expose the conserve_flux flag..." exposes in the top-level API an option that was already implemented in deforest.pyx, which uses the Jacobian to re-scale output pixel values in a flux-conserving way, as described in the DeForest (2004) paper.
• "Remove extra interpolation step..." fixes a problem in which this implementation was actually doing too much work, which turns out to have been detrimental. Again, here's a notebook discussing this in detail, demonstrating the fix, and showing that once again the test cases didn't suffer too much.
• "Add Gaussian filter kernel..." adds a Gaussian kernel as an alternative to the Hann kernel. This kernel trades some blurring of the output image for significantly-improved anti-aliasing and better flux-conserving behavior. (A later commit updates the code demo in docs/celestial.rst to show both kernels.)
• "Expand adaptive resampling tests" uses the conserve_flux mode to run tests over a grid of affine transformations, where the fact that flux should be conserved makes it really easy to quantify if the output is correct. The grid is 7D and so has to be pretty coarse to not slow down the testing too much, but all outputs show flux conservation to within 0.4%.
• "Documentation updates" does exactly that
• (EDIT: this commit has been moved to Changed default filter kernel and boundary mode in reproject_adaptive, and removed order argument. #291 ) "Make Gaussian kernel default..." uses the Gaussian kernel for reprojections with this algorithm, with the Hann left as an option. This is the most significant change in this PR series, as anyone using reproject_adaptive in their code will see very different results after updating to a version of reproject with this change. Their output images will be a bit blurrier, but better anti-aliased. (To be sure, every output image will also change due to the bug fixes, though those changes will hopefully be much smaller.) I'll argue that making this kernel the default is the right choice, though, because it strengthens the anti-aliasing quite a bit, and that anti-aliasing is what sets this algorithm apart from the interpolating and exact algorithms. Using the Gaussian by default offers the best anti-aliasing by default when a user comes to this function for anti-aliasing. I'm happy to discuss this point further and offer plots and examples.

[Cadair]

I'm wondering if separating things into three PRs that build on each other was the right strategy after all...

It was! You will have to update them as we merge them in order, but that's fine. Thanks so much for putting so much effort into making them easier to review!

[astrofrog]

@svank - feel free to rebase this now that #275 is merged! I can then try and review this :)

[svank]

@astrofrog Done. Thanks for reviewing!

[Cadair]

I think I just caused a conflict on this by merging #283 sorry

[svank]

@Cadair It was a pretty easy merge conflict to fix, no worries. I've rebased this and the other PR, and I think the CI failures are unrelated.

[Cadair]

Thanks, I shall try and find some time this week to address the CI issues and review your remaining PRs.

[Cadair]

"Make Gaussian kernel default..." uses the Gaussian kernel for reprojections with this algorithm, with the Hann left as an option. This is the most significant change in this PR series, as anyone using reproject_adaptive in their code will see very different results after updating to a version of reproject with this change.

This change does worry me. I am not really sure on what @astrofrog's policy on this kind of thing has been historically. I wonder if we should instead of just changing it in the next release issue a deprecation warning saying it will be changed in the following release? Making changes to things which affect peoples outputs without some kind of warning worries me.

[Cadair]

@svank hope you don't mind but I rebased this for you. Mainly so the CI runs properly and so that we get a doc build preview 🙂

[codecov]

Codecov Report

Merging #276 (cf40c3c) into main (3cfc7a8) will decrease coverage by 0.02%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main     #276      +/-   ##
==========================================
- Coverage   94.32%   94.30%   -0.03%     
==========================================
  Files          23       23              
  Lines         723      720       -3     
==========================================
- Hits          682      679       -3     
  Misses         41       41              

Impacted Files	Coverage Δ
reproject/adaptive/core.py	93.10% <100.00%> (ø)
reproject/adaptive/high_level.py	100.00% <100.00%> (ø)

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[astrofrog]

Regarding the change to the default - maybe remove that from this PR and we can discuss in a separate issue/PR and see what makes sense? I think having a deprecation warning would be sensible and then changing the default one major release later.

[svank]

Yeah, an abrupt change might not be the best. I'm on travel right now, so I'll separate out that one commit as soon as I can (no later than Monday). I think the other remaining PR had another default-behavior change around image boundaries that might need to likewise be pulled out.

[svank]

I've removed the defaults change from this PR and will submit it in another PR

[Cadair]

I am not going to pretend I have gone through all the cython changes and understand them, but I think this looks good.

The only other worry I have is the removal of the order keyword.

[Cadair]

The removal of this order keyword is also a breaking change. What was it being used for before this PR? Can we deprecate it?

[svank]

That's an inconvenient but good point. The order argument was controlling a linear interpolation step being applied before the weighted averaging that's at the core of this algorithm. The linear step isn't present in the paper describing this algorithm, and it turns out to have subtle effects that reduce the quality of the reprojection, as described and demonstrated here. It's an extra step with no upside that I'm aware of, so removing it completely makes sense to me.

One way to make it a non-breaking change might be to keep the order argument and just ignore it (and also raise a deprecation warning?). Outputs will also change slightly because of the removal of the linear interpolation step (not nearly as much as with a change to using the Gaussian kernel by default), but I'd put those changes in the "bug-fix" category. This PR and #274 both have other commits fixing bugs that will change outputs slightly.

[Cadair]

I would be happy with the kwarg still being there and just being ignored.

[svank]

Sounds good. I assume ignoring + raising a DeprecationWarning is the right path, so we can excise it completely later on. I'll wait to see if @astrofrog has thoughts or other review comments, though.

[Cadair]

I think raise a deprecation warning if order is anything other than the default value, as an approximation of user changed this, and then we can pull it in the subsequent release.

[svank]

I used one of the astropy decorators, which seems to trigger if any input value is explicitly provided:

In [10]: reproject.reproject_adaptive((data, w), w, data.shape, order=None)
WARNING: AstropyDeprecationWarning: "order" was deprecated in version 0.9 and will be removed in a future version.  [warnings]
Out[10]: 

[Cadair]

hey @svank sorry for leaving these to languish for a while. Is there anything I can do to help them along?

[svank]

@Cadair Thanks for checking. I just added a commit to address the deprecation of the 'order' argument we discussed earlier. Not sure if you're otherwise ready to merge this or if we're waiting for further review.

[astrofrog]

Looks good to me now with the deprecated order argument - thanks!