reslice.py

"""
==========================================================
Load CT slices and plot axial, sagittal and coronal images
==========================================================

This example illustrates loading multiple files, sorting them by slice
location, building a 3D image and reslicing it in different planes.

.. usage:

   reslice.py <glob>
   where <glob> refers to a set of DICOM image files.

   Example: python reslice.py "*.dcm". The quotes are needed to protect
   the glob from your system and leave it for the script.

.. note:

   Uses numpy and matplotlib.

   Tested using series 2 from here
   http://www.pcir.org/researchers/54879843_20060101.html
"""

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

#获取根目录文件
def get_root_file(root_file):
    list_file = os.listdir(root_file)
    list_file.sort(key=lambda x: str(x.split('.')[0]))
    return list_file

root_dicom="E:\\bigorgs\\Datasets\\1\\"
list_dicom = get_root_file(root_dicom)
num_dicom = len(list_dicom)
files = []
for i in range(num_dicom):
    file_root = root_dicom + list_dicom[i]
    dcm = pydicom.dcmread(file_root,force=True)
    dcm.file_meta.TransferSyntaxUID = pydicom.uid.ImplicitVRLittleEndian
    # print(dcm)
    print("loading: {}".format(file_root))
    # files.append(pydicom.dcmread(file_root,force=True))
    files.append(dcm)



# load the DICOM files
# files = []
# print('glob: {}'.format(sys.argv[1]))
# for fname in glob.glob(sys.argv[1], recursive=False):
#     print("loading: {}".format(fname))
#     files.append(pydicom.dcmread(fname))

print("file count: {}".format(len(files)))

# skip files with no SliceLocation (eg scout views)
slices = []
skipcount = 0
for f in files:
    if hasattr(f, 'SliceLocation'):
        slices.append(f)
    else:
        skipcount = skipcount + 1

print("skipped, no SliceLocation: {}".format(skipcount))

# ensure they are in the correct order
slices = sorted(slices, key=lambda s: s.SliceLocation)

# pixel aspects, assuming all slices are the same
ps = slices[0].PixelSpacing
ss = slices[0].SliceThickness
ax_aspect = ps[1]/ps[0]
sag_aspect = ps[1]/ss
cor_aspect = ss/ps[0]
# print(slices[0].pixel_array.shape)
# create 3D array
img_shape = list(slices[0].pixel_array.shape)
img_shape.append(len(slices))
img3d = np.zeros(img_shape)

# fill 3D array with the images from the files
for i, s in enumerate(slices):
    img2d = s.pixel_array
    img3d[:, :, i] = img2d

# plot 3 orthogonal slices
a1 = plt.subplot(1, 3, 1)
ax=img3d[:, :, img_shape[2]//2]
print("ax:{}".format(ax))
plt.imshow(img3d[:, :, img_shape[2]//2])
a1.set_aspect(ax_aspect)
# a1.set_aspect('auto')

a3 = plt.subplot(1, 3, 2)
cor=img3d[img_shape[0]//2, :, :]
print("cor:{}".format(cor))
plt.imshow(img3d[img_shape[0]//2, :, :])
a3.set_aspect(cor_aspect)
# a3.set_aspect('auto')

a2 = plt.subplot(1, 3, 3)
sag=img3d[:, img_shape[1]//2, :].T
print("sag:{}".format(sag))
plt.imshow(img3d[:, img_shape[1]//2, :].T)
a2.set_aspect(sag_aspect)             #坐标轴比例
# a2.set_aspect('auto')



plt.show()