今天是因为可以用py而高兴的一天。
昨天老板淡淡地回了一句，sierpinski地毯画得挺好的。
我思考了五秒钟之后，想起来作业其实是sierpinski垫片。
三角垫片比地毯难做多了。
因为你用的像素点，其实是矩形，你可以把像素点当矩形用。
三角形，随机算法可以，你要我用传统算法来解那个三角垫片，我只想给你以窝jio。除非你给我一个三角网格。
所以在jpg里玩，没tai意nan思le。你说你解完，反不反走样。。。是吧。。。反走样很蠢，走样了很丑。
所以，我决定用.ply来存这玩意。。。话说，既然这么用，干脆搞个大新闻好了。
不然你说，在三维空间里存了一堆平面上的三角，哈哈哈哈哈次藕爆了。
可是道理我都懂，为啥底是矩形呢。。。[图片来自百度，侵删]
所以，我们要整个更好看的哈哈哈哈哈。
算法逻辑：
每个四面体生成四个小四面体并多出六个顶点。
最终是要输出三角面片，可是在这之前我一直是体在拆解啊。
所以我中间过程存储体的信息，最后一步拆成四个面就好了。
每一轮，四面体信息矩阵扩大四倍。
每个四面体拆解都会成四个小四面体存入对应四行。
顶点通过取中点取到六个，信息补上6个。
import numpy as np
import cv2 def beiup(x1,x2,x3,x4):
x=np.zeros((6,3),dtype=float)
for j in range(3):
x[0][j]=(x1[j]+x2[j])/2
x[1][j]=(x2[j]+x3[j])/2
x[2][j]=(x3[j]+x1[j])/2
x[3][j]=(x4[j]+x1[j])/2
x[4][j]=(x4[j]+x2[j])/2
x[5][j]=(x4[j]+x3[j])/2
return x
def yeiup(yei,jsq):
new=np.zeros((4,4),dtype=int)
new[0][0]=yei[0]
new[0][1]=jsq
new[0][2]=jsq+2
new[0][3]=jsq+3
new[1][0]=jsq
new[1][1]=yei[1]
new[1][2]=jsq+1
new[1][3]=jsq+4
new[2][0]=jsq+2
new[2][1]=jsq+1
new[2][2]=yei[2]
new[2][3]=jsq+5
new[3][0]=jsq+3
new[3][1]=jsq+4
new[3][2]=jsq+5
new[3][3]=yei[3]
jsq+=6
return new,jsq
def feiwrite(yeimat,times):
k=4**times
feimat=np.zeros((k*4,3),dtype=int)
for i in range(k):
feimat[i*4+0][0]=yeimat[i][0]
feimat[i*4+0][1]=yeimat[i][2]
feimat[i*4+0][2]=yeimat[i][1]
feimat[i*4+1][0]=yeimat[i][0]
feimat[i*4+1][1]=yeimat[i][1]
feimat[i*4+1][2]=yeimat[i][3]
feimat[i*4+2][0]=yeimat[i][0]
feimat[i*4+2][1]=yeimat[i][3]
feimat[i*4+2][2]=yeimat[i][2]
feimat[i*4+3][0]=yeimat[i][1]
feimat[i*4+3][1]=yeimat[i][2]
feimat[i*4+3][2]=yeimat[i][3]
return feimat def plywrite(fname,beimat,feimat,potnumber,trinumber):
potpiplist=np.zeros((3),dtype=float)
tripiplist=np.zeros((3),dtype=int)
fply = open(fname,'w')
fply.write('ply\n')
fply.write('format ascii 1.0\n')
fply.write('comment VCGLIB generated\n')
fply.write('element vertex '+str(potnumber)+'\n')
fply.write('property float x\n')
fply.write('property float y\n')
fply.write('property float z\n')
fply.write('element face '+str(trinumber)+'\n')
fply.write('property list uchar int vertex_indices\n')
fply.write('end_header\n')
for i in range(potnumber):
potpiplist[0]=beimat[i][0]
potpiplist[1]=beimat[i][1]
potpiplist[2]=beimat[i][2]
fply.write('%0.4f %0.4f %0.4f\n' %(potpiplist[0],potpiplist[1],potpiplist[2]))
for i in range(trinumber):
tripiplist[0]=feimat[i][0]
tripiplist[1]=feimat[i][1]
tripiplist[2]=feimat[i][2]
fply.write('%d %d %d %d\n' %(3,tripiplist[0],tripiplist[1],tripiplist[2]))
fply.close()
return fname
解释一下变量名。
我最近在捣腾一个TPS的大项目。
ply文件有两个主要部分
一个是点坐标信息，beimat 顶点数*3, float
一个是面片顶点信息，feimat 面片数*3, int 对应beimat的第i个。
以前没有四面体一说，yeimat 体数*4 , int
bei是我的名字，fei是我姐上的名。
口哼哼，我不管，反正项目我最后封装掉的。
beiup解决了如何拆出六个中点；
yeiup解决了如何拆出四个体；
feiwrite解决了如何从一个体拆除四个面片。
import numpy as np
import cv2
import tps
import sie times=1
if times<4:
beimat=np.zeros((600,3),dtype=float)
if times<6:
beimat=np.zeros((2100,3),dtype=float)
if times<10:
beimat=np.zeros((9000,3),dtype=float) yeimat=np.zeros((1,4),dtype=int)
'''
0,0,0; 1,0,0; 0.5,g3/2,0 ;0.5 g3/6,g6/3)
'''
beimat[1][0]=1
beimat[2][0]=0.5
beimat[2][1]=(3**0.5)/2
beimat[3][0]=0.5
beimat[3][1]=(3**0.5)/6
beimat[3][2]=(6**0.5)/3 yeimat[0][1]=1
yeimat[0][2]=2
yeimat[0][3]=3
###end for setting jsq=4
for k in range(times):
yeinew=np.zeros((4**k*4,4),dtype=int)
for n in range(4**k):
beimat[jsq:jsq+6,:]=sie.beiup(beimat[int(yeimat[n][0])],beimat[int(yeimat[n][1])],beimat[int(yeimat[n][2])],beimat[int(yeimat[n][3])])
[yeinew[4*n:4*n+4,:],jsq]=sie.yeiup(yeimat[n],jsq)
yeimat=yeinew feimat=sie.feiwrite(yeimat,times) fname=str("sie"+str(int(times))+".ply")
print(jsq) fname=sie.plywrite(fname,beimat[0:jsq,:],feimat,jsq,4**times*4)
主程序有三部分：
决定beimat的大小
初始化第一个四面体（顶点坐标+编号）
主循环体+输出。
你说我不能算出来beimat应该有多大么？
毕竟是学树学的，这如果算不出来实在是无能。。。问题是，我有必要算么？
如果你要算的话，就是x=x*4-6的一个递归数列。。。可以，但没必要。
我jsq最后会告诉我有多少个坐标点的。。。我只要，一开始不把beimat设太小就好了。
这里四个是1-4次迭代；下方几个都是5/6次迭代。
嘿嘿嘿！py万岁！