python实现投影变换

def warp_perspective(src, M, width, height, origin_x, origin_y, flags=cv2.INTER_NEAREST, borderMode=cv2.BORDER_CONSTANT,
                     borderValue=0, dst=None):

    B_SIZE = 32

    if dst == None:
        dst = np.zeros((height, width, 3), dtype=src.dtype)

    # Set interpolation mode.
    interpolation = flags & cv2.INTER_MAX
    if interpolation == cv2.INTER_AREA:
        raise Exception('Area interpolation is not supported!')

    # Prepare matrix.    
    M = M.astype(np.float64)
    if not(flags & cv2.WARP_INVERSE_MAP):
        M = cv2.invert(M)[1]
    M = M.flatten()

    x_dst = y_dst = 0
    for y in range(-origin_y, height, B_SIZE):
        for x in range(-origin_x, width, B_SIZE):
            # Block dimensions.
            bw = min(B_SIZE, width - x_dst)
            bh = min(B_SIZE, height - y_dst)

            # To avoid dimension errors.
            if bw <= 0 or bh <= 0:
                break

            # View of the destination array.
            dpart = dst[y_dst:y_dst+bh, x_dst:x_dst+bw]

            # Original code used view of array here, but we're using numpy array's.
            XY = np.zeros((bh, bw, 2), dtype=np.int16)
            A = np.zeros((bh, bw), dtype=np.uint16)

            for y1 in range(bh):
                X0 = M[0]*x + M[1]*(y + y1) + M[2]
                Y0 = M[3]*x + M[4]*(y + y1) + M[5]
                W0 = M[6]*x + M[7]*(y + y1) + M[8]

                if interpolation == cv2.INTER_NEAREST:
                    for x1 in range(bw):
                        W = np.float64(W0 + M[6]*x1);
                        if W != 0:
                            W = np.float64(1.0)/W

                        X = np.int32((X0 + M[0]*x1)*W)
                        Y = np.int32((Y0 + M[3]*x1)*W)
                        XY[y1, x1][0] = np.int16(X)
                        XY[y1, x1][1] = np.int16(Y)
                else:
                    for x1 in range(bw):
                        W = np.float64(W0 + M[6]*x1);
                        if W != 0:
                            W = cv2.INTER_TAB_SIZE/W

                        X = np.int32((X0 + M[0]*x1)*W)
                        Y = np.int32((Y0 + M[3]*x1)*W)
                        XY[y1, x1][0] = np.int16((X >> cv2.INTER_BITS) + origin_x)
                        XY[y1, x1][1] = np.int16((Y >> cv2.INTER_BITS) + origin_y)
                        A[y1, x1] = np.int16(((Y & (cv2.INTER_TAB_SIZE-1))*cv2.INTER_TAB_SIZE + (X & (cv2.INTER_TAB_SIZE-1))))

            if interpolation == cv2.INTER_NEAREST:
                cv2.remap(src, XY, None, interpolation, dst=dpart,
                          borderMode=borderMode, borderValue=borderValue)
            else:
                cv2.remap(src, XY, A, interpolation, dst=dpart,
                          borderMode=borderMode, borderValue=borderValue)

            x_dst += B_SIZE
        x_dst = 0
        y_dst += B_SIZE

    return dst