numpyを使って複数の矩形のIoUを一度に高速に計算する

pythonのnumpyライブラリを使って、1つの矩形と複数の矩形とのIoU (Intersection over Union)を一度に高速に計算する方法を紹介します。計算時間の計測結果も記載し、1つずつIoUを計算した場合に比べてどのくらい高速化できるのか比較も行います。

この記事でできること

pythonのnumpyライブラリを使って、ある1つの矩形に対して複数の矩形のIoU (Intersection over Union)を一気に計算する方法を紹介します。IoUはコンピュータビジョンやニューラルネットワーク(ディープラーニング)等の機械学習を用いた物体検出AIで重要となる概念です。

pythonのnumpyでブロードキャスト演算を使うと、for文を使わずに複数の演算を一度に記述することができます。ブロードキャスト演算は記述が簡潔になるだけでなく、高速化も図ることができます。

この記事では、IoU計算をブロードキャスト演算で記述し複数のIoUを一度に計算する実装例を紹介します。さらに、for文とpythonのリスト(list、配列)を使って1つずつIoUを計算した場合との計算時間の比較を行います。

結論としては、(かなり偏ったケースであるものの)約26倍高速化できました。

IoU(Intersection over Union)について

IoUとその実装例については次の過去記事で紹介しました。

この過去記事では、numpyを使わずにpythonのリストを使って1つずつIoUを計算する関数”iou()”を掲載しています。

本記事では、すぐ次の節でnumpyを使った複数IoUを一度に計算するサンプルコードを紹介しますが、その計算時間と過去記事の1つずつIoUを計算する関数”iou()”との計算時間を、このページの下の方で比較したいと思います。

複数の矩形のIoUを一度に計算するnunpyの実装例

numpyを使った複数IoUを一度に計算する関数のサンプルコードは次の通りです。1つの矩形を表すnumpy配列”a”と、複数の矩形が格納されたnumpy配列”b”とを引数として入力し、矩形”a”と”b”に含まれる全ての矩形とのIoUをfor文無しで計算します。

import numpy as np

# 矩形aと、複数の矩形bのIoUを計算
def iou_np(a, b):
    # aは1つの矩形を表すshape=(4,)のnumpy配列
    # array([xmin, ymin, xmax, ymax])
    # bは任意のN個の矩形を表すshape=(N, 4)のnumpy配列
    # 2次元目の4は、array([xmin, ymin, xmax, ymax])
    
    # 矩形aの面積a_areaを計算
    a_area = (a[2] - a[0] + 1) \
             * (a[3] - a[1] + 1)
    # bに含まれる矩形のそれぞれの面積b_areaを計算
    # shape=(N,)のnumpy配列。Nは矩形の数
    b_area = (b[:,2] - b[:,0] + 1) \
             * (b[:,3] - b[:,1] + 1)
    
    # aとbの矩形の共通部分(intersection)の面積を計算するために、
    # N個のbについて、aとの共通部分のxmin, ymin, xmax, ymaxを一気に計算
    abx_mn = np.maximum(a[0], b[:,0]) # xmin
    aby_mn = np.maximum(a[1], b[:,1]) # ymin
    abx_mx = np.minimum(a[2], b[:,2]) # xmax
    aby_mx = np.minimum(a[3], b[:,3]) # ymax
    # 共通部分の矩形の幅を計算。共通部分が無ければ0
    w = np.maximum(0, abx_mx - abx_mn + 1)
    # 共通部分の矩形の高さを計算。共通部分が無ければ0
    h = np.maximum(0, aby_mx - aby_mn + 1)
    # 共通部分の面積を計算。共通部分が無ければ0
    intersect = w*h
    
    # N個のbについて、aとのIoUを一気に計算
    iou = intersect / (a_area + b_area - intersect)
    return iou

演算内容は普通のIoU計算と同じです(当然)。numpyの関数を使用することで複数の演算(加算や乗算、max/min)を同時に行っています。

このコードで使用している矩形の座標(xmin, ymin, xmax, ymax)の意味については、過去記事「pythonで2つの矩形(長方形)のIoUを計算する」で説明していますので参考にしてください。物体検出AIでよく使われる表現方法となっています。

IoUの計算時間を計測

上で紹介した、複数IoUを一度に計算する関数”iou_np(a, b)”と、過去記事の1つずつIoUを計算する関数”iou()”の処理時間をそれぞれ計測します。

Windowsでプログラムの処理時間を計測する際、OSタイマーの割込時間に左右されて正確に時間を計測できないことが多々あります。本記事では以前紹介したOSタイマーの変更方法(この段落のすぐ下のリンク)を使用して、OSタイマーをデフォルト(15.625ミリ秒?)から1ミリ秒に設定し、より正確に計測していきます。

後述のサンプルコード内の関数”windll.winmm.timeBeginPeriod(1)”や、”windll.winmm.timeEndPeriod(1)”の意味や使い方を説明した記事になります。

ただし、結果として気休め程度の処理となりました。時間計測をする上で必須なわけではありません。

複数のIoUをnumpyで一気に計算する場合の処理時間

次のサンプルコードで、複数IoUを一度に計算する関数”iou_np(a, b)”の計算時間を計測します。

if __name__ == '__main__':
    import numpy as np
    import time
    
    from ctypes import windll #new
    # タイマー精度を1msec単位にする
    windll.winmm.timeBeginPeriod(1)
    

    np.random.seed(314)
    
    ious = [] # 計算したIoUを格納するリスト
    
    # IoUを計算する任意の矩形
    a = [100, 500, 200, 600]
    a_np = np.array(a)
        
    time_list = [] # IoU処理の時間を格納する配列(リスト)の初期化
    for dummy in range(100): # 100回繰り返し
    
    
        bboxes_np = np.empty((0, 4)) # 乱数で生成した格納する配列(リスト)の初期化
        for i in range(1000): # 1000個の矩形を乱数で生成
            # 矩形aの座標を0～100範囲でずらした矩形を乱数で生成
            rnd_bbox = a_np + (np.random.rand(4) - 0.5)*200
            rnd_bbox = rnd_bbox.reshape(1, -1)
            bboxes_np = np.append(bboxes_np, rnd_bbox, axis=0)

        start = time.perf_counter() # 時間計測スタート
        _iou = iou_np(a_np, bboxes_np)
        end = time.perf_counter() # 時間計測終わり
         # 矩形aと1000個の矩形のIoUの処理時間
        time_sec = end - start # スタートから終わりまでの時間[秒]
        
        ious.append(_iou)
        
        time_list.append(time_sec*1000) # ミリ秒で配列time_listに追加
    
    print("Average time: ", sum(time_list)/len(time_list), " [msec]")
    
    print("bboxes_np", bboxes_np[0][0:10])
    
    # タイマー精度を戻す
    windll.winmm.timeEndPeriod(1)
    
    
    import matplotlib.pyplot as plt
    
    # ヒストグラムの横軸。0.1刻み0～3の範囲の配列(リスト)
    bin_list1 = [i*0.1 for i in range(31)]
    #print("bin_list1", bin_list1)
    
    plt.hist(time_list, bins=bin_list1) # ヒストグラム生成
    plt.savefig("histogram_time1_fast.png") # 画像として出力
    plt.clf() # 図をクリア # 計算したIoUを格納するリスト

“np.random.seed(314)”は乱数シードの設定です。乱数シードについては、過去記事「pythonで乱数シード(seed)を設定する３種の方法」で紹介しました。

このサンプルコードでは、矩形”a= [100, 500, 200, 600]”を一つ定義し、矩形”a”の座標に基づいてランダムに1000個の矩形を生成してnumpy配列”b”に格納し、”a”と”b”のIoUを計算します。関数”iou_np(a, b)”で計算しているのは、矩形”a”と1000個の矩形が格納された”b”とのそれぞれのIoUです。

上のコードでは、この演算を100回繰り返した計算時間の平均をprint文で出力しています。私の環境では次のように約0.09ミリ秒(msec)となりました。

Average time:  0.09082099999999205  [msec]

また、このサンプルコードの後半では、pythonのmatplotlibライブラリを使用し、100回繰り返した計算時間それぞれをヒストグラムとして画像出力しています。

ヒストグラムの生成方法や実装例については次の過去記事で紹介しましたので、参照ください。

出力されたヒストグラムは次のようになりました。横軸が処理時間[ミリ秒(msec)]で縦軸が度数です。

0.0～0.2ミリ秒の間に固まったヒストグラムとなっています。平均が約0.09ミリ秒ですので、妥当なヒストグラムかなと思います。

1つずつIoUを計算する場合の処理時間

次に、pythonのリストを使って1つずつIoUを計算する関数”iou()”の計算時間を計測します。

    import random
    import time
    random.seed(314)
    
    from ctypes import windll #new
    # タイマー精度を1msec単位にする
    windll.winmm.timeBeginPeriod(1)
    
    
    # IoUを計算する任意の矩形
    a = [100, 500, 200, 600]
    
    ious = [] # 計算したIoUを格納するリスト
    
    time_list = [] # IoU処理の時間を格納する配列(リスト)の初期化    
    for dummy in range(100): # 100回実行
    
    
        bboxes = [] # 乱数で生成した格納する配列(リスト)の初期化
        for i in range(1000): # 1000個の矩形を乱数で生成
            # 矩形aの座標を0～100範囲でずらした矩形を乱数で生成
            rnd_bbox = [e + (-100 + 200*random.random()) for e in a]
            bboxes.append(rnd_bbox)

        time_sec = 0 # 矩形aと1000個の矩形のIoUの処理時間初期化
        # 1000回iou(a, box)を実行
        for box in bboxes:
            start = time.perf_counter() # 時間計測スタート
            _iou = iou(a, box)
            end = time.perf_counter() # 時間計測終わり
            time_sec += end - start # スタートから終わりまでの時間[秒]
            
            ious.append(_iou)
        
        
        time_list.append(time_sec*1000) # ミリ秒で配列time_listに追加
    
    print("Average time: ", sum(time_list)/len(time_list), " [msec]")
    
    print("bboxes", bboxes[0:10])
    
    # タイマー精度を戻す
    windll.winmm.timeEndPeriod(1)
    
    
    import matplotlib.pyplot as plt
    
    # ヒストグラムの横軸。0.1刻み0～3の範囲の配列(リスト)
    bin_list1 = [i*0.1 for i in range(31)]
    print("bin_list1", bin_list1)
    
    plt.hist(time_list, bins=bin_list1) # ヒストグラム生成
    plt.savefig("histogram_time1.png") # 画像として出力
    plt.clf() # 図をクリア    

このコードもnumpyでIoU計算した時と同様に、矩形”a= [100, 500, 200, 600]”を一つ定義し、矩形”a”の座標に基づいてランダムに1000個の矩形を生成しています。矩形”a”とリスト”b”内の1000個の矩形とのIoUそれぞれを、for文と関数”iou(a, box)”を使って1つずつ計算していきます。

1000個のIoU計算を100回繰り返した計算時間の平均がprint文で出力されます。私の環境では次のように約2.44ミリ秒(msec)となりました。

Average time:  2.4406220000000287  [msec]

一つ上の節で紹介した、「複数のIoUをnumpyで一気に計算する場合」は平均約0.09ミリ秒でしたので、約26倍高速化できたことになります。

ヒストグラムは次のように出力されました。横軸が処理時間[ミリ秒(msec)]、縦軸が度数です。

平均約2.44ミリ秒に対し、2.0～2.3ミリ秒の間に固まったヒストグラムとなりました。仕方がないことですが、処理時間にばらつきがあるようです。

おわりに

pythonのnumpyライブラリを使って、1つの矩形と複数の矩形のIoUを一気に、高速に計算する方法を紹介しました。

1000個のIoUを計算するのに、平均で約2.44ミリ秒から約0.09ミリ秒約26倍高速化できました。（1ケースの高速化例ですので、実際はここまで高速化できないケースもあるかと思います）

IoUはコンピュータビジョンやニューラルネットワーク(ディープラーニング)等の機械学習を用いた物体検出AIで重要となる概念で、検出された矩形が正解であるかの判定に使われたり、検出された複数の矩形を統合する際の指標に使われたりしています。

検出された複数の矩形を統合するアルゴリズムであるNMS (Non-Maximum Suppression)や、その発展形であるsoft-NMSについては、以前実装例などを紹介しましたので、下にリンクを貼っておきます。

この記事で掲載した高速なIoUの計算の発展形として、numpyライブラリを使ったNMSの高速な実装例も次に紹介しています。20倍?くらい高速化できました。

他に、NMSの改良手法のsoft-NMSの高速化についても今後記載したいと思います。