PythonでCTのDICOMデータをサーフェスレンダリングする方法（骨をobjで出力まで）

import pydicom
import numpy as np
import matplotlib.pyplot as plt
import sys
import glob

# DICOM ファイルを読み込み
files = []
print('DICOMファイルの場所: {}'.format("C:\\Users\\patient_name/*"))
for fname in glob.glob("C:\\Users\\patient_name/*", recursive=False):
    files.append(pydicom.dcmread(fname))

# slicesというリストを作成。ファイルに、slicelocationという属性があれば追加していく。
slices = []
skipcount = 0
for f in files:
    if hasattr(f, 'SliceLocation'):
            if f.SeriesNumber==5: #5回CTを撮影し、5回目のデータのみを抽出
                slices.append(f)
    else:
        skipcount = skipcount + 1

print("該当ファイル数: {}".format(len(slices)))
print("スライス位置がないファイル: {}".format(skipcount))

# スライスの順番をそろえる
slices = sorted(slices, key=lambda s: s.SliceLocation) # .SliceLocationで場所情報を取得
# len(slices) = 270

# アスペクト比を計算する
ps = slices[0].PixelSpacing    # ps = [0.571, 0.571] 1ピクセルの [y, x] 長さ
ss = slices[0].SliceThickness  # ss = "3.0" スライスの厚み
ax_aspect = ps[0]/ps[1]        # yの長さ/xの長さ =  1
sag_aspect = ss/ps[0]          # スライスの厚み/y軸方向への１ピクセルの長さ = 3/0.571 = 5.25
cor_aspect = ss/ps[1]          # スライスの厚み/x軸方向への１ピクセルの長さ = 3/0.571 = 5.25

# 空の3Dのnumpy arrayを作成する 
img_shape = list(slices[0].pixel_array.shape) # img_shape = [512, 512]
img_shape.insert(0,len(slices)) # img_shape = [270, 512, 512]  
img3d = np.zeros(img_shape) # 空のarrayを作る　np.zeros([270, 512, 512])

# 3Dのnumpy arrayを作る
for i, s in enumerate(slices):
    img2d = s.pixel_array    # img2d.shape = (512, 512) 
    img3d[i,:, :] = img2d    # img3d.shape = (270, 512, 512)

# プロットする
middle0 = img_shape[0]//2  #3d array の(z軸)頭尾軸中央　270/2 = 135
middle1 = img_shape[1]//2  #3d array の(y軸)上下軸中央　512/2 = 256
middle2 = img_shape[2]//2  #3d array の(x軸)左右軸中央　512/2 = 256

a1 = plt.subplot(1, 3, 1)
plt.imshow(img3d[middle0,:, :])
a1.set_aspect(ax_aspect)  # ax_aspect = 1

a2 = plt.subplot(1, 3, 2)
plt.imshow(img3d[:, middle1, :])
a2.set_aspect(cor_aspect) # cor_aspect = 5.25

a3 = plt.subplot(1, 3, 3)
plt.imshow(img3d[:, :, middle2]) # 
a3.set_aspect(sag_aspect) # sag_aspect = 5.25

plt.show()

np.amin(img3d)

np.amax(img3d)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import measure

# マーチングキュープ法を実施
verts, faces, normals, values = measure.marching_cubes(img3d, 
                                                     level=299,
                                                     spacing=(3,0.57, 0.57))

# Matplotlibを使って表示します。
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')

mesh = Poly3DCollection(verts[faces])
mesh.set_edgecolor('k')
ax.add_collection3d(mesh)

ax.set_xlabel("x-axis")
ax.set_ylabel("y-axis")
ax.set_zlabel("z-axis") 

ax.set_xlim(0, 800)  
ax.set_ylim(0, 800)  
ax.set_zlim(0, 800)  

plt.tight_layout()
plt.show()

import copy

#　DICOMのデータが入ったimg3dをコピーし、名前をimg3d_thresholdとします。
img3d_threshold = copy.deepcopy(img3d)
#　img3dと同じ形の0だけのnumpy arrayを作ります。dtype=boolとしてください。
img3d_boolean = np.zeros(img_shape, dtype=bool)

#　CT値が300未満の場合
indice_under = img3d_threshold < 300
# img3d_thresholdに0を入れます。
img3d_threshold[indice_under] = 0
# img3d_booleanにfalseを入れます。
img3d_boolean[indice_under] = False

#　CT値が300以上の場合
indice_over = img3d_threshold >= 300
# img3d_thresholdに2を入れます。
img3d_threshold[indice_over] = 2
# img3d_booleanにtrueを入れます。
img3d_boolean[indice_over] = True

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import measure

# マーチングキュープ法を実施
verts, faces, normals, values = measure.marching_cubes(img3d_threshold, 
                                                      level=1,
                                                      spacing=(3,0.57, 0.57),
                                                      mask=img3d_boolean)

# グラフを作成
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')

mesh = Poly3DCollection(verts[faces])
mesh.set_edgecolor('k')
ax.add_collection3d(mesh)

ax.set_xlabel("x-axis")
ax.set_ylabel("y-axis")
ax.set_zlabel("z-axis")

ax.set_xlim(0, 800)  
ax.set_ylim(0, 800)  
ax.set_zlim(0, 800)  

plt.tight_layout()
plt.show()

faces = faces+1

thefile = open('DICOM_bone.obj', 'w')

for item in verts:
  thefile.write("v {0} {1} {2}\n".format(item[0],item[1],item[2]))

for item in normals:
  thefile.write("vn {0} {1} {2}\n".format(item[0],item[1],item[2]))

for item in faces:
  thefile.write("f {0}//{0} {1}//{1} {2}//{2}\n".format(item[0],item[1],item[2]))  

thefile.close()