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.

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