Raytracing diamonds

This comes from a question posed to me and others in 2019 by Stefan Langerman, a computational geometer who also happens to be a diamond dealer. The classic brilliant cut shape of diamonds has a flat “table” face on the top (visible) side of the gem, and many other smaller facets on the sides (“crown”) and back (“pavilion”), with angles chosen so that most of the light that enters through the table is reflected back out through the table, making the diamond look bright. Can we simulate this reflection process efficiently, without having to trace all the reflections of all (or a sufficiently dense sample of all) possible light rays?

This appears to be an interesting problem even for a very simple abstraction of the diamond’s shape: a two-dimensional isosceles triangle, with the base of the triangle as its table and the two equal sides as its pavilion. What you need to know about light reflections in gems is the critical angle, determined from the high refractive index of diamonds and other gems. Inside the gem, rays that hit a face at an angle shallower than the critical angle; rays that hit more steeply than the critical angle will escape. Any light ray hitting the table will always enter the gem, and will be refracted so that it enters along a ray steeper than the critical angle. It will then bounce around within the gem until it happens to hit a face steeply enough to escape. The question is to predict where this escape ray escapes, more quickly than the obvious algorithm of calculating explicitly the entire sequence of reflections that it performs. Or if we can’t, can we at least provide good bounds on the lengths of these reflection sequences?

Three of the simplest cases involve the isosceles triangles for which the triangle tiles the plane by copies reflected across its edges, the equilateral triangle, isosceles right triangle, and 30-30-120 isosceles triangle. The paths of a reflected light ray within such a triangle are combinatorially equivalent to the paths of a straight line through this triangulation: they pass through the same types of edges, at the same angles.

For all three of these triangle shapes, a line that enters the gem through any edge can only pass through a constant number of tiles in the tiling, before passing through another edge with the same angle, at which it must necessarily leave the gem (if it has not done so already). For instance, the red line shown, if it enters at its bottom left endpoint on the table of a green triangle, will exit after at most three reflections, through the parallel table of a blue triangle. The constant in the number of reflections here is absolute, not depending on the critical angle.

The next simplest case is the isosceles triangle with apex angle \(2\pi/5\). This triangle does not tile the plane, but we can still form a partial tiling by “unrolling” the path of reflections through the triangle into a line through a sequence of reflected triangles. In order to avoid exiting through one of the sides (by passing through it more perpendicularly than the entrance angle) and in order to miss the closest two copies of the table that are parallel to the entrance, there is only one passage that can be followed in the unrolled triangulation, the one shown by the red ray in the figure. Because this passage is perpendicular to one of the reflected sides of the triangle, a ray that is more perpendicular to the table than to this side can only continue for \(O(1)\) steps through this passage before hitting an exit. Here, though, the critical angle comes into play: critical angles that are closer to a right angle will allow larger (but constant) numbers of reflections.

Isosceles triangles with apex \(\pi/4\) are even messier. When a ray entering the table exits by a side, or exits by the table without first bouncing off the table, it always does so within $O(1)$ steps. The trickier cases are when a ray entering through the table misses the opposite table in the octagon formed by reflected copies of the triangle, and passes through a sequence of diagonal or horizontal reflections into other octagons before escaping. In the case that there are no horizontal reflections, only diagonal ones, we just need to divide the vertical distance that the ray needs to travel to hit the next row of vertical edges by the slope of the ray, and round. The case that there is at least one horizontal reflection (shown) is a little more compicated, but again appears to involve a finite calculation. So it seems that with a little more calculation it is still possible to determine algorithmically what happens in this case, in constant time rather than having to simulate the unbounded number of bounces triangle-by-triangle.

All of these triangles have angles that are rational multiples of \(2\pi\). For reflections within polygons of this type, the reflected line will only pass through a constant number of slopes. The rays of a single slope that hit a given line of the polygon can be describes as an interval, and the mapping on rays that takes each ray starting on one of the sides with one of these slopes to the start of its next reflection is a special kind of mapping studied in the theory of dynamical systems, an interval exchange transformation. So from a geometric question on ray-tracing, we are led to a more abstract question: given an interval exchange transformation, and a starting point in that transformation, how easy is it to determine the number of times one would need to repeat the transformation to return back to the same interval?

I have no real progress to report on any of these questions, in part because the continuous nature of these problems lead to difficult numerical issues. Instead, I have a new arXiv preprint “The Complexity of Iterated Reversible Computation” (arXiv:2112.11607 ) on related but more combinatorial problems. Suppose one has an invertible mapping \(f\) on a finite (but exponentially large) set of states, that is easy to compute, such as an interval exchange transformation on integers rather than real numbers. Suppose also that you want to repeat this mapping for a given initial state \(x\) and a large number \(n\) of steps. Is there anything you can do to compute the result \(f^{(n)}(x)\) of this iterated mapping, more efficiently than the obvious method of just computing \(x\), \(f(x)\), \(f(f(x))\), one step at a time? One property that appears highly relevant here is that these mappings are reversible: they are one-to-one, and more strongly the inverse mapping is also an interval exchange transformation. In contrast, it is easy to find mappings that are not reversible and for which computing iterated values is \(\mathsf{PSPACE}\)-complete: consider, for instance, the state transition map of a space-bounded Turing machine.

For some similar systems, we do know of effective methods in practice if not in theory. For Conway’s Game of Life, for instance, it appears to be \(\mathsf{PSPACE}\)-complete to compute what happens to an initial state after \(n\) steps (although I’m not sure this exact claim has been published in the literature), but hashlife (as implemented in Golly) can run many complicated and large patterns in an amount of runtime substantially less than the number of generated steps. The Game of Life is not reversible, but some other interesting cellular automata with similar behavior such as Critters are. Can we say anything about the computational complexity of simulating Critters or other reversible rules? My new paper addresses these questions by formulating a complexity class that includes computing \(f^{(n)}(x)\) for reversible cellular automaton simulation, integer interval exchange transformations, and others. It sandwiches this class between \(\mathsf{P}^{\mathsf{PH}}\) and \(\mathsf{PSPACE}\), and finds complete problems for it including the simulation of reversible cellular automata and of piecewise linear tranformations. The complete cellular automata include the two-dimensional billiard ball model, already known for its ability to simulate reversible circuits, and a new one-dimensional automata that can simulate some two-dimensional automata including the billiard ball model. Because Critters can in turn simulate the billiard ball model, it is also complete in the same sense. But I wasn’t able to prove completeness for integer interval exchange transformations, let alone say anything useful about the non-integer transformations needed to understand diamond reflections.

(Discuss on Mastodon or YC)

PS, added 2021-12-27 to appease the request on the YC discussion for an actual ray-traced diamond image: Here’s one from my paper “Isometric diamond subgraphs”, otherwise unrelated to this latest post:

Generated by POVray from a scene generated by this Python script, which all ran ok in 2008 when I wrote the paper (haven’t re-tested it recently).