diff --git a/README.md b/README.md index 15f800e..47ee6d4 100644 --- a/README.md +++ b/README.md @@ -15,8 +15,8 @@ You get: time-varying acceleration). - **`library/`** — one well-documented class per engineering domain (dynamics, statics, thermo, fluids, materials, gears, electrical, - mechatronics, modal/vibration, controls) plus a `Utility` class for files - and dialogs. + mechatronics, modal/vibration, controls), a set of radiology / + medical-imaging classes, and a `Utility` class for files and dialogs. Every method uses keyword-only arguments and carries a docstring, so the code itself doubles as a reference you can read and learn from. @@ -269,6 +269,68 @@ modal.fit_damped_response(time=t, displacement=x, guess_frequency_hz=3.0) That simulate-then-fit round-trip is checked in the test suite: the fit must return the parameters the simulation started with. +## Radiology & medical imaging + +Six classes cover the physics and image handling a radiology student +(radiologic technologist, resident, or medical physics student) works with +every day. Run the guided tour with `uv run radiology_demo.py`. It prints one +worked example per modality and saves `radiology.png`, which shows a real CT +slice in abdomen, lung, and bone windows. + +| Class | Module | What it covers | +| --- | --- | --- | +| `Radiation` | `library/radiation.py` | Inverse square law, attenuation (mu, HVL, TVL, shielding), photon energy, radioactive and effective half-life, equivalent/effective dose with ICRP 103 weighting factors | +| `Radiography` | `library/radiography.py` | mAs, 15% rule, exposure-maintenance (SID) formula, grid conversion, heat units, magnification, geometric unsharpness, detector Nyquist limit, deviation index | +| `CT` | `library/ct.py` | Hounsfield units, DICOM rescale, window width/level with presets, rough HU tissue lookup, DLP and effective dose estimate | +| `MRI` | `library/mri.py` | Larmor frequency, T1 recovery, T2 decay, spin-echo signal, T1/T2/PD weighting from TR/TE | +| `Ultrasound` | `library/ultrasound.py` | Wavelength, echo depth, acoustic impedance and reflection, attenuation, axial resolution, PRF depth limit, Doppler shift and velocity | +| `Dicom` | `library/dicom.py` | Read DICOM with pydicom, summarize technical tags, pixels in HU, pixel spacing, basic de-identification, PNG export | + +```python +from pydicom.data import get_testdata_file # sample images, no patient data +from library import CT, Dicom, Radiography, Radiation + +radiation = Radiation() +radiation.effective_dose(organ_equivalent_doses={"thyroid": 10.0}) # -> 0.4 mSv + +Radiography().exposure_maintenance_mas(original_mas=10, original_sid=100, new_sid=180) # 32.4 + +dicom, ct = Dicom(), CT() +ds = dicom.read(path=get_testdata_file("CT_small.dcm")) +hu = dicom.pixel_array_in_units(dataset=ds) # stored values -> Hounsfield units +lung = ct.apply_window(hounsfield=hu, preset="lung") # whole slice windowed in one call +safe = dicom.anonymize(dataset=ds) # blank PHI before sharing +``` + +> **Teaching use only.** These are textbook calculations for learning. They +> are not validated for clinical decisions, patient dosimetry, or shielding +> design. Use your department's protocols and a qualified medical physicist +> for those. `Dicom.anonymize` removes common direct identifiers only, which +> is *not* full de-identification under DICOM PS3.15 or HIPAA. + +### Python packages for medical imaging + +`pydicom` is already installed. Add the others as you need them with `uv add `. + +| Package | Use it for | +| --- | --- | +| [pydicom](https://pydicom.github.io/) | Read, edit, and write DICOM files (installed) | +| [SimpleITK](https://simpleitk.org/) | 3-D volumes, resampling, registration, segmentation; reads whole DICOM series | +| [nibabel](https://nipy.org/nibabel/) | NIfTI and other neuroimaging formats (research MRI) | +| [scikit-image](https://scikit-image.org/) | Filtering, thresholding, edge detection, Radon transform / CT reconstruction demos | +| [highdicom](https://highdicom.readthedocs.io/) | Create DICOM segmentations, structured reports, and annotations | +| [pynetdicom](https://pydicom.github.io/pynetdicom/) | Talk to a PACS (C-FIND / C-MOVE / C-STORE) | +| [MONAI](https://monai.io/) | Deep learning for medical imaging (built on PyTorch) | +| [TotalSegmentator](https://github.com/wasserth/TotalSegmentator) | Pretrained automatic segmentation of 100+ structures on CT/MR | +| [pyradiomics](https://pyradiomics.readthedocs.io/) | Extract quantitative radiomic features from images + masks | + +**Learning resources:** [Radiopaedia](https://radiopaedia.org/) (cases and +physics articles), [pydicom user guide](https://pydicom.github.io/pydicom/stable/tutorials/index.html), +[NIST XCOM attenuation data](https://physics.nist.gov/PhysRefData/Xcom/html/xcom1.html), +[ICRP Publication 103](https://www.icrp.org/publication.asp?id=ICRP%20Publication%20103), +[The Cancer Imaging Archive](https://www.cancerimagingarchive.net/) (free +de-identified datasets). + ## Tests The library methods are pure functions with known textbook answers, which makes @@ -342,6 +404,7 @@ Add any of them with `uv add `. first-principles-code/ ├── main.py # displacement demo (start here) ├── gui.py # MATLAB-style console GUI (uv run gui.py) +├── radiology_demo.py # medical-imaging tour (uv run radiology_demo.py) ├── library/ # one class per engineering domain │ ├── __init__.py # re-exports every class │ ├── units.py # shared Pint unit registry (ureg, Q_) @@ -356,7 +419,13 @@ first-principles-code/ │ ├── gear.py │ ├── electrical.py │ ├── mechatronics.py -│ └── controls.py +│ ├── controls.py +│ ├── radiation.py # radiology: physics, decay, dose +│ ├── radiography.py # radiology: x-ray technique & geometry +│ ├── ct.py # radiology: Hounsfield units & windowing +│ ├── mri.py # radiology: relaxation & contrast +│ ├── ultrasound.py # radiology: pulse-echo & Doppler +│ └── dicom.py # radiology: DICOM I/O with pydicom ├── tests/ # pytest suite (uv run pytest) ├── pyproject.toml # project metadata & dependencies ├── SPEC.md # full project design spec diff --git a/SPEC.md b/SPEC.md index dbf15fd..651a918 100644 --- a/SPEC.md +++ b/SPEC.md @@ -36,6 +36,8 @@ Initial class list: - a class for classical vibration and modal calculations - controls.py - a class for P, PD, PID, and other control technology calculations +- Radiology / medical imaging (added): radiation.py, radiography.py, ct.py, + mri.py, ultrasound.py, dicom.py (uses pydicom) - "please suggest other classes" ## Convntions diff --git a/gui.py b/gui.py index 92574a7..a4a9b83 100644 --- a/gui.py +++ b/gui.py @@ -38,6 +38,7 @@ from library import ( Dynamics, Statics, Modal, Thermo, Fluids, Material, Gear, Electrical, Mechatronics, Controls, Utility, + Radiation, Radiography, CT, MRI, Ultrasound, Dicom, ) @@ -66,6 +67,13 @@ def __init__(self): "electrical": Electrical(), "mechatronics": Mechatronics(), "controls": Controls(), + # Radiology / medical imaging. + "radiation": Radiation(), + "radiography": Radiography(), + "ct": CT(), + "mri": MRI(), + "ultrasound": Ultrasound(), + "dicom": Dicom(), } # Recorded (label, command, result-as-text) for the workspace + export. self.records = [] diff --git a/library/__init__.py b/library/__init__.py index dc10dfa..3161c79 100644 --- a/library/__init__.py +++ b/library/__init__.py @@ -34,6 +34,14 @@ from library.controls import Controls from library.utility import Utility +# Medical imaging / radiology classes. +from library.radiation import Radiation +from library.radiography import Radiography +from library.ct import CT +from library.mri import MRI +from library.ultrasound import Ultrasound +from library.dicom import Dicom + # Unit-aware modelling: the shared Pint registry and the Beam dataclass. from library.units import ureg, Q_ from library.beam import Beam @@ -50,6 +58,12 @@ "Mechatronics", "Controls", "Utility", + "Radiation", + "Radiography", + "CT", + "MRI", + "Ultrasound", + "Dicom", "Beam", "ureg", "Q_", diff --git a/library/ct.py b/library/ct.py new file mode 100644 index 0000000..8787949 --- /dev/null +++ b/library/ct.py @@ -0,0 +1,185 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2026 Evan Gress + +"""ct.py - computed tomography: Hounsfield units, windowing, and dose metrics. + +A CT scanner measures linear attenuation coefficients; the Hounsfield scale +re-expresses them relative to water so that water = 0 HU and air = -1000 HU. +A CT image spans ~4000 HU, but a monitor shows ~256 grey levels and the eye +far fewer, so radiologists choose a WINDOW (a slice of the HU range) to +display. This module does both steps with NumPy, so it works on a single +number or on a whole 512 x 512 image in one call (vectorization). +""" + +import numpy as np + + +class CT: + """Computed-tomography image and dose calculations.""" + + # Common window presets as (width, level) in HU. Values follow widely + # taught defaults (e.g. Radiopaedia); departments tweak these locally. + WINDOW_PRESETS = { + "brain": (80, 40), + "subdural": (200, 75), + "stroke": (40, 40), + "lung": (1500, -600), + "mediastinum": (350, 50), + "abdomen": (400, 50), + "liver": (150, 30), + "bone": (1800, 400), + } + + # Coarse teaching ranges as (upper bound in HU, label), checked in order. + # Real tissues overlap, and contrast, beam hardening, and partial-volume + # effects all shift values - treat these as a study aid, not a diagnosis. + HU_REFERENCE_RANGES = [ + (-900, "air"), + (-200, "lung / aerated tissue"), + (-20, "fat"), + (20, "water / simple fluid"), + (100, "soft tissue / blood"), + (300, "contrast-enhanced tissue"), + (float("inf"), "bone / calcification / metal"), + ] + + def hounsfield_units(self, *, attenuation_coefficient, water_attenuation, air_attenuation=0.0): + """Convert a linear attenuation coefficient to Hounsfield units. + + HU = 1000 * (mu - mu_water) / (mu_water - mu_air) + + Works element-wise on NumPy arrays as well as single numbers. + + :param attenuation_coefficient: tissue mu [1/cm] (number or array) + :param water_attenuation: mu of water at the same beam energy [1/cm] + :param air_attenuation: mu of air [1/cm]; ~0 so it defaults to 0 + :returns: CT number [HU] + """ + mu = np.asarray(attenuation_coefficient, dtype=float) + return 1000.0 * (mu - water_attenuation) / (water_attenuation - air_attenuation) + + def rescale_to_hounsfield(self, *, stored_values, rescale_slope=1.0, rescale_intercept=-1024.0): + """Convert raw stored pixel values to HU: HU = value * slope + intercept. + + DICOM files usually store unsigned integers; the RescaleSlope and + RescaleIntercept tags (0028,1053 / 0028,1052) say how to recover HU. + Forgetting this step is a classic beginner bug: every value is off by + ~1024. See :class:`library.dicom.Dicom` for reading these from a file. + + :param stored_values: raw pixel values (number or array) + :param rescale_slope: DICOM RescaleSlope + :param rescale_intercept: DICOM RescaleIntercept + :returns: values in HU (float array) + """ + return np.asarray(stored_values, dtype=float) * rescale_slope + rescale_intercept + + def window_bounds(self, *, window_width=None, window_level=None, preset=None): + """Lowest and highest HU shown by a window. + + Anything below ``level - width/2`` displays black, anything above + ``level + width/2`` displays white. Pass either width AND level, or a + ``preset`` name from :attr:`WINDOW_PRESETS`. + + :param window_width: window width W [HU] + :param window_level: window level / centre L [HU] + :param preset: optional preset name, e.g. "lung" + :returns: tuple (lower HU, upper HU) + """ + width, level = self._resolve_window(window_width, window_level, preset) + return level - width / 2.0, level + width / 2.0 + + def apply_window(self, *, hounsfield, window_width=None, window_level=None, preset=None, + output_max=255): + """Map HU values to display grey levels for a chosen window. + + Values are clipped to the window, then scaled linearly to 0..output_max. + A narrow window (brain, W=80) spreads a small HU range over all grey + levels -> high contrast; a wide window (bone, W=1800) shows everything + but with low contrast. + + :param hounsfield: HU value or array (e.g. a full CT slice) + :param window_width: window width W [HU] + :param window_level: window level L [HU] + :param preset: optional preset name instead of width/level + :param output_max: brightest grey level (255 for 8-bit display) + :returns: numpy array of grey levels (uint8 when output_max == 255) + """ + lower, upper = self.window_bounds( + window_width=window_width, window_level=window_level, preset=preset, + ) + hu = np.asarray(hounsfield, dtype=float) + # np.clip and the arithmetic below run over the whole image at once - + # no Python loop over pixels needed. + scaled = (np.clip(hu, lower, upper) - lower) / (upper - lower) * output_max + if output_max == 255: + return np.round(scaled).astype(np.uint8) + return scaled + + def classify_hounsfield(self, *, hounsfield): + """Return a rough tissue label for a single HU value (teaching aid). + + :param hounsfield: CT number [HU] + :returns: label string from :attr:`HU_REFERENCE_RANGES` + """ + for upper_bound, label in self.HU_REFERENCE_RANGES: + if hounsfield < upper_bound: + return label + return self.HU_REFERENCE_RANGES[-1][1] + + def dose_length_product(self, *, ctdi_vol_mgy, scan_length_cm): + """Dose-length product DLP = CTDIvol * scan length [mGy*cm]. + + CTDIvol is the average dose within one slice of a standard phantom; + DLP multiplies by how much of the patient was scanned. Both appear on + every scanner's dose report. + + :param ctdi_vol_mgy: volume CT dose index [mGy] + :param scan_length_cm: irradiated scan length [cm] + :returns: DLP [mGy*cm] + """ + return ctdi_vol_mgy * scan_length_cm + + # Adult conversion coefficients k [mSv / (mGy*cm)] from AAPM Report 96 / + # the European Guidelines (EUR 16262). They turn DLP into an estimate of + # effective dose for a "standard" adult. + DLP_TO_EFFECTIVE_DOSE = { + "head": 0.0021, + "neck": 0.0059, + "chest": 0.014, + "abdomen_pelvis": 0.015, + "chest_abdomen_pelvis": 0.015, + } + + def effective_dose_from_dlp(self, *, dlp_mgy_cm, region): + """Estimate adult effective dose: E ~= k * DLP. + + Example: a routine head CT with DLP ~ 1000 mGy*cm gives ~2 mSv. + + :param dlp_mgy_cm: dose-length product [mGy*cm] + :param region: key of :attr:`DLP_TO_EFFECTIVE_DOSE` + :returns: effective dose estimate [mSv] + """ + if region not in self.DLP_TO_EFFECTIVE_DOSE: + raise ValueError(f"unknown region {region!r}; use one of: " + f"{', '.join(self.DLP_TO_EFFECTIVE_DOSE)}") + return self.DLP_TO_EFFECTIVE_DOSE[region] * dlp_mgy_cm + + # ------------------------------------------------------------------ + # Helpers + # ------------------------------------------------------------------ + def _resolve_window(self, width, level, preset): + """Return (width, level) from either a preset name or explicit values. + + The leading underscore marks this as an internal helper - a Python + convention meaning "not part of the public API". + """ + if preset is not None: + if preset not in self.WINDOW_PRESETS: + raise ValueError(f"unknown preset {preset!r}; use one of: " + f"{', '.join(self.WINDOW_PRESETS)}") + return self.WINDOW_PRESETS[preset] + if width is None or level is None: + raise ValueError("give window_width and window_level, or a preset") + if width <= 0: + raise ValueError("window_width must be positive") + return width, level diff --git a/library/dicom.py b/library/dicom.py new file mode 100644 index 0000000..cd4e8c1 --- /dev/null +++ b/library/dicom.py @@ -0,0 +1,160 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2026 Evan Gress + +"""dicom.py - reading and handling DICOM medical images with pydicom. + +DICOM is the file format (and network protocol) used by every CT, MR, x-ray, +and ultrasound system. A DICOM file holds two things: + +* a header of tagged metadata - modality, pixel spacing, rescale values, and + patient identifiers (Protected Health Information, PHI); +* the pixel data, which pydicom hands back as a NumPy array. + +This class wraps the handful of pydicom calls a student needs: load a file, +get the image in real units (HU for CT), summarize the header, and strip PHI +before sharing an image. pydicom ships sample files, so you can practise +without any patient data:: + + from pydicom.data import get_testdata_file + path = get_testdata_file("CT_small.dcm") +""" + +import copy + +import numpy as np +import pydicom + + +class Dicom: + """Load, inspect, convert, and de-identify DICOM images.""" + + # Header keywords copied into :meth:`summary`. Only technical fields - no + # patient identifiers - so the summary is safe to print or log. + SUMMARY_KEYWORDS = ( + "Modality", + "BodyPartExamined", + "Manufacturer", + "Rows", + "Columns", + "PixelSpacing", + "SliceThickness", + "KVP", + "XRayTubeCurrent", + "ExposureTime", + "RescaleSlope", + "RescaleIntercept", + "WindowCenter", + "WindowWidth", + ) + + # Direct identifiers to blank when de-identifying. This is a teaching + # subset: full de-identification follows DICOM PS3.15 Annex E and also + # needs dates, UIDs, private tags, and text burned into the pixels handled. + PHI_KEYWORDS = ( + "PatientName", + "PatientID", + "PatientBirthDate", + "PatientAddress", + "PatientTelephoneNumbers", + "OtherPatientIDs", + "ReferringPhysicianName", + "InstitutionName", + "InstitutionAddress", + "AccessionNumber", + ) + + def read(self, *, path): + """Read a DICOM file from disk. + + :param path: file path (str or pathlib.Path) + :returns: a ``pydicom.Dataset`` - access tags as attributes, + e.g. ``ds.Modality``, and pixels as ``ds.pixel_array`` + """ + return pydicom.dcmread(path) + + def summary(self, *, dataset): + """Collect the key technical header fields into a plain dict. + + Uses ``dataset.get(keyword)`` rather than ``dataset.keyword`` so a + missing tag returns None instead of raising - many tags are optional + and vary by modality. + + :param dataset: a pydicom Dataset + :returns: dict {keyword: value or None} + """ + # A dict comprehension: build the whole dictionary in one expression. + return {keyword: dataset.get(keyword) for keyword in self.SUMMARY_KEYWORDS} + + def pixel_array_in_units(self, *, dataset): + """Pixel data converted to real-world values (HU for CT). + + Applies value * RescaleSlope + RescaleIntercept. If the tags are + absent (common for plain radiographs) the stored values are returned + unchanged as floats. + + :param dataset: a pydicom Dataset containing pixel data + :returns: float numpy array + """ + slope = float(dataset.get("RescaleSlope", 1.0)) + intercept = float(dataset.get("RescaleIntercept", 0.0)) + return dataset.pixel_array.astype(float) * slope + intercept + + def pixel_spacing_mm(self, *, dataset): + """Physical size of one pixel as (row spacing, column spacing) in mm. + + Needed for any real measurement: a length in pixels times the spacing + gives millimetres. + + :param dataset: a pydicom Dataset + :returns: tuple (row_mm, col_mm), or None if the tag is missing + """ + spacing = dataset.get("PixelSpacing") + if spacing is None: + return None + return float(spacing[0]), float(spacing[1]) + + def anonymize(self, *, dataset, replacement="ANONYMIZED"): + """Blank the direct patient identifiers and drop private tags. + + Works on a DEEP copy so the original dataset is left untouched - + mutating a caller's object behind their back is a common source of + bugs. (``dataset.copy()`` is only a shallow copy: the data elements are + shared, so editing one would silently edit the original too.) + + :param dataset: a pydicom Dataset + :param replacement: text placed in each blanked field + :returns: a new, de-identified Dataset + """ + cleaned = copy.deepcopy(dataset) + for keyword in self.PHI_KEYWORDS: + if keyword in cleaned: + # Dates must stay valid DICOM dates, so blank them instead. + is_date = cleaned[keyword].VR == "DA" + setattr(cleaned, keyword, "" if is_date else replacement) + # Vendors hide all sorts of things (sometimes PHI) in private tags. + cleaned.remove_private_tags() + return cleaned + + def save(self, *, dataset, path): + """Write a Dataset back to disk as a DICOM file. + + :param dataset: a pydicom Dataset + :param path: output file path + """ + dataset.save_as(path) + + def to_png(self, *, image, path): + """Save a 2-D greyscale array (e.g. a windowed CT slice) as a PNG. + + Uses matplotlib's ``imsave`` so no extra imaging library is needed. + + :param image: 2-D numpy array + :param path: output .png path + """ + import matplotlib.pyplot as plt # imported here: only needed for export + + image = np.asarray(image) + # For 8-bit display data, pin the grey scale to 0..255 so matplotlib + # does not re-stretch the contrast you chose with your window. + limits = {"vmin": 0, "vmax": 255} if image.dtype == np.uint8 else {} + plt.imsave(path, image, cmap="gray", **limits) diff --git a/library/mri.py b/library/mri.py new file mode 100644 index 0000000..f78352a --- /dev/null +++ b/library/mri.py @@ -0,0 +1,131 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2026 Evan Gress + +"""mri.py - magnetic resonance imaging signal physics. + +MRI contrast comes from how fast hydrogen nuclei (protons) recover after an RF +pulse: T1 is longitudinal recovery, T2 is transverse decay. By choosing the +repetition time (TR) and echo time (TE) the operator decides which of these +differences dominates the image - T1-weighted, T2-weighted, or proton-density +(PD) weighted. These methods let you compute that contrast instead of just +memorizing "fat is bright on T1". +""" + +import math + + +class MRI: + """MRI precession, relaxation, and spin-echo contrast calculations.""" + + # Gyromagnetic ratio over 2*pi, in MHz per tesla. + GYROMAGNETIC_RATIO_MHZ_PER_T = { + "1H": 42.577, + "13C": 10.708, + "19F": 40.078, + "23Na": 11.262, + "31P": 17.235, + } + + # Approximate relaxation times at 1.5 T in milliseconds (T1, T2), plus a + # relative proton density. Literature values vary by +/- 20% or more with + # field strength and measurement method - use them to compare tissues, + # not as reference standards. + TISSUE_PROPERTIES_1_5T = { + "fat": {"t1": 260, "t2": 80, "proton_density": 0.9}, + "white_matter": {"t1": 790, "t2": 90, "proton_density": 0.7}, + "grey_matter": {"t1": 920, "t2": 100, "proton_density": 0.8}, + "csf": {"t1": 4000, "t2": 2000, "proton_density": 1.0}, + "muscle": {"t1": 870, "t2": 50, "proton_density": 0.75}, + "liver": {"t1": 500, "t2": 40, "proton_density": 0.7}, + } + + def larmor_frequency(self, *, field_strength_t, nucleus="1H"): + """Precession (Larmor) frequency f = gamma_bar * B0. + + The RF coil must transmit and receive at exactly this frequency; the + scanner's gradients make it vary with position, which is how MRI knows + WHERE the signal came from. + + :param field_strength_t: main magnetic field B0 [T] + :param nucleus: key of :attr:`GYROMAGNETIC_RATIO_MHZ_PER_T` + :returns: frequency [MHz] (1H at 1.5 T -> ~63.9 MHz) + """ + if nucleus not in self.GYROMAGNETIC_RATIO_MHZ_PER_T: + raise ValueError(f"unknown nucleus {nucleus!r}") + return self.GYROMAGNETIC_RATIO_MHZ_PER_T[nucleus] * field_strength_t + + def t1_recovery(self, *, time_ms, t1_ms, equilibrium_magnetization=1.0): + """Longitudinal recovery after a 90 deg pulse: Mz = M0 * (1 - exp(-t/T1)). + + After one T1, magnetization is back to ~63% of equilibrium. + + :param time_ms: time since the 90 deg pulse [ms] + :param t1_ms: tissue T1 [ms] + :param equilibrium_magnetization: M0 (defaults to 1, i.e. a fraction) + :returns: Mz, same unit as M0 + """ + return equilibrium_magnetization * (1.0 - math.exp(-time_ms / t1_ms)) + + def t2_decay(self, *, time_ms, t2_ms, initial_magnetization=1.0): + """Transverse decay: Mxy = M0 * exp(-t/T2). + + After one T2, only ~37% of the transverse signal remains. + + :param time_ms: time since excitation [ms] + :param t2_ms: tissue T2 [ms] + :param initial_magnetization: Mxy at t = 0 (defaults to 1) + :returns: Mxy, same unit as the initial magnetization + """ + return initial_magnetization * math.exp(-time_ms / t2_ms) + + def spin_echo_signal(self, *, repetition_time_ms, echo_time_ms, t1_ms=None, t2_ms=None, + proton_density=1.0, tissue=None): + """Relative spin-echo signal S = PD * (1 - exp(-TR/T1)) * exp(-TE/T2). + + TR controls how much T1 recovery happens between pulses; TE controls + how much T2 decay happens before the echo is read. Pass explicit + T1/T2/PD, or a ``tissue`` name from :attr:`TISSUE_PROPERTIES_1_5T`. + + :param repetition_time_ms: TR [ms] + :param echo_time_ms: TE [ms] + :param t1_ms: tissue T1 [ms] + :param t2_ms: tissue T2 [ms] + :param proton_density: relative proton density [-] + :param tissue: optional tissue name that supplies T1, T2, and PD + :returns: relative signal intensity [-] + """ + if tissue is not None: + props = self.TISSUE_PROPERTIES_1_5T[tissue] + t1_ms, t2_ms, proton_density = props["t1"], props["t2"], props["proton_density"] + if t1_ms is None or t2_ms is None: + raise ValueError("give t1_ms and t2_ms, or a tissue name") + recovered = 1.0 - math.exp(-repetition_time_ms / t1_ms) + decayed = math.exp(-echo_time_ms / t2_ms) + return proton_density * recovered * decayed + + def sequence_weighting(self, *, repetition_time_ms, echo_time_ms): + """Classify a spin-echo TR/TE pair as T1-, T2-, or PD-weighted. + + Rule of thumb at 1.5 T: + short TR (< 800 ms) + short TE (< 30 ms) -> T1-weighted + long TR (> 2000 ms) + long TE (> 80 ms) -> T2-weighted + long TR + short TE -> PD-weighted + + :param repetition_time_ms: TR [ms] + :param echo_time_ms: TE [ms] + :returns: "T1-weighted", "T2-weighted", "PD-weighted", or "mixed" + """ + short_tr = repetition_time_ms < 800 + long_tr = repetition_time_ms > 2000 + short_te = echo_time_ms < 30 + long_te = echo_time_ms > 80 + + if short_tr and short_te: + return "T1-weighted" + if long_tr and long_te: + return "T2-weighted" + if long_tr and short_te: + return "PD-weighted" + # Short TR + long TE gives little signal and poor contrast; anything + # in the in-between zone mixes several contrast mechanisms. + return "mixed" diff --git a/library/radiation.py b/library/radiation.py new file mode 100644 index 0000000..552102a --- /dev/null +++ b/library/radiation.py @@ -0,0 +1,254 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2026 Evan Gress + +"""radiation.py - radiation physics and radiation protection. + +The physics every radiology student meets in the first semester: how x-ray +intensity falls off with distance, how a beam is attenuated by tissue and +shielding, how a radionuclide decays, and how absorbed dose [Gy] is turned +into the protection quantities equivalent dose and effective dose [Sv]. + +All the protection weighting factors come from ICRP Publication 103 (2007), +which is the basis of current regulations in most countries. + +These are teaching calculations. They are NOT a substitute for a qualified +medical physicist, a validated dosimetry system, or your local regulations. +""" + +import math + + +class Radiation: + """Radiation physics, attenuation, decay, and dosimetry calculations.""" + + # ICRP 103 radiation weighting factors w_R (dimensionless). Neutrons are + # missing on purpose: their w_R depends on energy, so they get a method of + # their own (see :meth:`neutron_weighting_factor`). + RADIATION_WEIGHTING_FACTORS = { + "photon": 1.0, # x-rays and gamma rays + "electron": 1.0, # beta particles, muons + "proton": 2.0, + "alpha": 20.0, # also fission fragments and heavy ions + } + + # ICRP 103 tissue weighting factors w_T. They describe what fraction of the + # whole-body stochastic risk each organ carries, so they sum to exactly 1. + TISSUE_WEIGHTING_FACTORS = { + "red_bone_marrow": 0.12, + "colon": 0.12, + "lung": 0.12, + "stomach": 0.12, + "breast": 0.12, + "remainder": 0.12, # 13 remaining tissues, averaged together + "gonads": 0.08, + "bladder": 0.04, + "oesophagus": 0.04, + "liver": 0.04, + "thyroid": 0.04, + "bone_surface": 0.01, + "brain": 0.01, + "salivary_glands": 0.01, + "skin": 0.01, + } + + # Planck constant times speed of light, expressed in keV * nm. Lets us go + # between photon energy and wavelength without carrying SI exponents. + HC_KEV_NM = 1.239842 + + # ------------------------------------------------------------------ + # Distance + # ------------------------------------------------------------------ + def inverse_square_law(self, *, initial_intensity, initial_distance, new_distance): + """Intensity at a new distance from a point source. + + Radiation spreads over a sphere whose area grows with r^2, so intensity + falls as 1/r^2: I2 = I1 * (d1 / d2)^2. Doubling your distance from the + tube cuts your exposure to one quarter - the cheapest protection there is. + + :param initial_intensity: intensity (or exposure / dose rate) at d1, any unit + :param initial_distance: original distance d1 (any length unit) + :param new_distance: new distance d2 (same unit as d1) + :returns: intensity at d2, in the same unit as ``initial_intensity`` + """ + if initial_distance <= 0 or new_distance <= 0: + raise ValueError("distances must be positive") + return initial_intensity * (initial_distance / new_distance) ** 2 + + # ------------------------------------------------------------------ + # Attenuation and shielding + # ------------------------------------------------------------------ + def transmitted_intensity(self, *, initial_intensity, attenuation_coefficient, thickness): + """Narrow-beam (Beer-Lambert) attenuation: I = I0 * exp(-mu * x). + + Each thin slab of material removes the same FRACTION of the photons + that reach it, which is why the curve is exponential, not linear. + + :param initial_intensity: incident intensity I0 (any unit) + :param attenuation_coefficient: linear attenuation coefficient mu [1/cm] + :param thickness: material thickness x [cm] + :returns: transmitted intensity, same unit as I0 + """ + return initial_intensity * math.exp(-attenuation_coefficient * thickness) + + def linear_attenuation_coefficient(self, *, mass_attenuation_coefficient, density): + """Linear coefficient mu [1/cm] from mass coefficient mu/rho [cm^2/g]. + + Tables (e.g. NIST XCOM) list mu/rho because it does not depend on the + material's physical state; multiply by density to get mu. + + :param mass_attenuation_coefficient: mu/rho [cm^2/g] + :param density: rho [g/cm^3] + :returns: mu [1/cm] + """ + return mass_attenuation_coefficient * density + + def half_value_layer(self, *, attenuation_coefficient): + """Half-value layer HVL = ln(2) / mu. + + The thickness that cuts the beam intensity in half. Clinically the HVL + (in mm of aluminium) is how beam QUALITY is specified. + + :param attenuation_coefficient: mu [1/cm] + :returns: HVL [cm] + """ + return math.log(2) / attenuation_coefficient + + def tenth_value_layer(self, *, attenuation_coefficient): + """Tenth-value layer TVL = ln(10) / mu (~3.32 HVLs). + + Shielding designers prefer TVLs: each one cuts intensity by 10x. + + :param attenuation_coefficient: mu [1/cm] + :returns: TVL [cm] + """ + return math.log(10) / attenuation_coefficient + + def shielding_thickness(self, *, attenuation_coefficient, transmission): + """Thickness needed to let only a given fraction of the beam through. + + Solve I/I0 = exp(-mu*x) for x: x = -ln(T) / mu. + + Note: real shielding has scatter "build-up", so true broad-beam barriers + must be thicker than this narrow-beam estimate (see NCRP Report 147). + + :param attenuation_coefficient: mu [1/cm] + :param transmission: desired transmitted fraction T, 0 < T <= 1 + :returns: required thickness [cm] + """ + if not 0 < transmission <= 1: + raise ValueError("transmission must be in (0, 1]") + return -math.log(transmission) / attenuation_coefficient + + # ------------------------------------------------------------------ + # Photon energy + # ------------------------------------------------------------------ + def photon_wavelength(self, *, energy_kev): + """Photon wavelength from energy: lambda = h*c / E. + + Applied to the tube voltage this gives the Duane-Hunt minimum + wavelength: a 100 kVp beam contains no photon above 100 keV, so no + wavelength shorter than ~0.0124 nm. + + :param energy_kev: photon energy [keV] + :returns: wavelength [nm] + """ + return self.HC_KEV_NM / energy_kev + + # ------------------------------------------------------------------ + # Radioactive decay (nuclear medicine) + # ------------------------------------------------------------------ + def decay_constant(self, *, half_life): + """Decay constant lambda = ln(2) / T_half. + + :param half_life: physical half-life (any time unit) + :returns: decay constant [1 / that time unit] + """ + return math.log(2) / half_life + + def remaining_activity(self, *, initial_activity, half_life, elapsed_time): + """Activity left after decay: A = A0 * exp(-lambda * t). + + Used daily in nuclear medicine to work out how much of a dose is left + at injection time (e.g. Tc-99m, T_half = 6.01 h). + + :param initial_activity: A0 (MBq, mCi, ...) + :param half_life: half-life (any time unit) + :param elapsed_time: elapsed time (same unit as half_life) + :returns: remaining activity, same unit as A0 + """ + decay = self.decay_constant(half_life=half_life) + return initial_activity * math.exp(-decay * elapsed_time) + + def effective_half_life(self, *, physical_half_life, biological_half_life): + """Effective half-life in the body: 1/Te = 1/Tp + 1/Tb. + + A radiopharmaceutical leaves the patient two ways at once - it decays + (physical) AND is excreted (biological) - so the rates add, and the + effective half-life is always shorter than either one. + + :param physical_half_life: Tp (any time unit) + :param biological_half_life: Tb (same unit) + :returns: effective half-life Te (same unit) + """ + return (physical_half_life * biological_half_life) / ( + physical_half_life + biological_half_life + ) + + # ------------------------------------------------------------------ + # Dosimetry: Gy -> Sv + # ------------------------------------------------------------------ + def neutron_weighting_factor(self, *, energy_mev): + """ICRP 103 continuous radiation weighting factor for neutrons. + + :param energy_mev: neutron energy [MeV] + :returns: w_R [-] (peaks at ~20 around 1 MeV) + """ + if energy_mev <= 0: + raise ValueError("neutron energy must be positive") + # Three branches straight from ICRP 103, eq. (4.3). + if energy_mev < 1.0: + return 2.5 + 18.2 * math.exp(-math.log(energy_mev) ** 2 / 6) + if energy_mev <= 50.0: + return 5.0 + 17.0 * math.exp(-math.log(2 * energy_mev) ** 2 / 6) + return 2.5 + 3.25 * math.exp(-math.log(0.04 * energy_mev) ** 2 / 6) + + def equivalent_dose(self, *, absorbed_dose, radiation_type): + """Equivalent dose H_T = w_R * D_T. + + Absorbed dose (Gy = J/kg) says how much energy was deposited; the + weighting factor accounts for how biologically damaging that type of + radiation is. 1 Gy of alpha harms tissue ~20x more than 1 Gy of x-rays. + + :param absorbed_dose: absorbed dose D [Gy] + :param radiation_type: key of :attr:`RADIATION_WEIGHTING_FACTORS` + :returns: equivalent dose [Sv] + """ + try: + weight = self.RADIATION_WEIGHTING_FACTORS[radiation_type] + except KeyError: + valid = ", ".join(self.RADIATION_WEIGHTING_FACTORS) + raise ValueError(f"unknown radiation_type {radiation_type!r}; use one of: {valid}") + return weight * absorbed_dose + + def effective_dose(self, *, organ_equivalent_doses): + """Effective dose E = sum over organs of w_T * H_T. + + Lets you compare a chest x-ray with a head CT on one risk-related + scale. It is a POPULATION protection quantity - not an individual + patient's risk estimate. + + Organs you leave out contribute nothing, so this naturally handles a + partial-body exposure. + + :param organ_equivalent_doses: dict {organ: H_T [Sv]}, organs named as + in :attr:`TISSUE_WEIGHTING_FACTORS` + :returns: effective dose [Sv] + """ + unknown = set(organ_equivalent_doses) - set(self.TISSUE_WEIGHTING_FACTORS) + if unknown: + raise ValueError(f"unknown organ(s): {sorted(unknown)}") + # A generator expression inside sum(): weight each organ, then add. + return sum( + self.TISSUE_WEIGHTING_FACTORS[organ] * dose + for organ, dose in organ_equivalent_doses.items() + ) diff --git a/library/radiography.py b/library/radiography.py new file mode 100644 index 0000000..bc28ff0 --- /dev/null +++ b/library/radiography.py @@ -0,0 +1,195 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2026 Evan Gress + +"""radiography.py - projection x-ray technique and image geometry. + +The "technique chart" arithmetic a radiographer does at the console: setting +mAs, trading kVp for mAs with the 15% rule, compensating for a change in +distance or grid, and predicting magnification and blur from the geometry of +tube, patient, and image receptor (IR). + +Abbreviations used throughout: + SID source-to-image distance (tube focal spot -> detector) + SOD source-to-object distance (focal spot -> anatomy) + OID object-to-image distance (anatomy -> detector); SID = SOD + OID +""" + +import math + + +class Radiography: + """Exposure-technique and geometric calculations for projection radiography.""" + + # Grid conversion factors (Bucky factors): how much the mAs must rise to + # keep the same IR exposure once a grid is soaking up scatter. + # Values are the standard textbook set (Bushong / Carlton & Adler). + GRID_CONVERSION_FACTORS = { + "none": 1, + "5:1": 2, + "6:1": 3, + "8:1": 4, + "12:1": 5, + "16:1": 6, + } + + # Heat-unit rectification constants per generator type. Texts differ + # slightly for high-frequency units (1.40 - 1.45); 1.45 follows Bushong. + GENERATOR_HEAT_FACTORS = { + "single_phase": 1.00, + "three_phase_6_pulse": 1.35, + "three_phase_12_pulse": 1.41, + "high_frequency": 1.45, + } + + # ------------------------------------------------------------------ + # Exposure factors + # ------------------------------------------------------------------ + def mas(self, *, tube_current_ma, exposure_time_s): + """Tube current-time product mAs = mA * s. + + mAs sets the QUANTITY of x-ray photons. Reciprocity says 200 mA x 0.1 s + and 400 mA x 0.05 s give the same 20 mAs - but the shorter time freezes + patient motion better. + + :param tube_current_ma: tube current [mA] + :param exposure_time_s: exposure time [s] + :returns: mAs [mA*s] + """ + return tube_current_ma * exposure_time_s + + def fifteen_percent_rule(self, *, kvp, mas, steps=1): + """Trade kVp for mAs while keeping the same exposure to the IR. + + Raising kVp by 15% roughly doubles the exposure reaching the detector, + so mAs can be halved. Each "step" applies the rule once; negative steps + go the other way (lower kVp 15%, double mAs). Raising kVp this way + lowers patient dose but also lowers subject contrast. + + :param kvp: starting peak tube voltage [kVp] + :param mas: starting mAs + :param steps: number of 15% adjustments (+ raises kVp, - lowers it) + :returns: tuple (new kVp, new mAs) + """ + new_kvp = kvp * 1.15 ** steps + new_mas = mas / 2 ** steps + return new_kvp, new_mas + + def exposure_maintenance_mas(self, *, original_mas, original_sid, new_sid): + """Density/exposure maintenance formula: mAs2 = mAs1 * (SID2 / SID1)^2. + + The inverse square law in reverse: moving the tube further away spreads + the beam out, so the mAs must rise by the square of the distance ratio + to keep the same detector exposure. + + :param original_mas: mAs used at the original SID + :param original_sid: original source-to-image distance (any length unit) + :param new_sid: new SID (same unit) + :returns: new mAs + """ + return original_mas * (new_sid / original_sid) ** 2 + + def grid_mas(self, *, original_mas, original_grid, new_grid): + """mAs compensation when adding, removing, or changing a grid. + + mAs2 = mAs1 * GCF2 / GCF1, with GCF from :attr:`GRID_CONVERSION_FACTORS`. + + :param original_mas: mAs used with the original grid + :param original_grid: grid ratio key, e.g. "none" or "12:1" + :param new_grid: grid ratio key for the new setup + :returns: new mAs + """ + factors = self.GRID_CONVERSION_FACTORS + for grid in (original_grid, new_grid): + if grid not in factors: + raise ValueError(f"unknown grid {grid!r}; use one of: {', '.join(factors)}") + return original_mas * factors[new_grid] / factors[original_grid] + + def heat_units(self, *, kvp, tube_current_ma, exposure_time_s, generator="high_frequency", + exposures=1): + """Anode heat load in heat units: HU = kVp * mA * s * C * exposures. + + Used with the tube rating chart to check a series of exposures (e.g. a + fluoro run or CT scan) will not overheat the anode. + + :param kvp: peak tube voltage [kVp] + :param tube_current_ma: tube current [mA] + :param exposure_time_s: time per exposure [s] + :param generator: key of :attr:`GENERATOR_HEAT_FACTORS` + :param exposures: number of exposures in the series + :returns: heat units [HU] + """ + if generator not in self.GENERATOR_HEAT_FACTORS: + raise ValueError(f"unknown generator {generator!r}") + factor = self.GENERATOR_HEAT_FACTORS[generator] + return kvp * tube_current_ma * exposure_time_s * factor * exposures + + # ------------------------------------------------------------------ + # Geometry + # ------------------------------------------------------------------ + def magnification_factor(self, *, sid, sod): + """Magnification factor MF = SID / SOD (always >= 1). + + Anatomy far from the detector is projected larger - which is why a PA + chest makes the heart look smaller than an AP chest. + + :param sid: source-to-image distance (any length unit) + :param sod: source-to-object distance (same unit) + :returns: magnification factor [-] + """ + if sod <= 0 or sod > sid: + raise ValueError("need 0 < SOD <= SID (the object sits between tube and detector)") + return sid / sod + + def object_size(self, *, image_size, sid, sod): + """True anatomical size from a measured image: object = image / MF. + + :param image_size: size measured on the image (any length unit) + :param sid: source-to-image distance + :param sod: source-to-object distance + :returns: true object size, same unit as ``image_size`` + """ + return image_size / self.magnification_factor(sid=sid, sod=sod) + + def geometric_unsharpness(self, *, focal_spot_size, sid, oid): + """Penumbra (edge blur) width: Ug = focal spot * OID / SOD. + + A real focal spot is not a point, so every edge gets a blurred fringe. + Minimize it with a small focal spot, long SID, and anatomy close to the + detector (small OID). + + :param focal_spot_size: effective focal spot size [mm] + :param sid: source-to-image distance (any length unit) + :param oid: object-to-image distance (same unit as SID) + :returns: geometric unsharpness [mm] + """ + sod = sid - oid + if sod <= 0: + raise ValueError("OID must be smaller than SID") + return focal_spot_size * oid / sod + + # ------------------------------------------------------------------ + # Digital detectors + # ------------------------------------------------------------------ + def nyquist_limit(self, *, pixel_pitch_mm): + """Highest resolvable spatial frequency of a digital detector. + + One line pair (a dark + a bright line) needs at least two pixels, so + f_Nyquist = 1 / (2 * pixel pitch). Finer detail than this aliases. + + :param pixel_pitch_mm: detector element spacing [mm] + :returns: limiting spatial resolution [line pairs / mm] + """ + return 1.0 / (2.0 * pixel_pitch_mm) + + def exposure_index_deviation(self, *, exposure_index, target_exposure_index): + """Deviation index DI = 10 * log10(EI / EI_target) (IEC 62494-1). + + A vendor-neutral check on whether a digital image was over- or + under-exposed. DI = 0 is on target; +1 is ~26% too much, +3 is double. + Rule of thumb: |DI| <= 1 is acceptable; DI > +3 is a dose problem. + + :param exposure_index: EI reported by the detector + :param target_exposure_index: EI_T set for this exam/body part + :returns: deviation index [-] + """ + return 10.0 * math.log10(exposure_index / target_exposure_index) diff --git a/library/ultrasound.py b/library/ultrasound.py new file mode 100644 index 0000000..9e26eb1 --- /dev/null +++ b/library/ultrasound.py @@ -0,0 +1,159 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2026 Evan Gress + +"""ultrasound.py - diagnostic ultrasound (sonography) physics. + +An ultrasound scanner sends a short pulse into the body and times the echoes. +It assumes sound always travels at 1540 m/s in soft tissue - that single +assumption explains depth calculation, many artefacts, and the trade-off +between frequency, resolution, and penetration. Doppler adds blood-flow +velocity from the frequency shift of moving red cells. +""" + +import math + + +class Ultrasound: + """Pulse-echo, attenuation, and Doppler calculations for sonography.""" + + # Speed of sound every scanner assumes for soft tissue [m/s]. + SOFT_TISSUE_SPEED = 1540.0 + + # Typical acoustic properties: density [kg/m^3] and speed of sound [m/s]. + MEDIA = { + "air": {"density": 1.2, "speed": 330.0}, + "fat": {"density": 925.0, "speed": 1450.0}, + "water": {"density": 1000.0, "speed": 1480.0}, + "soft_tissue": {"density": 1050.0, "speed": 1540.0}, + "liver": {"density": 1060.0, "speed": 1550.0}, + "muscle": {"density": 1075.0, "speed": 1580.0}, + "bone": {"density": 1900.0, "speed": 4080.0}, + } + + def wavelength(self, *, frequency_mhz, speed=SOFT_TISSUE_SPEED): + """Wavelength lambda = c / f. + + :param frequency_mhz: transducer frequency [MHz] + :param speed: speed of sound [m/s] (default: soft tissue) + :returns: wavelength [mm] (5 MHz in soft tissue -> 0.308 mm) + """ + # Work in SI (m/s divided by Hz gives metres), then convert m -> mm. + return speed / (frequency_mhz * 1e6) * 1000.0 + + def echo_depth(self, *, round_trip_time_us, speed=SOFT_TISSUE_SPEED): + """Reflector depth from echo time: d = c * t / 2. + + Divide by two because the pulse travels down AND back. This "13 us per + cm" rule is hard-wired into every scanner. + + :param round_trip_time_us: time from pulse to echo [microseconds] + :param speed: speed of sound [m/s] + :returns: depth [cm] + """ + return speed * round_trip_time_us * 1e-6 / 2.0 * 100.0 + + def acoustic_impedance(self, *, density, speed): + """Acoustic impedance Z = rho * c [rayl = kg/(m^2*s)]. + + :param density: density [kg/m^3] + :param speed: speed of sound [m/s] + :returns: impedance [rayl] + """ + return density * speed + + def intensity_reflection_coefficient(self, *, impedance_1, impedance_2): + """Fraction of intensity reflected at a boundary: R = ((Z2-Z1)/(Z2+Z1))^2. + + Soft-tissue/soft-tissue boundaries reflect < 1% (good: the beam carries + on to deeper structures). Tissue/air reflects ~99.9% - which is why we + need coupling gel and why bowel gas blocks the view. + + :param impedance_1: Z of the first medium [rayl] + :param impedance_2: Z of the second medium [rayl] + :returns: reflected fraction R [0..1] + """ + return ((impedance_2 - impedance_1) / (impedance_2 + impedance_1)) ** 2 + + def reflection_between(self, *, medium_1, medium_2): + """Reflection coefficient between two named media from :attr:`MEDIA`. + + A convenience wrapper that chains the two methods above - a small + example of building higher-level helpers from simple ones. + + :param medium_1: name of the first medium, e.g. "soft_tissue" + :param medium_2: name of the second medium, e.g. "air" + :returns: reflected intensity fraction R [0..1] + """ + z1 = self.acoustic_impedance(**self.MEDIA[medium_1]) + z2 = self.acoustic_impedance(**self.MEDIA[medium_2]) + return self.intensity_reflection_coefficient(impedance_1=z1, impedance_2=z2) + + def attenuation_db(self, *, frequency_mhz, path_length_cm, coefficient_db_per_cm_mhz=0.5): + """Total attenuation A = alpha * f * L [dB]. + + Soft tissue loses ~0.5 dB per cm per MHz. Higher frequency gives better + resolution but less penetration - the central trade-off when choosing + a probe. For a pulse-echo, the path length is TWICE the depth. + + :param frequency_mhz: frequency [MHz] + :param path_length_cm: distance travelled [cm] + :param coefficient_db_per_cm_mhz: attenuation coefficient [dB/(cm*MHz)] + :returns: attenuation [dB] + """ + return coefficient_db_per_cm_mhz * frequency_mhz * path_length_cm + + def axial_resolution(self, *, frequency_mhz, cycles_per_pulse=2, speed=SOFT_TISSUE_SPEED): + """Axial resolution = spatial pulse length / 2 = n * lambda / 2. + + Two reflectors along the beam closer than this blur into one echo. + + :param frequency_mhz: frequency [MHz] + :param cycles_per_pulse: number of cycles in each pulse + :param speed: speed of sound [m/s] + :returns: axial resolution [mm] + """ + return cycles_per_pulse * self.wavelength(frequency_mhz=frequency_mhz, speed=speed) / 2.0 + + def maximum_depth(self, *, pulse_repetition_frequency_hz, speed=SOFT_TISSUE_SPEED): + """Deepest unambiguous echo for a given PRF: d_max = c / (2 * PRF). + + Each echo must return before the next pulse is fired, so imaging deeper + forces a lower PRF (and lower frame rate). + + :param pulse_repetition_frequency_hz: PRF [Hz] + :param speed: speed of sound [m/s] + :returns: maximum depth [cm] + """ + return speed / (2.0 * pulse_repetition_frequency_hz) * 100.0 + + def doppler_shift(self, *, transmitted_frequency_mhz, velocity, angle_deg, + speed=SOFT_TISSUE_SPEED): + """Doppler frequency shift df = 2 * f0 * v * cos(theta) / c. + + :param transmitted_frequency_mhz: transmitted frequency f0 [MHz] + :param velocity: reflector (blood) speed [m/s] + :param angle_deg: angle between beam and flow [degrees] + :param speed: speed of sound [m/s] + :returns: Doppler shift [Hz] (conveniently in the audible range) + """ + f0_hz = transmitted_frequency_mhz * 1e6 + return 2.0 * f0_hz * velocity * math.cos(math.radians(angle_deg)) / speed + + def doppler_velocity(self, *, transmitted_frequency_mhz, doppler_shift_hz, angle_deg, + speed=SOFT_TISSUE_SPEED): + """Blood velocity from a measured Doppler shift (the inverse problem). + + v = df * c / (2 * f0 * cos(theta)). Because cos(theta) sits in the + denominator, angle errors explode near 90 deg - keep the insonation + angle at or below 60 deg in clinical practice. + + :param transmitted_frequency_mhz: f0 [MHz] + :param doppler_shift_hz: measured shift [Hz] + :param angle_deg: insonation angle [degrees] + :param speed: speed of sound [m/s] + :returns: velocity [m/s] + """ + cosine = math.cos(math.radians(angle_deg)) + if abs(cosine) < 1e-6: + raise ValueError("at 90 degrees there is no Doppler shift to measure") + return doppler_shift_hz * speed / (2.0 * transmitted_frequency_mhz * 1e6 * cosine) diff --git a/pyproject.toml b/pyproject.toml index 2cbd735..d4a2789 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,6 +11,7 @@ dependencies = [ "matplotlib>=3.10.9", "numpy>=2.4.6", "pint>=0.25.3", + "pydicom>=3.0.2", "scipy>=1.17.1", "scipy-stubs~=1.17.1", ] diff --git a/radiology_demo.py b/radiology_demo.py new file mode 100644 index 0000000..e1e4fd3 --- /dev/null +++ b/radiology_demo.py @@ -0,0 +1,130 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2026 Evan Gress + +"""radiology_demo.py - a worked tour of the medical-imaging classes. + +Run it with: uv run radiology_demo.py + +It walks through one calculation per modality, then loads a real (sample, +anonymous) CT slice with pydicom, converts it to Hounsfield units, and shows +it in three different windows - the same slice, three very different images. +The figure is saved to radiology.png. +""" + +import matplotlib.pyplot as plt +import numpy as np +from pydicom.data import get_testdata_file + +from library import CT, MRI, Dicom, Radiation, Radiography, Ultrasound + + +def radiation_protection(): + """Distance and shielding: the two cheapest ways to cut dose.""" + radiation = Radiation() + print("Radiation protection") + print("-" * 48) + # A loop over distances shows the inverse square law in action. + for distance_m in (1, 2, 3, 4): + rate = radiation.inverse_square_law( + initial_intensity=100.0, initial_distance=1.0, new_distance=distance_m, + ) + print(f" {distance_m} m from source : {rate:6.2f} uGy/h") + + # Lead at ~100 keV: mu/rho ~5.55 cm^2/g (NIST XCOM), rho = 11.35 g/cm^3. + mu_lead = radiation.linear_attenuation_coefficient( + mass_attenuation_coefficient=5.55, density=11.35, + ) + hvl_mm = radiation.half_value_layer(attenuation_coefficient=mu_lead) * 10 + print(f" HVL of lead at 100 keV : {hvl_mm:.2f} mm\n") + + +def radiography_technique(): + """Adjust a technique chart for a new SID and a grid.""" + radiography = Radiography() + mas = radiography.mas(tube_current_ma=200, exposure_time_s=0.05) + new_mas = radiography.exposure_maintenance_mas(original_mas=mas, original_sid=100, + new_sid=180) + gridded = radiography.grid_mas(original_mas=new_mas, original_grid="none", new_grid="12:1") + print("Radiography technique") + print("-" * 48) + print(f" 10 mAs at 100 cm -> {new_mas:.1f} mAs at 180 cm -> {gridded:.0f} mAs with 12:1 grid\n") + + +def mri_contrast(): + """Why fat is bright on T1 and CSF is bright on T2.""" + mri = MRI() + # Each sequence is a dict of keyword arguments - unpacked with ** below. + sequences = { + "T1 (TR 500/TE 15)": {"repetition_time_ms": 500, "echo_time_ms": 15}, + "T2 (TR 4000/TE 100)": {"repetition_time_ms": 4000, "echo_time_ms": 100}, + } + print("MRI spin-echo signal (relative)") + print("-" * 48) + for name, timing in sequences.items(): + signals = {tissue: mri.spin_echo_signal(**timing, tissue=tissue) + for tissue in ("fat", "white_matter", "csf")} + brightest = max(signals, key=lambda tissue: signals[tissue]) + row = " ".join(f"{t}={s:.2f}" for t, s in signals.items()) + print(f" {name:20s} {row} -> brightest: {brightest}") + print() + + +def ultrasound_doppler(): + """Doppler velocity, and why the insonation angle matters.""" + ultrasound = Ultrasound() + print("Ultrasound Doppler (5 MHz, true flow 0.5 m/s, shift measured at 60 deg)") + print("-" * 48) + shift = ultrasound.doppler_shift(transmitted_frequency_mhz=5, velocity=0.5, angle_deg=60) + # If the operator mis-sets the angle cursor, the velocity estimate is wrong. + for assumed_angle in (50, 60, 70): + v = ultrasound.doppler_velocity(transmitted_frequency_mhz=5, doppler_shift_hz=shift, + angle_deg=assumed_angle) + print(f" angle cursor at {assumed_angle} deg -> {v:.2f} m/s") + print() + + +def ct_windowing(): + """Load a sample CT slice and display it in three windows.""" + dicom, ct = Dicom(), CT() + dataset = dicom.read(path=get_testdata_file("CT_small.dcm")) + hu = dicom.pixel_array_in_units(dataset=dataset) + + print("CT slice (pydicom sample - no patient data)") + print("-" * 48) + for key, value in dicom.summary(dataset=dataset).items(): + if value is not None: + print(f" {key:18s}: {value}") + print(f" HU range : {hu.min():.0f} to {hu.max():.0f}") + print(f" centre pixel : {hu[64, 64]:.0f} HU " + f"({ct.classify_hounsfield(hounsfield=hu[64, 64])})") + + presets = ["abdomen", "lung", "bone"] + figure, axes = plt.subplots(1, len(presets) + 1, figsize=(14, 4)) + + # Panel 1: narrow-beam attenuation curves for a few materials. + thickness_cm = np.linspace(0, 10, 200) + for label, mu in {"soft tissue (0.18/cm)": 0.18, "bone (0.48/cm)": 0.48}.items(): + axes[0].plot(thickness_cm, np.exp(-mu * thickness_cm) * 100, label=label) + axes[0].set(title="Attenuation at ~60 keV", xlabel="thickness [cm]", + ylabel="transmitted [%]") + axes[0].legend() + + # Panels 2-4: the same HU data through different windows. + for axis, preset in zip(axes[1:], presets): + width, level = CT.WINDOW_PRESETS[preset] + axis.imshow(ct.apply_window(hounsfield=hu, preset=preset), cmap="gray", + vmin=0, vmax=255) + axis.set_title(f"{preset} window (W{width} / L{level})") + axis.axis("off") + + figure.tight_layout() + figure.savefig("radiology.png", dpi=120) + print("\n Saved figure -> radiology.png") + + +if __name__ == "__main__": + radiation_protection() + radiography_technique() + mri_contrast() + ultrasound_doppler() + ct_windowing() diff --git a/tests/test_radiology.py b/tests/test_radiology.py new file mode 100644 index 0000000..a79048f --- /dev/null +++ b/tests/test_radiology.py @@ -0,0 +1,308 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2026 Evan Gress + +"""Tests for the radiology / medical-imaging classes. + +Each expected value is a textbook answer a radiology student could check by +hand, so these tests double as worked examples. +""" + +import math + +import numpy as np +import pytest +from pydicom.data import get_testdata_file + +from library import CT, MRI, Dicom, Radiation, Radiography, Ultrasound + +radiation = Radiation() +radiography = Radiography() +ct = CT() +mri = MRI() +ultrasound = Ultrasound() +dicom = Dicom() + + +# --- Radiation --------------------------------------------------------------- + +def test_inverse_square_doubling_distance_quarters_intensity(): + assert radiation.inverse_square_law( + initial_intensity=100.0, initial_distance=1.0, new_distance=2.0, + ) == pytest.approx(25.0) + + +def test_one_hvl_halves_the_beam(): + mu = 0.2 # 1/cm + hvl = radiation.half_value_layer(attenuation_coefficient=mu) + assert hvl == pytest.approx(math.log(2) / mu) + assert radiation.transmitted_intensity( + initial_intensity=1.0, attenuation_coefficient=mu, thickness=hvl, + ) == pytest.approx(0.5) + + +def test_tvl_is_about_3_32_hvls(): + tvl = radiation.tenth_value_layer(attenuation_coefficient=1.0) + hvl = radiation.half_value_layer(attenuation_coefficient=1.0) + assert tvl / hvl == pytest.approx(3.3219, rel=1e-4) + + +def test_shielding_thickness_inverts_attenuation(): + x = radiation.shielding_thickness(attenuation_coefficient=0.5, transmission=0.01) + assert radiation.transmitted_intensity( + initial_intensity=1.0, attenuation_coefficient=0.5, thickness=x, + ) == pytest.approx(0.01) + + +def test_duane_hunt_minimum_wavelength_at_100_kvp(): + assert radiation.photon_wavelength(energy_kev=100) == pytest.approx(0.0124, rel=1e-3) + + +def test_tc99m_after_two_half_lives(): + """Tc-99m (6.01 h): after 12.02 h a quarter of the activity remains.""" + assert radiation.remaining_activity( + initial_activity=800.0, half_life=6.01, elapsed_time=12.02, + ) == pytest.approx(200.0) + + +def test_effective_half_life_is_shorter_than_both(): + te = radiation.effective_half_life(physical_half_life=6.0, biological_half_life=24.0) + assert te == pytest.approx(4.8) + + +@pytest.mark.parametrize("radiation_type, expected_sv", [ + ("photon", 2.0), ("electron", 2.0), ("proton", 4.0), ("alpha", 40.0), +]) +def test_equivalent_dose(radiation_type, expected_sv): + assert radiation.equivalent_dose( + absorbed_dose=2.0, radiation_type=radiation_type, + ) == pytest.approx(expected_sv) + + +def test_unknown_radiation_type_raises(): + with pytest.raises(ValueError): + radiation.equivalent_dose(absorbed_dose=1.0, radiation_type="gamma-ish") + + +def test_neutron_weighting_peaks_near_20_at_1_mev(): + assert radiation.neutron_weighting_factor(energy_mev=1.0) == pytest.approx(20.7, abs=0.1) + + +def test_tissue_weighting_factors_sum_to_one(): + assert sum(Radiation.TISSUE_WEIGHTING_FACTORS.values()) == pytest.approx(1.0) + + +def test_uniform_whole_body_dose_equals_effective_dose(): + """If every organ gets 1 mSv, effective dose must also be 1 mSv.""" + uniform = {organ: 1.0 for organ in Radiation.TISSUE_WEIGHTING_FACTORS} + assert radiation.effective_dose(organ_equivalent_doses=uniform) == pytest.approx(1.0) + + +def test_partial_body_effective_dose(): + assert radiation.effective_dose( + organ_equivalent_doses={"thyroid": 10.0}, + ) == pytest.approx(0.4) + + +# --- Radiography ------------------------------------------------------------- + +def test_mas_reciprocity(): + assert radiography.mas(tube_current_ma=200, exposure_time_s=0.1) == pytest.approx( + radiography.mas(tube_current_ma=400, exposure_time_s=0.05) + ) + + +def test_fifteen_percent_rule_round_trip(): + kvp, mas = radiography.fifteen_percent_rule(kvp=70, mas=20, steps=1) + assert (kvp, mas) == (pytest.approx(80.5), pytest.approx(10.0)) + back = radiography.fifteen_percent_rule(kvp=kvp, mas=mas, steps=-1) + assert back == (pytest.approx(70), pytest.approx(20)) + + +def test_exposure_maintenance_100_to_180_cm(): + assert radiography.exposure_maintenance_mas( + original_mas=10, original_sid=100, new_sid=180, + ) == pytest.approx(32.4) + + +def test_grid_mas_adding_12_to_1_grid(): + assert radiography.grid_mas(original_mas=5, original_grid="none", new_grid="12:1") == 25 + + +def test_heat_units_high_frequency(): + assert radiography.heat_units( + kvp=80, tube_current_ma=200, exposure_time_s=0.1, generator="high_frequency", + ) == pytest.approx(2320.0) + + +def test_magnification_and_true_size(): + assert radiography.magnification_factor(sid=100, sod=80) == pytest.approx(1.25) + assert radiography.object_size(image_size=50, sid=100, sod=80) == pytest.approx(40) + + +def test_geometric_unsharpness(): + assert radiography.geometric_unsharpness( + focal_spot_size=1.0, sid=100, oid=20, + ) == pytest.approx(0.25) + + +def test_nyquist_limit_for_100_um_detector(): + assert radiography.nyquist_limit(pixel_pitch_mm=0.1) == pytest.approx(5.0) + + +def test_deviation_index_doubling_is_plus_three(): + assert radiography.exposure_index_deviation( + exposure_index=800, target_exposure_index=400, + ) == pytest.approx(3.01, abs=0.01) + + +# --- CT ---------------------------------------------------------------------- + +def test_hounsfield_reference_points(): + hu = ct.hounsfield_units(attenuation_coefficient=[0.0, 0.19, 0.38], water_attenuation=0.19) + np.testing.assert_allclose(hu, [-1000.0, 0.0, 1000.0]) + + +def test_rescale_to_hounsfield(): + np.testing.assert_allclose(ct.rescale_to_hounsfield(stored_values=[0, 1024]), [-1024, 0]) + + +def test_window_bounds_preset(): + assert ct.window_bounds(preset="brain") == (0.0, 80.0) + assert ct.window_bounds(preset="lung") == (-1350.0, 150.0) + + +def test_apply_window_maps_ends_to_black_and_white(): + grey = ct.apply_window(hounsfield=[-500, 0, 40, 80, 1000], preset="brain") + assert grey.dtype == np.uint8 + assert list(grey) == [0, 0, 128, 255, 255] + + +def test_window_requires_width_and_level(): + with pytest.raises(ValueError): + ct.window_bounds(window_width=400) + + +@pytest.mark.parametrize("hu, label", [ + (-1000, "air"), (-700, "lung / aerated tissue"), (-100, "fat"), + (0, "water / simple fluid"), (45, "soft tissue / blood"), + (1000, "bone / calcification / metal"), +]) +def test_classify_hounsfield(hu, label): + assert ct.classify_hounsfield(hounsfield=hu) == label + + +def test_head_ct_effective_dose(): + dlp = ct.dose_length_product(ctdi_vol_mgy=60, scan_length_cm=16) + assert dlp == pytest.approx(960) + assert ct.effective_dose_from_dlp(dlp_mgy_cm=dlp, region="head") == pytest.approx(2.016) + + +# --- MRI --------------------------------------------------------------------- + +def test_larmor_frequency_at_1_5_t(): + assert mri.larmor_frequency(field_strength_t=1.5) == pytest.approx(63.87, abs=0.01) + + +def test_relaxation_after_one_time_constant(): + assert mri.t1_recovery(time_ms=500, t1_ms=500) == pytest.approx(1 - math.exp(-1)) + assert mri.t2_decay(time_ms=80, t2_ms=80) == pytest.approx(math.exp(-1)) + + +def test_fat_brighter_than_csf_on_t1(): + tr, te = 500, 15 + fat = mri.spin_echo_signal(repetition_time_ms=tr, echo_time_ms=te, tissue="fat") + csf = mri.spin_echo_signal(repetition_time_ms=tr, echo_time_ms=te, tissue="csf") + assert fat > csf + + +def test_csf_brighter_than_white_matter_on_t2(): + tr, te = 4000, 100 + csf = mri.spin_echo_signal(repetition_time_ms=tr, echo_time_ms=te, tissue="csf") + wm = mri.spin_echo_signal(repetition_time_ms=tr, echo_time_ms=te, tissue="white_matter") + assert csf > wm + + +@pytest.mark.parametrize("tr, te, weighting", [ + (500, 15, "T1-weighted"), (4000, 100, "T2-weighted"), + (3000, 15, "PD-weighted"), (1500, 50, "mixed"), +]) +def test_sequence_weighting(tr, te, weighting): + assert mri.sequence_weighting(repetition_time_ms=tr, echo_time_ms=te) == weighting + + +# --- Ultrasound -------------------------------------------------------------- + +def test_wavelength_5_mhz_soft_tissue(): + assert ultrasound.wavelength(frequency_mhz=5) == pytest.approx(0.308) + + +def test_13_microseconds_per_cm_rule(): + assert ultrasound.echo_depth(round_trip_time_us=13) == pytest.approx(1.0, rel=0.01) + + +def test_tissue_air_reflects_almost_everything(): + assert ultrasound.reflection_between(medium_1="soft_tissue", medium_2="air") > 0.99 + assert ultrasound.reflection_between(medium_1="soft_tissue", medium_2="fat") < 0.01 + + +def test_attenuation_round_trip(): + # 3 MHz, 10 cm deep -> 20 cm path -> 30 dB + assert ultrasound.attenuation_db(frequency_mhz=3, path_length_cm=20) == pytest.approx(30) + + +def test_maximum_depth_for_prf(): + assert ultrasound.maximum_depth(pulse_repetition_frequency_hz=5000) == pytest.approx(15.4) + + +def test_doppler_shift_and_inverse_agree(): + shift = ultrasound.doppler_shift(transmitted_frequency_mhz=5, velocity=0.5, angle_deg=60) + assert shift == pytest.approx(1623.4, rel=1e-3) + velocity = ultrasound.doppler_velocity( + transmitted_frequency_mhz=5, doppler_shift_hz=shift, angle_deg=60, + ) + assert velocity == pytest.approx(0.5) + + +def test_doppler_at_90_degrees_raises(): + with pytest.raises(ValueError): + ultrasound.doppler_velocity(transmitted_frequency_mhz=5, doppler_shift_hz=100, + angle_deg=90) + + +# --- DICOM (uses pydicom's bundled sample CT - no patient data) -------------- + +@pytest.fixture +def ct_dataset(): + return dicom.read(path=get_testdata_file("CT_small.dcm")) + + +def test_dicom_summary_has_modality(ct_dataset): + summary = dicom.summary(dataset=ct_dataset) + assert summary["Modality"] == "CT" + assert "PatientName" not in summary + + +def test_dicom_pixels_are_converted_to_hu(ct_dataset): + raw = ct_dataset.pixel_array.astype(float) + hu = dicom.pixel_array_in_units(dataset=ct_dataset) + np.testing.assert_allclose(hu, raw * float(ct_dataset.RescaleSlope) + + float(ct_dataset.RescaleIntercept)) + + +def test_dicom_anonymize_copies_and_strips_phi(ct_dataset, tmp_path): + original_name = str(ct_dataset.PatientName) + clean = dicom.anonymize(dataset=ct_dataset) + assert str(clean.PatientName) == "ANONYMIZED" + assert str(ct_dataset.PatientName) == original_name # original untouched + + out = tmp_path / "anon.dcm" + dicom.save(dataset=clean, path=out) + assert str(dicom.read(path=out).PatientName) == "ANONYMIZED" + + +def test_windowed_slice_exports_to_png(ct_dataset, tmp_path): + hu = dicom.pixel_array_in_units(dataset=ct_dataset) + image = ct.apply_window(hounsfield=hu, preset="abdomen") + out = tmp_path / "slice.png" + dicom.to_png(image=image, path=out) + assert out.stat().st_size > 0 diff --git a/uv.lock b/uv.lock index 502023a..277d8cd 100644 --- a/uv.lock +++ b/uv.lock @@ -61,6 +61,7 @@ dependencies = [ { name = "matplotlib" }, { name = "numpy" }, { name = "pint" }, + { name = "pydicom" }, { name = "scipy" }, { name = "scipy-stubs" }, ] @@ -75,6 +76,7 @@ requires-dist = [ { name = "matplotlib", specifier = ">=3.10.9" }, { name = "numpy", specifier = ">=2.4.6" }, { name = "pint", specifier = ">=0.25.3" }, + { name = "pydicom", specifier = ">=3.0.2" }, { name = "scipy", specifier = ">=1.17.1" }, { name = "scipy-stubs", specifier = "~=1.17.1" }, ] @@ -342,6 +344,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/54/20/4d324d65cc6d9205fabedc306948156824eb9f0ee1633355a8f7ec5c66bf/pluggy-1.6.0-py3-none-any.whl", hash = "sha256:e920276dd6813095e9377c0bc5566d94c932c33b27a3e3945d8389c374dd4746", size = 20538, upload-time = "2025-05-15T12:30:06.134Z" }, ] +[[package]] +name = "pydicom" +version = "3.0.2" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/7a/de/52aaf905f1f0ae7aba85996e2592ea2c1fe49157f3cfbcd1871965bdb51d/pydicom-3.0.2.tar.gz", hash = "sha256:5942bfc2d72c6fa4b3b5b62c527f54b7f2355f21d6f5d296df6bb30188df6a4f", size = 2886792, upload-time = "2026-03-19T21:46:20.935Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/46/e0/60466c6d712dad2cf807df315e39863e91609ffd1064ecb835994460bbda/pydicom-3.0.2-py3-none-any.whl", hash = "sha256:abf971a5440f84dbaf42c4b6758e30e62480902584f8b270b9a5d146e278a07b", size = 2376822, upload-time = "2026-03-19T21:46:19.042Z" }, +] + [[package]] name = "pygments" version = "2.20.0"