UVa12303 Composite Transformations

题目链接

  UVa12303 Composite Transformations

题意

  空间中有n个点和m个平面，你的任务是按顺序向它们施加t个变换，输出每个点的最终位置和每个平面的最终方程。一共有3种变换，如表下表所示。

变换	说明
TRANSLATE a b c	点(x,y,z)变成(x+a,y+b,z+c)
ROTATE a b c theta	把每个点旋转theta 度。旋转轴是向量(a,b,c)，并且穿过原点。当旋转轴指向你时，旋转看上去是逆时针的
SCALE a b c	点(x,y,z)变成(ax,by,cz)

输入格式

  输入只有一组数据。第一行是3个整数n, m, t（1≤n，m≤50 000，1≤t≤1 000）；以下n行每行包含一个点的坐标；用以下m行每行4个整数a, b, c, d，描述一个平面ax+by+cz+d=0（a, b, c 不全为0）；以下t行描述各个操作。操作中的参数a, b, c, d 是绝对值不超过10 的实数，最多保留小数点后两位。参数theta 是0~359 之间的整数。

输出格式

  对于每个点，输出一行，即变换后的坐标。对于每个平面，输出一行，即变换后平面方程的参数a, b, c 和d。输出应保证$a^2+b^2+c^2=1$，但如果有多种表示方法，任意一种均可。输出应保留两位小数。本题允许一定的浮点误差。

分析

  先说几个坑点：1、题面交代输入用4个整数a, b, c, d描述一个平面，实际数据会有浮点数，所以应该用double接收，否则WA；2、题面说输出接受0.05的误差，可能是相对误差（即5%），因为生成测试数据后Udebug输出的结果和lrj仓库代码输出的结果相差挺大的，所以果断全程double来算就行了。

  三个仿射变换矩阵中，平移、缩放矩阵好写，难的是旋转矩阵，具体可参考三维仿射变换矩阵。
  绕原点且单位法向量为$(a,b,c)$的轴旋转的变换矩阵为：
$$\left[\begin{array}{cccc}\cos \theta +a^2(1-\cos \theta )& ab(1-\cos \theta )-c\sin \theta & ac(1-\cos \theta )+b\sin \theta & 0\\ ab(1-\cos \theta )+c\sin \theta & \cos \theta +b^2(1-\cos \theta )& bc(1-cos\theta )-a\sin \theta & 0\\ ac(1-\cos \theta )-b\sin \theta & bc(1-\cos \theta )+a\sin \theta & \cos \theta +c^2(1-\cos \theta )& 0\\ 0& 0& 0& 1\end{array}\right]$$
  《训练指南》建议用三点式求解平面的变换结果，那么需要解决如下问题：已知平面方程，如何让找出平面上不共线的三个点？具体可参考平面的四种方程及一些应用。
  已知平面的一般式方程$\mathbf{ax}+\mathbf{by}+\mathbf{cz}+\mathbf{d}=0$，由于$\mathbf{a},\mathbf{b},\mathbf{c}不全为0$，则：
  $\mathbf{a}\ne 0$时，点$(-\frac{\mathbf{d}}{\mathbf{a}},0,0),(-\frac{\mathbf{c}+\mathbf{d}}{\mathbf{a}},0,1),(-\frac{\mathbf{b}+\mathbf{d}}{\mathbf{a}},1,0)$在平面上且三点不共线；
  $\mathbf{b}\ne 0$时，点$(0,-\frac{\mathbf{d}}{\mathbf{b}},0),(0,-\frac{\mathbf{c}+\mathbf{d}}{\mathbf{b}},1),(1,-\frac{\mathbf{a}+\mathbf{d}}{\mathbf{b}},0)$在平面上且三点不共线；
  $\mathbf{c}\ne 0$时，点$(0,0,-\frac{\mathbf{d}}{\mathbf{c}}),(0,1,-\frac{\mathbf{b}+\mathbf{d}}{\mathbf{c}}),(1,0,-\frac{\mathbf{a}+\mathbf{d}}{\mathbf{c}})$在平面上且三点不共线。
  编程时推荐选择$\mathbf{a},\mathbf{b},\mathbf{c}$中绝对值最大那个来构造点。

AC 代码

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

#define N 50010
double pl[N][4], p[N][3], l[3][4], r[4][4], v[3][4], e = M_PI/180;

double mul(const double (&r)[4], const double *v3) {
    return r[0]*v3[0] + r[1]*v3[1] + r[2]*v3[2] +  + r[3];
}

void load_trans(double a, double b, double c) {
    l[0][0] = 1.; l[0][1] = 0.; l[0][2] = 0.; l[0][3] = a;
    l[1][0] = 0.; l[1][1] = 1.; l[1][2] = 0.; l[1][3] = b;
    l[2][0] = 0.; l[2][1] = 0.; l[2][2] = 1.; l[2][3] = c;
}

void load_rot(double a, double b, double c, double theta) {
    double si = sin(theta), co = cos(theta), s = sqrt(a*a+b*b+c*c); a /= s; b /= s; c /= s;
    l[0][0] = co+a*a*(1-co); l[0][1] = a*b*(1-co)-c*si; l[0][2] = a*c*(1-co)+b*si; l[0][3] = 0.;
    l[1][0] = a*b*(1-co)+c*si; l[1][1] = co+b*b*(1-co); l[1][2] = b*c*(1-co)-a*si; l[1][3] = 0.;
    l[2][0] = a*c*(1-co)-b*si; l[2][1] = b*c*(1-co)+a*si; l[2][2] = co+c*c*(1-co); l[2][3] = 0.;
}

void load_scale(double a, double b, double c) {
    l[0][0] = a; l[0][1] = 0.; l[0][2] = 0.; l[0][3] = 0.;
    l[1][0] = 0.; l[1][1] = b; l[1][2] = 0.; l[1][3] = 0.;
    l[2][0] = 0.; l[2][1] = 0.; l[2][2] = c; l[2][3] = 0.;
}

void trans_plane(double a, double b, double c, double d) {
    double x = max(max(abs(a), abs(b)), abs(c));
    if (abs(a) == x) {
        v[0][0] = -d/a; v[0][1] = 0.; v[0][2] = .0;
        v[1][0] = -(c+d)/a; v[1][1] = 0.; v[1][2] = 1.;
        v[2][0] = -(b+d)/a; v[2][1] = 1.; v[2][2] = 0.;
    } else if (abs(b) == x) {
        v[0][0] = 0.; v[0][1] = -d/b;  v[0][2] = 0.;
        v[1][0] = 0.; v[1][1] = -(c+d)/b; v[1][2] = 1.;
        v[2][0] = 1.; v[2][1] = -(a+d)/b; v[2][2] = 0.;
    } else {
        v[0][0] = 0.; v[0][1] = 0.; v[0][2] = -d/c;
        v[1][0] = 0.; v[1][1] = 1.; v[1][2] = -(b+d)/c;
        v[2][0] = 1.; v[2][1] = 0.; v[2][2] = -(a+d)/c;
    }
    for (int i=0; i<3; ++i)
        l[i][0] = mul(r[0], v[i]), l[i][1] = mul(r[1], v[i]), l[i][2] = mul(r[2], v[i]);
    double x1 = l[1][0]-l[0][0], y1 = l[1][1]-l[0][1], z1 = l[1][2]-l[0][2],
        x2 = l[2][0]-l[0][0], y2 = l[2][1]-l[0][1], z2 = l[2][2]-l[0][2],
        e = y1*z2 - z1*y2, f = x2*z1 - x1*z2, g = x1*y2 - x2*y1, s = sqrt(e*e + f*f + g*g);
        e /= s; f /= s; g /= s;
    cout << e << ' ' << f << ' ' << g << ' ' << -(e*l[0][0] + f*l[0][1] + g*l[0][2]) << endl;
}

void solve() {
    int m, n, t; char s[10]; double a, b, c, d; cin >> n >> m >> t;
    for (int i=0; i<n; ++i) cin >> p[i][0] >> p[i][1] >> p[i][2];
    for (int i=0; i<m; ++i) cin >> pl[i][0] >> pl[i][1] >> pl[i][2] >> pl[i][3];
    for (int i=0; i<4; ++i) for (int j=0; j<4; ++j) r[i][j] = i==j ? 1. : 0.;
    while (t--) {
        cin >> s >> a >> b >> c;
        if (s[0] == 'T') load_trans(a, b, c);
        else if (s[0] == 'R') cin >> d, load_rot(a, b, c, d*e);
        else load_scale(a, b, c);
        for (int i=0; i<3; ++i) for (int j=0; j<4; ++j)
            v[i][j] = l[i][0]*r[0][j] + l[i][1]*r[1][j] + l[i][2]*r[2][j] + l[i][3]*r[3][j];
        for (int i=0; i<3; ++i) for (int j=0; j<4; ++j) r[i][j] = v[i][j];
    }
    for (int i=0; i<n; ++i) cout << mul(r[0], p[i]) << ' ' << mul(r[1], p[i]) << ' ' << mul(r[2], p[i]) << endl;
    for (int i=0; i<m; ++i) trans_plane(pl[i][0], pl[i][1], pl[i][2], pl[i][3]);
}

int main() {
    ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
    cout << fixed << setprecision(2);
    solve();
    return 0;
}