diff --git a/docs/source/examples/data/016/test016.dcm b/docs/source/examples/data/016/test016.dcm deleted file mode 100644 index 0eb6d78..0000000 Binary files a/docs/source/examples/data/016/test016.dcm and /dev/null differ diff --git a/docs/source/examples/data/test.dcm b/docs/source/examples/data/test.dcm deleted file mode 100644 index cd9948f..0000000 Binary files a/docs/source/examples/data/test.dcm and /dev/null differ diff --git a/docs/source/examples/data/test0034_2.dcm b/docs/source/examples/data/test0034_2.dcm deleted file mode 100644 index 68a9b1a..0000000 Binary files a/docs/source/examples/data/test0034_2.dcm and /dev/null differ diff --git a/docs/source/examples/data/test016.dcm b/docs/source/examples/data/test016.dcm deleted file mode 100644 index 6d5f808..0000000 Binary files a/docs/source/examples/data/test016.dcm and /dev/null differ diff --git a/docs/source/examples/data/testimages/0034.dcm b/docs/source/examples/data/testimages/0034.dcm deleted file mode 100644 index 0ac8dec..0000000 Binary files a/docs/source/examples/data/testimages/0034.dcm and /dev/null differ diff --git a/docs/source/examples/data/testimages/spect_counts.dcm b/docs/source/examples/data/testimages/spect_counts.dcm deleted file mode 100644 index 21ad36d..0000000 Binary files a/docs/source/examples/data/testimages/spect_counts.dcm and /dev/null differ diff --git a/docs/source/examples/data/testimages/spect_counts_out.dcm b/docs/source/examples/data/testimages/spect_counts_out.dcm deleted file mode 100644 index 05896e5..0000000 Binary files a/docs/source/examples/data/testimages/spect_counts_out.dcm and /dev/null differ diff --git a/docs/source/tutorials/SPECT2SUV/SPECT2SUV.ipynb b/docs/source/tutorials/SPECT2SUV/SPECT2SUV.ipynb deleted file mode 100644 index 478a918..0000000 --- a/docs/source/tutorials/SPECT2SUV/SPECT2SUV.ipynb +++ /dev/null @@ -1,183 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "6d072b76-26c2-4789-9d4a-7be2f0170f5c", - "metadata": {}, - "outputs": [], - "source": [ - "from doodle.dicomtools.dicomtools import DicomModify\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "from pydicom.dataset import Dataset\n", - "from pydicom.uid import generate_uid\n", - "\n", - "%load_ext autoreload\n", - "%autoreload 2\n", - "%matplotlib inline" - ] - }, - { - "cell_type": "markdown", - "id": "0c6f9dbd-0afa-4f6a-87ec-a81120df22ec", - "metadata": {}, - "source": [ - "\n", - "### Point to the SPECT image in counts (the one that you want to make quantitative)\n", - "\n", - "### and set the output path (the location where you want the QSPECT image to be saved at the end)" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "154c4542-32ab-4da3-89d1-170e10a539db", - "metadata": {}, - "outputs": [], - "source": [ - "from pathlib import Path\n", - "\n", - "def _find_examples_dir() -> Path:\n", - " candidates = [\n", - " Path.cwd() / \"examples\" / \"data\",\n", - " Path.cwd() / \"docs\" / \"source\" / \"examples\" / \"data\",\n", - " ]\n", - " for candidate in candidates:\n", - " if candidate.exists():\n", - " return candidate\n", - " raise FileNotFoundError(\"Could not locate docs example data directory\")\n", - "\n", - "EXAMPLES_DIR = _find_examples_dir()\n", - "spect_counts = str(EXAMPLES_DIR / \"testimages\" / \"016.dcm\")\n", - "output_path = str(Path.cwd() / \"test016.dcm\")\n" - ] - }, - { - "cell_type": "markdown", - "id": "82c8ff82-32de-43f5-936a-f63597c35bc5", - "metadata": {}, - "source": [ - "### Set the calibration factor for the camera of the centre that you're using" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "14c487d1-3d64-4794-a7e2-02e0ec98ba67", - "metadata": {}, - "outputs": [], - "source": [ - "CF = 0.10800584442987242\n", - "img=DicomModify(spect_counts,CF)" - ] - }, - { - "cell_type": "markdown", - "id": "15ff709a-66de-45be-960c-c7b0ab1da03d", - "metadata": {}, - "source": [ - "### Specify the following information from the filled in form from the injection" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "id": "3a7b4268-5548-4100-8173-95eae051ec75", - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
injected_activity_MBqinjection_datetime
07395.7448952022-06-16 09:18:00
\n", - "
" - ], - "text/plain": [ - " injected_activity_MBq injection_datetime\n", - "0 7395.744895 2022-06-16 09:18:00" - ] - }, - "execution_count": 12, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "weight = 113.4\n", - "height = 1.785\n", - "injection_date = '20220616'\n", - "pre_inj_activity = 7450\n", - "pre_inj_time = '0804'\n", - "post_inj_activity = 14.4\n", - "post_inj_time = '0955'\n", - "injection_time = '0918'\n", - "\n", - "#The activity meteer scale factor is a factor to multiply activity values if the calibration setting of the dose calibrator has changed\n", - "activity_meter_scale_factor = 1\n", - "\n", - "\n", - "inj_df = img.make_bqml_suv(weight=weight,height=height,injection_date=injection_date,pre_inj_activity=pre_inj_activity,pre_inj_time=pre_inj_time,post_inj_activity=post_inj_activity,post_inj_time=post_inj_time,injection_time=injection_time,activity_meter_scale_factor=activity_meter_scale_factor)\n", - "img.ds.save_as(output_path)\n", - "inj_df" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c4cbee90-0ee5-4915-8a4b-fd875a69d804", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "DOODLE", - "language": "python", - "name": "doodle" - }, - "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.9.7" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} \ No newline at end of file diff --git a/docs/source/tutorials/index.rst b/docs/source/tutorials/index.rst index e55dd70..165954e 100644 --- a/docs/source/tutorials/index.rst +++ b/docs/source/tutorials/index.rst @@ -10,5 +10,5 @@ Hands-on walkthroughs that demonstrate common PyTheranostics workflows. dicom_organization/dicom_organization Data_Ingestion_Examples/Data_Ingestion_Examples segmentation/total_segmentator_tutorial - SPECT2SUV/SPECT2SUV + qSPECT/qSPECT_Tutorial ROI_Mapping_Tutorial/ROI_Mapping_Tutorial diff --git a/docs/source/tutorials/qSPECT/qSPECT_Tutorial.ipynb b/docs/source/tutorials/qSPECT/qSPECT_Tutorial.ipynb new file mode 100644 index 0000000..c47967f --- /dev/null +++ b/docs/source/tutorials/qSPECT/qSPECT_Tutorial.ipynb @@ -0,0 +1,320 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "title", + "metadata": {}, + "source": [ + "# qSPECT Tutorial\n", + "\n", + "This tutorial converts a reconstructed SPECT DICOM image from scanner counts into a quantitative SPECT (qSPECT) image in Bq/mL using `DicomModify.make_bqml_suv`. It also adds the patient and radiopharmaceutical metadata needed by software that calculates standardized uptake values (SUVs).\n", + "\n", + "> **Before you begin:** the input must be a reconstructed SPECT image whose stored pixels represent counts. The method currently supports Siemens and GE DICOM data and expects acquisition timing, projection, voxel-size, and manufacturer metadata to be present. Do not use an image that the scanner has already calibrated into activity concentration." + ] + }, + { + "cell_type": "markdown", + "id": "conversion-overview", + "metadata": {}, + "source": [ + "## What the conversion does\n", + "\n", + "For each voxel, `make_bqml_suv` converts scanner counts to Bq/mL using the acquisition frame duration, the number of projections, a camera calibration factor, and the voxel volume:\n", + "\n", + "$$A_{voxel} [Bq/mL] = \\frac{counts}{frame\\ duration \\times projections} \\times CF \\times \\frac{10^6}{voxel\\ volume [mL]}$$\n", + "\n", + "Here, `CF` is in MBq/(counts/s). The method rescales the result into signed 16-bit DICOM pixels, records the Bq/mL real-world value mapping, adds injection and patient metadata, applies decay-correction metadata, and assigns new SOP and Series Instance UIDs. The returned table is a summary of the injection and scan timing; the modified image remains in `image.ds` until it is saved." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "imports", + "metadata": {}, + "outputs": [], + "source": [ + "from pytheranostics.dicomtools.dicomtools import DicomModify" + ] + }, + { + "cell_type": "markdown", + "id": "paths-heading", + "metadata": {}, + "source": [ + "## 1. Select the reconstructed counts image\n", + "\n", + "Enter the path to the reconstructed SPECT DICOM image exported by the scanner in counts. The output path is created alongside the input with `_qspect` appended to its filename, leaving the source image unchanged. You may set `output_path` to a different location if preferred." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "paths", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Input: ../../examples/data/testimages/016.dcm\n", + "Output: ./016_qspect.dcm\n" + ] + } + ], + "source": [ + "# Replace this with the path to your reconstructed SPECT DICOM in counts.\n", + "spect_counts_path = \"../../examples/data/testimages/016.dcm\" #Path(\"/path/to/reconstructed_spect_counts.dcm\")\n", + "\n", + "output_path = \"./016_qspect.dcm\"\n", + "\n", + "print(f\"Input: {spect_counts_path}\")\n", + "print(f\"Output: {output_path}\")" + ] + }, + { + "cell_type": "markdown", + "id": "calibration-heading", + "metadata": {}, + "source": [ + "## 2. Supply the camera calibration factor\n", + "\n", + "Use the site-, camera-, collimator-, energy-window-, isotope-, and reconstruction-specific calibration factor established by your quantitative imaging protocol. Its required unit is **MBq/(counts/s)**. The value below belongs only to the example dataset and must not be reused for clinical data without validation." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "calibration", + "metadata": {}, + "outputs": [], + "source": [ + "calibration_factor = 0.10800584442987242 # MBq/(counts/s); example only\n", + "image = DicomModify(str(spect_counts_path), CF=calibration_factor)" + ] + }, + { + "cell_type": "markdown", + "id": "injection-heading", + "metadata": {}, + "source": [ + "## 3. Enter patient and injection information\n", + "\n", + "Transcribe these values from the administration record. Activities are in MBq, weight is in kg, and **height is in cm**. Dates use `YYYYMMDD`; times use 24-hour `HHMM`.\n", + "\n", + "The method decay-corrects both syringe measurements to the injection time and subtracts the residual activity to estimate the net injected activity. The default half-life is 574,300 seconds (Lu-177), and the default radiopharmaceutical label is `Lutetium-PSMA-617`. Pass different values when processing another radionuclide or agent." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "injection-inputs", + "metadata": {}, + "outputs": [], + "source": [ + "weight_kg = 113.4\n", + "height_cm = 178.5\n", + "\n", + "injection_date = \"20220616\"\n", + "pre_inj_activity_mbq = 7450.0\n", + "pre_inj_time = \"0804\"\n", + "post_inj_activity_mbq = 14.4\n", + "post_inj_time = \"0955\"\n", + "injection_time = \"0918\"\n", + "\n", + "# Use a validated correction here only if the dose calibrator requires one.\n", + "activity_meter_scale_factor = 1.0" + ] + }, + { + "cell_type": "markdown", + "id": "convert-heading", + "metadata": {}, + "source": [ + "## 4. Convert the image and inspect the summary\n", + "\n", + "This call modifies the in-memory DICOM dataset. For Siemens data, `n_detectors` defaults to 2 and is used with `NumberOfFramesInRotation` to recover the total projection count. Set it explicitly if the acquisition used a different detector configuration." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "convert", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/esquinas/projects/venvs/.dosimetry/lib/python3.14/site-packages/pydicom/valuerep.py:440: UserWarning: The value length (17) exceeds the maximum length of 16 allowed for VR DS. Please see for allowed values for each VR.\n", + " warn_and_log(msg)\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
patient_idweight_kgheight_cmpre_inj_activity_MBqpre_inj_datetimepost_inj_activity_MBqpost_inj_datetimeinjected_activity_MBqinjection_datetimescan_datetimedelta_t_days
0PR21-CAVA-0016113.4178.57450.02022-06-16 08:04:0014.42022-06-16 09:55:007395.7448952022-06-16 09:18:002022-06-17 11:37:04.4781.09658
\n", + "
" + ], + "text/plain": [ + " patient_id weight_kg height_cm pre_inj_activity_MBq \\\n", + "0 PR21-CAVA-0016 113.4 178.5 7450.0 \n", + "\n", + " pre_inj_datetime post_inj_activity_MBq post_inj_datetime \\\n", + "0 2022-06-16 08:04:00 14.4 2022-06-16 09:55:00 \n", + "\n", + " injected_activity_MBq injection_datetime scan_datetime \\\n", + "0 7395.744895 2022-06-16 09:18:00 2022-06-17 11:37:04.478 \n", + "\n", + " delta_t_days \n", + "0 1.09658 " + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "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", + ")\n", + "\n", + "injection_summary" + ] + }, + { + "cell_type": "markdown", + "id": "save-heading", + "metadata": {}, + "source": [ + "## 5. Save and verify the qSPECT DICOM\n", + "\n", + "Save to the explicit output path only after the conversion succeeds. Then reopen the file and verify the series description and real-world mapping. Quantitative use should also include your site's usual image- and metadata-QC checks." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "save-verify", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Saved: ./016_qspect.dcm\n", + "Series description: QSPECT_Lu PSMA 617,6/17/2022 [WB Recon - AC ]\n", + "Mapped units: Bq/ml\n", + "Real-world slope: 51.68902587890625\n", + "Real-world intercept: 0.0\n" + ] + } + ], + "source": [ + "image.ds.save_as(output_path)\n", + "\n", + "mapping = image.ds.RealWorldValueMappingSequence[0]\n", + "print(f\"Saved: {output_path}\")\n", + "print(f\"Series description: {image.ds.SeriesDescription}\")\n", + "print(f\"Mapped units: {mapping.MeasurementUnitsCodeSequence[0].CodeValue}\")\n", + "print(f\"Real-world slope: {mapping.RealWorldValueSlope}\")\n", + "print(f\"Real-world intercept: {mapping.RealWorldValueIntercept}\")" + ] + }, + { + "cell_type": "markdown", + "id": "qc", + "metadata": {}, + "source": [ + "## Quality-control checklist\n", + "\n", + "Before using the result for quantification, confirm that:\n", + "\n", + "- the input pixels were counts and the reconstruction matches the calibration-factor protocol;\n", + "- the DICOM manufacturer, frame duration, number of projections, and voxel dimensions are correct;\n", + "- patient weight and height use kg and cm, respectively;\n", + "- syringe activities, measurement times, injection time, half-life, and radiopharmaceutical are correct;\n", + "- the calculated injected activity and scan-to-injection interval in `injection_summary` are plausible; and\n", + "- a known phantom or other site-approved test confirms the output values in Bq/mL.\n", + "\n", + "The DICOM stores integer pixels. Consumers must apply the `RealWorldValueSlope` and `RealWorldValueIntercept` in `RealWorldValueMappingSequence` to recover Bq/mL values." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".dosimetry (3.14.4.final.0)", + "language": "python", + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pytheranostics/dicomtools/dicomtools.py b/pytheranostics/dicomtools/dicomtools.py index 209d5c9..1931c46 100644 --- a/pytheranostics/dicomtools/dicomtools.py +++ b/pytheranostics/dicomtools/dicomtools.py @@ -1,94 +1,190 @@ """Utility functions for reading and modifying nuclear medicine DICOM files.""" +import re import time +import warnings from datetime import datetime from pathlib import Path -from typing import Any, List, Tuple +from typing import Any, List, Optional, Tuple import numpy as np import pandas as pd import pydicom import SimpleITK -from pydicom.dataset import Dataset +from pydicom.dataset import Dataset, FileDataset from pydicom.uid import generate_uid from pytheranostics.shared.radioactive_decay import get_activity_at_injection +_MANUFACTURER_ALIASES = { + "siemens": "siemens", + "siemens healthineers": "siemens", + "siemens medical solutions": "siemens", + "siemens medical systems": "siemens", + "siemens nm": "siemens", + "ge": "ge", + "ge healthcare": "ge", + "ge medical systems": "ge", +} + class DicomModify: - """Helper that edits DICOM headers/pixel data for quantitative SPECT studies.""" + """Edit DICOM headers and pixel data for quantitative SPECT studies. + + Parameters + ---------- + fname : str + Path to the source DICOM file. + CF : float + Camera calibration factor used to convert count-rate data to activity + concentration. + """ - def __init__(self, fname, CF): - """Load the DICOM file and store calibration info.""" - self.ds = pydicom.dcmread(fname) - self.CF = CF - self.fname = fname + def __init__(self, fname: str, CF: float) -> None: + """Load a DICOM file and store calibration information. + + Parameters + ---------- + fname : str + Path to the source DICOM file. + CF : float + Camera calibration factor used to convert count-rate data to + activity concentration, in units of MBq/(counts/s). + """ + self.ds: FileDataset = pydicom.dcmread(fname) + self.CF: float = CF + self.fname: str = fname def make_bqml_suv( self, - weight, - height, - injection_date, - pre_inj_activity, - pre_inj_time, - post_inj_activity, - post_inj_time, - injection_time, - activity_meter_scale_factor, - half_life=574300, - radiopharmaceutical="Lutetium-PSMA-617", - n_detectors=2, - ): - """Convert raw counts to BQML/SUV units and update the header accordingly.""" + weight: float, + height: float, + injection_date: str, + pre_inj_activity: float, + pre_inj_time: str, + post_inj_activity: float, + post_inj_time: str, + injection_time: str, + activity_meter_scale_factor: float = 1.0, + half_life: float = 574300, + radiopharmaceutical: str = "Lutetium-PSMA-617", + n_detectors: int = 2, + siemens_frame_duration_fallback: Optional[float] = None, + ) -> pd.DataFrame: + """Convert raw counts to BQML/SUV units and update DICOM metadata. + + Parameters + ---------- + weight : float + Patient weight in kilograms. + height : float + Patient height in centimeters. + injection_date : str + Injection date formatted as ``YYYYMMDD``. + pre_inj_activity : float + Pre-injection syringe activity in MBq. + pre_inj_time : str + Time of the pre-injection syringe measurement formatted as + ``HHMM``. + post_inj_activity : float + Post-injection syringe activity in MBq. + post_inj_time : str + Time of the post-injection syringe measurement formatted as + ``HHMM``. + injection_time : str + Injection time formatted as ``HHMM``. + activity_meter_scale_factor : float + Scale factor applied to the calculated injected activity. Default is 1.0. + half_life : float, optional + Radionuclide half-life in seconds. The default is 574300 for Lu-177. + radiopharmaceutical : str, optional + Radiopharmaceutical name written to the DICOM metadata. The default + is ``"Lutetium-PSMA-617"``. + n_detectors : int, optional + Number of detectors used to correct Siemens projection counts. The + default is 2. + siemens_frame_duration_fallback : float, optional + Frame duration in seconds to use for Siemens data when + ``ActualFrameDuration`` is missing or invalid. By default no + fallback is used and invalid DICOM metadata raises ``ValueError``. + + Returns + ------- + pandas.DataFrame + One-row summary of patient, injection, and scan timing metadata. + + Notes + ----- + This method mutates ``self.ds`` in place. Call :meth:`save` to write the + modified dataset to disk. + """ # Half-life is in seconds + if n_detectors <= 0: + raise ValueError(f"n_detectors must be positive; got {n_detectors}.") + + manufacturer = _normalize_dicom_manufacturer(self.ds) + # Siemens has an issue setting up the times. We are using the Acquisition time which is the time of the start of the last bed to harmonize. - if "siemens" in self.ds.Manufacturer.lower(): - self.ds.SeriesTime = self.ds.AcquisitionTime - self.ds.ContentTime = self.ds.AcquisitionTime - elif "ge" in self.ds.Manufacturer.lower(): # i think it applies to ge as well - self.ds.SeriesTime = self.ds.AcquisitionTime - self.ds.ContentTime = self.ds.AcquisitionTime + if manufacturer == "siemens": + acquisition_time = _require_dicom_text(self.ds, "AcquisitionTime") + self.ds.SeriesTime = acquisition_time + self.ds.ContentTime = acquisition_time + elif manufacturer == "ge": # i think it applies to ge as well + acquisition_time = _require_dicom_text(self.ds, "AcquisitionTime") + self.ds.SeriesTime = acquisition_time + self.ds.ContentTime = acquisition_time # Get the frame duration in seconds - frame_duration = ( - self.ds.RotationInformationSequence[0].ActualFrameDuration / 1000 + rotation_info = _require_dicom_sequence_item( + self.ds, "RotationInformationSequence" ) + frame_duration = _get_frame_duration_seconds( + rotation_info, + manufacturer, + siemens_frame_duration_fallback, + ) + # get number of projections because manufacturers scale by this in the dicomfile - if "siemens" in self.ds.Manufacturer.lower(): - n_proj = ( - self.ds.RotationInformationSequence[0].NumberOfFramesInRotation - * n_detectors + frames_in_rotation = _require_positive_dicom_int( + rotation_info, + "NumberOfFramesInRotation", + context="RotationInformationSequence[0]", + ) + if manufacturer == "siemens": + n_proj = frames_in_rotation * n_detectors + elif manufacturer == "ge": + n_proj = frames_in_rotation + else: + raise ValueError( + "Unsupported DICOM manufacturer for projection scaling: " + f"{_require_dicom_text(self.ds, 'Manufacturer')}" ) - elif "ge" in self.ds.Manufacturer.lower(): - n_proj = self.ds.RotationInformationSequence[0].NumberOfFramesInRotation + # get voxel volume in ml - vox_vol = np.append( - np.asarray(self.ds.PixelSpacing), float(self.ds.SliceThickness) - ) - vox_vol = np.prod(vox_vol / 10) + vox_vol = _get_voxel_volume_ml(self.ds) # Get image in Bq/ml - A = self.ds.pixel_array - A.astype(np.float64) + A = self.ds.pixel_array.astype(np.float64) A = A / (frame_duration * n_proj) * self.CF * 1e6 / vox_vol slope, intercept = dicom_slope_intercept(A) # update the PixelData - A = np.int16((A - intercept) / slope) # GE dicom is signed so np.int16 + # GE DICOM is signed, so use int16. + A = ((A - intercept) / slope).astype(np.int16) # bring the new image to the pixel bytes self.ds.PixelData = A.tobytes() - - self.ds.PixelData + _update_int16_pixel_metadata(self.ds, A) # update DICOM tags # self.ds.Units = 'BQML' - self.ds.SeriesDescription = "QSPECT_" + self.ds.SeriesDescription + series_description = _require_dicom_text(self.ds, "SeriesDescription") + self.ds.SeriesDescription = "QSPECT_" + series_description # add the RealWorldValueMappingSequence tag [0040,9096] - self.ds.add_new([0x0040, 0x9096], "SQ", []) + self.ds.add_new((0x0040, 0x9096), "SQ", []) self.ds.RealWorldValueMappingSequence += [Dataset(), Dataset()] for i in range(2): @@ -102,7 +198,7 @@ def make_bqml_suv( ) self.ds.RealWorldValueMappingSequence[i].LUTLabel = "BQML" - self.ds.RealWorldValueMappingSequence[i].add_new([0x0040, 0x08EA], "SQ", []) + self.ds.RealWorldValueMappingSequence[i].add_new((0x0040, 0x08EA), "SQ", []) self.ds.RealWorldValueMappingSequence[i].MeasurementUnitsCodeSequence += [ Dataset() ] @@ -115,9 +211,9 @@ def make_bqml_suv( self.ds.PatientSize = str(height / 100) # in m self.ds.DecayCorrection = "START" - self.ds.CorrectedImage.insert(0, "DECY") + _prepend_corrected_image_value(self.ds, "DECY") - self.ds.add_new([0x0054, 0x0016], "SQ", []) + self.ds.add_new((0x0054, 0x0016), "SQ", []) self.ds.RadiopharmaceuticalInformationSequence += [Dataset()] # values for net injected activity and injection date and time @@ -132,8 +228,12 @@ def make_bqml_suv( ) total_injected_activity = total_injected_activity * activity_meter_scale_factor - scan_datetime = datetime.strptime( - self.ds.SeriesDate + self.ds.SeriesTime, "%Y%m%d%H%M%S.%f" + series_date = _require_dicom_text(self.ds, "SeriesDate") + series_time = _require_dicom_text(self.ds, "SeriesTime") + scan_datetime = _parse_dicom_date_time( + series_date, + series_time, + context="SeriesDate/SeriesTime", ) delta_scan_inj = (scan_datetime - start_datetime).total_seconds() / ( 60 * 60 * 24 @@ -147,9 +247,9 @@ def make_bqml_suv( ) inj_dic = { - "patient_id": [self.ds.PatientID], + "patient_id": [_require_dicom_text(self.ds, "PatientID")], "weight_kg": [weight], - "height_m": [height], + "height_cm": [height], "pre_inj_activity_MBq": [pre_inj_activity], "pre_inj_datetime": [pre_inj_datetime], "post_inj_activity_MBq": [post_inj_activity], @@ -182,32 +282,49 @@ def make_bqml_suv( ].RadiopharmaceuticalStartDateTime = start_datetime.strftime("%Y%m%d%H%M%S.%f") # for storing as new series data - sop_ins_uid = self.ds.SOPInstanceUID - b = sop_ins_uid.split(".") - b[-1] = str(int(b[-1]) + 1) - self.ds.SOPInstanceUID = ".".join(b) - - ser_ins_uid = self.ds.SeriesInstanceUID - b = ser_ins_uid.split(".") - b.pop() - prefix = ".".join(b) + "." + sop_ins_uid = _require_dicom_text(self.ds, "SOPInstanceUID") + self.ds.SOPInstanceUID = _increment_uid_suffix(sop_ins_uid, "SOPInstanceUID") + + 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.MediaStorageSOPInstaceUID return inj_df - def save(self): - """Persist the modified dataset alongside the original file.""" - self.ds.save_as(f"{self.fname.split('.dcm')[0]}_out.dcm") + def save(self) -> None: + """Persist the modified dataset alongside the original file. + + Notes + ----- + The output path is generated by appending ``"_out"`` to the input file + stem and preserving the original suffix. + """ + path = Path(self.fname) + output_path = path.with_name(f"{path.stem}_out{path.suffix}") + self.ds.save_as(output_path) + + +def dicom_slope_intercept(img: np.ndarray) -> Tuple[float, float]: + """Calculate GE-style slope and intercept values. + Parameters + ---------- + img : numpy.ndarray + Floating-point image array to be quantized into signed 16-bit DICOM + pixel data. -def dicom_slope_intercept(img): - """Calculate GE-style slope/intercept for converting floats to signed int16. + Returns + ------- + tuple of float + Slope and intercept used to recover real-world voxel values from the + stored integer pixel data. - GE PET images are stored as signed int16 values with magnitude limited to - 32767. The computed slope ensures the largest absolute voxel value in the - floating-point array (e.g., MBq/mL) maps to this range once quantized, while - the intercept remains zero (GE convention). + Notes + ----- + GE PET images are stored as signed ``int16`` values with magnitude limited + to 32767. The computed slope maps the largest absolute voxel value to this + range once quantized, while the intercept remains zero. """ max_val = np.max(img) min_val = np.min(img) @@ -218,6 +335,334 @@ def dicom_slope_intercept(img): return float(slope), float(intercept) +def _update_int16_pixel_metadata( + ds: FileDataset, + pixel_array: np.ndarray, +) -> None: + """Update DICOM pixel metadata for signed 16-bit stored pixel data. + + Parameters + ---------- + ds : pydicom.dataset.FileDataset + Dataset whose ``PixelData`` has been rewritten. + pixel_array : numpy.ndarray + Signed 16-bit stored-value pixel array written to ``PixelData``. + + Notes + ----- + Quantitative scaling is stored in ``RealWorldValueMappingSequence``. The + conventional rescale transformation is therefore removed to prevent stale + source metadata or readers from applying an additional scale. + """ + ds.BitsAllocated = 16 + ds.BitsStored = 16 + ds.HighBit = 15 + ds.PixelRepresentation = 1 + ds.SmallestImagePixelValue = int(pixel_array.min()) + ds.LargestImagePixelValue = int(pixel_array.max()) + for keyword in ("RescaleSlope", "RescaleIntercept", "RescaleType"): + if keyword in ds: + del ds[keyword] + + +def _require_dicom_value( + ds: Dataset, + keyword: str, + context: str = "DICOM dataset", +) -> Any: + """Return a required DICOM value or raise a clear validation error.""" + if not hasattr(ds, keyword): + raise ValueError(f"Required DICOM tag '{keyword}' is missing from {context}.") + + value = getattr(ds, keyword) + if value is None: + raise ValueError(f"Required DICOM tag '{keyword}' is empty in {context}.") + if isinstance(value, str) and not value.strip(): + raise ValueError(f"Required DICOM tag '{keyword}' is empty in {context}.") + + return value + + +def _get_frame_duration_seconds( + rotation_info: Dataset, + manufacturer: str, + siemens_fallback: Optional[float] = None, +) -> float: + """Return the DICOM frame duration, with an opt-in Siemens fallback. + + ``ActualFrameDuration`` is stored in milliseconds, while the fallback is + supplied in seconds to match the value used by the conversion calculation. + """ + if siemens_fallback is not None: + try: + siemens_fallback = float(siemens_fallback) + except (TypeError, ValueError) as exc: + raise ValueError( + "siemens_frame_duration_fallback must be a positive number " + f"of seconds; got {siemens_fallback!r}." + ) from exc + if not np.isfinite(siemens_fallback) or siemens_fallback <= 0: + raise ValueError( + "siemens_frame_duration_fallback must be a positive number " + f"of seconds; got {siemens_fallback!r}." + ) + + try: + duration_ms = _require_positive_dicom_float( + rotation_info, + "ActualFrameDuration", + context="RotationInformationSequence[0]", + ) + except ValueError: + if manufacturer != "siemens" or siemens_fallback is None: + raise + warnings.warn( + "Siemens ActualFrameDuration is missing or invalid; using " + f"siemens_frame_duration_fallback={siemens_fallback} seconds.", + RuntimeWarning, + stacklevel=2, + ) + return siemens_fallback + + return duration_ms / 1000 + + +def _require_dicom_text( + ds: Dataset, + keyword: str, + context: str = "DICOM dataset", +) -> str: + """Return a required DICOM value as stripped text.""" + value = _require_dicom_value(ds, keyword, context) + text = str(value).strip() + if not text: + raise ValueError(f"Required DICOM tag '{keyword}' is empty in {context}.") + + return text + + +def _normalize_dicom_manufacturer(ds: Dataset) -> str: + """Return the supported canonical manufacturer name for a DICOM dataset.""" + manufacturer = _require_dicom_text(ds, "Manufacturer") + normalized = " ".join(manufacturer.casefold().replace("_", " ").split()) + + if normalized in _MANUFACTURER_ALIASES: + return _MANUFACTURER_ALIASES[normalized] + + # Manufacturer values commonly include a modality, division, or model name. + if "siemens" in normalized: + return "siemens" + if re.search(r"(?:^|[^a-z0-9])ge(?:$|[^a-z0-9])", normalized): + return "ge" + + supported = ", ".join(sorted(set(_MANUFACTURER_ALIASES.values()))) + raise ValueError( + f"Unsupported DICOM manufacturer for QSPECT conversion: " + f"{manufacturer!r}. Supported manufacturers: {supported}." + ) + + +def _require_dicom_float( + ds: Dataset, + keyword: str, + context: str = "DICOM dataset", +) -> float: + """Return a required DICOM value as a float.""" + value = _require_dicom_value(ds, keyword, context) + try: + return float(value) + except (TypeError, ValueError) as exc: + raise ValueError( + f"Required DICOM tag '{keyword}' in {context} must be numeric; " + f"got {value!r}." + ) from exc + + +def _require_dicom_int( + ds: Dataset, + keyword: str, + context: str = "DICOM dataset", +) -> int: + """Return a required DICOM value as an integer.""" + value = _require_dicom_value(ds, keyword, context) + try: + return int(value) + except (TypeError, ValueError) as exc: + raise ValueError( + f"Required DICOM tag '{keyword}' in {context} must be an integer; " + f"got {value!r}." + ) from exc + + +def _require_positive_dicom_float( + ds: Dataset, + keyword: str, + context: str = "DICOM dataset", +) -> float: + """Return a required DICOM value as a positive float.""" + value = _require_dicom_float(ds, keyword, context) + if value <= 0: + raise ValueError( + f"Required DICOM tag '{keyword}' in {context} must be positive; " + f"got {value}." + ) + + return value + + +def _require_positive_dicom_int( + ds: Dataset, + keyword: str, + context: str = "DICOM dataset", +) -> int: + """Return a required DICOM value as a positive integer.""" + value = _require_dicom_int(ds, keyword, context) + if value <= 0: + raise ValueError( + f"Required DICOM tag '{keyword}' in {context} must be positive; " + f"got {value}." + ) + + return value + + +def _require_dicom_sequence_item( + ds: Dataset, + keyword: str, + context: str = "DICOM dataset", + index: int = 0, +) -> Dataset: + """Return a required item from a DICOM sequence.""" + sequence = _require_dicom_value(ds, keyword, context) + try: + item = sequence[index] + except (IndexError, TypeError) as exc: + raise ValueError( + f"Required DICOM sequence '{keyword}' in {context} must contain " + f"item {index}." + ) from exc + + if not isinstance(item, Dataset): + raise ValueError( + f"Item {index} of DICOM sequence '{keyword}' in {context} must be " + f"a Dataset; got {type(item).__name__}." + ) + + return item + + +def _get_voxel_volume_ml(ds: Dataset) -> float: + """Return voxel volume in milliliters from required spacing tags.""" + pixel_spacing = _require_dicom_value(ds, "PixelSpacing") + try: + spacing_mm = np.asarray(pixel_spacing, dtype=float) + except (TypeError, ValueError) as exc: + raise ValueError( + "Required DICOM tag 'PixelSpacing' must contain numeric row and " + f"column spacing values; got {pixel_spacing!r}." + ) from exc + + if spacing_mm.size < 2: + raise ValueError( + "Required DICOM tag 'PixelSpacing' must contain row and column " + f"spacing values; got {pixel_spacing!r}." + ) + + if np.any(spacing_mm[:2] <= 0): + raise ValueError( + "Required DICOM tag 'PixelSpacing' must contain positive row and " + f"column spacing values; got {pixel_spacing!r}." + ) + + slice_thickness_mm = _require_positive_dicom_float(ds, "SliceThickness") + voxel_dimensions_mm = np.append(spacing_mm[:2], slice_thickness_mm) + voxel_volume_ml = float(np.prod(voxel_dimensions_mm / 10)) + if voxel_volume_ml <= 0: + raise ValueError( + "Voxel volume calculated from PixelSpacing and SliceThickness must " + f"be positive; got {voxel_volume_ml} ml." + ) + + return voxel_volume_ml + + +def _parse_dicom_date_time(date: str, time_value: str, context: str) -> datetime: + """Parse DICOM DA/TM text into a datetime.""" + for time_format in ("%H%M%S.%f", "%H%M%S", "%H%M"): + try: + return datetime.strptime(date + time_value, "%Y%m%d" + time_format) + except ValueError: + continue + + raise ValueError( + f"DICOM {context} must be formatted as YYYYMMDD and HHMM[SS[.ffffff]]; " + f"got {date!r} and {time_value!r}." + ) + + +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 _uid_prefix_for_generated_uid(uid: str, keyword: str) -> str: + """Return a UID prefix suitable for pydicom.uid.generate_uid.""" + uid_parts = uid.split(".") + if len(uid_parts) < 2 or "" in uid_parts: + raise ValueError( + f"Required DICOM tag '{keyword}' must be a dot-separated UID with " + f"at least two components; got {uid!r}." + ) + if not all(part.isdigit() for part in uid_parts): + raise ValueError( + f"Required DICOM tag '{keyword}' must contain only numeric UID " + f"components; got {uid!r}." + ) + + prefix = ".".join(uid_parts[:-1]) + "." + if len(prefix) > 54: + raise ValueError( + f"UID prefix derived from DICOM tag '{keyword}' is too long for " + f"generate_uid; got {len(prefix)} characters." + ) + + return prefix + + +def _prepend_corrected_image_value(ds: Dataset, value: str) -> None: + """Prepend a CorrectedImage value, creating the tag if needed.""" + corrected_image = getattr(ds, "CorrectedImage", []) + if corrected_image is None: + values = [] + elif isinstance(corrected_image, str): + values = [corrected_image] + else: + try: + values = list(corrected_image) + except TypeError as exc: + raise ValueError( + "DICOM tag 'CorrectedImage' must be text or a multi-value " + f"sequence; got {corrected_image!r}." + ) from exc + + if value not in values: + values.insert(0, value) + + ds.CorrectedImage = values + + def generate_basic_dcm_tags( img_size: Tuple[int, int, int], slice_thickness: float, diff --git a/pytheranostics/dosimetry/dvk.py b/pytheranostics/dosimetry/dvk.py index d7bb4ef..b8a8171 100644 --- a/pytheranostics/dosimetry/dvk.py +++ b/pytheranostics/dosimetry/dvk.py @@ -3,9 +3,10 @@ from __future__ import annotations import csv +import hashlib +import json import logging import re -import shutil import tempfile import urllib.error import urllib.request @@ -23,6 +24,9 @@ _DATA_DIR = Path(__file__).resolve().parent.parent / "data" / "voxel_kernels" _ZENODO_RECORD_ID = "7596345" +_ZENODO_DOWNLOAD_ATTEMPTS = 3 +_ZENODO_DOWNLOAD_TIMEOUT_S = 120 +_ZENODO_USER_AGENT = "PyTheranostics dose-kernel downloader" _DEFAULT_CROP_SIZE_MM = 25.0 _CSV_PATTERN = re.compile( r"^(?P[0-9]+[A-Za-z]+)" @@ -222,30 +226,96 @@ def _download_kernels_from_zenodo(isotope: str, kernel_dir: Path) -> None: If the archive download fails for a network-related reason. """ zenodo_name = _zenodo_isotope_name(isotope) - url = ( - f"https://zenodo.org/records/{_ZENODO_RECORD_ID}/files/" - f"{zenodo_name}.zip?download=1" - ) + archive_name = f"{zenodo_name}.zip" + record_url = f"https://zenodo.org/api/records/{_ZENODO_RECORD_ID}" kernel_dir.mkdir(parents=True, exist_ok=True) try: - with urllib.request.urlopen(url, timeout=120) as response: - with tempfile.NamedTemporaryFile(suffix=".zip", delete=False) as tmp_file: - shutil.copyfileobj(response, tmp_file) - archive_path = Path(tmp_file.name) + request = urllib.request.Request( + record_url, headers={"User-Agent": _ZENODO_USER_AGENT} + ) + with urllib.request.urlopen( + request, timeout=_ZENODO_DOWNLOAD_TIMEOUT_S + ) as response: + record = json.load(response) except urllib.error.HTTPError as exc: if exc.code == 404: logger.error( - "No dose voxel kernel archive is available on Zenodo for isotope %s.", - isotope, + "Dose voxel kernel record %s is unavailable on Zenodo.", + _ZENODO_RECORD_ID, ) raise NotImplementedError( f"Dose voxel kernels for isotope '{isotope}' are not available on Zenodo." ) from exc raise + file_record = next( + ( + item + for item in record.get("files", []) + if str(item.get("key", "")).casefold() == archive_name.casefold() + ), + None, + ) + if file_record is None: + raise NotImplementedError( + f"Dose voxel kernels for isotope '{isotope}' are not available on Zenodo." + ) + + links = file_record.get("links", {}) + download_url = links.get("content") or links.get("download") or links.get("self") + if not download_url: + raise FileNotFoundError( + f"Zenodo file '{archive_name}' does not include a download link." + ) + + expected_size = int(file_record["size"]) + checksum_algorithm, expected_checksum = str(file_record["checksum"]).split( + ":", maxsplit=1 + ) + if checksum_algorithm != "md5": + raise ValueError( + f"Unsupported Zenodo checksum algorithm: {checksum_algorithm!r}." + ) + + with tempfile.NamedTemporaryFile(suffix=".zip", delete=False) as tmp_file: + archive_path = Path(tmp_file.name) + try: + failures = [] + for attempt in range(1, _ZENODO_DOWNLOAD_ATTEMPTS + 1): + request = urllib.request.Request( + download_url, headers={"User-Agent": _ZENODO_USER_AGENT} + ) + digest = hashlib.md5() + downloaded_size = 0 + with urllib.request.urlopen( + request, timeout=_ZENODO_DOWNLOAD_TIMEOUT_S + ) as response, archive_path.open("wb") as archive_file: + while chunk := response.read(1024 * 1024): + archive_file.write(chunk) + digest.update(chunk) + downloaded_size += len(chunk) + + downloaded_checksum = digest.hexdigest() + if ( + downloaded_size == expected_size + and downloaded_checksum == expected_checksum + ): + break + + failures.append( + f"attempt {attempt}: size={downloaded_size}/{expected_size}, " + f"md5={downloaded_checksum}/{expected_checksum}" + ) + else: + raise OSError( + f"Downloaded kernel archive '{archive_name}' failed Zenodo " + f"integrity checks after {_ZENODO_DOWNLOAD_ATTEMPTS} attempts: " + + "; ".join(failures) + ) + with zipfile.ZipFile(archive_path) as archive: csv_members = [ member @@ -262,7 +332,7 @@ def _download_kernels_from_zenodo(isotope: str, kernel_dir: Path) -> None: finally: archive_path.unlink(missing_ok=True) - logger.info("Downloaded dose voxel kernels for %s from %s", isotope, url) + logger.info("Downloaded dose voxel kernels for %s from %s", isotope, download_url) def _ensure_kernel_files_available( diff --git a/pytheranostics/imaging_tools/tools.py b/pytheranostics/imaging_tools/tools.py index 8988ee8..e26abdf 100644 --- a/pytheranostics/imaging_tools/tools.py +++ b/pytheranostics/imaging_tools/tools.py @@ -313,7 +313,8 @@ def apply_qspect_dcm_scaling( dcm_data.RealWorldValueMappingSequence[0].RealWorldValueIntercept ) logger.info( - "Applying DICOM quantitative scaling for NM: slope=%s intercept=%s.", + "Applying NM quantitative scaling obtained from DICOM " + "RealWorldValueMappingSequence: slope=%s intercept=%s.", slope, intercept, ) diff --git a/pytheranostics/shared/radioactive_decay.py b/pytheranostics/shared/radioactive_decay.py index d259ec8..2982390 100644 --- a/pytheranostics/shared/radioactive_decay.py +++ b/pytheranostics/shared/radioactive_decay.py @@ -6,11 +6,13 @@ def decay_act(a_initial, delta_t, half_life): - """Return decayed activity after `delta_t` given the half-life.""" + """Return activity corrected by ``delta_t`` given the half-life. + + Positive time intervals decay activity forward in time. Negative time + intervals extrapolate a later measurement backward to an earlier time. + """ if np.any(np.asarray(a_initial) < 0): raise ValueError("a_initial must be positive") - if np.any(np.asarray(delta_t) < 0): - raise ValueError("delta_t must be positive") if np.any(np.asarray(half_life) < 0): raise ValueError("half_life must be positive") @@ -40,6 +42,19 @@ def get_activity_at_injection( injection_date + injection_time + "00.00", "%Y%m%d%H%M%S.%f" ) + if pre_datetime > inj_datetime: + raise ValueError( + "pre_inj_time must be at or before injection_time; " + f"got pre_inj_time={pre_inj_time!r} and " + f"injection_time={injection_time!r}." + ) + if post_datetime < inj_datetime: + raise ValueError( + "post_inj_time must be at or after injection_time; " + f"got post_inj_time={post_inj_time!r} and " + f"injection_time={injection_time!r}." + ) + delta_inj_pre = (inj_datetime - pre_datetime).total_seconds() delta_post_inj = (inj_datetime - post_datetime).total_seconds() diff --git a/tests/data/voxel_dosimetry_validation/README.md b/tests/data/voxel_dosimetry_validation/README.md index 55770a7..1ea557e 100644 --- a/tests/data/voxel_dosimetry_validation/README.md +++ b/tests/data/voxel_dosimetry_validation/README.md @@ -9,7 +9,7 @@ validation test: - `dosimetry_fit_defaults.json` (optional): custom ROI fit defaults for the validation run The precomputed RT-STRUCT DICOM files are downloaded on demand from Zenodo -record `21893683`: https://zenodo.org/records/21893683. +record `22036568`: https://zenodo.org/records/22036568. The validation test is implemented in `tests/test_voxel_dosimetry_validation.py`. GitHub Actions may skip this test when the external SNMMI Deep Blue dataset diff --git a/tests/test_dicomtools.py b/tests/test_dicomtools.py new file mode 100644 index 0000000..65e4d84 --- /dev/null +++ b/tests/test_dicomtools.py @@ -0,0 +1,94 @@ +import numpy as np +import pytest +from pydicom.dataset import Dataset + +from pytheranostics.dicomtools.dicomtools import ( + _get_frame_duration_seconds, + _normalize_dicom_manufacturer, + _update_int16_pixel_metadata, +) + + +@pytest.mark.parametrize( + ("manufacturer", "expected"), + [ + ("SIEMENS", "siemens"), + ("SIEMENS NM", "siemens"), + ("Siemens Healthineers", "siemens"), + ("ACME / SIEMENS Symbia", "siemens"), + ("GE MEDICAL SYSTEMS", "ge"), + ("GE Healthcare NM 870", "ge"), + ], +) +def test_normalize_dicom_manufacturer_aliases(manufacturer, expected): + dataset = Dataset() + dataset.Manufacturer = manufacturer + + assert _normalize_dicom_manufacturer(dataset) == expected + + +def test_normalize_dicom_manufacturer_rejects_unknown_value(): + dataset = Dataset() + dataset.Manufacturer = "Unknown Scanner Company" + + with pytest.raises(ValueError, match="Unsupported DICOM manufacturer"): + _normalize_dicom_manufacturer(dataset) + + +def test_normalize_dicom_manufacturer_does_not_match_ge_inside_a_word(): + dataset = Dataset() + dataset.Manufacturer = "Edge Imaging" + + with pytest.raises(ValueError, match="Unsupported DICOM manufacturer"): + _normalize_dicom_manufacturer(dataset) + + +def test_get_frame_duration_seconds_prefers_valid_dicom_value(): + rotation_info = Dataset() + rotation_info.ActualFrameDuration = 12000 + + assert _get_frame_duration_seconds(rotation_info, "siemens", 15) == 12 + + +def test_get_frame_duration_seconds_uses_opt_in_siemens_fallback(): + rotation_info = Dataset() + + with pytest.warns(RuntimeWarning, match="using siemens_frame_duration_fallback"): + duration = _get_frame_duration_seconds(rotation_info, "siemens", 15) + + assert duration == 15 + + +def test_get_frame_duration_seconds_remains_strict_without_fallback(): + with pytest.raises(ValueError, match="ActualFrameDuration"): + _get_frame_duration_seconds(Dataset(), "siemens") + + +def test_get_frame_duration_seconds_does_not_apply_siemens_fallback_to_ge(): + with pytest.raises(ValueError, match="ActualFrameDuration"): + _get_frame_duration_seconds(Dataset(), "ge", 15) + + +@pytest.mark.parametrize("fallback", [0, -1, float("nan"), "invalid"]) +def test_get_frame_duration_seconds_rejects_invalid_fallback(fallback): + rotation_info = Dataset() + rotation_info.ActualFrameDuration = 12000 + + with pytest.raises(ValueError, match="must be a positive number"): + _get_frame_duration_seconds(rotation_info, "siemens", fallback) + + +def test_update_int16_pixel_metadata_removes_conventional_rescale_mapping(): + dataset = Dataset() + dataset.RescaleSlope = "42.0" + dataset.RescaleIntercept = "7.0" + dataset.RescaleType = "US" + pixels = np.array([[-2, 0, 3]], dtype=np.int16) + + _update_int16_pixel_metadata(dataset, pixels) + + assert "RescaleSlope" not in dataset + assert "RescaleIntercept" not in dataset + assert "RescaleType" not in dataset + assert dataset.SmallestImagePixelValue == -2 + assert dataset.LargestImagePixelValue == 3 diff --git a/tests/test_dvk.py b/tests/test_dvk.py index cd90dd2..7112fc1 100644 --- a/tests/test_dvk.py +++ b/tests/test_dvk.py @@ -1,5 +1,9 @@ """Tests for DoseVoxelKernel CSV loading and selection.""" +import hashlib +import io +import json +import zipfile from pathlib import Path import numpy as np @@ -95,3 +99,44 @@ def _fake_download(isotope: str, kernel_dir_arg: Path) -> None: assert created["called"] is True assert kernel.kernel.shape == (3, 3, 3) + + +def test_kernel_download_retries_truncated_archive( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +): + """A truncated Zenodo response should be retried before ZIP extraction.""" + archive_buffer = io.BytesIO() + with zipfile.ZipFile(archive_buffer, mode="w") as archive: + archive.writestr( + "177Lu/177LU_6.647D_99%_5mm_3x3x3.csv", + "1,2\n3,4\n", + ) + archive_bytes = archive_buffer.getvalue() + checksum = hashlib.md5(archive_bytes).hexdigest() + record = { + "files": [ + { + "key": "177Lu.zip", + "size": len(archive_bytes), + "checksum": f"md5:{checksum}", + "links": {"self": "https://example.test/177Lu.zip"}, + } + ] + } + responses = [ + json.dumps(record).encode(), + archive_bytes[: len(archive_bytes) // 2], + archive_bytes, + ] + + def fake_urlopen(request, timeout): + assert timeout == dvk._ZENODO_DOWNLOAD_TIMEOUT_S + return io.BytesIO(responses.pop(0)) + + monkeypatch.setattr(dvk.urllib.request, "urlopen", fake_urlopen) + + kernel_dir = tmp_path / "voxel_kernels" + dvk._download_kernels_from_zenodo("Lu177", kernel_dir) + + assert not responses + assert (kernel_dir / "177Lu" / "177LU_6.647D_99%_5mm_3x3x3.csv").is_file() diff --git a/tests/test_imaging_tools.py b/tests/test_imaging_tools.py index e2a55c1..9ffd471 100644 --- a/tests/test_imaging_tools.py +++ b/tests/test_imaging_tools.py @@ -1,14 +1,41 @@ """Tests for imaging tools utilities.""" +import logging import shutil import numpy as np import pytest import SimpleITK +from pydicom.dataset import Dataset +from pydicom.sequence import Sequence from pytheranostics.imaging_tools import tools +def test_apply_qspect_scaling_logs_real_world_value_mapping_source( + tmp_path, monkeypatch, caplog +): + dicom_path = tmp_path / "qspect.dcm" + dicom_path.touch() + + mapping = Dataset() + mapping.RealWorldValueSlope = 2.0 + mapping.RealWorldValueIntercept = 3.0 + dataset = Dataset() + dataset.Modality = "NM" + dataset.DecayCorrection = "START" + dataset.RealWorldValueMappingSequence = Sequence([mapping]) + monkeypatch.setattr(tools.pydicom, "dcmread", lambda *args, **kwargs: dataset) + + input_array = np.arange(8, dtype=np.float32).reshape(2, 2, 2) + image = SimpleITK.GetImageFromArray(input_array) + with caplog.at_level(logging.INFO, logger=tools.__name__): + scaled = tools.apply_qspect_dcm_scaling(image, str(tmp_path)) + + np.testing.assert_allclose(SimpleITK.GetArrayFromImage(scaled), input_array * 2 + 3) + assert "obtained from DICOM RealWorldValueMappingSequence" in caplog.text + + def test_load_metadata_from_sample_spect_folder(spect_example_dir, tmp_path): """Ensure metadata extraction works on bundled DICOM samples.""" single_case_dir = tmp_path / "spect_case" diff --git a/tests/test_radioactive_decay.py b/tests/test_radioactive_decay.py index 45f34ad..5ab9d03 100644 --- a/tests/test_radioactive_decay.py +++ b/tests/test_radioactive_decay.py @@ -3,7 +3,7 @@ import numpy as np import pytest -from pytheranostics.shared.radioactive_decay import decay_act +from pytheranostics.shared.radioactive_decay import decay_act, get_activity_at_injection def test_decay_act(): @@ -31,28 +31,66 @@ def test_decay_act_array(): assert np.allclose(result, expected, rtol=1e-10) +def test_decay_act_negative_time_extrapolates_backward(): + """A negative interval back-corrects a later activity measurement.""" + measured_activity = 500.0 + half_life = 6.0 + delta_t = -2.0 + + expected = measured_activity * np.exp(-np.log(2) * delta_t / half_life) + + assert np.isclose( + decay_act(measured_activity, delta_t, half_life), expected, rtol=1e-10 + ) + + def test_invalid_inputs(): """Test that invalid inputs raise appropriate errors.""" with pytest.raises(ValueError): decay_act(-1000, 6.0, 12.0) # Negative activity with pytest.raises(ValueError): - decay_act(1000, -6.0, 12.0) # Negative half-life - - with pytest.raises(ValueError): - decay_act(1000, 6.0, -12.0) # Negative time - - -# def test_get_activity_at_injection(): -# """Test the get_activity_at_injection function.""" -# current_activity = 500 # MBq -# half_life = 6.0 # hours -# time_since_injection = 12.0 # hours - -# expected = current_activity * np.exp(np.log(2) * time_since_injection / half_life) -# result = decay_act(current_activity, time_since_injection, half_life) -# result = get_activity_at_injection( -# current_activity, half_life, time_since_injection -# ) - -# assert np.isclose(result, expected, rtol=1e-10) + decay_act(1000, 6.0, -12.0) # Negative half-life + + +def test_get_activity_at_injection_with_post_measurement_after_injection(): + """Correct pre/post syringe readings to the common injection time.""" + half_life = 574300.0 + + injection_datetime, injected_activity = get_activity_at_injection( + injection_date="20220616", + pre_inj_activity=7450.0, + pre_inj_time="0804", + post_inj_activity=14.4, + post_inj_time="0955", + injection_time="0918", + half_life=half_life, + ) + + expected_pre = 7450.0 * np.exp(-np.log(2) * (74 * 60) / half_life) + expected_post = 14.4 * np.exp(np.log(2) * (37 * 60) / half_life) + + assert injection_datetime.isoformat() == "2022-06-16T09:18:00" + assert np.isclose(injected_activity, expected_pre - expected_post, rtol=1e-10) + + +@pytest.mark.parametrize( + ("pre_inj_time", "post_inj_time", "error_message"), + [ + ("0930", "0955", "pre_inj_time must be at or before injection_time"), + ("0804", "0900", "post_inj_time must be at or after injection_time"), + ], +) +def test_get_activity_at_injection_rejects_invalid_chronology( + pre_inj_time, post_inj_time, error_message +): + with pytest.raises(ValueError, match=error_message): + get_activity_at_injection( + injection_date="20220616", + pre_inj_activity=7450.0, + pre_inj_time=pre_inj_time, + post_inj_activity=14.4, + post_inj_time=post_inj_time, + injection_time="0918", + half_life=574300.0, + ) diff --git a/tests/test_voxel_dosimetry_validation.py b/tests/test_voxel_dosimetry_validation.py index a630c13..291b750 100644 --- a/tests/test_voxel_dosimetry_validation.py +++ b/tests/test_voxel_dosimetry_validation.py @@ -2,6 +2,7 @@ from __future__ import annotations +import hashlib import json import os import re @@ -29,11 +30,12 @@ ) _EXPECTED_RESULTS_CSV: Final[Path] = _VALIDATION_ASSETS / "expected_results.csv" _EXPECTED_DF_AD_CSV: Final[Path] = _VALIDATION_ASSETS / "expected_df_ad.csv" -_RTSTRUCT_ZENODO_RECORD_ID: Final[str] = "21893683" +_RTSTRUCT_ZENODO_RECORD_ID: Final[str] = "22036568" _RTSTRUCT_ZENODO_RECORD_API_URL: Final[str] = ( f"https://zenodo.org/api/records/{_RTSTRUCT_ZENODO_RECORD_ID}" ) _RTSTRUCT_DOWNLOAD_TIMEOUT_S: Final[int] = 120 +_RTSTRUCT_DOWNLOAD_ATTEMPTS: Final[int] = 3 _ZENODO_USER_AGENT: Final[str] = "PyTheranostics test suite" _RTSTRUCT_SOP_CLASS_UID: Final[str] = "1.2.840.10008.5.1.4.1.1.481.3" _EXPECTED_RTSTRUCT_FILENAMES: Final[tuple[str, ...]] = ( @@ -91,12 +93,47 @@ def _copy_optional_validation_configs(project_base: Path) -> None: shutil.copy2(src, project_base / config_name) -def _download_url_to_file(url: str, output_path: Path) -> None: - request = urllib.request.Request(url, headers={"User-Agent": _ZENODO_USER_AGENT}) - with urllib.request.urlopen( - request, timeout=_RTSTRUCT_DOWNLOAD_TIMEOUT_S - ) as response, output_path.open("wb") as output_file: - shutil.copyfileobj(response, output_file) +def _download_url_to_file( + url: str, + output_path: Path, + *, + expected_size: int, + expected_checksum: str, +) -> None: + checksum_algorithm, checksum_value = expected_checksum.split(":", maxsplit=1) + if checksum_algorithm != "md5": + raise ValueError( + f"Unsupported Zenodo checksum algorithm: {checksum_algorithm!r}." + ) + + failures = [] + for attempt in range(1, _RTSTRUCT_DOWNLOAD_ATTEMPTS + 1): + request = urllib.request.Request( + url, headers={"User-Agent": _ZENODO_USER_AGENT} + ) + digest = hashlib.md5() + downloaded_size = 0 + with urllib.request.urlopen( + request, timeout=_RTSTRUCT_DOWNLOAD_TIMEOUT_S + ) as response, output_path.open("wb") as output_file: + while chunk := response.read(1024 * 1024): + output_file.write(chunk) + digest.update(chunk) + downloaded_size += len(chunk) + + downloaded_checksum = digest.hexdigest() + if downloaded_size == expected_size and downloaded_checksum == checksum_value: + return + + failures.append( + f"attempt {attempt}: size={downloaded_size}/{expected_size}, " + f"md5={downloaded_checksum}/{checksum_value}" + ) + + raise OSError( + f"Downloaded file '{output_path.name}' failed Zenodo integrity checks " + f"after {_RTSTRUCT_DOWNLOAD_ATTEMPTS} attempts: " + "; ".join(failures) + ) def _validate_rtstruct_file(rtstruct_path: Path) -> None: @@ -170,7 +207,13 @@ def _download_rtstruct_assets(project_base: Path) -> list[Path]: f"Zenodo file '{filename}' does not include a download link." ) rtstruct_path = target_dir / filename - _download_url_to_file(download_url, rtstruct_path) + file_record = files_by_name[filename] + _download_url_to_file( + download_url, + rtstruct_path, + expected_size=int(file_record["size"]), + expected_checksum=str(file_record["checksum"]), + ) _validate_rtstruct_file(rtstruct_path) return [target_dir / filename for filename in _EXPECTED_RTSTRUCT_FILENAMES]