pcshowを用いて​図を示しています．し​かし，綺麗に図が表示​されないため他に地形​を表示することが可能​でありそうな構文,又​は色合いの工夫などは​ありますでしょうか？

ある実験で示された地形を図としてプロットするために2mm間隔で座標点群を抽出しマットラボ上でプロット図として起こしました．地形の写真上の赤色で囲んだ部分をpcshowを用いて図で綺麗に出したいのですが，どうしてもぼやけてしまい上手く地形が綺麗に出てきません．

カラーマップエディターでカラーのMAXとMINの値を変えてみたり,単色の場合だったらマーカの色をより細かく刻んで設定してみたりしたのですがなかなか綺麗に表示ができません．

今はpcshowを用いて図を示しています．色は単色(出来ればグレー)かjetで図を示したいのですが，他にこのような地形を表示することが可能でありそうな構文,又は色合いの工夫などはありますでしょうか？

clear all

close all

filename ='case1 6h 2mm①-④.txt'
delimiterIn = ',';                
headerlinesIn = 3;               

A = importdata( filename, delimiterIn, headerlinesIn );

figure;
  
pcshow([A.data(:,1), A.data(:,2), A.data(:,3)]);
colormap gray;

ay = -0.00099;
el = 90.00099;
view(ay, el);
grid on;

% colormap('Direction','reverse')
%colorbar('Fontsize',12,'Direction','reverse')     

xlabel('\sl Longshore Distance (m)');
ax.XAxisLocation='bottom';
xticks([0.0525 0.5 1 1.5 2 2.5 3 3.5 4 4.5 4.94750])
xticklabels({'0','0.5','1','1.5','2','2.5','3','3.5','4','4.5','5'})

ylabel('\sl Crossshore Distance (m)');
ax.YAxisLocation = 'right'
yticks([0.00250 0.5 1 1.5 2 2.5 2.9975])
yticklabels({'0','0.5','1','1.5','2','2.5','3'})
%set(gca,'YDir','reverse')

   c = colorbar;
  c.Label.String = '\sl water depth(m)';
  c.Limits = [-0.1 0.1];

Answer 1

[Tohru Kikawada]

pcshow ですと 'MarkerSize' オプションで指定したドットの大きさによって、画像としてみた場合にはぼやけてしまいます。

もともとの点群の投影平面上の位置情報(organizedとも呼ばれます)は残っていますでしょうか。Kinectやステレオカメラ、LiDARで取得したデータであれば投影平面上の位置情報が残っていると思います。その場合、深度方向(Z軸)のデータだけ取り出すことができれば imagesc できれいに表示ができると思います。

一方、投影平面上の位置情報が欠落している場合には点群をメッシュ状にして表示することになると思います。その場合は delaunay 関数などが使えるかもしれません。

元のデータがあればもう少し具体的な提案ができるかもしれません。 ご参考になれば幸いです。

[Taka]

ありがとうございます．

点群の座標についてですが，Canon PowerShotG12というカメラを使って地形の周りを撮影し，その写真をPhotoScanという3Dモデルを作成するソフトを使って解析を行い，完成したモデルから座標を抽出しています． カメラには，位置情報を付けていないので，地形の4隅にマーカーを設置し，基準点を設けてそこからの距離(m)を測定し，それを座標として設定しています．

もとのデータですが，ファイルがかなり重いので，一部ですが載せました．

[Tohru Kikawada]

PhotoScanですとStructure from Motionによる点群生成と思いますので、投影面での位置情報を保持していないunorganizedな点群データですね。そうなると点群を補間してメッシュデータにする必要があります。

griddata を使って点群をグリッドデータとして補間する例を作成してみました。 surf 関数を使えば3次元の表面データとして可視化も可能です。ご参考になれば幸いです。

clear all
close all

filename ='case1 6h 2mm①-④.txt'
if ~exist(filename,'file')
    websave(filename,'https://jp.mathworks.com/matlabcentral/answers/uploaded_files/98621/case1%206h%20sample.txt');
end
delimiterIn = ',';
headerlinesIn = 3;

A = importdata( filename, delimiterIn, headerlinesIn );

figure;
pcshow([A.data(:,1), A.data(:,2), A.data(:,3)]);
colormap gray;

ay = -0.00099;
el = 90.00099;
view(ay, el);
grid on;
ax = gca;

% colormap('Direction','reverse')
%colorbar('Fontsize',12,'Direction','reverse')

xlabel('\sl Longshore Distance (m)');
ax.XAxisLocation='bottom';
xticks([0.0525 0.5 1 1.5 2 2.5 3 3.5 4 4.5 4.94750])
xticklabels({'0','0.5','1','1.5','2','2.5','3','3.5','4','4.5','5'})

ylabel('\sl Crossshore Distance (m)');
ax.YAxisLocation = 'right'
yticks([0.00250 0.5 1 1.5 2 2.5 2.9975])
yticklabels({'0','0.5','1','1.5','2','2.5','3'})
%set(gca,'YDir','reverse')

c = colorbar;
c.Label.String = '\sl water depth(m)';
c.Limits = [-0.1 0.1];

%%Point cloud to grid data
% Create labels
x = A.data(:,1);
y = A.data(:,2);
v = A.data(:,3);
xqlabel = linspace(min(x(:)),max(x(:)),1000);
yqlabel = linspace(min(y(:)),max(y(:)),1000);
[xq,yq] = meshgrid(xqlabel,yqlabel);

% Point cloud to grid
vq = griddata(x,y,v,xq,yq,'natural');

%%Image visualization
figure;
imagesc(xqlabel,yqlabel,vq);
ax = gca;
xlabel('\sl Longshore Distance (m)');
ax.XAxisLocation='bottom';
xticks([0.0525 0.5 1 1.5 2 2.5 3 3.5 4 4.5 4.94750])
xticklabels({'0','0.5','1','1.5','2','2.5','3','3.5','4','4.5','5'})
ylabel('\sl Crossshore Distance (m)');
ax.YAxisLocation = 'right';
yticks([0.00250 0.5 1 1.5 2 2.5 2.9975])
yticklabels({'0','0.5','1','1.5','2','2.5','3'})
colormap gray;
axis image;
c = colorbar;
c.Label.String = '\sl water depth(m)';
c.Limits = [-0.1 0.1];

%%Mesh visualization
figure, surf(xq,yq,vq,...
    'FaceColor',[1 1 1],...
    'EdgeColor','none',...
    'FaceLighting','gouraud');
light('Position',[3 3 1],'Style','local')
axis equal

[Taka]

ありがとうございます．

このような地形になりました．この地形ですが色合いは単色しかできないでしょうか？ 位置によって色合いを変更可能でしょうか？ 自分でやってみたのですが普通のプロット図になってしまい上手くできませんでした・・・・ 色々質問してしまいすいません．

[Tohru Kikawada]

位置によって色合いを変えるというのは具体的にどのように変えたいのでしょうか。たとえば、 jet でLongshoreに方向に色を変えるという感じでしょうか。

[Taka]

縦軸のCrossshore方向にjetのように色を変えたいのですが可能でしょうか？ あと光が写り込んでいて地形が見づらいので消したいのですが，消すことは可能ですか？

[Tohru Kikawada]

surf のオプションに色情報を渡すことで、色付けすることができますよ。また、 light オブジェクトで光源の位置や種類を設定できます。

clear all
close all

filename ='case1 6h 2mm①-④.txt'
if ~exist(filename,'file')
    websave(filename,'https://jp.mathworks.com/matlabcentral/answers/uploaded_files/98621/case1%206h%20sample.txt');
end
delimiterIn = ',';
headerlinesIn = 3;

A = importdata( filename, delimiterIn, headerlinesIn );

figure;
pcshow([A.data(:,1), A.data(:,2), A.data(:,3)]);
colormap gray;

ay = -0.00099;
el = 90.00099;
view(ay, el);
grid on;
ax = gca;

% colormap('Direction','reverse')
%colorbar('Fontsize',12,'Direction','reverse')

xlabel('\sl Longshore Distance (m)');
ax.XAxisLocation='bottom';
xticks([0.0525 0.5 1 1.5 2 2.5 3 3.5 4 4.5 4.94750])
xticklabels({'0','0.5','1','1.5','2','2.5','3','3.5','4','4.5','5'})

ylabel('\sl Crossshore Distance (m)');
ax.YAxisLocation = 'right'
yticks([0.00250 0.5 1 1.5 2 2.5 2.9975])
yticklabels({'0','0.5','1','1.5','2','2.5','3'})
%set(gca,'YDir','reverse')

c = colorbar;
c.Label.String = '\sl water depth(m)';
c.Limits = [-0.1 0.1];

%%Point cloud to grid data
% Create labels
x = A.data(:,1);
y = A.data(:,2);
v = A.data(:,3);
xqlabel = linspace(min(x(:)),max(x(:)),1000);
yqlabel = linspace(min(y(:)),max(y(:)),1000);
[xq,yq] = meshgrid(xqlabel,yqlabel);

% Point cloud to grid
vq = griddata(x,y,v,xq,yq,'natural');

%%Image visualization
figure;
imagesc(xqlabel,yqlabel,vq);
ax = gca;
xlabel('\sl Longshore Distance (m)');
ax.XAxisLocation='bottom';
xticks([0.0525 0.5 1 1.5 2 2.5 3 3.5 4 4.5 4.94750])
xticklabels({'0','0.5','1','1.5','2','2.5','3','3.5','4','4.5','5'})
ylabel('\sl Crossshore Distance (m)');
ax.YAxisLocation = 'right';
yticks([0.00250 0.5 1 1.5 2 2.5 2.9975])
yticklabels({'0','0.5','1','1.5','2','2.5','3'})
colormap gray;
axis image;
c = colorbar;
c.Label.String = '\sl water depth(m)';
c.Limits = [-0.1 0.1];

%%Mesh visualization

% Generate colormap along with the crossshore axis
colors = repmat(permute(jet(numel(xqlabel)),[1 3 2]),[1 numel(yqlabel) 1]);
figure1 = figure;

axes1 = axes('Parent',figure1);
hold(axes1,'on');

surf(xq,yq,vq,colors,'Parent',axes1,...
    'EdgeColor','none',...
    'FaceLighting','gouraud');

% Create a light
light('Parent',axes1,'Position',[1 -6.12323399573677e-17 1]);

% Create ylabel
ylabel('\sl Crossshore Distance (m)');

% Create xlabel
xlabel('\sl Longshore Distance (m)');

grid(axes1,'on');
axis(axes1,'tight');

% Configure the rest of the parameters
set(axes1,'CameraPosition',...
    [-21.8298827263687 0.583298749520218 1.76864432604658],'CameraUpVector',...
    [0.069435328776402 0.00668586171916293 0.997564050259824],'CameraViewAngle',...
    1.45621839490614,'DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',...
    [100.198662635289 2.9867302459682 1],'XTick',...
    [0.0525 0.5 1 1.5 2 2.5 3 3.5 4 4.5 4.9475],'XTickLabel',...
    {'0','0.5','1','1.5','2','2.5','3','3.5','4','4.5','5'},'YAxisLocation',...
    'right','YTick',[0.0025 0.5 1 1.5 2 2.5 2.9975],'YTickLabel',...
    {'0','0.5','1','1.5','2','2.5','3'});

[Taka]

質問に丁寧に対応していただき本当にありがとうございます！