pku1263 光线反射 at GCCFeli.cn

pku1263 光线反射

评论30次阅读2008.08.26 11:00; 作者：Felicia　

这个题是 sgu110 的二维版本，我直接把 sgu110 那题的`z`分量设成0的，嘿嘿：）
光线在球面的反射可以用向量法求解，不需要解方程。

/**********************************************************************
Author: WHU_GCC
Created Time: 2008年08月04日星期一 10时24分48秒
File Name: pku1263.cpp
Description:
**********************************************************************/
#include <iostream>
#include <cmath>
using namespace std;
#define out(x) (cout << #x << ": " << x << endl)
const int maxint = 0x7FFFFFFF;
template <class T> void get_max(T &a, const T &b) {b > a ? a = b : 1;}
template <class T> void get_min(T &a, const T &b) {b < a ? a = b : 1;}

const double eps = 1e-9;
const double pi = acos(-1);

struct point {
   double x, y, z;
   point() {}
   point(double _x, double _y, double _z) : x(_x), y(_y), z(_z) {}
};

struct vec {
   double x, y, z;
   vec() {}
   vec(double _x, double _y, double _z) : x(_x), y(_y), z(_z) {}
   vec(point p) : x(p.x), y(p.y), z(p.z) {}
   vec(point s, point e) : x(e.x - s.x), y(e.y - s.y), z(e.z - s.z) {}
};

int sgn(double d) {
   return (d > eps) - (d < -eps);
}

vec cmul(const vec &a, const vec &b) {
   return vec(a.y * b.z - b.y * a.z, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
}

double dmul(const vec &a, const vec &b) {
   return a.x * b.x + a.y * b.y + a.z * b.z;
}

double sqr(const double &x) {
   return x * x;
}

double len(const vec &a) {
   return sqrt(sqr(a.x) + sqr(a.y) + sqr(a.z));
}

double len2(const vec &a) {
   return sqr(a.x) + sqr(a.y) + sqr(a.z);
}

double dist2_point_point(const point &a, const point &b) {
   return sqr(a.x - b.x) + sqr(a.y - b.y) + sqr(a.z - b.z);
}

double dist2_point_line(const point &p, const point &a, const point &b) {
   return abs(len2(cmul(vec(p, a), vec(p, b)))) / len2(vec(a, b));
}

int inter_circle_line(const point &c, double r, const point &p0, const vec &v0, point &p1, vec &v1) {
   double h2 = dist2_point_line(c, p0, point(p0.x + v0.x, p0.y + v0.y, p0.z + v0.z));
   int flag = sgn(h2 - sqr(r));
   if (flag > 0) return 0;

   double sh = dmul(v0, vec(p0, c));
   if (sh <= 0) return 0;

   double l2 = dist2_point_point(p0, c);
   double t = sqrt(l2 - h2) - sqrt(sqr(r) - h2);
   double lv0 = len(v0);
   vec vt(v0.x * t / lv0, v0.y * t / lv0, v0.z * t / lv0);
   p1 = point(p0.x + vt.x, p0.y + vt.y, p0.z + vt.z);

   double t2 = 2 * (dmul(vt, vec(p1, c)) / r);
   vec vtt((p1.x - c.x) * t2 / r, (p1.y - c.y) * t2 / r, (p1.z - c.z) * t2 / r);

   v1 = vec(vt.x + vtt.x, vt.y + vtt.y, vt.z + vtt.z);

   return 1;
}

int n;

point c[100];
double r[100];

point p;
vec v;

int main() {
   int T = 1;
   while (scanf("%d", &n), n != 0) {
   printf("Scene %dn", T++);
   for (int i = 0; i < n; i++) {
   scanf("%lf%lf%lf", &c[i].x, &c[i].y, &r[i]);
   c[i].z = 0;
   }
   point tp;
   scanf("%lf%lf%lf%lf", &p.x, &p.y, &v.x, &v.y);
   p.z = 0;
   v.z = 0;

   int cnt = 0;
   while (cnt <= 10) {
   double mind = 1e100;
   int mini = -1;
   point mp;
   vec mv;
   for (int i = 0; i < n; i++) {
   point np;
   vec nv;
   if (inter_circle_line(c[i], r[i], p, v, np, nv)) {
   double td = len2(vec(p, np));
   if (td < mind) {
   mind = td;
   mini = i;
   mp = np;
   mv = nv;
   }
   }
   }
   if (mini == -1) {
   printf("inf");
   break;
   }
   if (cnt == 10) {
   printf("...");
   break;
   }
   printf("%d ", mini + 1);
   cnt++;
   p = mp;
   v = mv;
   }
   printf("nn");
   }
   return 0;
}

标签云: 计算几何 | 日志分类: 计算几何

如需转载烦请添加以下信息: 本文转自Felicia@GCCFeli.cn

评论 (0)
引用通告 (0)
发表评论引用通告

暂无评论.

暂无引用通告

发表回复