wsi_dicom_annotations/annotations/parametric_map/
document.rs

1use std::path::Path;
2
3use sha2::{Digest, Sha256};
4
5use crate::{
6    DerivedObjectProducer, DicomAnnotationContext, Error, InteroperabilityDiagnostic, Result,
7};
8
9use super::super::semantic_digest::{finish, update_algorithm, update_code, update_text};
10use super::frame::{encode_frame, FrameKey, PixelPrecision};
11use super::geometry::RasterGeometry;
12use super::profile::{RasterDType, RasterOutputPrecision};
13use super::{NormalizedRaster, RasterChannelSelection, RasterDescriptor, RasterProfile};
14
15pub struct ParametricMapDocument {
16    pub(super) source: DicomAnnotationContext,
17    pub(super) profile: RasterProfile,
18    pub(super) producer: DerivedObjectProducer,
19    pub(super) raster: NormalizedRaster,
20    pub(super) selected_channels: Vec<usize>,
21    pub(super) descriptor: RasterDescriptor,
22    pub(super) geometry: RasterGeometry,
23    pub(super) precision: PixelPrecision,
24    pub(super) frame_count: u32,
25    pub(super) dense: bool,
26    pub(super) channel_ranges: Vec<ValueRange>,
27    pub(super) pixel_digest: [u8; 32],
28    semantic_digest: String,
29    diagnostics: Vec<InteroperabilityDiagnostic>,
30}
31
32#[derive(Debug, Clone, Copy)]
33pub(super) struct ValueRange {
34    pub(super) minimum: f64,
35    pub(super) maximum: f64,
36}
37
38impl ParametricMapDocument {
39    pub fn open(
40        source: DicomAnnotationContext,
41        canonical_source: DicomAnnotationContext,
42        profile: RasterProfile,
43        raster_path: impl AsRef<Path>,
44        channel_selection: RasterChannelSelection,
45    ) -> Result<Self> {
46        let selected_channels = profile.select_channels(channel_selection)?;
47        let raster_path = raster_path.as_ref().to_path_buf();
48        let raster = NormalizedRaster::open(&profile, &raster_path)?;
49        let descriptor = raster.descriptor();
50        validate_descriptor(descriptor)?;
51        let precision = output_precision(&profile);
52        let geometry = RasterGeometry::resolve(&source, &canonical_source, &profile, descriptor)?;
53        let mut document = Self {
54            source,
55            profile,
56            producer: DerivedObjectProducer::library_default(9401, "WSI parametric maps"),
57            raster,
58            selected_channels,
59            descriptor,
60            geometry,
61            precision,
62            frame_count: 0,
63            dense: false,
64            channel_ranges: Vec::new(),
65            pixel_digest: [0; 32],
66            semantic_digest: String::new(),
67            diagnostics: Vec::new(),
68        };
69        document.scan()?;
70        Ok(document)
71    }
72
73    /// Replaces the neutral library identity with caller-owned producer metadata.
74    #[must_use]
75    pub fn with_producer(mut self, producer: DerivedObjectProducer) -> Self {
76        self.producer = producer;
77        self
78    }
79
80    #[must_use]
81    pub fn producer(&self) -> &DerivedObjectProducer {
82        &self.producer
83    }
84
85    #[must_use]
86    pub fn frame_count(&self) -> u32 {
87        self.frame_count
88    }
89
90    #[must_use]
91    pub fn semantic_digest(&self) -> &str {
92        &self.semantic_digest
93    }
94
95    #[must_use]
96    pub fn selected_channel_count(&self) -> usize {
97        self.selected_channels.len()
98    }
99
100    #[must_use]
101    pub fn diagnostics(&self) -> &[InteroperabilityDiagnostic] {
102        &self.diagnostics
103    }
104
105    fn scan(&mut self) -> Result<()> {
106        let mut digest = Sha256::new();
107        let mut ranges = vec![MutableRange::default(); self.selected_channels.len()];
108        let mut frame_count = 0_u64;
109        let mut missing_frame_count = 0_u64;
110        let mut missing_sample_count = 0_u64;
111        let mut padded_sample_count = 0_u64;
112        let tiles_across = self.descriptor.width.div_ceil(self.descriptor.tile_width);
113        let tiles_down = self.descriptor.height.div_ceil(self.descriptor.tile_height);
114        for (channel_ordinal, channel) in self.selected_channels.iter().copied().enumerate() {
115            let channel_ordinal = u32::try_from(channel_ordinal + 1)
116                .map_err(|_| Error::InvalidInput("raster channel count exceeds DICOM UL".into()))?;
117            for tile_row in 0..tiles_down {
118                for tile_column in 0..tiles_across {
119                    let key = FrameKey {
120                        channel,
121                        channel_ordinal,
122                        tile_row,
123                        tile_column,
124                    };
125                    let frame = encode_frame(&self.raster, self.descriptor, key)?;
126                    validate_encoded_precision(self.precision, self.descriptor, frame.bytes.len())?;
127                    missing_sample_count = missing_sample_count
128                        .checked_add(frame.missing_sample_count)
129                        .ok_or_else(|| {
130                            Error::InvalidInput("missing raster sample count overflows".into())
131                        })?;
132                    padded_sample_count = padded_sample_count
133                        .checked_add(frame.padded_sample_count)
134                        .ok_or_else(|| {
135                            Error::InvalidInput("padded raster sample count overflows".into())
136                        })?;
137                    if frame.all_missing {
138                        missing_frame_count += 1;
139                        continue;
140                    }
141                    digest.update(&frame.bytes);
142                    ranges[channel_ordinal as usize - 1]
143                        .observe(frame.finite_minimum, frame.finite_maximum);
144                    frame_count += 1;
145                }
146            }
147        }
148        if frame_count == 0 {
149            return Err(Error::InvalidInput(
150                "selected raster channels contain no finite samples".into(),
151            ));
152        }
153        self.frame_count = u32::try_from(frame_count).map_err(|_| {
154            Error::InvalidInput(
155                "selected raster contains more frames than Concatenation Frame Offset Number can represent"
156                    .into(),
157            )
158        })?;
159        self.dense = self.selected_channels.len() == 1 && missing_frame_count == 0;
160        self.channel_ranges = ranges
161            .into_iter()
162            .map(MutableRange::finish)
163            .collect::<Result<Vec<_>>>()?;
164        self.pixel_digest = digest.finalize().into();
165        self.semantic_digest = semantic_digest(self);
166        self.diagnostics = build_diagnostics(
167            &self.profile,
168            self.raster.normalizations(),
169            missing_sample_count,
170            padded_sample_count,
171            missing_frame_count,
172        );
173        Ok(())
174    }
175
176    pub(super) fn output_frames(&self) -> OutputFrameCursor<'_> {
177        OutputFrameCursor {
178            document: self,
179            channel_ordinal: 0,
180            tile_index: 0,
181        }
182    }
183}
184
185fn build_diagnostics(
186    profile: &RasterProfile,
187    source_normalizations: &[InteroperabilityDiagnostic],
188    missing_samples: u64,
189    padded_samples: u64,
190    omitted_frames: u64,
191) -> Vec<InteroperabilityDiagnostic> {
192    let mut diagnostics = source_normalizations.to_vec();
193    if profile.dtype.is_integer() {
194        diagnostics.push(InteroperabilityDiagnostic::normalized(
195            "RASTER_INTEGER_SCALED",
196            "$.raster.integer_scaling",
197            "integer samples were transformed with the explicitly declared slope, intercept, missing sentinel, and output precision",
198        ));
199    }
200    if missing_samples > 0 {
201        diagnostics.push(InteroperabilityDiagnostic::normalized(
202            "RASTER_NAN_CANONICALIZED",
203            "$.raster.values",
204            format!(
205                "canonicalized {missing_samples} missing samples to the selected quiet-NaN payload"
206            ),
207        ));
208    }
209    if padded_samples > 0 {
210        diagnostics.push(InteroperabilityDiagnostic::normalized(
211            "PM_EDGE_TILE_PADDED",
212            "$.raster.edge_tiles",
213            format!("padded {padded_samples} edge-tile samples with the canonical quiet NaN"),
214        ));
215    }
216    if omitted_frames > 0 {
217        diagnostics.push(InteroperabilityDiagnostic::normalized(
218            "PM_ALL_MISSING_TILES_OMITTED",
219            "$.raster.tiles",
220            format!("omitted {omitted_frames} all-missing tiles and encoded TILED_SPARSE"),
221        ));
222    }
223    diagnostics
224}
225
226pub(super) struct OutputFrameCursor<'a> {
227    document: &'a ParametricMapDocument,
228    channel_ordinal: usize,
229    tile_index: u64,
230}
231
232impl OutputFrameCursor<'_> {
233    pub(super) fn next(&mut self) -> Result<Option<super::frame::EncodedFrame>> {
234        let descriptor = self.document.descriptor;
235        let tiles_across = descriptor.width.div_ceil(descriptor.tile_width);
236        let tiles_down = descriptor.height.div_ceil(descriptor.tile_height);
237        let tiles_per_channel = u64::from(tiles_across) * u64::from(tiles_down);
238        while self.channel_ordinal < self.document.selected_channels.len() {
239            if self.tile_index == tiles_per_channel {
240                self.channel_ordinal += 1;
241                self.tile_index = 0;
242                continue;
243            }
244            let key = FrameKey {
245                channel: self.document.selected_channels[self.channel_ordinal],
246                channel_ordinal: self.channel_ordinal as u32 + 1,
247                tile_row: (self.tile_index / u64::from(tiles_across)) as u32,
248                tile_column: (self.tile_index % u64::from(tiles_across)) as u32,
249            };
250            self.tile_index += 1;
251            let frame = encode_frame(&self.document.raster, descriptor, key)?;
252            if !frame.all_missing {
253                return Ok(Some(frame));
254            }
255        }
256        Ok(None)
257    }
258}
259
260#[derive(Debug, Clone, Copy, Default)]
261struct MutableRange {
262    minimum: Option<f64>,
263    maximum: Option<f64>,
264}
265
266impl MutableRange {
267    fn observe(&mut self, minimum: Option<f64>, maximum: Option<f64>) {
268        if let Some(value) = minimum {
269            self.minimum = Some(self.minimum.map_or(value, |current| current.min(value)));
270        }
271        if let Some(value) = maximum {
272            self.maximum = Some(self.maximum.map_or(value, |current| current.max(value)));
273        }
274    }
275
276    fn finish(self) -> Result<ValueRange> {
277        match (self.minimum, self.maximum) {
278            (Some(minimum), Some(maximum)) => Ok(ValueRange { minimum, maximum }),
279            _ => Err(Error::InvalidInput(
280                "a selected raster channel contains no finite samples".into(),
281            )),
282        }
283    }
284}
285
286fn validate_descriptor(descriptor: RasterDescriptor) -> Result<()> {
287    if descriptor.tile_height == 0
288        || descriptor.tile_width == 0
289        || descriptor.tile_height > u32::from(u16::MAX)
290        || descriptor.tile_width > u32::from(u16::MAX)
291    {
292        return Err(Error::InvalidInput(
293            "PM tile rows and columns must each fit a nonzero DICOM US value".into(),
294        ));
295    }
296    Ok(())
297}
298
299fn output_precision(profile: &RasterProfile) -> PixelPrecision {
300    match profile.dtype {
301        RasterDType::Float32 => PixelPrecision::Float32,
302        RasterDType::Float64 => PixelPrecision::Float64,
303        _ => match profile
304            .integer_scaling
305            .as_ref()
306            .map(|scaling| scaling.output_precision)
307        {
308            Some(RasterOutputPrecision::Float64) => PixelPrecision::Float64,
309            Some(RasterOutputPrecision::Float32) | None => PixelPrecision::Float32,
310        },
311    }
312}
313
314fn validate_encoded_precision(
315    precision: PixelPrecision,
316    descriptor: RasterDescriptor,
317    actual: usize,
318) -> Result<()> {
319    let expected = super::frame::frame_byte_length(descriptor, precision)?;
320    if actual != expected {
321        return Err(Error::InvalidInput(format!(
322            "raster adapter produced {actual} bytes per frame, expected {expected} for the profile precision"
323        )));
324    }
325    Ok(())
326}
327
328fn semantic_digest(document: &ParametricMapDocument) -> String {
329    let mut digest = Sha256::new();
330    digest.update(b"dicom-parametric-map-v1\0");
331    update_text(&mut digest, document.source.study_instance_uid());
332    update_text(&mut digest, document.source.sop_instance_uid());
333    update_text(
334        &mut digest,
335        document.source.frame_of_reference_uid().unwrap_or(""),
336    );
337    for value in [
338        document.descriptor.height,
339        document.descriptor.width,
340        document.descriptor.tile_height,
341        document.descriptor.tile_width,
342    ] {
343        digest.update(value.to_le_bytes());
344    }
345    for value in document
346        .geometry
347        .origin
348        .into_iter()
349        .chain(document.geometry.orientation)
350        .chain(document.geometry.pixel_spacing)
351    {
352        digest.update(value.to_bits().to_le_bytes());
353    }
354    for channel in &document.selected_channels {
355        let channel = &document.profile.channels[*channel];
356        update_text(&mut digest, &channel.name);
357        update_code(&mut digest, &channel.quantity);
358        update_code(&mut digest, &channel.unit);
359    }
360    update_algorithm(&mut digest, &document.profile.algorithm);
361    digest.update(document.pixel_digest);
362    finish(digest)
363}
wsi_dicom_annotations/annotations/parametric_map/document.rs

wsi_dicom_annotations/annotations/parametric_map/
document.rs
