yamaken1343’s blog

技術ブログもどき

pythonで疎な二次元データをバイリニア法で補完する

何かいい感じのライブラリがあったら教えて欲しい

データについて

例えば, 以下のようなデータが対象になります.

f:id:yamaken1343:20180628190305p:plain

気象庁のサイトから持ってきた気温のデータですが, 気温が示されていない部分を補完し, もっともらしい温度を出力します.

本来は緯度経度等から正しい位置を持って来るべきですが, とりあえず適当に配置します. 0は未定義の値とします.*1

[[20.9  0.   0.   0.   0.   0.   0.  21.8]
 [ 0.  21.7  0.  21.9  0.   0.  21.2  0. ]
 [ 0.  21.7  0.  21.9  0.   0.   0.  20.6]
 [ 0.  22.3  0.   0.   0.   0.   0.   0. ]
 [22.   0.  22.8  0.  20.7  0.   0.  21.2]
 [ 0.  22.6  0.   0.   0.   0.   0.   0. ]
 [ 0.   0.   0.   0.   0.   0.   0.   0. ]
 [20.2  0.  23.2  0.  23.2 22.6  0.  22. ]]

色で表すと以下のようになります. 暗い領域は未定義を示します. f:id:yamaken1343:20180628190145p:plain

結果

こうなります.

[[20.9 21.2 21.6 21.8 21.6 21.5 21.6 21.8]
 [21.3 21.7 21.8 21.9 21.7 21.4 21.2 21.4]
 [21.6 21.7 21.8 21.9 21.6 21.4 21.  20.6]
 [22.1 22.3 22.4 21.8 21.  21.2 20.9 20.9]
 [22.  22.3 22.8 21.8 20.7 20.9 21.  21.2]
 [21.9 22.6 22.8 22.2 21.5 21.7 21.6 21.5]
 [21.  22.2 23.  22.6 22.4 22.1 21.8 21.7]
 [20.2 21.7 23.2 23.2 23.2 22.6 22.3 22. ]]

f:id:yamaken1343:20180628190127p:plain

コード

import numpy as np


def search_near(src, x, y):
    """
    データのある近傍4点を検索し, 座標を返す
    :param src: 検索する配列
    :param x: 基準点のx座標
    :param y: 基準点のy座標
    :return: 検索した4点を2*4のnumpy.arrayで返す
    """

    def min_pear(a, b, x, y):
        """
        基準点とリスト内の2点間の距離が最も小さいペアを返す
        :param a: x座標のリスト
        :param b: y座標のリスト
        :param x: 基準点のx
        :param y: 基準点のy
        :return: リスト内のペア
        """
        # 近傍に点が4点以下でもあるだけの点を用いて補完を行う
        # if len(a) == 0:
        #     return x, y

        add = np.power((a - x), 2) + np.power((b - y), 2)
        return a[np.argmin(add)], b[np.argmin(add)]

    # 処理簡便化のためデータの有無を2値で持つ
    src = np.array(src != 0)

    # 処理対象の点にデータがあるとき, その点を4近傍として返す
    if src[x, y]:
        return np.array([[x, y], [x, y], [x, y], [x, y]])

    try:
        # スライスされるため, indexに入る値はsrcのインデックスと互換性がないことに注意する
        # 左上
        index_x, index_y = np.where(src[0:x, 0:y])
        min_x, min_y = min_pear(index_x, index_y, x, y)  # 基準点は右下
        upper_left = (min_x, min_y)
        # 右上
        index_x, index_y = np.where(src[x:, 0:y])
        min_x, min_y = min_pear(index_x, index_y, 0, y)
        upper_right = (min_x + x, min_y)
        # 左下
        index_x, index_y = np.where(src[0:x, y:])
        min_x, min_y = min_pear(index_x, index_y, x, 0)
        lower_left = (min_x, min_y + y)
        # 左下
        index_x, index_y = np.where(src[x:, y:])
        min_x, min_y = min_pear(index_x, index_y, 0, 0)
        lower_right = (min_x + x, min_y + y)

        near_4_points = (upper_left, upper_right, lower_left, lower_right)

        return np.array(near_4_points)

    # データが見つからなくても止まらないようにする
    except:
        return None


def bilinear(src):
    """
    疎なデータを持つ配列に対しバイリニア補完を行う
    :param src: 入力配列
    :return: 補完後の配列
    """
    it = np.nditer(src, flags=['multi_index'])
    dst = np.zeros(src.shape)
    # 配列の各要素をなめる
    while not it.finished:
        idx = it.multi_index  # 現在参照する要素のインデックス
        it.iternext()

        # near four points 近傍四点を取得
        n4p = search_near(src, idx[0], idx[1])
        # 4点のどれかにデータがない場合の処理
        if n4p is None:
            continue

        # 参照する要素と近傍4点の距離を取得
        near_4_points_dist = np.array([np.linalg.norm(n4p[0] - idx), np.linalg.norm(n4p[1] - idx),
                                       np.linalg.norm(n4p[2] - idx), np.linalg.norm(n4p[3] - idx)])

        # 近傍点と参照点が同じ場合と近傍点がないときnanが入るので, 変換する
        near_4_points_dist[np.isnan(near_4_points_dist)] = 0

        # 参照要素とデータの点が重なった時の処理
        if near_4_points_dist.sum() == 0:
            dst[idx] = src[idx]
            continue

        # 近傍点と参照点が同じ場合無視したいので次の処理でゼロにするために代入
        near_4_points_dist[near_4_points_dist == 0] = np.inf

        # 小さいほうが値が大きくなるように逆数を取る(距離なので近いほうが良いというふうにしたい)
        near_4_points_score = 1 / near_4_points_dist

        # 和が1になるように正規化
        near_4_points_score /= near_4_points_score.sum()

        # 結果を格納する配列に代入
        for i in range(len(n4p)):
            a, b = n4p[i]
            dst[idx] += near_4_points_score[i] * src[a, b]

    return dst

def main():
    np.set_printoptions(precision=1)
    in = np.array([[20.9, 0, 0, 0, 0, 0, 0, 21.8],
                  [0, 21.7, 0, 21.9, 0, 0, 21.2, 0],
                  [0, 21.7, 0, 21.9, 0, 0, 0, 20.6],
                  [0, 22.3, 0, 0, 0, 0, 0, 0],
                  [22.0, 0, 22.8, 0, 20.7, 0, 0, 21.2],
                  [0, 22.6, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0],
                  [20.2, 0, 23.2, 0, 23.2, 22.6, 0, 22.0]])
    print(in)
    out = bilinear(in)
    print(out)


if __name__ == '__main__':
    main()

解説

7/20変更 左上領域右上領域左下領域右下領域が重ならないように範囲を設定しています
min_pearのすぐ下のコメントを外すと4点未満でもあるだけの点を用いて補完を行うようになります.
ヒートマップを作成するコードは省いています.

*1:本来はNoneとかを使うべき