A trick for dealing with periodic coordinates, not discussed in the article, is to convert each periodic coordinate (eg, an angle theta) into a pair of cartesian coordinates (x,y) on the circle, and then compute the distance in cartesian space. Ie, convert to x = cos(theta), y = sin(theta).<p>In the 2d example in the post, you would rescale the coordinates so the unit cell is from (-pi,pi) in both dimensions, and then the distance formula would be<p><pre><code> sqrt(((cos(x1) - cos(x2))**2 + ((sin(x1) - sin(x1))**2 + ((cos(y1) - cos(y2))**2 + ((sin(y1) - sin(y1))**2)
</code></pre>
This works well for doing things like determining the "nearest" point to another point, and similar operations, and has some other nice properties.<p>(I forget the name of this system, but it is commonly used for calculations involving dihedral angles in MD simulations)
I really enjoy this kind of post. It's short, simple and something most people could figure out by looking at the problem. Even then, it's still something that makes you go 'huh, how <i>would</i> you do that?'. Until the intuition about toroidal space kicks in, that is.
The may be a simpler way yet. If the coordinates are normalized from [0..1) and represented as 16bit integers with all bits fractional, you can simply take the difference x1-x2 and treat the result as signed. The result will be in the range [-0.5 .. 0.5). If you then square that number to get a 32bit result it will be correct for that component of the distance formula with no conditional code at all. Representation can be everything.
Isn't the topological manifold of the "wrap around" space described in the post spherical? A toroidal manifold would more complicated since the two generalized coordinates do not commute (order matters).
This article essentially describes a modulo (%) operator for floats, which is well-defined and completely analogous to the integer modulo operation. I always wondered why these aren't part of the FPU command set, or at least part of modern programming languages, e.g.:<p><pre><code> 1.0 % 5.0 = 1.0
6.0 % 5.0 = 1.0
7.0 % 5.0 = 2.0
7.4 % 5.0 = 2.4
7.4 % 1.0 = 0.4
</code></pre>
With that in place, the torodial distance becomes a one-liner without any case distinctions (i.e. without ifs).
This is one of my favorite things to do on the train while on the way to work. I select two random passing objects out the window, and use a general assumption about the train's speed. Then I visualize these objects on a toroidal space, and after a few minutes, I draw the rest of the owl.
This type of calculation is used on a very regular basis in molecular dynamics simulation. The wrap-around space is implemented using periodic boundary conditions. Calculating the distance between points is most often done using what is called the minimum image convention. In one dimension for a "box" of length L the minimum image distance between particles i and j: xij = xj - xi can be calculated through xij <- xij - L*nint(xij/L). This avoids any if statements.