• A Mathematica Optimization Problem
  http://xahlee.org/UnixResource_dir/writ/Mathematica_optimization.html

norm = Function[#/Sqrt@(Plus @@ (#^2))];

delta = Sqrt[$MachineEpsilon];

myInfinity = 100000. ;

Raytrace code

Clear[RaySphere];
RaySphere = 
        Compile[{o1, o2, o3, d1, d2, d3, c1, c2, c3, r}, 
          Block[{v = {c1 - o1, c2 - o2, c3 - o3}, 
              b = d1*(c1 - o1) + d2*(c2 - o2) + d3*(c3 - o3), 
              discriminant = -(c1 - o1)^2 - (c2 - 
                        o2)^2 + (d1*(c1 - o1) + d2*(c2 - o2) + 
                        d3*(c3 - o3))^2 - (c3 - o3)^2 + r^2, disc, 
              t1, t2}, 
            If[discriminant < 0. , myInfinity, disc = Sqrt[discriminant]; 
              If[(t1 = b - disc) > 0. , t1, 
                If[(t2 = b + disc) ≤ 0. , myInfinity, 
                  t2]]]]];

Remove[Intersect];
Intersect[{o1_, o2_, o3_}, {d1_, d2_, d3_}][{lambda_, n_}, 
        Sphere[{c1_, c2_, c3_}, r_]] := 
      Block[{lambda2 = RaySphere[o1, o2, o3, d1, d2, d3, c1, c2, c3, r]}, 
        If[lambda2 ≥ lambda, {lambda, n}, {lambda2, 
            norm[{o1, o2, o3} + 
                lambda2 *{d1, d2, d3} - {c1, c2, c3}]}]]

Intersect[{o1_, o2_, o3_}, {d1_, d2_, d3_}][{lambda_, n_}, 
        Bound[{c1_, c2_, c3_}, r_, s_]] := 
      Block[{lambda2 = RaySphere[o1, o2, o3, d1, d2, d3, c1, c2, c3, r]}, 
        If[lambda2 ≥ lambda, {lambda, n}, 
          Fold[Intersect[{o1, o2, o3}, {d1, d2, d3}], {lambda, n}, 
            s]]]

Clear[neglight, nohit]
neglight = N@norm[{1, 3, -2}];
nohit = {myInfinity, {0. , 0. , 0. }};

Clear[RayTrace];
RayTrace[o_, d_, scene_] := 
      Block[{lambda, n, g, 
          p}, {lambda, n} = Intersect[o, d][nohit, 
            scene]; 
If[lambda == myInfinity, 0, 
          g = n . neglight; 
If[g ≤ 0, 
            0, {lambda, 
                n} = Intersect[o + lambda d + delta n, neglight][nohit, 
                scene]; 
If[lambda < myInfinity, 0, 
              g]]]]

Clear[Create];
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]}]
]
]]

Main[level_, n_, ss_] := 
      With[{scene = Create[level, {0. , -1. , 4. }, 1. ]}, 
        Table[Sum[
              RayTrace[{0, 0, 0}, 
                N@norm[{(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}]]

\(\(\(Timing[
      Export["image.pgm", Graphics@\(Raster@Main[3, 200, N@3]\)]]\)\(
\)
    \)\)