Doubts about correctness of calculate_reliability! implementation in Unwrap

[wheeheee]

As mentioned in #576, copying here for convenience. It would be nice if any past contributors could shed some light on the algorithm, and even if there's nothing wrong a bit more documentation would be nicer.

For non-border regions, reliability is calculated over a cube around the pixel (3^N pixels, described in the comments as "nearest neighbours"), although there are only 2N Edges involving one pixel in populate_edges!, which are the nearest neighbours in the context of lattices.
But if there are circular dims / faces are connected, then pixel_shifts_border contains duplicate shifts. This makes more sense if only 2N neighbours should have been included in the reliability calculation, since then there would only be 1 shift per border pixel that wraps around.
Also, the reliability of pixels on the edges, including corners, are just rand(rng) and are not calculated.
I don't have access to the cited publications, so I can't check whether the current implementation is actually correct.