[AstroPy] Inverting SIP distortions

[Jacob Vanderplas]

Hello,
I'm working on a problem that requires the inversion of SIP
corrections: that is, I have points in world coordinates that I want
to project onto an image plane.

My FITS header contains both the (A_p_q, B_p_q) coefficients for the
forward transformation, as well as the (AP_p_q, BP_p_q) required for
the backward transformation, according to the Shupe et al 2005
reference.

The astropy WCS object has methods sip_foc2pix() and sip_pix2foc(),
which appeared to be what I needed, but these don't seem to actually
invert each other. That is:

>>> from astropy.wcs import WCS
>>> wcs = WCS(header)
>>> Xpix = [[1, 1]]
>>> Xfoc = wcs.sip_pix2foc(Xpix, 0)
>>> print wcs.sip_foc2pix(Xfoc, 0)
[[ 1.2522118   1.74039233]]

I would have expected this to recover the input pixels, but it's
clearly not the case.  Is there any functionality in astropy that will
utilize the available AP/BP coefficients in my header to perform the
inverse SIP distortion?  Is there anything I'm doing wrong or missing
here?
Thanks,
   Jake