阅读 123

Rust实现计算几何

本文介绍如何使用Rust程序设计语言实现常见的计算几何算法,本文假设所有的问题都是基于二维世界的。

数据结构

点/向量

使用坐标可以表示点或者向量。

#[derive(Copy,Clone,PartialEq)]
struct Point {
    x: f64,
    y: f64
}

type Vector = Point;
复制代码

线段/直线

使用两个不相同的点可以表示线段或者直线。

#[derive(Copy,Clone)]
struct Segment {
    start: Point,
    end: Point
}

type Line = Segment;
复制代码

圆可以使用一个点和半径来表示。

#[derive(Copy,Clone)]
struct Circle {
    center: Point,
    radius: f64
}
复制代码

多边形

多边形可以使用一组点来表示。

#[derive(Clone)]
struct Polygen {
    points: Vec<Point>
}
复制代码

算法

向量的基本运算

在数据结构的基础上,可以定义出向量加、向量减、长度、标量乘等基本运算。

impl Point {

    // 向量加
    fn add(&self, other: &Self) -> Point {
        Point {
            x: self.x + other.x,
            y: self.y + other.y
        }
    }

    // 向量减
    fn sub(&self, other: &Self) -> Point {
        Point {
            x: self.x - other.x,
            y: self.y - other.y
        }
    }

    // 向量长度
    fn abs(&self) -> f64 {
        (self.x * self.x + self.y * self.y).sqrt()
    }

    // 标量乘
    fn scale(&self, x: f64) -> Point {
        Point { x: self.x * x, y: self.y * x }
    }
}
复制代码

向量内积

设向量ab的夹角为\theta0\le\theta\le \pi),那么ab 的内积为a\cdot b=|a||b|\cos\theta

根据余弦定理

\cos\theta=\frac{|a|^2+|b|^2-|b-a|^2}{2|a||b|}

因此

a\cdot b= \frac{|a|^2+|b|^2-|b-a|^2}{2}=x_ax_b + y_ay_b
impl Point {
    fn dot(&self, other: &Self) -> f64 {
        self.x * other.x + self.y * other.y
    }
}
复制代码

向量外积

设向量ab的夹角为\theta0\le\theta\le \pi),那么ab 的外积为a\cdot b=|a||b|\sin\theta

根据三维空间的外积公式

|a\times b|=(y_az_b-z_ay_b,x_az_b-x_az_b,x_ay_b-y_ax_b)

外积大小为

|a\times b|=|a||b|\sin\theta=x_ay_b-y_ax_b
impl Point {
    fn cross(&self, other: &Self) -> f64 {
        self.x * other.y - self.y * other.x
    }
}
复制代码

判断直线垂直

如果两条直线垂直,那么他们对应的向量内积为零。

impl Segment {
    // 返回线段对应的向量
    fn vector(&self) -> Vector {
        self.end.sub(&self.start)
    }
    
    fn is_orthogonal(&self, other: &Self) -> bool {
        self.vector().dot(&other.vector()).abs() < std::f64::EPSILON
    }    
}
复制代码

判断直线水平

如果两条直线水平,那么他们对应的向量外积为零。

impl Segment {
    fn is_parallel(&self, other: &Self) -> bool {
        self.vector().cross(&other.vector()).abs() < std::f64::EPSILON
    }
}
复制代码

计算投影

计算点p在直线\overrightarrow {p_1p_2}上的投影的方法如下:

  1. 计算出p和p1构成的向量hypo=p-p_1和直线对应的向量base=p_2-p_1
  2. 利用余弦可以计算出投影点x到p_1的长度t=|hypo|\cos\theta,其中\cos\theta=\frac{hypo\cdot base}{|hypo||base|}
  3. 通过和线段p_1p_2方向相同的单位向量e,即可计算出x=p_1+te
impl Point {
    // 计算向量之间余弦
    fn cosine(&self, other: &Self) -> f64 {
        self.dot(other) / self.abs() / other.abs()
    }
    
    // 返回和向量在同一方向的单位向量
    fn unit(&self) -> Point {
        let abs = self.abs();
        Point { x: self.x / abs, y: self.y / abs }
    }
}
 
impl Segment {   
    // 计算点p在直线上的投影
    fn projection(&self, p: &Point) -> Point {
        if *p == self.start {
            *p
        } else {
            let hypo = p.sub(&self.start);
            let base = self.vector();
            let t = base.cosine(&hypo) * hypo.abs();
            self.start.add(&self.vector().unit().scale(t))
        }
    }
}
复制代码

计算映像

通过投影点p',可以很方便计算出映像点x=p+2(p'-p)

impl Segment {
    // 点p在直线上的映像
    fn reflection(&self, p: &Point) -> Point {
        p.add(&self.projection(&p).sub(&p).scale(2.0))
    }
}
复制代码

计算距离

点到点的距离

点到点的距离就是两点构成的线段长度。

impl Point {
    fn distance_to_point(&self, other: &Self) -> f64 {
        self.sub(other).abs()
    }
}
复制代码

点到直线的距离

向量外积就是两个向量为边构成的平行四边形面积,点p'到直线\overrightarrow{p_0p_1}的距离也就是\overrightarrow{p_0p'}\overrightarrow{p_0p_1}构成的平行四边形的高,因此距离就是\frac{\overrightarrow{p_0p_1}\times\overrightarrow{p_0p'}}{|\overrightarrow{p_0p_1}|}

impl Point {
    fn distance_to_line(&self, line: &Line) -> f64 {
        self.sub(&line.start).cross(&line.vector()).abs() / line.vector().abs()
    }
}
复制代码

点到线段的距离

  1. 向量\overrightarrow{p_2p_1}与向量\overrightarrow{pp_1}的夹角\theta大于90 度(或小于-90度时,距离为点p到点p_1的距离。
  2. 向量\overrightarrow{p_1p_2}与向量\overrightarrow{pp_2} 的夹角\theta大于90度(或小于-90度)时,距离为点p到点p_2的距离。
  3. 除上述两种情况外,直接求点p到直线p_1p_2的距离。
impl Point {
    fn distance_to_segment(&self, s: &Segment) -> f64 {
        if self.sub(&s.start).dot(&s.end.sub(&s.start)) < std::f64::EPSILON {
            self.distance_to_point(&s.start)
        } else if self.sub(&s.end).dot(&s.start.sub(&s.end)) < std::f64::EPSILON {
            self.distance_to_point(&s.end)
        } else {
            self.distance_to_line(s)
        }
    }
}
复制代码

线段到线段的距离

如果两条线段相交,则它们的距离为0,否则线段s_1与线段s_2 的距离为以下四个距离中最小的一个:

  1. 线段s_1与线段s_2的端点s_2.p_1的距离。
  2. 线段s_1与线段s_2的端点s_2.p_2的距离。
  3. 线段s_2与线段s_1的端点s_1.p_1的距离。
  4. 线段s_2与线段s_1的端点s_1.p_2的距离。
impl Segment {
    fn distance(&self, other: &Self) -> f64 {
        if self.intersect(other) {
            0.0
        } else {
            self.start.distance_to_segment(other)
                .min(self.end.distance_to_segment(other))
                .min(other.start.distance_to_segment(self))
                .min(other.end.distance_to_segment(self))
        }
    }
}
复制代码

计算向量相对方向

  1. 外积的大小\overrightarrow{p_0p_1}\times \overrightarrow{p_0p_2}为正时,可确定\overrightarrow{p_0p_2}\overrightarrow{p_0p_1}的逆时针位置。
  2. 外积的大小\overrightarrow{p_0p_1}\times \overrightarrow{p_0p_2}为负时, 可确定\overrightarrow{p_0p_2}\overrightarrow{p_0p_1}的顺时针位置。
  3. 上述情况都不符合时,表示p_2位于直线\overrightarrow{p_0p_1}上。\cos\theta\theta大于90度或小于-90度时为负,因此内积心\overrightarrow{p_0p_1}\cdot\overrightarrow{p_0p_2}为负时,可确定p_2位于线 段\overrightarrow{p_0p_1}后方,即顺序为p_2\rightarrow p_0 \rightarrow p_1
  4. 如果|\overrightarrow{p_0p_2}|>|\overrightarrow{p_0p_1}|,可以确定p_0\rightarrow p_1 \rightarrow p_2
  5. 不符合4时,可确定p_2位于线段\overrightarrow{p_0p_1}上。
#[derive(PartialEq)]
enum ClockDirection {
    OnlineBack = 2,
    CounterClockwise = 1,
    Clockwise = -1,
    OnlineFront = -2,
    OnSegment = 0
}

impl Point {
    fn ccw(&self, p1: &Self, p2: &Self) -> ClockDirection {
        let v1 = p1.sub(self);
        let v2 = p2.sub(self);
        let c = v1.cross(&v2);
        let d = v1.dot(&v2);
        if c > std::f64::EPSILON {
            ClockDirection::CounterClockwise
        } else if c < -std::f64::EPSILON {
            ClockDirection::Clockwise
        } else if d < -std::f64::EPSILON {
            ClockDirection::OnlineBack
        } else if v1.abs() < v2.abs() {
            ClockDirection::OnlineFront
        } else {
            ClockDirection::OnSegment
        }
    }
}
复制代码

判断线段相交

分别检查两条线段,如果双方都符合“另一条线段的两个端点分别位于当前线段的顺时针方向和逆时针方向”,则两条线段相交。

impl Segment {
    fn intersect(&self, other: &Self) -> bool {
        (self.start.ccw(&self.end, &other.start) as i32) * (self.start.ccw(&self.end, &other.end) as i32) <= 0 &&
        (other.start.ccw(&other.end, &self.start) as i32) * (other.start.ccw(&other.end, &self.end) as i32) <= 0
    }
}
复制代码

线段的交点

利用点到直线距离的方法,我们可以得到线段s_1的两个断点到线段s_2的距离d_1d_2,通过它们的比例可以计算出交点在线段s_1的相对位置,从而得到端点位置。

impl Segment {
    fn cross_point(&self, s: &Segment) -> Point {
        let base = self.vector();
        let d1 = base.cross(&s.start.sub(&self.start)).abs() / base.abs();
        let d2 = base.cross(&s.end.sub(&self.start)).abs() / base.abs();
        let t = d1 / (d1 + d2);
        s.start.add(&s.vector().scale(t))
    }
}
复制代码

圆是直线的交点

求圆是直线的交点过程如下:

  1. 求出圆心到直线的投影点pr
  2. 利用勾股定理,可以求出投影点到两个交点的距离base
  3. 利用直线相同方向的单位向量e,可以求出两个交点分别为pr-base\cdot epr+base\cdot e
impl Circle {
    fn cross_points_with_line(&self, line: &Segment) -> (Point, Point) {
        let pr = line.projection(&self.center);
        let e = line.vector().unit();
        let pr_center_length = pr.sub(&self.center).abs();
        let base = (self.radius * self.radius - pr_center_length * pr_center_length).sqrt();
        (pr.add(&e.scale(-base)), pr.add(&e.scale(base)))
    }
}
复制代码

圆与圆的交点

求两个圆交点的方法如下:

  1. 求出两个圆的圆心距d,这个圆心距就是c1圆心到c2圆心的向量;
  2. 由两圆圆心以及其中一个交点所组成的三角形的三条边分别为cl.r、c2.r和d。 根据余弦定理可以求岀向量c2.c-cl.c与cl.c到某交点的向量的夹角a。然后再求出 c2.c-cl.c与x轴的夹角t。
  3. 交点就是以圆心cl.c为起点,大小为cl.r,角度为t+a和t-a的两个向量。
impl Point {
    // 计算向量角度
    fn atan2(&self) -> f64 {
        self.y.atan2(self.x)
    }
    
    // 极坐标生成向量
    fn polar(a: f64, r: f64) -> Point {
        Point { x: a.cos() * r, y: a.sin() * r }
    }
}

impl Circle {
    fn cross_points_with_circle(&self, circle: &Self) -> (Point, Point) {
        let d = self.center.sub(&circle.center).abs();
        let a = ((self.radius * self.radius + d * d - circle.radius * circle.radius) / (2.0 * self.radius * d)).acos();
        let t = circle.center.sub(&self.center).atan2();
        (self.center.add(&Point::polar(t + a, self.radius)), self.center.add(&Point::polar(t - a, self.radius)))
    }
}
复制代码

点的内包

对于由点g_1,\cdots,g_n构成的多边形,我们需要判断点p是否包含在多边形中。 对于构成多边形各边的线段,设g_i-p与g_{i+1}-p分别为向量a和向量b。

  • 如果a和b外积大小为0且内积小于等于0,则点p位于g_ig_{i+1}上。
  • 至于射线与线段g_ig_{i+1}是否相交,可以通过a和b构成的平行四边 形的面积正负,即a和b外积的大小来判断。首先我们调整向量a和b,将y值较小的向量定为a。在这个状态下,如果a和b的外积大小为正(a到b为逆时针)且a和b位于射线两侧, 则可确定射线与边交叉。相交次数为奇数时表示“ 内包”,为偶数时表示“ 不内包”。
#[derive(PartialEq)]
enum Containment {
    Out = 0,
    On = 1,
    In = 2
}

impl Polygen {
    fn from(points: Vec<Point>) -> Polygen {
        Polygen { points: points }
    }

    fn contains(&self, point: &Point) -> Containment {
        let mut cnt = 0;
        for i in 0..self.points.len() {
            let mut a = self.points[i].sub(point);
            let mut b = self.points[(i+1)%self.points.len()].sub(point);
            if a.cross(&b).abs() < std::f64::EPSILON && a.dot(&b) < std::f64::EPSILON {
                return Containment::On;
            }
            if a.y > b.y {
                let temp = a;
                a = b;
                b = temp;
            }
            if a.y < std::f64::EPSILON && b.y > std::f64::EPSILON && a.cross(&b) > std::f64::EPSILON {
                cnt += 1;
            }
        }
        if cnt % 2 == 0 {
            Containment::Out
        } else {
            Containment::In
        }
    }
}
复制代码

凸包

凸包可以使用安德鲁算法求出:

  1. 将给定集合的点按x坐标升序排列,x相同的按y坐标升序排列。
  2. 按下述流程创建凸包的上部:将排序后的点按照x坐标从小到大的顺序加入凸包,如果新加人的点使得不再是凸多边形,则逆序删除之前加入的点,直至重新成为凸多边形为止。
  3. 按下述流程创建凸包的下部:将排序后的点按照x坐标从大到小的顺序加入凸包。 如果新加入的点使得凸包不再是凸多边形,则逆序删除之前加人的点,直至重新成为凸多边形为止。

是否构成凸多边形可以使用向量相对方向判断。

impl Polygen {
    fn convex_hull(_points: &Vec<Point>) -> Polygen {
        let mut points = _points.clone();
        points.sort_by(|a, b| a.x.partial_cmp(&b.x).unwrap().then(a.y.partial_cmp(&b.y).unwrap()));
        // 计算上边缘
        let mut hull: Vec<Point> = Vec::new();
        for &point in points.iter() {
            while hull.len() > 1 {
                let p0 = hull[hull.len()-1];
                let p1 = hull[hull.len()-2];
                let direction = p0.ccw(&p1, &point);
                if direction == ClockDirection::CounterClockwise || direction == ClockDirection::OnlineBack {
                    break;
                }
                hull.pop();
            }
            hull.push(point);
        }
        // 计算下边缘
        points.sort_by(|a, b| a.x.partial_cmp(&b.x).unwrap().reverse().then(a.y.partial_cmp(&b.y).unwrap().reverse()));
        let n_up_bound = hull.len();
        for &point in points.iter() {
            while hull.len() > n_up_bound {
                let p0 = hull[hull.len()-1];
                let p1 = hull[hull.len()-2];
                let direction = p0.ccw(&p1, &point);
                if direction == ClockDirection::CounterClockwise || direction == ClockDirection::OnlineBack {
                    break;
                }
                hull.pop();
            }
            hull.push(point);
        }
        hull.reverse();
        hull.pop();
        Polygen { points: hull }
    }
}
复制代码

完整实现

完整的实现和相应的测试可见GitHub仓库

关注下面的标签,发现更多相似文章
评论