Add multi-timepoint qSPECT workflow and improve GE DICOM quantification

pytheranostics/dicomtools/dicomtools.py

@@ -26,6 +26,7 @@
     "ge healthcare": "ge",
     "ge medical systems": "ge",
 }
+_GE_PIXEL_SCALE_TAG = (0x0011, 0x103B)
 
 
 class DicomModify:
@@ -116,7 +117,12 @@ def make_bqml_suv(
         Notes
         -----
         This method mutates ``self.ds`` in place. Call :meth:`save` to write the
-        modified dataset to disk.
+        modified dataset to disk. For GE data, the stored pixels are multiplied
+        by private DICOM tag ``(0011,103B) Pixel Scale`` when it is present.
+        The camera calibration factor must independently match the GE
+        reconstruction protocol, including the factor-of-four count scaling
+        caused by ``projections multiplication`` when resolution recovery is
+        enabled.
         """
         # Half-life is in seconds
 
@@ -166,6 +172,10 @@ def make_bqml_suv(
 
         # Get image in Bq/ml
         A = self.ds.pixel_array.astype(np.float64)
+        ge_pixel_scale = 1.0
+        if manufacturer == "ge":
+            ge_pixel_scale = _get_ge_pixel_scale(self.ds)
+            A *= ge_pixel_scale
         A = A / (frame_duration * n_proj) * self.CF * 1e6 / vox_vol
 
         slope, intercept = dicom_slope_intercept(A)
@@ -250,6 +260,7 @@ def make_bqml_suv(
             "patient_id": [_require_dicom_text(self.ds, "PatientID")],
             "weight_kg": [weight],
             "height_cm": [height],
+            "ge_pixel_scale": [ge_pixel_scale],
             "pre_inj_activity_MBq": [pre_inj_activity],
             "pre_inj_datetime": [pre_inj_datetime],
             "post_inj_activity_MBq": [post_inj_activity],
@@ -283,11 +294,15 @@ def make_bqml_suv(
 
         # for storing as new series data
         sop_ins_uid = _require_dicom_text(self.ds, "SOPInstanceUID")
-        self.ds.SOPInstanceUID = _increment_uid_suffix(sop_ins_uid, "SOPInstanceUID")
+        new_sop_ins_uid = _generate_new_uid_from_source(sop_ins_uid, "SOPInstanceUID")
+        self.ds.SOPInstanceUID = new_sop_ins_uid
+        if getattr(self.ds, "file_meta", None) is not None:
+            self.ds.file_meta.MediaStorageSOPInstanceUID = new_sop_ins_uid
 
         ser_ins_uid = _require_dicom_text(self.ds, "SeriesInstanceUID")
-        prefix = _uid_prefix_for_generated_uid(ser_ins_uid, "SeriesInstanceUID")
-        self.ds.SeriesInstanceUID = generate_uid(prefix=prefix)
+        self.ds.SeriesInstanceUID = _generate_new_uid_from_source(
+            ser_ins_uid, "SeriesInstanceUID"
+        )
 
         # self.ds.MediaStorageSOPInstaceUID
         return inj_df
@@ -365,6 +380,67 @@ def _update_int16_pixel_metadata(
             del ds[keyword]
 
 
+def _get_ge_pixel_scale(ds: Dataset) -> float:
+    """Return the GE private Pixel Scale correction, defaulting to one.
+
+    GE may reduce reconstructed stored counts to remain within the 16-bit pixel
+    range. Private tag ``(0011,103B)`` records the multiplier required to
+    recover the quantitative count values.
+
+    Parameters
+    ----------
+    ds : pydicom.dataset.Dataset
+        GE DICOM dataset containing the optional private Pixel Scale tag.
+
+    Returns
+    -------
+    float
+        Positive multiplier to apply to the stored pixel array.
+
+    Raises
+    ------
+    ValueError
+        If the tag is present but is not a finite positive number.
+    """
+    element = ds.get(_GE_PIXEL_SCALE_TAG)
+    if element is None:
+        return 1.0
+
+    value = element.value
+    if isinstance(value, bytes):
+        try:
+            value = value.decode("ascii").strip(" \x00")
+        except UnicodeDecodeError as exc:
+            raise ValueError(
+                "GE Pixel Scale DICOM tag (0011,103B) must contain a positive "
+                f"numeric value; got non-ASCII bytes {value!r}."
+            ) from exc
+
+    try:
+        pixel_scale = float(value)
+    except (TypeError, ValueError) as exc:
+        raise ValueError(
+            "GE Pixel Scale DICOM tag (0011,103B) must contain a positive "
+            f"numeric value; got {value!r}."
+        ) from exc
+
+    if not np.isfinite(pixel_scale) or pixel_scale <= 0:
+        raise ValueError(
+            "GE Pixel Scale DICOM tag (0011,103B) must contain a finite positive "
+            f"value; got {value!r}."
+        )
+
+    if not np.isclose(pixel_scale, 1.0):
+        warnings.warn(
+            "Applying GE Pixel Scale from DICOM tag (0011,103B): "
+            f"stored counts will be multiplied by {pixel_scale}.",
+            RuntimeWarning,
+            stacklevel=2,
+        )
+
+    return pixel_scale
+
+
 def _require_dicom_value(
     ds: Dataset,
     keyword: str,
@@ -600,25 +676,12 @@ def _parse_dicom_date_time(date: str, time_value: str, context: str) -> datetime
     )
 
 
-def _increment_uid_suffix(uid: str, keyword: str) -> str:
-    """Increment the final numeric component of a DICOM UID."""
-    uid_parts = uid.split(".")
-    if (
-        not uid_parts
-        or "" in uid_parts
-        or not all(part.isdigit() for part in uid_parts)
-    ):
-        raise ValueError(
-            f"Required DICOM tag '{keyword}' must be a dot-separated numeric "
-            f"UID ending in a numeric component; got {uid!r}."
-        )
-
-    uid_parts[-1] = str(int(uid_parts[-1]) + 1)
-    return ".".join(uid_parts)
-
+def _generate_new_uid_from_source(uid: str, keyword: str) -> str:
+    """Generate a new UID, preserving the source root when it safely fits.
 
-def _uid_prefix_for_generated_uid(uid: str, keyword: str) -> str:
-    """Return a UID prefix suitable for pydicom.uid.generate_uid."""
+    If removing the final component still leaves a prefix longer than pydicom's
+    54-character limit, a UID using pydicom's default root is generated instead.
+    """
     uid_parts = uid.split(".")
     if len(uid_parts) < 2 or "" in uid_parts:
         raise ValueError(
@@ -633,12 +696,16 @@ def _uid_prefix_for_generated_uid(uid: str, keyword: str) -> str:
 
     prefix = ".".join(uid_parts[:-1]) + "."
     if len(prefix) > 54:
-        raise ValueError(
+        warnings.warn(
             f"UID prefix derived from DICOM tag '{keyword}' is too long for "
-            f"generate_uid; got {len(prefix)} characters."
+            f"generate_uid ({len(prefix)} characters); generating a new UID "
+            "with pydicom's default root.",
+            RuntimeWarning,
+            stacklevel=2,
         )
+        return generate_uid()
 
-    return prefix
+    return generate_uid(prefix=prefix)
 
 
 def _prepend_corrected_image_value(ds: Dataset, value: str) -> None:

tests/test_dicomtools.py

@@ -1,14 +1,61 @@
 import numpy as np
 import pytest
-from pydicom.dataset import Dataset
+from pydicom.dataset import Dataset, FileDataset, FileMetaDataset
+from pydicom.sequence import Sequence
+from pydicom.uid import UID, ExplicitVRLittleEndian
 
 from pytheranostics.dicomtools.dicomtools import (
+    DicomModify,
     _get_frame_duration_seconds,
+    _get_ge_pixel_scale,
     _normalize_dicom_manufacturer,
     _update_int16_pixel_metadata,
 )
 
 
+def _write_minimal_ge_spect(
+    path,
+    pixel_scale=None,
+    sop_instance_uid="1.2.826.0.1.3680043.8.498.1",
+    series_instance_uid="1.2.826.0.1.3680043.8.498.2",
+):
+    file_meta = FileMetaDataset()
+    file_meta.TransferSyntaxUID = ExplicitVRLittleEndian
+    file_meta.MediaStorageSOPInstanceUID = sop_instance_uid
+    dataset = FileDataset(path, {}, file_meta=file_meta, preamble=b"\0" * 128)
+    dataset.Manufacturer = "GE MEDICAL SYSTEMS"
+    dataset.Modality = "NM"
+    dataset.PatientID = "TEST001"
+    dataset.SeriesDate = "20220101"
+    dataset.SeriesTime = "090000"
+    dataset.AcquisitionTime = "090000"
+    dataset.SeriesDescription = "RECON_COUNTS"
+    dataset.SOPInstanceUID = sop_instance_uid
+    dataset.SeriesInstanceUID = series_instance_uid
+    dataset.Rows = 1
+    dataset.Columns = 1
+    dataset.NumberOfFrames = 1
+    dataset.SamplesPerPixel = 1
+    dataset.PhotometricInterpretation = "MONOCHROME2"
+    dataset.BitsAllocated = 16
+    dataset.BitsStored = 16
+    dataset.HighBit = 15
+    dataset.PixelRepresentation = 1
+    dataset.PixelSpacing = [10, 10]
+    dataset.SliceThickness = 10
+    dataset.CorrectedImage = ["ATTN"]
+    dataset.PixelData = np.array([100], dtype=np.int16).tobytes()
+
+    rotation = Dataset()
+    rotation.ActualFrameDuration = 1000
+    rotation.NumberOfFramesInRotation = 1
+    dataset.RotationInformationSequence = Sequence([rotation])
+    if pixel_scale is not None:
+        dataset.add_new((0x0011, 0x103B), "DS", str(pixel_scale))
+
+    dataset.save_as(path)
+
+
 @pytest.mark.parametrize(
     ("manufacturer", "expected"),
     [
@@ -92,3 +139,94 @@ def test_update_int16_pixel_metadata_removes_conventional_rescale_mapping():
     assert "RescaleType" not in dataset
     assert dataset.SmallestImagePixelValue == -2
     assert dataset.LargestImagePixelValue == 3
+
+
+def test_get_ge_pixel_scale_defaults_to_one_when_tag_is_absent():
+    assert _get_ge_pixel_scale(Dataset()) == 1.0
+
+
+@pytest.mark.parametrize(("vr", "value"), [("DS", "2.5"), ("UN", b"4.0\x00")])
+def test_get_ge_pixel_scale_reads_numeric_private_tag(vr, value):
+    dataset = Dataset()
+    dataset.add_new((0x0011, 0x103B), vr, value)
+
+    with pytest.warns(RuntimeWarning, match="Applying GE Pixel Scale"):
+        pixel_scale = _get_ge_pixel_scale(dataset)
+
+    assert pixel_scale == pytest.approx(
+        float(value.rstrip(b"\x00")) if isinstance(value, bytes) else float(value)
+    )
+
+
+@pytest.mark.parametrize("value", ["invalid", "0", "-2", "nan", b"\xff"])
+def test_get_ge_pixel_scale_rejects_invalid_values(value):
+    dataset = Dataset()
+    vr = "UN" if isinstance(value, bytes) else "LO"
+    dataset.add_new((0x0011, 0x103B), vr, value)
+
+    with pytest.raises(ValueError, match="GE Pixel Scale DICOM tag"):
+        _get_ge_pixel_scale(dataset)
+
+
+def test_make_bqml_suv_applies_ge_pixel_scale(tmp_path):
+    unscaled_path = tmp_path / "unscaled.dcm"
+    scaled_path = tmp_path / "scaled.dcm"
+    _write_minimal_ge_spect(unscaled_path)
+    _write_minimal_ge_spect(scaled_path, pixel_scale=4)
+
+    conversion_kwargs = {
+        "weight": 70,
+        "height": 170,
+        "injection_date": "20220101",
+        "pre_inj_activity": 1000,
+        "pre_inj_time": "0800",
+        "post_inj_activity": 10,
+        "post_inj_time": "0820",
+        "injection_time": "0810",
+    }
+    unscaled = DicomModify(str(unscaled_path), CF=1.0)
+    unscaled_summary = unscaled.make_bqml_suv(**conversion_kwargs)
+    scaled = DicomModify(str(scaled_path), CF=1.0)
+    with pytest.warns(RuntimeWarning, match="Applying GE Pixel Scale"):
+        scaled_summary = scaled.make_bqml_suv(**conversion_kwargs)
+
+    unscaled_slope = float(
+        unscaled.ds.RealWorldValueMappingSequence[0].RealWorldValueSlope
+    )
+    scaled_slope = float(scaled.ds.RealWorldValueMappingSequence[0].RealWorldValueSlope)
+    assert scaled_slope == pytest.approx(unscaled_slope * 4)
+    assert unscaled_summary.loc[0, "ge_pixel_scale"] == 1.0
+    assert scaled_summary.loc[0, "ge_pixel_scale"] == 4.0
+
+
+def test_make_bqml_suv_falls_back_for_long_uid_prefix(tmp_path):
+    source_path = tmp_path / "long_uid.dcm"
+    source_sop_uid = "2.16.840.1.114362.1.12164994.27025353676.690008536.501.7081"
+    source_series_uid = "2.16.840.1.114362.1.12164994.27025353676.690008536.501.7082"
+    _write_minimal_ge_spect(
+        source_path,
+        sop_instance_uid=source_sop_uid,
+        series_instance_uid=source_series_uid,
+    )
+
+    image = DicomModify(str(source_path), CF=1.0)
+    with pytest.warns(RuntimeWarning, match="too long for generate_uid"):
+        image.make_bqml_suv(
+            weight=70,
+            height=170,
+            injection_date="20220101",
+            pre_inj_activity=1000,
+            pre_inj_time="0800",
+            post_inj_activity=10,
+            post_inj_time="0820",
+            injection_time="0810",
+        )
+
+    new_sop_uid = str(image.ds.SOPInstanceUID)
+    new_series_uid = str(image.ds.SeriesInstanceUID)
+    assert UID(new_sop_uid).is_valid
+    assert UID(new_series_uid).is_valid
+    assert new_sop_uid != source_sop_uid
+    assert new_series_uid != source_series_uid
+    assert new_sop_uid != new_series_uid
+    assert str(image.ds.file_meta.MediaStorageSOPInstanceUID) == new_sop_uid

workflows/README.md

@@ -0,0 +1,20 @@
+# PyTheranostics workflows
+
+This directory contains task-oriented notebooks for running common processing
+steps. Unlike the notebooks under `docs/source/tutorials`, these workflows are
+intended to be configured and executed on local data rather than read as guided
+lessons.
+
+## Available workflows
+
+- `qSPECT/counts_to_bqml.ipynb`: convert reconstructed SPECT DICOM images from
+  scanner counts to quantitative Bq/mL DICOM files.
+
+## Usage conventions
+
+- Review and edit each notebook's configuration cells before execution.
+- Keep patient data, credentials, and generated outputs outside the repository.
+- Do not commit executed cell outputs containing patient information.
+- Move reusable processing logic into `pytheranostics`; notebooks should focus
+  on configuration, orchestration, review, and quality control.
+