Quick fix for drizzle cases with negative values in the pixmap transformation

[gbrammer]

I discovered some problematic behavior when trying to use drizzle to map to an output wcs that is contained within the input wcs, where the result was an empty image. That was a symptom, for which I think the more general cause is that the pixmap computed by drizzle.calc_pixmap.calc_pixmap(input_wcs, output_wcs) has negative values that trip pixel range checks within the tdriz C code, which then faults out without doing anything.

This would be fixed by specifying xmin/xmax/ymin/ymax values at runtime with drizzle.dodrizzle, but it would preferable to make the desired behavior more transparent by default. This PR computes updates to xmin/xmax/ymin/ymax for cases where the derived pixmap has negative values. I think ideally the fix would be implemented within the tdriz C code, since the fix here could result in expensive tests on large index matrices in cases of large input image dimensions.

[gbrammer]

OK, my fix doesn't seem to work for rotations between the input/output WCS where the simple min/max test isn't valid.

[jdavies-st]

Interesting issue. If you're getting negative values in the pixmap, it means that the output frame is too small for the input image(s). Though I see that is explicitly what you're trying to do.

That means for your use case, that the size of the pixmap is wrong. It should be sized to the output frame, not the input frame. It is currently sized to the input frame, as usually users intend to resample a bunch of smaller images onto a larger mosaic, so the pixmap is sized to the input frame(s). That's the assumption right now.

One solution is to have a check to see which is bigger, input frame or output frame in calc_pixmap(), and then trim to either or if there are negative values in one or the other, though I'm not sure what logic one would use, as rotations could make that difficult to figure out. Or have a switch in that function that allows the user to specify the behavior?

[jdavies-st]

You can see here

https://github.com/spacetelescope/drizzle/blob/master/drizzle/calc_pixmap.py#L39

that it sizes the grid to be the size of the first wcs object, the input frame size. But you want it to be sized to the second wcs object, the drizzled output frame.

[gbrammer]

Thanks for the feedback @jdavies-st. Whatever is different about drizzlepac.astrodrizzle.adrizzle.do_driz compared to drizzle.dodrizzle.dodrizzle, the former has my expected/desired behavior in that it appears to be agnostic about the relative size of the input/output pixel grids.

It doesn't seem like changing pixmap to have the size of the output grid will work, as there is a test that its dimensions match that of the input image here: https://github.com/spacetelescope/drizzle/blob/master/drizzle/src/cdrizzleapi.c#L221

FYI, my motivation here is having the drizzle capabilities in a more compact package for running on, e.g., a virtual machine. drizzlepac brings a lot of things with it, and even so, I was unsuccessful installing it on a compact "amazonlinux" Docker machine.

[jdavies-st]

Ah interesting @gbrammer. Of course the test could be changed too. 😉

But given that the version in drizzlepac works as you expect, it suggests that the code here should work the same way. Not sure if this belongs in the python layer or the C layer. @bernie-simon might have thoughts on that.

Could you provide a simple programatic test case that illustrates the difference in behavior between the drizzlepac version and this version?

Agree that a lightweight version of this is very desirable!

[gbrammer]

@jdavies-st: here's a demo script that shows the issue: test_drizzle.py.txt.

I found a fix in the c code that seems to work in the most recent commit above, skipping the "clip_bounds" step. I'm going to close the PR since my fixes are now hacks that shouldn't be in the main repo.

If you look in detail in the plot produced by the demo, you'll see that the footprints are a bit different between astrodrizzle/drizzle in the reverse "2->1" transformation that at least works to produce some output. I don't understand where that comes from since I would expect the wcs mapping for this simple square WCS to be pretty robust to any differences in the mapping implementation.