A routine to calculate whether two points are adjacent in the Delaunay Triangulation of a space.

DelaunayAdjacency.java

import com.badlogic.gdx.math.Vector2;
import com.badlogic.gdx.utils.Array;

/**
 * This class computes whether two points are adjacent in the Delaunay Triangulation of a set of points.
 * Throughout the implementation, it makes the assumption that points can be compared with reference equality for uniqueness.
 */
public abstract class DelaunayAdjacency {

	private static final float EPSILON = 0.000001f;

	/**
	 * This function computes just enough of the Delaunay Triangulation
	 * of the space to calculate if points are neighbors.
	 *
	 * It hinges on the following principle:
	 * Given three points, A, B, and C,
	 * Let the area inside the circumcircle of line AB be known as S.
	 * If C is inside S, then The angle ACB is obtuse.
	 * If C is outside S, then the angle ACB is acute.
	 *
	 * Let S1 and S2 be the semicircular parts of S on either side of AB.
	 * Further, let P1 be the point in S1 such that the circumcircle of triangle P1AB contains no other points in P1.
	 * Let P2 be the similar point for S2.
	 * Note that this is equivalent to finding the point which maximizes the area of the circumcircle formed by itself and AB.
	 * Note that all angles of a chord in the same segment of a circle are equal.
	 * Therefore P1 is the point in S1 that makes the maximum angle AS1B.
	 *
	 * If neither P1 nor P2 exist, then the points that form the quadrilateral containing AB must both have acute angles on them.
	 * TODO: Prove that there must be a quadrilateral (or edge triangle) containing AB in this case.
	 * Therefore there is an edge between AB.
	 *
	 * Note that the circumcircle of the triangle formed by A, B, and any point in S1 must contain all of S2, and vice versa.
	 *
	 * If both P1 and P2 exist, then they must be contained within the circumcircle of the other with AB.
	 * Therefore AB cannot be connected.
	 *
	 * If only one of P1 and P2 exists, WLOG let that be PX.
	 * TODO: Prove that if AB are connected, PXAB must be a triangle in the triangulation.
	 * If there is any other point within the circumcircle of triangle PXAB, then AB cannot be connected.
	 * Otherwise they are connected.
	 *
	 * Math is awesome.
	 */
	public boolean arePointsNeighbors(Vector2 pointA, Vector2 pointB) {
		float x1 = pointA.x;
		float y1 = pointA.y;
		float x2 = pointB.x;
		float y2 = pointB.y;

		float midX = (x1 + x2) * 0.5f;
		float midY = (y1 + y2) * 0.5f;
		float rx = x1 - midX;
		float ry = y1 - midY;
		float radius = (float) Math.sqrt(rx * rx + ry * ry);

		Array<Vector2> pointsNearAB = getPointsInCircle(midX, midY, radius);

		Vector2 bestSideA = null;
		Vector2 bestSideB = null;
		float bestAngle = -Float.MAX_VALUE; // we only have to track the best angle on one side, because if we have angles on both sides we can early exit.
		float x3=x2, y3=y2; // Java forces these to be initialized.  This initialization guarantees that if the impossible happens, false will be returned.

		float abx = x2 - x1;
		float aby = y2 - y1;

		for (int c = 0, n = pointsNearAB.size; c < n; c++) {
			Vector2 point = pointsNearAB.get(c);
			if (point == pointA || point == pointB) continue;

			float x = point.x;
			float y = point.y;

			// Calculate difference vectors for the edges
			float dx1 = x - x1;
			float dy1 = y - y1;

			// Figure out if the point is on side A or side B.
			// From http://stackoverflow.com/questions/1560492/how-to-tell-whether-a-point-is-to-the-right-or-left-side-of-a-line
			// 'Should you find yourself in a situation where rounding error on this test
			//  is causing you problems, you will want to look up Jon Shewchuk's
			//  "Fast Robust Predicates for Computational Geometry"'
			float cross = abx * dy1 - aby * dx1;
			// Let side A be the side on which cross > 0.

			// if we have both a side A and side B, we can exit now.  They can't be connected.
			if (cross > 0) {
				if (bestSideB != null) return false;
			} else {
				if (bestSideA != null) return false;
			}

			// Otherwise we actually have to figure out whether this point has a better circumcircle than the existing one.
			float dx2 = x - x2;
			float dy2 = y - y2;

			// Normalize the vectors, then take their dot product to get the angle.
			float lenA = (float)Math.sqrt(dx1 * dx1 + dy1 * dy1);
			if (lenA != 0) {
				dx1 /= lenA;
				dy1 /= lenA;
			}
			float lenB = (float)Math.sqrt(dx2 * dx2 + dy2 * dy2);
			if (lenB != 0) {
				dx2 /= lenB;
				dy2 /= lenB;
			}

			float dot =  dx1 * dx2 + dy1 * dy2; // da dot db
			float angle = -dot; // way faster than acos, and equivalent for comparison.

			if (angle > bestAngle) {
				bestAngle = angle;
				x3 = x;
				y3 = y;
				if (cross > 0) {
					bestSideA = point;
				} else {
					bestSideB = point;
				}
			}
		}

		// At this point we cannot have both bestSideA and bestSideB, because of the early exit condition in the loop.
		Vector2 best;
		if (bestSideA == null) {
			if (bestSideB == null) return true; // If there are no other points in the circumcircle of AB, they must be connected.
			best = bestSideB;
		} else {
			best = bestSideA;
		}

		// We have three points, construct the circumcircle
		// This code copied from LibGDX's DelaunayTriangulator class.
		float xc, yc;
		float y1y2 = Math.abs(y1 - y2);
		float y2y3 = Math.abs(y2 - y3);
		if (y1y2 < EPSILON) {
			if (y2y3 < EPSILON) return false; // points are essentially collinear on y.  The circle is therefore infinite and must contain another point.
			float m2 = -(x3 - x2) / (y3 - y2);
			float mx2 = (x2 + x3) / 2f;
			float my2 = (y2 + y3) / 2f;
			xc = (x2 + x1) / 2f;
			yc = m2 * (xc - mx2) + my2;
		} else {
			float m1 = -(x2 - x1) / (y2 - y1);
			float mx1 = (x1 + x2) / 2f;
			float my1 = (y1 + y2) / 2f;
			if (y2y3 < EPSILON) {
				xc = (x3 + x2) / 2f;
				yc = m1 * (xc - mx1) + my1;
			} else {
				float m2 = -(x3 - x2) / (y3 - y2);
				float mx2 = (x2 + x3) / 2f;
				float my2 = (y2 + y3) / 2f;
				xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
				yc = m1 * (xc - mx1) + my1;
			}
		}

		float dx = x2 - xc;
		float dy = y2 - yc;
		float r = (float) Math.sqrt(dx * dx + dy * dy);

		// Otherwise, if there are no other points in the circumcircle, the points are connected.
		return !arePointsInCircle(xc, yc, r, pointA, pointB, best);
	}

    /** Returns an Array of all points within the given circle.
     * This class uses reference comparison to check for point equality, so if
     * this list contains the vectors passed as pointA and pointB to arePointsNeighbors,
     * it must return exactly those objects. */
    protected abstract Array<Vector2> getPointsInCircle(float xc, float yc, float radius);

    /** Returns true if there are points in the given circle which are not one of the `ignore` parameters.
     * Note that the circle passed in may be extremely large, to the point where even listing all the points in it is unreasonable. */
    protected abstract boolean arePointsInCircle(float xc, float yc, float radius, Vector2 ignore0, Vector2 ignore1, Vector2 ignore2);
}