- /
-
3DFractal/Surf_Cube_out
on 22 Nov 2023
- 8
- 2
- 0
- 0
- 1977
drawframe(1);
Write your drawframe function below
function drawframe(f1)
[X,Y] = meshgrid(-5:1:5);
R = sqrt(X.^2 + Y.^2) + eps;
Z = 5*sin(R)./R;
[f,v,~]=surf2patch(X,Y,Z);
TR.Points=v;
TR.ConnectivityList=f;
m=1;
plus_minus=1;
% [v,f] = createCube;
% TR.Points=v;
% TR.ConnectivityList=f;
index=round(linspace(0,2,48));
i=index(f1);
TR_out=Fractal_Kochsnow3D_square(TR,m,i,plus_minus);
c=jet(size(TR_out.Points,1));
% trisurf(TR_out,'FaceColor','blue')
patch('Faces',TR_out.ConnectivityList,'Vertices',TR_out.Points,'FaceVertexCData',c,'EdgeColor','k','FaceColor','interp','LineWidth',0.5);
axis equal off;
view(30+2*f1,30);
function TR_out = Fractal_Kochsnow3D_square(TR,m,n,plus_minus,varargin)
%% 缺省值处理。。 默认为向外扩展,可选为向内或者向外
if isempty(varargin)
inorout = 'out';
else
inorout = varargin{1};
end
%% 循环获取结果
%inorout1={'out';'in'};
if n==0
TR_out = TR;
else
for i = 1:n
P=[];
T=[];
for j = 1:size(TR.ConnectivityList,1)
%TR1 = triangulation([1,2,3],TR.Points(TR.ConnectivityList(j,:),:));
TR1.Points=TR.Points(TR.ConnectivityList(j,:),:);
TR1.ConnectivityList=TR.ConnectivityList(j,:);
Temp = KochSeed_square(TR1,m,plus_minus,inorout);
if j==1
P=[P;Temp.Points];
T=[T;Temp.ConnectivityList];
else
T=[T;Temp.ConnectivityList+size(P,1)];
P=[P;Temp.Points];
end
end
TR_out.Points = P;
TR_out.ConnectivityList = T;
TR=TR_out;
end
end
end
function TR_out=KochSeed_square(TR,m,plus_minus,varargin)
%% 缺省值处理。。 默认为向外扩展,可选为向内或者向外
if isempty(varargin)
inorout = 'out';
else
inorout = varargin{1};
end
% if strcmpi(inorout,'in')
% inorout = 1;
% else
% inorout = -1;
% end
%% 处理不同的inorout
if strcmpi(inorout,'in')
TR.Points=TR.Points([1,4,3,2],:);
end
%正方形的边长dl
dl=norm(TR.Points(1,:)-TR.Points(2,:));
%以dl/2构建正四面体,那么它的高h
h=dl/3*m;
%% points5-12
nn=[2:4,1];
points=[];
for i = 1:4
temp = NequalSegments_3Dline(TR.Points(i,:),TR.Points(nn(i),:),3);
points=[points;temp];
end
%% points12-16
nn1=[5,6]-4;
nn2=[10,9]-4;
for i = 1:2
temp = NequalSegments_3Dline(points(nn1(i),:),points(nn2(i),:),3);
points=[points;temp];
end
%% points17-20
TR_squ = triangulation([1,2,3],TR.Points(1:3,:));
faceNormal = TR_squ.faceNormal;
points=[points;points(end-3:end,:)+plus_minus*faceNormal*h];
%% points合并
points=[TR.Points;points];
%% ConnectivityList
ConnectivityList=[[1,5,13,12];
[5,6,15,13];
[6,2,7,15];
[7,8,16,15];
[8,3,9,16];
[9,10,14,16];
[10,4,11,14];
[11,12,13,14];
[13,15,19,17];
[15,16,20,19];
[16,14,18,20];
[14,13,17,18];
[17,19,20,18]
];
TR_out.Points=points;
TR_out.ConnectivityList=ConnectivityList;
end
function points = NequalSegments_3Dline(posA,posB,n)
points=[];
%% 方法1 采取从左到右进行等分,每次等分都是1:的形式。不断以新的点代替posA,这样会又累计误差
% for i = n-1:-1:1
% k=1/i;
% temp=(posA+k*posB)/(1+k);
% posA = temp;
% points=[points;temp];
% end
%% 方法2 按照公式,始终参照AB两个点,个人觉得该方法好.
up=1:n-1;
down=n-1:-1:1;
for i = 1:n-1
k=up(i)/down(i);
temp=(posA+k*posB)/(1+k);
points=[points;temp];
end
end
end