From 8015e0f7f2350d7c6b82c587681035c900c594e8 Mon Sep 17 00:00:00 2001 From: "Pedro L. Esquinas" Date: Wed, 26 Aug 2026 16:39:22 -0700 Subject: [PATCH 1/8] Apply GE DICOM Pixel Scale during qSPECT conversion --- pytheranostics/dicomtools/dicomtools.py | 74 ++++++++++++++++++++++++- 1 file changed, 73 insertions(+), 1 deletion(-) diff --git a/pytheranostics/dicomtools/dicomtools.py b/pytheranostics/dicomtools/dicomtools.py index 1931c46..a82c7d2 100644 --- a/pytheranostics/dicomtools/dicomtools.py +++ b/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], @@ -365,6 +376,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, From a2acbdeef505460cf2760b2998c801796dec452a Mon Sep 17 00:00:00 2001 From: "Pedro L. Esquinas" Date: Wed, 26 Aug 2026 16:39:32 -0700 Subject: [PATCH 2/8] Add tests for GE Pixel Scale handling --- tests/test_dicomtools.py | 101 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 100 insertions(+), 1 deletion(-) diff --git a/tests/test_dicomtools.py b/tests/test_dicomtools.py index 65e4d84..b89f56d 100644 --- a/tests/test_dicomtools.py +++ b/tests/test_dicomtools.py @@ -1,14 +1,55 @@ 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 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): + file_meta = FileMetaDataset() + file_meta.TransferSyntaxUID = ExplicitVRLittleEndian + 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 = "1.2.826.0.1.3680043.8.498.1" + dataset.SeriesInstanceUID = "1.2.826.0.1.3680043.8.498.2" + 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 +133,61 @@ 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 From ae3d300a08929de6cf522a31f5aa123bffee1f18 Mon Sep 17 00:00:00 2001 From: "Pedro L. Esquinas" Date: Wed, 26 Aug 2026 16:39:41 -0700 Subject: [PATCH 3/8] Document task-oriented processing workflows --- workflows/README.md | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 workflows/README.md diff --git a/workflows/README.md b/workflows/README.md new file mode 100644 index 0000000..38fc8b6 --- /dev/null +++ b/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. + From a5b8980acde483ce00144f54f63946174a160591 Mon Sep 17 00:00:00 2001 From: "Pedro L. Esquinas" Date: Wed, 26 Aug 2026 16:39:51 -0700 Subject: [PATCH 4/8] Add multi-timepoint counts-to-Bq/mL qSPECT workflow --- workflows/qSPECT/counts_to_bqml.ipynb | 447 ++++++++++++++++++++++++++ 1 file changed, 447 insertions(+) create mode 100644 workflows/qSPECT/counts_to_bqml.ipynb diff --git a/workflows/qSPECT/counts_to_bqml.ipynb b/workflows/qSPECT/counts_to_bqml.ipynb new file mode 100644 index 0000000..3afeb52 --- /dev/null +++ b/workflows/qSPECT/counts_to_bqml.ipynb @@ -0,0 +1,447 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "workflow-introduction", + "metadata": {}, + "source": [ + "# Convert reconstructed SPECT counts to Bq/mL\n", + "\n", + "This workflow finds reconstructed SPECT DICOM files below a user-selected directory, converts each image from scanner counts to quantitative Bq/mL with `DicomModify.make_bqml_suv`, and writes a new file beside the source using the suffix `_qSPECT.dcm`.\n", + "\n", + "One run represents **one radiopharmaceutical administration**. All time points therefore share the calibration factor, injection record, patient weight and height. Run the notebook separately for data from different administrations. Source files are never modified." + ] + }, + { + "cell_type": "markdown", + "id": "imports-heading", + "metadata": {}, + "source": [ + "## 1. Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "imports", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pydicom\n", + "from IPython.display import display\n", + "\n", + "from pytheranostics.dicomtools.dicomtools import DicomModify\n", + "from pytheranostics.shared.radioactive_decay import get_activity_at_injection" + ] + }, + { + "cell_type": "markdown", + "id": "calibration-heading", + "metadata": {}, + "source": [ + "## 2. Configure the input and calibration\n", + "\n", + "Set `input_directory` to the directory containing the reconstructed SPECT time points. Subdirectories are searched recursively. The calibration factor must be validated for the site, camera, collimator, energy windows, isotope, and reconstruction protocol used for these images. Its required unit is **MBq/(counts/s)**.\n", + "\n", + "> **GE SPECT/CT warning:** GE reconstruction can enable **projections multiplication**. When it is enabled together with resolution recovery, reconstructed counts may be scaled by a factor of 4. The calibration factor entered here must have been established using the same reconstruction setting and must therefore account for this factor. Do not apply a separate manual factor of 4 in this workflow. If reconstructed counts are reduced to remain within the 16-bit integer range, `DicomModify.make_bqml_suv` automatically reads GE private DICOM tag `(0011,103B) Pixel Scale` and multiplies the stored pixels by that value before quantification." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "input-configuration", + "metadata": {}, + "outputs": [], + "source": [ + "input_directory = Path(\"/path/to/reconstructed_spect_timepoints\")\n", + "calibration_factor = 0.0 # MBq/(counts/s); replace with your validated value\n", + "\n", + "# Discovery and output behavior\n", + "search_recursively = True\n", + "allowed_modalities = {\"NM\"} # Add \"PT\" only if your SPECT data uses that modality.\n", + "overwrite_existing = False" + ] + }, + { + "cell_type": "markdown", + "id": "administration-heading", + "metadata": {}, + "source": [ + "## 3. Enter the administration and patient information\n", + "\n", + "Transcribe these values from the administration record. Activities are in MBq, weight is in kg, height is in cm, dates use `YYYYMMDD`, and times use 24-hour `HHMM`. The pre-injection measurement must be at or before injection, and the post-injection measurement must be at or after injection." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "administration-configuration", + "metadata": {}, + "outputs": [], + "source": [ + "weight_kg = 0.0\n", + "height_cm = 0.0\n", + "\n", + "injection_date = \"YYYYMMDD\"\n", + "pre_inj_activity_mbq = 0.0\n", + "pre_inj_time = \"HHMM\"\n", + "post_inj_activity_mbq = 0.0\n", + "post_inj_time = \"HHMM\"\n", + "injection_time = \"HHMM\"\n", + "\n", + "radiopharmaceutical = \"Lutetium-PSMA-617\"\n", + "half_life_seconds = 574300.0 # Lu-177\n", + "activity_meter_scale_factor = 1.0" + ] + }, + { + "cell_type": "markdown", + "id": "scanner-heading", + "metadata": {}, + "source": [ + "## 4. Configure scanner-specific options\n", + "\n", + "For Siemens data, `n_detectors` is multiplied by `NumberOfFramesInRotation` to recover the total projection count. Leave the frame-duration fallback as `None` unless the DICOM `ActualFrameDuration` is missing or invalid and a validated duration is known. The fallback is only applied to Siemens data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "scanner-configuration", + "metadata": {}, + "outputs": [], + "source": [ + "n_detectors = 2\n", + "siemens_frame_duration_fallback = None # seconds, for example 15.0" + ] + }, + { + "cell_type": "markdown", + "id": "discovery-heading", + "metadata": {}, + "source": [ + "## 5. Discover reconstructed SPECT DICOM files\n", + "\n", + "The discovery step reads headers only, selects the configured modalities, and excludes files already ending in `_qSPECT.dcm`. Review the resulting table carefully before conversion. In particular, confirm that every candidate is a reconstructed counts image produced with the calibration-matched protocol." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "discover-reconstructions", + "metadata": {}, + "outputs": [], + "source": [ + "def discover_spect_reconstructions(root, recursive=True):\n", + " root = Path(root).expanduser().resolve()\n", + " if not root.is_dir():\n", + " raise NotADirectoryError(f\"Input directory does not exist: {root}\")\n", + "\n", + " candidates = root.rglob(\"*\") if recursive else root.glob(\"*\")\n", + " records = []\n", + " for path in sorted(candidate for candidate in candidates if candidate.is_file()):\n", + " if path.name.casefold().endswith(\"_qspect.dcm\"):\n", + " continue\n", + " try:\n", + " ds = pydicom.dcmread(path, stop_before_pixels=True)\n", + " except Exception:\n", + " continue\n", + "\n", + " modality = str(getattr(ds, \"Modality\", \"\")).upper()\n", + " if modality not in allowed_modalities:\n", + " continue\n", + "\n", + " pixel_scale_element = ds.get((0x0011, 0x103B))\n", + " records.append(\n", + " {\n", + " \"path\": path,\n", + " \"patient_id\": str(getattr(ds, \"PatientID\", \"\")),\n", + " \"modality\": modality,\n", + " \"manufacturer\": str(getattr(ds, \"Manufacturer\", \"\")),\n", + " \"series_date\": str(getattr(ds, \"SeriesDate\", \"\")),\n", + " \"series_time\": str(getattr(ds, \"SeriesTime\", \"\")),\n", + " \"series_description\": str(getattr(ds, \"SeriesDescription\", \"\")),\n", + " \"units\": str(getattr(ds, \"Units\", \"\")),\n", + " \"ge_pixel_scale\": (\n", + " str(pixel_scale_element.value) if pixel_scale_element else \"\"\n", + " ),\n", + " }\n", + " )\n", + "\n", + " return pd.DataFrame.from_records(records)\n", + "\n", + "\n", + "spect_reconstructions = discover_spect_reconstructions(\n", + " input_directory, recursive=search_recursively\n", + ")\n", + "if spect_reconstructions.empty:\n", + " raise FileNotFoundError(\n", + " f\"No SPECT DICOM files with modalities {sorted(allowed_modalities)} \"\n", + " f\"were found below {input_directory}.\"\n", + " )\n", + "\n", + "display(spect_reconstructions)" + ] + }, + { + "cell_type": "markdown", + "id": "validation-heading", + "metadata": {}, + "source": [ + "## 6. Validate the shared configuration\n", + "\n", + "These checks catch unset placeholders before any output is written. They do not replace verification of the administration record or calibration protocol." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "validate-configuration", + "metadata": {}, + "outputs": [], + "source": [ + "if calibration_factor <= 0:\n", + " raise ValueError(\"Set calibration_factor to a positive validated value.\")\n", + "if weight_kg <= 0:\n", + " raise ValueError(\"Set weight_kg to a positive value.\")\n", + "if height_cm <= 0:\n", + " raise ValueError(\"Set height_cm to a positive value.\")\n", + "if len(injection_date) != 8 or not injection_date.isdigit():\n", + " raise ValueError(\"injection_date must use YYYYMMDD.\")\n", + "for name, value in {\n", + " \"pre_inj_time\": pre_inj_time,\n", + " \"post_inj_time\": post_inj_time,\n", + " \"injection_time\": injection_time,\n", + "}.items():\n", + " if len(value) != 4 or not value.isdigit():\n", + " raise ValueError(f\"{name} must use HHMM.\")\n", + "if pre_inj_activity_mbq <= 0 or post_inj_activity_mbq < 0:\n", + " raise ValueError(\"Pre-injection activity must be positive and post-injection activity cannot be negative.\")\n", + "if half_life_seconds <= 0:\n", + " raise ValueError(\"half_life_seconds must be positive.\")\n", + "patient_ids = {value for value in spect_reconstructions[\"patient_id\"] if value}\n", + "if len(patient_ids) != 1:\n", + " raise ValueError(\n", + " \"All discovered time points must belong to exactly one patient; \"\n", + " f\"found patient IDs: {sorted(patient_ids) or 'none'}.\"\n", + " )\n", + "\n", + "print(f\"Configuration validated for {len(spect_reconstructions)} time point(s).\")" + ] + }, + { + "cell_type": "markdown", + "id": "conversion-heading", + "metadata": {}, + "source": [ + "## 7. Convert all discovered time points\n", + "\n", + "Set `run_conversion = True` only after reviewing the candidate table and configuration. Existing outputs are skipped unless `overwrite_existing` is enabled. An error in one file is recorded in the summary without preventing the remaining time points from being attempted." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "convert-timepoints", + "metadata": {}, + "outputs": [], + "source": [ + "run_conversion = False\n", + "\n", + "if not run_conversion:\n", + " raise RuntimeError(\n", + " \"Review the candidate table and configuration, then set \"\n", + " \"run_conversion = True to create qSPECT files.\"\n", + " )\n", + "\n", + "results = []\n", + "for source_path in spect_reconstructions[\"path\"]:\n", + " output_path = source_path.with_name(f\"{source_path.stem}_qSPECT.dcm\")\n", + " if output_path.exists() and not overwrite_existing:\n", + " results.append(\n", + " {\n", + " \"source_path\": source_path,\n", + " \"output_path\": output_path,\n", + " \"status\": \"skipped\",\n", + " \"message\": \"Output already exists\",\n", + " }\n", + " )\n", + " continue\n", + "\n", + " try:\n", + " image = DicomModify(str(source_path), CF=calibration_factor)\n", + " injection_summary = image.make_bqml_suv(\n", + " weight=weight_kg,\n", + " height=height_cm,\n", + " injection_date=injection_date,\n", + " pre_inj_activity=pre_inj_activity_mbq,\n", + " pre_inj_time=pre_inj_time,\n", + " post_inj_activity=post_inj_activity_mbq,\n", + " post_inj_time=post_inj_time,\n", + " injection_time=injection_time,\n", + " activity_meter_scale_factor=activity_meter_scale_factor,\n", + " half_life=half_life_seconds,\n", + " radiopharmaceutical=radiopharmaceutical,\n", + " n_detectors=n_detectors,\n", + " siemens_frame_duration_fallback=siemens_frame_duration_fallback,\n", + " )\n", + " image.ds.save_as(output_path)\n", + "\n", + " result = injection_summary.iloc[0].to_dict()\n", + " result.update(\n", + " {\n", + " \"source_path\": source_path,\n", + " \"output_path\": output_path,\n", + " \"status\": \"converted\",\n", + " \"message\": \"\",\n", + " }\n", + " )\n", + " results.append(result)\n", + " except Exception as exc:\n", + " results.append(\n", + " {\n", + " \"source_path\": source_path,\n", + " \"output_path\": output_path,\n", + " \"status\": \"failed\",\n", + " \"message\": f\"{type(exc).__name__}: {exc}\",\n", + " }\n", + " )\n", + "\n", + "conversion_summary = pd.DataFrame.from_records(results)\n", + "display(conversion_summary)" + ] + }, + { + "cell_type": "markdown", + "id": "fov-activity-heading", + "metadata": {}, + "source": [ + "## 8. Review total activity in the field of view\n", + "\n", + "As a quantitative sanity check, integrate the calibrated Bq/mL values over the entire image volume and compare the result with the activity administered at injection. The total activity in the SPECT field of view would not normally be expected to exceed the injected activity. A value above 100% is flagged for review.\n", + "\n", + "This check is not a substitute for quantitative QC. Unexpected values can indicate an incorrect calibration factor, GE reconstruction scaling mismatch, Pixel Scale problem, acquisition-duration or projection-count error, voxel-geometry error, decay-correction mismatch, contamination, or an incorrect administration record. Conversely, activity below the injected value is expected because of physical decay, biological clearance, and activity outside the field of view." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "report-fov-activity", + "metadata": {}, + "outputs": [], + "source": [ + "def calculate_fov_activity_mbq(dicom_path):\n", + " ds = pydicom.dcmread(dicom_path)\n", + " mapping_sequence = getattr(ds, \"RealWorldValueMappingSequence\", None)\n", + " if not mapping_sequence:\n", + " raise ValueError(\n", + " f\"RealWorldValueMappingSequence is missing from {dicom_path}.\"\n", + " )\n", + "\n", + " mapping = mapping_sequence[0]\n", + " slope = float(mapping.RealWorldValueSlope)\n", + " intercept = float(mapping.RealWorldValueIntercept)\n", + " activity_concentration_bqml = (\n", + " ds.pixel_array.astype(np.float64) * slope + intercept\n", + " )\n", + "\n", + " pixel_spacing_mm = np.asarray(ds.PixelSpacing, dtype=float)\n", + " slice_thickness_mm = float(ds.SliceThickness)\n", + " voxel_volume_ml = float(\n", + " np.prod([pixel_spacing_mm[0], pixel_spacing_mm[1], slice_thickness_mm])\n", + " / 1000.0\n", + " )\n", + " return float(activity_concentration_bqml.sum() * voxel_volume_ml / 1e6)\n", + "\n", + "\n", + "_, injected_activity_mbq = get_activity_at_injection(\n", + " injection_date=injection_date,\n", + " pre_inj_activity=pre_inj_activity_mbq,\n", + " pre_inj_time=pre_inj_time,\n", + " post_inj_activity=post_inj_activity_mbq,\n", + " post_inj_time=post_inj_time,\n", + " injection_time=injection_time,\n", + " half_life=half_life_seconds,\n", + ")\n", + "injected_activity_mbq *= activity_meter_scale_factor\n", + "\n", + "fov_activity_records = []\n", + "successful_outputs = conversion_summary[\n", + " conversion_summary[\"status\"].isin([\"converted\", \"skipped\"])\n", + "]\n", + "for row in successful_outputs.itertuples(index=False):\n", + " fov_activity_mbq = calculate_fov_activity_mbq(row.output_path)\n", + " percent_injected = 100.0 * fov_activity_mbq / injected_activity_mbq\n", + " fov_activity_records.append(\n", + " {\n", + " \"source_path\": row.source_path,\n", + " \"output_path\": row.output_path,\n", + " \"fov_activity_MBq\": fov_activity_mbq,\n", + " \"injected_activity_MBq\": injected_activity_mbq,\n", + " \"percent_of_injected_activity\": percent_injected,\n", + " \"exceeds_injected_activity\": fov_activity_mbq > injected_activity_mbq,\n", + " }\n", + " )\n", + "\n", + "fov_activity_summary = pd.DataFrame.from_records(fov_activity_records)\n", + "display(fov_activity_summary)\n", + "\n", + "if fov_activity_summary.empty:\n", + " print(\"No converted qSPECT outputs were available for the FOV activity check.\")\n", + "elif fov_activity_summary[\"exceeds_injected_activity\"].any():\n", + " print(\n", + " \"WARNING: At least one qSPECT image contains more total FOV activity \"\n", + " \"than the calculated injected activity. Review the highlighted values \"\n", + " \"and quantitative conversion inputs before use.\"\n", + " )\n", + "else:\n", + " print(\"PASS: No qSPECT image exceeds the calculated injected activity.\")" + ] + }, + { + "cell_type": "markdown", + "id": "quality-control-heading", + "metadata": {}, + "source": [ + "## 9. Confirm completion and perform quality control\n", + "\n", + "The cell below fails if any candidate could not be converted. Before quantitative use, independently confirm the calibration protocol, acquisition timing, projection count, voxel geometry, decay correction, injected activity, quantitative scaling, and image appearance for every output." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "confirm-completion", + "metadata": {}, + "outputs": [], + "source": [ + "failed = conversion_summary[conversion_summary[\"status\"] == \"failed\"]\n", + "if not failed.empty:\n", + " raise RuntimeError(\n", + " f\"{len(failed)} qSPECT conversion(s) failed. Review conversion_summary.\"\n", + " )\n", + "\n", + "converted_count = int((conversion_summary[\"status\"] == \"converted\").sum())\n", + "skipped_count = int((conversion_summary[\"status\"] == \"skipped\").sum())\n", + "print(f\"Completed: {converted_count} converted, {skipped_count} skipped.\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".dosimetry (3.14.4.final.0)", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.14.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From c898e03bd1745ac77aacee97b859238862773e65 Mon Sep 17 00:00:00 2001 From: "Pedro L. Esquinas" Date: Thu, 27 Aug 2026 08:18:20 -0700 Subject: [PATCH 5/8] Add multi-timepoint counts-to-Bq/mL qSPECT workflow --- workflows/qSPECT/counts_to_bqml.ipynb | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/workflows/qSPECT/counts_to_bqml.ipynb b/workflows/qSPECT/counts_to_bqml.ipynb index 3afeb52..9d3f90d 100644 --- a/workflows/qSPECT/counts_to_bqml.ipynb +++ b/workflows/qSPECT/counts_to_bqml.ipynb @@ -22,7 +22,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "imports", "metadata": {}, "outputs": [], @@ -52,13 +52,13 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "input-configuration", "metadata": {}, "outputs": [], "source": [ - "input_directory = Path(\"/path/to/reconstructed_spect_timepoints\")\n", - "calibration_factor = 0.0 # MBq/(counts/s); replace with your validated value\n", + "input_directory = Path(\"/home/esquinas/projects/DosimetryPR21/CAGH-0005/Cycle01\")\n", + "calibration_factor = 0.235579647832098 # MBq/(counts/s); replace with your validated value\n", "\n", "# Discovery and output behavior\n", "search_recursively = True\n", @@ -438,7 +438,15 @@ "name": "python3" }, "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", "version": "3.14.4" } }, From 86f334d7ebd781e3d7f0b7946ed7a6b90e94abcc Mon Sep 17 00:00:00 2001 From: "Pedro L. Esquinas" Date: Thu, 27 Aug 2026 08:18:57 -0700 Subject: [PATCH 6/8] Revert "Add multi-timepoint counts-to-Bq/mL qSPECT workflow" This reverts commit c898e03bd1745ac77aacee97b859238862773e65. --- workflows/qSPECT/counts_to_bqml.ipynb | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/workflows/qSPECT/counts_to_bqml.ipynb b/workflows/qSPECT/counts_to_bqml.ipynb index 9d3f90d..3afeb52 100644 --- a/workflows/qSPECT/counts_to_bqml.ipynb +++ b/workflows/qSPECT/counts_to_bqml.ipynb @@ -22,7 +22,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "imports", "metadata": {}, "outputs": [], @@ -52,13 +52,13 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "input-configuration", "metadata": {}, "outputs": [], "source": [ - "input_directory = Path(\"/home/esquinas/projects/DosimetryPR21/CAGH-0005/Cycle01\")\n", - "calibration_factor = 0.235579647832098 # MBq/(counts/s); replace with your validated value\n", + "input_directory = Path(\"/path/to/reconstructed_spect_timepoints\")\n", + "calibration_factor = 0.0 # MBq/(counts/s); replace with your validated value\n", "\n", "# Discovery and output behavior\n", "search_recursively = True\n", @@ -438,15 +438,7 @@ "name": "python3" }, "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", "version": "3.14.4" } }, From bd533e7ec46e535b899f3c2b09876c0b9c88a693 Mon Sep 17 00:00:00 2001 From: "Pedro L. Esquinas" Date: Thu, 27 Aug 2026 09:17:04 -0700 Subject: [PATCH 7/8] Fix qSPECT UID generation for long DICOM prefixes --- pytheranostics/dicomtools/dicomtools.py | 43 ++++++++++------------ tests/test_dicomtools.py | 47 ++++++++++++++++++++++--- 2 files changed, 62 insertions(+), 28 deletions(-) diff --git a/pytheranostics/dicomtools/dicomtools.py b/pytheranostics/dicomtools/dicomtools.py index a82c7d2..fd98bc3 100644 --- a/pytheranostics/dicomtools/dicomtools.py +++ b/pytheranostics/dicomtools/dicomtools.py @@ -294,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 @@ -672,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( @@ -705,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: diff --git a/tests/test_dicomtools.py b/tests/test_dicomtools.py index b89f56d..00688fc 100644 --- a/tests/test_dicomtools.py +++ b/tests/test_dicomtools.py @@ -2,7 +2,7 @@ import pytest from pydicom.dataset import Dataset, FileDataset, FileMetaDataset from pydicom.sequence import Sequence -from pydicom.uid import ExplicitVRLittleEndian +from pydicom.uid import UID, ExplicitVRLittleEndian from pytheranostics.dicomtools.dicomtools import ( DicomModify, @@ -13,9 +13,15 @@ ) -def _write_minimal_ge_spect(path, pixel_scale=None): +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" @@ -24,8 +30,8 @@ def _write_minimal_ge_spect(path, pixel_scale=None): dataset.SeriesTime = "090000" dataset.AcquisitionTime = "090000" dataset.SeriesDescription = "RECON_COUNTS" - dataset.SOPInstanceUID = "1.2.826.0.1.3680043.8.498.1" - dataset.SeriesInstanceUID = "1.2.826.0.1.3680043.8.498.2" + dataset.SOPInstanceUID = sop_instance_uid + dataset.SeriesInstanceUID = series_instance_uid dataset.Rows = 1 dataset.Columns = 1 dataset.NumberOfFrames = 1 @@ -191,3 +197,36 @@ def test_make_bqml_suv_applies_ge_pixel_scale(tmp_path): 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 From 6a24ff89ca2fcd8bff01b5b7986b37d230d5f9a0 Mon Sep 17 00:00:00 2001 From: "Pedro L. Esquinas" Date: Thu, 27 Aug 2026 10:58:02 -0700 Subject: [PATCH 8/8] Update explanation for GE scaling and sanity checks on activity. --- workflows/qSPECT/counts_to_bqml.ipynb | 177 +++++++++++++++++++++++--- 1 file changed, 159 insertions(+), 18 deletions(-) diff --git a/workflows/qSPECT/counts_to_bqml.ipynb b/workflows/qSPECT/counts_to_bqml.ipynb index 3afeb52..db74ed4 100644 --- a/workflows/qSPECT/counts_to_bqml.ipynb +++ b/workflows/qSPECT/counts_to_bqml.ipynb @@ -29,6 +29,7 @@ "source": [ "from pathlib import Path\n", "\n", + "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import pydicom\n", @@ -47,7 +48,7 @@ "\n", "Set `input_directory` to the directory containing the reconstructed SPECT time points. Subdirectories are searched recursively. The calibration factor must be validated for the site, camera, collimator, energy windows, isotope, and reconstruction protocol used for these images. Its required unit is **MBq/(counts/s)**.\n", "\n", - "> **GE SPECT/CT warning:** GE reconstruction can enable **projections multiplication**. When it is enabled together with resolution recovery, reconstructed counts may be scaled by a factor of 4. The calibration factor entered here must have been established using the same reconstruction setting and must therefore account for this factor. Do not apply a separate manual factor of 4 in this workflow. If reconstructed counts are reduced to remain within the 16-bit integer range, `DicomModify.make_bqml_suv` automatically reads GE private DICOM tag `(0011,103B) Pixel Scale` and multiplies the stored pixels by that value before quantification." + "> **GE SPECT/CT warning:** GE **projections multiplication** and the private DICOM **Pixel Scale** are separate scaling mechanisms. When projections multiplication is enabled together with resolution recovery, reconstructed counts may be scaled by a factor of 4. Tag `(0011,103B) Pixel Scale`, on the other hand, records a multiplier used when GE reduces reconstructed values to fit the 16-bit stored-pixel range. `DicomModify.make_bqml_suv` automatically applies Pixel Scale when the tag is present; when it is absent, the workflow assumes a Pixel Scale of 1.0. Absence of this tag does **not** show that projections multiplication was disabled, and the x4 setting cannot be determined reliably from the reconstructed counts alone. Verify it from the saved reconstruction protocol, scanner console, protocol screenshots, or vendor metadata. If necessary, reconstruct the same phantom or projection data with the option enabled and disabled and compare integrated counts. The calibration factor must have been established with exactly the same reconstruction settings as the patient images and must already account for projections multiplication. Do not apply a separate manual factor of 4 in this workflow. The activity sanity check in block 8 may reveal a gross scaling mismatch, but it cannot prove whether projections multiplication was enabled because biological clearance and incomplete FOV coverage can conceal such a mismatch." ] }, { @@ -57,13 +58,13 @@ "metadata": {}, "outputs": [], "source": [ - "input_directory = Path(\"/path/to/reconstructed_spect_timepoints\")\n", + "input_directory = Path(\"/path/to/spect/data/\")\n", "calibration_factor = 0.0 # MBq/(counts/s); replace with your validated value\n", "\n", "# Discovery and output behavior\n", "search_recursively = True\n", "allowed_modalities = {\"NM\"} # Add \"PT\" only if your SPECT data uses that modality.\n", - "overwrite_existing = False" + "overwrite_existing = True" ] }, { @@ -83,13 +84,13 @@ "metadata": {}, "outputs": [], "source": [ - "weight_kg = 0.0\n", - "height_cm = 0.0\n", + "weight_kg = 0\n", + "height_cm = 0\n", "\n", "injection_date = \"YYYYMMDD\"\n", - "pre_inj_activity_mbq = 0.0\n", + "pre_inj_activity_mbq = 0\n", "pre_inj_time = \"HHMM\"\n", - "post_inj_activity_mbq = 0.0\n", + "post_inj_activity_mbq = 0\n", "post_inj_time = \"HHMM\"\n", "injection_time = \"HHMM\"\n", "\n", @@ -250,7 +251,7 @@ "metadata": {}, "outputs": [], "source": [ - "run_conversion = False\n", + "run_conversion = True\n", "\n", "if not run_conversion:\n", " raise RuntimeError(\n", @@ -322,7 +323,7 @@ "source": [ "## 8. Review total activity in the field of view\n", "\n", - "As a quantitative sanity check, integrate the calibrated Bq/mL values over the entire image volume and compare the result with the activity administered at injection. The total activity in the SPECT field of view would not normally be expected to exceed the injected activity. A value above 100% is flagged for review.\n", + "As a quantitative sanity check, integrate the calibrated Bq/mL values over the entire image volume. At each acquisition time, the measured FOV activity is compared with the net injected activity physically decayed to that time using the supplied half-life and a 10% tolerance. The earliest measured time point is also extrapolated back to injection time for visualization; no fit is performed.\n", "\n", "This check is not a substitute for quantitative QC. Unexpected values can indicate an incorrect calibration factor, GE reconstruction scaling mismatch, Pixel Scale problem, acquisition-duration or projection-count error, voxel-geometry error, decay-correction mismatch, contamination, or an incorrect administration record. Conversely, activity below the injected value is expected because of physical decay, biological clearance, and activity outside the field of view." ] @@ -334,6 +335,27 @@ "metadata": {}, "outputs": [], "source": [ + "def get_acquisition_datetime(ds):\n", + " acquisition_date = str(\n", + " getattr(ds, \"AcquisitionDate\", None)\n", + " or getattr(ds, \"SeriesDate\", None)\n", + " or getattr(ds, \"StudyDate\", \"\")\n", + " )\n", + " acquisition_time = str(\n", + " getattr(ds, \"AcquisitionTime\", None)\n", + " or getattr(ds, \"SeriesTime\", None)\n", + " or getattr(ds, \"StudyTime\", \"\")\n", + " )\n", + " if not acquisition_date or not acquisition_time:\n", + " raise ValueError(\"The DICOM acquisition date or time is missing.\")\n", + "\n", + " # DICOM TM may contain fractional seconds; pad shorter HH or HHMM values.\n", + " acquisition_time = acquisition_time.split(\".\")[0].ljust(6, \"0\")[:6]\n", + " return pd.to_datetime(\n", + " acquisition_date + acquisition_time, format=\"%Y%m%d%H%M%S\"\n", + " )\n", + "\n", + "\n", "def calculate_fov_activity_mbq(dicom_path):\n", " ds = pydicom.dcmread(dicom_path)\n", " mapping_sequence = getattr(ds, \"RealWorldValueMappingSequence\", None)\n", @@ -369,37 +391,140 @@ ")\n", "injected_activity_mbq *= activity_meter_scale_factor\n", "\n", + "injection_datetime = pd.to_datetime(\n", + " injection_date + injection_time, format=\"%Y%m%d%H%M\"\n", + ")\n", + "\n", + "activity_tolerance_fraction = 0.10\n", + "half_life_hours = half_life_seconds / 3600.0\n", + "decay_constant_per_hour = np.log(2.0) / half_life_hours\n", + "\n", "fov_activity_records = []\n", "successful_outputs = conversion_summary[\n", " conversion_summary[\"status\"].isin([\"converted\", \"skipped\"])\n", "]\n", "for row in successful_outputs.itertuples(index=False):\n", + " ds = pydicom.dcmread(row.output_path, stop_before_pixels=True)\n", + " acquisition_datetime = get_acquisition_datetime(ds)\n", + " time_hours = (acquisition_datetime - injection_datetime).total_seconds() / 3600.0\n", " fov_activity_mbq = calculate_fov_activity_mbq(row.output_path)\n", - " percent_injected = 100.0 * fov_activity_mbq / injected_activity_mbq\n", + " decay_corrected_injected_activity_mbq = injected_activity_mbq * np.exp(\n", + " -decay_constant_per_hour * time_hours\n", + " )\n", + " activity_tolerance_limit_mbq = decay_corrected_injected_activity_mbq * (\n", + " 1.0 + activity_tolerance_fraction\n", + " )\n", + " percent_of_decay_corrected_injected = (\n", + " 100.0 * fov_activity_mbq / decay_corrected_injected_activity_mbq\n", + " )\n", " fov_activity_records.append(\n", " {\n", " \"source_path\": row.source_path,\n", " \"output_path\": row.output_path,\n", + " \"acquisition_datetime\": acquisition_datetime,\n", + " \"time_hours\": time_hours,\n", " \"fov_activity_MBq\": fov_activity_mbq,\n", - " \"injected_activity_MBq\": injected_activity_mbq,\n", - " \"percent_of_injected_activity\": percent_injected,\n", - " \"exceeds_injected_activity\": fov_activity_mbq > injected_activity_mbq,\n", + " \"decay_corrected_injected_activity_MBq\": decay_corrected_injected_activity_mbq,\n", + " \"injected_plus_10_percent_MBq\": activity_tolerance_limit_mbq,\n", + " \"percent_of_decay_corrected_injected_activity\": percent_of_decay_corrected_injected,\n", + " \"within_10_percent_tolerance\": fov_activity_mbq <= activity_tolerance_limit_mbq,\n", " }\n", " )\n", "\n", "fov_activity_summary = pd.DataFrame.from_records(fov_activity_records)\n", + "if not fov_activity_summary.empty:\n", + " fov_activity_summary = fov_activity_summary.sort_values(\"time_hours\").reset_index(drop=True)\n", "display(fov_activity_summary)\n", "\n", + "if not fov_activity_summary.empty:\n", + " if (fov_activity_summary[\"time_hours\"] < 0).any():\n", + " raise ValueError(\"At least one acquisition occurred before injection.\")\n", + "\n", + " first_time_hours = float(fov_activity_summary.iloc[0][\"time_hours\"])\n", + " first_activity_mbq = float(fov_activity_summary.iloc[0][\"fov_activity_MBq\"])\n", + "\n", + " # Back-extrapolate only the earliest measurement; no curve fitting is used.\n", + " extrapolated_activity_t0_mbq = first_activity_mbq * np.exp(\n", + " decay_constant_per_hour * first_time_hours\n", + " )\n", + " activity_t0_tolerance_limit_mbq = (\n", + " injected_activity_mbq * (1.0 + activity_tolerance_fraction)\n", + " )\n", + "\n", + " extrapolation_summary = pd.DataFrame(\n", + " {\n", + " \"extrapolated_activity_t0_MBq\": [extrapolated_activity_t0_mbq],\n", + " \"injected_activity_MBq\": [injected_activity_mbq],\n", + " \"tolerance_percent\": [100.0 * activity_tolerance_fraction],\n", + " \"injected_plus_tolerance_MBq\": [activity_t0_tolerance_limit_mbq],\n", + " \"within_tolerance\": [\n", + " extrapolated_activity_t0_mbq <= activity_t0_tolerance_limit_mbq\n", + " ],\n", + " }\n", + " )\n", + " display(extrapolation_summary)\n", + "\n", + " plot_time_hours = np.linspace(0.0, first_time_hours, 250)\n", + " fitted_activity_mbq = extrapolated_activity_t0_mbq * np.exp(\n", + " -decay_constant_per_hour * plot_time_hours\n", + " )\n", + " injected_curve_time_hours = np.linspace(\n", + " 0.0, fov_activity_summary[\"time_hours\"].max(), 250\n", + " )\n", + " injected_tolerance_curve_mbq = (\n", + " injected_activity_mbq\n", + " * (1.0 + activity_tolerance_fraction)\n", + " * np.exp(-decay_constant_per_hour * injected_curve_time_hours)\n", + " )\n", + "\n", + " fig, ax = plt.subplots(figsize=(8, 5))\n", + " ax.scatter(\n", + " fov_activity_summary[\"time_hours\"],\n", + " fov_activity_summary[\"fov_activity_MBq\"],\n", + " label=\"Measured FOV activity\",\n", + " zorder=3,\n", + " )\n", + " ax.plot(\n", + " plot_time_hours,\n", + " fitted_activity_mbq,\n", + " label=\"Physical decay from first time point\",\n", + " )\n", + " ax.scatter(\n", + " [0.0],\n", + " [extrapolated_activity_t0_mbq],\n", + " marker=\"D\",\n", + " s=60,\n", + " label=f\"Extrapolated t=0: {extrapolated_activity_t0_mbq:.1f} MBq\",\n", + " zorder=4,\n", + " )\n", + " ax.plot(\n", + " injected_curve_time_hours,\n", + " injected_tolerance_curve_mbq,\n", + " color=\"tab:red\",\n", + " linestyle=\"--\",\n", + " label=\"Decay-corrected injected activity +10%\",\n", + " )\n", + " ax.set_xlabel(\"Time after injection (hours)\")\n", + " ax.set_ylabel(\"Activity in field of view (MBq)\")\n", + " ax.set_title(\"qSPECT activity and back-extrapolation from first time point\")\n", + " ax.grid(True, alpha=0.3)\n", + " ax.legend()\n", + " fig.tight_layout()\n", + " plt.show()\n", + "\n", "if fov_activity_summary.empty:\n", " print(\"No converted qSPECT outputs were available for the FOV activity check.\")\n", - "elif fov_activity_summary[\"exceeds_injected_activity\"].any():\n", + "elif not fov_activity_summary[\"within_10_percent_tolerance\"].all():\n", " print(\n", - " \"WARNING: At least one qSPECT image contains more total FOV activity \"\n", - " \"than the calculated injected activity. Review the highlighted values \"\n", - " \"and quantitative conversion inputs before use.\"\n", + " \"WARNING: At least one measured FOV activity exceeds the net injected \"\n", + " \"activity physically decayed to that time point plus 10%. Review the \"\n", + " \"time-point results before quantitative use.\"\n", " )\n", "else:\n", - " print(\"PASS: No qSPECT image exceeds the calculated injected activity.\")" + " print(\n", + " \"PASS: Every measured FOV activity is within the net injected activity \"\n", + " \"physically decayed to its acquisition time plus 10%.\"\n", + " )" ] }, { @@ -429,6 +554,14 @@ "skipped_count = int((conversion_summary[\"status\"] == \"skipped\").sum())\n", "print(f\"Completed: {converted_count} converted, {skipped_count} skipped.\")" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45005496", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -438,7 +571,15 @@ "name": "python3" }, "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", "version": "3.14.4" } },