クロフトンの公式でボリュームデータの表面積を計算する

はじめに

この記事では、CTスキャンなどのボリュームデータから、セグメンテーションされたラベルデータの表面積を正確に計算するための「クロフトンの公式」について解説します。3次元データの表面積計算は、医用画像処理をはじめとする様々な分野で重要となります。本記事では、クロフトンの公式の理論的背景から、Python/NumPyを用いた実装、さらには解像度と精度の関係までを詳しく説明します。

この記事を書いたひと

デジタルリアクタ合同会社 代表
望月 優輝 / Yuki Mochizuki   

機械学習・統計、数値計算などの領域を軸としたソフトウェアエンジニアリングを専門としています。スタートアップからグローバル企業まで、さまざまなスケールの企業にて、事業価値に直結する計算システムを設計・構築してきました。主に機械学習の応用分野において、出版・発表、特許取得等の実績があり、また、IT戦略やデータベース技術等の情報処理に関する専門領域の国家資格を複数保有しています。九州大学理学部卒業。日本ITストラテジスト協会正会員。

対象読者:

• 画像処理、特に医用画像処理に携わる研究者やエンジニア
• 3次元データの表面積計算に関心のある方
• Python/NumPyを用いた数値計算に興味のある方

記事のポイント:

• クロフトンの公式の概要と、幾何学的確率論におけるその位置づけ
• 2次元での直感的な理解から、3次元への拡張方法
• Python/NumPyによる具体的な実装例と、高速化のためのテクニック
• 解像度と計算精度の関係性に関する評価

クロフトンの公式とは

クロフトンの公式は、ボリュームデータ（例えばCTスキャン）内のセグメンテーションされたラベルデータの表面積を計算する際に非常に有効な手法です。この公式を利用することで、マーチングキューブ法などの他の手法と比較して、より正確な表面積を算出することが可能となります。

クロフトンの公式は、幾何学的確率論における重要な定理であり、物体の幾何学的性質（長さ、面積、体積など）を、ランダムな直線や平面との交差回数と関連付けます。

2次元での直感的な理解

3次元で考えると複雑になるため、まずは2次元の場合で、クロフトンの公式を直感的に理解してみましょう。例として、単位円に対して、X軸に平行な線を等間隔に多数交差させる状況を考えます。

同様に、他の角度からも線を引いていきます。図形が複雑な形状であるほど、または周長が長くなるほど、交差回数が増加することが予想されます。重要なのは、各方向から引いた線と図形の交差回数の平均が、図形の周長に比例するという点です。この性質は円に限らず、任意の2次元図形に適用でき、3次元の表面積計算にも拡張できます。

クロフトンの公式の詳細は、Wikipediaの解説を参照することをお勧めします。

3次元への拡張

2次元での考え方を3次元に拡張します。理論的には、無限の方向から平行な線を生成する必要がありますが、離散的に表現されたデータの表面積を計算する際には、特定の方向からの交差回数を数えることで近似できます。

具体的には、2次元グリッドデータの場合は4方向（縦、横、斜め2方向）、3次元グリッドデータの場合は13方向からの交差を考えます。これらの方向は、それぞれ8近傍、26近傍を対称性から2で割ったものです。

Python/NumPyによる実装

このセクションでは、PythonとNumPyを用いて、クロフトンの公式を実装し、3次元データの表面積を計算する具体的な方法を説明します。

まず、サンプルデータとして、3次元のNumPy配列を生成し、半球内部を1でマークします。

import numpy as np
import matplotlib.pyplot as plt

data = np.zeros((101,101,101), dtype=np.uint8)

for x in range(data.shape[0]):
    for y in range(data.shape[1]):
        for z in range(data.shape[2]):
            if (x-50)**2+(y-50)**2+(z-50)**2 < 50**2 and z < 50.5:
                data[x,y,z] = 1

plt.imshow(data[:, 50, :])
plt.show()

このコードは、中心が(50, 50, 50)で半径が50の半球をdata配列に作成し、y=50での断面を表示します。

次に、クロフトンの公式に基づき、13方向の交差回数をカウントして表面積を計算します。交差回数をカウントする際には、方向に応じて、1, \sqrt{2}, \sqrt{3} の重みを付けます。これらの重みは、交差方向における交差点の平均間隔の比率を表しています。ボクセルサイズが異方性である（XYZでボクセルサイズが異なる）場合、これらの係数は調整する必要があります。

参考: ボクセルサイズが異なる場合の重み付け

ボクセルサイズが異なる場合の重み付けについて、参考まで。

等方的なボクセルサイズの場合、カウントは1, \sqrt{2}, \sqrt{3}で割られます。これは、各交差方向における交差点の平均間隔の比率に対応します。ボクセルサイズが異なる場合、例えば各辺の長さが1, 2, 3である場合、方向に応じて1, 2, 3, \sqrt{1^2+2^2}, \sqrt{1^2+3^2}, \sqrt{3^2+2^2}, \sqrt{1^2+2^2+3^2} で割ることで、比率を適切に表現できます。

さらに、ボクセルサイズが大きくなると、1つの交差が表す表面積も増加します。そのため、上記の比率に各辺の長さの積（この例では 1 \times 2 \times 3 = 6）を掛けます。これにより、次元が（長さ^-1）*（体積）=（面積）となり、物理的な整合性が保たれます。

異方性ボクセルデータを前提とした実装方法の詳細については、"On the Analysis of Spatial Binary Images"を参照するとよいでしょう。今回は計算が複雑になるので、ボクセルサイズが(1, 1, 1)の場合として進めます。

NumPyによる高速化

計算を高速化するために、斜め方向にグリッド値を一つずつ取得するのではなく、np.roll関数を用いてまとめてずらす方法を採用します。np.rollは、配列の要素を指定した軸方向にシフトさせる関数です。これにより、NumPyの軸に沿った並列化された高速な演算を活用できます。np.rollによって画素が斜めにずれる様子を以下に示します。

np.rollでずらした際の境界条件の影響を避けるため、事前にnp.pad関数を用いて周囲に0のパディングを追加します。

padded = np.pad(data, 1, 'constant', constant_values=0)

np.padは、配列の周囲に指定した幅で定数値を埋める関数です。

軸方向に沿って0/1の反転回数を計算するには、np.diff関数を使用します。np.diffは、隣接する要素間の差分を計算する関数で、0→1の反転では+1、1→0の反転では-1を返します。

diff = np.abs(np.sign(np.diff(rolled, axis=axis1)))

この結果の絶対値を取り、合計することで、指定された軸方向の交差回数が求まります。

cross_section = np.sum(diff) * (2**-0.5)

実装例

上記のテクニックを組み合わせ、クロフトンの公式を実装した完全なコード例を以下に示します。

from copy import deepcopy
import numpy as np

data = np.zeros((101,101,101), dtype=np.uint8)
for x in range(data.shape[0]):
    for y in range(data.shape[1]):
        for z in range(data.shape[2]):
            if (x-50)**2+(y-50)**2+(z-50)**2 < 50**2 and z < 50.5:
                data[x,y,z] = 1

def crofton_surface(data):
    padded = np.pad(data, 1, 'constant', constant_values=0)

    cross_section_sum = 0

    # 軸 X, Y, Z に沿った交差
    for axis in [0, 1, 2]:
        diff = np.abs(np.sign(np.diff(padded,axis=axis)))
        cross_section = np.sum(diff)
        cross_section_sum += cross_section

    # 軸 XY, YZ, ZX, -XY, -YZ, -ZX に沿った交差
    for axis1 in [0, 1, 2]:
        for axis2 in [0, 1, 2]:
            if axis1 >= axis2:
                continue
            for sgn in [1, -1]:
                rolled = deepcopy(padded)
                for i in range(padded.shape[axis1]):
                    slicer = tuple([(slice(None) if j != axis1 else [i]) for j in range(3)])
                    rolled[slicer] = np.roll(rolled[slicer], i*sgn, axis=axis2)

                diff = np.abs(np.sign(np.diff(rolled, axis=axis1)))
                cross_section = np.sum(diff) * (2**-0.5)
                cross_section_sum += cross_section

    # 軸 XYZ, XY-Z, X-Y-Z, X-YZ に沿った交差
    rollsgn = [
        [1, 1],
        [1, -1],
        [-1, 1],
        [-1, -1]
    ]
    for sgns in rollsgn:
        rolled = deepcopy(padded)
        for i in range(padded.shape[0]):
            rolled[[i], : , :] = np.roll(rolled[[i], : , :], i*sgns[0], axis=1)
            rolled[[i], : , :] = np.roll(rolled[[i], : , :], i*sgns[1], axis=2)

        diff = np.abs(np.sign(np.diff(rolled, axis=0)))
        cross_section = np.sum(diff) * (3**-0.5)
        cross_section_sum += cross_section

    crofton_s = cross_section_sum / 13 * 2
    return crofton_s

print(crofton_surface(data))

このコードを実行すると、表面積として約23134.45が得られます。半球の理論値は約23561.94なので、誤差は約-1.8%です。この誤差は、主にグリッドデータによる離散化に起因します。

精度の評価

クロフトンの公式による表面積計算の精度を評価するために、全球を対象とし、解像度と理論値との誤差の関係を調べます。

for resolution in [3, 6, 12, 25, 50, 100, 200]:

    s = resolution * 2 + 1
    data = np.zeros((s,s,s), dtype=np.uint8)
    for x in range(data.shape[0]):
        for y in range(data.shape[1]):
            for z in range(data.shape[2]):
                if (x-resolution)**2+(y-resolution)**2+(z-resolution)**2 < resolution**2:
                    data[x,y,z] = 1

    ref_r = (np.sum(data) / (4 / 3 * np.pi))**(1/3)
    ref_val = ref_r * ref_r * 4 * np.pi

    crofton_val = crofton_surface(data)
    print(resolution, crofton_val / ref_val - 1)

このコードは、さまざまな解像度で球を作成し、クロフトンの公式で計算された表面積と、体積から逆算した理論値との誤差を比較します。

解像度が低い場合、ボクセルグリッド表現では球を正確に表現できないため、誤差が大きくなります。そのため、生成されたグリッド表現の体積から半径を逆算し、その半径を用いて表面積の理論値を計算しています。

半径50のグリッドデータ表現の場合、誤差は約-0.09%程度であることがわかります。十分に高い精度であるといえます。ただし、クロフトンの公式は、表面積を一貫して小さめに計算する傾向があります。これは、球表面の曲率が厳密には1/rであるのに対し、クロフトンの公式では13方向から交点を計算するため、ボクセルレベルまで拡大すると表面の曲率がゼロに近づくためです。

また、この結果から、小さなボクセル集合で形成されるボリュームに関しては、表面積の正確な計算が本質的に難しいことがわかります。例えば、CTデータにおいて、1ボクセルで観測された輝点が球なのか、立方体なのか、あるいは他の形状なのかを判断することは、情報が限られているため、そもそも難しい問題であることがわかるでしょう。

関連投稿:

Graph Neural Networksを用いた異常検知入門 確率的最適化アルゴリズムの実装と比較：遺伝的アルゴリズムと勾配降下法 因果推論入門：基本概念とPythonによる実践 モンテカルロ法における次元の呪いの緩和方法：重点サンプリングによる効率化