4次元分布における複数重心点座標の取得について

添付図のように4次元分布を作成しました。
4次元分布上に5個の丸い分布があります。
この5個の分布それぞれの重心点を取得したいです。
分布上から1点の重心座標は下記のコードで取得できると思いますが、1つの分布上から任意の個数だけ重心点を取得する方法がありましたら
教えて頂きたいです。
v = squeeze(v);
BW = true(size(v));
W = regionprops3(BW,v,'WeightedCentroid');
WC = W{1,{'WeightedCentroid'}};
1.jpg

댓글 수: 9

Kenta
Kenta 2019년 10월 18일
こんにちは、5つの丸い分布それぞれから重心を取り出す、つまり、合計5つの点を取得したい、ということで正しいでしょうか。また、そのデータを添付いただけると幸いです。
ryo tanaka
ryo tanaka 2019년 10월 18일
コメントありがとうございます。その通りです。
この分布から重心点を5個取得したいです。
4次元分布に用いているデータvは49×26×180 doubleなのですが、どのように添付すればよいのかわかりません。
どうしたらよいでしょうか?
Kenta
Kenta 2019년 10월 18일
편집: Kenta 2019년 10월 18일
ご返信ありがとうございます。
データvをmatlabファイルとして保存し、このコメント欄または質問欄のINSERTのクリップのボタンを押せば添付できます。save model vなどとコマンドを打てば、そのパスに保存することができます。
また、5つの塊をそれぞれどのように見分けるのか基準があれば望ましいです。
なにか条件などあれば教えてください。
ryo tanaka
ryo tanaka 2019년 10월 18일
ご返信ありがとうございます。
添付の仕方を教えて頂きありがとうございます。
コードは下記に記載しました。
しかし、mat fileの添付は下記の理由でできないようです。
初めてこのような文がでてきたので、よくわかりませんが、明日以降にファイルの添付をすることになるかと思います。申し訳ありません。
また、見分ける基準は各塊のピーク値(付近)に重心点が置ければいいかなと思います。
(あくまで、塊の中心ではなく塊のピーク点です。)
次の理由でこのファイルを添付できません:
  • アップロードは 1 日 10 件に制限されています。追加ファイルをアップロードする必要がある場合、今すぐ 1 つ以上のファイルを削除するか、または 24 時間が経過してからアップロードしてください。
load X.mat
load Y.mat
load Z.mat
load V.mat
xslice = [-1/2 12];
yslice = 15;
zslice = [519/2 349];
figure
slice(X, Y, Z, v, xslice, yslice, zslice) % display the slices
shading faceted
ylim([ymin-1 ymax+1])
view(-34,24)
% set(gcf, 'Color', 'none');
cb = colorbar; % create and label the colorbar
cb.Label.String = 'Temperature, C';
% 重心位置の決定
v = squeeze(v);
BW = true(size(v));
W = regionprops3(BW,v,'WeightedCentroid');
WC = W{1,{'WeightedCentroid'}};
Kenta
Kenta 2019년 10월 18일
こんにちは、コードを教えていただきありがとうございます。
本日、他に10個ほどのファイルをアップロードされたのでしょうか。
もし、いまからアップロードするものがたくさんあるのでしたら、
1つのファイルとしてまとめて、.matファイルにすればよいかもしれません。
ryo tanaka
ryo tanaka 2019년 10월 18일
返信ありがとうございます。
ファイルのアップロードは本日していません。
(というか、ファイルのアップロードはこれまでにやったことがありません。。)
1つのファイルさえもアップロードできない状況です。
ご迷惑おかけして申し訳ありません。
Kenta
Kenta 2019년 10월 18일
詳しく教えていただき、ありがとうございます。
ご自身の質問欄のINSERTの上にマウスを置けば、attachmentsと出るはずで、そこをクリックすればできると思いました。そこでエラーが出てしまうのでは少し厳しいですね。
もう少し試して難しければスタッフの方に連絡すれば良いと思います。
ryo tanaka
ryo tanaka 2019년 10월 18일
attachmentsをクリックしファイルを選択してもやはりエラーが出てきてしまいます。。
スタッフの方にご連絡したいと思います。
a.png
ryo tanaka
ryo tanaka 2019년 10월 19일
すみません、遅くなりました。
ようやくファイルを添付できましたので、お送りいたします。

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

 채택된 답변

Kenta
Kenta 2019년 10월 19일

6 개 추천

こんにちは。データを送っていただきありがとうございます。
X,Y,Z,vとありますが、X-Zは、vのボクセル座標を定めるためのグリッドで、それぞれのグリッド内の値が
そのX(やY・Z)座標に対応している、そして、vが実際にそのプロットの持つ強度である、ということでしょうか。
コメントにて、コードも添付いただきありがとうございました。申し訳ございませんが、こちらで改めてコードを書きました。
質問者様は、上のようにボクセルを定義して、各重心を求める方針だったと読み取りましたが、
下では、ポイントクラウドで処理しています。下のように書けば重心を求めるまでかなり短く書くことができます。
重心は一般的な3D上でのものを計算しています。本来はこの過程でvの強度を利用して、重みつきの重心を定義するのですかね?また、4次元、というキーワードがありますが、これを4D画像としてとらえるなら、vの値を利用せず重心を求めているのでもう少し変更する必要があるとは思います。適宜変更いただければと思います。
詳細は適宜コード内のコメントをご参照ください。よろしくお願いいたします。
clear;clc;close all
% data loading
load V.mat
% find out the points with non-zero intensity
idx=find(v~=0);
% extract the XYZ coordinate
[subx, suby, subz]=ind2sub(size(v),idx);
% acquire the intensity, it was not used, though
int = v(idx);
% convert the xyz into pointcloud variable
ptCloud=pointCloud([subx, suby, subz]);
% figure;pcshow(ptCloud,'MarkerSize',50)
% closest pair of points whose euclidean distance is less than "minDistance"
% was assigned into the same cluster
minDistance = 3;
[labels,numClusters] = pcsegdist(ptCloud,minDistance);
% note that each cluster could be segmented based on the distance
figure;pcshow(ptCloud.Location,labels+1,'MarkerSize',50)
colormap([hsv(numClusters);[0 0 0]])
title('Segmented Cloud Clusters');hold on
% the centroid of each cluster was caclulated and plotted on the 3D space
% note that any weighting based on the intensity was not applied
s = countcats(categorical(labels));
[~, c_label] = sort(s,'descend');
for i=1:5
roi_idx=find(labels==c_label(i));
roi = ptCloud.Location(roi_idx,:);
xyz=mean(roi,1);
plot3(xyz(1),xyz(2),xyz(3),'o','MarkerSize',20,'MarkerFaceColor','w');hold on
end
result.PNG

댓글 수: 6

ryo tanaka
ryo tanaka 2019년 10월 20일
回答ありがとうございます。
>X,Y,Z,vとありますが、X-Zは、vのボクセル座標を定めるためのグリッドで、それぞれグ>リッド内の値がそのX(やY・Z)座標に対応している、そして、vが実際にそのプロットの>持つ強度である、ということでしょうか。
その通りですね。Vがプロットのもつ強度になります。
>質問者様は、上のようにボクセルを定義して、各重心を求める方針だったと読み取りました>が、下では、ポイントクラウドで処理しています。下のように書けば重心を求めるまでかな>り短く書くことができます。
強度をもつ3次元座標(X,Y,Z)をボクセル上に投影させて、座標点にボリュームをもたせて、そのボクセル上から重心点を求めようとしていました。
ポイントクラウドというのは、ボクセルではなく、各座標点をそのまま点として格納する方法でしょうか?
>重心は一般的な3D上でのものを計算しています。本来はこの過程でvの強度を利用して、>重みつきの重心を定義するのですかね?
そうですね。強度値から重心を求めようとしていました。
>また、4次元、というキーワードがありますが、これを4D画像としてとらえるなら、vの>値を利用せず重心を求めているのでもう少し変更する必要があるとは思います。適宜変更い>ただければと思います。
自分がやってる方法だとvの値は利用されていなかったのですね、勉強不足でした。
参照コードの添付ありがとうございます。図に示してある白い丸が重心点でしょうか?
この5つの丸の座標はどのように出すのでしょうか?
Kenta
Kenta 2019년 10월 20일
편집: Kenta 2019년 10월 20일
以上、丁寧に返信ありがとうございます。またコードをもとにいろいろと試していただけると幸いです。
>>ポイントクラウドというのは、ボクセルではなく、各座標点をそのまま点として格納する>>方法でしょうか?
はい、その通りです。
はい、白い丸が重心です。その座標は以下のように書き換えると取得できます。
こちらではxyzの算術平均で計算しています。
コード実行後にxyzに値が格納されています。
よろしくお願いいたします。
xyz(i,:)=mean(roi,1);
plot3(xyz(i,1),xyz(i,2),xyz(i,3),'o','MarkerSize',20,'MarkerFaceColor','w');hold on
ryo tanaka
ryo tanaka 2019년 10월 20일
返信ありがとうございます。再度添付頂いたコードで実行しましたら、5つの重心点の座標を取得することができました。ありがとうございます。
>こちらではxyzの算術平均で計算しています
添付していただいたコードはxyzのmeanですので、各分布の中心位置を重心位置としているということですね。
できれば強度の値から重心位置を取得したいと思いますので、最後の重心位置を求めるコード部分を変更しないといけませんね。
また、2つ程お聞きしたいのですが、下記のコードでlabels+1と[0 0 0]がどういった意味をしめしているのかわからないので、よろしければ教えて頂きたいです。勉強不足で申し訳ありません。
figure;pcshow(ptCloud.Location,labels+1,'MarkerSize',50)
colormap([hsv(numClusters);[0 0 0]])
Kenta
Kenta 2019년 10월 20일
こんばんは、ご丁寧に返信ありがとうございます。
重心点を得られたとのことでよかったです。
>>できれば強度の値から重心位置を取得したいと思いますので、最後の重心位置を求めるコー>>ド部分を変更しないといけませんね
はい、そうですね。またいろいろと試したうえで不明点があれば新たに質問いただければと思います。
>>labels+1と[0 0 0]がどういった意味をしめしているのかわからない
詳しく目を通していただきありがとうございます。特に両者とも特に意味はありません。
labelsの値が1以上になってほしかったので、こうしましたが、labelsはすべて1以上の値が入っています。特に意図せず書いていました。
[0 0 0]のほうですが、[0 0 0]だけ独立してあるのではなくて、
hsv(numclusters)と縦方向に積み重ねています。; マークが縦に積むということを意味しています。つまり、カラーマップの終点を[0 0 0](黒)にしています。こうするとFigure中の点の色が一部変更されます。特に大きな意味はないです。下で、ほかのいろのパターンも探せます。
ryo tanaka
ryo tanaka 2019년 10월 21일
ご返信ありがとうございます。
>はい、そうですね。またいろいろと試したうえで不明点があれば新たに質問いただければ
>と思います。
ありがとうございます。とても助かります。
質問に回答していただきありがとうございます。
2つとも特に気にすることはないようで安心しました。
また、何か不明な点がありましたら、ご連絡させて頂きます。
Kenta
Kenta 2019년 10월 21일
はい、承知いたしました。解決してよかったです。

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

추가 답변 (0개)

태그

질문:

2019년 10월 17일

댓글:

2019년 10월 21일

Community Treasure Hunt

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

Start Hunting!