Skip to content

Commit

Permalink
DICOM write works - added anonymization
Browse files Browse the repository at this point in the history
  • Loading branch information
mcencini authored and mcencini committed Jan 23, 2024
1 parent 0abbc74 commit f74e438
Show file tree
Hide file tree
Showing 9 changed files with 206 additions and 98 deletions.
44 changes: 28 additions & 16 deletions src/deepmr/io/header/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -103,29 +103,41 @@ def read_acquisition_header(filepath, device="cpu", verbose=False, *args):
print(f"Trajectory shape: (ncontrasts={head.traj.shape[0]}, nviews={head.traj.shape[1]}, nsamples={head.traj.shape[2]}, ndim={head.traj.shape[-1]})")
if head.dcf is not None:
print(f"DCF shape: (ncontrasts={head.dcf.shape[0]}, nviews={head.dcf.shape[1]}, nsamples={head.dcf.shape[2]})")
if head.FA is not None:
if len(np.unique(head.FA)) > 1:
if head.FA is not None:
if len(np.unique(head.FA)) > 1:
if verbose == 2:
print(f"Flip Angle train length: {len(head.FA)}")
else:
FA = float(np.unique(head.FA)[0])
else:
FA = float(np.unique(head.FA)[0])
head.FA = FA
if verbose == 2:
print(f"Constant FA: {round(abs(FA), 2)} deg")
if head.TR is not None:
if len(np.unique(head.TR)) > 1:
if head.TR is not None:
if len(np.unique(head.TR)) > 1:
if verbose == 2:
print(f"TR train length: {len(head.TR)}")
else:
TR = float(np.unique(head.TR)[0])
else:
TR = float(np.unique(head.TR)[0])
head.TR = TR
if verbose == 2:
print(f"Constant TR: {round(TR, 2)} ms")
if head.TE is not None:
if len(np.unique(head.TE)) > 1:
if head.TE is not None:
if len(np.unique(head.TE)) > 1:
if verbose == 2:
print(f"Echo train length: {len(head.TE)}")
else:
TE = float(np.unique(head.TE)[0])
else:
TE = float(np.unique(head.TE)[0])
head.TE = TE
if verbose == 2:
print(f"Constant TE: {round(TE, 2)} ms")
if head.TI is not None and np.allclose(head.TI, 0.0) is False:
if len(np.unique(head.TI)) > 1:
if head.TI is not None and np.allclose(head.TI, 0.0) is False:
if len(np.unique(head.TI)) > 1:
if verbose == 2:
print(f"Inversion train length: {len(head.TI)}")
else:
TI = float(np.unique(head.TI)[0])
else:
TI = float(np.unique(head.TI)[0])
head.TI = TI
if verbose == 2:
print(f"Constant TI: {round(TI, 2)} ms")

tend = time.time()
Expand Down
44 changes: 28 additions & 16 deletions src/deepmr/io/image/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -89,29 +89,41 @@ def read_image(filepath, acqheader=None, device="cpu", verbose=0):
print(f"Trajectory shape: (ncontrasts={head.traj.shape[0]}, nviews={head.traj.shape[1]}, nsamples={head.traj.shape[2]}, ndim={head.traj.shape[-1]})")
if head.dcf is not None:
print(f"DCF shape: (ncontrasts={head.dcf.shape[0]}, nviews={head.dcf.shape[1]}, nsamples={head.dcf.shape[2]})")
if head.FA is not None:
if len(np.unique(head.FA)) > 1:
if head.FA is not None:
if len(np.unique(head.FA)) > 1:
if verbose == 2:
print(f"Flip Angle train length: {len(head.FA)}")
else:
FA = float(np.unique(head.FA)[0])
else:
FA = float(np.unique(head.FA)[0])
head.FA = FA
if verbose == 2:
print(f"Constant FA: {round(abs(FA), 2)} deg")
if head.TR is not None:
if len(np.unique(head.TR)) > 1:
if head.TR is not None:
if len(np.unique(head.TR)) > 1:
if verbose == 2:
print(f"TR train length: {len(head.TR)}")
else:
TR = float(np.unique(head.TR)[0])
else:
TR = float(np.unique(head.TR)[0])
head.TR = TR
if verbose == 2:
print(f"Constant TR: {round(TR, 2)} ms")
if head.TE is not None:
if len(np.unique(head.TE)) > 1:
if head.TE is not None:
if len(np.unique(head.TE)) > 1:
if verbose == 2:
print(f"Echo train length: {len(head.TE)}")
else:
TE = float(np.unique(head.TE)[0])
else:
TE = float(np.unique(head.TE)[0])
head.TE = TE
if verbose == 2:
print(f"Constant TE: {round(TE, 2)} ms")
if head.TI is not None and np.allclose(head.TI, 0.0) is False:
if len(np.unique(head.TI)) > 1:
if head.TI is not None and np.allclose(head.TI, 0.0) is False:
if len(np.unique(head.TI)) > 1:
if verbose == 2:
print(f"Inversion train length: {len(head.TI)}")
else:
TI = float(np.unique(head.TI)[0])
else:
TI = float(np.unique(head.TI)[0])
head.TI = TI
if verbose == 2:
print(f"Constant TI: {round(TI, 2)} ms")

# cast
Expand Down
3 changes: 2 additions & 1 deletion src/deepmr/io/image/common.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,4 +54,5 @@ def _anonymize(head):
head.ref_dicom.PatientName = ""
head.ref_dicom.PatientBirthDate = ""
head.ref_dicom.PatientWeight = ""
head.ref_dicom.PatientID = "Anon"
head.ref_dicom.PatientID = "Anon"
return head
94 changes: 79 additions & 15 deletions src/deepmr/io/image/dicom.py
Original file line number Diff line number Diff line change
Expand Up @@ -73,20 +73,20 @@ def read_dicom(filepath):

# unpack sequence
TI, TE, EC, TR, FA = uContrasts.transpose()

# squeeze
if sorted_image.shape[0] == 1:
sorted_image = sorted_image[0]

# initialize header
header = Header.from_dicom(dsets, firstVolumeIdx)

# update header
header.FA = FA
header.TI = TI
header.TE = TE
header.TR = TR

return sorted_image, header

def write_dicom(filename, image, filepath="./", head=None, series_description="", series_number_offset=0, series_number_scale=1000, rescale=False, anonymize=False, verbose=False):
Expand Down Expand Up @@ -162,6 +162,9 @@ def write_dicom(filename, image, filepath="./", head=None, series_description=""
if head is not None:
transpose = head.transpose
flip = head.flip
else:
transpose = None
flip = None

# cast image to numpy
image, windowRange = _prepare_image(image, transpose, flip, rescale)
Expand Down Expand Up @@ -198,17 +201,47 @@ def write_dicom(filename, image, filepath="./", head=None, series_description=""

# remove constant parameters from header
if head.FA is not None:
if len(np.unique(head.FA)) == 1:
if head.FA.size == 1:
head.ref_dicom.FlipAngle = float(abs(head.FA))
head.FA = None
if head.TR is not None:
if len(np.unique(head.TR)) == 1:
head.TR = None
if head.TE is not None:
if len(np.unique(head.TE)) == 1:
elif len(np.unique(head.FA)) == 1:
head.ref_dicom.FlipAngle = float(abs(head.FA[0]))
head.FA = None
else:
head.FA = list(abs(head.FA).astype(float))
if head.TE is not None and not(np.isinf(np.sum(head.TE))):
if head.TE.size == 1:
head.ref_dicom.EchoTime = float(head.TE)
head.TE = None
elif len(np.unique(head.TE)) == 1:
head.ref_dicom.EchoTime = float(head.TE[0])
head.TE = None
if head.TI is not None:
if len(np.unique(head.TI)) == 1:
else:
head.TE = list(head.TE.astype(float))
else:
head.TE = None
if head.TR is not None and not(np.isinf(np.sum(head.TR))):
if head.TR.size == 1:
head.ref_dicom.RepetitionTime = float(head.TR)
head.TR = None
elif len(np.unique(head.TR)) == 1:
head.ref_dicom.RepetitionTime = float(head.TR[0])
head.TR = None
else:
head.TR= list(head.TR.astype(float))
else:
head.TR = None
if head.TI is not None and not(np.isinf(np.sum(head.TI))):
if head.TI.size == 1:
head.ref_dicom.InversionTime = float(head.TI)
head.TI = None
elif len(np.unique(head.TI)) == 1:
head.ref_dicom.InversionTime = float(head.TI[0])
head.TI = None
else:
head.TI = list(head.TI.astype(float))
else:
head.TI = None

# generate position and slice location
pos, zloc = dicom._make_geometry_tags(affine, shape, resolution)
Expand Down Expand Up @@ -402,19 +435,24 @@ def _cast_to_complex(dsets_in):

if vendor == "Siemens":
return _cast_to_complex_siemens(dsets_in)

if vendor == "DeepMR":
return _cast_to_complex_deepmr(dsets_in)

def _get_vendor(dset):
"""
Get vendor from DICOM header.
"""
if dset.Manufacturer == "GE MEDICAL SYSTEMS":
return "GE"
if "GE MEDICAL SYSTEMS" in dset.Manufacturer.upper():
return "GEHC"

if dset.Manufacturer == "Philips Medical Systems":
if "PHILIPS" in dset.Manufacturer.upper():
return "Philips"

if dset.Manufacturer == "SIEMENS":
if "SIEMENS" in dset.Manufacturer.upper():
return "Siemens"

return "DeepMR"

def _cast_to_complex_ge(dsets_in):
"""
Expand Down Expand Up @@ -557,6 +595,32 @@ def _cast_to_complex_siemens(dsets_in):
# count number of instances
ninstances = img.shape[0]

# assign to pixel array
for n in range(ninstances):
dsets_out[n].pixel_array[:] = 0.0

return img, dsets_out

def _cast_to_complex_deepmr(dsets_in):
"""
Attempt to get magnitude from DeepMR and cast to complex.
"""
# initialize
magnitude = []

# allocate template out
dsets_out = []

# loop over dataset
for dset in dsets_in:
magnitude.append(dset.pixel_array)
dsets_out.append(dset)

img = np.stack(magnitude, axis=0).astype(np.complex64)

# count number of instances
ninstances = img.shape[0]

# assign to pixel array
for n in range(ninstances):
dsets_out[n].pixel_array[:] = 0.0
Expand Down
26 changes: 19 additions & 7 deletions src/deepmr/io/image/nifti.py
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,7 @@ def write_nifti(filename, image, filepath="./", head=None, series_description=""
if head is not None:
head = copy.deepcopy(head)
head.numpy()
# print(head)

# anonymize
if head is not None and anonymize:
Expand Down Expand Up @@ -147,6 +148,9 @@ def write_nifti(filename, image, filepath="./", head=None, series_description=""
if head is not None:
transpose = head.transpose
flip = head.flip
else:
transpose = None
flip = None

# cast image to numpy
image, windowRange = _prepare_image(image, transpose, flip, rescale)
Expand Down Expand Up @@ -177,22 +181,30 @@ def write_nifti(filename, image, filepath="./", head=None, series_description=""
json_dict["SliceThickness"] = str(head.ref_dicom.SliceThickness)
json_dict["EchoNumber"] = [str(n) for n in range(ncontrasts)]
if head.FA is not None:
if len(np.unique(head.FA)) == 1:
if head.FA.size == 1:
json_dict["FlipAngle"] = float(abs(head.FA))
elif len(np.unique(head.FA)) == 1:
json_dict["FlipAngle"] = float(abs(head.FA[0]))
else:
json_dict["FlipAngle"] = list(abs(head.FA).astype(float))
if head.TE is not None and not(np.isinf(head.TE).any()):
if len(np.unique(head.TE)) == 1:
if head.TE is not None and not(np.isinf(np.sum(head.TE))):
if head.TE.size == 1:
json_dict["EchoTime"] = float(head.TE) * 1e-3
elif len(np.unique(head.TE)) == 1:
json_dict["EchoTime"] = float(head.TE[0]) * 1e-3
else:
json_dict["EchoTime"] = list(head.TE.astype(float) * 1e-3)
if head.TR is not None and not(np.isinf(head.TR).any()):
if len(np.unique(head.TR)) == 1:
if head.TR is not None and not(np.isinf(np.sum(head.TR))):
if head.TR.size == 1:
json_dict["RepetitionTime"] = float(head.TR) * 1e-3
elif len(np.unique(head.TR)) == 1:
json_dict["RepetitionTime"] = float(head.TR[0]) * 1e-3
else:
json_dict["RepetitionTime"] = list(head.TR.astype(float) * 1e-3)
if head.TI is not None and not(np.isinf(head.TI).any()):
if len(np.unique(head.TI)) == 1:
if head.TI is not None and not(np.isinf(np.sum(head.TI))):
if head.TI.size == 1:
json_dict["InversionTime"] = float(head.TI) * 1e-3
elif len(np.unique(head.TI)) == 1:
json_dict["InversionTime"] = float(head.TI[0]) * 1e-3
else:
json_dict["InversionTime"] = list(head.TI.astype(float) * 1e-3)
Expand Down
44 changes: 28 additions & 16 deletions src/deepmr/io/kspace/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -179,29 +179,41 @@ def read_rawdata(filepath, acqheader=None, device="cpu", verbose=0):
print(f"Trajectory shape: (ncontrasts={head.traj.shape[0]}, nviews={head.traj.shape[1]}, nsamples={head.traj.shape[2]}, ndim={head.traj.shape[-1]})")
if head.dcf is not None:
print(f"DCF shape: (ncontrasts={head.dcf.shape[0]}, nviews={head.dcf.shape[1]}, nsamples={head.dcf.shape[2]})")
if head.FA is not None:
if len(np.unique(head.FA)) > 1:
if head.FA is not None:
if len(np.unique(head.FA)) > 1:
if verbose == 2:
print(f"Flip Angle train length: {len(head.FA)}")
else:
FA = float(np.unique(head.FA)[0])
else:
FA = float(np.unique(head.FA)[0])
head.FA = FA
if verbose == 2:
print(f"Constant FA: {round(abs(FA), 2)} deg")
if head.TR is not None:
if len(np.unique(head.TR)) > 1:
if head.TR is not None:
if len(np.unique(head.TR)) > 1:
if verbose == 2:
print(f"TR train length: {len(head.TR)}")
else:
TR = float(np.unique(head.TR)[0])
else:
TR = float(np.unique(head.TR)[0])
head.TR = TR
if verbose == 2:
print(f"Constant TR: {round(TR, 2)} ms")
if head.TE is not None:
if len(np.unique(head.TE)) > 1:
if head.TE is not None:
if len(np.unique(head.TE)) > 1:
if verbose == 2:
print(f"Echo train length: {len(head.TE)}")
else:
TE = float(np.unique(head.TE)[0])
else:
TE = float(np.unique(head.TE)[0])
head.TE = TE
if verbose == 2:
print(f"Constant TE: {round(TE, 2)} ms")
if head.TI is not None and np.allclose(head.TI, 0.0) is False:
if len(np.unique(head.TI)) > 1:
if head.TI is not None and np.allclose(head.TI, 0.0) is False:
if len(np.unique(head.TI)) > 1:
if verbose == 2:
print(f"Inversion train length: {len(head.TI)}")
else:
TI = float(np.unique(head.TI)[0])
else:
TI = float(np.unique(head.TI)[0])
head.TI = TI
if verbose == 2:
print(f"Constant TI: {round(TI, 2)} ms")

# cast
Expand Down
Loading

0 comments on commit f74e438

Please sign in to comment.