2013-team5/linetoline
从 Trac 迁移的文章
这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。
原文章内容如下：
{{{
void linetoline(const point3 &a, const point3 &b, const point3 &c, const point3 &d, point3 &p1, point3 &p2){
point3 ab=b-a, cd=d-c, ac=c-a, s = ab*cd;
double t1 = ((ac*cd)^s) / s.sqrlen();
double t2 = ((ac*ab)^s) / s.sqrlen();
p1 = a+ab*t1;
p2 = c+cd*t2;
}
//精度更高
}}}
{{{
double line_to_line(point3 u1,point3 u2,point3 v1,point3 v2,point3& p1,point3& p2){
point3 s0=(u2-u1)*(v2-v1); //两直线法向量
point3 tmp=s0*(s0^(v1-u1))/s0.sqrlen();
p2=intersection(u1+tmp,u2+tmp,v1,v2); //将l1通过法向量方向投影到l2所在平面,得到l3. p2即为l2与l3交点
p1=p2-tmp; //将p2通过法向量方向投影到l1上,得到p1
return (p1-p2).length(); //返回p1,p2之间的距离
}
//u1,u2为l1上不同的两点,v1,v2为l2上不同的两点
//p1,p2为所求线段的端点,其中p1在l1上,p2在l2上,返回值为该线段长度.
}}}
{{{
///2013杭州赛区网络赛D题，hdu4741
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const int MAXN = 1000;
const double eps = 1e-8;
const double PI = atan2(0.0, -1.0);
inline double sqr(double x){ return x * x; }
inline bool zero(double x){ return (x > 0 ? x : -x) < eps; }
inline int sgn(double x){ return (x > eps ? 1 : (x + eps < 0 ? -1 : 0)); }
struct point3{
double x, y, z;
point3(double x, double y, double z):x(x), y(y), z(z){}
point3() {}
bool operator == (const point3 & a) const{ return sgn(x - a.x) == 0 && sgn(y - a.y) == 0 && sgn(z - a.z) == 0; }
bool operator != (const point3 & a) const{ return sgn(x - a.x) != 0 || sgn(y - a.y) != 0 || sgn(z - a.z) != 0; }
bool operator < (const point3 & a) const{ return sgn(x - a.x) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) == 0 && sgn(z - a.z) < 0; }
point3 operator + (const point3 & a) const{ return point3(x + a.x, y + a.y, z + a.z); }
point3 operator - (const point3 & a) const{ return point3(x - a.x, y - a.y, z - a.z); }
point3 operator * (const double & a) const{ return point3(x * a, y * a, z * a); }
point3 operator / (const double & a) const{ return point3(x / a, y / a, z / a); }
point3 operator * (const point3 & a) const{ return point3(y * a.z - z * a.y, z * a.x - x * a.z, x * a.y - y * a.x); }//xmult
double operator ^ (const point3 & a) const{ return x * a.x + y * a.y + z * a.z; }//dmult
double sqrlen() const{ return sqr(x) + sqr(y) + sqr(z); }
double length() const{ return sqrt(sqrlen()); }
point3 trunc(double a) const{ return (*this) * (a / length()); }
point3 rotate(const point3 & a, const point3 & b, const point3 & c) const{ return point3(a ^ (*this), b ^ (*this), c ^ (*this)); }//abc正交且模为1
};
point3 intersection(const point3 & u1, const point3 & u2, const point3 & v1, const point3 & v2){
point3 p0 = (u1 - v1) * (v1 - v2), p1 = (u1 - u2) * (v1 - v2);
return u1 + (u2 - u1) * (sgn(p0 ^ p1) * sqrt(p0.sqrlen() / p1.sqrlen()));
}
double line_to_line(point3 u1,point3 u2,point3 v1,point3 v2,point3& p1,point3& p2){
point3 s0=(u2-u1)*(v2-v1);
point3 tmp=s0*(s0^(v1-u1))/s0.sqrlen();
p2=intersection(u1+tmp,u2+tmp,v1,v2);
p1=p2-tmp;
return (p1-p2).length();
}
int main(){
int T;
scanf("%d",&T);
while (T--){
point3 p[2][2],a1,a2;
int i,j,k;
for (i=0; i<2; i++)
for (j=0; j<2; j++)
scanf("%lf%lf%lf",&p[i][j].x,&p[i][j].y,&p[i][j].z);
double res=line_to_line(p[0][0],p[0][1],p[1][0],p[1][1],a1,a2);
printf("%.6lf\n",res);
printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf\n",a1.x,a1.y,a1.z,a2.x,a2.y,a2.z);
}
return 0;
}
}}}void linetoline(const point3 &a, const point3 &b, const point3 &c, const point3 &d, point3 &p1, point3 &p2){
point3 ab=b-a, cd=d-c, ac=c-a, s = ab*cd;
double t1 = ((ac*cd)^s) / s.sqrlen();
double t2 = ((ac*ab)^s) / s.sqrlen();
p1 = a+ab*t1;
p2 = c+cd*t2;
}
//精度更高
double line_to_line(point3 u1,point3 u2,point3 v1,point3 v2,point3& p1,point3& p2){
point3 s0=(u2-u1)*(v2-v1); //两直线法向量
point3 tmp=s0*(s0^(v1-u1))/s0.sqrlen();
p2=intersection(u1+tmp,u2+tmp,v1,v2); //将l1通过法向量方向投影到l2所在平面,得到l3. p2即为l2与l3交点
p1=p2-tmp; //将p2通过法向量方向投影到l1上,得到p1
return (p1-p2).length(); //返回p1,p2之间的距离
}
//u1,u2为l1上不同的两点,v1,v2为l2上不同的两点
//p1,p2为所求线段的端点,其中p1在l1上,p2在l2上,返回值为该线段长度.
///2013杭州赛区网络赛D题，hdu4741
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const int MAXN = 1000;
const double eps = 1e-8;
const double PI = atan2(0.0, -1.0);
inline double sqr(double x){ return x * x; }
inline bool zero(double x){ return (x > 0 ? x : -x) < eps; }
inline int sgn(double x){ return (x > eps ? 1 : (x + eps < 0 ? -1 : 0)); }
struct point3{
double x, y, z;
point3(double x, double y, double z):x(x), y(y), z(z){}
point3() {}
bool operator == (const point3 & a) const{ return sgn(x - a.x) == 0 && sgn(y - a.y) == 0 && sgn(z - a.z) == 0; }
bool operator != (const point3 & a) const{ return sgn(x - a.x) != 0 || sgn(y - a.y) != 0 || sgn(z - a.z) != 0; }
bool operator < (const point3 & a) const{ return sgn(x - a.x) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) == 0 && sgn(z - a.z) < 0; }
point3 operator + (const point3 & a) const{ return point3(x + a.x, y + a.y, z + a.z); }
point3 operator - (const point3 & a) const{ return point3(x - a.x, y - a.y, z - a.z); }
point3 operator * (const double & a) const{ return point3(x * a, y * a, z * a); }
point3 operator / (const double & a) const{ return point3(x / a, y / a, z / a); }
point3 operator * (const point3 & a) const{ return point3(y * a.z - z * a.y, z * a.x - x * a.z, x * a.y - y * a.x); }//xmult
double operator ^ (const point3 & a) const{ return x * a.x + y * a.y + z * a.z; }//dmult
double sqrlen() const{ return sqr(x) + sqr(y) + sqr(z); }
double length() const{ return sqrt(sqrlen()); }
point3 trunc(double a) const{ return (*this) * (a / length()); }
point3 rotate(const point3 & a, const point3 & b, const point3 & c) const{ return point3(a ^ (*this), b ^ (*this), c ^ (*this)); }//abc正交且模为1
};
point3 intersection(const point3 & u1, const point3 & u2, const point3 & v1, const point3 & v2){
point3 p0 = (u1 - v1) * (v1 - v2), p1 = (u1 - u2) * (v1 - v2);
return u1 + (u2 - u1) * (sgn(p0 ^ p1) * sqrt(p0.sqrlen() / p1.sqrlen()));
}
double line_to_line(point3 u1,point3 u2,point3 v1,point3 v2,point3& p1,point3& p2){
point3 s0=(u2-u1)*(v2-v1);
point3 tmp=s0*(s0^(v1-u1))/s0.sqrlen();
p2=intersection(u1+tmp,u2+tmp,v1,v2);
p1=p2-tmp;
return (p1-p2).length();
}
int main(){
int T;
scanf("%d",&T);
while (T--){
point3 p[2][2],a1,a2;
int i,j,k;
for (i=0; i<2; i++)
for (j=0; j<2; j++)
scanf("%lf%lf%lf",&p[i][j].x,&p[i][j].y,&p[i][j].z);
double res=line_to_line(p[0][0],p[0][1],p[1][0],p[1][1],a1,a2);
printf("%.6lf\n",res);
printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf\n",a1.x,a1.y,a1.z,a2.x,a2.y,a2.z);
}
return 0;
}