def rectify2d(image2d, aij, bij, resampling,
              naxis1out=None, naxis2out=None,
              ioffx=None, ioffy=None,
              debugplot=0):
    """Rectify image applying the provided 2D transformation.

    The rectified image correspond to the transformation given by:
        u = sum[i=0:order]( sum[j=0:j]( a_ij * x**(i - j) * y**j ))
        v = sum[i=0:order]( sum[j=0:j]( b_ij * x**(i - j) * y**j ))

    Parameters
    ----------
    image2d : 2d numpy array
        Initial image.
    aij : 1d numpy array
        Coefficients a_ij of the transformation.
    bij : 1d numpy array
        Coefficients b_ij of the transformation.
    resampling : int
        1: nearest neighbour, 2: flux preserving interpolation.
    naxis1out : int or None
        X-axis dimension of output image.
    naxis2out : int or None
        Y-axis dimension of output image.
    ioffx : int
        Integer offset in the X direction.
    ioffy : int
        Integer offset in the Y direction.
    debugplot : int
        Determines whether intermediate computations and/or plots
        are displayed. The valid codes are defined in
        numina.array.display.pause_debugplot

    Returns
    -------
    image2d_rect : 2d numpy array
        Rectified image.

    """

    # protections
    ncoef = len(aij)
    if len(bij) != ncoef:
        raise ValueError("aij and bij lengths are different!")

    # order of the polynomial transformation
    order = order_fmap(ncoef)
    if abs(debugplot) >= 10:
        print('--> rectification order:', order)

    # initial image dimension
    naxis2, naxis1 = image2d.shape

    # output image dimension
    if naxis1out is None:
        naxis1out = naxis1
    if naxis2out is None:
        naxis2out = naxis2
    if ioffx is None:
        ioffx = 0
    if ioffy is None:
        ioffy = 0

    if resampling == 1:
        # pixel coordinates (rectified image); since the fmap function
        # below requires floats, these arrays must use dtype=np.float
        j = np.arange(0, naxis1out, dtype=np.float) - ioffx
        i = np.arange(0, naxis2out, dtype=np.float) - ioffy
        # the cartesian product of the previous 1D arrays could be stored
        # as np.transpose([xx,yy]), where xx and yy are computed as follows
        xx = np.tile(j, (len(i),))
        yy = np.repeat(i, len(j))
        # compute pixel coordinates in original (distorted) image
        xxx, yyy = fmap(order, aij, bij, xx, yy)
        # round to nearest integer and cast to integer; note that the
        # rounding still provides a float, so the casting is required
        ixxx = np.rint(xxx).astype(np.int)
        iyyy = np.rint(yyy).astype(np.int)
        # determine pixel coordinates within available image
        lxxx = np.logical_and(ixxx >= 0, ixxx < naxis1)
        lyyy = np.logical_and(iyyy >= 0, iyyy < naxis2)
        lok = np.logical_and(lxxx, lyyy)
        # assign pixel values to rectified image
        ixx = xx.astype(np.int)[lok]
        iyy = yy.astype(np.int)[lok]
        ixxx = ixxx[lok]
        iyyy = iyyy[lok]
        # initialize result
        image2d_rect = np.zeros((naxis2out, naxis1out), dtype=np.float)
        # rectified image
        image2d_rect[iyy + ioffy, ixx + ioffx] = image2d[iyyy, ixxx]
    elif resampling == 2:
        # coordinates (rectified image) of the four corners, sorted in
        # anticlockwise order, of every pixel
        j = np.array(
            [[k - 0.5 - ioffx, k + 0.5 - ioffx,
              k + 0.5 - ioffx, k - 0.5 - ioffx] for k in range(naxis1out)]
        )
        i = np.array(
            [[k - 0.5 - ioffy, k - 0.5 - ioffy,
              k + 0.5 - ioffy, k + 0.5 - ioffy] for k in range(naxis2out)]
        )
        xx = np.reshape(np.tile(j, naxis2out), naxis1out * naxis2out * 4)
        yy = np.concatenate([np.reshape(i, naxis2out * 4)] * naxis1out)
        # compute pixel coordinates in original (distorted) image
        xxx, yyy = fmap(order, aij, bij, xx, yy)
        # indices of pixels in the rectified image
        ixx = np.repeat(np.arange(naxis1out), naxis2out)
        iyy = np.tile(np.arange(naxis2out), (naxis1out,))
        # rectified image (using cython function)
        image2d_rect = _resample(image2d, xxx, yyy, ixx, iyy,
                                 naxis1out, naxis2out)
    else:
        raise ValueError("Sorry, resampling method must be 1 or 2")

    # return result
    return image2d_rect
def rectify2d ( image2d , aij , bij , resampling , naxis1out = None , naxis2out = None , ioffx = None , ioffy = None , debugplot = 0 ) : ncoef = len ( aij ) if len ( bij ) != ncoef : raise ValueError ( "aij and bij lengths are different!" ) order = order_fmap ( ncoef ) if abs ( debugplot ) >= 10 : print ( '--> rectification order:' , order ) naxis2 , naxis1 = image2d . shape if naxis1out is None : naxis1out = naxis1 if naxis2out is None : naxis2out = naxis2 if ioffx is None : ioffx = 0 if ioffy is None : ioffy = 0 if resampling == 1 : j = np . arange ( 0 , naxis1out , dtype = np . float ) - ioffx i = np . arange ( 0 , naxis2out , dtype = np . float ) - ioffy xx = np . tile ( j , ( len ( i ) , ) ) yy = np . repeat ( i , len ( j ) ) xxx , yyy = fmap ( order , aij , bij , xx , yy ) ixxx = np . rint ( xxx ) . astype ( np . int ) iyyy = np . rint ( yyy ) . astype ( np . int ) lxxx = np . logical_and ( ixxx >= 0 , ixxx < naxis1 ) lyyy = np . logical_and ( iyyy >= 0 , iyyy < naxis2 ) lok = np . logical_and ( lxxx , lyyy ) ixx = xx . astype ( np . int ) [ lok ] iyy = yy . astype ( np . int ) [ lok ] ixxx = ixxx [ lok ] iyyy = iyyy [ lok ] image2d_rect = np . zeros ( ( naxis2out , naxis1out ) , dtype = np . float ) image2d_rect [ iyy + ioffy , ixx + ioffx ] = image2d [ iyyy , ixxx ] elif resampling == 2 : j = np . array ( [ [ k - 0.5 - ioffx , k + 0.5 - ioffx , k + 0.5 - ioffx , k - 0.5 - ioffx ] for k in range ( naxis1out ) ] ) i = np . array ( [ [ k - 0.5 - ioffy , k - 0.5 - ioffy , k + 0.5 - ioffy , k + 0.5 - ioffy ] for k in range ( naxis2out ) ] ) xx = np . reshape ( np . tile ( j , naxis2out ) , naxis1out * naxis2out * 4 ) yy = np . concatenate ( [ np . reshape ( i , naxis2out * 4 ) ] * naxis1out ) xxx , yyy = fmap ( order , aij , bij , xx , yy ) ixx = np . repeat ( np . arange ( naxis1out ) , naxis2out ) iyy = np . tile ( np . arange ( naxis2out ) , ( naxis1out , ) ) image2d_rect = _resample ( image2d , xxx , yyy , ixx , iyy , naxis1out , naxis2out ) else : raise ValueError ( "Sorry, resampling method must be 1 or 2" ) return image2d_rect
Rectify image applying the provided 2D transformation.

    The rectified image correspond to the transformation given by:
        u = sum[i=0:order]( sum[j=0:j]( a_ij * x**(i - j) * y**j ))
        v = sum[i=0:order]( sum[j=0:j]( b_ij * x**(i - j) * y**j ))

    Parameters
    ----------
    image2d : 2d numpy array
        Initial image.
    aij : 1d numpy array
        Coefficients a_ij of the transformation.
    bij : 1d numpy array
        Coefficients b_ij of the transformation.
    resampling : int
        1: nearest neighbour, 2: flux preserving interpolation.
    naxis1out : int or None
        X-axis dimension of output image.
    naxis2out : int or None
        Y-axis dimension of output image.
    ioffx : int
        Integer offset in the X direction.
    ioffy : int
        Integer offset in the Y direction.
    debugplot : int
        Determines whether intermediate computations and/or plots
        are displayed. The valid codes are defined in
        numina.array.display.pause_debugplot

    Returns
    -------
    image2d_rect : 2d numpy array
        Rectified image.