TSDF文件mhd,raw格式转点云的可视化，python实现

最近获得了一个TSDF表示的三维模型，文件格式是mhd和raw文件，我看很多医学图像也是用的这个格式，但是网上却找不到可视化的方法，于是在研究了TSDF的原理之后，自己用Python进行了TSDF到点云的简单转换，有一定的误差，但是效果还是可以的。

TSDF原理：

TSDF其实与体素表示提相近的，但是又有一定的差别。TSDF是固定了一个三维空间，如图1所示，这个三维空间有点类似于体素，空间设定了固定的长宽高，有一个固定的体积，和体素空间非常相似，但是与体素不同的是，这个空间里面每一个体素位置都会有一个值，这个值是每一个体素点到这个空间里面的对象的表面的距离，而且TSDF只进行这个三维空间的场景建模，如果建模过程中建模的范围超出了这个三维空间，那么超出的部分是不会记录到模型中的。如图2所示，是一个人脸的切面，也就是TSDF的其中一层的演示。因为只涉及到到表面的距离，在表面前面的值是正的，如果到表面后面了这个值就是负的，默认值是零，如果值的绝对值大于1，则用正负整数1来表示，从图2中可以看出来。

图1 TSDF的三维空间

图2 TSDF的一个切片的演示

那么如果要可视化TSDF，通常是使用光线射影raycasting，就是从你想要查看的视角射出光线，如果穿过了表面，那么就记录，如果没有就不记录值。

但是我目前还没有找到很好用的，于是我打算把TSDF转换成点云文件，然后进行查看。

我的做法是记录每一个边缘位置的体素的位置为点云坐标，比如图2里面展示的，从正值到负值的转换边界记录成点云坐标，这与实际的场景模型的坐标值是不对应的，所以目前只适合进行可视化，如果要得到正确的点云值还需要进行转换，我之后会考虑实现，目前先记录我把TSDF转换成可视化点云的过程。

我操作的是7scenes数据集里面提供的TSDF模型：RGB-D Dataset 7-Scenes - Microsoft Research

使用的是chess场景的raw文件。

我首先通过二进制方式读取了raw文件，并将二进制按照有符号short型进行转换，然后转变为numpy数组，并转为了512x512x512的数组。

接下来我对这些数组进行遍历，记录每一个维度中发生了符号变化的所有点边界，并存储到点云文件中。

最后读取点云文件，使用open3开源3D库进行可视化。

使用需安装：

pip install open3d
pip install numpy

下面是代码：

import numpy as np
import struct

TSDFfilename='chess.raw' # TSDF 二进制文件路径
PointCloudfilename="chessPointCloud.xyz" # 输出点云文件
h,w,d=512,512,512 # TSDF 的高度 宽度和深度

# 加载测试数据
f = open(TSDFfilename,'rb')
# 文档中包含的数字个数，而一个short数占2个字节
data_raw = struct.unpack('h'*(h*w*d),f.read(2*h*w*d)) # h: short H: usight short
f.close()
verify_data =  np.asarray(data_raw).reshape(h,w,d)

pointCloud = []
f_pc = open(PointCloudfilename, "w")

# 转换为点云的遍历过程
for i in range(h):
    for j in range(w):
        for k in range(1,d):
            if (verify_data[i][j][k]>0 and verify_data[i][j][k-1]<0):
                # print(i,j,k)
                pointCloud.append([i,j,k])
                f_pc.write(str(i)+" "+str(j)+" "+str(k)+" \n")

for i in range(h):
    for k in range(d):
        for j in range(1,w):
            if (verify_data[i][j][k]>0 and verify_data[i][j-1][k]<0):
                # print(i,j,k)
                if [i,j,k] in pointCloud: # 避免记录同样的点
                    pass
                else:
                    f_pc.write(str(i)+" "+str(j)+" "+str(k)+" \n")

for k in range(d):
    for j in range(w):
        for i in range(1,h):
            if (verify_data[i][j][k]>0 and verify_data[i-1][j][k]<0):
                # print(i,j,k)
                if [i,j,k] in pointCloud: # 避免记录同样的点
                    pass
                else:
                    f_pc.write(str(i)+" "+str(j)+" "+str(k)+" \n")
f_pc.close()

#下面部分的代码可以单独拿出去用作可视化，改一下文件名那里的参数就可以了
import open3d as o3d
 
print("Load a ply point cloud, print it, and render it")
pcd = o3d.io.read_point_cloud(PointCloudfilename)#可以替换成需要查看的点云文件
o3d.visualization.draw_geometries([pcd])

这是生成的点云的可视化效果（代码运行的时间挺长的）：

By the way, to open .mhd/.raw images you can use multiple software:
3D Slicer image computing platform | 3D Slicer
vv - Rayonnement, Images, Oncologie