Solución:
Aquí hay una versión funcional del código que necesitaba escribir. Escribirá una imagen DICOM en escala de grises de 16 bits a partir de un 2D dado array de píxeles. Según el estándar DICOM, los UID deben ser únicos para cada imagen y serie, por lo que este código no se preocupa, porque no sé qué hacen realmente los UID. Si alguien más lo hace, estaré feliz de agregarlo.
import dicom, dicom.UID
from dicom.dataset import Dataset, FileDataset
import numpy as np
import datetime, time
def write_dicom(pixel_array,filename):
"""
INPUTS:
pixel_array: 2D numpy ndarray. If pixel_array is larger than 2D, errors.
filename: string name for the output file.
"""
## This code block was taken from the output of a MATLAB secondary
## capture. I do not know what the long dotted UIDs mean, but
## this code works.
file_meta = Dataset()
file_meta.MediaStorageSOPClassUID = 'Secondary Capture Image Storage'
file_meta.MediaStorageSOPInstanceUID = '1.3.6.1.4.1.9590.100.1.1.111165684411017669021768385720736873780'
file_meta.ImplementationClassUID = '1.3.6.1.4.1.9590.100.1.0.100.4.0'
ds = FileDataset(filename, ,file_meta = file_meta,preamble=" "*128)
ds.Modality = 'WSD'
ds.ContentDate = str(datetime.date.today()).replace('-','')
ds.ContentTime = str(time.time()) #milliseconds since the epoch
ds.StudyInstanceUID = '1.3.6.1.4.1.9590.100.1.1.124313977412360175234271287472804872093'
ds.SeriesInstanceUID = '1.3.6.1.4.1.9590.100.1.1.369231118011061003403421859172643143649'
ds.SOPInstanceUID = '1.3.6.1.4.1.9590.100.1.1.111165684411017669021768385720736873780'
ds.SOPClassUID = 'Secondary Capture Image Storage'
ds.SecondaryCaptureDeviceManufctur = 'Python 2.7.3'
## These are the necessary imaging components of the FileDataset object.
ds.SamplesPerPixel = 1
ds.PhotometricInterpretation = "MONOCHROME2"
ds.PixelRepresentation = 0
ds.HighBit = 15
ds.BitsStored = 16
ds.BitsAllocated = 16
ds.SmallestImagePixelValue = '\x00\x00'
ds.LargestImagePixelValue = '\xff\xff'
ds.Columns = pixel_array.shape[0]
ds.Rows = pixel_array.shape[1]
if pixel_array.dtype != np.uint16:
pixel_array = pixel_array.astype(np.uint16)
ds.PixelData = pixel_array.tostring()
ds.save_as(filename)
return
if __name__ == "__main__":
# pixel_array = np.arange(256*256).reshape(256,256)
# pixel_array = np.tile(np.arange(256).reshape(16,16),(16,16))
x = np.arange(16).reshape(16,1)
pixel_array = (x + x.T) * 32
pixel_array = np.tile(pixel_array,(16,16))
write_dicom(pixel_array,'pretty.dcm')
El ejemplo anterior funciona, pero hace que muchas herramientas se quejen de los DICOM y ni siquiera se puedan leer usando itk / SimpleITK como una pila. La mejor manera que he encontrado para hacer DICOM a partir de numpy es utilizando las herramientas SimpleITK y generando los DICOM rebanada por rebanada. Un ejemplo básico (https://github.com/zivy/SimpleITK/blob/8e94451e4c0e90bcc6a1ffdd7bc3d56c81f58d80/Examples/DicomSeriesReadModifyWrite/DicomSeriesReadModifySeriesWrite.py) muestra cómo cargar, y luego volver a guardar archivos, realizar una transformación fácilmente. modificado utilizando el
import SimpleITK as sitk
filtered_image = sitk.GetImageFromArray(my_numpy_array)
El número de etiquetas en última instancia en la imagen de salida es bastante grande y, por lo tanto, crearlas manualmente es tedioso. Además, SimpleITK admite imágenes de 8, 16 y 32 bits, así como RGB, por lo que es mucho más fácil que hacerlas en pydicom.
(0008, 0008) Image Type CS: ['DERIVED', 'SECONDARY']
(0008, 0016) SOP Class UID UI: Secondary Capture Image Storage
(0008, 0018) SOP Instance UID UI: 1.2.826.0.1.3680043.2.1125.1.35596048796922805578234000521866725
(0008, 0020) Study Date DA: '20170803'
(0008, 0021) Series Date DA: '20170803'
(0008, 0023) Content Date DA: 0
(0008, 0030) Study Time TM: '080429.171808'
(0008, 0031) Series Time TM: '080429'
(0008, 0033) Content Time TM: 0
(0008, 0050) Accession Number SH: ''
(0008, 0060) Modality CS: 'OT'
(0008, 0064) Conversion Type CS: 'WSD'
(0008, 0090) Referring Physician's Name PN: ''
(0010, 0010) Patient's Name PN: ''
(0010, 0020) Patient ID LO: ''
(0010, 0030) Patient's Birth Date DA: ''
(0010, 0040) Patient's Sex CS: ''
(0018, 2010) Nominal Scanned Pixel Spacing DS: ['1', '3']
(0020, 000d) Study Instance UID UI: 1.2.826.0.1.3680043.2.1125.1.33389357207068897066210100430826006
(0020, 000e) Series Instance UID UI: 1.2.826.0.1.3680043.2.1125.1.51488923827429438625199681257282809
(0020, 0010) Study ID SH: ''
(0020, 0011) Series Number IS: ''
(0020, 0013) Instance Number IS: ''
(0020, 0020) Patient Orientation CS: ''
(0020, 0052) Frame of Reference UID UI: 1.2.826.0.1.3680043.2.1125.1.35696880630664441938326682384062489
(0028, 0002) Samples per Pixel US: 1
(0028, 0004) Photometric Interpretation CS: 'MONOCHROME2'
(0028, 0010) Rows US: 40
(0028, 0011) Columns US: 50
(0028, 0100) Bits Allocated US: 32
(0028, 0101) Bits Stored US: 32
(0028, 0102) High Bit US: 31
(0028, 0103) Pixel Representation US: 1
(0028, 1052) Rescale Intercept DS: "0"
(0028, 1053) Rescale Slope DS: "1"
(0028, 1054) Rescale Type LO: 'US'
(7fe0, 0010) Pixel Data OW: Array of 8000 bytes
Actualización 2020 🙂
Ninguna de estas respuestas funcionó para mí. Esto es lo que terminé con para guardar un segmento MR monocromático válido de 16 bpp que se muestra correctamente al menos en Slicer, Radiant y MicroDicom:
import pydicom
from pydicom.dataset import Dataset, FileDataset
from pydicom.uid import ExplicitVRLittleEndian
import pydicom._storage_sopclass_uids
image2d = image2d.astype(np.uint16)
print("Setting file meta information...")
# Populate required values for file meta information
meta = pydicom.Dataset()
meta.MediaStorageSOPClassUID = pydicom._storage_sopclass_uids.MRImageStorage
meta.MediaStorageSOPInstanceUID = pydicom.uid.generate_uid()
meta.TransferSyntaxUID = pydicom.uid.ExplicitVRLittleEndian
ds = Dataset()
ds.file_meta = meta
ds.is_little_endian = True
ds.is_implicit_VR = False
ds.SOPClassUID = pydicom._storage_sopclass_uids.MRImageStorage
ds.PatientName = "Test^Firstname"
ds.PatientID = "123456"
ds.Modality = "MR"
ds.SeriesInstanceUID = pydicom.uid.generate_uid()
ds.StudyInstanceUID = pydicom.uid.generate_uid()
ds.FrameOfReferenceUID = pydicom.uid.generate_uid()
ds.BitsStored = 16
ds.BitsAllocated = 16
ds.SamplesPerPixel = 1
ds.HighBit = 15
ds.ImagesInAcquisition = "1"
ds.Rows = image2d.shape[0]
ds.Columns = image2d.shape[1]
ds.InstanceNumber = 1
ds.ImagePositionPatient = r"0 1"
ds.ImageOrientationPatient = r"1 -1 "
ds.ImageType = r"ORIGINALPRIMARYAXIAL"
ds.RescaleIntercept = "0"
ds.RescaleSlope = "1"
ds.PixelSpacing = r"11"
ds.PhotometricInterpretation = "MONOCHROME2"
ds.PixelRepresentation = 1
pydicom.dataset.validate_file_meta(ds.file_meta, enforce_standard=True)
print("Setting pixel data...")
ds.PixelData = image2d.tobytes()
ds.save_as(r"out.dcm")
Tenga en cuenta lo siguiente:
- Pasar por el constructor FileDataset como sugieren los documentos de PyDicom no pudo crear un encabezado válido para mí
- validate_file_meta creará algunos elementos faltantes en el encabezado para usted (versión)
- Necesita especificar endianness y VR explícita / implícita dos veces: /
- Este método le permitirá crear un volumen válido y siempre que actualice ImagePositionPatient y InstanceNumber para cada segmento en consecuencia
- Asegúrate de que tu numpy array se convierte en formato de datos que tiene la misma cantidad de bits que su BitsStored
Te mostramos las comentarios y valoraciones de los lectores
Te invitamos a corroborar nuestra publicación añadiendo un comentario o puntuándolo te estamos eternamente agradecidos.