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)