This is part of what we discussed privately but can be relevant:

Most existing interpolation methods (for example those in Scipy, that Raul initially referenced) are interpolation-based: assuming each pixel as a point in its center on the coordinate grid. Then finding the new grid values based on the positions of the centers of the pixels of the new grid.

But so far, Warp doesn't use this type of interpolation (even though it is so much easier to implement!) because it is terrible for low signal-to-noise images and generally interpolation always adds so many systematics and adds noise.

What Warp currently does is to put the new pixel over the old one, find the polygon within each input pixel that overlaps with each output pixel. calculate that polygon's area, and multiply the fractional area by the pixel value (to get the fractional value) and finally to put that fractional value into the value of the output pixel. There is thus no interpolation involved and everything is robust and accurate and the signal-to-noise increases, its like a small convolution.

But I am not against using a point-based interpolation method of rebinning if someone else is interested to add it (I just don't have time to implement it myself)!