【SGU】110. Dungeon 计算几何——线对球的反射射线

传送门：【SGU】110. Dungeon

题目分析：本题重点就在于对给定入射射线AB计算反射射线BC。

对于空间上的入射射线的起点A，射线与球的交点B，球的圆心C。

因为向量AB总可以表示成|AB|*e，其中e为向量AB的方向向量，|AB|为向量AB的模，所以|AO - AB| = |BO| = R

将A，O，R带进去可以得到一个关于|AB|的一元二次方程，解方程求得x1，x2，取其中正的较小的一个。（可能存在无解，莫忘判断）

得到|AB|后，AB = e*|AB|，同时B = AB + Ａ，BO = O - B，此时BC = AB - AB·BO / |BO| * 2 * ( BO / |BO| )，做一个图就可以发现这个公式。

我们就是利用上面的公式，然后将求得的满足|AB|长度最小的交点B，以及相交的球的球心O带入来求得反射射线BC（新的入射射线AB），反射射线的起点B（新的入射射线的起点A）。

需要注意的是，可能起点在一个球上，这时候我们需要判断使得一条射线不能连续在同一个球上反射。

代码如下：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;

typedef long long LL ;

#define rep( i , a , b ) for ( int i = ( a ) ; i <  ( b ) ; ++ i )
#define For( i , a , b ) for ( int i = ( a ) ; i <= ( b ) ; ++ i )
#define rev( i , a , b ) for ( int i = ( a ) ; i >= ( b ) ; -- i )
#define travel( e , H , u ) for ( Edge* e = H[u] ; e ; e = e -> next )
#define clr( a , x ) memset ( a , x , sizeof a )

const int MAXN = 55 ;
const double INF = 1e10 ;
const double eps = 1e-5 ;

struct Point {
	double x , y , z ;
	Point () {}
	Point ( double x , double y , double z ) : x ( x ) , y ( y ) , z ( z ) {}
	Point operator + ( const Point& p ) const {
		return Point ( x + p.x , y + p.y , z + p.z ) ;
	}
	Point operator - ( const Point& p ) const {
		return Point ( x - p.x , y - p.y , z - p.z ) ;
	}
	Point operator * ( const double& b ) const {
		return Point ( x * b , y * b , z * b ) ;
	}
	Point operator / ( const double& b ) const {
		return Point ( x / b , y / b , z / b ) ;
	}
	void input () {
		scanf ( "%lf%lf%lf" , &x , &y , &z ) ;
	}
	void show () {
		printf ( "( %.5f , %.5f , %.5f )\n" , x , y , z ) ;
	}
} ;

typedef Point Vector ;

struct Sphere {
	Point p ;
	double r ;
	void input () {
		p.input () ;
		scanf ( "%lf" , &r ) ;
	}
} ;

Sphere ball[MAXN] ;
int n ;

int dcmp ( double x ) {
	return ( x > eps ) - ( x < -eps ) ;
}

double Length ( const Vector& p ) {
	return sqrt ( p.x * p.x + p.y * p.y + p.z * p.z ) ;
}

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

void solve () {
	Point pos ;
	Vector dir ;
	For ( i , 1 , n ) ball[i].input () ;
	pos.input () ;
	dir.input () ;
	dir = dir - pos ;
	int now = 0 ;
	For ( i , 1 , 11 ) {
		int idx = 0 ;
		double len = INF ;
		Point A , B , O ;
		A = pos ;
		For ( j , 1 , n ) {
			if ( now == j ) continue ;
			Vector AO = ball[j].p - A ;
			double a = dir.x * dir.x + dir.y * dir.y + dir.z * dir.z ;
			double b = - 2 * ( AO.x * dir.x + AO.y * dir.y + AO.z * dir.z ) ;
			double c = AO.x * AO.x + AO.y * AO.y + AO.z * AO.z - ball[j].r * ball[j].r ;
			double delta = b * b - 4 * a * c ;
			if ( dcmp ( delta ) < 0 ) continue ;
			double x1 = ( - b - sqrt ( delta ) ) / ( 2 * a ) ;
			double x2 = ( - b + sqrt ( delta ) ) / ( 2 * a ) ;
			//printf ( "%.5f %.5f\n" , x1 , x2 ) ;
			double x = INF ;
			if ( dcmp ( x1 ) > 0 ) x = min ( x , x1 ) ;
			if ( dcmp ( x2 ) > 0 ) x = min ( x , x2 ) ;
			if ( dcmp ( x - INF ) == 0 ) continue ;
			else if ( dcmp ( len - x ) > 0 ) {
				O = ball[j].p ;
				idx = j ;
				len = x ;
			}
		}
		if ( !idx ) break ;
		if ( i == 11 ) {
			printf ( " etc." ) ;
			break ;
		}
		B = dir * len + pos ;
		Vector BO = O - B ;
		//dir.show () ;
		dir = dir - BO * ( Dot ( dir , BO ) / pow ( Length ( BO ) , 2.0 ) * 2.0 ) ;
		//dir.show () ;
		pos = B ;
		if ( i == 1 ) printf ( "%d" , idx ) ;
		else printf ( " %d" , idx ) ;
		now = idx ;
	}
	printf ( "\n" ) ;
}

int main () {
	while ( ~scanf ( "%d" , &n ) ) solve () ;
	return 0 ;
}