Calculating the Distance Between Points in Wrap Around (Toroidal) Space | Hacker News Reader

Calculating the Distance Between Points in Wrap Around (Toroidal) Space

blog.demofox.org

106 pointsAtrix25639 comments

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).

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

    sqrt(((cos(x1) - cos(x2))**2 + ((sin(x1) - sin(x1))**2 + ((cos(y1) - cos(y2))**2 + ((sin(y1) - sin(y1))**2)

This works well for doing things like determining the "nearest" point to another point, and similar operations, and has some other nice properties.

(I forget the name of this system, but it is commonly used for calculations involving dihedral angles in MD simulations)

It depends on how you define "distance". My initial thought when reading your comment was that you couldn't possibly use the Euclidean distance because you'd be drawing a path between two points on the circle that involves points _not_ on the circle.

Using the Euclidean distance (Cartesian distance? I dunno) works because it's a monotonic underestimate of the distance along the circle. So you know if cartesian_dist(A, B) > cartesian_dist(A, C) then circle_dist(A, B) > circle_dist(A, C).

The original problem, though, is about computing distance, not about comparing distances. Comparing distances is just used as an intermediate step in the initial method given.

For anyone curious, I’m guessing “MD simulation” means https://en.wikipedia.org/wiki/Molecular_dynamics

Folks who use this kind of distance measure might note that it also works well as a metric between angles (i.e. rotations) or between points on the projective line. For instance, see https://www.youtube.com/watch?v=oJAn--vsAzc (unfortunately not too self-contained a source)

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 would 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.

x1 = 0.9

x2 = 0

x1-x2 = 0.9 => -0.4

I don't see how this works.

Add 0.5 to the result to convert the result in [-0.5 .. 0.5) back to [0..1)

    -0.4 + 0.5 => 0.1

No, don't subtract 0.5. Just reinterpret (through a cast) the MSB as a sign bit: 0.9 => -0.1

You could also possibly use BAMS and ignore a negative space. You'd also likely save a few bytes on your location data as well (unsigned short instead of float).

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).

When wrapping rectangular sections, you create a torus. Spherical wraparound requires trapezoids.

I think you're mistaken. Torii come from non-rectangular regions, too, though they'll always be parallelograms. Spheres cannot arise like this from any shape.

Mathematically, you'd say that a sphere is simply connected and hence is its own universal cover and hence is not the quotient of a plane. Meanwhile, the universal cover of a torus is a plane and the fundamental domain of that covering can always be taken to be a parallelogram. The study of all possible such parallelograms is the standard example of a Teichmuller space [1] which is well-known recently due to both the use of the idea in "inter-universal Teichmuller theory" developed for the famous ABC conjecture work in 2012 as well as for being the subject of study of a winner of the most recent Fields medal [2] who is especially notable for being the first female recipient and was in the news again recently sadly for having died of cancer at a young age.

[1] https://en.wikipedia.org/wiki/Teichm%C3%BCller_space [2] https://en.wikipedia.org/wiki/Maryam_Mirzakhani

I'm really confused.

Sperical wraparound can be viewed as a disk with it's perimeter identified to a point.

Maybe I'm just missing something obvious in the discussion here? Also, how is Teichmüller Theory relevant? Is there some deep insight I'm just missing?

Yes, that works. But it loses the geometry of the plane when you do that. Or in terms of topology it's not a covering space of the sphere. When you represent the torus as the plane cut up into squares, any small piece of the torus looks just like a small piece of the plane. But with your example, the North Pole (say) comes from this extra point you add in. So it's a topological quotient of the disk plus a point, but that's not a covering space. So it's much less useful, particularly if you want to compute distances. Distances (and angles and other geometric properties) behave well under coverings but not under arbitrary topological quotients like you have. In particular, there's no flat sphere so you can't possibly have it covered by the flat plane.

Teichmuller theory isn't actually important here, just a connection to a piece of mathematics that's been in the news.

Thanks.

The discussion just says "torus" leaving everyone to infer what category we are working in. It seems like people are inferring different things, and that's what confused me.

I'm not super familiar with homotopy theory but do know that universal covers pull back a lot of nice structure from the base space. Is there a theorem that makes this precise?

Anyway, since the original article here discusses distances, I'd have thought that we'd be wanting to hold ourselves to isometries or something, which universal covers usually aren't, e.g. R^2/Z^2 isn't isometric to a doughnut.

Of course, you're right that what mathematical structures are the equivalent of the everyday terms we're talking about here depends on exactly what we're trying to accomplish. (Though note that the definition of trapezoid versus a square inherently makes use of lengths or angles hence implies we need to be working on geometric and not purely topological structures to be able to make any statement about trapezoids that are meaningful. I.e. a trapezoid can be the domain of a covering of a torus but only when we are interested only in the topology and not the geometry and hence only when "trapezoid" is indistinguishable from "square.")

Yeah coverings can't be isometries (unless they're trivial), but they can always be local isometries which is what we'd most likely want here. R^2/Z^2 is (isometric to) a flat torus, which is, I'd argue, the nicest (though not the tastiest) torus. (Though the question of which flat torus is nicest is another matter and many would say it's not the torus that comes from the square!)

I don't know of any general statement for you about what pulls back by coverings but in general anything you can describe as "local" pulls back. That's pretty vague though, sorry!

After some thought (and before I read your comment), I am indeed mistaken. Take the long/lat mapping we put onto Earth - it contains trapezoids.

However, the sections on a torus are not always parallelograms. Take your donut and lay it on the table. It has an equator of sorts where the oil didn't cook it as dark. Imagine little rectangles around the donut, long sides resting on this outer equator, short sides touching adjacent rectangles. We continue stacking rectangles over the donut and toward the inner equator. If we maintain the same number of rectangles circling around the central axis (the one that passes through the center, perpendicular to the equatorial plane, but never intersects the donut), then the long sides of the rectangles must get shorter as we approach the inner equator. If we look at one of these rectangles, we can see that the long side toward the outer equator is longer than the side toward the inner equator. And that the short sides are no longer parallel. Thus, these sections are not parallelograms.

You're right that if you take a trapezoid and identify the sides, you get something that is topologically a torus. However, we were trying to measure distances on the torus and hence need it to be geometrically a torus - i.e. we need distances to preserved at least locally when we map from the domain (the square or trapezoid) to the torus. You can see that this isn't done with the trapezoid: you have to identify two edges that have different lengths, hence lengths must get messed up! Or in your example, each trapezoid on the donut has different size and shape so it's not useful for measuring distances. Though getting this to work for the donut is impossible, so you had a difficult task. (If your domain is in the flat plane, then the resulting torus must be flat as well. Which means it's not the kind of torus you can eat since it can't fit in 3 dimensional (flat) space along side us. But Pac Man and friends live on flat torii.)

Nonetheless, you're right that the domains don't always have to be parallelograms. Each covering has a parallelogram domain but in general coverings have an infinite number of possible domains. They satisfy some conditions though: essentially they're a parallelogram where the edges don't have to be straight lines so long as the left and the right edges (and the top and bottom edges) are the same curve just translated. So generally the parallelogram domain is easiest to work with (and in fact there are multiple parallelogram domains you can choose so one usually chooses the one with the smallest edge lengths).

That manifold is a torus - you can walk in circles by walking West-East or North-South, each of which repeats every 1 distance unit. If you walk along the diagonal you get a repeating path every 1.4 distance units. All "great circle" paths (straight lines on a manifold) on a sphere repeat every 1 unit of distance.

I'm not sure what you mean, could you elaborate?

The post says that for a torus, which is a product space of two (or more?) circles, we can analyze the the distance "along each circle" from a shared origin for a chosen basis. One circle is our "x" and the other is our "y". (Contrast with a plane and two lines.)

This seems reasonable if the objects we're looking to measure our distance to are expressed in terms of those circles as well, since this allows us to simply put a measure on each circle and get a distance formula for any two points that behaves the way we would like.

I'm out of my depth but I believe it's a torus because it's 2 cylinders, not circles?

A cylinder is a circle and a line. The circle goes around and the line is the height.

A torus is a circle and a circle. One circle goes from the outside to the inside and back around, and the other goes around the outside, if we're talking about a donut. We can see that any point on the surface of the donut is some combination of "go around the outside" and "walk around towards the center". You can also imagine taking a cylinder and bending it around to make a torus -- this changes the 'height' line into a second circle.

We're talking about the product space, so what that means is you can describe the space by (two) coordinates, one drawn from each "shape" (actually, manifold). A "torus" then differs from a "plane" in that you can wrap around in the x-direction and wrap around in the y-direction, independently, since the x-direction on a torus is a circle instead of a line and the y-direction on a torus is a circle instead of a line. By contrast, a cylinder you can only wrap around one way, because you can only wrap in the direction that has a circle for its coordinate space.

You cannot take a piece of paper (of any size) and wrap it around part of a donut.

I'm not sure what you're trying to say, could you elaborate?

you're thinking of all boundaries mapping to a point.

in your interpretation, the red circle moving to the boundary would look like this https://cdn.pixabay.com/photo/2013/03/08/05/53/frame-91452_9...

Instead, look at how the arrow types match here: https://upload.wikimedia.org/wikipedia/commons/3/37/Projecti...

the side you choose to glue first matters (compare before/after) https://upload.wikimedia.org/wikipedia/commons/b/ba/Inside-o...

and the orientation of the arrows https://upload.wikimedia.org/wikipedia/commons/6/64/Klein_Bo... also matters https://upload.wikimedia.org/wikipedia/commons/5/5c/Klein_bo...

This example is not the surface of a real donut. It's an abstract surface.

In the surface of a real donut you have a lot of problems due to the curvature and the different sizes, so the calculation is very complicated. For example, in a real donut the circle at the top and at the bottom of the donut have different length that the two circles at the equator. In the torus in this page all the horizontal circles have the same length.

This is a magical torus, where all the horizontal and vertical circles have the same length. It also has no Gaussian curvature. You can't construct one of this magical torus as a surface in the 3D space, but you can construct one as a surface in the 4D space.

A flat torus can be isometrically embedded in 3D space with a continuously differentiable embedding. Visualising this embedding using computer graphics was the aim of a French project, completed in 2012. Turned out a flat torus in 3D has cool looking fractal structure. http://www2.cnrs.fr/en/2027.htm

Nice result!

But, is this construction injective? I.E. do the little waves not produce any fake crossings in the surface?

IIRC from my Differential Geometry classes, immersions and embeddings are more general that I expected, and they include some strange cases. I don't remember the details, but looking at the Wikipedia page they include the embedding of the Klein bottle in R^3.

It's immersions that include some edge cases. Embeddings effectively identify a manifold with a submanifold of another manifold.

Klein bottles can't be embedded into R^3.

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.:

   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

With that in place, the torodial distance becomes a one-liner without any case distinctions (i.e. without ifs).

> I always wondered why these aren't part of the FPU command set, or at least part of modern programming languages

They are a part of modern (and not-so-modern) programming languages (% in Python, Rust, etc.). Even C from quarter of a century ago supports them (fmod). There are also at least two CPU instructions on x86, fprem and fprem1.

I couldn't even find an integer modulus operation when I was porting my fractal code to WebGL. https://www.shadertoy.com/view/MdSXWc

Did the % operator not work? I just tested 1 % 10 and it seemed to compile. I used the % operator in other shaders. Remember that WebGL is ultimately a sub and superset of C.

It does work! I wonder why I didn't think of it while transliterating from my older Python codebase.

It's possible it used to not work with an older version of WebGL/GLSL, and now it does.

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.

Reply on news.ycombinator.com