Extract DICOM conversion to util module and add proper affine matrix (#396)

arhowe00 · arhowe00 · commit a7eae350dc45 · 2025-11-20T11:54:03.000-05:00
- Move is_dicom_file_or_directory() and convert_dicom_to_nifti() from
  database.py to new util/volume_conversion.py
- Implement extract_affine_from_dicom() using ImageOrientationPatient,
  ImagePositionPatient, and PixelSpacing from DICOM headers
- Create tests/test_volume_conversion.py with unit tests for conversion
  functions
- Remove conversion function tests from test_database.py, keep
  integration tests for write_volume
- Update imports in database.py and test_database.py
diff --git a/src/openlifu/db/database.py b/src/openlifu/db/database.py
@@ -12,13 +12,15 @@
 
 import h5py
 import nibabel as nib
-import numpy as np
-import pydicom
 
 from openlifu.nav.photoscan import Photoscan, load_data_from_photoscan
 from openlifu.plan import Protocol, Run, Solution
 from openlifu.util.json import PYFUSEncoder
 from openlifu.util.types import PathLike
+from openlifu.util.volume_conversion import (
+    convert_dicom_to_nifti,
+    is_dicom_file_or_directory,
+)
 from openlifu.xdc import Transducer, TransducerArray
 from openlifu.xdc.util import load_transducer_from_file
 
@@ -29,94 +31,6 @@
 OnConflictOpts = Enum('OnConflictOpts', ['ERROR', 'OVERWRITE', 'SKIP'])
 
 
-def is_dicom_file_or_directory(path: PathLike) -> bool:
-    """
-    Check if a path is a DICOM file or directory containing DICOM files.
-
-    Args:
-        path: Path to check
-
-    Returns:
-        True if path is a DICOM file or directory with DICOM files, False otherwise
-    """
-    path = Path(path)
-
-    if path.is_file():
-        # check for 'DICM' magic bytes at offset 128
-        try:
-            with open(path, 'rb') as f:
-                f.seek(128)
-                return f.read(4) == b'DICM'
-        except OSError:
-            return False
-
-    elif path.is_dir():
-        for file in path.iterdir():
-            if file.is_file():
-                try:
-                    with open(file, 'rb') as f:
-                        f.seek(128)
-                        if f.read(4) == b'DICM':
-                            return True
-                except OSError:
-                    continue
-
-    return False
-
-
-def convert_dicom_to_nifti(input_path: PathLike, output_filepath: PathLike) -> None:
-    """
-    Convert DICOM file(s) to NIfTI format using pydicom and nibabel.
-
-    Args:
-        input_path: Path to either a DICOM file or directory containing DICOM files
-        output_filepath: Path where the output NIfTI file should be saved
-
-    Raises:
-        RuntimeError: If the conversion fails
-    """
-    input_path = Path(input_path)
-    output_filepath = Path(output_filepath)
-
-    try:
-        if input_path.is_file():
-            dicom_files = [input_path]
-        else:
-            # dicom files may not have .dcm extension
-            dicom_files = [f for f in input_path.iterdir() if f.is_file()]
-
-        if not dicom_files:
-            raise RuntimeError("No DICOM files found")
-
-        slices = []
-        for dcm_file in dicom_files:
-            try:
-                ds = pydicom.dcmread(dcm_file)
-                slices.append((ds.get('InstanceNumber', 0), ds.pixel_array))
-            except Exception:
-                # skip files that aren't valid dicom
-                continue
-
-        if not slices:
-            raise RuntimeError("No valid DICOM files found")
-
-        # sort by instance number - this is the slice order in the series
-        # so we reconstruct the 3D volume in the right order
-        slices.sort(key=lambda x: x[0])
-
-        # stack into 3D volume (handles both single and multiple slices)
-        volume = np.stack([s[1] for s in slices], axis=-1)
-
-        # identity affine for now - could extract from dicom headers in the future
-        affine = np.eye(4)
-
-        nifti_img = nib.Nifti1Image(volume, affine)
-        nib.save(nifti_img, str(output_filepath))
-
-    except Exception as e:
-        raise RuntimeError(f"DICOM to NIfTI conversion failed: {e}") from e
-
-
 class Database:
     def __init__(self, path: str | None = None):
         if path is None:
@@ -1096,7 +1010,6 @@ def load_volume(self, subject, volume_id):
         if not volume_data_filepath.exists() or not volume_data_filepath.is_file():
             self.logger.error(f"Volume data file not found for volume {volume_id}, subject {subject.id}")
             raise FileNotFoundError(f"Volume data file not found for volume {volume_id}, subject {subject.id}")
-        import nibabel as nib
         # Load the volume data using nibabel
         volume_data = nib.load(volume_data_filepath)
         self.logger.info(f"Loaded volume {volume_id} for subject {subject.id}")
diff --git a/src/openlifu/util/volume_conversion.py b/src/openlifu/util/volume_conversion.py
@@ -0,0 +1,162 @@
+"""Utilities for converting between different medical imaging formats."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+import nibabel as nib
+import numpy as np
+import pydicom
+
+from openlifu.util.types import PathLike
+
+
+def is_dicom_file_or_directory(path: PathLike) -> bool:
+    """
+    Check if a path is a DICOM file or directory containing DICOM files.
+
+    Args:
+        path: Path to check
+
+    Returns:
+        True if path is a DICOM file or directory with DICOM files, False otherwise
+    """
+    path = Path(path)
+
+    if path.is_file():
+        # check for 'DICM' magic bytes at offset 128
+        try:
+            with open(path, 'rb') as f:
+                f.seek(128)
+                return f.read(4) == b'DICM'
+        except OSError:
+            return False
+
+    elif path.is_dir():
+        for file in path.iterdir():
+            if file.is_file():
+                try:
+                    with open(file, 'rb') as f:
+                        f.seek(128)
+                        if f.read(4) == b'DICM':
+                            return True
+                except OSError:
+                    continue
+
+    return False
+
+
+def extract_affine_from_dicom(dicom_slices: list) -> np.ndarray:
+    """
+    Extract the affine transformation matrix from DICOM header information.
+
+    Args:
+        dicom_slices: List of tuples (instance_number, pixel_array, dicom_dataset)
+
+    Returns:
+        4x4 affine transformation matrix mapping voxel coordinates to world coordinates
+
+    Raises:
+        RuntimeError: If required DICOM tags are missing
+    """
+    # use the first slice to extract most parameters
+    first_ds = dicom_slices[0][2]
+
+    try:
+        # ImageOrientationPatient (0020,0037): direction cosines for row and column
+        orientation = np.array(first_ds.ImageOrientationPatient, dtype=float)
+        row_cosine = orientation[:3]  # direction cosines for rows
+        col_cosine = orientation[3:]  # direction cosines for columns
+
+        # ImagePositionPatient (0020,0032): position of the upper-left voxel
+        position = np.array(first_ds.ImagePositionPatient, dtype=float)
+
+        # PixelSpacing (0028,0030): spacing between pixels [row_spacing, col_spacing]
+        pixel_spacing = np.array(first_ds.PixelSpacing, dtype=float)
+        row_spacing = pixel_spacing[0]
+        col_spacing = pixel_spacing[1]
+
+    except AttributeError as e:
+        raise RuntimeError(
+            f"Missing required DICOM tag for affine calculation: {e}"
+        ) from e
+
+    # calculate slice direction as cross product of row and column directions
+    slice_cosine = np.cross(row_cosine, col_cosine)
+
+    # calculate slice spacing
+    if len(dicom_slices) > 1:
+        # calculate from the distance between first two slices
+        first_pos = np.array(dicom_slices[0][2].ImagePositionPatient, dtype=float)
+        second_pos = np.array(dicom_slices[1][2].ImagePositionPatient, dtype=float)
+        slice_spacing = np.linalg.norm(second_pos - first_pos)
+    else:
+        # single slice - try to get from SliceThickness or default to 1.0
+        slice_spacing = float(getattr(first_ds, 'SliceThickness', 1.0))
+
+    # construct the affine matrix
+    # The affine maps from voxel indices (i, j, k) to physical coordinates (x, y, z)
+    affine = np.eye(4)
+    affine[:3, 0] = row_cosine * row_spacing
+    affine[:3, 1] = col_cosine * col_spacing
+    affine[:3, 2] = slice_cosine * slice_spacing
+    affine[:3, 3] = position
+
+    return affine
+
+
+def convert_dicom_to_nifti(input_path: PathLike, output_filepath: PathLike) -> None:
+    """
+    Convert DICOM file(s) to NIfTI format using pydicom and nibabel.
+
+    The affine transformation matrix is extracted from DICOM headers to properly
+    map voxel coordinates to scanner coordinates in the NIfTI output.
+
+    Args:
+        input_path: Path to either a DICOM file or directory containing DICOM files
+        output_filepath: Path where the output NIfTI file should be saved
+
+    Raises:
+        RuntimeError: If the conversion fails
+    """
+    input_path = Path(input_path)
+    output_filepath = Path(output_filepath)
+
+    try:
+        if input_path.is_file():
+            dicom_files = [input_path]
+        else:
+            # dicom files may not have .dcm extension
+            dicom_files = [f for f in input_path.iterdir() if f.is_file()]
+
+        if not dicom_files:
+            raise RuntimeError("No DICOM files found")
+
+        slices = []
+        for dcm_file in dicom_files:
+            try:
+                ds = pydicom.dcmread(dcm_file)
+                # store instance number, pixel array, and dataset for affine extraction
+                slices.append((ds.get('InstanceNumber', 0), ds.pixel_array, ds))
+            except Exception:
+                # skip files that aren't valid dicom
+                continue
+
+        if not slices:
+            raise RuntimeError("No valid DICOM files found")
+
+        # sort by instance number - this is the slice order in the series
+        # so we reconstruct the 3D volume in the right order
+        slices.sort(key=lambda x: x[0])
+
+        # stack into 3D volume (handles both single and multiple slices)
+        volume = np.stack([s[1] for s in slices], axis=-1)
+
+        # extract affine from DICOM headers
+        affine = extract_affine_from_dicom(slices)
+
+        nifti_img = nib.Nifti1Image(volume, affine)
+        nib.save(nifti_img, str(output_filepath))
+
+    except Exception as e:
+        raise RuntimeError(f"DICOM to NIfTI conversion failed: {e}") from e
diff --git a/tests/test_database.py b/tests/test_database.py
@@ -15,11 +15,12 @@
 
 from openlifu import Solution
 from openlifu.db import Session, Subject, User
-from openlifu.db.database import Database, OnConflictOpts, is_dicom_file_or_directory
+from openlifu.db.database import Database, OnConflictOpts
 from openlifu.db.session import TransducerTrackingResult
 from openlifu.geo import ArrayTransform, Point
 from openlifu.nav.photoscan import Photoscan
 from openlifu.plan import Protocol, Run
+from openlifu.util.volume_conversion import is_dicom_file_or_directory
 from openlifu.xdc import Transducer
 
 
@@ -744,27 +745,6 @@ def test_get_transducer_absolute_filepaths(example_database, tmp_path: Path, reg
     else:
         assert absolute_file_paths["transducer_body_abspath"] is None
 
-def test_is_dicom_file_or_directory(tmp_path: Path):
-    test_dicom_file = Path(__file__).parent / "resources" / "CT_small.dcm"
-    assert test_dicom_file.exists(), "CT_small.dcm test file should exist"
-    assert is_dicom_file_or_directory(test_dicom_file)
-
-    non_dicom_file = tmp_path / "test.nii"
-    non_dicom_file.write_bytes(b'not a dicom file')
-    assert not is_dicom_file_or_directory(non_dicom_file)
-
-    dicom_dir = tmp_path / "dicom_series"
-    dicom_dir.mkdir()
-    shutil.copy(test_dicom_file, dicom_dir / "CT_small.dcm")
-    assert is_dicom_file_or_directory(dicom_dir)
-
-    non_dicom_dir = tmp_path / "non_dicom"
-    non_dicom_dir.mkdir()
-    (non_dicom_dir / "file.txt").write_text("not dicom")
-    assert not is_dicom_file_or_directory(non_dicom_dir)
-
-    assert not is_dicom_file_or_directory(tmp_path / "nonexistent.dcm")
-
 def test_write_volume_dicom(example_database: Database):
     """Test writing a volume from DICOM file - conversion to NIfTI and storage"""
     subject_id = "example_subject"
diff --git a/tests/test_volume_conversion.py b/tests/test_volume_conversion.py
