feat: Add radiology and medical-imaging classes

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>`.
+
+| 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 <package>`.
 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

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

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 = []

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_",

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