计算几何之基础总结（未完）

略颓废啊，用了两个半星期才搞完了计算几何的基础，而且还没做多少题，觉得晕呼呼的，还是在进入数据结构之前先总结一下吧，不然就忘掉了。

0 - Preface

题目列表大概就是 http://blog.renren.com/share/347994200/7898523786 ，做的很少。

1 - Vector

叉积

a × b = x1y2 - x2y1

叉积主要是用来判断向量的位置关系，如果 **a** × **b** > 0 说明 b 在 a 的顺时针方向。

1
2
3
4
5
6

inline double CrossProduct(const Point& A, const Point& B, const Point& C, const Point& D)
{
  const double x1 = B.x-A.x, y1 = B.y-A.y;
  const double x2 = D.x-C.x, y2 = D.y-C.y;
  return x1*y2 - x2*y1;
}

点积

a ・ b = x1x2 + y1y2

点积觉得除了数学书上写的一般也就没什么特殊用处了，无非就是判断共线，不然就是求角（求余弦）。

直线和线段相交

只要做半次跨立试验就可以了，线段P1P2 和 直线AB 相交，只需要考察 (**AP1** × **AB**) * (**AP2** × **AB**) < 0 即可，也就是 直线AB 跨过 线段P1P2

线段交点 ( Method 1 )

这个问题比较棘手，如果是解析几何的话非常愉悦，但是如果用向量方法来解决就有点不舒服了。我一开始是在 http://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect 这里看到了一个方法，然后打了POJ1269，这个方法很好理解，运用叉积的分配率，但是缺点是Point类必须使用浮点，而且有向量相加。

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

// POJ1269
#include <cmath>
#include <cstdio>

struct Point
{
  double x, y;
  Point(): x(.0), y(.0) { }
  Point(const double a, const double b): x(a), y(b) { }
} p1, p2, p3, p4, p;
inline Point operator+(const Point& a, const Point& b) { return Point(a.x+b.x, a.y+b.y); }
inline Point operator-(const Point& a, const Point& b) { return Point(a.x-b.x, a.y-b.y); }
inline Point operator*(const Point& a, const double b) { return Point(a.x*b, a.y*b); }
inline double operator^(const Point& a, const Point& b) { return a.x*b.y - a.y*b.x; }

int Intersection(const Point& p1, const Point& p2, const Point& p3, const Point& p4, Point& p)
{
  const Point r = p2-p1, s = p4-p3;
  const double cross = r^s;
  if (std::abs(cross) < 1e-7)
    return std::abs((p2-p1)^(p3-p1)) < 1e-7 ? 2 : 0;
  const double t = ((p3-p1)^s) / cross;
  p = p1 + r*t;
  return 1;
}
int T;
int main()
{
  printf("INTERSECTING LINES OUTPUTn");
  for (scanf("%d", &T); T--; )
  {
    scanf("%lf%lf%lf%lf%lf%lf%lf%lf", &p1.x, &p1.y, &p2.x, &p2.y, &p3.x, &p3.y, &p4.x, &p4.y);
    const int ret = Intersection(p1, p2, p3, p4, p);
    if (ret == 1)
      printf("POINT %.2f %.2fn", p.x, p.y);
    else
      printf(ret == 0 ? "NONEn" : "LINEn");
  }
  printf("END OF OUTPUT");
}

点在线段上

这个没什么好说的了，就是点在线段上，看起来很简单，写起来略长。

1
2
3
4
5

inline bool OnSegment(const Point& p, const Point& a, const Point& b)
{
  return std::min(a.x, b.x) <= p.x && p.x <= std::max(a.x, b.x)
      && std::min(a.y, b.y) <= p.y && p.y <= std::max(a.y, b.y);
}

* 线段和线段相交 ( Method 1 )

对与规范相交，需要做两次跨立试验。需要注意的是，如果两条线平行但是不相交，跨立试验会出错，需要特判。网上很多说要快速排斥试验的，觉得没什么必要，或者说我没有遇到需要快速排斥试验的题目。

最棘手的是非规范相交，有两种情况：

  • 点在线段上
  • 线段和线段部分重叠

点P 在 线段AB 上的判断方法就是： **AP** × **AB** = 0 并且 **PA** ・ **PB** < 0 。第一句的意思是面积为零，也就是 点P 在 直线AB 上，第二句表示 <**PA**, **PB**> > pi/2，又因为点在直线上了，因此夹角只能是 pi ，也就是两条向量反向，即点在线段上。

线段和线段部分重叠的判断方法是：首先有个面积=0，接下来就是某个点在线段上了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

int Intersection(const Point& p1, const Point& p2, const Point& p3, const Point& p4)
{
  const int c1 = CrossProduct(p1, p3, p1, p2);
  const int c2 = CrossProduct(p1, p4, p1, p2);
  const int c3 = CrossProduct(p3, p1, p3, p4);
  const int c4 = CrossProduct(p3, p2, p3, p4);
  if (c1*c2 < 0 && c3*c4 < 0) return 5;                     // at a point
  if (c1 == 0 && DotProduct(p3, p1, p3, p2) < 0) return 3;  // at p3
  if (c2 == 0 && DotProduct(p4, p1, p4, p2) < 0) return 4;  // at p4
  if (c3 == 0 && DotProduct(p1, p3, p1, p4) < 0) return 1;  // at p1
  if (c4 == 0 && DotProduct(p2, p3, p2, p4) < 0) return 2;  // at p2
  if (c1*c2*c3*c4 == 0 && (
      OnSegment(p3, p1, p2) || OnSegment(p4, p1, p2) ||
      OnSegment(p1, p3, p4) || OnSegment(p2, p3, p4)))
    return 6;                                               // on segment
  return 0;                                                 // none
}

! 线段和线段相交 ( Method 2 )

上面那个方法和明显就能发现是我自己意淫的，虽然说是对的，但是求不了交点，而且也很麻烦，正解应该是这样的：

需要用到定比分点的知识，参见：http://baike.baidu.com/view/900146.htm

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

int Intersection(const Point& p1, const Point& p2, const Point& p3, const Point& p4, Point& p)
{
  const double s1 = CrossProduct(p1, p3, p1, p4);
  const double s2 = CrossProduct(p2, p3, p2, p4);
  const double s3 = CrossProduct(p3, p1, p3, p2);
  const double s4 = CrossProduct(p4, p1, p4, p2);
  if (s1*s2 < 0 && s3*s4 < 0)
  {
    p.x = (s1*p2.x - s2*p1.x) / (s1 - s2);
    p.y = (s1*p2.y - s2*p1.y) / (s1 - s2);
    return 5;
  }
  if (FloatEqual(s1, 0) && FloatLessOrEqual(DotProduct(p1, p3, p1, p4), 0)) { p = p1; return 1; }
  if (FloatEqual(s2, 0) && FloatLessOrEqual(DotProduct(p2, p3, p2, p4), 0)) { p = p2; return 2; }
  if (FloatEqual(s3, 0) && FloatLessOrEqual(DotProduct(p3, p1, p3, p2), 0)) { p = p3; return 3; }
  if (FloatEqual(s4, 0) && FloatLessOrEqual(DotProduct(p4, p1, p4, p2), 0)) { p = p4; return 4; }
  return 0;
}

多边形面积

对任何多边形都适用。

1/2 | (OP1 × OP2) * (OP2 × OP3) * … * (OPn × OP1) |

2 - Convex Hull

我只会两种方法，一种 Graham ，一种单调链。

凸包不会单独考，都是用来奠基的，正常打随便打都可以，但是有些脑残题目会要你保留凸包上每个点（包括共线的），比如POJ1228，这个时候就要好好考虑了。据我所知 Graham 需要一个小麻烦的特判，所以这题我打了单调链。

Graham

我尽量封装的好点，但是我发现比较函数不能放到 class 里面，然后就跑去 StackOverflow 上面问了下，除非 C++ 11 ，否则无解，那就只能破坏封装了。 http://stackoverflow.com/questions/13801311/compare-function-in-a-class-for-stdsort-cant-compile

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

Point base;
inline bool cmp(const Point& a, const Point& b)
{
  static int tmp;
  return (tmp = CrossProduct(base, a, base, b)) < 0 || (tmp == 0 && Dist(base, a) < Dist(base, b));
}
class ConvexHull
{
  Point p[MAXN];
public:
  int size;
  void Create(Point (&)[MAXN], const int);
  double Perimeter() const
  {
    double m = .0;
    for (int i = 1; i < size; ++i)
      m += std::sqrt(Dist(p[i-1], p[i]));
    return m;
  }
};
void ConvexHull::Create(Point (&a)[MAXN], const int n)
{
  base = a[0];
  for (int i = 1; i < n; ++i)
    if (a[i].x < base.x || (a[i].x == base.x && a[i].y < base.y))
      base = a[i];
  std::sort(a, a+n, cmp);
  p[0] = a[0], p[1] = a[1];
  size = 2;
  for (int i = 2; i < n; ++i)
  {
    while (size >= 2 && CrossProduct(a[i], p[size-1], a[i], p[size-2]) <= 0) --size;
    p[size++] = a[i];
  }
  p[size++] = p[0];
}

单调链

去年我一直打单调链，但是看到的标程都是 Graham ，于是今年决定全部打 Graham。但是郁闷的事情是，我的 Graham 被POJ1228卡掉了，只好老老实实回来打单调链。其实单调链除了速度慢了点，好处还是很明显的：非常好打。下面的程序是POJ1228的。

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

inline bool cmp(const Point& lhs, const Point& rhs) { return lhs.x < rhs.x || (lhs.x == rhs.x && lhs.y < rhs.y); }
class ConvexHull
{
  Point p[MAXN+1];
  void add(const Point& x)
  {
    while (size >= 2 && CrossProduct(p[size-1], x, p[size-2], x) > 0) --size;
    p[size++] = x;
  }
public:
  int size;
  void Create(Point (&a)[MAXN], const int n)
  {
    std::sort(a, a+n, cmp);
    p[0] = a[0], p[1] = a[1];
    size = 2;
    for (int i = 2; i < n; ++i) add(a[i]);
    for (int i = n-2; i >= 0; --i) add(a[i]);
  }
  bool Unique() const
  {
    for (int i = 0, j = 1; j < size; i = j-1)
    {
      while (j < size && CrossProduct(p[i], p[i+1], p[i], p[j]) == 0) ++j;
      if (j-i < 3)
        return false;
    }
    return true;
  }
};

Comments