计算几何：两条直线的位置关系

       想学习一点计算几何的知识，估计需要重新把中学的解析几何捡起来才行。丢了那么多年了，想捡并不容易。晚饭后和老婆出去散步，碰到个老同事，原本搞软件开发的，做得相当有成就。闲聊中得知，人家已经走上管理岗位，早不再写程序了。我都有点不好意思开口说，自己为了教学，还在从零开始，学一些以前不会或者忘了的东东。

        虽然差不多是从头再来，但多少有些印象，知道什么是简单的，好入门的。于是就选择了这么一个最容易的的问题：给定两条直线，判断它们的位置关系。如果相交，则求出交点。

        由高中数学知识可知：两条直线有位置关系有三种：相交、平行、重合。

        （1）若两条直线都存在斜率，则当斜率不相等时，两直线相交；斜率相等，则可根据两直线在Y轴的截距是否相等判断是平行还是重合：截距相等则重合。斜率相等时，直线重合也可根据两条直线是否有公共点来确定。

        （2）若两条直线都没有斜率，则不重合时，就平行。

        （3）若一条直线有斜率，而另一条没有斜率，则相交。

        PKU1269 Intersecting Lines 正是一道直白的题目。

        题目大意：给出N组直线，每组两条直线，每条直线给出2个点的坐标。判断这两条直线的位置关系。

        分析：由给定2个点的坐标，可以判定并求出这条直线的斜率。之后，依据两条直线斜率的有无，分成4种情况分别处理。为了判定重合，需要计算出直线在Y轴的斜率。我用的是y = kx + b斜截式方程的变形：b = y - kx。计算两条直线的交点，相当于联立求解一个二元一次方程。

        根据前面的介绍，可以写出一个基本的，可以AC的程序。参考程序：

/*
Problem: 1269		
Memory: 208K		Time: 0MS
Language: C++		Result: Accepted
*/
#include<cstdio>
#include<cmath>
#include<string>
using namespace std;
#define eps	 1E-8
#define zero(x) (((x)>0?(x):-(x)) < eps) 

struct point {
	double x, y;
};
struct line {
	point a, b;
};

bool get_k(line u, double & k)		//计算直线的斜率k，若不存在，返回false
{
	if(u.a.x == u.b.x) return false;
	k = (u.a.y - u.b.y) / (u.a.x - u.b.x);
	return true;
}

double get_b(line u, double k)		//计算直线在Y轴的截距
{
	double b = u.a.y - k * u.a.x;
	return b;
}

point intersection(line u, line v)
{
	point ret = u.a;
	double t = ((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))
			 / ((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));
	ret.x += (u.b.x-u.a.x)*t;
	ret.y += (u.b.y-u.a.y)*t; 
	return ret;
}

int main()
{
	int n;
	line u, v;
	point res;
	double ku, kv, bu, bv;
	bool fu, fv;
	string ans;
	
	puts("INTERSECTING LINES OUTPUT");
	
	scanf("%d", &n);
	while(n-- > 0)
	{
		scanf("%lf%lf%lf%lf", &u.a.x, &u.a.y, &u.b.x, &u.b.y);
		scanf("%lf%lf%lf%lf", &v.a.x, &v.a.y, &v.b.x, &v.b.y);
		fu = get_k(u, ku);
		fv = get_k(v, kv);
		if(fu && fv)					//两条直线都有斜率
		{
			if( zero(ku - kv ) )		//斜率相等
			{
				bu = get_b(u, ku);
				bv = get_b(v, kv);
			
				if ( zero( bu - bv ))	//截距相等
					ans = "LINE";
				else
					ans = "NONE";
			}
			else
			{
				res = intersection(u, v);
				ans = "POINT";
			}
		}
		else if( !fu && ! fv)			//都没有斜率
		{
			if( zero(u.a.x - v.a.x) )	//如果X轴截距相等，则重合
				ans = "LINE";
			else
				ans = "NONE";
		}
		else if ( !fu )					//u没有斜率
		{
			bv = get_b(v, kv);
			res.x = u.a.x;
			res.y = kv * u.a.x + bv;
			ans = "POINT";
		}
		else							//v没有斜率
		{
			bu = get_b(u, ku);
			res.x = v.a.x;
			res.y = ku * v.a.x + bu;
			ans = "POINT";
		}
		if(ans == "POINT")
			printf("POINT %.2lf %.2lf\n", res.x, res.y);
		else
			printf("%s\n", ans.c_str());
	}
	puts("END OF OUTPUT");
	return 0;
}

        写好程序，通过样例后，提交到POJ，却WA掉了。检查程序，并对照别人AC的程序，没发现区别。心里很不是滋味。难道我第一次争战计算几何就出师不利？无奈中打开版上的DISCUSS，却发现WA声一片。有人叫道：“poj的G++到底是咋回事？每次WA后换成C++就能过.....”我按其提示把语言换成C++，果然AC。不过为何G++过不了还是不得而知。是double型的输入输出函数的格式控制符%lf要换成%f吗？我试过，还是不行。等下次我试试用cin, cout，看行不？

        更好的，更漂亮的方法是利用向量的叉积作判断。其背后的原理是什么，我现在不懂。留个链接：

http://www.cppblog.com/Felicia/archive/2007/08/21/30535.html

 

        上面那个求两条直线交点的代码是照抄ACM模板的，我没有亲手推导一下。在国内OI比赛中是不允许带任何资料进场的，所以要准确记住这么一大坨代码好难。我就想，抽空推导一下吧，兴许记得住。然后我从直线的两点式方程开始，变形，再变形，.......，最后好不容易把x坐标的解表示出来了，而y坐标的解实在不知该如何化简了。后来又想，可能先利用两点式方程求出解析式，即把一般方程ax + by + c = 0的a, b, c三个系数先求出来，再联立解二元一次方程就容易了。推导了一下，果然成功。最后突然认识到，既然已经求出了斜率和截距，利用最简单的斜截式方程，不是很容易表示出交点的解吗？ 实际编码试了一下，果然可行。代码如下：

point intersection(double k1, double b1, double k2, double b2)	//斜截方程求交点
{
	point ret;
	ret.x = (b2 - b1) / (k1 - k2);
	ret.y = (b2 * k1 - b1 * k2) / (k1 - k2);
	return ret;
}