openbnct_core/
external_dose.rs

1// SPDX-License-Identifier: Apache-2.0
2
3//! External photon/hadron dose import (`openbnct.external-dose/0.1.0`) and
4//! dose-field resampling.
5//!
6//! This is the single-scalar counterpart of the component-dose interchange:
7//! an external pipeline (treatment-planning export, MCNP/PHITS/Geant4 dose
8//! tallies, a research code) contributes one absolute absorbed-dose field on
9//! a declared grid plus the fractionation the dose was delivered in. The
10//! importer validates it into an [`ExternalDoseBundle`] whose provenance
11//! binds the document hash (`external-dose:<system>:sha256:<hash>`).
12//!
13//! Fractionation is declared, never guessed: `uniform` splits the total into
14//! `count` equal per-fraction doses; `explicit` carries the full per-fraction
15//! dose arrays. The biological layer (`openbnct-bio`) turns either into a
16//! BED/EQD2 field; combining that field with a BNCT biological bundle is a
17//! separate explicit step with its own compatibility gates.
18
19use serde::{Deserialize, Serialize};
20use thiserror::Error;
21
22use crate::{ExternalProducer, GridGeometry, ValidationError, grid_geometry_equivalent};
23
24/// Current schema token for external dose documents.
25pub const EXTERNAL_DOSE_SCHEMA: &str = "openbnct.external-dose/0.1.0";
26
27/// What the imported dose field physically is. `physical` is absorbed dose;
28/// `rbe_weighted` declares the producer already applied an RBE model — the
29/// basis is carried through BED conversion so a weighted field can never be
30/// silently read as physical dose downstream.
31#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
32#[serde(rename_all = "snake_case")]
33pub enum ExternalDoseQuantity {
34    Physical,
35    RbeWeighted,
36}
37
38/// How the total dose divides into fractions.
39#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
40#[serde(tag = "kind", rename_all = "snake_case", deny_unknown_fields)]
41pub enum ExternalFractionation {
42    /// `count` identical fractions; per-fraction dose is `values / count`.
43    Uniform { count: u32 },
44    /// Explicit per-fraction dose arrays, each voxel-aligned to `geometry`.
45    Explicit { doses: Vec<Vec<f64>> },
46}
47
48/// A `openbnct.external-dose/0.1.0` document as emitted by the producer.
49#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
50#[serde(deny_unknown_fields)]
51pub struct ExternalDoseDocument {
52    #[serde(deserialize_with = "crate::deserialize_contract_id")]
53    pub schema_version: String,
54    pub case_id: String,
55    /// DICOM frame-of-reference UID tying the grid to a patient space.
56    #[serde(default, skip_serializing_if = "Option::is_none")]
57    pub frame_of_reference_uid: Option<String>,
58    pub geometry: GridGeometry,
59    pub producer: ExternalProducer,
60    /// What the dose physically is; `rbe_weighted` producers must say so.
61    pub quantity: ExternalDoseQuantity,
62    /// Absolute absorbed dose in gray over the whole course.
63    pub values: Vec<f64>,
64    /// Optional per-voxel 1-sigma absolute uncertainty, same unit.
65    #[serde(default, skip_serializing_if = "Option::is_none")]
66    pub absolute_standard_uncertainty: Option<Vec<f64>>,
67    /// The fractionation the course was delivered in.
68    pub fractionation: ExternalFractionation,
69}
70
71/// The validated, provenance-bound form produced by [`import_external_dose`].
72#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
73#[serde(deny_unknown_fields)]
74pub struct ExternalDoseBundle {
75    #[serde(deserialize_with = "crate::deserialize_contract_id")]
76    pub schema_version: String,
77    pub case_id: String,
78    #[serde(default, skip_serializing_if = "Option::is_none")]
79    pub frame_of_reference_uid: Option<String>,
80    pub geometry: GridGeometry,
81    pub producer: ExternalProducer,
82    pub quantity: ExternalDoseQuantity,
83    pub values: Vec<f64>,
84    #[serde(default, skip_serializing_if = "Option::is_none")]
85    pub absolute_standard_uncertainty: Option<Vec<f64>>,
86    pub fractionation: ExternalFractionation,
87    /// `external-dose:<system>:sha256:<document hash>`.
88    pub provenance_id: String,
89}
90
91impl ExternalDoseBundle {
92    /// Number of fractions declared by the producer.
93    #[must_use]
94    pub fn fraction_count(&self) -> usize {
95        match &self.fractionation {
96            ExternalFractionation::Uniform { count } => *count as usize,
97            ExternalFractionation::Explicit { doses } => doses.len(),
98        }
99    }
100
101    /// Per-fraction dose at `voxel`, expanding a uniform split.
102    pub fn fraction_dose(&self, fraction: usize, voxel: usize) -> f64 {
103        match &self.fractionation {
104            ExternalFractionation::Uniform { count } => self.values[voxel] / f64::from(*count),
105            ExternalFractionation::Explicit { doses } => doses[fraction][voxel],
106        }
107    }
108}
109
110#[derive(Debug, Error)]
111pub enum ExternalDoseError {
112    #[error("unsupported schema_version {0:?}")]
113    UnsupportedSchema(String),
114    #[error("field {0} must be nonempty")]
115    EmptyField(&'static str),
116    #[error("invalid geometry: {0}")]
117    InvalidGeometry(#[from] ValidationError),
118    #[error("dose values are not finite non-negative voxel values")]
119    InvalidDoseValues,
120    #[error("uncertainty values are not finite non-negative voxel values")]
121    InvalidUncertainty,
122    #[error("fractionation declares {0} doses; each must cover the grid")]
123    InvalidFractionDoses(usize),
124    #[error("fraction count must be at least one")]
125    ZeroFractions,
126}
127
128/// Validate an external-dose document into a provenance-bound bundle.
129///
130/// `document_sha256` is the SHA-256 of the document's bytes — the caller
131/// hashes the exact file so provenance binds content, not a name.
132pub fn import_external_dose(
133    document: &ExternalDoseDocument,
134    document_sha256: &str,
135) -> Result<ExternalDoseBundle, ExternalDoseError> {
136    if !crate::schema_matches(&document.schema_version, EXTERNAL_DOSE_SCHEMA) {
137        return Err(ExternalDoseError::UnsupportedSchema(
138            document.schema_version.clone(),
139        ));
140    }
141    if document.case_id.trim().is_empty() {
142        return Err(ExternalDoseError::EmptyField("case_id"));
143    }
144    if document.producer.system.trim().is_empty() {
145        return Err(ExternalDoseError::EmptyField("producer.system"));
146    }
147    if document.producer.version.trim().is_empty() {
148        return Err(ExternalDoseError::EmptyField("producer.version"));
149    }
150    if document.producer.normalization.trim().is_empty() {
151        return Err(ExternalDoseError::EmptyField("producer.normalization"));
152    }
153    if let Some(uid) = &document.frame_of_reference_uid
154        && uid.trim().is_empty()
155    {
156        return Err(ExternalDoseError::EmptyField("frame_of_reference_uid"));
157    }
158    let voxel_count = document.geometry.voxel_count()?;
159    if document.values.len() != voxel_count
160        || document.values.iter().any(|v| !v.is_finite() || *v < 0.0)
161    {
162        return Err(ExternalDoseError::InvalidDoseValues);
163    }
164    if let Some(sigma) = &document.absolute_standard_uncertainty
165        && (sigma.len() != voxel_count || sigma.iter().any(|v| !v.is_finite() || *v < 0.0))
166    {
167        return Err(ExternalDoseError::InvalidUncertainty);
168    }
169    match &document.fractionation {
170        ExternalFractionation::Uniform { count } if *count == 0 => {
171            return Err(ExternalDoseError::ZeroFractions);
172        }
173        ExternalFractionation::Explicit { doses } => {
174            if doses.is_empty() {
175                return Err(ExternalDoseError::ZeroFractions);
176            }
177            for dose in doses {
178                if dose.len() != voxel_count || dose.iter().any(|v| !v.is_finite() || *v < 0.0) {
179                    return Err(ExternalDoseError::InvalidFractionDoses(doses.len()));
180                }
181            }
182            // Explicit fractions must sum to the declared total — otherwise
183            // the two dose statements silently disagree.
184            for voxel in 0..voxel_count {
185                let sum: f64 = doses.iter().map(|d| d[voxel]).sum();
186                let total = document.values[voxel];
187                if (sum - total).abs() > 1e-6 * total.abs().max(1e-12) {
188                    return Err(ExternalDoseError::InvalidFractionDoses(doses.len()));
189                }
190            }
191        }
192        ExternalFractionation::Uniform { .. } => {}
193    }
194
195    Ok(ExternalDoseBundle {
196        schema_version: document.schema_version.clone(),
197        case_id: document.case_id.clone(),
198        frame_of_reference_uid: document.frame_of_reference_uid.clone(),
199        geometry: document.geometry.clone(),
200        producer: document.producer.clone(),
201        quantity: document.quantity,
202        values: document.values.clone(),
203        absolute_standard_uncertainty: document.absolute_standard_uncertainty.clone(),
204        fractionation: document.fractionation.clone(),
205        provenance_id: format!(
206            "external-dose:{}:sha256:{document_sha256}",
207            document.producer.system
208        ),
209    })
210}
211
212/// Dose-field resampling method. Trilinear interpolation at voxel centers is
213/// the only supported scheme; the choice is recorded in whatever record the
214/// resampled field lands in.
215#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
216#[serde(rename_all = "snake_case")]
217pub enum ResampleMethod {
218    Trilinear,
219}
220
221#[derive(Debug, Error, PartialEq)]
222pub enum ResampleError {
223    #[error("resampling requires axis-aligned (identity-direction) grids")]
224    NonAxisAligned,
225    #[error("source field has {0} values; grid holds {1} voxels")]
226    LengthMismatch(usize, usize),
227    #[error("target grid is not covered by the source field on axis {axis}")]
228    UncoveredTarget { axis: usize },
229    #[error("invalid geometry: {0}")]
230    InvalidGeometry(#[from] ValidationError),
231}
232
233/// Resample a voxel field onto a target grid by trilinear interpolation at
234/// voxel centers.
235///
236/// Strict co-registration: both grids must be axis-aligned and every target
237/// voxel center must lie inside the source field's outer bounds (source
238/// centers at the edge extrapolate flatly to the half-spacing face). A
239/// target point outside the source extent is a hard error — never a silent
240/// zero.
241pub fn resample_trilinear(
242    values: &[f64],
243    source: &GridGeometry,
244    target: &GridGeometry,
245) -> Result<Vec<f64>, ResampleError> {
246    let identity = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0];
247    if source.direction != identity || target.direction != identity {
248        return Err(ResampleError::NonAxisAligned);
249    }
250    let src_nx = source.shape[0] as usize;
251    let src_ny = source.shape[1] as usize;
252    let src_nz = source.shape[2] as usize;
253    let src_count = src_nx * src_ny * src_nz;
254    if values.len() != src_count {
255        return Err(ResampleError::LengthMismatch(values.len(), src_count));
256    }
257    let target_count = target.voxel_count()?;
258    // Reuse the shared equivalent-grid check: an identical grid resamples to
259    // the same values with no interpolation cost.
260    if grid_geometry_equivalent(source, target) {
261        return Ok(values.to_vec());
262    }
263
264    let mut out = Vec::with_capacity(target_count);
265    let [tnx, tny, tnz] = target.shape.map(|d| d as usize);
266    for k in 0..tnz {
267        for j in 0..tny {
268            for i in 0..tnx {
269                let mut value = 0.0;
270                let mut weights = [(0_usize, 0.0_f64); 8];
271                for axis in 0..3 {
272                    let p = match axis {
273                        0 => target.origin_mm[0] + i as f64 * target.spacing_mm[0],
274                        1 => target.origin_mm[1] + j as f64 * target.spacing_mm[1],
275                        _ => target.origin_mm[2] + k as f64 * target.spacing_mm[2],
276                    };
277                    // Fractional source index of the target center; a target
278                    // center may sit at most half a spacing outside the outer
279                    // centers (flat extrapolation to the face), never beyond.
280                    let f = (p - source.origin_mm[axis]) / source.spacing_mm[axis];
281                    let n = source.shape[axis] as f64;
282                    if !(-0.5..=n - 0.5).contains(&f) {
283                        return Err(ResampleError::UncoveredTarget { axis });
284                    }
285                    let f = f.clamp(0.0, n - 1.0);
286                    let lo = f.floor() as usize;
287                    let hi = (lo + 1).min(source.shape[axis] as usize - 1);
288                    weights[axis * 2] = (lo, 1.0 - (f - lo as f64));
289                    weights[axis * 2 + 1] = (hi, f - lo as f64);
290                }
291                for (bi, wi) in [weights[0], weights[1]] {
292                    for (bj, wj) in [weights[2], weights[3]] {
293                        for (bk, wk) in [weights[4], weights[5]] {
294                            value += values[bi + src_nx * bj + src_nx * src_ny * bk] * wi * wj * wk;
295                        }
296                    }
297                }
298                out.push(value);
299            }
300        }
301    }
302    Ok(out)
303}
304
305#[cfg(test)]
306mod tests {
307    use super::*;
308
309    fn grid() -> GridGeometry {
310        GridGeometry {
311            shape: [2, 1, 1],
312            spacing_mm: [5.0, 5.0, 5.0],
313            origin_mm: [0.0, 0.0, 0.0],
314            direction: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
315        }
316    }
317
318    fn producer() -> ExternalProducer {
319        ExternalProducer {
320            system: "course-sim".into(),
321            version: "1".into(),
322            normalization: "absolute gray".into(),
323        }
324    }
325
326    fn document() -> ExternalDoseDocument {
327        ExternalDoseDocument {
328            schema_version: EXTERNAL_DOSE_SCHEMA.into(),
329            case_id: "case".into(),
330            frame_of_reference_uid: None,
331            geometry: grid(),
332            producer: producer(),
333            quantity: ExternalDoseQuantity::Physical,
334            values: vec![60.0, 60.0],
335            absolute_standard_uncertainty: Some(vec![0.6, 0.6]),
336            fractionation: ExternalFractionation::Uniform { count: 30 },
337        }
338    }
339
340    #[test]
341    fn valid_document_binds_provenance() {
342        let bundle = import_external_dose(&document(), "ab12").unwrap();
343        assert_eq!(bundle.provenance_id, "external-dose:course-sim:sha256:ab12");
344        assert_eq!(bundle.fraction_count(), 30);
345        assert!((bundle.fraction_dose(0, 0) - 2.0).abs() < 1e-12);
346    }
347
348    #[test]
349    fn rejects_wrong_schema_and_bad_values() {
350        let mut bad = document();
351        bad.schema_version = "openbnct.external-dose/0.0.0".into();
352        assert!(import_external_dose(&bad, "x").is_err());
353        let mut bad = document();
354        bad.values.push(1.0); // 3 values on a 2-voxel grid
355        assert!(import_external_dose(&bad, "x").is_err());
356        let mut bad = document();
357        bad.values[0] = f64::NAN;
358        assert!(import_external_dose(&bad, "x").is_err());
359        let mut bad = document();
360        bad.fractionation = ExternalFractionation::Uniform { count: 0 };
361        assert!(import_external_dose(&bad, "x").is_err());
362    }
363
364    #[test]
365    fn explicit_fractions_must_sum_to_total() {
366        let mut doc = document();
367        doc.fractionation = ExternalFractionation::Explicit {
368            doses: vec![vec![30.0, 30.0], vec![30.0, 30.0]],
369        };
370        assert!(import_external_dose(&doc, "x").is_ok());
371        doc.fractionation = ExternalFractionation::Explicit {
372            doses: vec![vec![30.0, 30.0], vec![29.0, 30.0]],
373        };
374        assert!(import_external_dose(&doc, "x").is_err());
375    }
376
377    #[test]
378    fn resample_identity_is_exact() {
379        let values = vec![1.0, 2.0];
380        assert_eq!(
381            resample_trilinear(&values, &grid(), &grid()).unwrap(),
382            values
383        );
384    }
385
386    #[test]
387    fn resample_trilinear_interpolates_centers() {
388        // Source: 2 voxels at x = 0, 5 mm with values 0 and 10.
389        // Target: 3 voxels at x = 0, 2.5, 5 mm → 0, 5, 10.
390        let target = GridGeometry {
391            shape: [3, 1, 1],
392            spacing_mm: [2.5, 5.0, 5.0],
393            ..grid()
394        };
395        let out = resample_trilinear(&[0.0, 10.0], &grid(), &target).unwrap();
396        assert!((out[0] - 0.0).abs() < 1e-12);
397        assert!((out[1] - 5.0).abs() < 1e-12);
398        assert!((out[2] - 10.0).abs() < 1e-12);
399    }
400
401    #[test]
402    fn resample_rejects_uncovered_target() {
403        // Target centered beyond the source extent (x = 10 mm past the edge).
404        let target = GridGeometry {
405            shape: [1, 1, 1],
406            origin_mm: [20.0, 0.0, 0.0],
407            ..grid()
408        };
409        assert_eq!(
410            resample_trilinear(&[1.0, 2.0], &grid(), &target),
411            Err(ResampleError::UncoveredTarget { axis: 0 })
412        );
413    }
414
415    #[test]
416    fn resample_rejects_rotated_grids() {
417        let mut rotated = grid();
418        rotated.direction = [0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0];
419        assert_eq!(
420            resample_trilinear(&[1.0, 2.0], &rotated, &grid()),
421            Err(ResampleError::NonAxisAligned)
422        );
423    }
424}
openbnct_core/external_dose.rs

openbnct_core/
external_dose.rs
