classic_to_enhanced_Philips_fMRI.py

#!/usr/bin/env python
#
# Convert multiple classic DICOMs to a single enhanced DICOM. Works only for
# Philips fMRIs with a single stack that can be organized with dimensions
# stack / slice / time. Required other files were generated by
# dicom_standard_to_yaml.py:
#    enhanced-mr-image-root.yaml  
#    shared_tags_fmri.yaml   (manually edited for fMRI from enhanced-mr-image.yaml)
#
# Probably requires python 3.
# Quite slow: ~1 hour for a typical 25k frames multiband fMRI.

import argparse
import copy
import difflib
import functools
import numpy
import os
import pandas
import pydicom
import uuid
import datetime
import yaml
import sys

# DICOM UID root. To generate compliant DICOMs, you must use one that is uniquely assigned to you,
# and maintain uniqueness within your subset of the numbering space
# https://www.medicalconnections.co.uk/FreeUID/
UIDROOT = 'EDIT_ME'

# Implementation class info for this program. It needs its own UID. Version name just needs to be
# some short informative string like 'MySite cl2enh'
IMPLEMENTATION_CLASS_UID = UIDROOT + '.1'
IMPLEMENTATION_VERSION_NAME = 'EDIT_ME'

# Strange value that shows up in Philips diffusion tags sometimes instead of zero. We'll replace 
# it with zero
WACKY_DIFFUSION_VALUES = [
    float(7230000000000000073493057058404159389631884489254582023887810025787823751168),
    ]

# All tags OTHER than these are required to be the same for all classic files (frames).
# DiffusionBValue and DiffusionGradientOrientation are left out here because diffusion
# support is not implemented. This list is fMRI specific
VARYING_TAGS = [
    'MediaStorageSOPInstanceUID',
    'SOPInstanceUID',
    'InstanceCreationTime',
    'InstanceNumber',
    'InStackPositionNumber',
    'ImagePositionPatient',
    'TemporalPositionIdentifier',
    'AcquisitionTime',
    'ContentTime',
    'DimensionIndexValues',
    'WindowCenter',
    'WindowWidth',
    ]


def parse_args() -> tuple:
    """ Get command line arguments """
    parser = argparse.ArgumentParser()
    parser.add_argument('--dicom_dir', required=True) 
    parser.add_argument('--out_dir', default='.', required=False)
    args = parser.parse_args()
    assert (os.path.isdir(args.dicom_dir)), 'dicom_dir does not exist'
    assert (os.path.isdir(args.out_dir)), 'out_dir does not exist'
    return args.dicom_dir, args.out_dir


def get_file_list(dicom_dir: str) -> list:
    """ Recursively find files in a directory """
    files = list()
    for r,d,f in os.walk(dicom_dir):
        files = files + [os.path.join(r,x) for x in f]
    return files


def read_dicoms(files:list) -> tuple:
    """
    Load a list of DICOM files, return a tuple containing a list of filenames
    for those that are actually DICOM and a matching list of pydicom.dataset.Dataset
    """
    dicom_files = list()
    datasets = list()
    for f in files:
        if not pydicom.misc.is_dicom(f):
            print(f'Skipping non-DICOM file {f}')
            continue
        ds = pydicom.dcmread(f,stop_before_pixels=False)
        datasets.append(ds)
        dicom_files.append(f)
    return dicom_files, datasets
    

def organize_dicoms_fmri(dicom_files:list, datasets:list) -> tuple:
    """
    Specific to fMRI with stack/slice/time dimension organization and 
    a single stack
    """

    # Clean up the datasets a bit
    for ds in datasets:
        
        # Remove private tags. But, retain the scaling value in one specific
        # Philips field 0x2005100e
        # Scale slope is value 0e in the Philips MR Imaging DD 001 private block
        tg = ds.private_block(0x2005, 'Philips MR Imaging DD 001').get_tag(0x0e)
        scale_slope = ds[tg].value
        
        ds.remove_private_tags()
        
        # To create a private block with a specific number. Philips classic has
        # this at the root level in Philips MR Imaging DD 001 (0x10). For the enhanced,
        # we will put it in SharedFunctionalGroupsSequence as Philips MR Imaging DD 005
        # (0x14) as 2005,140f -> 2005,100e. So we will need to create both of 001 and
        # 005. private_block method seems to just create them in order from 0x10
        bk1 = ds.private_block(0x2005, 'Philips MR Imaging DD 001', create=True)        
        bk1.add_new(0x0e, 'FL', scale_slope)
        
        # Some non-diffusion frames have weird values in diffusion fields
        fix_diffusion_values(ds)
        
    # Find or compute fmri-specific frame info
    temporal_position = list()
    slice_position = list()
    instance_uid = list()
    image_position_patient = list()
    for ds in datasets:
        
        # We track the original instance UID, but don't use it
        instance_uid.append(str(ds.SOPInstanceUID))

        # Grab the original temporal position info
        temporal_position.append(int(ds.TemporalPositionIdentifier))
        
        # Spatial position
        image_position_patient.append(ds.ImagePositionPatient)
        
        # Compute the signed slice position on the through-slice axis
        # https://nipy.org/nibabel/dicom/dicom_orientation.html
        # https://itk.org/pipermail/insight-users/2003-September/004762.html
        iop = numpy.array(ds.ImageOrientationPatient)
        ipp = numpy.array(ds.ImagePositionPatient)
        slice_normal = numpy.cross(iop[0:3], iop[3:6])
        slice_position.append(numpy.dot(slice_normal, ipp))
        
    # Check counts
    utemporal_position = numpy.sort(numpy.unique(temporal_position))
    num_times = len(utemporal_position)
    uslice_position = numpy.sort(numpy.unique(slice_position))
    num_slices = len(uslice_position)
    numF = len(datasets)
    print(f'Found {num_times} time points, {num_slices} slices, {numF} frames')
    if len(utemporal_position)<2:
        raise Exception('Only 1 temporal position - not an fmri series?')
    if numF != num_times*num_slices:
        raise Exception(f'Time and slice not completely sampled - expected {num_times*num_slices} frames')
    
    # Compute InStackPositionNumber
    slice_position_idx = [numpy.where(uslice_position==x)[0][0]+1 for x in slice_position]
    
    # File info dataframe
    dicom_info = pandas.DataFrame({
        'dicom_file': dicom_files,
        'image_position_patient': image_position_patient,
        'temporal_position': temporal_position,
        'instance_uid': instance_uid,
        'slice_position': slice_position,
        'slice_position_idx': slice_position_idx,
        })
    
    # Resort info and datasets by time, then slice
    dicom_info.sort_values(['temporal_position','slice_position'],inplace=True)
    datasets = [datasets[x] for x in dicom_info.index.values]
   
    # Create PerFrameFunctionalGroupsSequence items (fmri specific)
    frames = list()
    for r in range(len(datasets)):
        
        frame = pydicom.dataset.Dataset()
        
        # Create the PlanePositionSequence 
        pps = pydicom.dataset.Dataset()
        pps.ImagePositionPatient = datasets[r].ImagePositionPatient
        frame.PlanePositionSequence = [pps]
        
        # Create the FrameContentSequence
        fcs = pydicom.dataset.Dataset()
        fcs.FrameAcquisitionDateTime = pydicom.valuerep.DT(datasets[r].AcquisitionDate)
        fcs.FrameReferenceDateTime = fcs.FrameAcquisitionDateTime
        fcs.FrameAcquisitionDuration = datasets[r].AcquisitionDuration
        fcs.StackID = '1'
        fcs.InStackPositionNumber = dicom_info.slice_position_idx.iloc[r]
        fcs.DimensionIndexValues = [
            1,fcs.InStackPositionNumber,datasets[r].TemporalPositionIdentifier
            ]
        fcs.TemporalPositionIndex = datasets[r].TemporalPositionIdentifier
        frame.FrameContentSequence = [fcs]
        
        # Create the InstanceNumber
        frame.InstanceNumber =  pydicom.valuerep.IS(r + 1)
        
        # FrameVOILUTSequence
        fvs = pydicom.dataset.Dataset()
        fvs.WindowCenter = datasets[r].WindowCenter
        fvs.WindowWidth = datasets[r].WindowWidth
        frame.FrameVOILUTSequence = [fvs]
        
        # TemporalPositionSequence - we'll just leave this one out
        
        # Add new frame to the list
        frames.append(frame)      
        
    return dicom_info, datasets, frames


def callback_tags_delete(ds, el):
    """
    Callback function for verify_matching_tags(). Works on the original
    dataset so no need for a return value
    """
    if el.tag in VARYING_TAGS:
        del ds[el.tag]


def verify_matching_tags(dicom_files:list, datasets:list):
    """ 
    Verify all tags match for a list of pydicom.dataset.Dataset, except
    the exceptions in VARYING_TAGS.
     """
    shared_tags = copy.deepcopy(datasets[0])
    shared_tags.walk(callback_tags_delete)
    shared_tags.file_meta.walk(callback_tags_delete)
    strrep_shared = str(shared_tags).split("\n")
        
    strreps = list()
    for ds in datasets:
        thisds = copy.deepcopy(ds)
        thisds.walk(callback_tags_delete)
        thisds.file_meta.walk(callback_tags_delete)
        lines = str(thisds).split("\n")
        strreps.append(lines)

    # https://pydicom.github.io/pydicom/stable/auto_examples/plot_dicom_difference.html
    for s,strrep in enumerate(strreps):
        diff_g = difflib.context_diff(
            strrep,strrep_shared,n=0,fromfile=dicom_files[0],tofile=dicom_files[s]
            )
        diffs = [x for x in diff_g]
        if diffs:
            print('\n-------------------------------------------------------------------')
            for line in diffs:
                print(line)
            print('-------------------------------------------------------------------\n')    
            raise Exception(f'Value mismatch in DICOM files')

    
def fix_diffusion_values(ds:pydicom.dataset.Dataset):
    """
    Some Philips non-diffusion scans get a strange value in diffusion-related
    fields. We'll reset it to zero if we find it
    """
    if ds.DiffusionBValue in WACKY_DIFFUSION_VALUES:
        ds.DiffusionBValue = 0.0;
    for d in range(3):
        if ds.DiffusionGradientOrientation[d] in WACKY_DIFFUSION_VALUES:
            ds.DiffusionGradientOrientation[d] = 0.0


def make_instance_uid(prefix:str=UIDROOT) -> str:
    """ Use current time to microseconds and random number to generate unique ID """
    timestamp = datetime.datetime.strftime(datetime.datetime.now(),"%Y%H%M%S%f")
    prefix = f'{prefix}.{timestamp}.'
    uid = pydicom.uid.generate_uid(prefix)
    return uid
    
    
def make_enhanced_root(root_tags_file:str, datasets:list) -> pydicom.dataset.Dataset:
    """
    Create and populate the root level fields of an enhanced DICOM
    based on info in a classic. Update several fields to be compatible
    with enhanced format. Return the enhanced dataset. Assumes fMRI with
    dimensions stack/slice/time.
    """
    with open(root_tags_file,'r') as f:
        root_tags = yaml.safe_load(f)
    
    enhanced = copy.deepcopy(datasets[0])
    enhanced.remove_private_tags()
    enhanced.clear()
    
    # Get as much as we can from the classic. Make sure we don't get functional
    # group sequences (though they shouldn't be there)
    root_tags.remove('SharedFunctionalGroupsSequence')
    root_tags.remove('PerFrameFunctionalGroupsSequence')
    for tag in root_tags:
        if datasets[0].get(tag):
            enhanced[tag] = datasets[0][tag]
        
    # Add the .1 that indicates Enhanced to the SOP class
    enhanced.SOPClassUID = enhanced.SOPClassUID + '.1'
    
    # Create an instance UID for the file
    enhanced.SOPInstanceUID = make_instance_uid()
    
    # Add back InstanceNumber
    enhanced.InstanceNumber = 1
    
    # Update the meta info
    enhanced.file_meta.MediaStorageSOPClassUID = enhanced.SOPClassUID
    enhanced.file_meta.MediaStorageSOPInstanceUID = enhanced.SOPInstanceUID
    enhanced.file_meta.ImplementationClassUID = IMPLEMENTATION_CLASS_UID
    enhanced.file_meta.ImplementationVersionName = IMPLEMENTATION_VERSION_NAME
    
    # Number of frames
    enhanced.NumberOfFrames = len(datasets)
    
    # Estimate acquisition duration
    enhanced.AcquisitionDuration = enhanced.NumberOfFrames * enhanced.AcquisitionDuration
    
    # Dimension UID specific to this file
    DimensionOrganizationUID = make_instance_uid()
    
    # Dimension info. Only Stack / Slice / Time is supported
    DimensionOrganizationSequence = pydicom.dataset.Dataset()
    DimensionOrganizationSequence.DimensionOrganizationUID = DimensionOrganizationUID
    enhanced.DimensionOrganizationSequence = [DimensionOrganizationSequence]
    
    d_stack = pydicom.dataset.Dataset()
    d_stack.DimensionOrganizationUID = DimensionOrganizationUID
    d_stack.DimensionIndexPointer = pydicom.tag.Tag('StackID')
    d_stack.FunctionalGroupPointer = pydicom.tag.Tag('FrameContentSequence')
    d_stack.DimensionDescriptionLabel = 'Stack ID'

    d_slice = pydicom.dataset.Dataset()
    d_slice.DimensionOrganizationUID = DimensionOrganizationUID
    d_slice.DimensionIndexPointer = pydicom.tag.Tag('InStackPositionNumber')
    d_slice.FunctionalGroupPointer = pydicom.tag.Tag('FrameContentSequence')
    d_slice.DimensionDescriptionLabel = 'In Stack Position Number'
    
    d_time = pydicom.dataset.Dataset()
    d_time.DimensionOrganizationUID = DimensionOrganizationUID
    d_time.DimensionIndexPointer = pydicom.tag.Tag('TemporalPositionIndex')
    d_time.FunctionalGroupPointer = pydicom.tag.Tag('FrameContentSequence')
    d_time.DimensionDescriptionLabel = 'Temporal Position Index'
    
    enhanced.DimensionIndexSequence = [d_stack, d_slice, d_time]
    
    return enhanced


def insert_shared_tags(dataset:pydicom.dataset.Dataset, tags:dict):
    """
    Recursively create nested sequences from the tags list, and insert
    values from dataset when we find them
    """
    # These are equivalent and return the value:
    #    enhanced.SharedFunctionalGroupsSequence
    #    enhanced['SharedFunctionalGroupsSequence'].value
    # But this returns the data element:
    #    enhanced['SharedFunctionalGroupsSequence']
    ds = pydicom.dataset.Dataset()
    for tag, val in tags.items():    
        if val:  # If it's a sequence, we need to recurse
            ds.add_new(tag, 'SQ', [insert_shared_tags(dataset, val)])
        elif dataset.get(tag):  # Otherwise it's a value
            ds.add_new(tag, pydicom.datadict.dictionary_VR(tag), dataset[tag].value)
    return ds


def make_enhanced_shared(
        enhanced:pydicom.dataset.Dataset,
        dataset:pydicom.dataset.Dataset,
        shared_tags
        ):
    """
    Create the SharedFunctionalGroupsSequence in an enhanced DICOM. Add all sequences
    in shared_tags, and copy any elements that exist in dataset.
    """
    with open(shared_tags,'r') as f:
        shared_tags = yaml.safe_load(f)

    # Philips private tag that holds scale slope, in classic
    scale_slope = dataset.get_private_item(0x2005, 0x0e, 'Philips MR Imaging DD 001').value

    # Build the inner seq for enhanced and insert the scale slope value
    ssseq = pydicom.dataset.Dataset()
    sbk1 = ssseq.private_block(0x2005, 'Philips MR Imaging DD 001', create=True)        
    sbk1.add_new(0x0e, 'FL', scale_slope)
    
    # Build the outer blocks for enhanced. pydicom won't let us specify the block
    # number, so we have to make all five to get the 005 block numbered as 0x14 
    # as it is in Philips practice.
    ebk1 = enhanced.private_block(0x2005, 'Philips MR Imaging DD 001', create=True)
    ebk2 = enhanced.private_block(0x2005, 'Philips MR Imaging DD 002', create=True)
    ebk3 = enhanced.private_block(0x2005, 'Philips MR Imaging DD 003', create=True)
    ebk4 = enhanced.private_block(0x2005, 'Philips MR Imaging DD 004', create=True)
    ebk5 = enhanced.private_block(0x2005, 'Philips MR Imaging DD 005', create=True)
    ebk5.add_new(0x0f, 'SQ', [ssseq])
    
    # Public tags for the shared functional groups sequence
    ds = insert_shared_tags(dataset, shared_tags)
    enhanced.SharedFunctionalGroupsSequence = [ds]


def make_enhanced_frames(
        enhanced:pydicom.dataset.Dataset,
        datasets:list,
        frames:list
        ):
    """
    Add FrameContentSequence and PixelData to an enhanced DICOM, sourced from
    datasets and frames from organize_dicoms_fmri()
    """
    PixelData = datasets[0].PixelData
    PerFrameFunctionalGroupsSequence = frames[0:1]
    
    for r in range(1,len(datasets)):
        PixelData = PixelData + datasets[r].PixelData
        PerFrameFunctionalGroupsSequence.append(frames[r])

    enhanced.PixelData = PixelData
    enhanced.PerFrameFunctionalGroupsSequence = PerFrameFunctionalGroupsSequence


def main():
    dicom_dir, out_dir = parse_args()
    files = get_file_list(dicom_dir)
    print(f'Number of files: {len(files)}')
    dicom_files, datasets = read_dicoms(files)
    dicom_info, datasets, frames = organize_dicoms_fmri(dicom_files, datasets)
    verify_matching_tags(dicom_info.dicom_file, datasets)
    enhanced = make_enhanced_root('enhanced-mr-image-root.yaml', datasets)
    make_enhanced_shared(enhanced, datasets[0], 'shared_tags_fmri.yaml')
    make_enhanced_frames(enhanced, datasets, frames)
    enhanced.save_as(os.path.join(out_dir, enhanced.SOPInstanceUID))


if __name__ == '__main__':
    main()