matlab便便,豆粑粑 三维散点 根据点密度画等值曲面 contour3D

clear all

clc

close all

%% outline

% plot scatter points, surf, countour3D

%% main

% load data ---------------------------------------------------------------

dat=load('cv_357.txt');

n=length(dat(:,1));

x=dat(:,1);

y=dat(:,2);

z=dat(:,3);

% grid size ---------------------------------------------------------------

n_cv=30;

cv_g=linspace(-pi,pi,n_cv);

cv_int=0.5; % 2 calculate point denisty, half of cubic size

cv_x=cv_g;

cv_y=cv_g;

cv_z=cv_g;

[cv_X,cv_Y,cv_Z]=meshgrid(cv_g,cv_g,cv_g);

% points denisty ----------------------------------------------------------

for i=1:n_cv

ind_x= (x>=cv_x(i)-cv_int) & (x

for j=1:n_cv

ind_y= (y>=cv_y(j)-cv_int) & (y

for k=1:n_cv

ind_z= (z>=cv_z(k)-cv_int) & (z

p_n = sum(ind_x & ind_y & ind_z);

cv_V(j,i,k)= p_n;

end

end

end

% %% 4 conv

% save Mat_post_3D_contour.mat

%

% load Mat_post_3D_contour.mat

% 2 smooth the points density, change 1 to 3, or 5 or 7 according to your

% data

cv_V=smooth3(cv_V,'box',1);

%% original points --------------------------------------------------------

figure(1)

plot3(x,y,z,'.')

axis([-pi pi -pi pi -pi pi])

axis equal tight

view(20,20)

box on

ax=gca;

ax.BoxStyle = 'full';

%% surf version -----------------------------------------------------------

figure(2)

hold on

for i=1:3:n_cv

surf( cv_X(:,:,i),cv_Y(:,:,i),cv_Z(:,:,i),cv_V(:,:,i),...

'edgecolor','none','facealpha',0.05+0.01*i)

end

axis([-pi pi -pi pi -pi pi])

colormap(flipud(jet))

axis equal tight

view(20,20)

box on

ax=gca;

ax.BoxStyle = 'full';

%% contour3D --------------------------------------------------------------

figure(3)

n_layer=9; % layer number

col_mm=flipud(jet(n_layer)); % color map

lev_add=2.3.^[2:n_layer+1]; % which layer you want to plot

lev_beg=80;                 % begin layer

% contour 3D with each layer

for i = 1: n_layer

mm_lev=lev_beg+lev_add(i);

p1 = patch(isosurface(cv_V,mm_lev),'FaceColor',col_mm(i,:),...

'EdgeColor','none','FaceAlpha',0.1+0.015*i);

isonormals(cv_V,p1);

end

% other setting

axis equal

view(20,20);

axis vis3d equal tight

camlight headlight;

lighting phong

set(gca,'GridLineStyle',':')

box on

ax=gca;

ax.BoxStyle = 'full';

grid on

%% logs

% typed by : mm

% mod : 2017年 03月 16日 星期四 15:23:22 CST

% contact me : meatball1982@163.com

%