範囲の中から最大値や最小値を表示したいです。

18 views (last 30 days)
yuya nakashima
yuya nakashima on 24 Jul 2019
Commented: yuya nakashima on 3 Sep 2019
現在、以下のソースを用いて分布の表示をしているのですが、適当に範囲を設定してその中の最大値や最小値、
また放射線治療で用いられる線量体積ヒストグラム(Dose Volume Histgram)のようなものを表示したいと考えているのですが、
どうすればよろしいでしょうか。
clear ;close all;clc
warning off;
[Xtest,Ytest,Ztest] = meshgrid(-9:1:9); %Xtest,Ytest,Ztest座標に-9~9まで1刻みの箱を作成
V=zeros(19,19,19);
p=-0.0444*(sqrt(Xtest(:,:,4).^2+Ytest(:,:,4).^2)).^2 + 0.1857*sqrt(Xtest(:,:,4).^2+Ytest(:,:,4).^2) + 3.55;
for a=10:14;
p=p-0.5;
V(:,:,a)=p;
end
for b=6:10;
p=p+0.5;
V(:,:,b)=p;
end
xslice = [0];
yslice = [-9,0,9];
zslice = [0];
h = slice(Xtest,Ytest,Ztest,V,xslice,[],zslice);
xlabel('X');ylabel('Y');zlabel('Z');
caxis([min(V(:)) max(V(:))])
axis equal
set(h,'EdgeColor','none','FaceColor','interp',...
'FaceAlpha','interp')
alpha('color')
alphamap('increase',.3)
rotate3d

Accepted Answer

Shunichi Kusano
Shunichi Kusano on 25 Jul 2019
線量体積ヒストグラムについて知らなかったので調べてみたのですが、
問題とする領域の中で、被曝線量がx以上である領域の比率f(x)が線量体積ヒストグラムという認識でいいでしょうか。
ご提示のプログラムでいうところのVが、各ボクセルにおける被曝線量だと仮定しますと、問題とする領域の中での被ばく線量の(1-累積分布関数)が、線量体積ヒストグラムに相当するのかなと思います(累積分布関数はx"以下"である比率)。
という妄想の仮定のもとですが、下記のような形で計算できるのかなと思います。
%% 適当に範囲を選択。抜き出す範囲は図などを見ながらピックアップする必要があります。
% 2019aであればdrawcuboidという関数が便利かと思います。
Vpart = V(2:8,3:5,5:9); % とりあえず適当に範囲を指定しました。
% 最小値
min(Vpart(:)) % (:)が重要で、値を1次元にしています。つまりVpartに含まれる全ての値を対象に最小値を抽出する、という意味になります。
% 最大値
max(Vpart(:))
%% ヒストグラムから累積分布関数を求めて1から引く
[N, edges] = histcounts(Vpart(:), [-1.5:0.1:3]);
cumN = cumsum(N); % 累積カウント
sumN = sum(N); % 全voxel数
dvh = (sumN - cumN) / sumN; % 1 - 累積分布関数
figure, plot(edges(2:end), dvh, 'b.-');
%% statistics and machine learning toolboxをお持ちの場合、関数から直接計算可能
% データの累積分布関数の'function'を'survivor'に設定すると線量体積ヒストグラム
[dvh, x] = ecdf(Vpart(:),'function', 'survivor');
hold on;
plot(x, dvh, 'r.-');
私の理解が誤っているかもしれません。ご教示いただければ、それに応じて修正したプログラムに書き換えますので、コメントいただけましたら幸いです。
  1 Comment
yuya nakashima
yuya nakashima on 3 Sep 2019
私の思っていたとおりにうまくいきました。 ありがとうございます。

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!