目的

画像処理における座標変換では，変換後の画像の(i,j)のindexに該当する元画像の座標を逆変換を用いて求めて，Bilinear等の手法で補間するのが一般的です。

CやC++等では仕方なしに2重ループを回していたのですがNumpyとPythonではforループは極力使いたくないと思っていた所，効率の良さそうなコードをGithubのとあるページで発見したので自分なりに変えてみました。

言うまでも無いことですが，Affineだったり射影変換の場合はそれ専用のOpenCVのライブラリを使ったほうが早いです。

コードの流れ

１. 変換により参照される元画像の座標も求める

変換前の(i,j)の座標indexを記録した配列を作ります。
arangeを用いて配列を作成し，tileやrepeatを用いて重ね合わせればあっという間にできますね。

以下の例では画像中心cx,cyから差し引いて，正規化した後に球面極座標に変換するために円周率をかけてたりしますが同じことです。

xline = (np.arange(wid)-cx)*2*pi/wid
yline = -(np.arange(hei)-cy)*pi/hei
Xsline = np.tile(xline,hei)
Ysline = yline.repeat(wid,axis=0)

あとは，これらの配列に対して様々な非線形変換を当てはめればOKです。
正弦余弦や対数など簡単な関数ならnumpyに用意があるので配列にそのまま渡せば全て一様に処理してくれてfor文を書くよりも断然早いです。

２. 求めた座標を基に変換後の画像をBilinear補間で求める

先程求めた座標に対応するピクセル値を補間でもってきて最後に画像の形にReshapeしておしまいです。

以下のコードでは，ピクセルindexの最大値，最小値を超えた場合はループさせていますが，一般的には0などの値でうめる（Padding）するのが良いでしょう。

def _bilinear_interpolation_loop(frame, screen_coord):
    ''' Calculate 
    input: frame and its subpixel coordinate [x,y]
    '''
    frame_height,frame_width,frame_channel =frame.shape
    #uf = np.mod(screen_coord[0],1) * frame_width  # long - width
    #vf = np.mod(screen_coord[1],1) * frame_height  # lat - height

    x0 = np.floor(screen_coord[0]).astype(int)  # coord of pixel to bottom left
    y0 = np.floor(screen_coord[1]).astype(int)
    uf = screen_coord[0]  # long - width
    vf = screen_coord[1]  # lat - height
    x2 = np.add(x0, np.ones(uf.shape).astype(int))  # coords of pixel to top right
    y2 = np.add(y0, np.ones(vf.shape).astype(int))
    
    # Assume Loop 
    x2 = np.where(x2 > frame_width-1, 0, x2)
    y2 = np.where(y2 > frame_height-1, 0, y2)
    
    base_y0 = np.multiply(y0, frame_width)
    base_y2 = np.multiply(y2, frame_width)

    A_idx = np.add(base_y0, x0)
    B_idx = np.add(base_y2, x0)
    C_idx = np.add(base_y0, x2)
    D_idx = np.add(base_y2, x2)

    flat_img = np.reshape(frame, [-1, frame_channel])
    #print(flat_img.shape)
    A = np.take(flat_img, A_idx, axis=0)
    B = np.take(flat_img, B_idx, axis=0)
    C = np.take(flat_img, C_idx, axis=0)
    D = np.take(flat_img, D_idx, axis=0)
    
    wa = np.multiply(x2 - uf, y2 - vf)
    wb = np.multiply(x2 - uf, vf - y0)
    wc = np.multiply(uf - x0, y2 - vf)
    wd = np.multiply(uf - x0, vf - y0)
    #print(wa,wb,wc,wd)
    # interpolate
    AA = np.multiply(A.astype(np.float32), np.array([wa, wa, wa]).T)
    BB = np.multiply(B.astype(np.float32), np.array([wb, wb, wb]).T)
    CC = np.multiply(C.astype(np.float32), np.array([wc, wc, wc]).T)
    DD = np.multiply(D.astype(np.float32), np.array([wd, wd, wd]).T)
    out = np.reshape(np.round(AA + BB + CC + DD).astype(np.uint8), [frame_height, frame_width, 3])
    return out

実行例

以下のブログにて正距円筒図法の画像を元に全天球カメラを仮想的に回転させた際の画像を作成していました。