Python dicom.read_file() Examples

def read_images_info(path):
    for subdir, dirs, files in os.walk(path):
        dcms = glob.glob(os.path.join(subdir, '*.dcm'))
        if len(dcms) > 1:
            slices = [dicom.read_file(dcm) for dcm in dcms]
            slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
            images = np.stack([s.pixel_array for s in slices], axis=0).astype(np.float32)
            images = images + slices[0].RescaleIntercept
    
    orig_shape = images.shape
    
    inplane_scale = slices[0].PixelSpacing[0] / PIXEL_SPACING
    inplane_size = int(np.rint(inplane_scale * slices[0].Rows / 2) * 2)
    return orig_shape, inplane_size 

def wrapper_from_file(file_like, *args, **kwargs):
    ''' Create DICOM wrapper from `file_like` object

    Parameters
    ----------
    file_like : object
       filename string or file-like object, pointing to a valid DICOM
       file readable by ``pydicom``
    *args : positional
    **kwargs : keyword
        args to ``dicom.read_file`` command.  ``force=True`` might be a
        likely keyword argument.

    Returns
    -------
    dcm_w : ``dicomwrappers.Wrapper`` or subclass
       DICOM wrapper corresponding to DICOM data type
    '''
    import dicom
    fobj = allopen(file_like)
    dcm_data = dicom.read_file(fobj, *args, **kwargs)
    return wrapper_from_data(dcm_data) 

def load_scan(path):
    slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]
    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
    return slices 

def load_scan(path):
    slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]
    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
    return slices 

def load_ch4(self):
        ch4dirs = [d for d in os.listdir(self.directory) if '4ch' in d]
        max_val = 0
        max_dir = None
        self.ch4_images = None
        for d in ch4dirs:
            fn = [f for f in os.listdir(os.path.join(self.directory, d)) if 'dcm' in f][0]
            series = int(f.split('-')[1])
            if series > max_val:
                max_val = series
                max_dir = d
        if max_dir is not None:
            ch4_fns = [f for f in os.listdir(os.path.join(self.directory, max_dir))
                    if 'dcm' in f]
            ch4_fns = sorted(ch4_fns, key=lambda f: int(f.split('.')[0].split('-')[2]))
            ch4_ims = np.empty((len(ch4_fns), self.img_size, self.img_size))
            for i,fn in enumerate(ch4_fns):
                d = dicom.read_file(os.path.join(self.directory, max_dir, fn))
                ch4_ims[i] = crop_resize(d.pixel_array, self.img_size)
                if i == 0:
                    short_edge = min(d.pixel_array.shape)
                    self.ch4_line_mult = float(d.PixelSpacing[0])*short_edge*1./self.img_size
            self.ch4_images = ch4_ims 

def load_contour(contour, img_path):
    filename = "IM-%s-%04d.dcm" % (SAX_SERIES[contour.case], contour.img_no)
    full_path = os.path.join(img_path, contour.case, filename)
    f = dicom.read_file(full_path)
    ctrs = np.loadtxt(contour.ctr_path, delimiter=" ").astype(np.int)
    label = np.zeros(f.pixel_array.shape, dtype=np.uint8)
    cv2.fillPoly(label, [ctrs], 255)
    img,lab = getAlignImg(f,label);
    lx,ly = img.shape;
    assert(lx==ly);
    xm,ym = np.where(lab>127);
    if xm.size<30:
        xm,ym = lx//2,ly//2;
    xm = np.mean(xm);
    ym = np.mean(ym);
    delta = int(lx*0.62)//2;#cut middle 160x160 from 256x256 for sunny brook data
    assert(delta<xm and delta<ym);
    xm,ym,delta = int(xm),int(ym),int(delta);
    img = img[xm-delta:xm+delta,ym-delta:ym+delta];
    lab = lab[xm-delta:xm+delta,ym-delta:ym+delta];
    return cv2.resize(img, (SZ,SZ)), cv2.resize(lab, (SZ,SZ)) 

def load_ch4(self):
        ch4dirs = [d for d in os.listdir(self.directory) if '4ch' in d]
        max_val = 0
        max_dir = None
        self.ch4_images = None
        for d in ch4dirs:
            fn = [f for f in os.listdir(os.path.join(self.directory, d)) if 'dcm' in f][0]
            series = int(f.split('-')[1])
            if series > max_val:
                max_val = series
                max_dir = d
        if max_dir is not None:
            ch4_fns = [f for f in os.listdir(os.path.join(self.directory, max_dir))
                    if 'dcm' in f]
            ch4_fns = sorted(ch4_fns, key=lambda f: int(f.split('.')[0].split('-')[2]))
            ch4_ims = np.empty((len(ch4_fns), self.img_size, self.img_size))
            for i,fn in enumerate(ch4_fns):
                d = dicom.read_file(os.path.join(self.directory, max_dir, fn))
                ch4_ims[i] = crop_resize(d.pixel_array, self.img_size)
                if i == 0:
                    short_edge = min(d.pixel_array.shape)
                    self.ch4_line_mult = float(d.PixelSpacing[0])*short_edge*1./self.img_size
            self.ch4_images = ch4_ims 

def is_file_a_dicom(file):
    """
    Check whether a given file is of type DICOM

    :param file: path to the file to identify
     :type file: str

    :return: True if the file is DICOM, False otherwise
     :rtype: bool

    """

    try:
        dicom.read_file(file)
    except InvalidDicomError:
        return False
    return True 

def read_ct_scan(folder_name):
    # Read the slices from the dicom file
    slices = [dicom.read_file(folder_name + '/' + s) for s in os.listdir(folder_name)]

    # Sort the dicom slices in their respective order
    slices.sort(key=lambda x: int(x.InstanceNumber))

    # Get the pixel values for all the slices
    slices = np.stack([s.pixel_array for s in slices])
    slices[slices == -2000] = 0
    return slices 

def get_data(lst,preproc):
   data = []
   result = []
   for path in lst:
       f = dicom.read_file(path)
       img = preproc(f.pixel_array.astype(float) / np.max(f.pixel_array))
       dst_path = path.rsplit(".", 1)[0] + ".64x64.jpg"
       scipy.misc.imsave(dst_path, img)
       result.append(dst_path)
       data.append(img)
   data = np.array(data, dtype=np.uint8)
   data = data.reshape(data.size)
   data = np.array(data, dtype=np.str_)
   data = data.reshape(data.size)
   return [data,result] 

def animate(i):
        d = dicom.read_file(file_list[i])
        img = d.pixel_array.astype('int')
        im.set_data(img)
        return im, 

def read_dicom(filename):
    d = dicom.read_file(filename)
    metadata = {}
    for attr in dir(d):
        if attr[0].isupper() and attr != 'PixelData':
            try:
                metadata[attr] = getattr(d, attr)
            except AttributeError:
                pass
    return np.array(d.pixel_array), metadata 

def read_dicom(filename):
    d = dicom.read_file(filename)
    metadata = {}
    for attr in dir(d):
        if attr[0].isupper() and attr != 'PixelData':
            try:
                metadata[attr] = getattr(d, attr)
            except AttributeError:
                pass
    return np.array(d.pixel_array), metadata 

def load_contour(contour, img_path):
    filename = "IM-%s-%04d.dcm" % (SAX_SERIES[contour.case], contour.img_no)
    full_path = os.path.join(img_path, contour.case, filename)
    f = dicom.read_file(full_path)
    img = f.pixel_array.astype(np.int)
    ctrs = np.loadtxt(contour.ctr_path, delimiter=" ").astype(np.int32)
    label = np.zeros_like(img, dtype="uint8")
    cv2.fillPoly(label, [ctrs], 1)
    return img, label 

def read_dicom(filename):
    d = dicom.read_file(filename)
    metadata = {}
    for attr in dir(d):
        if attr[0].isupper() and attr != 'PixelData':
            try:
                metadata[attr] = getattr(d, attr)
            except AttributeError:
                pass
    return np.array(d.pixel_array), metadata 

def get_data(lst,preproc):
   data = []
   result = []
   for path in lst:
       f = dicom.read_file(path)
       img = preproc(f.pixel_array.astype(float) / np.max(f.pixel_array))
       dst_path = path.rsplit(".", 1)[0] + ".64x64.jpg"
       scipy.misc.imsave(dst_path, img)
       result.append(dst_path)
       data.append(img)
   data = np.array(data, dtype=np.uint8)
   data = data.reshape(data.size)
   data = np.array(data, dtype=np.str_)
   data = data.reshape(data.size)
   return [data,result] 

def dicom_datasets_from_zip(zip_file):
    datasets = []
    for entry in zip_file.namelist():
        if entry.endswith('/'):
            continue  # skip directories

        entry_pseudo_file = zip_file.open(entry)

        # the pseudo file does not support `seek`, which is required by
        # dicom's lazy loading mechanism; use temporary files to get around this;
        # relies on the temporary files not being removed until the temp
        # file is garbage collected, which should be the case because the
        # dicom datasets should retain a reference to the temp file
        temp_file = tempfile.TemporaryFile()
        temp_file.write(entry_pseudo_file.read())
        temp_file.flush()
        temp_file.seek(0)

        try:
            dataset = dicom.read_file(temp_file)
            datasets.append(dataset)
        except dicom.errors.InvalidDicomError as e:
            msg = 'Skipping invalid DICOM file "{}": {}'
            logger.info(msg.format(entry, e))

    if len(datasets) == 0:
        raise DicomImportException('Zipfile does not contain any valid DICOM files')

    return datasets 

def read_dicom_series(directory, filepattern="image_*"):
    """ Reads a DICOM Series files in the given directory.
    Only filesnames matching filepattern will be considered"""

    if not os.path.exists(directory) or not os.path.isdir(directory):
        raise ValueError("Given directory does not exist or is a file : " + str(directory))
    print
    '\tRead Dicom', directory
    lstFilesDCM = natsort.natsorted(glob.glob(os.path.join(directory, filepattern)))
    print
    '\tLength dicom series', len(lstFilesDCM)
    # Get ref file
    RefDs = dicom.read_file(lstFilesDCM[0])
    # Load dimensions based on the number of rows, columns, and slices (along the Z axis)
    ConstPixelDims = (int(RefDs.Rows), int(RefDs.Columns), len(lstFilesDCM))
    # The array is sized based on 'ConstPixelDims'
    ArrayDicom = np.zeros(ConstPixelDims, dtype=RefDs.pixel_array.dtype)

    # loop through all the DICOM files
    for filenameDCM in lstFilesDCM:
        # read the file
        ds = dicom.read_file(filenameDCM)
        # store the raw image data
        ArrayDicom[:, :, lstFilesDCM.index(filenameDCM)] = ds.pixel_array

    return ArrayDicom 

def load_scan(path):
    # 子文件夹
    slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]
    slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))
    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

    return slices 

def load_contour(contour, img_path):
    filename = 'IM-%s-%04d.dcm' % (sax_series_dict[contour.case], contour.img_no)
    full_path = os.path.join(img_path, contour.case, filename)
    f = dicom.read_file(full_path)
    img = f.pixel_array.astype(np.uint8)
    ctrs = np.loadtxt(contour.ctr_path, delimiter=' ').astype(np.int)
    label = np.zeros_like(img, dtype='uint8')
    cv2.fillPoly(label, [ctrs], 1)
    return cv2.resize(img, (img_size,img_size)), cv2.resize(label, (img_size,img_size)) 

def read_dicom_with_pydicom(dicom_file, dicom_fields):
    """
    Read DICOM file using PyDICOM python library.

    :param dicom_file: DICOM file to read
     :type dicom_file: str
    :param dicom_fields: Dictionary containing DICOM fields and values
     :type dicom_fields: dict

    :return: updated dictionary of DICOM fields and values
     :rtype : dict

    """

    # Read DICOM file
    dicom_dataset = dicom.read_file(dicom_file)

    # Grep information from DICOM header and store them
    # into dicom_fields dictionary under flag Value
    # Dictionnary of DICOM values to be returned
    for name in dicom_fields:
        try:
            description = dicom_fields[name]['Description']
            value = dicom_dataset.data_element(description).value
            dicom_fields[name]['Value'] = value
        except:
            continue

    return dicom_fields 

def pydicom_zapping(dicom_file, dicom_fields):
    """
    Actual zapping method for PyDICOM

    :param dicom_file: DICOM to de-identify
     :type dicom_file: str
    :param dicom_fields: Dictionary with DICOM fields & values to use
     :type dicom_fields: dict

    :return: None

    """

    dicom_dataset = dicom.read_file(dicom_file)

    for name in dicom_fields:
        new_val = ""
        if 'Value' in dicom_fields[name]:
            new_val = dicom_fields[name]['Value'].strip()

        if dicom_fields[name]['Editable'] is True:
            try:
                dicom_dataset.data_element(
                    dicom_fields[name]['Description']).value = new_val
            except:
                continue
        else:
            try:
                dicom_dataset.data_element(
                    dicom_fields[name]['Description']).value = ''
            except:
                continue
    dicom_dataset.save_as(dicom_file) 

def get_data(lst,preproc):
   data = []
   result = []
   for path in lst:
       f = dicom.read_file(path)
       img = preproc(f.pixel_array.astype(float) / np.max(f.pixel_array))
       dst_path = path.rsplit(".", 1)[0] + ".64x64.jpg"
       scipy.misc.imsave(dst_path, img)
       result.append(dst_path)
       data.append(img)
   data = np.array(data, dtype=np.uint8)
   data = data.reshape(data.size)
   data = np.array(data, dtype=np.str_)
   data = data.reshape(data.size)
   return [data,result] 

def get_3d_data(path):
    slices = [dicom.read_file(os.path.join(path, s)) for s in os.listdir(path)]
    slices.sort(key=lambda x: int(x.InstanceNumber))
    return np.stack([s.pixel_array for s in slices]) 

def read_images_labels(path):
    # Read the images and labels from a folder containing both dicom files
    for subdir, dirs, files in os.walk(path):
        dcms = glob.glob(os.path.join(subdir, '*.dcm'))
        if len(dcms) == 1:
            structure = dicom.read_file(dcms[0])
            contours = read_structure(structure)
        elif len(dcms) > 1:
            slices = [dicom.read_file(dcm) for dcm in dcms]
            slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
            images = np.stack([s.pixel_array for s in slices], axis=0).astype(np.float32)
            images = images + slices[0].RescaleIntercept
    
    images = normalize(images)
    labels = get_labels(contours, images.shape, slices)
    inplane_scale = slices[0].PixelSpacing[0] / PIXEL_SPACING
    inplane_size = int(np.rint(inplane_scale * slices[0].Rows / 2) * 2)
    z_scale = slices[0].SliceThickness / SLICE_THICKNESS
    z_size = int(np.rint(z_scale * images.shape[0]))
    
    if inplane_size != INPLANE_SIZE or z_scale != 1:
        images = resize(images, (z_size, inplane_size, inplane_size), mode='constant')
        new_labels = np.zeros_like(images, dtype=np.float32)
        for z in range(N_CLASSES):
            roi = resize((labels == z + 1).astype(np.float32), (z_size, inplane_size, inplane_size), mode='constant')
            new_labels[roi >= 0.5] = z + 1
        labels = new_labels
        if inplane_size != INPLANE_SIZE:
            if inplane_size > INPLANE_SIZE:
                crop = int((inplane_size - INPLANE_SIZE) / 2)
                images = images[:, crop : crop + INPLANE_SIZE, crop : crop + INPLANE_SIZE]
                labels = labels[:, crop : crop + INPLANE_SIZE, crop : crop + INPLANE_SIZE]
            else:
                pad = int((INPLANE_SIZE - new_size) / 2)
                images = np.pad(images, ((0, 0), (pad, pad), (pad, pad)), 'constant')
                labels = np.pad(labels, ((0, 0), (pad, pad), (pad, pad)), 'constant')
    
    return images, labels 

def read_images(path):
    for subdir, dirs, files in os.walk(path):
        dcms = glob.glob(os.path.join(subdir, '*.dcm'))
        if len(dcms) > 1:
            slices = [dicom.read_file(dcm) for dcm in dcms]
            slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
            images = np.stack([s.pixel_array for s in slices], axis=0).astype(np.float32)
            images = images + slices[0].RescaleIntercept
    images = normalize(images)
    
    inplane_scale = slices[0].PixelSpacing[0] / PIXEL_SPACING
    inplane_size = int(np.rint(inplane_scale * slices[0].Rows / 2) * 2)
    z_scale = slices[0].SliceThickness / SLICE_THICKNESS
    z_size = int(np.rint(z_scale * images.shape[0]))
    
    if inplane_size != INPLANE_SIZE or z_scale != 1:
        images = resize(images, (z_size, inplane_size, inplane_size), mode='constant')
        if inplane_size != INPLANE_SIZE:
            if inplane_size > INPLANE_SIZE:
                crop = int((inplane_size - INPLANE_SIZE) / 2)
                images = images[:, crop : crop + INPLANE_SIZE, crop : crop + INPLANE_SIZE]
            else:
                pad = int((INPLANE_SIZE - new_size) / 2)
                images = np.pad(images, ((0, 0), (pad, pad), (pad, pad)))
    
    return images 

def _read_projections(folder, indices):
    """Read mayo projections from a folder."""
    datasets = []

    # Get the relevant file names
    file_names = sorted([f for f in os.listdir(folder) if f.endswith(".dcm")])

    if len(file_names) == 0:
        raise ValueError('No DICOM files found in {}'.format(folder))

    file_names = file_names[indices]

    data_array = None

    for i, file_name in enumerate(tqdm.tqdm(file_names,
                                            'Loading projection data')):
        # read the file
        dataset = dicom.read_file(folder + '/' + file_name)

        # Get some required data
        rows = dataset.NumberofDetectorRows
        cols = dataset.NumberofDetectorColumns
        hu_factor = dataset.HUCalibrationFactor
        rescale_intercept = dataset.RescaleIntercept
        rescale_slope = dataset.RescaleSlope

        # Load the array as bytes
        proj_array = np.array(np.frombuffer(dataset.PixelData, 'H'),
                              dtype='float32')
        proj_array = proj_array.reshape([rows, cols], order='F').T

        # Rescale array
        proj_array *= rescale_slope
        proj_array += rescale_intercept
        proj_array /= hu_factor

        # Store results
        if data_array is None:
            # We need to load the first dataset before we know the shape
            data_array = np.empty((len(file_names), cols, rows),
                                  dtype='float32')

        data_array[i] = proj_array[:, ::-1]
        datasets.append(dataset)

    return datasets, data_array 

def _read_dicom_image(self, sv, t):
        filename = self._filename(sv[0],t,sv[1]);
        d = dicom.read_file(filename)
        img = getAlignImg(d);
        #img = d.pixel_array;
        res = np.array(img,dtype=np.float)
        midx = tuple(sv);

        shift  = np.array([0,0]);
        #crop the center
        if res.shape[0]>res.shape[1]:
            s = (res.shape[0]-res.shape[1])//2;
            res = res[s:s+res.shape[1],:];
            shift[1] = s;
        else:
            s = (res.shape[1]-res.shape[0])//2;
            res = res[:,s:s+res.shape[0]];
            shift[0] = s;

        #crop or stretch to the same size
        if self.img_L>0 and (res.shape[0] != self.img_L):
            #print("crop or fill",filename);
            if res.shape[0]>self.img_L:#because h=w after crop
                s = (res.shape[0]-self.img_L)//2;
                res = res[s:s+self.img_L,s:s+self.img_L];
                shift = shift + s;
            else:
                s = (self.img_L-res.shape[0])//2;
                res2 = np.zeros((self.img_L,self.img_L));
                res2[s:s+res.shape[0],s:s+res.shape[0]] = res;
                res = res2;
                shift = shift - s;

        if t==self.time[0]:
            m = self.M[midx];
            self.M[midx] = (m[0],m[1],getPos(m,shift[0],shift[1]));

        return res;

    #example code, check distance of two planes (even if one is rotated and shifted)
    #x =getPos(M1,50,50);
    #ii = getIdx(M2,M1,50,50);
    #y = getPos(M2,ii[0],ii[1]);
    #np.sqrt(np.sum((y-x)**2)); 

def convertimages():
    print("Converting image patient name, birthdate and id to match pinnacle\n")
    global patientname
    global pid
    global dob
    global FrameUID
    global imageslice
    global SeriesUID
    global StudyInstanceUID
    global imageuid
    global patientfolder
    global posrefind
    global imagesetnumber
    global image_orientation
    global flag_noimages

    if not os.path.exists("%s%s/ImageSet_%s.DICOM"%(Inputf,patientfolder, imagesetnumber)):
        #Image set folder not found, need to ignore patient
        #Will want to call a function to be written that will create image set files from the condensed pixel data file
        print("Image files do not exist. Creating image files")
        createimagefiles()
        return
    for file in os.listdir("%s%s/ImageSet_%s.DICOM"%(Inputf,patientfolder, imagesetnumber)):
        if file == '11026.1.img':
            continue
        imageds = dicom.read_file("%s%s/ImageSet_%s.DICOM/%s"%(Inputf, patientfolder, imagesetnumber, file), force=True)
        imageds.PatientsName = patientname
        imageds.PatientID = pid
        imageds.PatientsBirthDate = dob
        imageslice.append(imageds.SliceLocation)
        imageuid.append(imageds.SOPInstanceUID)
        image_orientation = imageds.ImageOrientationPatient
        tempinstuid = imageds.SOPInstanceUID
        posrefind = imageds.PositionReferenceIndicator
        imageds.SOPInstanceUID = tempinstuid
        imageds.FrameOfReferenceUID = FrameUID
        imageds.StudyInstanceUID = StudyInstanceUID
        imageds.SeriesInstanceUID = SeriesUID
        file_meta = Dataset()
        file_meta.TransferSyntaxUID = GTransferSyntaxUID
        file_meta.MediaStorageSOPClassUID = '1.2.840.10008.5.1.4.1.1.2'
        file_meta.MediaStorageSOPInstanceUID = tempinstuid
        file_meta.ImplementationClassUID = gImplementationClassUID
        imageds.file_meta = file_meta
        preamble = getattr(imageds, "preamble", None)
        if not preamble:
            preamble = b'\x00'*128
        currfile = DicomFile(Outputf+"%s/CT.%s.dcm"%(patientfolder, tempinstuid), 'wb')
        currfile.write(preamble)
        currfile.write(b'DICM')
        dicom.write_file(Outputf+"%s/CT.%s.dcm"%(patientfolder,tempinstuid), imageds, False)
        print("Current image: ", file)
####################################################################################################################################################
####################################################################################################################################################                      
      

####################################################################################################################################################
# Function: createimagefiles()
# This function will create dicom image files for each slice using the condensed pixel data from file ImageSet_%s.img
####################################################################################################################################################