Files
Bioxel/bioxelnodes/parse.py
Ma Nan 854af4c3bc feat: support more format
feat: support .mrc .map
feat: support .ome.tif
feat: improved user experience for importing
feat: better contrast algorithm
2024-07-22 08:51:12 +08:00

370 lines
13 KiB
Python

from pathlib import Path
import numpy as np
from .exceptions import CancelledByUser
from .utils import get_text_index_str
try:
import SimpleITK as sitk
from pyometiff import OMETIFFReader
import mrcfile
except:
...
SUPPORT_EXTS = ['', '.dcm', '.DCM', '.DICOM', '.ima', '.IMA',
'.bmp', '.BMP',
'.PIC', '.pic',
'.gipl', '.gipl.gz',
'.jpg', '.JPG', '.jpeg', '.JPEG',
'.lsm', '.LSM',
'.tif', '.TIF', '.tiff', '.TIFF',
'.mnc', '.MNC',
'.mrc', '.rec',
'.mha', '.mhd',
'.hdf', '.h4', '.hdf4', '.he2', '.h5', '.hdf5', '.he5',
'.nia', '.nii', '.nii.gz', '.hdr', '.img', '.img.gz',
'.nrrd', '.nhdr',
'.png', '.PNG',
'.vtk',
'.ome.tiff', '.ome.tif',
'.mrc', '.mrc.gz', '.map', '.map.gz']
OME_EXTS = ['.ome.tiff', '.ome.tif',
'.tif', '.TIF', '.tiff', '.TIFF']
MRC_EXTS = ['.mrc', '.mrc.gz', '.map', '.map.gz']
SEQUENCE_EXTS = ['.bmp', '.BMP',
'.jpg', '.JPG', '.jpeg', '.JPEG',
'.tif', '.TIF', '.tiff', '.TIFF',
'.png', '.PNG']
DICOM_EXTS = ['', '.dcm', '.DCM', '.DICOM']
FH_EXTS = ['', '.dcm', '.DCM', '.DICOM', '.ima', '.IMA',
'.gipl', '.gipl.gz',
'.mnc', '.MNC',
'.mrc', '.rec',
'.mha', '.mhd',
'.nia', '.nii', '.nii.gz', '.hdr', '.img', '.img.gz',
'.hdf', '.h4', '.hdf4', '.he2', '.h5', '.hdf5', '.he5',
'.nrrd', '.nhdr',
'.vtk',
'.gz']
def get_ext(filepath: str) -> str:
file_path = Path(filepath)
if file_path.name.endswith(".nii.gz"):
return ".nii.gz"
elif file_path.name.endswith(".img.gz"):
return ".img.gz"
elif file_path.name.endswith(".gipl.gz"):
return ".gipl.gz"
elif file_path.name.endswith(".ome.tiff"):
return ".ome.tiff"
elif file_path.name.endswith(".ome.tif"):
return ".ome.tif"
elif file_path.name.endswith(".mrc.gz"):
return ".mrc.gz"
elif file_path.name.endswith(".map.gz"):
return ".map.gz"
else:
return file_path.suffix
def get_sequence_name(filepath: str) -> str:
ext = get_ext(filepath)
filename = Path(filepath).name.removesuffix(ext)
index: str = get_text_index_str(filename)
return filename.removesuffix(index)
def get_sequence_index(filepath: str) -> int:
ext = get_ext(filepath)
filename = Path(filepath).name.removesuffix(ext)
index: str = get_text_index_str(filename)
return int(index) if index else 0
def collect_sequence(filepath: str):
file_path = Path(filepath).resolve()
files = list(file_path.parent.iterdir())
file_dict = {}
for file in files:
if file.is_file() \
and get_ext(file_path) == get_ext(file) \
and get_sequence_name(file_path) == get_sequence_name(file):
file_dict[get_sequence_index(file)] = file
for key in file_dict.copy().keys():
if not file_dict.get(key+1) \
and not file_dict.get(key-1):
del file_dict[key]
file_dict = dict(sorted(file_dict.items()))
sequence = [str(f) for f in file_dict.values()]
if len(sequence) == 0:
sequence = [str(file_path)]
return sequence
def parse_volumetric_data(filepath: str, series_id="", progress_callback=None):
"""Parse any volumetric data to numpy with shap (T,X,Y,Z,C)
Args:
filepath (str): file path
series_id (str, optional): DICOM series id. Defaults to "".
Returns:
_type_: _description_
"""
ext = get_ext(filepath)
if progress_callback:
progressing = progress_callback(0, "Reading the Data...")
if not progressing:
raise CancelledByUser
is_sequence = False
if ext in SEQUENCE_EXTS:
sequence = collect_sequence(filepath)
if len(sequence) > 1:
is_sequence = True
if ext in MRC_EXTS and not is_sequence:
# TODO: much to do with mrc
with mrcfile.open(filepath) as mrc:
volume = mrc.data
# mrc.print_header()
# print(volume.shape)
# print(mrc.voxel_size)
if mrc.is_single_image():
volume = np.expand_dims(volume, axis=0) # expend frame
volume = np.expand_dims(volume, axis=-1) # expend Z
volume = np.expand_dims(volume, axis=-1) # expend channel
elif mrc.is_image_stack():
volume = np.expand_dims(volume, axis=-1) # expend Z
volume = np.expand_dims(volume, axis=-1) # expend channel
elif mrc.is_volume():
volume = np.expand_dims(volume, axis=0) # expend frame
volume = np.expand_dims(volume, axis=-1) # expend channel
elif mrc.is_volume_stack():
volume = np.expand_dims(volume, axis=-1) # expend channel
name = Path(filepath).name.removesuffix(ext).replace(" ", "-")
shape = volume.shape[1:4]
spacing = (mrc.voxel_size.x, mrc.voxel_size.y, mrc.voxel_size.z)
meta = {
"name": name,
"description": "",
"shape": shape,
"spacing": spacing,
"origin": (0, 0, 0),
"direction": (1, 0, 0, 0, 1, 0, 0, 0, 1),
"frame_count": volume.shape[0],
"channel_count": volume.shape[-1],
"is_oriented": False
}
elif ext in OME_EXTS and not is_sequence:
reader = OMETIFFReader(fpath=filepath)
ome_volume, metadata, xml_metadata = reader.read()
if progress_callback:
progressing = progress_callback(0.5, "Transpose to 'TXYZC'...")
if not progressing:
raise CancelledByUser
try:
# print(ome_volume.shape)
# for key in metadata:
# print(f"{key},{metadata[key]}")
ome_order = metadata['DimOrder BF Array']
except:
ome_order = "TCZYX"
if ome_volume.ndim == 2:
ome_order = ome_order.replace("T", "")\
.replace("C", "").replace("Z", "")
bioxel_order = (ome_order.index('X'),
ome_order.index('Y'))
volume = np.transpose(ome_volume, bioxel_order)
volume = np.expand_dims(volume, axis=0) # expend frame
volume = np.expand_dims(volume, axis=-1) # expend Z
volume = np.expand_dims(volume, axis=-1) # expend channel
elif ome_volume.ndim == 3:
# -> XYZC
ome_order = ome_order.replace("T", "").replace("C", "")
bioxel_order = (ome_order.index('X'),
ome_order.index('Y'),
ome_order.index('Z'))
volume = np.transpose(ome_volume, bioxel_order)
volume = np.expand_dims(volume, axis=0) # expend frame
volume = np.expand_dims(volume, axis=-1) # expend channel
elif ome_volume.ndim == 4:
# -> XYZC
ome_order = ome_order.replace("T", "")
bioxel_order = (ome_order.index('X'),
ome_order.index('Y'),
ome_order.index('Z'),
ome_order.index('C'))
volume = np.transpose(ome_volume, bioxel_order)
volume = np.expand_dims(volume, axis=0) # expend frame
elif ome_volume.ndim == 5:
# -> TXYZC
bioxel_order = (ome_order.index('T'),
ome_order.index('X'),
ome_order.index('Y'),
ome_order.index('Z'),
ome_order.index('C'))
volume = np.transpose(ome_volume, bioxel_order)
shape = volume.shape[1:4]
try:
spacing = (metadata['PhysicalSizeX'],
metadata['PhysicalSizeY'],
metadata['PhysicalSizeZ'])
except:
spacing = (1, 1, 1)
name = Path(filepath).name.removesuffix(ext).replace(" ", "-")
meta = {
"name": name,
"description": "",
"shape": shape,
"spacing": spacing,
"origin": (0, 0, 0),
"direction": (1, 0, 0, 0, 1, 0, 0, 0, 1),
"frame_count": volume.shape[0],
"channel_count": volume.shape[-1],
"is_oriented": False
}
else:
if ext in DICOM_EXTS:
dir_path = Path(filepath).resolve().parent
reader = sitk.ImageSeriesReader()
reader.MetaDataDictionaryArrayUpdateOn()
reader.LoadPrivateTagsOn()
series_files = reader.GetGDCMSeriesFileNames(
str(dir_path), series_id)
reader.SetFileNames(series_files)
itk_volume = reader.Execute()
# for k in reader.GetMetaDataKeys(0):
# v = reader.GetMetaData(0, k)
# print(f'({k}) = = "{v}"')
def get_meta(key):
try:
stirng = reader.GetMetaData(0, key).removesuffix(" ")
if stirng in ["No study description",
"No series description"]:
return None
else:
return stirng
except:
return None
study_description = get_meta("0008|1030")
series_description = get_meta("0008|103e")
series_modality = get_meta("0008|0060")
name = study_description or dir_path.name
if series_description and series_modality:
description = f"{series_description}-{series_modality}"
elif series_description:
description = series_description
elif series_modality:
description = series_modality
else:
description = ""
name = name.replace(" ", "-")
description = description.replace(" ", "-")
elif ext in SEQUENCE_EXTS:
itk_volume = sitk.ReadImage(sequence)
name = get_sequence_name(filepath).replace(" ", "-")
description = ""
else:
itk_volume = sitk.ReadImage(filepath)
name = Path(filepath).name.removesuffix(ext).replace(" ", "-")
description = ""
# for key in itk_volume.GetMetaDataKeys():
# print(f"{key},{itk_volume.GetMetaData(key)}")
if progress_callback:
progressing = progress_callback(0.5, "Transpose to 'TXYZC'...")
if not progressing:
raise CancelledByUser
if itk_volume.GetDimension() == 3:
itk_volume = sitk.DICOMOrient(itk_volume, 'RAS')
volume = sitk.GetArrayFromImage(itk_volume)
# transpose ijk to kji
if volume.ndim == 4:
volume = np.transpose(volume, (2, 1, 0, 3))
else:
volume = np.transpose(volume)
volume = np.expand_dims(volume, axis=-1) # expend channel
volume = np.expand_dims(volume, axis=0) # expend frame
meta = {
"name": name,
"description": description,
"shape": tuple(itk_volume.GetSize()),
"spacing": tuple(itk_volume.GetSpacing()),
"origin": tuple(itk_volume.GetOrigin()),
"direction": tuple(itk_volume.GetDirection()),
"frame_count": volume.shape[0],
"channel_count": volume.shape[-1],
"is_oriented": True
}
if itk_volume.GetDimension() == 4:
# FIXME: not sure...
direction = np.array(itk_volume.GetDirection())
direction = direction.reshape(3, 3) if itk_volume.GetDimension() == 3 \
else direction.reshape(4, 4)
direction = direction[1:, 1:]
direction = tuple(direction.flatten())
volume = sitk.GetArrayFromImage(itk_volume)
if volume.ndim == 5:
volume = np.transpose(volume, (0, 3, 2, 1, 4))
else:
volume = np.transpose(volume, (0, 3, 2, 1))
volume = np.expand_dims(volume, axis=-1)
meta = {
"name": name,
"description": description,
"shape": tuple(itk_volume.GetSize()[:3]),
"spacing": tuple(itk_volume.GetSpacing()[:3]),
"origin": tuple(itk_volume.GetOrigin()[:3]),
"direction": direction,
"frame_count": volume.shape[0],
"channel_count": volume.shape[-1],
"is_oriented": False
}
return volume, meta