3D Hausdorff Distance (percentile付き)の算出方法(Python)

調べてもぱっと出てこなくて困ったので書いておきます。

(知らない人向け)そもそもHausdorff Distanceって何？

セグメンテーションの評価指標というと、DiceやIoU,PresicionやRecallなどがありますが、これに加えてハウスドルフ距離という評価指標も使用する場合があります。ハウスドルフ距離は、正解ラベル/予測したマスクのある点から、最も近い予測値したマスク/正解ラベルまでの距離の最大値で表されます [1] 。セグメンテーション対象の境界の精度が重要である医療画像の分野において評価指標として使用される場合が多いです [2] [3]。

必要となった経緯

CT画像のセグメンテーションの評価指標としてHausdorff距離を使用する際に、先行研究論文の多くで3次元ベースの90~95% Hausdorff Distanceを使用しており、これに習ってHDを算出したかったからです。
3次元ベースのHD算出コードが見つからない、パーセンタイルのHD算出コードが見つからないで結構探すのに苦労しました。

実装

環境はKaggle Notebookです(2021/09/27)。
実装Notebookはこちらです。

MONAIというライブラリを使用することで簡単に算出できます。
pytorchとnumpyが必要で、monaiからの戻り値はpytorchのTensor型となっています。

#https://docs.monai.io/en/latest/installation.html
#インストール
!pip install monai

必要なライブラリをインポートします。
後で3D描画を行なうので、色々入ってます。

import numpy as np
from monai.metrics import compute_hausdorff_distance
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
from skimage import measure
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

import warnings
warnings.filterwarnings('ignore')

さて、まずは2つの3次元の配列arr1,arr2のある2点についてHDを算出してみます。
実際にセグメンテーションのコードを動かすことを想定して、配列の次元はそれぞれ、(batch, n_classes, X, Y, Z)とすることにします。(というか、MONAIの想定する入力がこれです)
今回は、バッチサイズは1,クラス数は2,(X,Y,Z)方向に32ピクセルずつ存在する三次元配列を作成します。

次に、双方の各クラスの適当な場所に1点だけ1を格納します。この点同士の3次元上でのユークリッド距離を算出することになります。
ユークリッド距離を手動で算出した後に、MONAIのcompute_hausdorff_distance関数を動かして、結果が同じであることを確認します。引数についてはこちらのページ を参照してください。

#HD算出
n_classes = 2
#(batch,n_classes,X,Y,Z) array must be one-hot format
arr1 = np.zeros((1,n_classes,32,32,32),dtype=np.uint8) 
arr2 = np.zeros((1,n_classes,32,32,32),dtype=np.uint8)

ch = 0
arr1[0,ch,0,0,0] = 1
arr2[0,ch,8,8,8] = 1

ch = 1
arr1[0,ch,0,0,0] = 1
arr2[0,ch,16,16,16] = 1

print("manual : ")
# Euclidean distance between two points. 2点のユークリッド距離の手動計算
print(f"ch0:{round((((8-0)**2)*3)**0.5,4)}, ch1:{round((((16-0)**2)*3)**0.5,4)}") 
#-> ch0:13.8564, ch1:27.7128

print("MONAI : ")
HD = compute_hausdorff_distance(arr1,arr2,include_background=True)
print(f"ch0:{HD[0][0]:.4f}, ch1:{HD[0][1]:.4f}") 
#-> ch0:13.8564, ch1:27.7128

#変数HDの中身：
#tensor([[13.8564, 27.7128]], dtype=torch.float64)

compute_hausdorff_distanceの戻り値は(batch, class)の形となっています。
手動のユークリッド距離とMONAIで算出したHDの出力が等しく、正しく計算できていることが確認できました。
次は、2つの立方体について算出してみます。

n_classes = 1
#(batch,n_classes,X,Y,Z) array must be one-hot format
arr1 = np.zeros((1,n_classes,32,32,32),dtype=np.uint8)
arr2 = np.zeros((1,n_classes,32,32,32),dtype=np.uint8)

#２つの立方体についてHDを計算
arr1[0,0,5:15,5:15,5:15] = 1
arr2[0,0,20:30,20:30,20:30] = 1

図示するとこんな感じです。

3次元プロットの関数と実行

def plot_3d(arr1,arr2):
    
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    
    verts, faces,_,_ = measure.marching_cubes_lewiner(arr1)
    mesh = Poly3DCollection(verts[faces], alpha=0.2)
    mesh.set_facecolor("navy")
    ax.add_collection3d(mesh)
    
    verts, faces,_,_ = measure.marching_cubes_lewiner(arr2)
    mesh = Poly3DCollection(verts[faces], alpha=0.2)
    mesh.set_facecolor("red")
    ax.add_collection3d(mesh)

    ax.set_xlim(0, arr1.shape[0])
    ax.set_ylim(0, arr1.shape[1])
    ax.set_zlim(0, arr1.shape[2])

    plt.show()

plot_3d(arr1[0,0,:,:,:],arr2[0,0,:,:,:])

HD = compute_hausdorff_distance(arr1,arr2,include_background=True)
print(f"{HD[0][0]:.4f}") 
#->25.9808

立体に対しても、HDを計算することができました。ちなみに、このHDを手計算する場合、青いオブジェクトの(5,5,5)の位置から赤いオブジェクトの(20,20,20)の位置へのユークリッド距離ということになります。(今回の例は立方体なので、双方からの距離の最大値が等しく、赤の(30,30,30)から青の(15,15,15)でも同じです。)

最後に、外れ値を含めた配列を作成し、これに対する95%HDを算出して、外れ値を考慮したHDが算出できることを確認しましょう。

n_classes = 1
#(batch,n_classes,X,Y,Z) array must be one-hot format
arr1 = np.zeros((1,n_classes,32,32,32),dtype=np.uint8)
arr2 = np.zeros((1,n_classes,32,32,32),dtype=np.uint8)

arr1[0,0,5:15,5:15,5:15] = 1
arr1[0,0,1:3,1:3,1:3] = 1 #outlier 外れ値 
arr2[0,0,20:30,20:30,20:30] = 1

plot_3d(arr1[0,0,:,:,:],arr2[0,0,:,:,:])

HD_100 = compute_hausdorff_distance(arr1,arr2,include_background=True)
HD_95 = compute_hausdorff_distance(arr1,arr2,include_background=True,percentile=95)

#95% HD excludes the effect of outliers
#95%では外れ値の影響を除外できていることが分かる
print(f"HD 100% : {HD_100[0][0]:.4f}") 
#-> HD 100% : 32.9090
print(f"HD 95% : {HD_95[0][0]:.4f}") 
#-> HD 95% : 24.2899

100%HDでは小さい外れ値のオブジェクトからの距離がHDとして採用されています。
セグメンテーションではこのような小さなノイズが出現することがあるため、これでは評価指標として適切ではない場合があります。
そこで、95%のHDを採用してやることで、先の外れ値の存在しなかった場合のHD 25.9808と近い、24.2899という値を得ることができました。

実装コード

参考文献

[1] ハウスドルフ距離 (Hausdorff distance)
https://zellij.hatenablog.com/entry/20111206/p1
[2] Improving automatic delineation for head and neck organs at risk byDeep Learning Contouring
https://www.sciencedirect.com/science/article/pii/S0167814019331111?via%3Dihub
[3] Boundary loss for highly unbalanced segmentation
https://arxiv.org/abs/1812.07032