前言:在地球物理勘探,流体空间分布等多种场景中,定位空间点P(x,y,x)的物理属性值Q,并绘制三维空间分布图,对我们洞察空间场景有十分重要的意义。
1. 三维立体图的基本要件:
- 全空间网格化
- 网格节点的物理属性值
2.数据准备
数据不易贴,我放在了百度网盘:点击下载数据
大概如下形式:
TIP:
这里的数据矩阵为v(5276),可以看成一本27页纸,每页绘制了5*6的网格,然后27页纸叠在一起。当你理解本图绘制后,数据可以随意制作。
3.主要函数:slice、isosurface、patch
NOTE:想了解每个函数的用法可以在matlab command windows中 输入 doc slice即可获得slice的matlab详解
4.主要代码:
TIP:你可以将代码复制到一起运行,此处我是单独展示,分别绘图
1.单独切片,切法可自定义。便于展示局部细节
clc,clear
v=csvread(\'vdata.txt\');
v=reshape(v,5,27,6);
[x,y,z]=meshgrid(1:27,1:5,1:6);
h=figure(1);
set(h,\'name\',\'取单切片\')
subplot(221)
slice(x,y,z,v,[],[1],[]);
shading interp
set(gca,\'zdir\',\'reverse\');
axis equal
grid on
subplot(222)
slice(x,y,z,v,[],[2],[]);
shading interp
colormap(\'jet\')
set(gca,\'zdir\',\'reverse\');
axis equal
grid on
subplot(223)
slice(x,y,z,v,[],[3],[]);
shading interp
set(gca,\'zdir\',\'reverse\');
axis equal
grid on
subplot(224)
slice(x,y,z,v,[],[4],[]);
shading interp
set(gca,\'zdir\',\'reverse\');
axis equal
grid on
效果图:
2.全空间立体切片
h2=figure(2);
set(h2,\'name\',\'全空间切片\',\'MenuBar\',\'none\',\'ToolBar\',\'none\')
slice(x,y,z,v,[1:2:27],[2 3 4],[2 3 4 5])
shading interp
colorbar
colormap(\'jet\')
set(gca,\'zdir\',\'reverse\');
axis equal
grid on
box on
效果图:
3.立体包络图
h3=figure(3);
set(h3,\'name\',\'定值包络立体图\',\'MenuBar\',\'none\',\'ToolBar\',\'none\')
set(gcf,\'InvertHardcopy\',\'off\')
fw=350; %%此值为最外层包络面取值
fv=isosurface(x,y,z,v,fw);
p=patch(fv);
set(p,\'facecolor\',\'b\',\'edgecolor\',\'none\');
patch(isocaps(x,y,z,v, fw), \'FaceColor\', \'interp\', \'EdgeColor\', \'none\');
colorbar
colormap(\'jet\')
box on
daspect([1,1,1])
view(3)
set(gca,\'zdir\',\'reverse\',\'color\',[0.2,0.2,0.2]);
camlight
camproj perspective
lighting phong
axis equal
grid on
title([\'最外层表面的值为: \' , num2str(fw)]);
NOTE:上述代码中的fw是作主要可变参数,控制最外层包络面的值
如下演示多个fw下不同的效果图: