データ点を補完して3次元マップを作成したい

30 views (last 30 days)
kou sasai
kou sasai on 15 Jun 2021
Commented: kou sasai on 24 Jun 2021
Matlabを用いて下の図のような半円柱内の磁場分布を示す3次元マップを作成したいです。その過程の補完・マップ作成について教えていただきたいです。
高さをz方向として、原点を通るxz平面に関して格子点1cm間隔の磁場データを取得してあります。データはエクセルに記録し、位置x,y(=0),zとその地点の磁場強度Bを[x,0,z,B]として4× n の形で記録しています。(原点は円柱の中心にしています)
このデータを出発点に以下の操作をしたいです。
step1.xz平面の磁場データの補完(内挿)
   1cm間隔のデータを線形補完して滑らかな分布を示すデータにしたいです
   interp2を使用することを考えましたがどのように引数を設定すれば良いかわからず詰まっています。
step2.3次元への拡張
   原点からの距離と高さが等しい点は同じ磁場強度として扱うことで平面の磁場分布データを立体に拡張したいです
   (得られた平面のデータを、z軸(高さ方向の軸)を回転の軸としてまわすことで3次元に拡張するイメージです)
step3.マップ表示
   データから図のように色で磁場強度を示す半円柱の3次元マップを作成したいです
   またx軸上およびz軸上の磁場強度を表すグラフも一緒に作成したいです(まさに図のような形で)
どこか一部分だけでも構いませんのでご教授いただけると幸いです。

Accepted Answer

Hernia Baby
Hernia Baby on 15 Jun 2021
Step1
 interp1が適切かと思われます。
 図を参考にすると、x = 0:10:80; z = 0:10:150;とサイズが違うからです。
 interp1でxを補間し、meshgridで n×mこの行列を作るのがよいと思われます。
Step2
 z軸中心に回転するのであればcylinder関数が使えます。
Step3
 パラメータがないので、磁場強度を適当な関数で行っています。
 まずは円筒を作ります。
clear,clc,close all;
% 最大半径40mm,中心(40,40),高さ150mmの円筒
xlen = 80;
z = 0:5:150;
% 半径1~40の円をつくる
r = 1:40;
for i = 1:length(r)
[xx,yy,zz] = cylinder(r(i));
X(i,:) = xx(1,:);
Y(i,:) = yy(1,:);
end
clear i r xx yy zz
% 中心(0,0)から(40,40)にシフト
X = X + 40;
Y = Y + 40;
 次に磁場強度を作ります。
 ポイントはmeshgridでf(x,z)の関数にすることです。
cnt = 1;
while cnt <= length(X(:,1))
x = X(cnt,:);
% point f(x,z)の関数にする
[xt,zt] = meshgrid(x,z);
rxt = 110.*(1-sin(2*pi*xt./xlen/2)); % xについてはここをinterp1
rzt = 250.*(1-sin(2*pi*zt./max(z)/2)); % zについてはここをinterp1
mesh=rzt.*rxt;
str(:,:,cnt) = mesh';
cnt = cnt+1;
end
 最後に次元を合わせます。
 ここでは、(θ,r) -> (x,y)にしているので、21×半径の数×高さの次元に直してます。
 コードとしては正直汚いです。
Z = repmat(z,length(X(1,:)),1)';
Z = repmat(Z,1,1,length(X(:,1)));
X = X'; Y = Y';
X = permute(repmat(X,1,1,length(z)),[1 2 3]);
Y = permute(repmat(Y,1,1,length(z)),[1 2 3]);
Z = permute(Z,[2 3 1]);
str = permute(str,[1 3 2]);
clear x z xt zt rxt rzt mesh cnt
 今回はsurfで円盤を積み重ねます。
 x軸上&z軸上の磁場強度を見る場合はplot関数を使うと可能です。
 subplotを使うと並記できるのでお勧めです。
cnt = 1;
while cnt <= length(X(1,1,:))
surf(X(:,:,cnt),Y(:,:,cnt),Z(:,:,cnt),str(:,:,cnt)), hold on
xlabel('X[mm]'); ylabel('Y[mm]'); zlabel('Z[mm]');
cnt = cnt +1;
end
ylim([40 80]);
  3 Comments
Hernia Baby
Hernia Baby on 22 Jun 2021
xは等間隔と仮定した場合、reshape関数で磁場強度をmeshにしてinterp2をかけるとよさそうです。
reshapeまでの方法を書いておきます。
clc,clear;
x = repmat([1:20],1,5);
z = repelem([1:5],20);
M = randi(10,[1,100]);
table(x',z',M','VariableNames',{'x','z','B'})
ans = 100×3 table
x z B __ _ _ 1 1 4 2 1 4 3 1 8 4 1 5 5 1 6 6 1 2 7 1 6 8 1 9 9 1 1 10 1 5 11 1 2 12 1 4 13 1 5 14 1 9 15 1 5 16 1 9
meshの形にします
xu = unique(x);
zu = unique(z);
Mu = reshape(M,[],length(zu))
Mu = 20×5
4 3 9 6 8 4 10 2 5 2 8 5 6 9 7 5 7 4 9 9 6 4 5 2 5 2 1 4 4 1 6 4 6 1 1 9 6 5 2 4 1 1 5 5 4 5 2 6 8 5
surf(zu,xu,Mu,'FaceAlpha',0.5)
kou sasai
kou sasai on 24 Jun 2021
疑問全て解決しました。ご丁寧に対応してくださりありがとうございます。

Sign in to comment.

More Answers (0)

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!