Radiomics超入門：Neighbourhood gray tone difference matrix特徴

2023年9月24日 10:22

概要

[Amadasun1989]は、グレーレベル共起行列（GLCM）に代わるNeighbourhood gray tone difference matrix（近傍グレートーン差行列, NGTDM）を提案しました。

NGTDMは、離散化されたグレーレベル$${i}$$を持つボクセルのグレーレベルとチェビシェフ距離$${δ}$$内の近傍ボクセルの平均グレーレベルの差の合計を考慮します。

3D については、AmadasunとKingによるオリジナルの定義を拡張することができます。位置$${k=(kx,ky,kz)}$$におけるボクセルの離散化済みのグレーレベルを$${X_{d,k}}$$とします。

そして、$${(kx,ky,kz)}$$を中心として、$${(kx,ky,kz)}$$自身を除いた近傍領域内の平均グレーレベル$${\overline{X_k}}$$を求めます。

$$
\overline{X_k}=\frac{1}{W} \displaystyle\sum_{m_z=-\delta}^{\delta} \displaystyle\sum_{m_y=-\delta}^{\delta} \displaystyle\sum_{m_x=-\delta}^{\delta} X_d (k_x+m_x, k_y+m_y, k_z+m_z) {, }(m_x,m_y,m_z)\not=(0,0,0)
$$

$${m}$$は、$${k}$$からみたときの近傍ボクセル位置への方向ベクトルです。

$${W=(2δ+1)^3-1}$$は 3D 近傍領域のボクセル数（中心ボクセルを除く）です。2D の場合、$${W=(2δ+1)^2-1}$$になり、異なるスライス間の平均は計算されません。

グレーレベル$${i}$$の近傍グレートーン差（Neighbourhood gray tone difference）$${s_i}$$は次のように求められます。

$$
s_i=\displaystyle\sum_{k}^{N_v} \lvert {i-\overline{X_k}} \lvert { }[X_d(\bold k) = \text{i and k has a valid neighbourhood}]
$$

ここで、[...]はアイバーソン条件式であり、ボクセル $${k}$$ のグレーレベル $${X_{d,k}}$$ が $${i}$$ に等しく、かつ、そのボクセルが有効な近傍ボクセルをもつという条件が共に真であれば 1 （つまり、計算する）となり、そうでなければ 0 となります。$${N_v}$$ は ROI 信号強度マスクのボクセル数です。

2D の例を下図に示します。

2DにおけるNGTDMの例
(a) 元画像（離散化済みグレーレベル）とマスク領域
(b) NGTDM
Addapted from IBSI ref Fig.18

この例では、距離$${δ=1}$$、近傍は8近傍を用いています。まず、NGTDMのグレーレベル$${i=1}$$に着目してみます。$${n_1=0}$$となっています。これは、中心画素となるグレーレベル「1」の有効な画素が(a)の矩形マスク内にないため、「0」です。これに伴い、近傍グレートーン差$${s_1}$$も 0.0 です。

次に、グレーレベル$${i=2}$$は、マスク内に2つあります。よって、$${n_2=2}$$です。それぞれの近傍画素の平均値は (1+2+2+1+3+4+2+4)/8 と (1+2+3+4+4+4+1+2)/8です。したがって、$${ s_2=|2-19/8|+|2-21/8|=1}$$となります。ここで、この例では、近傍画素の平均を計算する際に、マスク内の画素に限定していないことに注意が必要です。同様に、$${ s3=|3-19/8|=0.625}$$、$${ s4=|4-17/8|=1.825}$$となります。

この例では、[Amadasun1989]によるオリジナルの定義に倣った計算を行っています。しかし、IBSIの定義はオリジナルの定義とは異なる定義を用いています。違いは、有効な近傍領域に対する考え方で、近傍領域画素が、ROI マスク内に存在する場合にのみ計算するという制約を設けています。

よって、先程の式を次のように書き換えて計算します。

$$
\overline{X_k}=\frac{1}{W_k} \displaystyle\sum_{m_z=-\delta}^{\delta} \sum_{m_y=-\delta}^{\delta} \sum_{m_x=-\delta}^{\delta} X_d (k+m) [ m \not=0{ and k+m in ROI}]
$$

近傍のボクセルをカウントする際にも、ROIの中にあるかを条件に加えます。

$$
{W_k}=\displaystyle\sum_{m_z=-\delta}^{\delta} \sum_{m_y=-\delta}^{\delta} \sum_{m_x=-\delta}^{\delta} [ m \not=0 \text{ and } k+m \text{ in ROI}]
$$

近傍グレートーン差は、近傍のボクセルが検出されたときに加算されます。

$$
s_i=\displaystyle\sum_{k}^{N_v} \lvert {i-\overline{X_k}} \lvert { }[X_d(k) =i \text{ and } W_k \not=0]
$$

多くの NGTDM 特徴は、ROI 信号強度マスク内のグレーレベルの数$${N_g}$$と、確率分布$${p_i=n_i/N_{v,c}}$$に依存しています。ここで、$${N_{v,c}=\sum n_i}$$は少なくとも1つの近傍を持つボクセルの総数です。

全てのボクセルが少なくとも1つの近傍を持つ場合、$${N_{v,c}=N_v}$$となります。さらに、$${N_{g,p}≦N_g}は、$${p_i＞0}$$であるグレーレベルの数を意味します。上の例では、$${N_g=4}$$、$${N_{g,p}=3}$$です。

特徴の集約

GLSZMと同様の考え方です。

距離による重みづけ

デフォルトの近傍領域は，チェビシェフ・ノルムを用いて定義されます。これに限るわけではなく、マンハッタン・ノルムやユークリッド・ノルムを使うこともできます。この場合、平均グレーレベル$${\overline{X_k}}$$をより一般的に定義する必要があります。

$$
\overline{X_k}=\frac{1}{W_k} \displaystyle\sum_{m_z=-\delta}^{\delta} \sum_{m_y=-\delta}^{\delta} \sum_{m_x=-\delta}^{\delta} X_d (k+m) [\| m \| \le \delta{ and }m\not=0 \text{ and }k+m \text{ in ROI}]
$$

$$
W_k=\displaystyle\sum_{m_z=-\delta}^{\delta} \displaystyle\sum_{m_y=-\delta}^{\delta} \displaystyle\sum_{m_x=-\delta}^{\delta} [\| m \| \le \delta{ and }m\not=0 \text{ and } k+m \text{ in ROI}]
$$

NGTDM の距離による重み付けは比較的簡単です。$${w}$$ を中心ボクセルから見た時の位置ベクトルの長さ$${m}$$ に依存する重みとして、例えば、$${ w=\| m\|^{-1}}$$ または $${w=exp(-\| m \| ^2)}$$ とします。すると、平均グレーレベルは次のように計算されます。

$$
\overline{X_k}=\frac{1}{W_k} \displaystyle\sum_{m_z=-\delta}^{\delta} \sum_{m_y=-\delta}^{\delta} \sum_{m_x=-\delta}^{\delta} w(m)X_d (k+m) [\| m \| \le \delta{ and }m\not=0 \text{ and }k+m \text{ in ROI}]
$$

$$
W_k=\displaystyle\sum_{m_z=-\delta}^{\delta} \displaystyle\sum_{m_y=-\delta}^{\delta} \displaystyle\sum_{m_x=-\delta}^{\delta} w(m)[\| m \| \le \delta{ and }m\not=0 \text{ and } k+m \text{ in ROI}]
$$

なお、IBSI では距離による重み付けは非推奨となっています。

NGTDM特徴

Coarseness

一般に、大きいスケールのパターンでは、粗い（大まかな）テクスチャのグレーレベルの差は小さくなります。差を合計すると、信号強度の空間的な変化率のレベルを示すことができます [Amadasun1989]。Coarseness（粗さ）は次のように定義されます。

$$
F_{ngt.coarseness}=\frac{1}{\sum_{i=1}^{N_g} {p_i}{s_i}}
$$

$${\sum_{i=1}^{N_g} {p_i}{s_i}}$$は 0 に近くなることがあるため、特徴量の最大値は任意の数、例えば、$${10^6}$$が設定されます。AmadasunとKingは元々、分母に不特定の小さな数$${\epsilon}$$を代入することでこの問題を回避していましたが、明示的に最大値を設定することで、より一貫性を持たせることができる可能性があります。

Contrast

Contrastは、グレーレベルのダイナミックレンジと信号強度変化の空間周波数に依存する特徴量です。次のように定義されます。

$$
F_{ngt.contrast} = \displaystyle \big ( \frac {1}{N_{g,p}(N_{g,p}-1)} \sum_{i_1=1}^{N_g}\sum_{i_2=1}^{N_g}p_{i_1}p_{i_2}(i_1-i_2)^2 \big) \big ( \frac {1}{N_{v,c}} \sum_{i=1}^{N_g} s_i \big )
$$

グレーレベル確率$${p_{i1}}$$と$${p_{i2}}$$は、$${p_i}$$のコピーで、それぞれ異なる繰り返しインデックスを持つことのみが違いです。式の第 1 項はグレーレベルのダイナミックレンジを考慮し、第 2 項はボリューム内の信号強度変化の尺度です。

$${N_{g,p}=1}$$のとき、$${F_{ngt.contrast}=0}$$です。

Busyness

隣接するボクセル間のグレーレベルの変化が大きいテクスチャは、"ビジー"であるとされています。Busynessは次のように定義されます。

$$
F_{ngt.busyness} = \frac {\sum_{i=1}^{N_g}p_i s_i} {\sum_{i_1=1}^{N_g}\sum_{i_2=1}^{N_g} i_1p_{i_1} - i_2p_{i_2}},{ }p_{i_1} \not=0 \text{ and }p_{i_2} \not=0
$$

しかし、この元の定義は、分母は常に0と評価されるという誤った定式化になっています。そこで、少し異なる定義 [Hatt2016]を用います。

$$
F_{ngt.busyness} = \frac {\sum_{i=1}^{N_g}p_i s_i} {\sum_{i_1=1}^{N_g}\sum_{i_2=1}^{N_g} | i_1p_{i_1} - i_2p_{i_2}|},{ }p_{i_1} \not=0 \text{ and }p_{i_2} \not=0
$$

$${N_{g,p}=1}$$のとき、$${F_{ngt.busyness}=0}$$です。

Complexity

複雑なテクスチャは一様ではなく、グレーレベルが大きく変化するのが一般的です。この複雑さを評価する指標がComplexityです。

$$
F_{ntg.complexity} = \displaystyle\frac {1} {N_{v,c}} \sum_{i_1=1}^{N_g}\sum_{i_2=1}^{N_g} |i_1 - i_2| \frac{p_{i_1}s_{i_1} + p_{i_2}s_{i_2}}{p_{i_1}+p_{i_2}},{ }p_{i_1} \not=0 \text{ and }p_{i_2} \not=0
$$

Strength

テクスチャの強さは次のように定義されています。

$$
F_{ngt.strength} = \frac {\sum_{i_1=1}^{N_g}\sum_{i_2=1}^{N_g} (p_{i_1}+p_{i_2})(i_1 - i_2)^2} {\sum_{i=1}^{N_g}s_i},{ }p_{i_1} \not=0 \text{ and }p_{i_2} \not=0
$$

もし、$${\sum_{i=1}^{N_g}s_i=0}$$のとき、$${F_{ngt.strength}=0}$$です。

実践

RadiomicsJを用いて実践していきます。

NGTDMの計算

上記の 2D サンプル例を用いて、NGTDMを出力してみます。

byte pixels[] = new byte[16];
byte r0[] = new byte[] { 1, 2, 2, 3 };
byte r1[] = new byte[] { 1, 2, 3, 3 };
byte r2[] = new byte[] { 4, 2, 4, 1 };
byte r3[] = new byte[] { 4, 1, 2, 3 };
// flatten to create ByteProcessor
int i= 0;
for(byte[] r: new byte[][] {r0,r1,r2,r3}) {
	for(byte v:r) {
		pixels[i++] = v;
	}
}
ImageProcessor bp = new ByteProcessor(4, 4, pixels);
ImagePlus imp = new ImagePlus("sample", bp);

byte roi_mask[] = new byte[16];
byte m0[] = new byte[] { 0, 0, 0, 0 };
byte m1[] = new byte[] { 0, 1, 1, 0 };
byte m2[] = new byte[] { 0, 1, 1, 0 };
byte m3[] = new byte[] { 0, 0, 0, 0 };
i= 0;
for(byte[] r: new byte[][] {m0,m1,m2,m3}) {
	for(byte v:r) {
		roi_mask[i++] = v;
	}
}
ImageProcessor bp2 = new ByteProcessor(4, 4, roi_mask);
ImagePlus mask = new ImagePlus("sample_mask", bp2);
		
RadiomicsJ.targetLabel = 1;
int delta=1;
NGTDMFeatures test = new NGTDMFeatures(imp/*descretized*/,mask,RadiomicsJ.targetLabel,delta);
// re-calculate
test.fillNGTDM(true);//amadasunAlgorithms
System.out.println(test.toString());

//出力
i | Ni | Pi | Si
1.0 0.0 0.0 0.0
2.0 2.0 0.5 1.0
3.0 1.0 0.25 0.625
4.0 1.0 0.25 1.875

特徴の計算

RadiomicsJを用いて、IBSIデジタルファントムから特徴を計算してみます。$${\delta = 1}$$です。

ImagePlus ds_pair[] = TestDataLoader.digital_phantom1();
ImagePlus imp = ds_pair[0];
ImagePlus mask = ds_pair[1];

RadiomicsJ.targetLabel = 1;

//ImagePlus img, ImagePlus mask, int label, Integer delta, boolean useBinCount, Integer nBins, Double binWidth
NGTDMFeatures f = new NGTDMFeatures(imp, mask ,RadiomicsJ.targetLabel, 1, true, 6 ,null);
System.out.println(f.toString());
		
System.out.println(NGTDMFeatureType.Coarseness+":"+f.calculate(NGTDMFeatureType.Coarseness.id()));//OK
System.out.println(NGTDMFeatureType.Contrast+":"+f.calculate(NGTDMFeatureType.Contrast.id()));//OK
System.out.println(NGTDMFeatureType.Busyness+":"+f.calculate(NGTDMFeatureType.Busyness.id()));//OK
System.out.println(NGTDMFeatureType.Complexity+":"+f.calculate(NGTDMFeatureType.Complexity.id()));//OK
System.out.println(NGTDMFeatureType.Strength+":"+f.calculate(NGTDMFeatureType.Strength.id()));//OK

//出力
i	Ni	Pi	Si
1	50	0.676	39.947
2	0	0.000	0.000
3	1	0.014	0.200
4	16	0.216	20.825
5	0	0.000	0.000
6	7	0.095	24.127

Coarseness:0.029604225700771966
Contrast:0.5837109155843375
Busyness:6.543568456913214
Complexity:13.539763601407849
Strength:0.7634954330521357

RadiomicsJの引用はこちら

Kobayashi, T. RadiomicsJ: a library to compute radiomic features. Radiol Phys Technol 15, 255–263 (2022). https://doi.org/10.1007/s12194-022-00664-4

RadiomicsJのリンク

https://github.com/tatsunidas/RadiomicsJ

Stay visionary

この記事が気に入ったらサポートをしてみませんか？