delta = Sqrt[$MachineEpsilon];
RaySphere[o_, d_, c_, r_] :=
Block[{v, b, disc, t1, t2},
v = c - o;
b = v.d;
disc = Sqrt[b^2 - v.v + r^2];
t2 = b + disc;
If[Im[disc] != 0 || t2 <= 0, \[Infinity],
t1 = b - disc;
If[t1 > 0, t1, t2]]
]
Intersect[o_, d_][{lambda_, n_}, Sphere[c_, r_]] :=
Block[{lambda2 = RaySphere[o, d, c, r]},
If[lambda2 >= lambda, {lambda, n}, {lambda2,
Normalize[o + lambda2 d - c]}]
]
Intersect[o_, d_][{lambda_, n_}, Bound[c_, r_, s_]] :=
Block[{lambda2 = RaySphere[o, d, c, r]},
If[lambda2 >= lambda, {lambda, n},
Fold[Intersect[o, d], {lambda, n}, s]]
]
neglight = N@Normalize[{1, 3, -2}];
nohit = {\[Infinity], {0, 0, 0}};
RayTrace[o_, d_, scene_] :=
Block[{lambda, n, g, p},
{lambda, n} = Intersect[o, d][nohit, scene];
If[lambda == \[Infinity], 0,
g = n.neglight;
If[g <= 0, 0,
{lambda, n} =
Intersect[o + lambda d + delta n, neglight][nohit, scene];
If[lambda < \[Infinity], 0, g]]]
]
Create[level_, c_, r_] :=
Block[{obj = Sphere[c, r]},
If[level == 1, obj,
Block[{a = 3*r/Sqrt[12], Aux},
Aux[x1_, z1_] := Create[level - 1, c + {x1, a, z1}, 0.5 r];
Bound[c,
3 r, {obj, Aux[-a, -a], Aux[a, -a], Aux[-a, a], Aux[a, a]}]]]]
scene = Create[1, {0, -1, 4}, 1];
Main[level_, n_, ss_] :=
Block[{scene = Create[level, {0, -1, 4}, 1]},
Table[
Sum[
RayTrace[{0, 0, 0},
N@Normalize[{(x + s/ss/ss)/n - 1/2, (y + Mod[s, ss]/ss)/n - 1/2,
1}], scene], {s, 0, ss^2 - 1}]/ss^2, {y, 0, n - 1},
{x, 0, n - 1}]]
AbsoluteTiming[Export["image.pgm", Graphics@Raster@Main[9, 512, 4]]]
