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:50]#抽取一部分
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起点id
b1qi_id = [b5]+b1 +[b2]
b6 = z2[s1[i1+1]-1]#x,y终点id
b1zg_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==1
assert len(a[1])%2==1
assert len(a[2])%2==1
for i in range(len(datatoudeng)):
a = datatoudeng[i]
assert len(a[0])%2==1
assert len(a[1])%2==1
assert len(a[2])%2==1
finaldata = datatoudeng +dataweideng#未插值
final = datatoudeng #所有分叉,头等分叉,尾等分叉
# print(final)
###################################################
#未插值前,先对长度不相等的列表进行旋转,以便后面求得最大最小值
#finalrotate保存了原来分叉的旋转若干角度的列表
# for file in finalrotate:
#这里的file类似于上面的final保存了每个分叉
#对于此处旋转通用
def rotate(angle,valuex,valuey):
rotatex = math.cos(angle)*valuex -math.sin(angle)*valuey
rotatey = math.cos(angle)*valuey + math.sin(angle)* valuex
#做转换为列表
rotatex = rotatex.tolist()
rotatey = rotatey.tolist()
xy = rotatex + rotatey
return xy
finalrotate = []
for i in range(len(final)):
rotatedata = []
for j in range(0,360,3):#每隔10,3,5度旋转一次
zhu = final[i][0]
zuo = final[i][1]
you = final[i][2]
zhu = np.array(zhu)
zuo = np.array(zuo)
you = np.array(you)
rotatezhu = rotate(j,zhu[0:len(zhu)//2],zhu[(len(zhu)//2):(len(zhu)-1)])
rotatezuo = rotate(j,zuo[0:len(zuo)//2],zuo[(len(zuo)//2):(len(zuo)-1)])
rotateyou = rotate(j,you[0:len(you)//2],you[(len(you)//2):(len(you)-1)])
rotatezhu.append(zhu[-1])
rotatezuo.append(zuo[-1])
rotateyou.append(you[-1])
rotatefencha = [rotatezhu]+[rotatezuo]+[rotateyou]
rotatedata.append(rotatefencha)
finalrotate.append(rotatedata)
##################################################
#插值方法:第一种策略,在头部相等分叉,主分支头部加入做小值,左右分支尾部加入最小值,最终使它们维度相等
final_list = []
for final in finalrotate:
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)]
#将一个分叉的三个分支的x整成一个列表,然后求它们的最大最小值,以便后面归一化
x_list = zhu_x.copy()
x_list.extend(zuo_x)
x_list.extend(you_x)
min_x = min(x_list)
max_x = max(x_list)
#将一个分叉的三个分支的y整成一个列表,然后求它们的最大最小值,以便后面归一化
y_list = zhu_y.copy()
y_list.extend(zuo_y)
y_list.extend(you_y)
min_y = min(y_list)
max_y = max(y_list)
#从这里开始两种策略出现不同,因为插入的值不同
#对于头部相等,主分支,在前面索引0插入n个x_list最小值,对于左右分支,在尾部-1,插入m个y_list最小值,最终使它们归一化后值为0
#对于final37希望每一个维度为60,所以,让他们每一个列表插入(30-len_list)个最小值
for i in range((30-len(zhu_x))):
zhu_x.insert(0,min_x)
zhu_y.insert(0,min_y)
#左分支在尾部插入最小值,左右分支头部要去掉一个,所以应比主分支多插入一个
for i in range((31-len(zuo_x))):
zuo_x.append(min_x)
zuo_y.append(min_y)
for i in range((31-len(you_x))):
you_x.append(min_x)
you_y.append(min_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_y
zuo_xy = zuo_x + zuo_y
you_xy = you_x + you_y
#这里在将坐标点与管径接起来
zhu = zhu_xy + zhu_diam
zuo = zuo_xy + zuo_diam
you = you_xy + you_diam
fencha = [zhu]+[zuo]+[you]
finaldata.append(fencha)
final = np.array(finaldata)
final_list.append(final)
final = np.array(final_list)
final = final.reshape(-1,3,61)
# print(final.shape)
#################################################################################
#归一化处理不去除管径信息
finalSubCAM = []
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,或归一化管径normd
finalSubCAM.append(normxy)
finaldata = np.array(finalSubCAM)
np.random.shuffle(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()
#############################################################
#遍历分叉,去除与原点相连的坐标
# finaldata = finaldata.tolist()
# for i in range(len(finaldata)):
# zhu = finaldata[i][0]
# zuo = finaldata[i][1]
# you = finaldata[i][2]
# zhu_x = zhu[0:30]
# zhu_y = zhu[30:60]
# zhu_dim =zhu[-1]
# zuo_x = zuo[0:30]
# zuo_y = zuo[30:60]
# zuo_dim = zuo[-1]
# you_x = you[0:30]
# you_y = you[30:60]
# you_dim = you[-1]
# zhu_pltx = []
# zhu_plty = []
# for j in range(len(zhu_x)):
# if zhu_x[j] != zhu_y[j]:
# zhu_pltx.append(zhu_x[j])
# zhu_plty.append(zhu_y[j])
# zuo_pltx = []
# zuo_plty = []
# for j in range(len(zuo_x)):
# if zuo_x[j] != zuo_y[j]:
# zuo_pltx.append(zuo_x[j])
# zuo_plty.append(zuo_y[j])
# you_pltx = []
# you_plty = []
# for j in range(len(you_x)):
# if you_x[j] != you_y[j]:
# you_pltx.append(you_x[j])
# you_plty.append(you_y[j])
# plt.figure(figsize=(128,128),dpi=1)
# width0 = np.log(zhu_dim)
# width1 = np.log(zuo_dim)
# width2 = np.log(you_dim)
# plt.plot(zhu_pltx,zhu_plty,color='red',linewidth=width0*100)
# plt.plot([zhu_pltx[-1]]+zuo_pltx,[zhu_plty[-1]]+zuo_plty,color='blue',linewidth=width1*100)
# plt.plot([zhu_pltx[-1]]+you_pltx,[zhu_plty[-1]]+you_plty,color='green',linewidth=width2*100)
# plt.xticks(np.arange(0,1,0.1))
# plt.yticks(np.arange(0,1,0.1))
# # plt.axis('off')
# 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()