[Shootout-list] Ray tracer developments

[Jon Harrop]

--Boundary-00=_7VQdCIbHtM9e9H1
Content-Type: text/plain;
  charset="iso-8859-1"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

On Sunday 01 May 2005 15:50, Jon Harrop wrote:
> On Sunday 01 May 2005 15:39, Sebastien Loisel wrote:
> > > The inheritance hierarchy is certainly one of the most verbose parts.
> > > But what would you use instead? You could use a tagged union, but that
> > > is much more common in C than in C++.
> >
> > I've got a lot of work so I haven't looked in detail, but I suspect
> > you can replace everything by a single "node" struct, with a list of
> > children. If the list is empty, it's a sphere, otherwise it's a group.
>
> An excellent idea! ...

The new C++ version is about 10 lines shorter, weighing in at 85 LOC, and ~10% 
slower on x86.

-- 
Dr Jon D Harrop, Flying Frog Consultancy Ltd.
Objective CAML for Scientists
http://www.ffconsultancy.com/products/ocaml_for_scientists

--Boundary-00=_7VQdCIbHtM9e9H1
Content-Type: text/x-c++src;
  charset="iso-8859-1";
  name="ray3.cpp"
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment;
	filename="ray3.cpp"

#include <vector>
#include <iostream>
#include <limits>
#include <cmath>
using namespace std;
numeric_limits<double> dbl;
double delta = sqrt(dbl.epsilon()), infinity = dbl.infinity(), pi = M_PI;
struct Vec { // 3D vector
  double x, y, z;
  Vec(double x2, double y2, double z2) : x(x2), y(y2), z(z2) {}
};
Vec operator+(Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y, a.z + b.z); }
Vec operator-(Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y, a.z - b.z); }
Vec operator*(double a, Vec b) { return Vec(a * b.x, a * b.y, a * b.z); }
double dot(Vec a, Vec b) { return a.x*b.x + a.y*b.y + a.z*b.z; }
Vec unitise(Vec a) { return (1 / sqrt(dot(a, a))) * a; }
struct Ray { Vec orig, dir; Ray(Vec o, Vec d) : orig(o), dir(d) {} };
struct Tree {
  vector<Tree> child;
  Vec center;
  double radius;
  Tree(Vec c, double r) : child(), center(c), radius(r) {}
};
double ray_sphere(const Ray ray, const Tree &tree) {
  Vec v = tree.center - ray.orig;
  double b = dot(v, ray.dir),
    disc = b*b - dot(v, v) + tree.radius * tree.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);
}
void intersect(double &lambda, Vec &normal, const Ray ray, const Tree &tree) {
  double l = ray_sphere(ray, tree);
  if (l >= lambda) return;
  if (tree.child.size() == 0) {
    lambda = l;
    normal = unitise(ray.orig + l * ray.dir - tree.center);
  } else
    for (vector<Tree>::const_iterator it=tree.child.begin();
	 it!=tree.child.end(); ++it)
      intersect(lambda, normal, ray, *it);
}
double ray_trace(const double weight, const Vec light, const Ray ray,
		 const Tree &tree) {
  double lambda = infinity;
  Vec normal(0, 0, 0);
  intersect(lambda, normal, ray, tree);
  if (lambda == infinity) return 0;
  Vec o = ray.orig + lambda * ray.dir + delta * normal;
  double g = -dot(normal, light), l = infinity;
  if (g <= 0) return 0.;
  intersect(l, normal, Ray(o, Vec(0, 0, 0) - light), tree);
  return (l == infinity ? g : 0);
}
Tree create(int level, double r, double x, double y, double z) {
  if (level == 1) return Tree(Vec(x, y, z), r);
  Tree group = Tree(Vec(x, y, z), 3*r);
  group.child.push_back(Tree(Vec(x, y, z), r));
  double rn = 3*r/sqrt(12.);
  for (int dz=-1; dz<=1; dz+=2)
    for (int dx=-1; dx<=1; dx+=2)
      group.child.push_back(create(level-1, r/2, x-dx*rn, y+rn, z-dz*rn));
  return group;
}
int main(int argc, char *argv[]) {
  int level = (argc==2 ? atoi(argv[1]) : 6), w = 512, h = 512, ss = 4;
  Tree tree = create(level, 1, 0, -1, 0);
  cout << "P2\n" << w << " " << h << "\n255\n";
  for (int y=h-1; y>=0; --y) {
    for (int x=0; x<w; ++x) {
      double g=0;
      for (int dx=0; dx<ss; ++dx)
	for (int dy=0; dy<ss; ++dy) {
	  Vec d(x+double(dx)/ss-w/2, y+double(dy)/ss-h/2, max(w, h));
	  g += ray_trace(1, unitise(Vec(-1, -3, 2)), Ray(Vec(0, 0, -4),
							 unitise(d)), tree);
	}
      cout << min(255, int(256. * g / (ss*ss))) << " ";
    }
    cout << endl;
  }
  return 0;
}

--Boundary-00=_7VQdCIbHtM9e9H1--