openbnct_core/
lib.rs

1// SPDX-License-Identifier: Apache-2.0
2
3#![forbid(unsafe_code)]
4
5mod exposure;
6mod external_dose;
7mod interchange;
8mod registration;
9mod stats;
10mod systematic;
11
12use std::collections::BTreeSet;
13
14use serde::{Deserialize, Serialize};
15use thiserror::Error;
16
17pub use exposure::{
18    BoundFileReference, EXPOSURE_PLAN_SCHEMA, Exposure, ExposureCovariance, ExposurePlan,
19    ExposurePlanError, WeightBasis, accumulate_exposures,
20};
21pub use external_dose::{
22    EXTERNAL_DOSE_SCHEMA, ExternalDoseBundle, ExternalDoseDocument, ExternalDoseError,
23    ExternalDoseQuantity, ExternalFractionation, ResampleError, ResampleMethod,
24    import_external_dose, resample_trilinear,
25};
26pub use interchange::{
27    COMPONENT_DOSE_INTERCHANGE_SCHEMA, ComponentDoseInterchange, ExternalProducer, ExternalTotal,
28    InterchangeError, grid_geometry_equivalent, import_component_dose,
29};
30pub use registration::{
31    LandmarkPair, REGISTRATION_SCHEMA, Registration, RegistrationError, RegistrationMethod,
32    RigidTransform, declared_registration, fit_landmark_transform, landmark_registration,
33};
34pub use stats::{
35    dose_covering_percent, equivalent_uniform_dose, masked_values, mean, volume_at_least,
36};
37pub use systematic::{
38    RegionUncertainty, SYSTEMATIC_UNCERTAINTY_QUALIFICATION, SYSTEMATIC_UNCERTAINTY_SCHEMA,
39    SourceSummary, SystematicError, SystematicUncertaintyReport, UncertaintySource,
40    boron_field_sigma, combine_total_sigma, combine_voxel_sigma, positioning_sigma,
41    region_uncertainty, relative_component_sigma, summarize_source,
42};
43
44/// A regular patient-coordinate voxel grid.
45///
46/// `direction` is row-major and maps voxel axes into the patient coordinate
47/// frame. Geometry importers must preserve the original DICOM frame of
48/// reference outside this numerical representation.
49#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
50#[serde(deny_unknown_fields)]
51pub struct GridGeometry {
52    pub shape: [u32; 3],
53    pub spacing_mm: [f64; 3],
54    pub origin_mm: [f64; 3],
55    pub direction: [f64; 9],
56}
57
58impl GridGeometry {
59    pub fn voxel_count(&self) -> Result<usize, ValidationError> {
60        if self.shape.contains(&0) {
61            return Err(ValidationError::EmptyGrid);
62        }
63        if self
64            .spacing_mm
65            .iter()
66            .any(|value| !value.is_finite() || *value <= 0.0)
67        {
68            return Err(ValidationError::InvalidSpacing);
69        }
70        if self.origin_mm.iter().any(|value| !value.is_finite())
71            || self.direction.iter().any(|value| !value.is_finite())
72        {
73            return Err(ValidationError::NonFiniteGeometry);
74        }
75        let axes = [
76            [self.direction[0], self.direction[3], self.direction[6]],
77            [self.direction[1], self.direction[4], self.direction[7]],
78            [self.direction[2], self.direction[5], self.direction[8]],
79        ];
80        let dot = |left: [f64; 3], right: [f64; 3]| {
81            left[0].mul_add(right[0], left[1].mul_add(right[1], left[2] * right[2]))
82        };
83        let determinant = self.direction[0]
84            * (self.direction[4] * self.direction[8] - self.direction[5] * self.direction[7])
85            - self.direction[1]
86                * (self.direction[3] * self.direction[8] - self.direction[5] * self.direction[6])
87            + self.direction[2]
88                * (self.direction[3] * self.direction[7] - self.direction[4] * self.direction[6]);
89        const TOLERANCE: f64 = 1.0e-6;
90        if axes
91            .iter()
92            .any(|axis| (dot(*axis, *axis) - 1.0).abs() > TOLERANCE)
93            || dot(axes[0], axes[1]).abs() > TOLERANCE
94            || dot(axes[0], axes[2]).abs() > TOLERANCE
95            || dot(axes[1], axes[2]).abs() > TOLERANCE
96            || (determinant - 1.0).abs() > TOLERANCE
97        {
98            return Err(ValidationError::InvalidDirection);
99        }
100
101        self.shape.iter().try_fold(1_usize, |count, extent| {
102            count
103                .checked_mul(*extent as usize)
104                .ok_or(ValidationError::GridTooLarge)
105        })
106    }
107
108    /// Return the patient-coordinate center of a voxel in millimetres.
109    pub fn voxel_center_lps_mm(&self, voxel: [u32; 3]) -> Result<[f64; 3], ValidationError> {
110        self.voxel_count()?;
111        if voxel
112            .into_iter()
113            .zip(self.shape)
114            .any(|(index, extent)| index >= extent)
115        {
116            return Err(ValidationError::VoxelOutOfBounds {
117                voxel,
118                shape: self.shape,
119            });
120        }
121        let local = [
122            f64::from(voxel[0]) * self.spacing_mm[0],
123            f64::from(voxel[1]) * self.spacing_mm[1],
124            f64::from(voxel[2]) * self.spacing_mm[2],
125        ];
126        Ok([
127            self.origin_mm[0]
128                + self.direction[0].mul_add(
129                    local[0],
130                    self.direction[1].mul_add(local[1], self.direction[2] * local[2]),
131                ),
132            self.origin_mm[1]
133                + self.direction[3].mul_add(
134                    local[0],
135                    self.direction[4].mul_add(local[1], self.direction[5] * local[2]),
136                ),
137            self.origin_mm[2]
138                + self.direction[6].mul_add(
139                    local[0],
140                    self.direction[7].mul_add(local[1], self.direction[8] * local[2]),
141                ),
142        ])
143    }
144
145    /// World-axis-aligned bounding box `(minimum, maximum)` of the grid's
146    /// voxel extents in LPS millimetres. `origin_mm` is voxel index
147    /// `[0,0,0]`'s center, so each face sits half a spacing beyond the
148    /// extreme centers along the (possibly rotated) voxel axes.
149    pub fn bounding_box_lps_mm(&self) -> Result<([f64; 3], [f64; 3]), ValidationError> {
150        self.voxel_count()?;
151        let mut minimum = self.origin_mm;
152        let mut maximum = self.origin_mm;
153        // Axis `axis`'s half-extent in voxel-index units is `shape/2` about
154        // the [0,0,0] center; walk the eight extreme-index corners.
155        for corner in 0..8 {
156            let local = [
157                (if corner & 1 == 0 {
158                    -0.5
159                } else {
160                    self.shape[0] as f64 - 0.5
161                }) * self.spacing_mm[0],
162                (if corner & 2 == 0 {
163                    -0.5
164                } else {
165                    self.shape[1] as f64 - 0.5
166                }) * self.spacing_mm[1],
167                (if corner & 4 == 0 {
168                    -0.5
169                } else {
170                    self.shape[2] as f64 - 0.5
171                }) * self.spacing_mm[2],
172            ];
173            for axis in 0..3 {
174                let value = self.origin_mm[axis]
175                    + self.direction[axis * 3].mul_add(
176                        local[0],
177                        self.direction[axis * 3 + 1]
178                            .mul_add(local[1], self.direction[axis * 3 + 2] * local[2]),
179                    );
180                minimum[axis] = minimum[axis].min(value);
181                maximum[axis] = maximum[axis].max(value);
182            }
183        }
184        Ok((minimum, maximum))
185    }
186}
187
188/// The four physical dose groups retained before biological weighting.
189#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Serialize, Deserialize)]
190#[serde(rename_all = "snake_case")]
191pub enum DoseComponent {
192    Boron,
193    Nitrogen,
194    Hydrogen,
195    Photon,
196}
197
198impl DoseComponent {
199    pub const REQUIRED: [Self; 4] = [Self::Boron, Self::Nitrogen, Self::Hydrogen, Self::Photon];
200}
201
202#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
203#[serde(rename_all = "snake_case")]
204pub enum DoseUnit {
205    Gray,
206    GrayPerSourceParticle,
207}
208
209/// Immutable identity of a scientific input artifact.
210#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
211#[serde(deny_unknown_fields)]
212pub struct ContentReference {
213    pub id: String,
214    pub sha256: String,
215}
216
217impl ContentReference {
218    pub fn validate(&self) -> Result<(), ContentReferenceError> {
219        if self.id.trim().is_empty() {
220            return Err(ContentReferenceError::EmptyId);
221        }
222        if !is_canonical_sha256(&self.sha256) {
223            return Err(ContentReferenceError::InvalidSha256);
224        }
225        Ok(())
226    }
227}
228
229/// Content reference used specifically for the component-definition profile.
230pub type ComponentProfileReference = ContentReference;
231
232#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
233#[serde(deny_unknown_fields)]
234pub struct DoseVolume {
235    pub component: DoseComponent,
236    pub unit: DoseUnit,
237    pub values: Vec<f64>,
238    /// One-sigma absolute standard uncertainty in the same unit as `values`.
239    pub absolute_standard_uncertainty: Option<Vec<f64>>,
240}
241
242impl DoseVolume {
243    /// Derive relative uncertainty for one voxel.
244    ///
245    /// Relative uncertainty is deliberately absent when no absolute
246    /// uncertainty exists or when the mean is zero.
247    pub fn relative_standard_uncertainty(
248        &self,
249        voxel_index: usize,
250    ) -> Result<Option<f64>, ValidationError> {
251        let mean = self
252            .values
253            .get(voxel_index)
254            .ok_or(ValidationError::DoseIndexOutOfBounds {
255                index: voxel_index,
256                length: self.values.len(),
257            })?;
258        let Some(uncertainty) = &self.absolute_standard_uncertainty else {
259            return Ok(None);
260        };
261        let absolute = uncertainty
262            .get(voxel_index)
263            .ok_or(ValidationError::UncertaintyLength {
264                component: self.component,
265                expected: self.values.len(),
266                actual: uncertainty.len(),
267            })?;
268        if *mean == 0.0 {
269            Ok(None)
270        } else {
271            Ok(Some(*absolute / *mean))
272        }
273    }
274}
275
276#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
277#[serde(rename_all = "snake_case")]
278pub enum TotalUncertaintyMethod {
279    /// Uncertainty comes from a dedicated physical-total estimator.
280    DedicatedEstimator,
281    /// Uncertainty was calculated from batch-level component covariance.
282    BatchCovariance,
283    /// No defensible total uncertainty is available.
284    Unavailable,
285}
286
287#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
288#[serde(deny_unknown_fields)]
289pub struct PhysicalTotalDoseVolume {
290    pub unit: DoseUnit,
291    pub values: Vec<f64>,
292    pub absolute_standard_uncertainty: Option<Vec<f64>>,
293    pub uncertainty_method: TotalUncertaintyMethod,
294}
295
296/// Named voxel mask in grid order (`i + nx*j + nx*ny*k`).
297#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
298#[serde(deny_unknown_fields)]
299pub struct RegionMask {
300    pub name: String,
301    pub voxels: Vec<bool>,
302}
303
304impl RegionMask {
305    /// Number of voxels included in the mask.
306    #[must_use]
307    pub fn included_voxel_count(&self) -> usize {
308        self.voxels.iter().filter(|voxel| **voxel).count()
309    }
310
311    /// `self` minus `other`, under `name`.
312    ///
313    /// Both masks must describe the same voxel count — masks do not carry
314    /// their grid, so the caller binds them to a shared case geometry. The
315    /// result must select at least one voxel.
316    pub fn subtract(
317        &self,
318        other: &RegionMask,
319        name: impl Into<String>,
320    ) -> Result<RegionMask, ValidationError> {
321        self.combine(other, name, |a, b| a && !b)
322    }
323
324    /// Union of `self` and `other`, under `name`.
325    pub fn union(
326        &self,
327        other: &RegionMask,
328        name: impl Into<String>,
329    ) -> Result<RegionMask, ValidationError> {
330        self.combine(other, name, |a, b| a || b)
331    }
332
333    /// Intersection of `self` and `other`, under `name`.
334    pub fn intersection(
335        &self,
336        other: &RegionMask,
337        name: impl Into<String>,
338    ) -> Result<RegionMask, ValidationError> {
339        self.combine(other, name, |a, b| a && b)
340    }
341
342    fn combine(
343        &self,
344        other: &RegionMask,
345        name: impl Into<String>,
346        op: impl Fn(bool, bool) -> bool,
347    ) -> Result<RegionMask, ValidationError> {
348        if self.voxels.len() != other.voxels.len() {
349            return Err(ValidationError::MaskVoxelCountMismatch {
350                name: self.name.clone(),
351                other: other.name.clone(),
352                expected: self.voxels.len(),
353                actual: other.voxels.len(),
354            });
355        }
356        let voxels = self
357            .voxels
358            .iter()
359            .copied()
360            .zip(other.voxels.iter().copied())
361            .map(|(a, b)| op(a, b))
362            .collect::<Vec<_>>();
363        let mask = RegionMask {
364            name: name.into(),
365            voxels,
366        };
367        if mask.included_voxel_count() == 0 {
368            return Err(ValidationError::EmptyMask(mask.name));
369        }
370        Ok(mask)
371    }
372
373    /// Mean and maximum of a per-voxel quantity inside this mask.
374    ///
375    /// `values` must be laid out on the same voxel count as the mask.
376    pub fn summarize(&self, values: &[f64]) -> Result<MaskDoseSummary, ValidationError> {
377        if values.len() != self.voxels.len() {
378            return Err(ValidationError::MaskValuesLength {
379                mask: self.name.clone(),
380                mask_voxels: self.voxels.len(),
381                values: values.len(),
382            });
383        }
384        let included = self
385            .voxels
386            .iter()
387            .copied()
388            .zip(values.iter().copied())
389            .filter_map(|(inside, value)| inside.then_some(value));
390        let mut sum = 0.0;
391        let mut maximum = f64::NEG_INFINITY;
392        let mut count = 0_usize;
393        for value in included {
394            if !value.is_finite() {
395                return Err(ValidationError::NonFiniteMaskValue {
396                    mask: self.name.clone(),
397                });
398            }
399            sum += value;
400            maximum = maximum.max(value);
401            count += 1;
402        }
403        if count == 0 {
404            return Err(ValidationError::EmptyMask(self.name.clone()));
405        }
406        Ok(MaskDoseSummary {
407            voxel_count: count,
408            mean: sum / count as f64,
409            maximum,
410        })
411    }
412
413    /// Centroid of the included voxels' centers in LPS millimetres.
414    ///
415    /// `geometry` must describe the grid this mask indexes; the mask must
416    /// select at least one voxel.
417    pub fn centroid_lps_mm(&self, geometry: &GridGeometry) -> Result<[f64; 3], ValidationError> {
418        let total = geometry.voxel_count()?;
419        if self.voxels.len() != total {
420            return Err(ValidationError::MaskVoxelCountMismatch {
421                name: self.name.clone(),
422                other: "geometry".into(),
423                expected: total,
424                actual: self.voxels.len(),
425            });
426        }
427        let nx = geometry.shape[0] as usize;
428        let ny = geometry.shape[1] as usize;
429        let mut centroid = [0.0_f64; 3];
430        let mut count = 0_usize;
431        for (index, inside) in self.voxels.iter().copied().enumerate() {
432            if !inside {
433                continue;
434            }
435            let voxel = [
436                (index % nx) as u32,
437                ((index % (nx * ny)) / nx) as u32,
438                (index / (nx * ny)) as u32,
439            ];
440            let center = geometry.voxel_center_lps_mm(voxel)?;
441            for axis in 0..3 {
442                centroid[axis] += center[axis];
443            }
444            count += 1;
445        }
446        if count == 0 {
447            return Err(ValidationError::EmptyMask(self.name.clone()));
448        }
449        for axis in &mut centroid {
450            *axis /= count as f64;
451        }
452        Ok(centroid)
453    }
454}
455
456/// Mean and maximum of a per-voxel quantity inside a `RegionMask`.
457#[derive(Debug, Clone, Copy, PartialEq)]
458pub struct MaskDoseSummary {
459    pub voxel_count: usize,
460    pub mean: f64,
461    pub maximum: f64,
462}
463
464/// Contract-id namespace emitted by current artifacts.
465pub const SCHEMA_PREFIX: &str = "openbnct.";
466
467/// Contract-id namespace used by artifacts written before the project was
468/// renamed from NCTForge to OpenBNCT. Such artifacts remain valid inputs and
469/// are normalized on read; new artifacts always emit [`SCHEMA_PREFIX`].
470pub const LEGACY_SCHEMA_PREFIX: &str = "nctforge.";
471
472/// Hyphenated tool/method-id namespace of the same pre-rename era (e.g.
473/// `nctforge-openmc-data-inspector/0.3.0`).
474pub const LEGACY_TOOL_PREFIX: &str = "nctforge-";
475
476/// Returns `id` with a legacy `nctforge.`/`nctforge-` contract or tool
477/// namespace replaced by the current `openbnct.`/`openbnct-` namespace. Ids
478/// without a legacy prefix are returned unchanged.
479pub fn normalize_contract_id(id: &str) -> String {
480    if let Some(rest) = id.strip_prefix(LEGACY_SCHEMA_PREFIX) {
481        return format!("{SCHEMA_PREFIX}{rest}");
482    }
483    if let Some(rest) = id.strip_prefix(LEGACY_TOOL_PREFIX) {
484        return format!("openbnct-{rest}");
485    }
486    id.to_string()
487}
488
489/// True when `actual` identifies the same contract as `expected`, tolerating
490/// the legacy `nctforge.` namespace on either side. Use this wherever an
491/// externally supplied artifact's `schema_version` is checked, so committed
492/// `nctforge.*` evidence remains readable.
493pub fn schema_matches(actual: &str, expected: &str) -> bool {
494    normalize_contract_id(actual) == normalize_contract_id(expected)
495}
496
497/// Serde `deserialize_with` for `schema_version` fields: accepts the legacy
498/// `nctforge.` namespace and normalizes it to `openbnct.` in memory, so
499/// pre-rename artifacts load into current contracts transparently.
500pub fn deserialize_contract_id<'de, D>(deserializer: D) -> Result<String, D::Error>
501where
502    D: serde::Deserializer<'de>,
503{
504    Ok(normalize_contract_id(&String::deserialize(deserializer)?))
505}
506
507/// Schema identifier carried by every `PhysicalDoseBundle`.
508pub const PHYSICAL_DOSE_BUNDLE_SCHEMA: &str = "openbnct.physical-dose-bundle/0.2.0";
509
510#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
511#[serde(deny_unknown_fields)]
512pub struct PhysicalDoseBundle {
513    #[serde(deserialize_with = "crate::deserialize_contract_id")]
514    pub schema_version: String,
515    pub case_id: String,
516    pub frame_of_reference_uid: Option<String>,
517    pub geometry: GridGeometry,
518    pub component_profile: ComponentProfileReference,
519    /// Material- and nuclear-data-specific neutron response curves.
520    pub response_set: ContentReference,
521    pub components: Vec<DoseVolume>,
522    /// A dedicated physical total, retained separately from component means.
523    pub physical_total: PhysicalTotalDoseVolume,
524    /// Identifier of the run manifest that binds inputs, engine, data, and logs.
525    pub provenance_id: String,
526}
527
528impl PhysicalDoseBundle {
529    pub fn validate(&self) -> Result<(), ValidationError> {
530        let voxel_count = self.geometry.voxel_count()?;
531        for (label, value) in [
532            ("schema_version", self.schema_version.as_str()),
533            ("case_id", self.case_id.as_str()),
534            ("provenance_id", self.provenance_id.as_str()),
535        ] {
536            if value.trim().is_empty() {
537                return Err(ValidationError::EmptyIdentifier(label));
538            }
539        }
540        self.component_profile
541            .validate()
542            .map_err(|_| ValidationError::InvalidContentReference("component_profile"))?;
543        self.response_set
544            .validate()
545            .map_err(|_| ValidationError::InvalidContentReference("response_set"))?;
546
547        let mut observed = BTreeSet::new();
548
549        for volume in &self.components {
550            if !observed.insert(volume.component) {
551                return Err(ValidationError::DuplicateComponent(volume.component));
552            }
553            if volume.values.len() != voxel_count {
554                return Err(ValidationError::DoseLength {
555                    component: volume.component,
556                    expected: voxel_count,
557                    actual: volume.values.len(),
558                });
559            }
560            if volume
561                .values
562                .iter()
563                .any(|value| !value.is_finite() || *value < 0.0)
564            {
565                return Err(ValidationError::InvalidDose(volume.component));
566            }
567            if volume.unit != self.physical_total.unit {
568                return Err(ValidationError::DoseUnitMismatch {
569                    component: volume.component,
570                    component_unit: volume.unit,
571                    total_unit: self.physical_total.unit,
572                });
573            }
574            if let Some(uncertainty) = &volume.absolute_standard_uncertainty {
575                if uncertainty.len() != voxel_count {
576                    return Err(ValidationError::UncertaintyLength {
577                        component: volume.component,
578                        expected: voxel_count,
579                        actual: uncertainty.len(),
580                    });
581                }
582                if uncertainty
583                    .iter()
584                    .any(|value| !value.is_finite() || *value < 0.0)
585                {
586                    return Err(ValidationError::InvalidUncertainty(volume.component));
587                }
588            }
589        }
590
591        for required in DoseComponent::REQUIRED {
592            if !observed.contains(&required) {
593                return Err(ValidationError::MissingComponent(required));
594            }
595        }
596
597        if self.physical_total.values.len() != voxel_count {
598            return Err(ValidationError::TotalDoseLength {
599                expected: voxel_count,
600                actual: self.physical_total.values.len(),
601            });
602        }
603        if self
604            .physical_total
605            .values
606            .iter()
607            .any(|value| !value.is_finite() || *value < 0.0)
608        {
609            return Err(ValidationError::InvalidTotalDose);
610        }
611        if let Some(uncertainty) = &self.physical_total.absolute_standard_uncertainty {
612            if uncertainty.len() != voxel_count {
613                return Err(ValidationError::TotalUncertaintyLength {
614                    expected: voxel_count,
615                    actual: uncertainty.len(),
616                });
617            }
618            if uncertainty
619                .iter()
620                .any(|value| !value.is_finite() || *value < 0.0)
621            {
622                return Err(ValidationError::InvalidTotalUncertainty);
623            }
624        }
625        match (
626            self.physical_total.absolute_standard_uncertainty.is_some(),
627            self.physical_total.uncertainty_method,
628        ) {
629            (false, TotalUncertaintyMethod::Unavailable)
630            | (true, TotalUncertaintyMethod::DedicatedEstimator)
631            | (true, TotalUncertaintyMethod::BatchCovariance) => {}
632            _ => return Err(ValidationError::InconsistentTotalUncertainty),
633        }
634
635        Ok(())
636    }
637}
638
639fn is_canonical_sha256(value: &str) -> bool {
640    value.len() == 64
641        && value
642            .bytes()
643            .all(|byte| byte.is_ascii_digit() || (b'a'..=b'f').contains(&byte))
644}
645
646#[derive(Debug, Error, PartialEq)]
647pub enum ValidationError {
648    #[error("voxel grid contains an empty dimension")]
649    EmptyGrid,
650    #[error("voxel spacing must be finite and greater than zero")]
651    InvalidSpacing,
652    #[error("voxel geometry contains a non-finite value")]
653    NonFiniteGeometry,
654    #[error("voxel direction matrix must be right-handed and orthonormal")]
655    InvalidDirection,
656    #[error("voxel count overflows the addressable platform size")]
657    GridTooLarge,
658    #[error("voxel {voxel:?} is outside grid shape {shape:?}")]
659    VoxelOutOfBounds { voxel: [u32; 3], shape: [u32; 3] },
660    #[error("physical dose component {0:?} is missing")]
661    MissingComponent(DoseComponent),
662    #[error("physical dose component {0:?} occurs more than once")]
663    DuplicateComponent(DoseComponent),
664    #[error("{component:?} has {actual} dose values; expected {expected}")]
665    DoseLength {
666        component: DoseComponent,
667        expected: usize,
668        actual: usize,
669    },
670    #[error("{0:?} contains a negative or non-finite physical dose")]
671    InvalidDose(DoseComponent),
672    #[error("{component:?} has {actual} uncertainty values; expected {expected}")]
673    UncertaintyLength {
674        component: DoseComponent,
675        expected: usize,
676        actual: usize,
677    },
678    #[error("{0:?} contains a negative or non-finite uncertainty")]
679    InvalidUncertainty(DoseComponent),
680    #[error("dose index {index} is outside volume length {length}")]
681    DoseIndexOutOfBounds { index: usize, length: usize },
682    #[error("required identifier {0} is empty")]
683    EmptyIdentifier(&'static str),
684    #[error("{0} must have a nonempty ID and canonical lowercase SHA-256 digest")]
685    InvalidContentReference(&'static str),
686    #[error("{component:?} uses {component_unit:?}, but the physical total uses {total_unit:?}")]
687    DoseUnitMismatch {
688        component: DoseComponent,
689        component_unit: DoseUnit,
690        total_unit: DoseUnit,
691    },
692    #[error("physical total has {actual} dose values; expected {expected}")]
693    TotalDoseLength { expected: usize, actual: usize },
694    #[error("physical total contains a negative or non-finite dose")]
695    InvalidTotalDose,
696    #[error("physical total has {actual} uncertainty values; expected {expected}")]
697    TotalUncertaintyLength { expected: usize, actual: usize },
698    #[error("physical total contains a negative or non-finite uncertainty")]
699    InvalidTotalUncertainty,
700    #[error("physical-total uncertainty and its method are inconsistent")]
701    InconsistentTotalUncertainty,
702    #[error(
703        "mask {name:?} has {actual} voxels but mask {other:?} has {expected}; masks must share one voxel count"
704    )]
705    MaskVoxelCountMismatch {
706        name: String,
707        other: String,
708        expected: usize,
709        actual: usize,
710    },
711    #[error("mask {0:?} selects no voxels")]
712    EmptyMask(String),
713    #[error("threshold window [{minimum}, {maximum}] must be finite with minimum <= maximum")]
714    InvalidThresholdWindow { minimum: f64, maximum: f64 },
715    #[error("mask {mask:?} covers {mask_voxels} voxels but the values array has {values} entries")]
716    MaskValuesLength {
717        mask: String,
718        mask_voxels: usize,
719        values: usize,
720    },
721    #[error("mask {mask:?} contains a non-finite value")]
722    NonFiniteMaskValue { mask: String },
723    #[error("dose selection {mask:?} contains a negative or non-finite value at voxel {index}")]
724    InvalidMaskedDose { mask: String, index: usize },
725    #[error("dose statistic {name} is invalid: {reason}")]
726    InvalidStatistic { name: &'static str, reason: String },
727}
728
729#[derive(Debug, Error, PartialEq, Eq)]
730pub enum ContentReferenceError {
731    #[error("content reference ID is empty")]
732    EmptyId,
733    #[error("content reference SHA-256 must be 64 lowercase hexadecimal characters")]
734    InvalidSha256,
735}
736
737#[cfg(test)]
738mod tests {
739    use super::*;
740
741    #[test]
742    fn legacy_nctforge_contract_ids_match_current_schemas() {
743        assert!(schema_matches(
744            "nctforge.physical-dose-bundle/0.2.0",
745            PHYSICAL_DOSE_BUNDLE_SCHEMA
746        ));
747        assert!(schema_matches(
748            PHYSICAL_DOSE_BUNDLE_SCHEMA,
749            "nctforge.physical-dose-bundle/0.2.0"
750        ));
751        assert!(!schema_matches(
752            "nctforge.physical-dose-bundle/0.1.0",
753            PHYSICAL_DOSE_BUNDLE_SCHEMA
754        ));
755        assert!(!schema_matches(
756            "other.physical-dose-bundle/0.2.0",
757            PHYSICAL_DOSE_BUNDLE_SCHEMA
758        ));
759        assert_eq!(
760            normalize_contract_id("nctforge.registration/0.1.0"),
761            "openbnct.registration/0.1.0"
762        );
763        assert_eq!(
764            normalize_contract_id("openbnct.registration/0.1.0"),
765            "openbnct.registration/0.1.0"
766        );
767    }
768
769    #[test]
770    fn legacy_schema_version_is_accepted_on_read() {
771        let mut bundle = valid_bundle();
772        bundle.schema_version = "nctforge.physical-dose-bundle/0.2.0".into();
773        assert!(bundle.validate().is_ok());
774    }
775
776    fn geometry() -> GridGeometry {
777        GridGeometry {
778            shape: [2, 2, 1],
779            spacing_mm: [1.0, 1.0, 2.0],
780            origin_mm: [0.0; 3],
781            direction: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
782        }
783    }
784
785    fn volume(component: DoseComponent) -> DoseVolume {
786        DoseVolume {
787            component,
788            unit: DoseUnit::GrayPerSourceParticle,
789            values: vec![1.0; 4],
790            absolute_standard_uncertainty: Some(vec![0.1; 4]),
791        }
792    }
793
794    fn valid_bundle() -> PhysicalDoseBundle {
795        PhysicalDoseBundle {
796            schema_version: PHYSICAL_DOSE_BUNDLE_SCHEMA.into(),
797            case_id: "synthetic".into(),
798            frame_of_reference_uid: None,
799            geometry: geometry(),
800            component_profile: ComponentProfileReference {
801                id: "openbnct.macroscopic-absorbed-dose.v1".into(),
802                sha256: "a".repeat(64),
803            },
804            response_set: ContentReference {
805                id: "openbnct.synthetic-response-set.v1".into(),
806                sha256: "b".repeat(64),
807            },
808            components: DoseComponent::REQUIRED.into_iter().map(volume).collect(),
809            physical_total: PhysicalTotalDoseVolume {
810                unit: DoseUnit::GrayPerSourceParticle,
811                values: vec![4.0; 4],
812                absolute_standard_uncertainty: Some(vec![0.2; 4]),
813                uncertainty_method: TotalUncertaintyMethod::DedicatedEstimator,
814            },
815            provenance_id: "manifest-sha256:synthetic".into(),
816        }
817    }
818
819    #[test]
820    fn rejects_incomplete_component_bundle() {
821        let mut bundle = valid_bundle();
822        bundle.components.clear();
823
824        assert_eq!(
825            bundle.validate(),
826            Err(ValidationError::MissingComponent(DoseComponent::Boron))
827        );
828    }
829
830    #[test]
831    fn rejects_non_orthonormal_grid_direction() {
832        let mut invalid = geometry();
833        invalid.direction[4] = 2.0;
834        assert_eq!(
835            invalid.voxel_count(),
836            Err(ValidationError::InvalidDirection)
837        );
838    }
839
840    #[test]
841    fn validates_complete_physical_dose_contract() {
842        assert_eq!(valid_bundle().validate(), Ok(()));
843    }
844
845    #[test]
846    fn serializes_only_canonical_component_names() {
847        assert_eq!(
848            serde_json::to_string(&DoseComponent::Hydrogen).unwrap(),
849            "\"hydrogen\""
850        );
851        assert!(serde_json::from_str::<DoseComponent>("\"hydrogen_recoil\"").is_err());
852    }
853
854    #[test]
855    fn derives_relative_uncertainty_but_not_for_zero_mean() {
856        let volume = DoseVolume {
857            component: DoseComponent::Boron,
858            unit: DoseUnit::Gray,
859            values: vec![0.0, 2.0],
860            absolute_standard_uncertainty: Some(vec![0.1, 0.2]),
861        };
862
863        assert_eq!(volume.relative_standard_uncertainty(0), Ok(None));
864        let relative = volume.relative_standard_uncertainty(1).unwrap().unwrap();
865        assert!((relative - 0.1).abs() < f64::EPSILON);
866    }
867
868    #[test]
869    fn rejects_mixed_component_and_total_units() {
870        let mut bundle = valid_bundle();
871        bundle.components[0].unit = DoseUnit::Gray;
872
873        assert_eq!(
874            bundle.validate(),
875            Err(ValidationError::DoseUnitMismatch {
876                component: DoseComponent::Boron,
877                component_unit: DoseUnit::Gray,
878                total_unit: DoseUnit::GrayPerSourceParticle,
879            })
880        );
881    }
882
883    #[test]
884    fn requires_consistent_total_uncertainty_state() {
885        let mut bundle = valid_bundle();
886        bundle.physical_total.uncertainty_method = TotalUncertaintyMethod::Unavailable;
887
888        assert_eq!(
889            bundle.validate(),
890            Err(ValidationError::InconsistentTotalUncertainty)
891        );
892    }
893
894    #[test]
895    fn rejects_noncanonical_component_profile_hash() {
896        let mut bundle = valid_bundle();
897        bundle.component_profile.sha256 = "A".repeat(64);
898
899        assert_eq!(
900            bundle.validate(),
901            Err(ValidationError::InvalidContentReference(
902                "component_profile"
903            ))
904        );
905    }
906
907    #[test]
908    fn rejects_noncanonical_response_set_hash() {
909        let mut bundle = valid_bundle();
910        bundle.response_set.sha256 = "B".repeat(64);
911
912        assert_eq!(
913            bundle.validate(),
914            Err(ValidationError::InvalidContentReference("response_set"))
915        );
916    }
917
918    fn mask(name: &str, voxels: &[bool]) -> RegionMask {
919        RegionMask {
920            name: name.into(),
921            voxels: voxels.to_vec(),
922        }
923    }
924
925    #[test]
926    fn mask_subtract_removes_shared_voxels() {
927        let organ = mask("ORGAN", &[true, true, true, true]);
928        let tumor = mask("TUMOR", &[false, true, true, false]);
929
930        let limited = organ.subtract(&tumor, "ORGAN-T").unwrap();
931
932        assert_eq!(limited.name, "ORGAN-T");
933        assert_eq!(limited.voxels, vec![true, false, false, true]);
934        assert_eq!(limited.included_voxel_count(), 2);
935    }
936
937    #[test]
938    fn mask_union_and_intersection() {
939        let a = mask("A", &[true, true, false, false]);
940        let b = mask("B", &[false, true, true, false]);
941
942        assert_eq!(
943            a.union(&b, "U").unwrap().voxels,
944            vec![true, true, true, false]
945        );
946        assert_eq!(
947            a.intersection(&b, "I").unwrap().voxels,
948            vec![false, true, false, false]
949        );
950    }
951
952    #[test]
953    fn mask_ops_reject_frame_mismatch() {
954        let a = mask("A", &[true; 4]);
955        let b = mask("B", &[true; 8]);
956
957        assert_eq!(
958            a.union(&b, "U"),
959            Err(ValidationError::MaskVoxelCountMismatch {
960                name: "A".into(),
961                other: "B".into(),
962                expected: 4,
963                actual: 8,
964            })
965        );
966    }
967
968    #[test]
969    fn mask_ops_reject_empty_result() {
970        let a = mask("A", &[true, true, false, false]);
971        let b = mask("B", &[false, false, true, true]);
972
973        assert_eq!(
974            a.intersection(&b, "I"),
975            Err(ValidationError::EmptyMask("I".into()))
976        );
977        assert_eq!(
978            a.subtract(&a.clone(), "E"),
979            Err(ValidationError::EmptyMask("E".into()))
980        );
981    }
982}
openbnct_core/lib.rs

openbnct_core/
lib.rs
