dicom_viewer/dcm/
roi.rs

1//! ROI statistics + measurement units.
2
3use crate::dcm::pixel::RawImage;
4use crate::dcm::study::Instance;
5
6/// First-moment statistics over a region of interest.
7///
8/// Produced by [`rect_stats`] and [`ellipse_stats`]. Values are in the
9/// same units as the underlying [`RawImage::values`] — Hounsfield Units
10/// on CT (because the modality LUT has been applied) and raw rescaled
11/// pixel values otherwise.
12#[derive(Debug, Clone, Copy)]
13pub struct RoiStats {
14    /// Arithmetic mean of pixels inside the ROI.
15    pub mean: f64,
16    /// Sample standard deviation (population formula with `count` in the
17    /// denominator — adequate for a non-diagnostic overlay).
18    pub std: f64,
19    /// Minimum pixel value inside the ROI.
20    pub min: f64,
21    /// Maximum pixel value inside the ROI.
22    pub max: f64,
23    /// Number of pixels summed.
24    pub count: usize,
25}
26
27impl RoiStats {
28    /// Format the stats for overlay display, suffixing `HU` on CT.
29    pub fn label(&self, modality: &str) -> String {
30        let unit = if modality.eq_ignore_ascii_case("CT") {
31            " HU"
32        } else {
33            ""
34        };
35        format!(
36            "n={} μ={:.1}{u} σ={:.1} min={:.0} max={:.0}",
37            self.count,
38            self.mean,
39            self.std,
40            self.min,
41            self.max,
42            u = unit
43        )
44    }
45}
46
47fn bbox(p1: [f32; 2], p2: [f32; 2], w: u32, h: u32) -> (u32, u32, u32, u32) {
48    let x0 = p1[0].min(p2[0]).max(0.0).round() as i32;
49    let y0 = p1[1].min(p2[1]).max(0.0).round() as i32;
50    let x1 = p1[0].max(p2[0]).round() as i32;
51    let y1 = p1[1].max(p2[1]).round() as i32;
52    let x0 = x0.clamp(0, w as i32 - 1) as u32;
53    let y0 = y0.clamp(0, h as i32 - 1) as u32;
54    let x1 = x1.clamp(0, w as i32 - 1) as u32;
55    let y1 = y1.clamp(0, h as i32 - 1) as u32;
56    (x0, y0, x1.max(x0), y1.max(y0))
57}
58
59/// Compute [`RoiStats`] over the axis-aligned rectangle spanned by `p1`
60/// and `p2` (in displayed-image-pixel coordinates).
61///
62/// Returns `None` when the rectangle degenerates after clamping to the
63/// image bounds.
64pub fn rect_stats(raw: &RawImage, p1: [f32; 2], p2: [f32; 2]) -> Option<RoiStats> {
65    let (x0, y0, x1, y1) = bbox(p1, p2, raw.width, raw.height);
66    if x1 <= x0 || y1 <= y0 {
67        return None;
68    }
69    let mut sum = 0.0f64;
70    let mut sumsq = 0.0f64;
71    let mut mn = f64::INFINITY;
72    let mut mx = f64::NEG_INFINITY;
73    let mut count = 0usize;
74    let stride = raw.width as usize;
75    for y in y0..=y1 {
76        for x in x0..=x1 {
77            let v = raw.values[y as usize * stride + x as usize] as f64;
78            sum += v;
79            sumsq += v * v;
80            if v < mn {
81                mn = v;
82            }
83            if v > mx {
84                mx = v;
85            }
86            count += 1;
87        }
88    }
89    if count == 0 {
90        return None;
91    }
92    let mean = sum / count as f64;
93    let var = (sumsq / count as f64 - mean * mean).max(0.0);
94    Some(RoiStats {
95        mean,
96        std: var.sqrt(),
97        min: mn,
98        max: mx,
99        count,
100    })
101}
102
103/// Compute [`RoiStats`] over the ellipse inscribed in the rectangle
104/// spanned by `p1` and `p2` (in displayed-image-pixel coordinates).
105///
106/// Returns `None` when the bounding rectangle degenerates after clamping
107/// to the image bounds.
108pub fn ellipse_stats(raw: &RawImage, p1: [f32; 2], p2: [f32; 2]) -> Option<RoiStats> {
109    let (x0, y0, x1, y1) = bbox(p1, p2, raw.width, raw.height);
110    if x1 <= x0 || y1 <= y0 {
111        return None;
112    }
113    let cx = (x0 as f64 + x1 as f64) / 2.0;
114    let cy = (y0 as f64 + y1 as f64) / 2.0;
115    let rx = ((x1 - x0) as f64 / 2.0).max(0.5);
116    let ry = ((y1 - y0) as f64 / 2.0).max(0.5);
117    let mut sum = 0.0f64;
118    let mut sumsq = 0.0f64;
119    let mut mn = f64::INFINITY;
120    let mut mx = f64::NEG_INFINITY;
121    let mut count = 0usize;
122    let stride = raw.width as usize;
123    for y in y0..=y1 {
124        for x in x0..=x1 {
125            let dx = (x as f64 - cx) / rx;
126            let dy = (y as f64 - cy) / ry;
127            if dx * dx + dy * dy > 1.0 {
128                continue;
129            }
130            let v = raw.values[y as usize * stride + x as usize] as f64;
131            sum += v;
132            sumsq += v * v;
133            if v < mn {
134                mn = v;
135            }
136            if v > mx {
137                mx = v;
138            }
139            count += 1;
140        }
141    }
142    if count == 0 {
143        return None;
144    }
145    let mean = sum / count as f64;
146    let var = (sumsq / count as f64 - mean * mean).max(0.0);
147    Some(RoiStats {
148        mean,
149        std: var.sqrt(),
150        min: mn,
151        max: mx,
152        count,
153    })
154}
155
156/// Distance between two displayed-pixel points in millimetres if pixel
157/// spacing is known, else pixels.
158pub fn length_label(p1: [f32; 2], p2: [f32; 2], inst: &Instance, raw: &RawImage) -> String {
159    let dx_px = (p2[0] - p1[0]) as f64;
160    let dy_px = (p2[1] - p1[1]) as f64;
161    let scale = raw.display_scale.max(1) as f64;
162    if let Some((row_sp, col_sp)) = inst.pixel_spacing {
163        let dx = dx_px * col_sp * scale;
164        let dy = dy_px * row_sp * scale;
165        let mm = (dx * dx + dy * dy).sqrt();
166        format!("{:.1} mm", mm)
167    } else {
168        let px = (dx_px * dx_px + dy_px * dy_px).sqrt();
169        format!("{:.0} px", px)
170    }
171}
172
173/// Angle at vertex `v` between rays to `p1` and `p2`, in degrees.
174///
175/// Returns `0.0` when either ray has near-zero length.
176pub fn angle_deg(p1: [f32; 2], v: [f32; 2], p2: [f32; 2]) -> f64 {
177    let a = (p1[0] - v[0], p1[1] - v[1]);
178    let b = (p2[0] - v[0], p2[1] - v[1]);
179    let dot = (a.0 * b.0 + a.1 * b.1) as f64;
180    let la = ((a.0 * a.0 + a.1 * a.1) as f64).sqrt();
181    let lb = ((b.0 * b.0 + b.1 * b.1) as f64).sqrt();
182    if la < 1e-6 || lb < 1e-6 {
183        return 0.0;
184    }
185    (dot / (la * lb)).clamp(-1.0, 1.0).acos().to_degrees()
186}
dicom_viewer/dcm/roi.rs

dicom_viewer/dcm/
roi.rs
