linestring 0.16.0

// SPDX-License-Identifier: MIT OR Apache-2.0
// Copyright (c) 2021,2023 lacklustr@protonmail.com https://github.com/eadf

// This file is part of the linestring crate.

/*
Copyright (c) 2021,2023 lacklustr@protonmail.com https://github.com/eadf

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

or

Copyright 2021,2023 lacklustr@protonmail.com https://github.com/eadf

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

mod impls;

use crate::{LinestringError, linestring_2d};
use itertools::Itertools;
use std::{collections, fs, io, io::Write, path::Path};
use vector_traits::{
    approx::{AbsDiffEq, UlpsEq},
    num_traits::{AsPrimitive, real::Real},
    prelude::{GenericVector3, *},
};

pub(crate) mod simplify;
#[cfg(test)]
mod tests;

/// A 3d line
#[derive(PartialEq, Copy, Clone)]
pub struct Line3<T: GenericVector3> {
    pub start: T,
    pub end: T,
}

impl<T: GenericVector3> Line3<T> {
    pub fn new(start: T, end: T) -> Self {
        Self { start, end }
    }

    /// returns (area of a triangle)²*4
    fn triangle_area_squared_times_4(p1: &T, p2: &T, p3: &T) -> T::Scalar {
        let v1_x = p1.x() - p2.x();
        let v1_y = p1.y() - p2.y();
        let v1_z = p1.z() - p2.z();

        let v2_x = p3.x() - p2.x();
        let v2_y = p3.y() - p2.y();
        let v2_z = p3.z() - p2.z();

        let x = v1_y * v2_z - v2_y * v1_z;
        let y = v1_x * v2_z - v2_x * v1_z;
        let z = v1_x * v2_y - v2_x * v1_y;
        x * x + y * y + z * z
    }

    /// Copy this Line3 into a Line2, populating the axes defined by 'plane'
    /// An axis will always try to keep its position (e.g. y goes to y if possible).
    /// That way the operation is reversible (in regard to axis positions).
    pub fn copy_to_2d(&self, plane: Plane) -> linestring_2d::Line2<T::Vector2> {
        linestring_2d::Line2::new(
            plane.point_to_2d::<T>(self.start),
            plane.point_to_2d::<T>(self.end),
        )
    }

    pub fn aabb(&self) -> T::Aabb {
        let mut aabb = T::Aabb::from_point(self.start);
        aabb.add_point(self.end);
        aabb
    }

    /// Apply an operation over each coordinate.
    /// Useful when you want to round the values of the coordinates.
    pub fn apply<F>(&mut self, f: &F)
    where
        F: Fn(T) -> T,
    {
        self.start = f(self.start);
        self.end = f(self.end);
    }

    /// Return a new line containing the transformed points
    pub fn apply_new<F: Fn(T) -> T>(self, f: F) -> Self {
        Self {
            start: f(self.start),
            end: f(self.end),
        }
    }
}

struct PriorityDistance<T: GenericVector3> {
    key: T::Scalar,
    index: usize,
}

#[doc(hidden)]
pub struct WindowIterator<'a, T: GenericVector3>(std::slice::Windows<'a, T>);

impl<T: GenericVector3> WindowIterator<'_, T> {
    // TODO: add is_empty once the unstable feature "exact_size_is_empty" has transitioned
    #[inline(always)]
    pub fn is_empty(&self) -> bool {
        self.0.len() == 0
    }

    #[inline(always)]
    pub fn len(&self) -> usize {
        self.0.len()
    }
}

#[doc(hidden)]
pub struct ChunkIterator<'a, T: GenericVector3>(std::slice::ChunksExact<'a, T>);

impl<'a, T: GenericVector3> ChunkIterator<'a, T> {
    // TODO: add is_empty once the unstable feature "exact_size_is_empty" has transitioned
    #[inline(always)]
    pub fn is_empty(&self) -> bool {
        self.0.len() == 0
    }

    #[inline(always)]
    pub fn len(&self) -> usize {
        self.0.len()
    }

    #[must_use]
    #[inline(always)]
    pub fn remainder(&self) -> &'a [T] {
        self.0.remainder()
    }
}

pub trait LineString3<T: GenericVector3> {
    /// Copy this linestring3 into a linestring2, keeping the axes defined by 'plane'
    /// An axis will always try to keep its position (e.g. y goes to y if possible).
    /// That way the operation is reversible (in regard to axis positions).
    fn copy_to_2d(&self, plane: Plane) -> Vec<T::Vector2>;

    /// Returns `true` if the last and first points of the collection are exactly the same,
    /// indicating a closed loop.
    /// Note that an empty linestring is also "closed" since first() and last() object are the
    /// same. I.e. None
    fn is_connected(&self) -> bool;

    //// Returns an iterator over consecutive pairs of points in the line string, forming continuous lines.
    /// This iterator is created using `.windows(2)` on the underlying point collection.
    #[must_use = "iterator adaptors are lazy and do nothing unless consumed"]
    fn window_iter(&self) -> WindowIterator<'_, T>;

    /// Returns an iterator over pairs of points in the line string, forming disconnected edges.
    /// This iterator is created using `.chunks_exact(2)` on the underlying point collection.
    #[must_use = "iterator adaptors are lazy and do nothing unless consumed"]
    fn chunk_iter(&self) -> ChunkIterator<'_, T>;

    /// Simplify using Ramer–Douglas–Peucker algorithm adapted for 3d
    fn simplify_rdp(&self, distance_predicate: T::Scalar) -> Self;

    /// Simplify using Visvalingam–Whyatt algorithm adapted for 3d.
    /// This algorithm will delete 'points_to_delete' of points from the polyline with the smallest
    /// area defined by one point and it's neighbours.
    fn simplify_vw(&self, points_to_delete: usize) -> Self;

    /// Apply an operation over each coordinate.
    /// Useful when you want to round the values of the coordinates.
    fn apply<F>(&mut self, f: &F)
    where
        F: Fn(T) -> T;

    /// Creates an iterator that generates intermediate points on a line segment
    /// with steps shorter or equal to the specified distance.
    ///
    /// The iterator ensures that the original points from the line segment are always included.
    /// Each step along the line segment is of equal or shorter length than the specified distance.
    ///
    /// # Parameters
    ///
    /// - `distance`: The desired maximum distance between generated points.
    ///
    /// # Returns
    ///
    /// Returns a `DiscreteLineSegmentIterator` that iterates over the generated points.
    fn discretize(&self, distance: T::Scalar) -> DiscreteLineSegmentIterator<'_, T>;
}

impl<T: GenericVector3> LineString3<T> for Vec<T> {
    /// Copy this linestring3 into a linestring2, keeping the axes defined by 'plane'
    /// An axis will always try to keep its position (e.g. y goes to y if possible).
    /// That way the operation is reversible (in regard to axis positions).
    fn copy_to_2d(&self, plane: Plane) -> Vec<T::Vector2> {
        plane.points_to_2d(self.iter()).collect()
    }

    /// Returns `true` if the last and first points of the collection are exactly the same,
    /// indicating a closed loop.
    /// Note that an empty linestring is also "closed" since first() and last() object are the
    /// same. I.e. None
    fn is_connected(&self) -> bool {
        self.is_empty() || self.first().unwrap() == self.last().unwrap()
    }

    //// Returns an iterator over consecutive pairs of points in the line string, forming continuous lines.
    /// This iterator is created using `.windows(2)` on the underlying point collection.
    fn window_iter(&self) -> WindowIterator<'_, T> {
        WindowIterator(self.windows(2))
    }

    /// Returns an iterator over pairs of points in the line string, forming disconnected edges.
    /// This iterator is created using `.chunks_exact(2)` on the underlying point collection.
    fn chunk_iter(&self) -> ChunkIterator<'_, T> {
        ChunkIterator(self.chunks_exact(2))
    }

    /// Simplify using Ramer–Douglas–Peucker algorithm adapted for 3d
    #[inline(always)]
    fn simplify_rdp(&self, distance_predicate: T::Scalar) -> Self {
        simplify::simplify_rdp(self, distance_predicate)
    }

    /// Simplify using Visvalingam–Whyatt algorithm adapted for 3d.
    /// This algorithm will delete 'points_to_delete' of points from the polyline with the smallest
    /// area defined by one point and it's neighbours.
    /// Simplify using Visvalingam–Whyatt algorithm. This algorithm will delete 'points_to_delete'
    /// of points from the polyline with the smallest area defined by one point and it's neighbours.
    #[inline(always)]
    fn simplify_vw(&self, points_to_delete: usize) -> Self {
        simplify::simplify_vw(self, points_to_delete)
    }

    fn apply<F: Fn(T) -> T>(&mut self, f: &F) {
        self.iter_mut().for_each(|v| *v = f(*v));
    }

    /// Creates an iterator that generates intermediate points on a line segment
    /// with steps shorter or equal to the specified distance.
    ///
    /// The iterator ensures that the original points from the line segment are always included.
    /// Each step along the line segment is of equal or shorter length than the specified distance.
    ///
    /// # Examples
    ///
    /// ```
    /// # use vector_traits::glam::Vec3;
    /// use linestring::prelude::LineString3;
    ///
    /// let points = vec![
    ///     Vec3::new(0.5, 0.0, 0.0),
    ///     Vec3::new(1.0, 0.0, 0.0),
    ///     Vec3::new(2.0, 0.0, 0.0),
    ///     Vec3::new(3.0, 0.0, 0.0),
    /// ];
    ///
    /// let distance = 0.8;
    /// let iterator = points.discretize(distance);
    ///
    /// let result: Vec<Vec3> = iterator.collect();
    ///
    /// // Ensure that the result contains the original points and intermediate steps
    /// assert_eq!(result, vec![
    ///     Vec3::new(0.5, 0.0, 0.0),
    ///     Vec3::new(1.0, 0.0, 0.0),
    ///     Vec3::new(1.5, 0.0, 0.0),
    ///     Vec3::new(2.0, 0.0, 0.0),
    ///     Vec3::new(2.5, 0.0, 0.0),
    ///     Vec3::new(3.0, 0.0, 0.0),
    /// ]);
    /// ```
    ///
    /// # Parameters
    ///
    /// - `distance`: The desired maximum distance between generated points.
    ///
    /// # Returns
    ///
    /// Returns a `DiscreteLineSegmentIterator` that iterates over the generated points.
    fn discretize(&self, distance: T::Scalar) -> DiscreteLineSegmentIterator<'_, T> {
        DiscreteLineSegmentIterator::new(self, distance)
    }
}

/// A workaround for Rust's limitations where external traits cannot be implemented for external types.
///
/// The `Approx` trait provides methods for performing approximate equality comparisons on types.
/// It serves as a workaround for Rust's limitations, allowing you to implement approximate
/// equality checks for types not originally designed with this capability.
///
/// This trait leverages the `approx` crate and its traits to perform approximate equality
/// comparisons. The methods in this trait are wrappers around the corresponding methods provided
/// by the `approx` crate.
///
pub trait Approx<T: GenericVector3> {
    /// Checks if two instances are nearly equal within a specified tolerance in ULPs (Units in the Last Place).
    ///
    /// This method delegates to the `approx::UlpsEq::ulps_eq` method, performing approximate equality checks
    /// one time per coordinate axis.
    fn ulps_eq(
        &self,
        other: &Self,
        epsilon: <T::Scalar as AbsDiffEq>::Epsilon,
        max_ulps: u32,
    ) -> bool;

    /// Checks if two instances are nearly equal within a specified absolute difference tolerance.
    ///
    /// This method delegates to the `approx::AbsDiffEq::abs_diff_eq` method, performing approximate equality checks
    /// one time per coordinate axis.
    fn abs_diff_eq(&self, other: &Self, epsilon: <T::Scalar as AbsDiffEq>::Epsilon) -> bool;
}

impl<T: GenericVector3> Approx<T> for Vec<T> {
    fn ulps_eq(
        &self,
        other: &Self,
        epsilon: <T::Scalar as AbsDiffEq>::Epsilon,
        max_ulps: u32,
    ) -> bool {
        self.len() == other.len()
            && self
                .iter()
                .zip(other.iter())
                .all(|(a, b)| a.is_ulps_eq(*b, epsilon, max_ulps))
    }

    fn abs_diff_eq(&self, other: &Self, epsilon: <T::Scalar as AbsDiffEq>::Epsilon) -> bool {
        self.len() == other.len()
            && self
                .iter()
                .zip(other.iter())
                .all(|(a, b)| a.is_abs_diff_eq(*b, epsilon))
    }
}

#[allow(clippy::many_single_char_names)]
#[inline(always)]
fn cross_abs_squared<T: GenericVector3>(a: T, b: T) -> T::Scalar {
    let x = a.y() * b.z() - a.z() * b.y();
    let y = a.z() * b.x() - a.x() * b.z();
    let z = a.x() * b.y() - a.y() * b.x();
    x * x + y * y + z * z
}

#[inline(always)]
/// The distance between the line a->b to the point p is the same as
/// distance = |(a-p)×(a-b)|/|a-b|
/// <https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Another_vector_formulation>
/// Make sure to *not* call this function with a-b==0
/// This function returns the distance²
pub fn distance_to_line_squared<T: GenericVector3>(l0: T, l1: T, p: T) -> T::Scalar {
    let l0_sub_l1 = l0 - l1;
    let l0_sub_p = l0 - p;
    cross_abs_squared(l0_sub_p, l0_sub_l1)
        / (l0_sub_l1.x() * l0_sub_l1.x()
            + l0_sub_l1.y() * l0_sub_l1.y()
            + l0_sub_l1.z() * l0_sub_l1.z())
}

/// Same as `distance_to_line_squared<T>` but it can be called when a-b might be 0.
/// It's a little slower because it does a==b test
#[inline(always)]
pub fn distance_to_line_squared_safe<T: GenericVector3>(l0: T, l1: T, p: T) -> T::Scalar {
    if l0.is_ulps_eq(
        l1,
        T::Scalar::default_epsilon(),
        T::Scalar::default_max_ulps(),
    ) {
        // give the point-to-point answer if the segment is a point
        l0.distance_sq(p)
    } else {
        let l0_sub_l1 = l0 - l1;
        let l0_sub_p = l0 - p;
        cross_abs_squared(l0_sub_p, l0_sub_l1)
            / (l0_sub_l1.x() * l0_sub_l1.x()
                + l0_sub_l1.y() * l0_sub_l1.y()
                + l0_sub_l1.z() * l0_sub_l1.z())
    }
}

/// Rudimentary save line strings to .obj file function
pub fn save_to_obj_file<T: GenericVector3>(
    filename: impl AsRef<Path>,
    object_name: &str,
    lines: &[Vec<Line3<T>>],
) -> Result<(), LinestringError> {
    let mut point_set = collections::HashMap::<String, usize>::new();
    // try to de-duplicate points
    for i in lines.iter() {
        for j in i.iter() {
            let a_str = format!("v {:+e} {:+e} {:+e}", j.start.x(), j.start.y(), j.start.z());
            let len = point_set.len();
            let _ = point_set.entry(a_str).or_insert(len);
            let a_str = format!("v {:+e} {:+e} {:+e}", j.end.x(), j.end.y(), j.end.z());
            let len = point_set.len();
            let _ = point_set.entry(a_str).or_insert(len);
        }
    }
    let mut file = io::BufWriter::new(fs::File::create(filename)?);
    writeln!(file, "o {object_name}")?;
    for (k, v) in point_set.iter().sorted_unstable_by(|a, b| a.1.cmp(b.1)) {
        writeln!(file, "{} #{}", k, v + 1)?;
    }
    for i in lines.iter() {
        for j in i.iter() {
            let str0 = format!("v {:+e} {:+e} {:+e}", j.start.x(), j.start.y(), j.start.z());
            if let Some(p0) = point_set.get(&str0) {
                let str1 = format!("v {:+e} {:+e} {:+e}", j.end.x(), j.end.y(), j.end.z());
                if let Some(p1) = point_set.get(&str1) {
                    writeln!(file, "l {} {}", p0 + 1, p1 + 1)?;
                } else {
                    return Err(LinestringError::InternalError(format!(
                        "HashMap can't find the key:{str1}"
                    )));
                }
            } else {
                return Err(LinestringError::InternalError(format!(
                    "HashMap can't find the key:{str0}"
                )));
            }
        }
    }
    file.flush()?;
    Ok(())
}

pub struct DiscreteLineSegmentIterator<'a, T: GenericVector3> {
    points: &'a Vec<T>,
    distance: T::Scalar,
    current_index: usize,
    substep_index: usize,
    num_substeps: usize,
    step: T,
}

impl<'a, T: GenericVector3> DiscreteLineSegmentIterator<'a, T> {
    fn new(points: &'a Vec<T>, distance: T::Scalar) -> Self {
        DiscreteLineSegmentIterator {
            points,
            distance,
            current_index: 0,
            substep_index: 0,
            num_substeps: 0,
            step: T::new_3d(T::Scalar::ZERO, T::Scalar::ZERO, T::Scalar::ZERO),
        }
    }

    fn initialize_step(&mut self) {
        let start_point = self.points[self.current_index];
        let end_point = self.points[self.current_index + 1];
        let direction = end_point - start_point;
        let segment_length = direction.magnitude();
        if segment_length > T::Scalar::ZERO {
            let num_substeps = (segment_length / self.distance).ceil();
            self.num_substeps = num_substeps.as_();
            self.step = direction / num_substeps;
        } else {
            // Handle the case when segment_length is exactly 0.0
            self.num_substeps = 1;
            self.step = T::new_3d(T::Scalar::ZERO, T::Scalar::ZERO, T::Scalar::ZERO);
        }
    }
}

impl<T: GenericVector3> Iterator for DiscreteLineSegmentIterator<'_, T>
where
    usize: AsPrimitive<<T as HasXY>::Scalar>,
{
    type Item = T;

    fn next(&mut self) -> Option<Self::Item> {
        if self.points.is_empty() {
            return None;
        }
        match self.points.len() - 1 {
            n if self.current_index < n => {
                if self.substep_index == 0 {
                    // Initialize step on a new line segment
                    self.initialize_step();
                }

                if self.substep_index < self.num_substeps {
                    let start_point = self.points[self.current_index];
                    let result = start_point + self.step * (self.substep_index.as_());
                    self.substep_index += 1;
                    return Some(result);
                }

                // Move to the next line segment
                self.current_index += 1;
                self.substep_index = 0;
                self.next()
            }
            n if self.current_index == n => {
                // return the very last point of the list
                self.current_index += 1;
                self.substep_index = 0;
                self.points.last().copied()
            }
            _ => None,
        }
    }
}