Incorrect Delauny triangulation returned for 2D boxes

[klywang]

Describe the bug Voronoi diagrams calculated in 2D boxes return the expected results, however, the Delauny triangulation returned is incorrect. Expected results can be obatined using scipy.spatial.Delauny or with 3D structural equivelants.

To Reproduce Using the following function:

import freud
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
import gsd.hoomd
import numpy as np

def voro(filename, is2D=False):
    frame = gsd.hoomd.open(filename, mode='rb')[0]
    pos = np.copy(frame.particles.position)

    # Crop
    pos = pos[abs(pos[:, 0]) < 20]
    pos = pos[abs(pos[:, 1]) < 20]
    box = freud.Box(Lx=40, Ly=40, Lz=frame.configuration.box[2])

    # Calculate Voronoi diagram
    voro = freud.locality.Voronoi().compute((box,  pos))
    if is2D is False:
        # Removes bonds between layers
        voro.nlist.filter(voro.nlist.distances > 0.1)

    # Plot
    fig = plt.figure(figsize=(20, 20))
    ax = fig.gca()
    if is2D:
        voro.plot(ax=ax)
    line_data = np.asarray([[pos[i],
                             pos[i] + box.wrap(pos[j] - pos[i])]
                            for i, j in voro.nlist])[:, :, :2]
    line_collection = LineCollection(line_data, colors='black', alpha=0.2)
    ax.add_collection(line_collection)
    ax.set_ylim([-20, 20])
    ax.set_xlim([-20, 20])
    plt.show()

Where plots from the 2D case were created using the line voro('2D_file.gsd', is2D=True and the flattened 3D case were created using the line voro('3D_file.gsd').

Scipy results were produced with the following code:

import matplotlib.pyplot as plt
import scipy.spatial
import numpy as np
import gsd.hoomd

def delauny(filename, is2D=False):
    plt.figure(figsize=(20, 20))
    frame = gsd.hoomd.open(filename, mode='rb')[0]
    pos = np.copy(frame.particles.position[:, 0:2])
    pos = pos[abs(pos[:, 0]) < 20]
    pos = pos[abs(pos[:, 1]) < 20]
    tri = scipy.spatial.Delaunay(pos)
    plt.triplot(pos[:, 0], pos[:, 1], tri.simplices.copy())
    plt.scatter(pos[:, 0], pos[:, 1], s=50, color='black')
    plt.ylim([-20, 20])
    plt.xlim([-20, 20])
    plt.show()

Where plots from the 2D case were created using the line delauny('2D_file.gsd', is2D=True and the flattened 3D case were created using the line delauny('3D_file.gsd').

Error output For the 2D case, the expected output, calculated from scipy.spatial.Delauny is: Freud returns the following result:

In the 3D case, you get the expected output. scipy.spatial.Delauny returns: Freud returns a similar results:

System configuration (please complete the following information):

• OS: macOS
• Version of Python [e.g. 3.7] 3.8.1
• Version of freud [e.g. 1.0] 2.1.0

[klywang]

Here are the simulation files. ideal_tiling.gsd is the true 2D case restart_sim.gsd is the 2D analogue simulations.zip