`add_array_2d` shape mismatch

[barentsen]

When adding a small 2d array to a larger one, the behaviour of add_array_2d seems ill-defined when the shapes do not align nicely. At the moment this crashes:

In [1]: from imageutils import add_array_2d
In [2]: import numpy as np
In [3]: add_array_2d(np.zeros((3,3)), np.ones((2,2)), (1, 1))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-111-5d52eedc3140> in <module>()
----> 1 imageutils.add_array_2d(np.zeros((3,3)), np.ones((2,2)), (1, 1))

/home/gb/bin/anaconda/lib/python2.7/site-packages/imageutils-0.0.dev106-py2.7-linux-x86_64.egg/imageutils/array_utils.pyc in add_array_2d(array_large, array_small, position)
    131         s_y, s_x, b_y, b_x = _get_slices(array_large.shape,
    132                                          array_small.shape, position)
--> 133         array_large[s_y, s_x] += array_small[b_y, b_x]
    134         return array_large
    135     else:

ValueError: operands could not be broadcast together with shapes (3,3) (2,2) (3,3) 

Whereas this works:

In [112]: imageutils.add_array_2d(np.zeros((3,3)), np.ones((2,2)), (0.4,0.4))
Out[112]: 
array([[ 1.,  1.,  0.],
       [ 1.,  1.,  0.],
       [ 0.,  0.,  0.]])

and this works:

In [113]: imageutils.add_array_2d(np.zeros((3,3)), np.ones((2,2)), (1.5, 1.5))
Out[113]: 
array([[ 0.,  0.,  0.],
       [ 0.,  1.,  1.],
       [ 0.,  1.,  1.]])

i.e. values between 0.5 and 1.4 fail in this example. What is the behaviour we ought to expect for position = (1, 1)?

[barentsen]

I'd like to add that I decided to review the behaviour of add_array_2d because I am experiencing the following problem with photutils.subtract_psf, which is related:

In [161]:  photutils.subtract_psf(np.zeros((10,10)), photutils.DiscretePRF(np.ones((7,7))), [[2,2.]], [1.])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-161-05cd08bf6117> in <module>()
----> 1 photutils.subtract_psf(np.zeros((10,10)), photutils.DiscretePRF(np.ones((7,7))), [[2,2.]], [1.])

/home/gb/dev/photutils/photutils/psf.pyc in subtract_psf(data, psf, positions, fluxes, mask)
    538         x = extract_array_2d(indices[1], psf.shape, position)
    539         psf_image = psf.eval(x, y, fluxes[i], position[0], position[1])
--> 540         data = add_array_2d(data, -psf_image, position)
    541     return data

/home/gb/dev/photutils/photutils/extern/imageutils/array_utils.pyc in add_array_2d(array_large, array_small, position)
    135         s_y, s_x, b_y, b_x = _get_slices(array_large.shape,
    136                                          array_small.shape, position)
--> 137         array_large[s_y, s_x] += array_small[b_y, b_x]
    138         return array_large
    139     else:

ValueError: operands could not be broadcast together with shapes (6,6) (5,5) (6,6) 

[cdeil]

@adonath I think you wrote this ... any thoughts?

Maybe the issue is that imageutils.add_array_2d uses the center as position ... probably it's better to use the edges instead (either via integers or slice objects), like scipy.ndimage and scikit-image do. I tried to think about a good API for image utils that operate on parts of images only at some points and started collecting some thoughts here but didn't get to a conclusion or good implementation ... now is the time to do better for astropy.image!

[adonath]

@barentsen Thank you for revealing this issue. I had quick look at the code and the problem is rather related to the fact that smaller array has an even shape. Originally this code was meant to be only used with odd shaped PSFs/PRFs, but the issue with subtract_psf shows, that there is a need to also handle even shaped cases properly, which can occur close to data boundaries, when only a fraction of the PSF/PRF overlaps with the data.

I agree with @cdeil to think about a general solution for handling image slices, bounding boxes, cutouts etc., because it is needed at so many places (only in photutils there are three different kind of implementations of similar functionality...).

[cdeil]

I now remember where we discussed this before ... in #4 I proposed to add a BoundingBox class.

I still think it would be a good idea to centralise the image slicing code (generating min / max from center and extent; handling edges and overlaps, ...) in a class. I'm not sure it should hold references to the large image array or the cutout image array ... I think only trying to rewrite the existing code that deals with slices in photutils, gammapy and ccdproc would reveal what is best to implement and for the end user.

[barentsen]

Should the generic problem of "adding arrays" perhaps be treated as a reprojection problem, rather than a slicing problem? At least for the purpose of photutils.subtract_psf, I would expect add_array_2d to behave as follows:

>>> imageutils.add_array_2d(np.zeros((4,4)), np.ones((3,3)), (1.5, 1.5))
array([[ 0.25,  0.5,  0.5, 0.25],
       [ 0.5,    1.,   1.,  0.5],
       [ 0.5,    1.,   1.,  0.5],
       [ 0.25,  0.5,  0.5, 0.25]])

... whereas right now the output is:

array([[ 0.,  0.,  0.,  0.],
       [ 0.,  1.,  1.,  1.],
       [ 0.,  1.,  1.,  1.],
       [ 0.,  1.,  1.,  1.]])

[adonath]

@barentsen The handling of subpixel positions, that you just described, is done internally in the photutils.DiscretePRF class, so there was no need for this functionality in subtract_psf and add_array_2d. Generally speaking you are totally right, that's the kind of behaviour you would expect for an add_array_2d function .

I have a quick fix for the issue with subtract_psf in https://github.com/adonath/imageutils/tree/fix_for_%2340, which should resolve the issue with even shaped PSFs/PRFs close to the data boundaries.

[astrofrog]

The original bug here is fixed by #42. @cdeil - I agree we may still want to find a better API, so feel free to open an issue for that if you have ideas.

[astrofrog]

This is now fixed in the Astropy version. @cdeil - feel free to open an issue on Astropy if you think the BoundingBox approach should be implemented.