一个基于坐标系的几何画图网页,对于计算几何,这个网页用来画图、打草稿啥的挺好用。

向量

点积(内积)

线性代数中的通用公式

对于n维的向量a_ib_i表示向量在第i个维度上的分量。

亦可化简为:

对于二维向量

  • 设向量的模分别为,其夹角为,则有

  • 向量自身的点积为其长度的平方。
  • 点积与方向的关系:
的方向
<0 钝角 反向
=0 直角 垂直
>0 锐角 同向

叉积(外积)

线性代数中的通用公式

涉及到行列式展开计算,这里用三维举例。

对于二维向量

  • 叉积可以用来求出法向量方向。
  • 叉积为向量围成的平行四边形的面积(带符号)。
  • ,则表示可通过逆时针方向旋转至方向,小于0则为顺时针,等于0则表示两个向量共线。

单位向量

对于向量,其单位向量

投影

对于向量上的投影为

例题

POJ 2007:Scrambled Polygon

代码模板

// 二维向量
struct Vector
{
    double x, y;

    Vector ()
    {}

    Vector (double x, double y) : x(x), y(y)
    {}

    // 点积
    double operator* (Vector v)
    {
        return x * v.x + y * v.y;
    }

    // 叉积
    double operator^ (Vector v)
    {
        return x * v.y - y * v.x;
    }

    // 乘法 
    Vector operator* (double d)
    {
        return Vector(x * d, y * d);
    }

    // 加法
    Vector operator+ (Vector v)
    {
        return Vector(x + v.x, y + v.y);
    }

    // 减法
    Vector operator- (Vector v)
    {
        return Vector(x - v.x, y - v.y);
    }

    // 模
    double length ()
    {
        return sqrt((*this) * (*this));
    }

    // 单位化
    Vector unit ()
    {
        return (*this) * (1.0 / length());
    }

};

点与向量的表示方法相同,都可以通过两个坐标来表示,直接重用向量即可,无需重复定义。

绕点旋转

点绕点旋转,其实可以直接理解成旋转到,其中点A为旋转的中心。

可以先将A平移至原点,旋转完成后,再平移回原来位置。

设:

  • 原向量与x轴夹角为,表示为(x,y)

  • 旋转后的向量与x轴夹角为,表示为

  • r

有如下推导:

代码模板

#define Point Vector

    // 作为Vector成员函数
    Vector rotate (double alpha)
    {
        return Vector(x + (x * cos(alpha) - y * sin(alpha)),
                      y + (x * sin(alpha) + y * cos(alpha)));
    }

例题

POJ 2624:4th Point

直线

两点连成线,用两个点相减即可表示。

点-线距离

对于点A、B所在直线和点P,P点到直线的距离有:

P的距离为

点-线位置

对于点A、B所在直线和点P,利用前面介绍的叉积与向量方向关系的性质,构造向量,通过来判断P和的位置关系:

过点做垂线

对于点A、B所在直线和点P:

垂线:将直线所在向量旋转,并平移至经过点P即可。

垂足:利用点积求出上投影的长度即为垂足的位置。

直线-直线的位置

对于点AB所在直线和C、D所在直线

如果二者不相交(平行或相等),则AB同一侧,其叉积值和符号都相同,相减为0。反之则不为零。

  • 如果二者不相交,则再通过A的距离来判断二者是否相等。距离为0则相等,否则平行。

  • 若二者相交,通过相似三角形的比例关系,可求得则其交点位置为:

代码

// 二维直线
struct Line
{
    Point x, y;

    // 通过两点构造直线
    Line (Point x, Point y) : x(x), y(y)
    {}

    // 点和方向向量来构造线
    static Line makeLine (Point x, Vector v)
    {
        return Line(x, x + v);
    }

    bool operator== (const Line &rhs) const
    {
        return x == rhs.x &&
               y == rhs.y;
    }

    bool operator!= (const Line &rhs) const
    {
        return !(rhs == *this);
    }

    // 线长度
    double length ()
    {
        return (y - x).length();
    }

    // 点到该直线的距离
    double dist (Point p)
    {
        return fabs((x - p) ^ (y - p)) / length();
    }

    // #define EPS 1e-6
    // 判断点和直线的位置
    int side (Point p)
    {
        double result = (y - x) ^(p - x);
        if (fabs(result) < EPS)
            return 0;   // 在线上
        else if (result > 0)
            return 1;   // 左侧
        else
            return -1;  // 右侧
    }

    // 过点做垂线
    Line vertical (Point p)
    {
        return makeLine(p, (y - x).rotate(PI / 2));
    }

    // 垂足
    Point foot (Point p)
    {
        Vector self = y - x;
        return x + self.unit() * self.project(p - x);
    }

    Point intersect (Line l)
    {
        double s1 = ((x-l.x) ^(l.y - l.x)) / 2 ;
        double s2 = ((l.y-l.x) ^(y - l.x)) / 2;
        if (fabs(s1 + s2) < EPS)
        {
            if (l.dist(x))
                return l.x;   // 重合
            else
                return l.y;   // 平行
        } else
            return x + (y - x) * (s1 / (s1 + s2)); // 交点
    }

};

线段

点-线段位置

对于点P和线段AB,构造向量,求叉积有:

对于共线的情况,再根据两个端点的坐标判断,P是否在线段AB上,即满足以下情况的,即为点P在线段AB上:

直线-线段位置

对于直线AB和线段CD,二者位置关系的判断:

  1. 构造向量,并求叉积
  2. 若叉积都不为0,且符号相同,则可认为点都在向量同一侧,或二者重合。
    • 如果一个端点到直线的距离为0,则二者重合。
    • 否则,二者不相交。
  3. 否则二者相交。交点位置求法同直线-直线交点位置。

线段-线段位置

规范相交:两条线段恰有一个不是端点的公共点。

对于直线AB和线段CD,利用叉积的符号判断是否为规范相交,即任一条线段的两个端点都分属另一条线段的两侧,满足以下条件:

如不满足以上条件,再判断两条线段的端点是否在另一条线段上,若是则为非规范相交,否则二者不相交。

例题

POJ 3304:Segments

POJ 1039:Pipe

代码

// 二维线段
struct Seg
{
    Point x, y;

    Seg (const Point &x, const Point &y) : x(x), y(y)
    {}

    // 判断点是否在线段内
    bool pointIn (Point p)
    {
        double cross = (x - p) ^(y - x);
        if (fabs(cross) < EPS)
            return 0;
        if (fabs(p.x - min(x.x, y.x)) < EPS &&
            fabs(p.y - min(x.y, y.y)) < EPS &&
            fabs(max(x.x, y.x) - p.x) < EPS &&
            fabs(max(x.y, y.y) - p.y) < EPS)
            return 1;
        return 0;
    }

    Point intersect (Line l)
    {
        int s1 = l.side(x);
        int s2 = l.side(y);
        if (s1 * s2 > 0)
        {
//            if (l.dist(x) < EPS)
//                return l.y; // 重合
//            else
            return Point(EPS / 10, EPS / 10);// 不相交
        } else if (s1 == 0 && s2 == 0)
        {
            return l.x;// 重合
        } else if (s1 == 0)   // 一个端点在直线上
        {
            return x;
        } else if (s2 == 0)  // 另一个端点在直线上
        {
            return y;
        } else
        {
            return Line(x, y).intersect(l); // 交点
        }
    }

    // 判断线段相交
    bool isIntersected (Seg s)
    {
        Point p1 = x;
        Point p2 = y;
        Point p3 = s.x;
        Point p4 = s.y;
        double s1 = (p3 - p1) ^(p4 - p1);
        double s2 = (p4 - p2) ^(p3 - p2);
        if (((p2 - p1) ^ (p3 - p1) * ((p2 - p1) ^ (p4 - p1))) < -EPS &&
            (((p4 - p3) ^ (p1 - p3)) * ((p4 - p3) ^ (p2 - p3))) < -EPS)
            return 1;   // 规范相交
        if (s.pointIn(p1) || s.pointIn(p2) || pointIn(p3) || pointIn(p4))
            return 1;   // 非规范相交

        return 0;
    }

    bool operator== (const Seg &rhs) const
    {
        return x == rhs.x &&
               y == rhs.y;
    }

    bool operator!= (const Seg &rhs) const
    {
        return !(rhs == *this);
    }

};

多边形

凸多边形

百度百科:凸多边形是一个内部为凸集的简单多边形。凸多边形(Convex Polygon)指如果把一个多边形的所有边中,任意一条边向两方无限延长成为一直线时,其他各边都在此直线的同旁,那么这个多边形就叫做凸多边形,其内角应该全不是优角,任意两个顶点间的线段位于多边形的内部或边上。

就是没有凹角的多边形,特点是所有边所在的直线都不会穿过其他边。

多边形面积

对多边形进行三角剖分:

  1. 设原点为O
  2. 按逆时针方向为个条边指定方向
  3. 对于每条边,累加

点-多边形位置

设多边形P的顶点序列为、待判定的点 q

射线法

算法思路

q做水平射线l

  • lP的边界不相交,则P的外部。

  • 若相交,则再分析交点数的奇偶性。

    交点数 点相对于多边形的位置
    奇数 内部
    偶数 外部
  • 需要注意的是,当多边形不是凸多边形时,需要考虑特殊情况:

    1. 射线恰好经过P的某个顶点
    2. 射线恰好与P的某条边重合
算法特点
  1. 运算速度快
  2. 精度高
  3. 特殊情况较多

转角法

算法思路

沿多边形走一圈,累计绕点旋转了多少角度。(需要保证边的方向一致)

角度 点相对于多边形的位置
0 外部
点在多边形的边上
内部

假设向量 其角度通过点积公式可以求出:

(仅适用于凸包)或使用累乘叉乘,判断点和边的位置关系:

  • 当点位于所有边向量的同一侧时(叉乘值同号),其在多边形内部。
  • 如果没有出现在同一侧,则为外部(叉乘值异号)。
  • 在多边形上的情况是可以直接判断出来的(叉乘为0)。
算法特点
  • 几乎没有特殊情况
  • 精度低、速度较慢(需要用到反三角函数、开方等)

求凸包(凸壳)

给定一个平面点集,要求找到一个最小的凸多边形,满足点集中所有点都在凸多边形内部或边上。

稳定凸包

对于已有的凸包中,无法通过在平面中添加一个点获取到更大的凸包,则该已有的凸包称为稳定凸包

特性:凸包的每条边上,都有大于等于3个的顶点。

具体可见例题【POJ 1228 Grandpa's Estate】。

半平面交

半平面

一条直线将一个平面分为2个半平面,直线是有向的,设直线方向左边的平面为我们研究的半平面(包含直线)。

半平面交

被(多个)半平面包含的点的集合。

性质:半平面交的结果是一个凸区域。

求交

对于每一条多边形的每一条边、表示半平面的直线

  1. 对多边形按照逆时针排序。

  2. 如果边的起点A位于半平面l中,则将点A加入结果集。

  3. 求边l的交点,并将结果加入结果集。(即使点A不在半平面中,也需要这步操作)

增量法

  1. 构造一个足够大的矩形。
  2. 依次用半平面和该矩形求交。

时间复杂度:

分治法(二分)

  1. n个半平面分成2个大小近似相等的集合。
  2. 递归地构造凸多边形区域A与。
  3. 合并A与。

时间复杂度:

例题

POJ 1569:Myacm Triangles
POJ 1113-Wall
POJ 1228-Grandpa's Estate

代码

// 二维多边形
struct Polygon
{
    // 点集
    vector<Point> points;

    Polygon (const vector<Point> &points) : points(points)
    {}

    Polygon ()
    {}

    /**
     * 根据已排序的点集来求半个凸壳(上半个或下半个)
     * @param pts 已排序的点集
     * @return 半个凸包中的点,有序,逆时针序
     */
    deque<Point> getHalfConvexHull (const vector<Point> &pts) const
    {
        deque<Point> deq;
        int si = pts.size();
        int i = 2;

        deq.pb(pts[0]);
        deq.pb(pts[1]);
        while (i < si)
        {
            Point nt = pts[i];
            ++i;
            while (deq.size() >= 2)
            {
                Point p2 = deq[deq.size() - 2], p1 = deq.back();
                Vector b = p2 - p1, a = p1 - nt;
                double d = (b ^ a);
                if (d <= 0)
                    deq.pop_back();
                else
                    break;
            }

            deq.pb(nt);
        }
        return deq;
    }

    /**
     * 水平序扫描法
     * 求该多边形的凸包
     * @return 凸包的点集,有序,逆时针序
     */
    vector<Point> getConvexHull ()
    {
        vector<Point> pts(points);
        vector<Point> con;


        // 2个点无法构成凸包
        if (pts.size() <= 2)
        {
            return con;
        }

        // 下凸壳
        sort(pts.begin(), pts.end());
        deque<Point> temp = getHalfConvexHull(pts);
        for (int i = 0; i < temp.size(); ++i)
        {
            con.pb(temp[i]);
        }
        temp.clear();


        // 上凸壳
        sort(pts.begin(), pts.end(), greater<Point>());
        temp = getHalfConvexHull(pts);
        for (int i = 1; i < temp.size() - 1; ++i)
        {
            con.pb(temp[i]);
        }

        return con;
    }
};