Add LPS -> RAS conversion for affine in dicom -> nifti pipeline (#396)

Author: arhowe00

src/openlifu/util/volume_conversion.py

@@ -49,12 +49,13 @@ def is_dicom_file_or_directory(path: PathLike) -> bool:
 def extract_affine_from_dicom(dicom_slices: list) -> np.ndarray:
     """
     Extract the affine transformation matrix from DICOM header information.
+    Converts from DICOM LPS (Left-Posterior-Superior) to NIfTI RAS.
 
     Args:
         dicom_slices: List of tuples (instance_number, pixel_array, dicom_dataset)
 
     Returns:
-        4x4 affine transformation matrix mapping voxel coordinates to world coordinates
+        4x4 affine transformation matrix mapping voxel coordinates to RAS world coordinates
 
     Raises:
         RuntimeError: If required DICOM tags are missing
@@ -71,47 +72,42 @@ def extract_affine_from_dicom(dicom_slices: list) -> np.ndarray:
         # 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]
+        # PixelSpacing is [row_spacing, col_spacing], so map to dy, dx
+        dy, dx = np.array(first_ds.PixelSpacing, dtype=float)
 
     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
+    # Compute Z direction and spacing (handling potential gantry tilt)
     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)
+        dz = np.dot(second_pos - first_pos, slice_cosine)
     else:
         # single slice - try to get from SliceThickness or default to 1.0
-        slice_spacing = float(getattr(first_ds, 'SliceThickness', 1.0))
+        dz = 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)
+    # Construct affine in DICOM LPS space
     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, 0] = row_cosine * dx
+    affine[:3, 1] = col_cosine * dy
+    affine[:3, 2] = slice_cosine * dz
     affine[:3, 3] = position
 
-    return affine
+    # Convert LPS to RAS by flipping X and Y axes
+    return np.diag([-1, -1, 1, 1]) @ 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
@@ -136,8 +132,8 @@ def convert_dicom_to_nifti(input_path: PathLike, output_filepath: PathLike) -> N
         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))
+                # Transpose to swap (Row, Col) -> (X, Y) for NIfTI
+                slices.append((ds.get('InstanceNumber', 0), ds.pixel_array.T, ds))
             except Exception:
                 # skip files that aren't valid dicom
                 continue
@@ -155,8 +151,7 @@ def convert_dicom_to_nifti(input_path: PathLike, output_filepath: PathLike) -> N
         # extract affine from DICOM headers
         affine = extract_affine_from_dicom(slices)
 
-        nifti_img = nib.Nifti1Image(volume, affine)
-        nib.save(nifti_img, str(output_filepath))
+        nib.save(nib.Nifti1Image(volume, affine), str(output_filepath))
 
     except Exception as e:
         raise RuntimeError(f"DICOM to NIfTI conversion failed: {e}") from e

tests/test_volume_conversion.py

@@ -45,8 +45,10 @@ def test_convert_dicom_to_nifti_single_file(tmp_path: Path):
 
     assert output_nifti.exists()
     img = nib.load(output_nifti)
-    assert img.shape[2] == 1  # single slice
-    assert not np.allclose(img.affine, np.eye(4))  # affine should not be identity
+    # Ensure 3D shape even for single slice
+    assert img.shape[2] == 1
+    # Affine should definitely not be identity for a valid medical image
+    assert not np.allclose(img.affine, np.eye(4))
 
 
 def test_convert_dicom_to_nifti_directory(tmp_path: Path):
@@ -57,24 +59,67 @@ def test_convert_dicom_to_nifti_directory(tmp_path: Path):
 
     assert output_nifti.exists()
     img = nib.load(output_nifti)
-    assert img.shape[2] == 2  # two slices in test series
-    assert not np.allclose(img.affine, np.eye(4))  # affine should not be identity
+    assert img.shape[2] == 2
+    assert not np.allclose(img.affine, np.eye(4))
 
 
-def test_extract_affine_from_dicom():
-    # mock DICOM dataset with minimal required tags
+def test_extract_affine_standard_orientation():
+    """
+    Test affine extraction for a standard axial slice.
+    Verifies proper mapping of Row/Col spacing and LPS->RAS conversion.
+    """
     mock_ds = Mock()
+    # Standard Axial: Row is X, Col is Y
     mock_ds.ImageOrientationPatient = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]
+    # Position in LPS (Left, Posterior, Superior)
     mock_ds.ImagePositionPatient = [10.0, 20.0, 30.0]
-    mock_ds.PixelSpacing = [0.5, 0.5]
+    # PixelSpacing: [RowSpacing (dy), ColSpacing (dx)]
+    mock_ds.PixelSpacing = [0.5, 0.8]
     mock_ds.SliceThickness = 2.0
 
     slices = [(1, None, mock_ds)]
     affine = extract_affine_from_dicom(slices)
 
-    assert affine.shape == (4, 4)
-    assert np.allclose(affine[3, :], [0, 0, 0, 1])  # homogeneous coordinate
-    assert np.allclose(affine[:3, 3], [10.0, 20.0, 30.0])  # position
+    # Expected RAS Matrix:
+    # 1. Spacing: X=0.8, Y=0.5, Z=2.0
+    # 2. Orientation: Standard Axial means X aligns with Right, Y with Anterior.
+    #    However, DICOM is LPS.
+    #    LPS X+ (Left) converts to RAS X- (Left).
+    #    LPS Y+ (Posterior) converts to RAS Y- (Posterior).
+    # 3. Origin: [10, 20, 30] LPS -> [-10, -20, 30] RAS
+
+    expected_affine = np.array([
+        [-0.8,  0.0,  0.0, -10.0],
+        [ 0.0, -0.5,  0.0, -20.0],
+        [ 0.0,  0.0,  2.0,  30.0],
+        [ 0.0,  0.0,  0.0,   1.0]
+    ])
+
+    np.testing.assert_allclose(affine, expected_affine)
+
+
+def test_extract_affine_multi_slice_spacing():
+    """Test that Z-spacing is calculated from slice positions, not SliceThickness."""
+    ds1 = Mock()
+    ds1.ImageOrientationPatient = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]
+    ds1.PixelSpacing = [1.0, 1.0]
+    ds1.ImagePositionPatient = [0.0, 0.0, 10.0]  # z = 10
+
+    ds2 = Mock()
+    ds2.ImagePositionPatient = [0.0, 0.0, 12.5]  # z = 12.5
+
+    # Create slice list with dummy pixel arrays
+    slices = [
+        (1, None, ds1),
+        (2, None, ds2)
+    ]
+
+    affine = extract_affine_from_dicom(slices)
+
+    # Z spacing should be 12.5 - 10.0 = 2.5
+    assert np.isclose(affine[2, 2], 2.5)
+    # Origin should be from the first slice, converted to RAS (Z is unchanged)
+    assert np.isclose(affine[2, 3], 10.0)
 
 
 def test_extract_affine_missing_tags():