计算几何学 | 投影 | Projection | C/C++实现

问题描述

对于给定的三个点p1、p2、p，从点p向通过p1、p2的直线引一条垂线，求垂足x的坐标。（点p在直线p1p2上的投影）

输入：
输入按以下格式给出：
$x_{p1}$xp1​ $y_{p1}$yp1​ $x_{p2}$xp2​ $y_{p2}$yp2​
$q$q
$x_{p0}$xp0​ $y_{p0}$yp0​
$x_{p1}$xp1​ $y_{p1}$yp1​
…
$x_{p_{q-1}}$xpq−1​​ $y_{p_{q-1}}$ypq−1​​
第1行给出p1、p2的坐标。接下来给出q个p的坐标用作问题。
输出：
根据各问题输出垂足x的坐标，每个问题占1行。输出允许误差不超过0.00000001。
限制：
1 ≤ q ≤ 1000
-10000 ≤ $x_i,y_i$xi​,yi​ ≤ 10000
p1、p2不是同一个点。

输入示例

0 0 3 4
1
2 5

输出示例

3.1200000000 4.1600000000

讲解

从点p向线段（或直线）s = p1p2引1条垂线，交点设为x。这个x就叫做点p的投影（projection)。

设s.p2 - s.p1为向量base，p - s.p1为向量hypo，点s.p1与点x的距离为t，hypo与x的距离为t，hypo与base的夹角为θ，则：
$t=|hypo|cosθ$t=∣hypo∣cosθ，$hypo·base=|hypo||base|cosθ$hypo⋅base=∣hypo∣∣base∣cosθ
于是有$t=hypo·base/|base|$t=hypo⋅base/∣base∣。
根据t与|base|的比例$r=t/|base|$r=t/∣base∣可得：
$x=s.p1+base*t/|base|=s.p1+base*(hypo·base)/|base|^2$x=s.p1+base∗t/∣base∣=s.p1+base∗(hypo⋅base)/∣base∣2
于是来看看如何用程序求点p在线段（直线）s上的投影。

点p在线段s上的投影：

Point project(Segment s, Point p) {
	Vector base = s.p2 - s.p1;
	double r = dot(p - s.p1, base) / norm(base);
	return s.p1 + base * r;
}

AC代码如下

class Point {//Point类，点 
	public:
		double x, y;
		
		Point(double x = 0, double y = 0): x(x), y(y) {}

		Point operator + (Point p) { return Point(x + p.x, y + p.y); }
		Point operator - (Point p) { return Point(x - p.x, y - p.y); }
		Point operator * (double a) { return Point(a * x, a * y); }
		Point operator / (double a) { return Point(x / a, y / a); }

		double abs() { return sqrt(norm()); }
		double norm() { return x * x + y * y; }
		
		bool operator < (const Point &p) const {
			return x != p.x ? x < p.x : y < p.y;
		}

		bool operator == (const Point &p) const {
			return fabs(x - p.x) < EPS && fabs(y - p.y) < EPS;
		}
};

typedef Point Vector;//Vector类，向量 

struct Segment{//Segment 线段 
	Point p1, p2;
};
double dot(Vector a, Vector b) {//内积 
	return a.x * b.x + a.y * b.y;
}

double cross(Vector a, Vector b) {//外积 
	return a.x*b.y - a.y*b.x;
}

bool isOrthogonal(Vector a, Vector b){//判断正交 
	return equals(dot(a, b), 0.0);
}

bool isOrthogonal(Point a1, Point a2, Point b1, Point b2){//判断正交 
	return isOrthogonal(a1 - a2, b1 - b2);
}

bool isOrthogonal(Segment s1, Segment s2){//判断正交 
	return equals(dot(s1.p2 - s1.p1, s2.p2 - s2.p1), 0.0);
}

bool isParallel(Vector a, Vector b){//判断平行 
	return equals(cross(a, b), 0.0);
}

bool isParallel(Point a1, Point a2, Point b1, Point b2){//判断平行 
	return isParallel(a1 - a2, b1 - b2);
}

bool isParallel(Segment s1, Segment s2){//判断平行 
	return equals(cross(s1.p2 - s1.p1, s2.p2 - s2.p1), 0.0);
}

Point project(Segment s, Point p) {//投影 对于给定的三个点p1、p2、p，从点p向通过
//p1、p2的直线引一条垂线，求垂足x的坐标。（点p在直线p1p2上的投影） 
	Vector base = s.p2 - s.p1;
	double r = dot(p - s.p1, base) / base.norm();
	return s.p1 + base * r;
} 

int main(){
	Point p1, p2, p;
	cin>>p1.x>>p1.y>>p2.x>>p2.y;
	
	Segment s;
	s.p1 = p1;
	s.p2 = p2;
	
	int q;
	cin>>q;
	
	while(q--){
		cin>>p.x>>p.y;
		p = project(s, p);
		printf("%.10f %.10f\n", p.x, p.y);
	}
}

对于代码中其他函数和类的详细解释，可参见：计算几何学