lsst.afw.math.binImage puts the corner of the input image in the corner of the output image. pfs.drp.stella.resampleKernelImage puts the center of the input image in the center of the output image. Consider the following example:

>>> import numpy as np
>>> from lsst.afw.image import ImageD
>>> from lsst.afw.math import binImage
>>> from lsst.geom import Box2I, Point2I, Point2D
>>> from pfs.drp.stella.images import getIndices, calculateCentroid
>>> from pfs.drp.stella import resampleKernelImage
>>> binning = 5
>>> box = Box2I(Point2I(-50, -50), Point2I(50, 50))
>>> image = ImageD(box)
>>> xx, yy = getIndices(image.getBBox())
>>> image.array[:] = np.exp(-0.5*(xx**2 + yy**2)/5.43**2)
>>> lsst = binImage(image, binning)
>>> pfs = resampleKernelImage(image, binning, Box2I(Point2I(-10, -10), Point2I(10, 10)))
>>> calculateCentroid(image)
namespace(x=-2.5493701622806684e-18, y=-1.1988050118056573e-17)
>>> calculateCentroid(lsst)
namespace(x=-40.39999999998449, y=-40.39999999998449)
>>> calculateCentroid(pfs)
namespace(x=-9.642385750780955e-18, y=-9.642172334948585e-18)

We created a PSF directly in the center of the input kernel, but lsst.afw.math.binImage moved the center so that it's not even centered on a pixel. However, pfs.drp.stella.resampleKernelImage puts the center exactly where it should be. In the figure below, the original oversampled image is on the left, the lsst image is in the middle, and the pfs image is on the right.

 Could you perhaps comment on the results of the followign example

from lsst.afw.image import ImageD
from lsst.afw.math import binImage
from lsst.geom import Box2I, Point2I, Point2D
from pfs.drp.stella.images import getIndices, calculateCentroid
from pfs.drp.stella import resampleKernelImagebinning = 10
box = Box2I(Point2I(-20, -20), Point2I(19, 19))
image = ImageD(box)xx, yy = getIndices(image.getBBox())
image.array[:] =np.ones((40,40))print('calculateCentroid(starting image): '+str(calculateCentroid(image)))lsst = binImage(image, binning)
pfs = resampleKernelImage(image, binning, Box2I(Point2I(-2, -2), Point2I(1, 1)))
print('calculateCentroid(after binImage + getBBox() minimum ): '+str(np.array([calculateCentroid(lsst).x,calculateCentroid(lsst).x])+20))
print('calculateCentroid(resampleKernelImage): '+str(calculateCentroid(pfs)))
plt.figure(figsize=(24,6))plt.subplot(131)
plt.title('via binImage')
plt.imshow(lsst.array,origin='lower')
plt.axhline(1.5,color='black')
plt.axvline(1.5,color='black')
cbar=plt.colorbar(fraction=0.046, pad=0.04)plt.subplot(132)
plt.title('via resampleKernelImage')
plt.imshow(pfs.array,origin='lower')
cbar=plt.colorbar(fraction=0.046, pad=0.04)
plt.axhline(1.6666666666666667,color='black')
plt.axvline(1.6666666666666667,color='black')

which creates the figure `example_constant_array.png` and output of 

calculateCentroid(starting image): namespace(x=-0.5, y=-0.5) 
calculateCentroid(after binImage + getBBox() minimum ): [1.5 1.5] 
calculateCentroid(resampleKernelImage): namespace(x=-0.33333333333333337, y=-0.33333333333333337)

In this example I create a constant array, 40x40 which I resample into 4x4, assuming oversampling of 10. This is somewhat more similar to my case in which I always dowsample by integer value (in the previous example 101x101 gets downsampled by factor of 5 into 21x21). I am quite confused by the result from the pfs algorithm, whilst understanding lsst output - I just get constant image again.
 

The more I look at the algorithm I think that there is a 0.5 pixel difference between the two which is creating all this confusion, i.e., standard what is the index of pixel question.

Because the image is constant, how the output image looks all comes down to how the output pixels overlap the input pixels. For the lsst image, all the output pixels overlap a full 10x10 pixels (the corner maps to the corner, and there's an exact multiple of the binning in the input image), and so the output image is constant. For the pfs image, not all output pixels are fully populated (the center maps to the center, which means 20,20 physical/LOCAL coordinate in the input maps to 2,2 physical/LOCAL in the output; therefore the column=2 physical/LOCAL output contains contributions from physical input columns 16 to 25 inclusive, so column=1 has 6 to 15, and column=0 has -4 to 5, but it is assumed there is no flux outside the physical bounds of the input image, which means -4 to -1 inclusive are zero, and only 0 through 5 inclusive are populated, hence it only gets 6/10 of the flux), and so the edges have 0.6 of the flux that the main image has, and the bottom-left corner has 0.6^2 = 0.36.

>>> (pfs.array[1,1]/pfs.array[1,0])**-1
0.6
>>> (pfs.array[1,1]/pfs.array[0,0])**-1
0.36

It occurred to me that your confusion might be due to being unaware of the pixel indexing conventions in use.

We find it useful to associate images with two related indexing conventions. The first is the LOCAL convention, which you can think of as physical coordinates (not strictly accurate, but good enough for now), where the bottom-left pixel is 0,0. For an image of size 4,5, the LOCAL coordinates run from 0,0 to 3,4 (inclusive). The second is the PARENT convention, or logical coordinates, where the bottom-left pixel is x0,y0 (often abbreviated as xy0). For an image of size 4,5 with an xy0 of 2,3, the PARENT coordinates run from 2,3 to 5,7 (inclusive).

The PARENT convention is so-named because it's what we use when we want to refer to the coordinate system of the full image from which we've taken a sub-image. We can have an image of size 100,100, and from it extract a sub-image, say of size 4,5 and with an xy0 of 2,3. When we operate on that sub-image (e.g., subImage.array[y,x]), we use LOCAL coordinates. When we want to talk about where something is in relation to the full image, we use the PARENT coordinates.

You can use image.getXY0() to get the xy0 of the image, and image.getBBox() to get the PARENT bounding box, and image.getBBox(lsst.afw.image.LOCAL) to get the LOCAL bounding box.

When we have a kernel image (as in a convolution kernel), where the concept of the center is important (because it corresponds to zero shift, or the center of a PSF), we set xy0 negative, so that 0,0 in the PARENT coordinate system is somewhere in the middle of our image (in the LOCAL coordinate system). So an image of size 4,5 with an xy0 of -2,-3 has PARENT coordinates that run from -2,-3 to 1,1 (inclusive), and the "center" (the 0,0 pixel in the PARENT frame) is at 2,3 in the LOCAL frame.

So when I say that "the center maps to the center", it means that the (center of the) 0,0 pixel in the PARENT frame of the input image maps to the (center of the) 0,0 pixel in the PARENT frame of the output image, and the borders of the output pixels on the input image follow from that. And "the corner maps to the corner" means that the (corner of the) 0,0 pixel in the LOCAL frame of the input image maps to the (corner of the) 0,0 pixel in the LOCAL frame of the output image, and the borders of the output pixels on the input image follow from that.

One subtlety is that when we're binning by an even factor, we can't put the center of the input image directly on the center of the output image without sub-pixel sampling, which we don't want to do in the rebinning operation. Instead for even binning factors, when we recenter ((recenterOversampledKernelImage) we move the center of the PSF to -0.5,-0.5 instead of 0,0 to account for this.

Code changes look fine. Minor comment made.

I apologize for not commenting here earlier. The changes are exactly as needed and described. Comment made by Paul was very helpful.

Merged.