三次元グレースケール​画像内で、指定した範​囲内のエントロピーの​計算の仕方

三次元行列を作成してそれを可視化した三次元イメージ画像があるのですが、自分で立体空間を指定して（100×100×100のイメージ内で1×1×1の範囲など）そのなかのエントロピー（黒が集まってれば集まってるほど高い）を計算したいです。 それか、範囲内の密度（黒であればあるほど高い）の計算の仕方なども教えてもらいたいです。

例として作成した三次元イメージを添付します。100×100×100doubleです

Answer 1

[michio]

グレースケールのエントロピーは entropy 関数. 範囲を指定しての局所的なエントロピーは entropyfilt 関数 で求められそうです。どちらも、Image Processing Toolbox の関数です。

entropyfilt 関数はデフォルトで、9行9列近傍のエントロピー値を計算しますが、3x3x3ですと、

 entropyfilt(I3D, true(3,3,3));

となります。(例えば I3D は 100x100x100 の配列とします) ただ、同じ値が集まっていればいるほど低い値となりますが、用途に合っていますか？

また、３次元画像の事前指定した範囲内の平均値を求める関数として imboxfilt3 関数 があります。グレースケール画像ですと局所的な平均輝度（白ほど高い値）が求まりますが、ryoya suzuki さんのおっしゃる密度の定義とは逆ですが、算出に応用できるかと思うので紹介しますね。

[thuru]

なるほど、エントロピーの関数は三次元でも利用できるんですね！ これは適用する範囲は指定できないのでしょうか？例えば立方体を4分割して、それぞれの範囲内でのエントロピーを計算できたらいいのですが、、 同じ値が集まっていればいるほど、というのは例え高い数値が集まっていようと低い数値が集まっているのと同じになってしまうということでしょうか？

平均値を求めるのは、そのまま使えそうですね！ありがとうございます

[michio]

上で紹介した関数では、分割してそれぞれのエントロピーを求めるのはできそうにないですね。明示的に４つに分割（手動で）して関数に与えるのがよいかと思います。

 entropy(I3D(1:25,:,:))
 entropy(I3D(26:50,:,:))
 entropy(I3D(51:75,:,:))
 entropy(I3D(76:100,:,:))

など、スクリプトを書いて対応できます。

>> 同じ値が集まっていればいるほど、というのは例え高い数値が集まっていようと低い数値が集まっているのと同じになってしまうということでしょうか？

関数が求めるものはおっしゃるとおりです。自身は、画像の専門家ではないですし、もし異なるエントロピーの定義の仕方が紹介されている参考資料などありましたら教えてください。対応する関数を探すことができるかもしれません。

[thuru]

なるほど、でしたらエントロピーを求めるよりかは分割して平均値をそれぞれ求めていったほうが分かりやすいですし実用的ですね。 平均値を分割して求めるやり方は上記のエントロピーの式のようにできるのでしょうか。

[michio]

そうですね。平均値を取るだけなら mean 関数で大丈夫です。 実際に試してみてうまくいかないようでしたら、実行内容と結果・望ましい結果を添えてまた投稿ください。

[thuru]

ありがとうございます。mean関数についてですが、調べたのですが三次元空間内で範囲を指定して平均値を出すというやり方が分からなくて詰まっている状態です。

[michio]

mean 関数自体には範囲を指定して～といった機能はないので、与える入力を工夫します。_I3D_ の一部分を取り出し、 mean 関数に与えるなどの方法はいかがでしょうか？

 subI = I3D(1:25,:,:); 
 mean(subI(:))

など。

[thuru]

こういったことが思いつかないのがだめですね、、ほんとにありがとうございます！！ 早速試してみたいと思います。 ちなみに上記のsubIの中のすべての和（積分？）の関数はintegralで可能でしょうか。

[michio]

和は sum でOKですよ。

[thuru]

試してみます。ありがとうございます。

[thuru]

clear all; close all; imtool close all;

A = (1:100)'; B = repmat(A,1,100);

I = mat2gray(B);

X = repmat(B,[1 1 100]);

subX1 = X(1:33,1:33,1:33);

subX2 = X(34:66,1:33,1:33);

subX3 = X(67:100,1:33,1:33);

subX4 = X(1:33,34:66,1:33);

subX5 = X(34:66,34:66,1:33);

subX6 = X(67:100,34:66,1:33);

subX7 = X(34:66,67:100,1:33);

subX8 = X(34:66,67:100,1:33);

subX9 = X(67:100,67:100,1:33);

subX10 = X(1:33,1:33,34:66);

subX11 = X(34:66,1:33,34:66);

subX12 = X(67:100,1:33,34:66);

subX13 = X(1:33,34:66,34:66);

subX14 = X(34:66,34:66,34:66);

subX15 = X(67:100,34:66,34:66);

subX16 = X(1:33,67:100,34:66);

subX17 = X(34:66,67:100,34:66);

subX18 = X(67:100,67:100,34:66);

subX19 = X(1:33,1:33,67:100);

subX20 = X(34:66,1:33,67:100);

subX21 = X(67:100,1:33,67:100);

subX22 = X(1:33,34:66,67:100);

subX23 = X(34:66,34:66,67:100);

subX24 = X(67:100,34:66,67:100);

subX25 = X(1:33,67:100,67:100);

subX26 = X(34:66,67:100,67:100);

subX27 = X(67:100,67:100,67:100);

A1=entropy(subX1(:)); A2=entropy(subX2(:)); A3=entropy(subX3(:)); A4=entropy(subX4(:)); A5=entropy(subX5(:)); A6=entropy(subX6(:)); A7=entropy(subX7(:)); A8=entropy(subX8(:)); A9=entropy(subX9(:)); A10=entropy(subX10(:)); A11=entropy(subX11(:)); A12=entropy(subX12(:)); A13=entropy(subX13(:)); A14=entropy(subX14(:)); A15=entropy(subX15(:)); A16=entropy(subX16(:)); A17=entropy(subX17(:)); A18=entropy(subX18(:)); A19=entropy(subX19(:)); A20=entropy(subX20(:)); A21=entropy(subX21(:)); A22=entropy(subX22(:)); A23=entropy(subX23(:)); A24=entropy(subX24(:)); A25=entropy(subX25(:)); A26=entropy(subX26(:)); A27=entropy(subX27(:));

上記のコードで計算すると、エントロピーの値がすべて０になります。どこかおかしいのでしょうか。

[michio]

 X = repmat(B,[1 1 100]);

を

 X = repmat(I,[1 1 100]);

に変更してみてください。

関数の詳細になりますが、 entropy関数 は入力変数として受け付けるデータタイプが logical, uint8, uint16, double 型とドキュメンテーションに記述されてます。

それぞれのデータは内部で im2uint8 と呼ばれる関数で uint8型へ キャストされ、 256階調でヒストグラムが作成されエントロピーが計算されます。 すなわち、それぞれの階調の画素の数 p を元に

 -sum(p.*log2(p))

とエントロピーが計算されます。

ここで、 im2uint8 の仕様ですが、入力変数のデータ型によって挙動が異なります。

• uint16型の場合、0 ～ 65535 までの値を 0 ～ 255 に量子化します。
• double型の場合は、0 ～ 1 までの実数値を 0 ～ 255 までの値へ量子化されることになります。

今回の場合

 X = repmat(B,[1 1 100]);

と実行されていたことから X が double 型で定義されてしまっていたため、1 以上のデータはすべて１として解釈され、結果すべての画素が 1 である画像として認識され、エントロピーは 0 と算出されていました。

参考： http://jp.mathworks.com/matlabcentral/answers/164622-estimation-measure-of-entropy

[thuru]

全ての値が0となった理由が分かりました。ありがとうございます。 Iで実行したところ、0～1までの値になってしまってたのですが、Iを256階調にして256でエントロピーの計算を行いたいのですが、uint8にすればいいのでしょうか？

[michio]

すいません、「Iを256階調にして256でエントロピーの計算を行いたい」とはどういうことでしょうか。もう少し詳しく教えてください。 現在の結果では要求を満たさない等、困る点を交えて頂けますと対策を考えることができるかもしれません。

[thuru]

今回自分が用いている三次元行列は、シミュレーションのために作成したものなのですが、今後エントロピーなどを計算したいと考えている画像はグレースケールの256階調のものなのです。ですからできればシミュレーションのときも0から255の中で計算したいのです。 上記のIは0から1の値なので0から255にしたいです。 説明下手ですみません。

[michio]

なるほどそうこうことですね。ありがとうございます。それには im2uint8 関数 が使えます。

 I = im2uint8(mat2gray(B));

といった形でしょうか。I は 0-255で定義されるグレースケール画像になります。

[thuru]

できました！ありがとうございます。

[michio]

:)