在端点处插入相等的值，主分支最好两个相等，然后，又去除了最后一个数，左右分支，都和主分支断开了
import numpy as np
import operator
import os
import copy
from matplotlib.font_manager import FontProperties
from scipy.interpolate import lagrange
import random
import matplotlib.pyplot as plt
import 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]#IDpart
z2 = data[:,8]#IDpart
diam = 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 = lxqi
lxzg.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 = []
#去掉起点id
for 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==1
for i in range(len(datatoudeng)):a = datatoudeng[i]assert len(a[0])%2==1assert len(a[1])%2==1assert len(a[2])%2==1
finaldata = 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,rotatey
rotatedata = []
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 = finaldata
for 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)
#####################################################
#final37主分支后两个数相等，那去掉主分支最后一个坐标
finaldata = finaldata.tolist()
final = []
for i in range(len(finaldata)):zhu = finaldata[i][0]zuo = finaldata[i][1]you = finaldata[i][2]#单独分开x,y,列表zhu_x = zhu[0:30]zhu_y = zhu[30:60]zhu_diam = [zhu[-1]]zuo_x = zuo[0:30]zuo_y = zuo[30:60]zuo_diam = [zuo[-1]]you_x = you[0:30]you_y = you[30:60]you_diam = [you[-1]]
############################################
#去掉主分支倒数后两个基本相等,所以决定去掉主分支最后一个，然后头部再补一个相等的数，使维度相等del zhu_x[-1]del zhu_y[-1]zhu_x.insert(0,zhu_x[0])zhu_y.insert(0,zhu_y[0])zhu_x.extend(zhu_y)zuo_x.extend(zuo_y)you_x.extend(you_y)zhu_x.extend(zhu_diam)zuo_x.extend(zuo_diam)you_x.extend(you_diam)fencha = [zhu_x] + [zuo_x] + [you_x]final.append(fencha)
finaldata = np.array(final)
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],color='red',linewidth=np.log(finaldata[i][0][-1]))
# plt.plot([finaldata[i][0][29]]+finaldata[i][1][0:30],[finaldata[i][0][59]]+finaldata[i][1][30:60],color='blue',linewidth=np.log(finaldata[i][1][-1]))
# plt.plot([finaldata[i][0][29]]+finaldata[i][2][0:30],[finaldata[i][0][59]]+finaldata[i][2][30:60],color='green',linewidth=np.log(finaldata[i][2][-1]))
# 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()依此判断两点之间的距离，然后在距离最大的两点之间插值
import numpy as np
import operator
import os
import copy
from matplotlib.font_manager import FontProperties
from scipy.interpolate import lagrange
import random
import matplotlib.pyplot as plt
import 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:1000]#抽取一部分
x1 = data[:,5]#x起点坐标
x2 = data[:,9]#x终点坐标
y1 = data[:,6]#y起
y2 = data[:,10]#y起
z1 = data[:,4]#IDpart
z2 = data[:,8]#IDpart
diam = 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 = lxqi
lxzg.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 = []
#去掉起点id
for 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==1
for i in range(len(datatoudeng)):a = datatoudeng[i]assert len(a[0])%2==1assert len(a[1])%2==1assert len(a[2])%2==1
finaldata = datatoudeng +dataweideng#未插值
final = datatoudeng #所有分叉，头等分叉，尾等分叉
##############################################################
#计算两点之间距离
def get_len(x1,x2,y1,y2):diff_x = (x1-x2)**2diff_y = (y1-y2)**2length = np.sqrt(diff_x+diff_y)return length
######################################################################
#插值方法三：在相邻两点坐标距离最大的地方插值
finaldata = []
for i in range(len(final)):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)]
# plt.plot(zhu_x,zhu_y,color='red')
# plt.plot(zuo_x,zuo_y,color='blue')
# plt.plot(you_x,you_y,color='green')while len(zhu_x)< 31:zhu_lin_list = []for j in range(1,len(zhu_x)):zhu_lin_len = get_len(zhu_x[j-1],zhu_x[j],zhu_y[j-1],zhu_y[j])zhu_lin_list.append(zhu_lin_len)zhu_max_index = zhu_lin_list.index(max(zhu_lin_list)) #j-1#不处理的话会出现nanif abs(zhu_x[zhu_max_index]-zhu_x[zhu_max_index+1])==0:zhu_x[zhu_max_index+1] = zhu_x[zhu_max_index+1] +1zhu_insert_x = np.linspace(zhu_x[zhu_max_index],zhu_x[zhu_max_index+1],3)#插入的点zhu_insert_x = zhu_insert_x[1]f_zhu = lagrange([zhu_x[zhu_max_index],zhu_x[zhu_max_index+1]],[zhu_y[zhu_max_index],zhu_y[zhu_max_index+1]])zhu_insert_y = f_zhu(zhu_insert_x)zhu_x.insert(zhu_max_index+1,zhu_insert_x)zhu_y.insert(zhu_max_index+1,zhu_insert_y)while len(zuo_x) < 31:zuo_lin_list = []for j in range(1,len(zuo_x)):zuo_lin_len = get_len(zuo_x[j-1],zuo_x[j],zuo_y[j-1],zuo_y[j])zuo_lin_list.append(zuo_lin_len)zuo_max_index = zuo_lin_list.index(max(zuo_lin_list)) #对应j-1if abs(zuo_x[zuo_max_index]-zuo_x[zuo_max_index+1])==0:zuo_x[zuo_max_index+1] = zuo_x[zuo_max_index+1] + 1zuo_insert_x = np.linspace(zuo_x[zuo_max_index],zuo_x[zuo_max_index+1],3)# #插入的点zuo_insert_x = zuo_insert_x[1]f_zuo = lagrange([zuo_x[zuo_max_index],zuo_x[zuo_max_index+1]],[zuo_y[zuo_max_index],zuo_y[zuo_max_index+1]])zuo_insert_y = f_zuo(zuo_insert_x)zuo_x.insert(zuo_max_index+1,zuo_insert_x)zuo_y.insert(zuo_max_index+1,zuo_insert_y)while len(you_x) < 31:you_lin_list = []for j in range(1,len(you_x)):#计算相邻两坐标的距离you_lin_len = get_len(you_x[j-1],you_x[j],you_y[j-1],you_y[j])#添加进列表中you_lin_list.append(you_lin_len)#计算距离最大的值对应的索引，对应x坐标的j-1和j之间的距离，最大you_max_index = you_lin_list.index(max(you_lin_list)) #对应j-1#然后，在两个最大点之间，平均插入一个数，作为插入x点if abs(you_x[you_max_index]-you_x[you_max_index+1]) == 0:you_x[you_max_index+1] = you_x[you_max_index+1] + 1you_insert_x = np.linspace(you_x[you_max_index],you_x[you_max_index+1],3)#插入的点you_insert_x = you_insert_x[1]#拉格朗日计算直线方程f_you = lagrange([you_x[you_max_index],you_x[you_max_index+1]],[you_y[you_max_index],you_y[you_max_index+1]])#插入的y点you_insert_y = f_you(you_insert_x)#将求得的x,y点插入对应位置you_x.insert(you_max_index+1,you_insert_x)you_y.insert(you_max_index+1,you_insert_y)
################################################
#可视化插值点
# plt.scatter(zhu_x,zhu_y,marker='*',color='red')
# plt.scatter(zuo_x,zuo_y,marker='*',color='blue')
# plt.scatter(you_x,you_y,marker='*',color='green')
# plt.show()
######################################################################处理头部相等,删除主分支最后一个，删除左右分支第一个del zhu_x[-1]del zhu_y[-1]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,rotatey
rotatedata = []
for i in range(0,360,360): #每隔3，30度旋转一次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)
################################################################
#归一化前可视化单张图片
# for i in range(len(finaldata)):
# plt.scatter(finaldata[i][0][0:30],finaldata[i][0][30:60],marker='*',color='red')
# plt.scatter(finaldata[i][1][0:30],finaldata[i][1][30:60],marker='*',color='blue')
# plt.scatter(finaldata[i][2][0:30],finaldata[i][2][30:60],marker='*',color='green')
# plt.show()
##################################################################
#归一化处理不去除管径信息
finalSubCAM = []
final = finaldata
for 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],color='red',linewidth=np.log(finaldata[i][0][-1]))
# plt.plot([finaldata[i][0][29]]+finaldata[i][1][0:30],[finaldata[i][0][59]]+finaldata[i][1][30:60],color='blue',linewidth=np.log(finaldata[i][1][-1]))
# plt.plot([finaldata[i][0][29]]+finaldata[i][2][0:30],[finaldata[i][0][59]]+finaldata[i][2][30:60],color='green',linewidth=np.log(finaldata[i][2][-1]))
# 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()