import numpy as npimport operatorimport osimport copyfrom matplotlib.font_manager import FontPropertiesfrom scipy.interpolate import lagrangeimport randomimport matplotlib.pyplot as pltimport math np.set_printoptions(suppress=True)# 把opt文件内的逗号变为空格#数据在我的百度云数据库txt文件,及opt文件np.set_printoptions(threshold=np.inf) #输出全部矩阵不带省略号random.seed(10)##########################################data = np.loadtxt('txt//final37.txt')# data = data[0:100]#抽取一部分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]!=-1 or z1[a1]!= z2[a1-1]]#id相同不等于0，或id不同# print(s1)lx = []#x1,x2相同的部分组成的列表lxqi = []lxzg = []for i1 in range(len(s1)-1):b1 = x1[s1[i1]:s1[i1+1]]b1 = b1.tolist()b2 = x2[s1[i1+1]-1]#s1[i1]相当于a1#b1 = b1 + [b2]#把与x2最后相连的一个数和x1拼接起来b5 = z1[s1[i1]]#x,y起点idb1qi_id = [b5]+b1 +[b2]b6 = z2[s1[i1+1]-1]#x,y终点idb1zg_id = [b6] + b1+[b2]lx.append(b1)lxqi.append(b1qi_id)lxzg.append(b1zg_id)###################################################ly = []#y坐标以及管径大小for i3 in range(len(s1)-1):b3 = y1[s1[i3]:s1[i3+1]]b3 = b3.tolist()b4 = y2[s1[i3+1]-1]#y最后一个不相等的数b3 = b3 + [b4]dm = diam[s1[i3+1]-1]b3 = b3 + [dm]#加上管径ly.append(b3)######################################################带有起点id的x坐标与y坐标合并for q1 in range(len(lxqi)):for q2 in range(len(ly[q1])):lxqi[q1].append(ly[q1][q2])#带有终点id的x坐标与y坐标合并for p1 in range(len(lxzg)):for p2 in range(len(ly[p1])):lxzg[p1].append(ly[p1][p2])lxqi.sort(key=operator.itemgetter(0))#排序，只按照第一个索引大小排序tou = lxqilxzg.sort(key=operator.itemgetter(0)) wei = lxzg # #########################################toudeng = []weideng = []for dwei in wei:for i in range(len(tou)-1):if dwei[0] ==tou[i][0] and dwei[0]==tou[i+1][0]:toud = [dwei,tou[i],tou[i+1]]toudeng.append(toud)for dtou in tou:for i in range(len(wei)-1):if dtou[0] == wei[i][0] and dtou[0]==wei[i+1][0]:weid = [wei[i],wei[i+1],dtou]weideng.append(weid)# ###################################################datatoudeng = []dataweideng = []#去掉起点idfor i in range(len(toudeng)):a = toudeng[i][0][1::]b = toudeng[i][1][1::]c = toudeng[i][2][1::]d = [a]+[b]+[c]datatoudeng.append(d)for i in range(len(weideng)):a1 = weideng[i][0][1::]b1 = weideng[i][1][1::]c1 = weideng[i][2][1::]d1 = [a1]+[b1]+[c1]dataweideng.append(d1)# print(dataweideng)#####################################################################判断管径信息是否加进列表，若未加进则只为x,y坐标，为偶数for i in range(len(dataweideng)):a = dataweideng[i]assert len(a[0])%2==1assert len(a[1])%2==1assert len(a[2])%2==1for i in range(len(datatoudeng)):a = datatoudeng[i]assert len(a[0])%2==1assert len(a[1])%2==1assert len(a[2])%2==1finaldata = datatoudeng +dataweideng#未插值final = datatoudeng #所有分叉，头等分叉，尾等分叉# print(final)#################################################################################插值方法：第二张策略，在端点处，补上与端点的值相等的数，然后使每个分支维度相等。finaldata = []for i in range(len(final)):#final[i]代表一个分叉，它有三个不同长度的分支zhu = final[i][0]zuo = final[i][1]you = final[i][2]zhu_diam = [zhu[-1]]zuo_diam = [zuo[-1]]you_diam = [you[-1]]zhu_x = zhu[0:len(zhu)//2]zuo_x = zuo[0:len(zuo)//2]you_x = you[0:len(you)//2]zhu_y = zhu[len(zhu)//2:(len(zhu)-1)]zuo_y = zuo[len(zuo)//2:(len(zuo)-1)]you_y = you[len(you)//2:(len(you)-1)]#从这里开始插值，对于头部相等的数主分支，在头部插值，左右分支分别在尾部端点插值，对于final37希望每一个维度为60，所以，让他们每一个列表插入（30-len(zhu_x))个值zhu_insert_x = zhu_x[0]#主分支的第一个坐标zhu_insert_y = zhu_y[0]zuo_insert_x = zuo_x[-1]#左分支的最后一个坐标zuo_insert_y = zuo_y[-1]you_insert_x = you_x[-1]#右分支的最后一个坐标you_insert_y = you_y[-1]for i in range(30-len(zhu_x)):zhu_x.insert(0,zhu_insert_x)zhu_y.insert(0,zhu_insert_y)#左右分支在尾部插入与尾部端点坐标相等的值，因为左右分支头部要去掉一个值，所以应比主分支多插入一个for i in range(31-len(zuo_x)):zuo_x.append(zuo_insert_x)zuo_y.append(zuo_insert_y)for i in range(31-len(you_x)):you_x.append(you_insert_x)you_y.append(you_insert_y)#从这里去掉左右分支开头与主分支相等的部分zuo_x = zuo_x[1::] you_x = you_x[1::]zuo_y = zuo_y[1::]you_y = you_y[1::]#这里将x和y列表再接起来zhu_xy = zhu_x + zhu_yzuo_xy = zuo_x + zuo_yyou_xy = you_x + you_y#这里再将坐标点与管径接起来zhu = zhu_xy + zhu_diamzuo = zuo_xy + zuo_diamyou = you_xy + you_diamfencha = [zhu]+[zuo]+[you]finaldata.append(fencha)final = np.array(finaldata) #数组维度(-1,3,61)print(final.shape)################################################################################## final = final[0:2,:,:]#选取一个分叉测试旋转效果x = final[:,:,0:30]y = final[:,:,30:60]diam = final[:,:,-1]diam = diam.reshape(-1,3,1)##########################################旋转def rotate(angle,valuex,valuey):rotatex = math.cos(angle)*valuex -math.sin(angle)*valueyrotatey = math.cos(angle)*valuey + math.sin(angle)* valuexreturn rotatex,rotateyrotatedata = []for i in range(0,360,3): #每隔3度旋转一次x1,y1 = rotate(i,x,y)rotate_final = np.concatenate((x1,y1,diam),axis=2)rotatedata.append(rotate_final)finaldata = []for file in rotatedata:for data in file:finaldata.append(data)finaldata = np.array(finaldata)max1 = np.max(finaldata)min1 = np.min(finaldata)print(max1)print(min1)print(finaldata.shape)###############################################################归一化前可视化单张图片# final = finaldata# for i in range(len(final)):# #plt.figure(figsize=(128,128),dpi=1)#plt.plot(final[i][0][0:30],final[i][0][30:60])#plt.plot(final[i][1][0:30],final[i][1][30:60])#plt.plot(final[i][2][0:30],final[i][2][30:60])# #plt.axis('off')#plt.show()# #plt.savefig('C:\\Users\\Administrator\\Desktop\\调整分辨率\\原始图\\resouce%d.jpg' %(i),dpi=1)# #plt.close()#################################################归一化处理不去除管径信息finalSubCAM = []final = finaldatafor i in range(len(final)):finalx = final[i][:,0:30]#(7,10,30)finaly = final[i][:,30:60]#(10,30,60)diameter = final[i][:,-1]diameter = diameter.reshape(3,1)Xmax = np.max(finalx)Xmin = np.min(finalx)Ymax = np.max(finaly)Ymin = np.min(finaly)Dmax = np.max(diameter)Dmin = np.min(diameter)normx = (finalx-Xmin)/(Xmax-Xmin)normy = (finaly-Ymin)/(Ymax-Ymin)normd = (diameter-Dmin)/(Dmax-Dmin)normxy = np.concatenate((normx,normy,diameter),axis=1) #加入原始管径diameter，或归一化管径normdfinalSubCAM.append(normxy)finaldata = np.array(finalSubCAM)np.random.shuffle(finaldata)print(np.min(finaldata))print(np.max(finaldata))print(finaldata.shape)#######################################################一种普通的可视化方法，此时画出来的图端点都连在了原点位置# finaldata = finaldata.tolist()# for i in range(len(finaldata)):#plt.plot(finaldata[i][0][0:30],finaldata[i][0][30:60])#plt.plot([finaldata[i][0][29]]+finaldata[i][1][0:30],[finaldata[i][0][59]]+finaldata[i][1][30:60])#plt.plot([finaldata[i][0][29]]+finaldata[i][2][0:30],[finaldata[i][0][59]]+finaldata[i][2][30:60])#plt.xticks(np.arange(0,1,0.1))#plt.yticks(np.arange(0,1,0.1))#plt.show()#############################################################np.save('C:\\Users\\Administrator\\Desktop\\重新整理血管网络\\final37在端点处插值.npy',finaldata)####################################################################每100张图片显示在一张图中# rows,cols = 10, 10# fig,axs = plt.subplots(rows,cols)# cnt = 0# for i in range(rows):#for j in range(cols):# xy = finaldata[cnt]#第n个分叉图，有三个分支，每个分支21个数# for k in range(len(xy)):# x = xy[k][0:30]# y = xy[k][30:60]# axs[i,j].plot(x,y,linewidth=2)# axs[i,j].axis('off')# cnt +=1# plt.show()
如果觉得《final37另一种插值方法：在端点处插入相等的值》对你有帮助，请点赞、收藏，并留下你的观点哦！
