1理论基础

使用DCM三维重建格式文件。DCM格式，可参考医疗CT断层图像DICOM标准和工业CT断层图像DICONDE标准。 在使用DCM三维重建格式文件时，实际上是从几个方面进行的DCM在文件中，获取三维结构信息。然后用面绘或体绘重建三维图像。

2 代码（python）

prepro.py

import matplotlib.pyplot as plt import pydicom.uid import sys from PyQt5 import QtGui import os import pydicom import glob from PIL import * import matplotlib.pyplot as plt from pylab import * from tkinter.filedialog import * import numpy import SimpleITK as sitk   def window(img):     win_min = 350     win_max = 50      for i in range(img.shape[0]):         img[i] = 255.0 * (img[i] - win_min) / (win_max - win_min)         min_index = img[i] < 0         img[i][min_index] = 0         max_index = img[i] > 255         img[i][max_index] = 255         img[i] = img[i] - img[i].min()         c = float(255) / img[i].max()         img[i] = img[i] * c      return img.astype(np.uint8)  def get_pixels_hu(ct_array,slope,intercept):     if slope != 1:         ct_array = slope*ct_array     ct_array  = intercet
    return ct_array

PathDicom = "D:/xiangmu/CT/data/CQ500_CT/CQ500CT181 CQ500CT181/Unknown Study/CT 0.625mm"
# PathDicom = r'D:\xiangmu\CT\data\CQ500_CT\CQ500CT0 CQ500CT0\Unknown Study\CT 4cc sec 150cc D3D on'
lstFilesDCM = []
SaveRawDicom = "SaveRaw/"

reader=sitk.ImageSeriesReader()
series_IDs=sitk.ImageSeriesReader.GetGDCMSeriesIDs(PathDicom)

nb_series=len(series_IDs)
# print(nb_series)

series_file_names=sitk.ImageSeriesReader.GetGDCMSeriesFileNames(PathDicom,series_IDs[0])
reader.SetFileNames(series_file_names)
image=reader.Execute()
# median 滤波
sitk_median = sitk.MedianImageFilter()
sitk_median.SetRadius(2)
image = sitk_median.Execute(image)
# sitk.WriteImage(sitk_median, 'sitk_median.mha')

#增强


print('before:',image)
# for root,dirs,files in os.walk(PathDicom):
# for file in files:
# print(os.path.join(root,file))
# image=sitk.ReadImage(os.path.join(root,file))
resample=sitk.ResampleImageFilter()
resample.SetInterpolator(sitk.sitkLinear)
# resample.SetInterpolator(sitk.sitkNearestNeighbor) #没有线性插值效果好
# resample.SetInterpolator(sitk.sitkBSpline) #效果不错
# resample.SetInterpolator(sitk.sitkGaussian) #效果不错，耗时较长
resample.SetDefaultPixelValue(0)
newsapcing=[0.488281,0.488281,0.625]  #每个维度像素的间距
# newsapcing=[0.488281,0.488281,0.488281]
resample.SetOutputSpacing(newsapcing)
resample.SetOutputOrigin(image.GetOrigin())     #源点
resample.SetOutputDirection(image.GetDirection())   #每个轴的方向

size=[512,512,330]
resample.SetSize(size)
resample.SetDefaultPixelValue(0)
new=resample.Execute(image)
print('after:',new)
# data=sitk.GetArrayFromImage(new)
sitk_hisequal = sitk.AdaptiveHistogramEqualizationImageFilter()
sitk_hisequal.SetAlpha(0.9)
sitk_hisequal.SetBeta(0.9)
sitk_hisequal.SetRadius(3)
sitk_hisequal = sitk_hisequal.Execute(new)
sitk.WriteImage(sitk_hisequal, os.path.join(SaveRawDicom, "sample_test" + ".mhd"))

# print(image)
# for dirName, subdirList, fileList in os.walk(PathDicom): #os.walk目录遍历器 根目录，文件夹，文件
# for filename in fileList:
# if ".dcm" in filename.lower(): # 判断文件是否为dicom文件
# lstFilesDCM.append(os.path.join(dirName, filename)) # 加入到列表中
#
# for i, filenameDCM in enumerate(lstFilesDCM):
# ds = sitk.ReadImage(filenameDCM)
# image_array = sitk.GetArrayFromImage(ds)
# # image_array = window(image_array)
# # ArrayDicom[:, :, i] = np.asarray(Image.fromarray(image_array[0]))
# # out = sitk.GetImageFromArray(image_array)
# ct_array=get_pixels_hu(image_array)
# import vtk
# def show(fileName):
# colors = vtk.vtkNamedColors()
#
# colors.SetColor("SkinColor", [255, 125, 64, 255])
# # colors.SetColor("SkinColor", [0, 255, 0, 200])
# colors.SetColor("BkgColor", [51, 77, 102, 255])
#
#
# aRenderer = vtk.vtkRenderer() #渲染器
# renWin = vtk.vtkRenderWindow() #渲染窗口
# renWin.AddRenderer(aRenderer)
# iren = vtk.vtkRenderWindowInteractor() #获取渲染窗口上发生的鼠标，键盘，事件事件。
# iren.SetRenderWindow(renWin)
#
#
# reader = vtk.vtkMetaImageReader() #读取、显示 .mha 图像
# reader.SetFileName(fileName)
#
#
# skinExtractor = vtk.vtkMarchingCubes()
# skinExtractor.SetInputConnection(reader.GetOutputPort())
# skinExtractor.SetValue(0, 500) #500的等值面对应于患者的皮肤，
#
# skinStripper = vtk.vtkStripper() #三角形剥离
# skinStripper.SetInputConnection(skinExtractor.GetOutputPort())
#
# skinMapper = vtk.vtkPolyDataMapper() #vtkPolyData数据显示时需要定义vtkPolyDataMapper对象，用来接受vtkPolyData数据以实现图形数据到渲染图元的转换。
# skinMapper.SetInputConnection(skinStripper.GetOutputPort())
# skinMapper.ScalarVisibilityOff() #阻止颜色映射
#
# skin = vtk.vtkActor()
# skin.SetMapper(skinMapper)
# skin.GetProperty().SetDiffuseColor(colors.GetColor3d("SkinColor"))
# skin.GetProperty().SetSpecular(.5) #光照
# skin.GetProperty().SetSpecularPower(10) #20
# skin.GetProperty().SetOpacity(.5) #设置透明度 .5
# # skin.GetProperty().SetColorsLevel(350) #加
# # skin.GetProperty().SetColorsWindow(50) #加
#
#
# boneExtractor = vtk.vtkMarchingCubes()
# boneExtractor.SetInputConnection(reader.GetOutputPort())
# boneExtractor.SetValue(0, 1150) #1150对应于患者的骨骼
# # boneExtractor.SetValue(0, 1000) # 1150对应于患者的骨骼
#
# boneStripper = vtk.vtkStripper()
# boneStripper.SetInputConnection(boneExtractor.GetOutputPort())
#
# boneMapper = vtk.vtkPolyDataMapper()
# boneMapper.SetInputConnection(boneStripper.GetOutputPort())
# boneMapper.ScalarVisibilityOff()
#
# bone = vtk.vtkActor()
# bone.SetMapper(boneMapper)
# bone.GetProperty().SetDiffuseColor(colors.GetColor3d("Ivory")) #象牙色
# bone.GetProperty().SetOpacity(.8) #加
# # bone.GetProperty().SetColorsLevel(2000) #加
# # bone.GetProperty().SetColorsWindow(500) #加
#
#
# outlineData = vtk.vtkOutlineFilter()
# outlineData.SetInputConnection(reader.GetOutputPort())
#
# mapOutline = vtk.vtkPolyDataMapper()
# mapOutline.SetInputConnection(outlineData.GetOutputPort())
#
# outline = vtk.vtkActor()
# outline.SetMapper(mapOutline)
# outline.GetProperty().SetColor(colors.GetColor3d("Black"))
#
# aCamera = vtk.vtkCamera()
# aCamera.SetViewUp(0, 0, -1)
# aCamera.SetPosition(0, -1, 0)
# aCamera.SetFocalPoint(0, 0, 0)
# aCamera.ComputeViewPlaneNormal()
# aCamera.Azimuth(30.0)
# aCamera.Elevation(30.0)
#
# aRenderer.AddActor(outline)
# aRenderer.AddActor(skin)
# aRenderer.AddActor(bone)
# aRenderer.SetActiveCamera(aCamera)
# aRenderer.ResetCamera()
# aCamera.Dolly(1.5)
#
# aRenderer.SetBackground(colors.GetColor3d("BkgColor"))
# renWin.SetSize(640, 480)
# # renWin.SetSize(240, 100)
#
# aRenderer.ResetCameraClippingRange()
#
# iren.Initialize()
# iren.Start()
#
# show('SaveRaw/sample.mhd')

先运行prepro.py，生成三维信息，然后保存。之后再分别运行prepro_mian.py（使用面绘制重建出三维图像）或运行prepro_ti.py（使用体绘制重建出三维图像）。两种方法根据需要选择。 prepro_mian.py

import vtk
def show(fileName):
    colors = vtk.vtkNamedColors()

    colors.SetColor("SkinColor", [255, 125, 64, 255])
    # colors.SetColor("SkinColor", [0, 255, 0, 200])
    colors.SetColor("BkgColor", [51, 77, 102, 255])


    aRenderer = vtk.vtkRenderer()           #渲染器
    renWin = vtk.vtkRenderWindow()          #渲染窗口
    renWin.AddRenderer(aRenderer)
    iren = vtk.vtkRenderWindowInteractor()          #获取渲染窗口上发生的鼠标，键盘，事件事件。
    iren.SetRenderWindow(renWin)


    reader = vtk.vtkMetaImageReader()   #读取、显示 .mha 图像
    reader.SetFileName(fileName)


    skinExtractor = vtk.vtkMarchingCubes()
    skinExtractor.SetInputConnection(reader.GetOutputPort())
    skinExtractor.SetValue(0, 500)              #500的等值面对应于患者的皮肤，

    skinStripper = vtk.vtkStripper()        #三角形剥离
    skinStripper.SetInputConnection(skinExtractor.GetOutputPort())

    skinMapper = vtk.vtkPolyDataMapper()    #vtkPolyData数据显示时需要定义vtkPolyDataMapper对象，用来接受vtkPolyData数据以实现图形数据到渲染图元的转换。
    skinMapper.SetInputConnection(skinStripper.GetOutputPort())
    skinMapper.ScalarVisibilityOff()    #阻止颜色映射

    skin = vtk.vtkActor()
    skin.SetMapper(skinMapper)
    skin.GetProperty().SetDiffuseColor(colors.GetColor3d("SkinColor"))
    skin.GetProperty().SetSpecular(.5)      #光照
    skin.GetProperty().SetSpecularPower(10)     #20
    skin.GetProperty().SetOpacity(.5)       #设置透明度 .5
    # skin.GetProperty().SetColorsLevel(350) #加
    # skin.GetProperty().SetColorsWindow(50) #加


    boneExtractor = vtk.vtkMarchingCubes()
    boneExtractor.SetInputConnection(reader.GetOutputPort())
    boneExtractor.SetValue(0, 1150)         #1150对应于患者的骨骼
    # boneExtractor.SetValue(0, 1000) # 1150对应于患者的骨骼

    boneStripper = vtk.vtkStripper()
    boneStripper.SetInputConnection(boneExtractor.GetOutputPort())

    boneMapper = vtk.vtkPolyDataMapper()
    boneMapper.SetInputConnection(boneStripper.GetOutputPort())
    boneMapper.ScalarVisibilityOff()

    bone = vtk.vtkActor()
    bone.SetMapper(boneMapper)
    bone.GetProperty().SetDiffuseColor(colors.GetColor3d("Ivory"))      #象牙色
    bone.GetProperty().SetOpacity(.8)               #加
    # bone.GetProperty().SetColorsLevel(2000) #加
    # bone.GetProperty().SetColorsWindow(500) #加


    outlineData = vtk.vtkOutlineFilter()
    outlineData.SetInputConnection(reader.GetOutputPort())

    mapOutline = vtk.vtkPolyDataMapper()
    mapOutline.SetInputConnection(outlineData.GetOutputPort())

    outline = vtk.vtkActor()
    outline.SetMapper(mapOutline)
    outline.GetProperty().SetColor(colors.GetColor3d("Black"))

    aCamera = vtk.vtkCamera()
    aCamera.SetViewUp(0, 0, -1)
    aCamera.SetPosition(0, -1, 0)
    aCamera.SetFocalPoint(0, 0, 0)
    aCamera.ComputeViewPlaneNormal()
    aCamera.Azimuth(30.0)
    aCamera.Elevation(30.0)

    aRenderer.AddActor(outline)
    aRenderer.AddActor(skin)
    aRenderer.AddActor(bone)
    aRenderer.SetActiveCamera(aCamera)
    aRenderer.ResetCamera()
    aCamera.Dolly(1.5)

    aRenderer.SetBackground(colors.GetColor3d("BkgColor"))
    renWin.SetSize(640, 480)
    # renWin.SetSize(240, 100)

    aRenderer.ResetCameraClippingRange()

    iren.Initialize()
    iren.Start()

show('SaveRaw/sample.mhd')

prepro_ti.py

import vtk

# file_path=r'D:/xiangmu/CT/data/CQ500_CT/CQ500CT181 CQ500CT181/Unknown Study/CT 0.625mm'
reader = vtk.vtkMetaImageReader()  # 读取、显示 .mha 图像
reader.SetFileName('SaveRaw/sample.mhd')

ren = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

volumeMapper = vtk.vtkGPUVolumeRayCastMapper()
volumeMapper.SetInputConnection(reader.GetOutputPort())
volumeMapper.SetBlendModeToComposite()

volumeColor = vtk.vtkColorTransferFunction()    #颜色
volumeColor.AddRGBPoint(0,    0.0, 0.0, 0.0)
volumeColor.AddRGBPoint(500,  1.0, 0.5, 0.3)
volumeColor.AddRGBPoint(1000, 1.0, 0.5, 0.3)
volumeColor.AddRGBPoint(1150, 1.0, 1.0, 0.9)

volumeScalarOpacity = vtk.vtkPiecewiseFunction()    #不透明度，越大越不透明
volumeScalarOpacity.AddPoint(0,    0.00)
volumeScalarOpacity.AddPoint(500,  0.15)
volumeScalarOpacity.AddPoint(1000, 0.15)
volumeScalarOpacity.AddPoint(1150, 0.85)

volumeGradientOpacity = vtk.vtkPiecewiseFunction()
volumeGradientOpacity.AddPoint(0,   0.0)
volumeGradientOpacity.AddPoint(90,  0.5)
volumeGradientOpacity.AddPoint(100, 1.0)

volumeProperty = vtk.vtkVolumeProperty()
volumeProperty.SetColor(volumeColor)
volumeProperty.SetScalarOpacity(volumeScalarOpacity)
# volumeProperty.SetGradientOpacity(volumeGradientOpacity)
volumeProperty.SetInterpolationTypeToLinear()
volumeProperty.ShadeOn()
volumeProperty.SetAmbient(0.9)  #设置环境
volumeProperty.SetDiffuse(0.9)  #设置漫反射
volumeProperty.SetSpecular(0.9) #设置镜面反射

volume = vtk.vtkVolume()
volume.SetMapper(volumeMapper)
volume.SetProperty(volumeProperty)

ren.AddViewProp(volume)

camera = ren.GetActiveCamera()
c = volume.GetCenter()
camera.SetFocalPoint(c[0], c[1], c[2])
# camera.SetPosition(c[0] + 400, c[1], c[2])
camera.SetPosition(c[0] + 500, c[1], c[2])
camera.SetViewUp(0, 0, -1)

renWin.SetSize(640, 480)
# renWin.SetSize(940, 650)

renWin.Render()
iren.Initialize()

iren.Start()

数据文件下载地址： python CT切片图像三维重建（数据和代码）.zip