1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
use std::{cmp::Ordering, ops::Deref};

use super::SweepPoint;
use crate::{
    line_intersection::line_intersection, Coord, GeoFloat, GeoNum, Kernel, Line, LineIntersection,
    Orientation,
};

/// Either a line segment or a point.
///
/// The coordinates are ordered (see [`SweepPoint`]) and a line
/// segment must have distinct points (use the `Point` variant if the
/// coordinates are the equal).
#[derive(Clone, Copy)]
pub enum LineOrPoint<T: GeoNum> {
    Point(SweepPoint<T>),
    Line {
        left: SweepPoint<T>,
        right: SweepPoint<T>,
    },
}

impl<T: GeoNum> std::fmt::Debug for LineOrPoint<T> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            LineOrPoint::Point(p) => f.debug_tuple("Pt").field(&p.x_y()).finish(),
            LineOrPoint::Line { left, right } => f
                .debug_tuple("LPt")
                .field(&left.x_y())
                .field(&right.x_y())
                .finish(),
        }
    }
}

impl<T: GeoNum> From<SweepPoint<T>> for LineOrPoint<T> {
    fn from(pt: SweepPoint<T>) -> Self {
        Self::Point(pt)
    }
}

impl<T: GeoNum> From<(SweepPoint<T>, SweepPoint<T>)> for LineOrPoint<T> {
    fn from((start, end): (SweepPoint<T>, SweepPoint<T>)) -> Self {
        match start.cmp(&end) {
            Ordering::Less => Self::Line {
                left: start,
                right: end,
            },
            Ordering::Equal => Self::Point(start),
            Ordering::Greater => Self::Line {
                left: end,
                right: start,
            },
        }
    }
}

/// Convert from a [`Line`] ensuring end point ordering.
impl<T: GeoNum> From<Line<T>> for LineOrPoint<T> {
    fn from(l: Line<T>) -> Self {
        let start: SweepPoint<T> = l.start.into();
        let end = l.end.into();
        (start, end).into()
    }
}

/// Convert from a [`Coord`]
impl<T: GeoNum> From<Coord<T>> for LineOrPoint<T> {
    fn from(c: Coord<T>) -> Self {
        Self::Point(c.into())
    }
}

impl<T: GeoNum> LineOrPoint<T> {
    /// Checks if the variant is a line.
    #[inline]
    pub fn is_line(&self) -> bool {
        matches!(self, Self::Line { .. })
    }

    /// Return a [`Line`] representation of self.
    #[inline]
    pub fn line(&self) -> Line<T> {
        match self {
            LineOrPoint::Point(p) => Line::new(**p, **p),
            LineOrPoint::Line { left, right } => Line::new(**left, **right),
        }
    }

    #[inline]
    pub fn left(&self) -> SweepPoint<T> {
        match self {
            LineOrPoint::Point(p) => *p,
            LineOrPoint::Line { left, .. } => *left,
        }
    }

    #[inline]
    pub fn right(&self) -> SweepPoint<T> {
        match self {
            LineOrPoint::Point(p) => *p,
            LineOrPoint::Line { right, .. } => *right,
        }
    }

    #[cfg(test)]
    pub fn coords_equal(&self, other: &LineOrPoint<T>) -> bool {
        self.is_line() == other.is_line() && self.end_points() == other.end_points()
    }

    #[inline]
    pub fn end_points(&self) -> (SweepPoint<T>, SweepPoint<T>) {
        match self {
            LineOrPoint::Point(p) => (*p, *p),
            LineOrPoint::Line { left, right } => (*left, *right),
        }
    }

    pub fn new(left: SweepPoint<T>, right: SweepPoint<T>) -> Self {
        if left == right {
            Self::Point(left)
        } else {
            Self::Line { left, right }
        }
    }

    pub fn orient2d(&self, other: Coord<T>) -> Orientation {
        let (left, right) = match self {
            LineOrPoint::Point(p) => (**p, **p),
            LineOrPoint::Line { left, right } => (**left, **right),
        };
        T::Ker::orient2d(left, right, other)
    }
}

/// Equality based on ordering defined for segments as per algorithm.
impl<T: GeoNum> PartialEq for LineOrPoint<T> {
    #[inline]
    fn eq(&self, other: &Self) -> bool {
        self.partial_cmp(other) == Some(Ordering::Equal)
    }
}

/// Ordering defined for segments as per algorithm.
///
/// Requires the following conditions:
///
/// 1. If comparing two lines, both the left ends must be strictly
/// smaller than both right ends.
///
/// 2. A point is treated as a infinitesimal small vertical segment
/// centered at its coordinates.
impl<T: GeoNum> PartialOrd for LineOrPoint<T> {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        match (self, other) {
            (LineOrPoint::Point(p), LineOrPoint::Point(o)) => {
                if p == o {
                    Some(Ordering::Equal)
                } else {
                    // Unequal points do not satisfy pre-condition and
                    // can't be ordered.
                    None
                }
            }
            (LineOrPoint::Point(_), LineOrPoint::Line { .. }) => {
                other.partial_cmp(self).map(Ordering::reverse)
            }
            (LineOrPoint::Line { left, right }, LineOrPoint::Point(p)) => {
                if p > right || left > p {
                    return None;
                }
                Some(
                    T::Ker::orient2d(**left, **right, **p)
                        .as_ordering()
                        .then(Ordering::Greater),
                )
            }
            (
                LineOrPoint::Line {
                    left: left_a,
                    right: right_a,
                },
                LineOrPoint::Line {
                    left: left_b,
                    right: right_b,
                },
            ) => {
                if left_a > left_b {
                    return other.partial_cmp(self).map(Ordering::reverse);
                }
                if left_a >= right_b || left_b >= right_a {
                    return None;
                }

                // Assertion: p1 <= p2
                // Assertion: pi < q_j
                Some(
                    T::Ker::orient2d(**left_a, **right_a, **left_b)
                        .as_ordering()
                        .then_with(|| {
                            T::Ker::orient2d(**left_a, **right_a, **right_b).as_ordering()
                        }),
                )
            }
        }
    }
}

impl<T: GeoFloat> LineOrPoint<T> {
    /// Intersect a line with self and return a point, a overlapping segment or `None`.
    ///
    /// The `other` argument must be a line variant (debug builds will panic otherwise).
    pub fn intersect_line(&self, other: &Self) -> Option<Self> {
        debug_assert!(other.is_line(), "tried to intersect with a point variant!");

        let line = other.line();
        match self {
            LineOrPoint::Point(p) => {
                use crate::Intersects;
                if line.intersects(&**p) {
                    Some(*self)
                } else {
                    None
                }
            }
            LineOrPoint::Line { left, right } => {
                line_intersection(self.line(), line).map(|l| match l {
                    LineIntersection::SinglePoint {
                        intersection,
                        is_proper,
                    } => {
                        let mut pt = intersection;
                        if is_proper && (&pt == left.deref()) {
                            if left.x == right.x {
                                pt.y = pt.y.next_after(T::infinity());
                            } else {
                                pt.x = pt.x.next_after(T::infinity());
                            }
                        }
                        pt.into()
                    }
                    LineIntersection::Collinear { intersection } => intersection.into(),
                })
            }
        }
    }

    pub fn intersect_line_ordered(&self, other: &Self) -> Option<Self> {
        let ord = self.partial_cmp(other);
        match self.intersect_line(other) {
            Some(Self::Point(p)) => {
                // NOTE: A key issue with using non-exact numbers (f64, etc.) in
                // this algo. is that line-intersection may return
                // counter-intuitive points.
                //
                // Specifically, this causes two issues:
                //
                // 1. The point of intersection r lies between the end-points in
                // the lexicographic ordering. However, with finite repr., the
                // line (1, 1) - (1 + eps, -1), where eps is ulp(1), does not
                // admit r that lies between the end-points. Further, the
                // end-points may be a very bad approximation to the actual
                // intersection points (eg. intersect with x-axis).
                //
                // We detect and force r to be greater than both end-points; the
                // other case is not easy to handle as the sweep has already
                // progressed to a p strictly > r already.
                //
                // 2. The more severe issue is that in general r may not lie
                // exactly on the line. Thus, with the segment stored on the
                // active-segments tree (B-Tree / Splay), this may in adverse
                // cases, cause the ordering between the segments to be
                // incorrect, hence invalidating the segments. This is not easy
                // to correct without a intrusive data-structure built
                // specifically for this algo., that can track the neighbors of
                // tree-nodes, and fix / report this issue. The crate
                // `btree-slab` seems like a great starting point.
                let (mut x, y) = p.x_y();

                let c = self.left();
                if x == c.x && y < c.y {
                    x = x.next_after(T::infinity());
                }

                let p = Coord { x, y }.into();
                debug_assert!(
                    p >= self.left(),
                    "line intersection before first line: {p:?}\n\tLine({lp1:?} - {lp2:?}) X Line({lp3:?} - {lp4:?})",
                    lp1 = self.left(),
                    lp2 = self.right(),
                    lp3 = other.left(),
                    lp4 = other.right(),
                );
                debug_assert!(
                    p >= other.left(),
                    "line intersection before second line: {p:?}\n\tLine({lp1:?} - {lp2:?}) X Line({lp3:?} - {lp4:?})",
                    lp1 = self.left(),
                    lp2 = self.right(),
                    lp3 = other.left(),
                    lp4 = other.right(),
                );

                if let Some(ord) = ord {
                    let l1 = LineOrPoint::from((self.left(), p));
                    let l2 = LineOrPoint::from((other.left(), p));
                    let cmp = l1.partial_cmp(&l2).unwrap();
                    if l1.is_line() && l2.is_line() && cmp.then(ord) != ord {
                        debug!(
                            "ordering changed by intersection: {l1:?} {ord:?} {l2:?}",
                            l1 = self,
                            l2 = other
                        );
                        debug!("\tparts: {l1:?}, {l2:?}");
                        debug!("\tintersection: {p:?} {cmp:?}");

                        // RM: This is a complicated intersection that is changing the ordering.
                        // Heuristic: approximate with a trivial intersection point that preserves the topology.
                        return Some(if self.left() > other.left() {
                            self.left().into()
                        } else {
                            other.left().into()
                        });
                    }
                }
                Some(Self::Point(p))
            }
            e => e,
        }
    }
}

#[cfg(test)]
mod tests {
    use std::cmp::Ordering;

    use geo_types::{Coord, LineString};
    use wkt::ToWkt;

    use crate::{GeoFloat, GeoNum, Kernel};

    use super::LineOrPoint;

    // Used for debugging sweep fp issues
    #[test]
    #[ignore]
    fn check_ordering() {
        let pt_7 = Coord::from((-32.57812499999999, 241.33427773853316));
        let pt_8 = Coord::from((-36.11348070978957, 237.7989220287436));
        let pt_13 = Coord::from((-25.507080078124993, 248.40532266040816));
        let pt_14 = Coord::from((-36.48784219165816, 237.424560546875));
        let _pt_15 = Coord::from((4.4929199218750036, 196.44379843334184));
        let pt_16 = Coord::from((-36.048578439260666, 237.8638242992725));
        let pt_17 = Coord::from((3.545624214480127, 197.39109414073673));

        fn check_isection<T: GeoFloat>(abcd: [Coord<T>; 4]) -> Option<LineOrPoint<T>> {
            let l1 = LineOrPoint::from((abcd[0].into(), abcd[1].into()));
            let l2 = LineOrPoint::from((abcd[2].into(), abcd[3].into()));
            l1.intersect_line_ordered(&l2)
        }
        fn check_lines<T: GeoNum>(abcd: [Coord<T>; 4]) -> Ordering {
            let l1 = LineOrPoint::from((abcd[0].into(), abcd[1].into()));
            let l2 = LineOrPoint::from((abcd[2].into(), abcd[3].into()));
            l1.partial_cmp(&l2).unwrap()
        }

        eprintln!(
            "(14-17) {cmp:?} (14-16)",
            cmp = check_lines([pt_14, pt_17, pt_14, pt_16])
        );
        eprintln!(
            "(8-16) {cmp:?} (14-16)",
            cmp = check_lines([pt_8, pt_16, pt_14, pt_16]),
        );
        eprintln!(
            "(8-7) {cmp:?} (14-16)",
            cmp = check_lines([pt_8, pt_7, pt_14, pt_16]),
        );
        eprintln!(
            "(8-7) {cmp:?} (14-13)",
            cmp = check_lines([pt_8, pt_7, pt_14, pt_13]),
        );
        eprintln!(
            "(8-7) {isect:?} (14-13)",
            isect = check_isection([pt_8, pt_7, pt_14, pt_13]),
        );
        let l87 = LineString::new(vec![pt_8, pt_16, pt_7]);
        let lo = LineString::new(vec![pt_14, pt_16, pt_13]);
        eprintln!("l1: {}", l87.to_wkt());
        eprintln!("lo: {}", lo.to_wkt());

        eprintln!(
            "pred: {:?}",
            <f64 as GeoNum>::Ker::orient2d(pt_8, pt_7, pt_17)
        );
        eprintln!(
            "pred: {:?}",
            <f64 as GeoNum>::Ker::orient2d(pt_8, pt_14, pt_16)
        );
    }
}