3D计算几何的旋转与求交

被faebdc的shadow坑惨了。在考场上才发现自己的3D计算几何知识基本为0。

旋转

已知点\(P=(x,y,z)\),把它沿着\(\overrightarrow a=(x_0,y_0,z_0)\)转动\(\theta^{\circ}\)(从与向量相对的角度看过去,点是逆时针旋转的) 首先可以给出绕三个轴的旋转矩阵 \[R_x(\theta)=\begin{bmatrix} 1&0&0\\ 0&\cos\theta&-\sin\theta\\ 0&\sin\theta&\cos\theta \end{bmatrix}\] \[R_y(\theta)=\begin{bmatrix} \cos\theta&0&\sin\theta\\ 0&1&0\\ -\sin\theta&0&\cos\theta \end{bmatrix}\] \(R_y\)有些不同,因为\(R_y\)的变化中\(z\)等价成了二维的\(x\)轴,\(x\)被等价成了\(y\)\[R_z(\theta)=\begin{bmatrix} \cos\theta&-\sin\theta&0\\ \sin\theta&\cos\theta&0\\ 0&0&1 \end{bmatrix}\]

对于这三个特殊的旋转轴我们已经有了解决方案了,那么对于任意的轴\(\overrightarrow a\)要怎么办呢?我们可以将这个旋转分解: 将整个坐标轴旋转,使得旋转轴\(\overrightarrow a\)\(z\)轴重合 再将点\(P\)\(z\)轴旋转\(\theta\)角 再将整个坐标轴旋转回原位

算了直接上结论。。Miskcoo \[\begin{bmatrix} \cos \theta + x^2(1 - \cos \theta) & xy(1 - \cos \theta) - z\sin \theta & xz(1 - \cos \theta) + y\sin \theta \\ yx(1 - \cos \theta) + z\sin \theta & \cos \theta + y^2(1 - \cos \theta) & yz(1 - \cos \theta) - x\sin \theta \\ zx(1 - \cos \theta) - y\sin \theta & zy(1 - \cos \theta) + x\sin \theta & \cos \theta + z^2(1 - \cos \theta) \end{bmatrix}\] 还是很优美的。

线面相交

设平面的法向量为\(\overrightarrow n\),平面上有一点\(P\)。有直线\(l=Q+t\overrightarrow v\)

设交点为\(C\),则

\[\left\{ \begin{aligned} &(C-P)\cdot \overrightarrow n=0\\ &C=Q+t'\overrightarrow v \end{aligned} \right.\] 解得\(t'=\dfrac {(P-Q)\cdot\overrightarrow n}{\overrightarrow v\cdot\overrightarrow n}\)

\(\overrightarrow v\cdot\overrightarrow n\)的时候,显然交点不存在。

这种法线的操作对于二维的线线相交依然适用,比面积法或者体积法显然得多,好记得多,而且数值稳定性并没有损失。

一道题

附上一道练手题目:Shadow

这道题目存在更简洁的做法,但是为了强行套模板,我选择了一种非常难受的转化:把坐标系强行转化成地面是坐标平面,转化成把\((0,0,0),B,C\)构成的三角形转到\(xOy\)面上。详见代码,实在是恶心得不想多说。

写完1A好评。

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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
#include <cstdio>
struct Istream {
template <class T>
Istream &operator >>(T &x) {
static char ch;static bool neg;
for(ch=neg=0;ch<'0' || '9'<ch;neg|=ch=='-',ch=getchar());
for(x=0;'0'<=ch && ch<='9';(x*=10)+=ch-'0',ch=getchar());
x=neg?-x:x;
return *this;
}
Istream &operator >>(double &x) {
scanf("%lf",&x);
return *this;
}
}is;
struct Ostream {
template <class T>
Ostream &operator <<(T x) {
x<0 && (putchar('-'),x=-x);
static char stack[233];static int top;
for(top=0;x;stack[++top]=x%10+'0',x/=10);
for(top==0 && (stack[top=1]='0');top;putchar(stack[top--]));
return *this;
}
Ostream &operator <<(double const &x) {
printf("%.2lf",x);
return *this;
}
Ostream &operator <<(char ch) {
putchar(ch);
return *this;
}
}os;
#include <cmath>
#include <cassert>
#include <algorithm>
const int MAXN=100+11;
const double eps=1e-7,pi=acos(-1);
inline int sgn(double const &x) { return (x>eps)-(x<-eps); }
#define SELF_OPER(Type,op) template <class T> Type &operator op##=(T const &x) { return *this=*this op x; }
namespace _3D {
typedef struct Point Vector;
struct Point {
double x,y,z;
Point() {}
Point(double const &x,double const &y,double const &z):x(x),y(y),z(z) {}
double Length() const { return sqrt(x*x+y*y+z*z); }
friend double Inner(Vector const &a,Vector const &b) { return a.x*b.x+a.y*b.y+a.z*b.z; }
friend Vector Outer(Vector const &a,Vector const &b) { return Vector(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x); }
friend Point operator +(Point const &lhs,Vector const &rhs) { return Point(lhs.x+rhs.x,lhs.y+rhs.y,lhs.z+rhs.z); }
friend Vector operator -(Point const &lhs,Point const &rhs) { return Vector(lhs.x-rhs.x,lhs.y-rhs.y,lhs.z-rhs.z); }
friend Vector operator *(Vector const &lhs,double const &rhs) { return Vector(lhs.x*rhs,lhs.y*rhs,lhs.z*rhs); }
friend Vector operator /(Vector const &lhs,double const &rhs) { return Vector(lhs.x/rhs,lhs.y/rhs,lhs.z/rhs); }
SELF_OPER(Point,+)
SELF_OPER(Point,-)
SELF_OPER(Point,*)
SELF_OPER(Point,/)
friend double Angle(Vector const &v1,Vector const &v2) { return asin(Outer(v1,v2).Length()/(v1.Length()*v2.Length())); }
friend Istream &operator >>(Istream &is,Point &p) {
is>>p.x>>p.y>>p.z;
return is;
}
friend Ostream &operator <<(Ostream &os,Point const &p) {
os<<'('<<p.x<<','<<p.y<<','<<p.z<<')';
return os;
}
}p[MAXN];
struct Line {
Point p;
Vector v;
Line() {}
Line(Point const &p,Vector const &v):p(p),v(v) {}
};
Point Rotate(Line const &ln,Point p,double const &rad) {
static double co,si;
static Vector a;
si=sin(rad),co=cos(rad);
a=ln.v/ln.v.Length();
double T[][3]={{co+a.x*a.x*(1-co),a.x*a.y*(1-co)-a.z*si,a.x*a.z*(1-co)+a.y*si},
{a.y*a.x*(1-co)+a.z*si,co+a.y*a.y*(1-co),a.y*a.z*(1-co)-a.x*si},
{a.z*a.x*(1-co)-a.y*si,a.z*a.y*(1-co)+a.x*si,co+a.z*a.z*(1-co)}};
p-=ln.p;
p=Point(p.x*T[0][0]+p.y*T[1][0]+p.z*T[2][0],
p.x*T[0][1]+p.y*T[1][1]+p.z*T[2][1],
p.x*T[0][2]+p.y*T[1][2]+p.z*T[2][2]);
return p+=ln.p;
}
struct Plane {
Point p;Vector nv;
Plane(Point const &p1,Point const &p2,Point const &p3):p(p1),nv(Outer(p2-p1,p3-p1)) {}
};
const Plane GROUND(Point(0,0,0),Point(1,0,0),Point(0,1,0));
Point Cross(Line ln,Plane pln) {
double d1=-Inner(pln.nv,ln.p-pln.p),d2=Inner(pln.nv,ln.v);//d1/d2
if(d2==0) {
throw;
} else {
return ln.p+ln.v*(d1/d2);
}
}
}
namespace _2D {
typedef struct Point Vector;
struct Point {
double x,y;
Point() {}
Point(double const &x,double const &y):x(x),y(y) {}
friend bool operator <(Point const &lhs,Point const &rhs) { return lhs.x==rhs.x?lhs.y<rhs.y:lhs.x<rhs.x; }
friend Vector operator +(Point const &lhs,Vector const &rhs) { return Vector(lhs.x+rhs.x,lhs.y+rhs.y); }
friend Vector operator -(Point const &lhs,Point const &rhs) { return Vector(lhs.x-rhs.x,lhs.y-rhs.y); }
friend Vector operator *(Vector const &lhs,double const &rhs) { return Vector(lhs.x*rhs,lhs.y*rhs); }
friend Vector operator /(Vector const &lhs,double const &rhs) { return Vector(lhs.x/rhs,lhs.y/rhs); }
friend double Inner(Vector const &a,Vector const &b) { return a.x*b.x+a.y*b.y; }
friend double Outer(Vector const &a,Vector const &b) { return a.x*b.y-a.y*b.x; }
}p[MAXN];
struct Line {
Point p;Vector v;
Line() {}
Line(Point const &p,Vector const &v):p(p),v(v) {}
};
Point Cross(Line const &ln1,Line const &ln2) {
Vector norm(-ln2.v.y,ln2.v.x);
double d1=-Inner(norm,ln1.p-ln2.p),d2=Inner(ln1.v,norm);
return ln1.p+ln1.v*(d1/d2);
}
double HullSize(Point p[],int n) {
std::sort(p+1,p+1+n);
static Point stack[MAXN];
int top=0;
for(int i=1;i<=n;++i) {
while(top>=2 && Outer(stack[top]-stack[top-1],p[i]-stack[top])<=0) {
top--;
}
stack[++top]=p[i];
}
int mid=top;
for(int i=n;i>=1;--i) {
while(top>=mid+1 && Outer(stack[top]-stack[top-1],p[i]-stack[top])<=0) {
top--;
}
stack[++top]=p[i];
}
if(top) {
top--;
}
double res=0;
for(int i=1;i<=top;++i) {
res+=Outer(stack[i],stack[i==top?1:i+1]);
}
return res/2;
}
}
using namespace _3D;
void Transform(Point g[],Point &light,Point p[],int n) {
Point &o=g[0];
for(int i=1;i<=n;++i) {
p[i]-=o;
}
light-=o;
g[2]-=o;
g[1]-=o;
g[0]-=o;
Line axis;double angle;
axis=Line(o,Vector(-g[1].y,g[1].x,0));
angle=sgn(g[1].x)==0 && sgn(g[1].y)==0?pi/2:Angle(g[1],Vector(g[1].x,g[1].y,0));
if(sgn(Rotate(axis,g[1],angle).z)!=0) {
angle*=-1;
}
for(int i=1;i<=n;++i) {
p[i]=Rotate(axis,p[i],angle);
}
light=Rotate(axis,light,angle);
g[1]=Rotate(axis,g[1],angle);
g[2]=Rotate(axis,g[2],angle);
assert(sgn(g[1].z)==0);
_2D::Point _crs=_2D::Cross(_2D::Line(_2D::Point(g[2].x,g[2].y),_2D::Vector(-g[1].y,g[1].x)),_2D::Line(_2D::Point(0,0),_2D::Point(g[1].x,g[1].y)));
Point crs=Point(_crs.x,_crs.y,0);
axis=Line(o,g[1]);
angle=sgn(g[2].x-crs.x)==0 && sgn(g[2].y-crs.y)==0?pi/2:Angle(g[2]-crs,Point(g[2].x,g[2].y,0)-crs);
if(sgn(Rotate(axis,g[2],angle).z)!=0) {
angle*=-1;
}
for(int i=1;i<=n;++i) {
p[i]=Rotate(axis,p[i],angle);
}
light=Rotate(axis,light,angle);
g[2]=Rotate(axis,g[2],angle);
assert(sgn(g[2].z)==0);
}
int main() {
#ifdef DEBUG
freopen("input","r",stdin);
#else
freopen("shadow.in","r",stdin);
freopen("shadow.out","w",stdout);
#endif
Point g[3];is>>g[0]>>g[1]>>g[2];
Point light;is>>light;
int n;is>>n;
for(int i=1;i<=n;++i) {
is>>p[i];
}
Transform(g,light,p,n);
Plane pln(g[0],g[1],g[2]);
for(int i=1;i<=n;++i) {
Point crs=Cross(Line(light,p[i]-light),pln);
_2D::p[i]=_2D::Point(crs.x,crs.y);
}
os<<_2D::HullSize(_2D::p,n)<<'\n';
return 0;
}

Powered by Hexo and Hexo-theme-hiker

Copyright © 2013 - 2017 Hey,Lucida here. All Rights Reserved.

Lucida Lu 保留所有权利