從單個文件，加入管徑後將數據處理成分叉血管，並畫出網路

import numpy as npimport operatorimport copyfrom matplotlib.font_manager import FontPropertiesfrom scipy.interpolate import lagrangeimport randomimport matplotlib.pyplot as plt# 把opt文件內的逗號變為空格np.set_printoptions(threshold=np.inf) #輸出全部矩陣不帶省略號# with open(opt\SubCAM.opt,r) as f:# file = f.read()# file = file.strip().split(
)# for data in file:# # print(data)# data = data.replace(,, )#replace不能用於列表# data = data.replace(-1,0)# with open(txt\SubCAM1.txt,a) as f:# f.write(data+
)##########################################data = np.loadtxt(txt\SubCAM.txt)x1 = data[:,5]#x起點坐標x2 = data[:,9]#x終點坐標y1 = data[:,6]#y起y2 = data[:,10]#y起z1 = data[:,4]#IDpartz2 = data[:,8]#IDpartdiam = data[:,12]s1 = [a1 for a1 in range(1,len(x1)-1) if z1[a1]==z2[a1-1]!=0 or z1[a1]!= z2[a1-1]]#id相同不等於0，或id不同lx = []#x1,x2相同的部分組成的列表kx = []#x2不同的單個數據組成的列表for i1 in range(len(s1)-1): b1 = x1[s1[i1]:s1[i1+1]] b1 = b1.tolist() lx.append(b1) b2 = x2[s1[i1]-1] kx.append(b2)kx = kx[1::]for i2 in range(len(kx)): lx[i2].append(kx[i2])# ##########################################ly = []ky = []km = []for i3 in range(len(s1)-1): b3 = y1[s1[i3]:s1[i3+1]] b3 = b3.tolist() ly.append(b3) b4 = y2[s1[i3]-1]#y最後一個不相等的數 dm = diam[s1[i3]-1] ky.append(b4) km.append(dm)ky = ky[1::]km = km[1::]for i4 in range(len(ky)): ly[i4].append(ky[i4]) ly[i4].append(km[i4])# ##########################################lz = []kz = []for i5 in range(len(s1)-1): b5 = z1[s1[i5]:s1[i5+1]] b5 = b5.tolist() lz.append(b5) b6 = z2[s1[i5]-1]#y最後一個不相等的數 kz.append(b6)kz = kz[1::]for i6 in range(len(kz)): lz[i6].append(kz[i6])# print(lz)#####################################################lxwei =copy.deepcopy(lx)#深拷貝，不改變原來的值lywei = copy.deepcopy(ly)for i7 in range(len(lz)): lx[i7].insert(0,lz[i7][0])#插入起點id# ly[i7].insert(0,lz[i7][0])#ly插不插都行for i8 in range(len(lz)): lxwei[i8].insert(0,lz[i8][-1])#插入終點id# lywei[i8].insert(0,lz[i8][-1])#ly插不插都行# print(lx)# print(ly)for q1 in range(len(lx)): for q2 in range(len(ly[q1])): lx[q1].append(ly[q1][q2])for p1 in range(len(lxwei)): for p2 in range(len(lywei[p1])): lxwei[p1].append(lywei[p1][p2])lx.sort(key=operator.itemgetter(0))#排序，只按照第一個索引大小排序tou = lxlxwei.sort(key=operator.itemgetter(0)) wei = lxwei # print(tou)# print(wei)# #########################################final = [] for i in range(len(tou)-1): for dwei in wei: if dwei[0]==tou[i][0] and dwei[0]==tou[i+1][0]: zuhe = [dwei,tou[i],tou[i+1]] final.append(zuhe)# print(final)# ###################################################finaldata= []for i in range(len(final)): a = final[i][0][1::] b = final[i][1][1::] c = final[i][2][1::] finaldata.append(a) finaldata.append(b) finaldata.append(c)# print(finaldata)############################################管徑沒插進去for i in range(len(finaldata)): if len(finaldata[i])%2==0: finaldata[i].append(60.76)# ####################################### #可視化，未插值# w1 = finaldata[3] # w2 = finaldata[4]# w3 = finaldata[5]# x1 = w1[0:(len(w1)-1)//2]# y1 = w1[((len(w1)-1)//2):len(w1)-1]# x2 = w2[0:(len(w2)-1)//2]# y2 = w2[((len(w2)-1)//2):(len(w2)-1)]# x3 = w3[0:((len(w3)-1)//2)]# y3 = w3[((len(w3)-1)//2):(len(w3)-1)]# plt.plot(x1,y1)# plt.plot(x2,y2)# plt.plot(x3,y3)# plt.show()# ###################################################### #插值，起點不同for i in range(len(finaldata)): lenth = len(finaldata[i])//2 a = finaldata[i][(lenth-2):lenth] if a[0]==a[1]: a[1] = a[1]+0.000001 b = finaldata[i][-3:-1] f = lagrange(a,b) if i % 3 != 0: fx = [random.uniform(a[0],a[1]) for _ in range(11-lenth)] else: fx = [random.uniform(a[0],a[1]) for _ in range(10-lenth)] fy = f(fx) for u in range(len(fx)-1,-1,-1): finaldata[i].insert((lenth-1),fx[u]) for n in range(len(fy)): finaldata[i].insert(-2,fy[n]) if i % 3 != 0: lenth = len(finaldata[i])//2 finaldata[i].remove(finaldata[i][lenth]) finaldata[i].remove(finaldata[i][0])# ###################################################### #保存分叉數據為txt形式for i in range(len(finaldata)): hang = finaldata[i] with open(C:\Users\Administrator\Desktop\分叉數據\a1.txt,a) as f: f.write(str(hang)+
)# ##################################################### 一個一個畫出分叉# for i in range(len(finaldata)//3):# w1 = finaldata[0+3*i] # w2 = finaldata[1+3*i]# w3 = finaldata[2+3*i]# x1 = w1[0:10]# y1 = w1[10:20]# x2 = [w1[9]]+w2[0:10]# y2 = [w1[19]]+w2[10:20]# x3 = [w1[9]]+w3[0:10]# y3 = [w1[19]]+w3[10:20]# plt.plot(x1,y1)# plt.plot(x2,y2)# plt.plot(x3,y3)# plt.show()# ################################################################# # # 畫出整個血管網路# lx = []# ly = []# for i in range(len(finaldata)):# lx.append(finaldata[i][0:10])# ly.append(finaldata[i][10:20])# for i in range(len(lx)):# plt.plot(lx[i],ly[i])# plt.xlabel(x axis,size=16)# plt.ylabel(y axis,size=16)# font = FontProperties(fname=rc:windowsfontssimsun.ttc,size=16)# plt.title(u血管分叉,fontproperties=font)# plt.show()

推薦閱讀：