function warp, im, slice, interp, UNWARP=unwarp
    w = (size(im))[1]
    ; the matrix [[a, b],[c,d]] is the warp transform
    a = slice.a & b = slice.b & c = slice.c & d = slice.d
    if KEYWORD_SET(unwarp) then begin
        m = INVERT([[a,b],[c,d]])
        a = m[0,0] & b = m[1,0] & c = m[0,1] & d = m[1,1]
    endif
    det = a*d-b*c
    ; we need to shift the coordinates by w/2 before and after the warp
    ; this corresponds to a post-warp shift of
    s = (1/det)*[[d,-b],[-c,a]]##[[-w/2.],[-w/2.]]+[[w/2.],[w/2.]]
    ; remap using the inverse transform as linear coordinate lookup polynomials
    P = [[s[0],-b/det],[d/det,0]]
    Q = [[s[1],a/det],[-c/det,0]]
    if interp eq 2 then $
        return, poly_2d(im, P, Q, 2, CUBIC=-.5) $
    else $
        return, poly_2d(im, P, Q, interp)
end

; perform a flux-preserving warp
; this is nominally smarter than using interpolation
; convolve by the shape of a pixel before sampling
function warp_high_quality, im, slice
    pixel = DBLARR(15,15)
    pixel[5:9,5:9] = 1d
    unwarped_pixel = warp(pixel, slice, 1, /UNWARP)
    im2 = CONVOL(im, REBIN(unwarped_pixel,3,3))
    return, warp(im2, slice, 0)
end

; perform a point-preserving warp
; this has the benefit that stars remain pointlike
; flux is also conserved, but this routine would be
; terrible for continuous images
function warp_points, im, slice
    w = (size(im))[1]
    pts = where(im ne 0d)
    if pts[0] lt 0 then return, im
    x = pts mod w
    y = pts / w
    vals = im[pts]

    xyw = [[x],[y],[DBLARR(N_ELEMENTS(pts))+1d]]
    T = [[1,0,w/2+1],[0,1,w/2+1]] ## [[slice.a,slice.b,0],[slice.c,slice.d,0],[0,0,1]] ## [[1,0,-w/2],[0,1,-w/2],[0,0,1]]
    uv = T ## xyw
    return, hist2d(uv[*,0], uv[*,1], vals, MIN1=1.,MIN2=1.,MAX1=w,MAX2=w)
end