Add block_reduce and block_replicate functions

[larrybradley]

The downsample and upsample (Cython) functions weren't merged as part of the large imageutils PR (#3201) because they could be accomplished using ndarray tricks (which are a bit slower than the Cython-based functions). This PR implements downsample and upsample using the ndarray tricks, which are likely too obscure for most users. Both functions can handle 1D or 2D data, and downsample will trim the input array if necessary.

[mhvk]

Two quick comments:

• why the restriction to 1 or 2D? It should be possible to generalize quite easily, possibly with different block sizes in each dimension.
• especially for upsample, I found usemean a confusing term. In other contexts, I have seen conserve_flux? (which implies that sum(data) is conserved)

[hamogu]

I think this is more restricted that necessary. It basically offers two function that can be used: sum and mean. At least for downsampling there could be a parameter that defaults tofunction=np.sum. For the mean you would set function=np.mean, but I've used the same trick with the standard deviation in the past and it could also be used with min or max.
This would also make the implementation simpler - see inline comment.

[hamogu]

Instead of sum this could be a parameter that the use could set to any jumpy ufunc.

[astrofrog]

+1 to @hamogu's suggestion

[larrybradley]

Thanks for the suggestions. I started to implement them, but then I realized that downsample would essentially be a duplication of scikit-image's block_reduce:

http://scikit-image.org/docs/dev/api/skimage.measure.html#block-reduce

There is no exact equivalent of upsample in scikit-image (but it does have an interpolation-based resize function).

So, should we not add a rebining/downsampling function to astropy and tell people to use skimage.measure.block_reduce? I don't think that's a good solution because most astronomers aren't even aware of scikit-image. Should we copy skimage.measure.block_reduce in astropy.extern?

[hamogu]

I suggest to advertise scikit-image. The original docs for imutils had a sentence that we would only implement what is not available in numpy, scipy, or scikit-image.
We interface with external packages at other places in astropy, e.g. we use beautifulsoup to read HTML tables.
However, we definitely need to have some documentation that says: you might expect to find a function here that performs the following task and this is why we did not implement it in astropy. Use the package found here instead.

That being said, as long as we only need one or two functions from scikit-image they could live in extern as long as we have a big warning sign that they are only there temporarily for convenience and will be removed if we see that we need more functionality from scikit-image so that it becomes impractical to maintain a separate copy.
I think many astropy users use anaconda and they have scikit-image installed anyway.

[embray]

I also agree that for functionality that is already available there scikit-image should be promoted.

[Cadair]

if scikit-image dosen't have the upscale functionality why not see if they would be interested in adding it there? As more and more people are using packaging systems like anaconda, worrying about peoples access to packages is, imo, becoming less of a concern.

[embray]

I think maybe as @hamogu suggested, it can't hurt to include such features provisionally in Astropy if there is a strong need. But they should also be offered to scikit-image as perhaps an eventual permanent home. As long as it's something that doesn't have deeply embedded astronomy-specific functionality, which it sounds like upsample doesn't.

[mwcraig]

@larrybradley:

• Sounds like people are receptive to updating the docs to point to scikit-image...do you want to update this PR go that route?
• If not, the fail on travis is because the axis argument to np.sum has to be an integer, not a tuple, in numpy prior to 1.7 (fail is in numpy 1.6)

[hamogu]

Even though I said it above, I'll answer here again for the record:

• downsample -> change docs to point to scikit-image (might be worth mentioning skimage.measure.local_sum and skimage.transform.downscale_local_mean, too)
• upsample -> scikit-image also has an upsample function. Change the docs to point to that function. This can be achieved with the right parameters to skimage.transform.resize, most importantly order=0. This should have an example along the lines of:

In [24]: a
Out[24]: 
array([[ 0.  ,  0.09,  0.18,  0.27],
       [ 0.36,  0.45,  0.55,  0.64],
       [ 0.73,  0.82,  0.91,  1.  ]])

In [25]: resize(a, (6,8), order=0)
Out[25]: 
array([[ 0.  ,  0.  ,  0.09,  0.09,  0.18,  0.18,  0.27,  0.27],
       [ 0.  ,  0.  ,  0.09,  0.09,  0.18,  0.18,  0.27,  0.27],
       [ 0.36,  0.36,  0.45,  0.45,  0.55,  0.55,  0.64,  0.64],
       [ 0.36,  0.36,  0.45,  0.45,  0.55,  0.55,  0.64,  0.64],
       [ 0.73,  0.73,  0.82,  0.82,  0.91,  0.91,  1.  ,  1.  ],
       [ 0.73,  0.73,  0.82,  0.82,  0.91,  0.91,  1.  ,  1.  ]])

Note: skimage.transform.resize will also normalize the resized image, unless preserve_range=True is set. However, preserve_range is a keyword that was only recently added, that's why my example above (that I wrote on a laptop with an older version) is normalized to 0.

[perrygreenfield]

Let me advocate a contrary view. There are actually a few issues in the background that need clarification (i.e., we need to reconsider the currently adopted principles).

1. Don't duplicate functionality elsewhere (in this case sickit-image)
2. Don't add required dependencies outside of numpy.

These are at odds with each other, of course, and I think we need to be clearer on how to resolve the conflict between them in general.

The two could be satisfied by copying the code from scikit-image into astropy. That has a drawback of imposing a requirement to keep tabs on the original code in the event there are improvements or bug fixes there. Or we can relax our dependencies requirement (yes, we can make it optional, but for such basic functionality, that doesn't really seem like a workable solution).

One solution I am not in favor of is simply directing people to scikit-image for that functionality. As developers that are familiar with the scientific python landscape, that doesn't seem like such a big deal. But for people who are new to this, having basic image operations coming from a number of different sources (and also having potentially different API's or conventions) is likely to be a great annoyance.

I'd argue that even if we do allow a direct dependence on scikit-image, that it makes a lot of sense to repackage these disparate tools into one package by effectively collecting them into one namespace even when underlying imports come from a number of different areas. It also allows us to present a more unified API for these tools by wrapping them with adapters.

Finally, if for some reason we absolutely don't want a scikit-image dependency in this case, then the above code is better than no code at all. This is a perfect illustration of the perfect being the enemy of the good. Just accept it, and if someone wants a more general implementation, add it yourself! Rejecting additional functionality because it isn't as general as we would like doesn't help to encourage contributions.

[astrofrog]

Just one note (which I've expanded on in the mailing list) - we do ultimately want our own downsample/upsample wrappers that deal with WCS, even if they have optional dependencies under-the-hood.

[embray]

@perrygreenfield

Don't add required dependencies outside of numpy

I've advocated in the recent past for gradually relaxing this requirement, since I think the installation issues have improved a lot in recent years. So I thought I'd just put that out there--don't want to derail the discussion though.

[mwcraig]

A couple comments, then a concrete (if not particularly imaginative) proposal that can hopefully wrap this up promptly.

• This PR is trying to address a very concrete issue -- photutils, ccdproc and perhaps other affiliated packages need this sort of functionality, and it has hopped around a bit among the affiliated packages and then imageutils, but didn't end up merged into 1.0 for the reasons @larrybradley described in the initial message.
• I agree that wherever possible we should avoid reinventing the wheel, but...
• ...it would be nice for astropy to make using things like scikit-image for astronomical data as easy and fool-proof as possible.

Proposal:

• Implement downsample and upsample as thin wrappers around the appropriate scikit-image , ideally decorated with @support_nddata to ease use with NDData/derivatives and facilitate eventual handling of WCS...though if @larrybradley doesn't want to do that in this PR it can be implemented later.
• Provide only minimal docstrings for each function, referring readers to the appropriate scikit-image docs for details of the implementation.

I do think there is some value astropy can add here by wrapping and decorating.

@larrybradley -- what are your thoughts?

[larrybradley]

My current plans for this PR are to:

• change the function names to block_reduce and block_replicate
• make block_reduce a wrapper around skimage.measure.block_reduce, but with the following changes:
  • in cases where the image shape along any axis is not perfectly divisible by the block size, then trim the array (from the end) instead of padding the array, which is what the scikit-image version does
  • add the @support_nddata decorator
• keep the current implementation of block_replicate (i.e. upsample), but also add the @support_nddata decorator and extend it to 3D arrays (or generalize it to n-dimensional arrays, but I didn't see an obvious way to do that with the current implementation -- suggestions welcome!). Simply wrapping skimage.transform.resize will not work here because resize requires a float array in the range from -1 to 1 (and int arrays don't work either) [in general one has to be careful about scikit-image input dtypes].

One other feature that I would like to add to these functions is handling masks, but I will likely add that in a separate PR.

I completely agree that these functions should eventually support WCS, but that will likely be easier once gWCS is ready.

[mwcraig]

Sounds great to me!

[mwcraig]

Also, if you could ping me when it is ready for review that would be helpful.

[larrybradley]

This is ready for review.

CC: @mwcraig

[hamogu]

Do we need to wrap this a try except block in case the dependency is not available?

[mwcraig]

Looks like no....in coordinates, where scipy is needed for some things, the import is not wrapped in the functions that use it but the tests are wrapped to avoid failures due to imports.

See:

astropy/astropy/coordinates/matching.py

Line 402 in 321a41c

from scipy import spatial

[hamogu]

I like the way this looks (assuming Travis passes).

[mwcraig]

@larrybradley -- all of the test failures are because scikit-image isn't installed.

Looks like the correct place to add the dependency for travis would be at: https://github.com/astropy/astropy/blob/master/.continuous-integration/travis/setup_dependencies_common.sh#L33

For appveyor please add it at:
https://github.com/astropy/astropy/blob/master/appveyor.yml#L44

[mwcraig]

Could you please add a scikit-image import here in a try/except, and if the import fails, skip the tests that require scipy?

An example of how to do that is at:

astropy/astropy/coordinates/tests/test_matching.py

Line 20 in 321a41c

try:

[mwcraig]

Forgot to comment on the actual code -- that part looks fine!

[larrybradley]

Can someone please restart travis test 9179.14? I think that must be a transient failure.

[embray]

Done.

[larrybradley]

Thanks, @embray. All tests pass now.

[mwcraig]

@larrybradley -- looks good! Just read the code over one more time, I think this is ready to merge.

[wkerzendorf]

This looks very good! Thanks @larrybradley - merging