irAEMM/dicomUtils/loadDicom.py

import pydicom
import numpy
import sys
import os
import nibabel
import json
import datetime
import re

# load the DICOM files
def load(dir):
    files = []
    print('Loading: {}'.format(dir))
    for fname in os.listdir(dir):
        print("loading: {}".format(fname))
        files.append(pydicom.dcmread(os.path.join(dir,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)
    return slices;

    # pixel aspects, assuming all slices are the same

    for s in slices:
        print("Pixel spacing: {}".format(s.PixelSpacing))

    return;

def convertToNIfTI(slices):

    # create 3D array
    img_shape = list(slices[0].pixel_array.shape)
    img_shape.append(len(slices))
    img3d = numpy.zeros(img_shape)

    # get position and orientation parameters

    x0=numpy.array([float(f) for f in slices[0].ImagePositionPatient])
    x1=numpy.array([float(f) for f in slices[-1].ImagePositionPatient])
    n=(x1-x0)/(len(slices)-1.);

    f=numpy.array([float(f) for f in slices[0].ImageOrientationPatient])

    s=numpy.array([float(f) for f in slices[0].PixelSpacing])
    
    #create affine 
    a=numpy.zeros((4,4))
    f1=numpy.zeros(4)
    f1[0:3]=f[0:3]
    #swap to account for row/column to (c,r) position mismatch
    a[:,1]=f1*s[1]
    f2=numpy.zeros(4)
    f2[0:3]=f[3:6]
    a[:,0]=f2*s[0]
    nn=numpy.zeros(4)
    nn[0:3]=n[0:3]
    a[:,2]=nn

    xn=numpy.ones(4)
    xn[0:3]=x0[0:3]
    a[:,3]=xn

     

    # fill 3D array with the images from the files
    for i, s in enumerate(slices):
        img2d = s.pixel_array
        img3d[:, :, i] = img2d
    #orientation and pixel spacing
    img = nibabel.Nifti1Image(img3d, a)
    img.header.set_slope_inter(float(slices[0].RescaleSlope),float(slices[0].RescaleIntercept))

    return img

def writeAnonymousSeries(slices,path,patientID, studyUID, fields):
    uid=uuid()
    seriesUID=uid.generateSeriesUUID('volume')
    i=0;
    for s in slices:
        outFile="DCM{:04d}.dcm".format(i)
        outFile=os.path.join(path,outFile)
        writeAnonymousDicomFile(s,patientID,studyUID,seriesUID,outFile,fields)
        i=i+1




def writeAnonymousDicomFile(slice, patientID, studyUID, seriesUID, outFile, fields):
    uid=uuid()
    instanceUID=uid.generateSOPInstanceUUID('volume')


    file_meta = pydicom.Dataset()
    file_meta.MediaStorageSOPClassUID = slice.file_meta.MediaStorageSOPClassUID
    file_meta.MediaStorageSOPInstanceUID = instanceUID
    file_meta.ImplementationClassUID = slice.file_meta.ImplementationClassUID
    file_meta.TransferSyntaxUID=pydicom.uid.ImplicitVRLittleEndian

    print("Setting dataset values...")
    # Create the FileDataset instance (initially no data elements, but file_meta
    # supplied)
    ds = pydicom.FileDataset(outFile, {},
                         file_meta=file_meta, preamble=b"\0" * 128)

    # Add the data elements -- not trying to set all required here. Check DICOM
    # standard
    ds.PatientName = patientID
    ds.PatientID = patientID
    ds.StudyInstanceUID=studyUID
    ds.SeriesInstanceUID=seriesUID
    ds.SOPInstanceUID=instanceUID

    # Set the transfer syntax
    ds.is_little_endian = True
    ds.is_implicit_VR = True

    # Set creation date/time
    #dt = datetime.datetime.now()
    #ds.ContentDate = dt.strftime('%Y%m%d')
    #timeStr = dt.strftime('%H%M%S.%f')  # long format with micro seconds
    ds.ContentTime = slice.ContentTime
    #add variables from Daniel's list (fields)
    for f in fields:
        try:
            ds[f]=slice[f];
        except:
            print("Field {} missing".format(f))

    ds.PixelData=slice.PixelData
    print("Writing test file", outFile)
    ds.save_as(outFile)
    print("File saved.")

class uuid:

    def __init__(self):
        fhome=os.path.expanduser('~')
        fconfig=os.path.join(fhome,'.uuid','config.json')
        self.basePath=os.path.join(fhome,'.uuid')
        with open(fconfig,'r') as f:
            self.config=json.load(f)
        self.labelUUID={'series':'1','study':'2','instance':'3','frameOfReference':4}
        self.dataUUID={'volume':'1','segmentation':'2','transformation':'3'}

    def getTimeCode(x):
        hour=x.strftime("%H")
        hour=re.sub('^0','',hour)
        return hour+x.strftime("%M%S")
 
    def generateStudyUUID(self,type):

        x=datetime.datetime.now()
        date=x.strftime("%Y%m%d")
        studyFile=os.path.join(self.basePath,'studyCount'+date+'.txt')

        try:
            f=open(studyFile,"r")
            id=int(f.readline())
            id=id+1
            f.close()
        except:
            id=0

        studyId="{}.{}.{}.{}.{}".format(self.config["baseUUID"],self.labelUUID['study'],
                        self.dataUUID[type],date,id)

        f=open(studyFile,"w")
        f.write("{}".format(id))
        f.close()
        return studyId


    def generateSOPInstanceUUID(self,type):

        baseUUID=self.config["baseUUID"]
        x=datetime.datetime.now()
        date=x.strftime("%Y%m%d")
        instanceFile=os.path.join(self.basePath,'instanceCount'+date+'.txt')

        try:
            f=open(instanceFile,"r")
            id=int(f.readline())
            id=id+1
            f.close()
        except:
            id=0

        instanceId="{}.{}.{}.{}.{}".format(baseUUID,self.labelUUID['instance'],
                        self.dataUUID[type],date,id)

        f=open(instanceFile,"w")
        f.write("{}".format(id))
        f.close()
        return instanceId


    def generateFrameOfReferenceUUID(self,type):

        baseUUID=self.config["baseUUID"]
        basePath=self.basePath
        x=datetime.datetime.now()
        date=x.strftime("%Y%m%d")
        forFile=os.path.join(basePath,'frameCount'+date+'.txt')

        try:
            f=open(studyFile,"r")
            id=int(f.readline())
            id=id+1
            f.close()
        except:
            id=0

        forId="{}.{}.{}.{}.{}".format(baseUUID,self.labelUUID['frameOfReference'],
                        self.dataUUID[type],date,id)

        f=open(forFile,"w")
        f.write("{}".format(id))
        f.close()
        return forId
  
    def generateSeriesUUID(self,type):

        baseUUID=self.config["baseUUID"]
        x=datetime.datetime.now()
        seriesInstanceUid=baseUUID+'.'+self.labelUUID['series']+'.'
        seriesInstanceUid+=self.dataUUID[type]+'.'+x.strftime("%Y%m%d")+'.'+uuid.getTimeCode(x)
        return seriesInstanceUid