SGU 110 – 3D Gemo Practice

/* Good Gemo Pratice
 * No reference, CSCI3260 is helpful!
 * Line - Sphere intersection
 * Point - Line distance
 * Reflection vector compute
 * WA on 6: should intersect the very first sphere in the direction
 */
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
 
#define EPS 1e-9
#define feq(x,y) (fabs((x)-(y))<EPS)
#define fgt(x,y) ((x)>((y)+EPS))
#define fge(x,y) (((x)+EPS)>(y))
 
typedef struct Pt{
    double x,y,z,r;
    Pt(){}
    Pt(double x,double y,double z):x(x),y(y),z(z){}
    Pt operator+(const Pt&p){ return Pt(x+p.x,y+p.y,z+p.z); }
    Pt operator-(const Pt&p){ return Pt(x-p.x,y-p.y,z-p.z); }
    Pt operator*(double k){ return Pt(x*k,y*k,z*k); }
    double operator*(const Pt&p){ return x*p.x+y*p.y+z*p.z; }
    Pt operator^(const Pt&p){ return Pt(y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x); }
    double len2(){ return x*x+y*y+z*z; }
    double len(){ return sqrt(len2()); }
    bool operator==(const Pt&p){ return feq(x,p.x)&&feq(y,p.y)&&feq(z,p.z); }
    void eat(){ scanf("%lf%lf%lf",&x,&y,&z); }
    void out(){ printf("(%lf,%lf,%lf)\n",x,y,z); }	
} Pt;
 
double P2Ldist(Pt P, Pt A, Pt B)
{
	return ((P-A)^(B-A)).len() / (B-A).len();
}
 
bool L2Sinter(Pt A, Pt B, Pt S, Pt& insP, double& tt)
{
	Pt norm = (B-A)*(1/(B-A).len());
	double a,b,c,delta,t1,t2;
	a=norm.len2();
	b=2*(norm*(A-S));
	c=(A-S).len2()-S.r*S.r;
	delta = b*b - 4*a*c;
	t1 = (-b-sqrt(delta))/(2*a);
	if(fgt(t1,0)) {
		insP = A+norm*t1;
		tt = t1;
		return true;
	} else return false;
}
 
Pt newDir(Pt V, Pt N)
{
	N = N*(1/N.len());
	return V-N*(2*(V*N));
}
 
int main()
{	
	int n,pos;
	double tt,mint;
	Pt s[55],A,B,insP,nDIR,recP;
 
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
	{
		s[i].eat();
		scanf("%lf",&s[i].r);
	}
	A.eat(); B.eat();
	bool never=false;
	int refcnt=0;
	while(refcnt <= 10)
	{		
		never=true;
		mint=1000000000;
		for(int i=1;i<=n;i++)
		{
			if(fge(s[i].r, P2Ldist(s[i],A,B)))
				if(L2Sinter(A,B,s[i],insP,tt))				 
				{	
					never=false;					
					if(fgt(mint, tt)) 
					{
						pos=i;	recP=insP;	mint=tt;
					}
				}
		}
		if(never) break; 
		else {
			refcnt++;
			if(refcnt>10) printf("etc.\n");
			else printf("%d ",pos);
			nDIR=newDir(B-A, recP-s[pos]);
			A=recP;
			B=A+nDIR;
		}
	}
	return 0;
}

Leave a Reply

Your email address will not be published.