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

조회 수: 18 (최근 30일)
yuya nakashima
yuya nakashima 2019년 7월 24일
댓글: yuya nakashima 2019년 9월 3일
現在、以下のソースを用いて分布の表示をしているのですが、適当に範囲を設定してその中の最大値や最小値、
また放射線治療で用いられる線量体積ヒストグラム(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

채택된 답변

Shunichi Kusano
Shunichi Kusano 2019년 7월 25일
線量体積ヒストグラムについて知らなかったので調べてみたのですが、
問題とする領域の中で、被曝線量が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
yuya nakashima
yuya nakashima 2019년 9월 3일
私の思っていたとおりにうまくいきました。 ありがとうございます。

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 確率分布에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!