#include #include #include #include using namespace std; numeric_limits real; double delta = sqrt(real.epsilon()), infinity = real.infinity(); struct Vec { // 3D vector double x, y, z; Vec(double x2, double y2, double z2) : x(x2), y(y2), z(z2) {} }; Vec operator+(const Vec &a, const Vec &b) { return Vec(a.x+b.x, a.y+b.y, a.z+b.z); } Vec operator-(const Vec &a, const Vec &b) { return Vec(a.x-b.x, a.y-b.y, a.z-b.z); } Vec operator*(double a, const Vec &b) { return Vec(a*b.x, a*b.y, a*b.z); } double dot(const Vec &a, const Vec &b) { return a.x*b.x + a.y*b.y + a.z*b.z; } double length(const Vec &a) { return sqrt(dot(a, a)); } Vec unitise(const Vec &a) { return (1 / sqrt(dot(a, a))) * a; } typedef pair Hit; struct Sphere; struct Scene { virtual ~Scene() {}; virtual void intersect(Hit &, const Vec &) const = 0; virtual bool intersect(const Vec &, const Vec &) const = 0; virtual Sphere bound(Sphere b) const = 0; }; struct Sphere : public Scene { Vec center; double radius; Sphere(Vec c, double r) : center(c), radius(r) {} ~Sphere() {} double ray_sphere(const Vec &dir) const { const Vec &v = center; double b = dot(v, dir), disc = b*b - dot(v, v) + radius * radius; if (disc < 0) return infinity; double d = sqrt(disc), t2 = b + d; if (t2 < 0) return infinity; double t1 = b - d; return (t1 > 0 ? t1 : t2); } bool sray_sphere(const Vec &orig, const Vec &dir) const { Vec v = center - orig; double b = dot(v, dir), disc = b*b - dot(v, v) + radius * radius; return (disc < 0 ? false : b + sqrt(disc) >= 0); } void intersect(Hit &hit, const Vec &dir) const { double l = ray_sphere(dir); if (l < hit.first) hit = Hit(l, unitise(l*dir - center)); } bool intersect(const Vec &orig, const Vec &dir) const { return sray_sphere(orig, dir); } Sphere bound(Sphere b) const { return Sphere(b.center, max(b.radius, length(center - b.center) + radius)); } }; typedef list Scenes; struct Group : public Scene { Sphere b; Scenes child; Group(Sphere b, Scenes c) : b(b), child(c) {} ~Group() { for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it) delete *it; } void intersect(Hit &hit, const Vec &dir) const { double l = b.ray_sphere(dir); if (l < hit.first) for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it) (*it)->intersect(hit, dir); } bool intersect(const Vec &orig, const Vec &dir) const { if (!b.sray_sphere(orig, dir)) return false; for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it) if ((*it)->intersect(orig, dir)) return true; return false; } Sphere bound(Sphere b) const { Sphere b2 = b; for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it) b2 = (*it)->bound(b2); return b2; } }; double ray_trace(const Vec &light, const Vec &dir, const Scene &s) { Hit hit(infinity, Vec(0, 0, 0)); s.intersect(hit, dir); if (hit.first == infinity) return 0; double g = dot(hit.second, light); if (g >= 0) return 0.; Vec sorig = hit.first*dir + delta*hit.second; return (s.intersect(sorig, -1. * light) ? 0 : -g); } Scene *create(int level, Vec c, double r) { Scene *s = new Sphere(c, r); if (level == 1) return s; Scenes child; child.push_back(s); double rn = 3*r/sqrt(12.); for (int dz=-1; dz<=1; dz+=2) for (int dx=-1; dx<=1; dx+=2) child.push_back(create(level-1, c + rn*Vec(dx, 1, dz), r/2)); Sphere b2(c + Vec(0, r, 0), r); for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it) b2 = (*it)->bound(b2); return new Group(b2, child); } int main(int argc, char *argv[]) { int level = (argc==3 ? atoi(argv[1]) : 9), n = (argc==3 ? atoi(argv[2]) : 512), ss = 4; Vec light = unitise(Vec(-1, -3, 2)); Scene *s(create(level, Vec(0, -1, 4), 1)); cout << "P5\n" << n << " " << n << "\n255\n"; for (int y=n-1; y>=0; --y) for (int x=0; x