wsi_dicom_annotations/annotations/model/
geometry.rs

1use serde::{Deserialize, Serialize};
2
3use crate::{Error, Result};
4
5const MIN_POLYGON_AREA: f64 = 0.5;
6
7#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
8#[serde(deny_unknown_fields)]
9pub struct Point2 {
10    pub x: f64,
11    pub y: f64,
12}
13
14#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
15#[serde(deny_unknown_fields)]
16pub struct Point3 {
17    pub x: f64,
18    pub y: f64,
19    pub z: f64,
20}
21
22impl Point3 {
23    #[must_use]
24    pub const fn new(x: f64, y: f64, z: f64) -> Self {
25        Self { x, y, z }
26    }
27}
28
29impl Point2 {
30    #[must_use]
31    pub const fn new(x: f64, y: f64) -> Self {
32        Self { x, y }
33    }
34}
35
36#[derive(Debug, Clone, Copy, PartialEq, Eq)]
37pub enum AnnotationGraphicType {
38    Point,
39    Polygon,
40    Polyline,
41    Ellipse,
42    Rectangle,
43}
44
45impl AnnotationGraphicType {
46    #[must_use]
47    pub const fn dicom_value(self) -> &'static str {
48        match self {
49            Self::Point => "POINT",
50            Self::Polygon => "POLYGON",
51            Self::Polyline => "POLYLINE",
52            Self::Ellipse => "ELLIPSE",
53            Self::Rectangle => "RECTANGLE",
54        }
55    }
56
57    pub(crate) fn from_dicom(value: &str) -> Result<Self> {
58        match value.trim() {
59            "POINT" => Ok(Self::Point),
60            "POLYGON" => Ok(Self::Polygon),
61            "POLYLINE" => Ok(Self::Polyline),
62            "ELLIPSE" => Ok(Self::Ellipse),
63            "RECTANGLE" => Ok(Self::Rectangle),
64            other => Err(Error::Unsupported(format!(
65                "ANN graphic type {other:?} is not recognized"
66            ))),
67        }
68    }
69}
70
71#[derive(Debug, Clone, PartialEq)]
72pub enum AnnotationGeometry {
73    Points(Vec<Point2>),
74    Polygons(Vec<Vec<Point2>>),
75    ReadOnly {
76        graphic_type: AnnotationGraphicType,
77        coordinates: Vec<f64>,
78        primitive_point_indices: Vec<u32>,
79        coordinate_dimensions: usize,
80    },
81}
82
83impl AnnotationGeometry {
84    #[must_use]
85    pub const fn graphic_type(&self) -> AnnotationGraphicType {
86        match self {
87            Self::Points(_) => AnnotationGraphicType::Point,
88            Self::Polygons(_) => AnnotationGraphicType::Polygon,
89            Self::ReadOnly { graphic_type, .. } => *graphic_type,
90        }
91    }
92
93    #[must_use]
94    pub const fn editable(&self) -> bool {
95        matches!(self, Self::Points(_) | Self::Polygons(_))
96    }
97
98    pub(crate) fn annotation_count(&self) -> usize {
99        match self {
100            Self::Points(points) => points.len(),
101            Self::Polygons(polygons) => polygons.len(),
102            Self::ReadOnly {
103                graphic_type,
104                coordinates,
105                primitive_point_indices,
106                coordinate_dimensions,
107            } => match graphic_type {
108                AnnotationGraphicType::Point => coordinates.len() / coordinate_dimensions,
109                AnnotationGraphicType::Polygon | AnnotationGraphicType::Polyline => {
110                    primitive_point_indices.len()
111                }
112                AnnotationGraphicType::Ellipse | AnnotationGraphicType::Rectangle => {
113                    coordinates.len() / (coordinate_dimensions * 4)
114                }
115            },
116        }
117    }
118}
119
120pub(super) fn validate_points(points: impl Iterator<Item = Point2>) -> Result<()> {
121    if points
122        .into_iter()
123        .any(|point| !point.x.is_finite() || !point.y.is_finite())
124    {
125        return Err(Error::InvalidInput(
126            "annotation coordinates must be finite".into(),
127        ));
128    }
129    Ok(())
130}
131
132/// Validates a finite, non-self-intersecting, implicitly closed polygon ring.
133pub fn validate_polygon(points: &[Point2]) -> Result<()> {
134    if points.len() < 3 {
135        return Err(Error::InvalidInput(
136            "an annotation polygon needs at least three vertices".into(),
137        ));
138    }
139    validate_points(points.iter().copied())?;
140    if points.first() == points.last() {
141        return Err(Error::InvalidInput(
142            "ANN polygons close implicitly and must not repeat the first vertex".into(),
143        ));
144    }
145    if polygon_signed_area(points).abs() < MIN_POLYGON_AREA {
146        return Err(Error::InvalidInput(
147            "an annotation polygon has zero area".into(),
148        ));
149    }
150    if polygon_self_intersects(points) {
151        return Err(Error::InvalidInput(
152            "an annotation polygon crosses itself".into(),
153        ));
154    }
155    Ok(())
156}
157
158pub(crate) fn clockwise_polygon(points: &[Point2]) -> impl Iterator<Item = Point2> + '_ {
159    let reverse = polygon_signed_area(points) < 0.0;
160    (0..points.len()).map(move |index| {
161        let index = if reverse {
162            points.len() - 1 - index
163        } else {
164            index
165        };
166        points[index]
167    })
168}
169
170/// Reports whether `point` is inside `polygon` using an even-odd ray cast.
171///
172/// Boundary points are not assigned special semantics. Callers that require
173/// strict containment must separately reject boundary intersections.
174#[must_use]
175pub fn polygon_contains_point(polygon: &[Point2], point: Point2) -> bool {
176    let mut inside = false;
177    for (a, b) in polygon
178        .iter()
179        .zip(polygon.iter().cycle().skip(1))
180        .take(polygon.len())
181    {
182        if (a.y > point.y) != (b.y > point.y) {
183            let crossing_x = (b.x - a.x) * (point.y - a.y) / (b.y - a.y) + a.x;
184            if point.x < crossing_x {
185                inside = !inside;
186            }
187        }
188    }
189    inside
190}
191
192/// Reports whether any edges from two implicitly closed polygons intersect.
193#[must_use]
194pub fn polygon_boundaries_intersect(first: &[Point2], second: &[Point2]) -> bool {
195    first
196        .iter()
197        .zip(first.iter().cycle().skip(1))
198        .take(first.len())
199        .any(|(first_start, first_end)| {
200            second
201                .iter()
202                .zip(second.iter().cycle().skip(1))
203                .take(second.len())
204                .any(|(second_start, second_end)| {
205                    segments_intersect(*first_start, *first_end, *second_start, *second_end)
206                })
207        })
208}
209
210/// Returns the signed area of an implicitly closed polygon.
211///
212/// Positive values use the coordinate system's counterclockwise winding.
213#[must_use]
214pub fn polygon_signed_area(points: &[Point2]) -> f64 {
215    points
216        .iter()
217        .zip(points.iter().cycle().skip(1))
218        .take(points.len())
219        .map(|(a, b)| a.x * b.y - b.x * a.y)
220        .sum::<f64>()
221        * 0.5
222}
223
224/// Reports whether non-adjacent edges of an implicitly closed polygon intersect.
225#[must_use]
226pub fn polygon_self_intersects(points: &[Point2]) -> bool {
227    for first in 0..points.len() {
228        let first_next = (first + 1) % points.len();
229        for second in (first + 1)..points.len() {
230            let second_next = (second + 1) % points.len();
231            if first_next == second || second_next == first {
232                continue;
233            }
234            if segments_intersect(
235                points[first],
236                points[first_next],
237                points[second],
238                points[second_next],
239            ) {
240                return true;
241            }
242        }
243    }
244    false
245}
246
247fn segments_intersect(a1: Point2, a2: Point2, b1: Point2, b2: Point2) -> bool {
248    let o1 = orientation(a1, a2, b1);
249    let o2 = orientation(a1, a2, b2);
250    let o3 = orientation(b1, b2, a1);
251    let o4 = orientation(b1, b2, a2);
252    if ((o1 > 0.0 && o2 < 0.0) || (o1 < 0.0 && o2 > 0.0))
253        && ((o3 > 0.0 && o4 < 0.0) || (o3 < 0.0 && o4 > 0.0))
254    {
255        return true;
256    }
257    (o1.abs() <= f64::EPSILON && point_on_segment(b1, a1, a2))
258        || (o2.abs() <= f64::EPSILON && point_on_segment(b2, a1, a2))
259        || (o3.abs() <= f64::EPSILON && point_on_segment(a1, b1, b2))
260        || (o4.abs() <= f64::EPSILON && point_on_segment(a2, b1, b2))
261}
262
263fn point_on_segment(point: Point2, a: Point2, b: Point2) -> bool {
264    point.x >= a.x.min(b.x)
265        && point.x <= a.x.max(b.x)
266        && point.y >= a.y.min(b.y)
267        && point.y <= a.y.max(b.y)
268}
269
270fn orientation(a: Point2, b: Point2, c: Point2) -> f64 {
271    (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)
272}
wsi_dicom_annotations/annotations/model/geometry.rs

wsi_dicom_annotations/annotations/model/
geometry.rs
