From 33d3b4f35cbb9cd24246e9f574bfca14766c8880 Mon Sep 17 00:00:00 2001
From: Julius Steiglechner <julius.steiglechner@tuebingen.mpg.de>
Date: Thu, 19 Dec 2024 11:57:54 +0100
Subject: [PATCH] Provide data functionality.
---
src/__init__.py | 0
src/data_flow/__init__.py | 0
src/data_flow/dicom_operations.py | 50 ++++++
src/data_flow/itk_operations.py | 188 ++++++++++++++++++++++
src/data_flow/itk_transform_operations.py | 141 ++++++++++++++++
src/data_flow/nifti_operations.py | 140 ++++++++++++++++
6 files changed, 519 insertions(+)
create mode 100644 src/__init__.py
create mode 100644 src/data_flow/__init__.py
create mode 100755 src/data_flow/dicom_operations.py
create mode 100755 src/data_flow/itk_operations.py
create mode 100644 src/data_flow/itk_transform_operations.py
create mode 100644 src/data_flow/nifti_operations.py
diff --git a/src/__init__.py b/src/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/data_flow/__init__.py b/src/data_flow/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/data_flow/dicom_operations.py b/src/data_flow/dicom_operations.py
new file mode 100755
index 0000000..e373a8b
--- /dev/null
+++ b/src/data_flow/dicom_operations.py
@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""Module provides functionality to deal with DICOM files."""
+
+import datetime as dt
+from pathlib import Path
+from typing import Tuple
+
+import pydicom
+
+
+def read_subjects_age_gender(dcm_path: Path) -> Tuple[int, str]:
+ """
+ Read age and gender for a subject from a DCM file.
+
+ Parameters
+ ----------
+ dcm_path : Path
+ Path to a DCM file.
+
+ Returns
+ -------
+ Tuple[int, str]
+ (age, gender).
+
+ """
+ ds = pydicom.dcmread(dcm_path)
+
+ try:
+ age = int(ds["PatientAge"].value)
+ except AssertionError:
+ birth_date = ds["PatientBirthDate"].value
+ birth_date = dt.date(year=birth_date[:4],
+ month=birth_date[4:6],
+ day=birth_date[6:8])
+ acq_date = ds["AcquisitionDate"].value
+ acq_date = dt.date(year=acq_date[:4],
+ month=acq_date[4:6],
+ day=acq_date[6:8])
+ age = acq_date.year - birth_date.year
+ if acq_date.month < birth_date.month:
+ age -= 1
+ elif (
+ acq_date.month == birth_date.month
+ ) and (acq_date.day < birth_date.day):
+ age -= 1
+
+ gender = ds["PatientSex"].value
+
+ return age, gender
diff --git a/src/data_flow/itk_operations.py b/src/data_flow/itk_operations.py
new file mode 100755
index 0000000..70752c0
--- /dev/null
+++ b/src/data_flow/itk_operations.py
@@ -0,0 +1,188 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This module provides functionality for handling ITK images.
+
+Mainly, there are conversions to and from them.
+"""
+
+from typing import List, Tuple
+
+import itk
+import numpy as np
+
+ITK_COORD_LPI: np.ndarray = np.array([-1, -1, 1], dtype=float)
+
+
+def convert_affine_to_itk_params(
+ affine: np.ndarray,
+) -> Tuple[List, List, itk.Matrix]:
+ """
+ Get itk conform information from affine matrix.
+
+ Parameters
+ ----------
+ affine : np.ndarray
+ Affine matrix.
+
+ Raises
+ ------
+ ValueError
+ If affine spatial dimensionality is not 3.
+
+ Returns
+ -------
+ Tuple[List, List, itk.Matrix]
+ Affine space with itk conform values: spacing, origin, direction.
+
+ """
+ if affine.shape[0] != 4:
+ raise ValueError(
+ f'Spatial dimensionality is not 3, but {affine.shape[0] - 1}.'
+ )
+
+ # Extract affine informations
+ affine_direction = affine[:3, :3]
+ affine_spacing: np.ndarray = np.linalg.norm(affine_direction, axis=0)
+ affine_origin = affine[:3, 3]
+
+ # Convert affine informations
+ affine_origin = affine_origin * ITK_COORD_LPI
+ affine_direction = affine_direction / affine_spacing
+ affine_direction = (affine_direction.T * ITK_COORD_LPI).T
+
+ # Return itk conform values
+ return (
+ affine_spacing.tolist(),
+ affine_origin.tolist(),
+ itk.matrix_from_array(affine_direction),
+ )
+
+
+def convert_itk_params_to_affine(itk_img: itk.Image) -> np.ndarray:
+ """
+ Extract an affine matrix from itk image.
+
+ Parameters
+ ----------
+ itk_img : itk.Image
+ ITK image.
+
+ Raises
+ ------
+ ValueError
+ If image dimensionality does not match.
+
+ Returns
+ -------
+ affine : np.ndarray
+ Affine matrix.
+
+ """
+ if itk_img.GetImageDimension() != 3:
+ raise ValueError(
+ f'Image dimensionality isnot 3, but {itk_img.GetImageDimension()}.'
+ )
+
+ # Extract image information
+ affine_origin = itk_img.GetOrigin()
+ affine_spacing = itk_img.GetSpacing()
+ affine_direction = itk.array_from_matrix(itk_img.GetDirection())
+
+ # Convert affine information
+ affine_origin = affine_origin * ITK_COORD_LPI
+ affine_direction = affine_direction * affine_spacing
+ affine_direction = (affine_direction.T * ITK_COORD_LPI).T
+
+ # Conclude affine information
+ affine = np.eye(4, dtype=float)
+ affine[:3, :3] = affine_direction
+ affine[:3, 3] = affine_origin
+
+ return affine
+
+
+def convert_nifti_array_to_itk(
+ img: np.ndarray,
+ affine: np.ndarray,
+ voxel_type: itk.itkCType = itk.F,
+) -> itk.Image:
+ """
+ Convert a physical 3D-image to an ITK Image.
+
+ Parameters
+ ----------
+ img : np.ndarray
+ Input image effectivly in F ordering.
+ affine : np.ndarray
+ Affine matrix defining the physical cartesian space of image grid.
+ voxel_type : itk.itkCType, optional
+ Data type of voxels. The default is itk.F.
+
+ Raises
+ ------
+ ValueError
+ If image dimensionality is not 3.
+ If affine spatial dimensionality is not 3.
+
+ Returns
+ -------
+ itk_img : itk.Image
+ Converted ITK image.
+
+ """
+ if len(img.shape) != 3:
+ raise ValueError(
+ f'Image dimensionality is not 3, but {len(img.shape)}.')
+
+ if img.flags['C_CONTIGUOUS']:
+ # Is set to F-Odering so that shapes of np and ITK image match.
+ # Note, that NIfTI1 images are effectivly in F-Ordering.
+ img = np.asfortranarray(img)
+
+ affine_spacing, affine_origin, affine_direction = \
+ convert_affine_to_itk_params(affine=affine)
+
+ # Convert to ITK, data is copied
+ itk_img = itk.image_from_array(img)
+ itk_img = itk_img.astype(voxel_type)
+
+ # Set image informations
+ itk_img.SetOrigin(affine_origin)
+ itk_img.SetSpacing(affine_spacing)
+ itk_img.SetDirection(affine_direction)
+
+ return itk_img
+
+
+def convert_itk_img_to_nifti_array(
+ itk_img: itk.Image,
+ dtype: np.dtype = float,
+) -> Tuple[np.ndarray, np.ndarray]:
+ """
+ Convert an ITK image to a Numpy array with a defined physical space.
+
+ Parameters
+ ----------
+ itk_img : itk.Image
+ Image in ITK format.
+ dtype : np.dtype, optional
+ Voxel dtype. The default is float.
+
+ Returns
+ -------
+ img : TYPE
+ Image arrage in F-ordering.
+ affine : TYPE
+ Corresponding affine matrix.
+
+ """
+ affine = convert_itk_params_to_affine(itk_img)
+
+ # Convert to array, data is copied.
+ # This will always be a C ordered array.
+ # Has to be transposed to match NIfTI1 convention of F-Ordering.
+ img = itk.array_from_image(itk_img).T
+ img = img.astype(dtype)
+
+ return img, affine
diff --git a/src/data_flow/itk_transform_operations.py b/src/data_flow/itk_transform_operations.py
new file mode 100644
index 0000000..b89eb4c
--- /dev/null
+++ b/src/data_flow/itk_transform_operations.py
@@ -0,0 +1,141 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Dec 19 11:44:28 2024
+
+@author: jsteiglechner
+"""
+
+from typing import Tuple
+
+import itk
+import numpy as np
+
+
+def extract_elastix_parameters(
+ transform_parameters: itk.ParameterObject,
+) -> Tuple[np.ndarray, np.ndarray]:
+ """
+ Extract parameters of elastix object to affine parameters.
+
+ Parameters
+ ----------
+ transform_parameters : itk.ParameterObject
+ Parameters of operation.
+
+ Returns
+ -------
+ center : TYPE
+ center of rotation point.
+ elx_parameters : TYPE
+ transform parameters.
+
+ """
+ parameter_map = transform_parameters.GetParameterMap(0)
+
+ center = np.array(
+ [float(p) for p in parameter_map['CenterOfRotationPoint']]
+ )
+ elx_parameters = np.array(
+ [float(p) for p in parameter_map['TransformParameters']]
+ )
+
+ return center, elx_parameters
+
+
+def get_transformation_matrix(
+ itk_transform: itk.Transform,
+) -> np.ndarray:
+ """
+ Extract relevant a linear itk transform to generate the affine matrix.
+
+ Parameters
+ ----------
+ itk_transform : itk.Transform
+ Linear transform.
+
+ Returns
+ -------
+ transformation_matrix : TYPE
+ Affine transformation matrix.
+
+ """
+ affine_direction = itk.array_from_matrix(itk_transform.GetMatrix())
+ affine_origin = np.asarray(itk_transform.GetTranslation())
+
+ # Conclude affine information
+ transformation_matrix = np.eye(4, dtype=float)
+ transformation_matrix[:3, :3] = affine_direction
+ transformation_matrix[:3, 3] = affine_origin
+
+ return transformation_matrix
+
+
+def build_euler_transformation(
+ center: np.ndarray,
+ transform: np.ndarray,
+) -> itk.Euler3DTransform:
+ """
+ Build an Euler transformation from elastix parameters.
+
+ Parameters
+ ----------
+ center : np.ndarray
+ Center of rotation point in shape (dim, ).
+ transform : np.ndarray
+ Transform parameters in shape (dim * (dim+1), ).
+
+ Returns
+ -------
+ euler_transform : TYPE
+ Applicable Euler transform.
+
+ """
+ euler_transform = itk.Euler3DTransform[itk.D].New()
+
+ fixed_parameters = itk.OptimizerParameters[itk.D](len(center))
+ for i, p in enumerate(center):
+ fixed_parameters[i] = p
+ euler_transform.SetFixedParameters(fixed_parameters)
+
+ itk_parameters = itk.OptimizerParameters[itk.D](len(transform))
+ for i, p in enumerate(transform):
+ itk_parameters[i] = p
+ euler_transform.SetParameters(itk_parameters)
+
+ return euler_transform
+
+
+def build_affine_transformation(
+ center: np.ndarray,
+ transform: np.ndarray,
+) -> itk.AffineTransform:
+ """
+ Build an Affine transformation from elastix parameters.
+
+ Parameters
+ ----------
+ center : np.ndarray
+ Center of rotation point in shape (dim, ).
+ transform : np.ndarray
+ Transform parameters in shape (dim * (dim+1), ).
+
+ Returns
+ -------
+ affine_transform : TYPE
+ Applicable affine transform.
+
+ """
+ affine_transform = itk.AffineTransform[itk.D, len(center)].New()
+
+ fixed_parameters = itk.OptimizerParameters[itk.D](len(center))
+ for i, p in enumerate(center):
+ fixed_parameters[i] = p
+ affine_transform.SetFixedParameters(fixed_parameters)
+
+ itk_parameters = itk.OptimizerParameters[itk.D](len(transform))
+ for i, p in enumerate(transform):
+ itk_parameters[i] = p
+ affine_transform.SetParameters(itk_parameters)
+
+ return affine_transform
diff --git a/src/data_flow/nifti_operations.py b/src/data_flow/nifti_operations.py
new file mode 100644
index 0000000..066f969
--- /dev/null
+++ b/src/data_flow/nifti_operations.py
@@ -0,0 +1,140 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This module provides functionality for handling NIfTI images with nibabel.
+
+Created on Fri Aug 4 14:19:10 2023
+
+@author: jsteiglechner
+"""
+
+from pathlib import Path
+from typing import Tuple, Union
+
+import nibabel as nib
+import numpy as np
+
+
+def read_nii(filepath: Union[str, Path]) -> Tuple[np.ndarray, np.ndarray]:
+ """
+ Read brain MRI and corresponding affine from .nii file.
+
+ Note that nibabel generates a fortran ordered array.
+
+ Parameters
+ ----------
+ filepath : Union[str, Path]
+ Path to MRI.
+
+ Returns
+ -------
+ image : TYPE
+ Image data.
+ affine_matrix : TYPE
+ Corresponsing affine matrix.
+
+ """
+ nifti = nib.load(filepath)
+ image = nifti.dataobj[:]
+ affine_matrix = nifti.affine
+
+ return image, affine_matrix
+
+
+def save_nii(
+ image: np.ndarray,
+ affine_matrix: np.ndarray,
+ filepath: Union[str, Path],
+ dtype: np.dtype = float,
+) -> None:
+ """
+ Save brain MRI with affine matrix in NIfTI file format.
+
+ Parameters
+ ----------
+ image : np.ndarray
+ Image data.
+ affine_matrix : np.ndarray
+ Corresponding affine matrix.
+ filepath : Union[str, Path]
+ Path where MRI should be saved.
+ dtype : np.dtype, optional
+ The type of the saved array. The default is float.
+
+ Returns
+ -------
+ None
+ DESCRIPTION.
+
+ """
+ nifti = nib.Nifti1Image(image, affine_matrix, dtype=dtype)
+ nib.save(nifti, filepath)
+
+
+def get_datatype(image: np.ndarray) -> np.dtype:
+ """
+ Get datatype of image data.
+
+ Parameters
+ ----------
+ image : np.ndarray
+ NIfTI image dataobj.
+
+ Returns
+ -------
+ datatype : np.dtype
+ numpy datatype.
+
+ """
+ if np.all(np.equal(np.mod(image, 1), 0)):
+ if np.min(image) >= 0:
+ if np.max(image) <= (2**8 - 1):
+ datatype = np.uint8
+ elif np.max(image) <= (2**16 - 1):
+ datatype = np.uint16
+ else:
+ datatype = np.uint32
+ elif (np.min(image) > -2**15) and (np.max(image) < 2**15):
+ datatype = np.int16
+ else:
+ datatype = np.int32
+ else:
+ datatype = np.dtype('float')
+
+ return datatype
+
+
+def convert_to_original_datatype(
+ image: np.ndarray,
+ datatype: np.dtype,
+) -> np.ndarray:
+ """
+ Convert an image to a certain original datatype.
+
+ Parameters
+ ----------
+ image : np.ndarray
+ Image dataobject.
+ datatype : np.dtye
+ original datatype.
+
+ Returns
+ -------
+ image : TYPE
+ Converted image.
+
+ """
+ if datatype == np.dtype('float'):
+ image = image.astype(float)
+ elif datatype == np.uint8:
+ image = nib.casting.float_to_int(image, 'uint8')
+ elif datatype == np.uint16:
+ image = nib.casting.float_to_int(image, 'uint16')
+ elif datatype == np.uint32:
+ image = nib.casting.float_to_int(image, 'uint32')
+ elif datatype == np.int16:
+ image = nib.casting.float_to_int(image, 'int16')
+ else:
+ image = nib.casting.float_to_int(image, 'int32')
+
+ return image
--
GitLab