Commit be3569f2 authored by Joe Lyons's avatar Joe Lyons 🚴
Browse files

Resolve "Add code to handle VASO fMRI sequence"

parent 0bc2d8a2
#!/usr/bin/python3
# New module for YNiC imaging archive service code
# Edited by ADG/JL Jan 2020 to handle Siemens VASO data
from copy import deepcopy
from os import chmod, walk, unlink
from os.path import join, basename
from shutil import copy2
from tempfile import mkdtemp
import tarfile
from .errors import YIASMixedTypeError, YIASNoSeriesError, YIASDuplicateError
from .utils import sanitise_string
import tarfile
import nibabel.nicom.dicomwrappers
# We filter UserWarning messages from nibabel.
# It warns that the DICOM readers are highly experimental, unstable
# and only work for Siemens time-series at the moment, which is fine.
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
try:
from pydicom import read_file
......@@ -24,8 +30,10 @@ try:
except ImportError: # pragma: nocover
from dicom.errors import InvalidDicomError
__all__ = []
def anonymise_dicom(dcm):
"""
Anonymise the relevant Patient fields in a DICOM file object
......@@ -38,8 +46,10 @@ def anonymise_dicom(dcm):
return dcm
__all__.append('anonymise_dicom')
def incoming_from_tarball(tarname, tmpdir):
"""
Routine which creates an incoming directory (inside tmpdir) from a
......@@ -56,8 +66,10 @@ def incoming_from_tarball(tarname, tmpdir):
return i
__all__.append('incoming_from_tarball')
class IncomingHandler(object):
"""
This class is responsible for parsing a directory containing DICOM
......@@ -68,7 +80,7 @@ class IncomingHandler(object):
"""Directory which we are handling"""
self.directory = directory
"A list of filenames which were could not be parsed as DICOM files"""
"""A list of filenames which were could not be parsed as DICOM files"""
self.bad_files = []
"""A list of filenames which were DICOM files but contained duplicate SOPInstanceUIDs"""
......@@ -80,6 +92,9 @@ class IncomingHandler(object):
# Parse all DICOM files in our incoming directory
self.parse_directory()
# Now split directory for VASO data if necessary
self.split_studies_VASO()
# Now split directory by echo times if necessary
self.split_studies_by_echo()
......@@ -117,6 +132,12 @@ class IncomingHandler(object):
for studyuid, study in list(self.studies.items()):
study.split_series_by_echo()
def split_studies_VASO(self):
"""
Iterate over all known series and split if VASO data
"""
for studyuid, study in list(self.studies.items()):
study.split_series_VASO()
###########################################################################
# Routines for extracting study information
......@@ -125,6 +146,7 @@ class IncomingHandler(object):
"""
Iterator over the studies which tries to go in a logical order
"""
study_keys = list(self.studies.keys())
study_numbers = [self.studies[key].StudyDateTime for key in study_keys]
sorted_indices = [x[0] for x in sorted(enumerate(study_numbers), key=lambda y:y[1])]
......@@ -132,6 +154,7 @@ class IncomingHandler(object):
for idx in sorted_indices:
yield self.studies[study_keys[idx]]
__all__.append('IncomingHandler')
......@@ -165,7 +188,6 @@ class DICOMStudy(object):
"""
self.study_description_override = None
def __str__(self):
return self.str_indent(0)
......@@ -224,6 +246,24 @@ class DICOMStudy(object):
self.series.update(old_series.split_series_by_echo_time())
def split_series_VASO(self):
"""
Walks through all series in the study and check if VASO data
split is necessary
"""
# Walk through our series and, where necessary, split by echo time
# We can't iterate over the dictionary as we may remove things as we go
for uid in list(self.series.keys()):
curr_series_name = self.series[uid].SeriesDescription
# Break out if series is not VASO
if "VASO" not in curr_series_name:
continue
# If we have multiple echos, split the series and remove the old one
old_series = self.series.pop(uid)
self.series.update(old_series.split_series_VASO_samples())
###########################################################################
# Routines for extracting series information
......@@ -355,8 +395,10 @@ class DICOMStudy(object):
return self.get_first_series().SanitisedStudyDescription
__all__.append('DICOMStudy')
class DICOMSeries(object):
"""
This class is responsible for representing a single DICOM series
......@@ -432,6 +474,7 @@ class DICOMSeries(object):
'ORIGINAL': self.IsOriginal,
'PHYSIO': self.IsPhysiological,
'NIFTI': self.ShouldCreateNIFTI,
'VASO': self.IsVASO,
'ISSS': self.IsISSS}
flags = ', '.join([ x for x in sorted(list(flags.keys())) if flags[x] ])
......@@ -489,6 +532,64 @@ class DICOMSeries(object):
return series_map
def split_series_VASO_samples(self):
"""
If the series has more than one sample at each slice (e.g. for VASO)
data, return multiple series objects in a dictionary (each with a
new UID key) where the data has been split out.
If there is only one sample per slice, nothing will be done and None will be
returned.
The grouping info is held in the ICE_dims flag of the CSA header, at index 3 and 4.
When the text string is split on underscore ( _ )
e.g. ICE_Dims: '1_1_1_2_292_1_1_1_1_1_1_1_357'
index 3 is the group and 4 is the volume number in that group
"""
# Create a holder to capture the grouping information for each volume
grouping_info = []
groups = []
for imageuid, image in list(self.images.items()):
imagename = image.filename
# Load data
img = nibabel.nicom.dicomwrappers.wrapper_from_file(imagename)
# The grouping info is held in the ICE_dims flag of the CSA header, at index 3 and 4
ice_dims = img.csa_header.get('tags').get('ICE_Dims').get('items')[0].split('_')[3:5]
grouping_info.append(ice_dims)
# Keep track of number of groups
if ice_dims[0] not in groups:
groups.append(ice_dims[0])
series_map = {}
sample_indices = {}
# groups is sorted, so we sort them this way
for group_idx, thisGroup in enumerate(groups):
new_uid = self.seriesUID + '.%d' % (group_idx)
series_map[new_uid] = DICOMSeries(new_uid)
sample_indices[thisGroup] = new_uid
# 1-index the echo descriptions to be more user-friendly
series_map[new_uid].extraSeriesNumber = 'S%d' % (group_idx+1)
series_map[new_uid].extraDescription = '_SAMPLE%d' % (group_idx+1)
## series_map[new_uid].multi_echo = (echo_idx+1, len(self.EchoTimes)) # not relevant?
pos_cnt = 0
for imageuid, image in list(self.images.items()):
group_to_join = str(grouping_info[pos_cnt][0])
series_map[ sample_indices[group_to_join] ].add_dicom(image)
pos_cnt += 1
return series_map
def deduplicate_images(self):
"""
Sorts through instance file names and checks whether there are any
......@@ -532,7 +633,6 @@ class DICOMSeries(object):
copy2(imagename, join(dirname, fname))
chmod(join(dirname, fname), mode)
def anonymise_series(self, dirname, mode=0o444):
"""
Anonymise series into a pre-created directory
......@@ -776,6 +876,20 @@ class DICOMSeries(object):
return True
@property
def IsVASO(self):
"""
Check whether this was acquired by one of the VASO sequences
Unfortunately, we have to do this based on the series name as
there's nothing in the header to tell us.
"""
try:
name = self.get_first_property('SeriesDescription')
return name.lower().startswith('vaso')
except ValueError:
return False
@property
def IsISSS(self):
"""
......@@ -949,7 +1063,6 @@ class DICOMSeries(object):
return self.get_first_property('SeriesDescription') + \
self.extraDescription
@property
def SanitisedSeriesDescription(self):
"""
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment