三维体素图绘制实验与教程

还是实验课内容，数据是Transunet的标签数据存成了npy格式，想把多器官分割结果用立体的方式表现出来，完整代码在最后，前面主要是搞的一些测试。体素的绘图具体参考了链接【1】，主要工具是python matplotlib包。matplotlib的一大主要问题就是大数据画图慢，因而下面体素的流程演示的数据采用【1】中的随机数生成。

xyzvalues = np.random.choice(range(0,10), size=(10,10,10))

首先结论是在大数据量时，散点图是绝对优于体素的，核心原因就是快，而在数据较少时，通常体素图画着比较快乐，散点图可能会分层。
另外勾线细一点比较好看。

三维器官效果展示

• 体素：
  
• 散点：
  
  

详细流程

1.体素绘图

• 核心命令：

fig = plt.figure(figsize=(7, 4.5))
ax = fig.gca(projection='3d')
ax.voxels(xyzvalues, facecolors=colorsvalues, edgecolor='k',shade=False,)

结果图如下，voxels的更详细用法见【4】：

edgecolor=None时效果如下：

• 颜色表选取

mycolormap = plt.get_cmap('coolwarm')
xyzmaxvalue=xyzvalues.max()
relativevalue=np.round(xyzvalues/xyzmaxvalue,1)

plt提供了多种颜色表，具体见链接【5】

另外注意用于生成颜色表的数据必须归一化，不然会变成这样：

• 不透明颜色表生成

colorsvalues=mycolormap(relativevalue)

• 透明颜色表生成

sizeall=xyzvalues.size
zt=np.reshape(relativevalue,(sizeall,))
alpha=0.5
colorsvalue=np.array([(mycolormap(i)[0],mycolormap(i)[1],mycolormap(i)[2],alpha) for i in zt])
colorsvalues=np.reshape(colorsvalue,(d,w,h,4))

增加透明度，减轻遮挡，效果如下：

2.散点绘图
因为体素画着太慢，在此基础上我们尝试散点绘图.散点绘图速度快了许多，而且能鼠标拖动转了，就是有些时候效果非常的诡异。

• 核心代码

ax.scatter(YL,XL,ZL,s=5, c=ZTL, cmap='coolwarm',edgecolors="k",linewidth=0.3,depthshade=True,alpha=0.3)

这个没什么能讲的，主要注意的是画黑边，大半径比较好看，透明度点多的时候其实没什么效果

• 坐标变换
  先把数据的深度提到第三维里：

lab=np.load("TransUNet/testre/case0001_predicyion.npy")
labb=np.transpose(lab,(2,1,0))

提取绘图范围：

xlim=[200,320]
ylim=[180,300]
zlim=[95,105]

xyzvalues=labb[xlim[0]:xlim[1],ylim[0]:ylim[1],zlim[0]:zlim[1]]

生成坐标：

x=np.linspace(xlim[0],xlim[1]-1,-xlim[0]+xlim[1])
y=np.linspace(ylim[0],ylim[1]-1,-ylim[0]+ylim[1])
z=np.linspace(zlim[0],zlim[1]-1,-zlim[0]+zlim[1])
ZALL=zlim[1]-zlim[0]

XX,YY=np.meshgrid(x,y)
X=np.repeat(XX, ZALL,axis=1)
X=np.reshape(X,(d,w,h))
Y=np.repeat(YY, ZALL,axis=1)
Y=np.reshape(Y,(d,w,h))
zzz=np.ones((d,w,h))
Z=zzz*z

X=np.reshape(X,(-1,))
Y=np.reshape(Y,(-1,))
Z=np.reshape(Z,(-1,))

与体素不同，散点需要把没有标签的数据剔除，不然也会上色。

XL=[];YL=[];ZL=[];ZTL=[]
for i in range(sizeall):
    if zt[i]>0:
        XL.append(X[i])
        YL.append(Y[i])
        ZL.append(Z[i])
        ZTL.append(zt[i])

这部分代码写得时候太困了，应该有简化版本，但是懒了.jpg。

• 不同大小的点的诡异效果：
  在选取的点数量较少的时候，散点图绘制会出现一些问题，首先看一眼体素效果：
  
  小半径散点：
  
  大半径散点：
  
  所以说这种画法存在一个问题，小半径z轴连不起来，很明显一层一层的，大半径又膨胀太厉害。应当是因为z轴的刻度大，xy刻度小引起的，理论上应当可以通过改变坐标轴显示范围和刻度间隔来缓解，或者干脆选点的时候选成立方体范围。

完整代码

1.体素式绘图

lab=np.load("TransUNet/testre/case0001_predicyion.npy")
labb=np.transpose(lab,(2,1,0))
xlim=[200,320]
ylim=[180,300]
zlim=[95,105]
xyzvalues=labb[xlim[0]:xlim[1],ylim[0]:ylim[1],zlim[0]:zlim[1]]

mycolormap = plt.get_cmap('viridis')
xyzminvalue=xyzvalues.min()
xyzmaxvalue=xyzvalues.max()

relativevalue=np.round(xyzvalues/xyzmaxvalue,1)
colorsvaluesr = np.empty(xyzvalues.shape, dtype=object)
alpha=0.5
d,w,h=xyzvalues.shape
sizeall=xyzvalues.size
zt=np.reshape(relativevalue,(sizeall,))



colorsvalue=np.array([(mycolormap(i)[0],mycolormap(i)[1],mycolormap(i)[2],alpha) for i in zt])
colorsvalues=np.reshape(colorsvalue,(d,w,h,4))


fig = plt.figure(figsize=(7, 4.5))
ax = fig.gca(projection='3d')
pos=ax.voxels(xyzvalues, facecolors=colorsvalues, edgecolor='k',linewidth=0.3,shade=False,)



ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Voxel Map')

plt.show()

2.散点式绘图

lab=np.load("TransUNet/testre/case0001_predicyion.npy")
labb=np.transpose(lab,(2,1,0))
xlim=[200,320]
ylim=[180,300]
zlim=[95,105]
xyzvalues=labb[xlim[0]:xlim[1],ylim[0]:ylim[1],zlim[0]:zlim[1]]

xyzminvalue=xyzvalues.min()
xyzmaxvalue=xyzvalues.max()

relativevalue=np.round(xyzvalues/xyzmaxvalue,1)
alpha=0.5
d,w,h=xyzvalues.shape
sizeall=xyzvalues.size
zt=np.reshape(relativevalue,(sizeall,))

fig = plt.figure(figsize=(7, 4.5))# Make a figure and axes with dimensions as desired.
ax = fig.gca(projection='3d')

x=np.linspace(xlim[0],xlim[1]-1,-xlim[0]+xlim[1])
y=np.linspace(ylim[0],ylim[1]-1,-ylim[0]+ylim[1])
z=np.linspace(zlim[0],zlim[1]-1,-zlim[0]+zlim[1])
ZALL=zlim[1]-zlim[0]

XX,YY=np.meshgrid(x,y)
X=np.repeat(XX, ZALL,axis=1)
X=np.reshape(X,(d,w,h))
Y=np.repeat(YY, ZALL,axis=1)
Y=np.reshape(Y,(d,w,h))
zzz=np.ones((d,w,h))
Z=zzz*z

X=np.reshape(X,(-1,))
Y=np.reshape(Y,(-1,))
Z=np.reshape(Z,(-1,))
XL=[];YL=[];ZL=[];ZTL=[]
for i in range(sizeall):
    if zt[i]>0:
        XL.append(X[i])
        YL.append(Y[i])
        ZL.append(Z[i])
        ZTL.append(zt[i])
ax.scatter(YL,XL,ZL,s=5, c=ZTL, cmap='coolwarm',edgecolors="k",linewidth=0.3,depthshade=True,alpha=0.3)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Voxel Map')

plt.show()

参考

【1】https://blog.youkuaiyun.com/weixin_39771351/article/details/111293632
【2】https://blog.youkuaiyun.com/qq_31225201/article/details/95619527
【3】https://blog.youkuaiyun.com/qq_37851620/article/details/100642566
【4】https://www.jianshu.com/p/03ed963eb51a
【5】https://blog.youkuaiyun.com/lly1122334/article/details/88535217
【6】https://blog.youkuaiyun.com/weixin_44520259/article/details/89917026