In this computational essay, I investigated what happened when a cue ball bounced around the interior of a polygon for an arbitrary number of bounces. The goal was to discover patterns in the distribution of the coverings created by these processes. This data was then combined with image processing techniques to extract insights about these coverings and the effect of initial conditions on them.
Introduction
This project was motivated by a theorem from measure theory: given a bouncing ball in a square, we can cover the entire square in infinite time will eventually cover the entire square by starting the ball at an angle with an irrational tangent. I aimed to generalize this to more types of shapes. My process began with a programmatic method of generating polygons. Initial tests were run on random polygons, though as I progressed it became rhombi, as rhombi are relatively simple geometrically. In order to implement a bouncing ball, I decided to use a raycasting algorithm to simulate bouncing the ball iteratively hundreds or thousands of times. Once we have a plot of those bounces, we can use the image data to approximate the distribution of paths within the shape.
Raycasting is an algorithm that originates in 3D computer graphics. In the algorithm, we take a large number of rays and cast them at a scene, find where they intersect parts of that scene, and then generate an image based on that. For this project we have the sides of a polygon instead of a scene and a bouncing point instead of a pixel. In Chapter 2 we will examine what happens when this procedure is done in rhombi, squares, and randomly generated n-gons. The image data is then used to plot a probability distribution of the data based on the intensity of each pixel in our polygons.

Polygons

Initialization

The first task in making Polypool simulations is setting up data. For our rhombi, all we need is a length and two angles. For more complex shapes, we would need hyper-complex parameterizations, so we randomly generate them instead.
Take an angle and generate a unit rhombus represented as a Polygon:
In[]:=
genRhombus[angle_]:=CanonicalizePolygon[Parallelogram[{0,0},{{1,0},AngleVector[angle]}]]

Utility Functions

Raycasting is fundamentally a geometric algorithm. The following functions are all utilities that are necessary to make the main ray algorithm work.
sideList
takes a polygon and finds the list of
Line
objects that make up the edges of that polygon:
In[]:=
sideList[poly_]:=Line/@Partition[Partition[Flatten[List@@RegionBoundary[poly]],2],2,1];
Converts a point in a polygon to the edge of the polygon it resides on:
In[]:=
getCurrentSide[poly_,point_]:=First[SortBy[sideList[poly],RegionDistance[#,point]&]][[1]];
Convert a line as a pair of coordinate pairs to vector form:
In[]:=
lineToVecRep[line_]:=line[[2]]-line[[1]]
This function takes a vector and a pair of coordinate pairs (representing a line) and creates a new vector representing the path of a ball after hitting the side (See Acknowledgements):
In[]:=
calculateNewDirection[originVector_,sideVec_]:=Normalize[(2*Projection[originVector,sideVec])-originVector]
calculateIntersection
takes a ray and a line segment and computes the point where they intersect (this is the core of the raycast algorithm, see Acknowledgements):
In[]:=
calculateIntersection[ray_,segment_]:=​​Quiet@Solve[t{segment[[1]][[1]],segment[[1]][[2]]}+(1-t){segment[[2]][[1]],segment[[2]][[2]]}=={ray[[1]][[1]],ray[[1]][[2]]}+u{ray[[2]][[1]],ray[[2]][[2]]},{t,u},​​Assumptions->0<=t<=1&&u>0​​]
Pool and Pretty Pictures

The Main Raycasting Algorithm

The actual raycasting is implemented in a few functions, notably
castRayPolygon
and
castRayPolygonNestable
and their quadrilateral equivalents. These functions take an abstract representation of a ray and figure out where it will collide with the polygon again. This can be naively iterated over to compute many raycasts, though it is very slow, so
castRayPolygonNestable
allows the algorithm to be executed in accordance with the functional paradigm and Wolfram Kernel’s optimizations. This section contains two sets of raycasting algorithms. These are
castRayPolygon
and
castRayQuad
, and their corresponding nest wrappers. The reason for this is simply that the numerical algorithm used in
castRayQuad
did not work on general polygons, and vise versa. These functions synthesize all of the previous utilities and data to produce arbitrarily iterated raycasts.

Raycasting Implementations

The collision detection subroutine used by all the other raycasting algorithms:
In[]:=
detectCollision[ray_,sides_]:=Module[{possibleCollisions,dist,pointOfCollision,direction},​​direction=ray[[2]];​​possibleCollisions=calculateIntersection[ray,#]&/@(Sequence@@@sides);​​dist=Max[(Flatten[u/.possibleCollisions])/.u->-Infinity];​​ray[[1]]+dist*direction​​]
This function is fundamental to raycasting. It takes one ray and a list of sides and figures out what the point of intersection between that ray and those line segments is. Since ray intersections are the cornerstone of raycasting, this function ends up being very important. Next we have the three casting functions,
castRayQuad
,
castRayPolygon
, and
generateInitialCast
. These all use the same underlying raycasting concepts I’ve discussed throughout this whole essay, but are slightly different.
In[]:=
castRayQuad[quad_,point_,originVector_]:=Module[{currentSide,ray,sides,collisions,dist,newDirection,possibleCollisions,pointOfCollision},​​sides=sideList[quad];​​currentSide=getCurrentSide[quad,point];​​newDirection=calculateNewDirection[originVector,lineToVecRep[currentSide]];​​ray=HalfLine[point,newDirection];​​pointOfCollision=detectCollision[ray,sides];​​{pointOfCollision,point,newDirection}​​]
castRayQuad
is a function that takes a quadrilateral, a point, and a vector.
point
represents a point on one line of the polygon and
originVector
represents the direction of the ray that got us to that point.
In[]:=
castRayPolygon[poly_,point_,originVector_]:=Module[{currentSide,ray,sides,collisions,pointOfCollision,oppRay,newDirection,possibleCollisions},​​sides=sideList[poly];​​currentSide=getCurrentSide[poly,point];​​newDirection=calculateNewDirection[originVector,lineToVecRep[currentSide]];​​ray=HalfLine[point,newDirection];​​possibleCollisions=Apply[List,RegionIntersection[ray,RegionBoundary[poly]]][[1]];​​pointOfCollision=Last[SortBy[EuclideanDistance[#,point]&][possibleCollisions]];​​{pointOfCollision,currentSide,point,newDirection}​​]
castRayPolygon
is a slightly different version of the above function that works in general polygons. Otherwise, it works the same way.
Cast a ray at a given angle at a certain length along one side of a polygon:
In[]:=
generateInitialCast[quad_,angle_,length_]:=Module[{sides,ray,possibleCollisions,dist,pointOfCollision,direction,startPoint},​​sides=sideList[quad];​​startPoint={sides[[1,1,1,1]]+length,0};​​direction=AngleVector[angle];​​ray=HalfLine[startPoint,direction];​​pointOfCollision=detectCollision[ray,sides];​​{pointOfCollision,pointOfCollision-startPoint}​​]
This function ties together all the initial conditions we set in Part One to start the raycasting process. It accompanies
castRayQuad
to allow us to set specific initial conditions for our raycasts. It returns the point along the polygon that would be the first collision given the initial conditions.

castRayPolygonNestable

This function, while important, is much simpler than
castRayPolygon
because it just abstracts the logic away from the previous function. It works with either the polygon or quadrilateral versions.
Allows the raycasting algorithm to be
Nest
ed:
In[]:=
castRayPolygonNestable[poly_,{point_,originVector_}]:=​​ castRayPolygon[poly,point,originVector][[{1,4}]]
In[]:=
castRayQuadNestable[poly_,{point_,originVector_}]:=​​ castRayQuad[poly,point,originVector][[{1,3}]]

And Finally...

We now have the utility to actually begin bouncing some cueballs. In this example, we will use a rhombus. Rhombi are interesting shapes to use as pool tables because while they are geometrically and algebraically simple, they create interesting lattice pattern that make them good test shapes. The following block of code-text walks through the raycast algorithm on a rhombus. Later on, this will be modularized and expanded into a succinct function block.
Initializes an empty list of points:
In[]:=
points={};
Creates data to represent our initial conditions:
In[]:=
quad=genRhombus[Pi/4];​​initialAngle=Pi/6;​​initialLength=0.5;
This takes our rhombus and a choice of initial angle and length along the rhombus edge:
In[]:=
{breakPoint,startDirection}=generateInitialCast[quad,initialAngle,initialLength];
Plots first raycast and sets up the points list:
In[]:=
castInfo=castRayQuad[quad,breakPoint,startDirection];​​AppendTo[points,castInfo[[1]]];
Iterates the raycast 1,000 times:
In[]:=
Nest[(c=castRayQuadNestable[quad,#];AppendTo[points,c[[1]]];c)&,{breakPoint,startDirection},1000];
Splits the points into overlapping pairs representing line segments:
In[]:=
individualLines=Line/@Partition[points,2,1];
Displays the resulting lines as an image:
In[]:=
rayPainting=Style[Graphics[{Opacity[0.05],individualLines}],Antialiasing->False]
Out[]=
As you can see, rhombi are simultaneously very interesting and very uninteresting at the same time. This specific rhombus has alternating areas of high density in the corners and edges and much lower density in the very middle and in a cross shape around it. We’ll analyze this more rigorously in the following section.

Image Processing 101

As you can see, we can get interesting patterns from even the humble rhombus. We will now quantify what interesting means through some convenient image processing. This manifests as a grayscale histogram. This histogram is created by converting to grayscale data and compensating for opacity.
​
Compensation for opacity is necessary because of the way we drew images. To draw pixels that get darker where more rays intersect, we use
Opacity
. However, in image processing, opacity is a mulitplicative value, not additive. This means that every time a ray crosses another ray, their color values are scaled, not added together. To undo this exponential effect, we take the logarithm. Once we have all that data, it can be very easily plotted onto a histogram as below.
Rasterizes the raycast image to get pixel data:
In[]:=
raster=Rasterize[rayPainting];
Takes the logarithm and negates the grayscale values of the raster:
Plots the resulting data as a histogram:

General N-gons

Modularization

We will now “show off” our new raycast generator by computing raycasts in a very thin rhombus:
Compute raycast data in a very thin rhombus:
Renders the result of rhombus casting:
Along with the polygon itself, we can also analyze the and histograms. With rhombi, the distributions are always unimodal. Notably, thin rhombi have more spread distributions compared to rhombi with larger interior angles.
Squares are a good baseline, as they almost always have the same lattice-like structure. The distributions of squares are almost always roughly uniform as well because rays bouncing inside squares generally don’t intersect themselves repeatedly.
As discussed earlier, we have a version of the full cast algorithm for general polygons with random conditions:
We will start off with a regular hexagon just to give an example, and then the code for a random polygon will be provided. As with rhombi, we can display the same sort of data. Notice that polygons with random initial conditions have a much greater variety in distributions and transforms, though regular polygons will give more uniform patterns.
Render hexagon cast:
The histogram distribution of regular polygons are roughly unimodal, though there are example that don’t display this behavior exactly. They are also usually symmetrical.
As you can see, this randomly generated septagon has a much different distribution than the the rhombus or square example. This generally holds for most random polygons, with a lot of them even having areas of almost zero density after several thousand bounces. Randomized polygons almost always have much more spread out intensity distributions than any of the other classes I examined.
Generate a histogram of pixel intensity:
Conclusions
As we can see, this algorithm can produce a huge variety of results in a wide array of possible polygons. Given even more computing power and higher resolution images, more complete data sets could be derived . Of course, could be explored more rigorously using pure mathematics and analysis, though that is out of the scope of this essay. Further directions include more optimizations of this algorithm and a more unified codebase (less fragmentation between classes of polygons) and better analysis of raycast data. I also would have loved to quantitatively analyze the histograms created by these distributions, though there wasn’t enough time to derive a method for doing so. One other possibility I discussed with a groupmate was employing Fourier analysis on the resulting images, but creating useful interpretations of that information was out of the scope of this project.

Acknowledgements

I’d like to thank my wonderful mentor Dr. Lyman Hurd for all his assistance with mathematics, program implementations, code reviews, and putting up with my insanity for 7 hours a day. I’d also like to thank program faculty Dugan Hammock for assistance with writing the intersection algorithms and fellow program student Christopher Gilbert for mathematical assistance with the ray direction calculations. I’d also like to thank my group mate Max Niederman and program directors Eryn Gillam and Megan Davis for being sources of endless helpful distractions during many a stressful coding session and director Rory Foulger for teaching me the Wolfram Language. Finally, I’d like to thank Dr. Stephen Wolfram for giving me this project and for telling me about the irrational angle covering theorem in squares.

References

◼
  • 1 . William, Simon Woods, & Matthias Odisio. (2013, June 1). Getting RGBA data from an image. Mathematica Stack Exchange. https://mathematica.stackexchange.com/questions/26770/getting-rgba-data-from-an-image
    • ◼
  • 2 . Points, lines, and planes. Point, Line, Plane. (n.d.). https://paulbourke.net/geometry/pointlineplane/
    • ◼
  • 3. Users of Mathematica Stack Exchange. (1959, February 1). Calculate the 2d Fourier transform of an image. Mathematica Stack Exchange. https://mathematica.stackexchange.com/questions/29203/calculate-the-2d-fourier-transform-of-an-image

    • CITE THIS NOTEBOOK

    Polypool: explorations into infinitely bouncing cueballs in geometric shapes​
    by Bertie Bennett​
    Wolfram Community, STAFF PICKS, July 12, 2024
    ​https://community.wolfram.com/groups/-/m/t/3215621