Извлечение сагиттальных и коронарных срезов с осевого обзора с использованием pydicom - PullRequest
Новая
       20

Извлечение сагиттальных и коронарных срезов с осевого обзора с использованием pydicom

0 голосов
/ 18 декабря 2018

Я пытаюсь прочитать серию файлов .dcm, которые по умолчанию показывают осевой вид.Ниже приведен код: 

import os
import numpy as np
import pydicom as dicom
from matplotlib import pyplot as plt

root_dir = 'mydcomDir'


def sortDcm():
        print('Given Path to the .dcm directory is: {}'.format(root_dir))
        slices = [dicom.read_file(root_dir + '/' + s) for s in os.listdir(root_dir)]
        slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
        pos1 = slices[int(len(slices)/2)].ImagePositionPatient[2]
        pos2 = slices[(int(len(slices)/2)) + 1].ImagePositionPatient[2]
        diff = pos2 - pos1
#        if diff > 0:
#            slices = np.flipud(slices)
        try:
            slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
        except:
            slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)

        for s in slices:
            s.SliceThickness = slice_thickness
#        print("from sorted dicom",len(slices))         
        return slices 


dcms = sortDcm()
ref_dicom = dcms[0]

d_array = np.zeros((ref_dicom.Columns,ref_dicom.Rows, len(dcms)), dtype=ref_dicom.pixel_array.dtype)

for dcm in dcms:
    d_array[:, :, dcms.index(dcm)] = dcm.pixel_array

#    fig = plt.figure(figsize=(12,12))
#    plt.subplot(1, 3, 1)
#    plt.title("Coronal")
#    plt.imshow(np.flipud(d_array[idx , :, :].T))
#    plt.subplot(1, 3, 2)
#    plt.title("Sagital")
#    plt.imshow(np.flipud(d_array[:, idy, :].T))
#    plt.subplot(1, 3, 3)
    plt.title("axial")
    plt.imshow(d_array[:, :, dcms.index(dcm)])
    plt.pause(0.001)

Как видно из кода, я не смог выяснить соответствующие idx и id для конкретного файла dcm.Поэтому мой вопрос заключается в том, как получить сагиттальные и корональные разрезы и нанести их на график, учитывая осевые разрезы?

Заранее спасибо.

Редактировать: @ColonelFazackerley ответил идеально.Я добавляю строку ниже, чтобы показать, как я ее использовал. 

# fill 3D array with the images from the files
for i, s in enumerate(slices):
    img2d = s.pixel_array
    img3d[:,:,i] = img2d
#then to view sagittal and coronal slices for each of the axial slice
for i, s in enumerate(slices):
    img2d = s.pixel_array
    img3d[:,:,i] = img2d
    corId = corId-1
    sagId = sagId-1
#    plot 3 orthogonal slices
    a1 = plt.subplot(1,3,1)
    plt.title('Axial')
    plt.imshow(img3d[:,:,i],'gray')
    a1.set_aspect(ax_aspect)

    a2 = plt.subplot(1,3,2)
    plt.title('Sagittal')
    plt.imshow(np.flipud(img3d[:,sagId,:].T),'gray')
    a2.set_aspect(sag_aspect)

    a3 = plt.subplot(1,3,3)
    plt.imshow(np.flipud(img3d[corId,:,:].T),'gray')
    a3.set_aspect(cor_aspect)
    plt.title('Coronal')
    plt.show()
    plt.pause(0.001)

1 Ответ

0 голосов
/ 26 декабря 2018 
"""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."""

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

# 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.read_file(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]

# 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(2,2,1)
plt.imshow(img3d[:,:,img_shape[2]//2])
a1.set_aspect(ax_aspect)

a2 = plt.subplot(2,2,2)
plt.imshow(img3d[:,img_shape[1]//2,:])
a2.set_aspect(sag_aspect)

a3 = plt.subplot(2,2,3)
plt.imshow(img3d[img_shape[0]//2,:,:].T)
a3.set_aspect(cor_aspect)

plt.show()

проверено на серии 2 из этого примера. Данные 3D CT http://www.pcir.org/researchers/54879843_20060101.html

редактировать примечание: принято в качестве примера в проект pydicom https://github.com/pydicom/pydicom/blob/master/examples/image_processing/reslice.py

...