Sorting issue?

[Av3r3tt]

Hi, I am trying to refactor your code into Swift, and am running into an issue. Is it possible your quicksort algorithm is not working as intended? I run into issues with the following set of blue noise points (but other sets of points also fail, either not sorting or not triangulating); your code triangulates successfully only 3 out of 4 times or so. I am still pulling my hair out wondering whether it is my code... could someone verify with this simple set of points?

coords = [12.0, 100.0, 181.0, 13.0, 72.0, 7.0, 274.0, 23.0, 133.0, 13.0, 153.0, 66.0, 248.0, 71.0, 56.0, 59.0, 97.0, 98.0, 207.0, 86.0, 496.0, 9.0, 546.0, 1.0, 427.0, 23.0, 384.0, 7.0, 332.0, 20.0, 507.0, 62.0, 450.0, 94.0, 578.0, 102.0, 364.0, 61.0, 317.0, 103.0, 400.0, 106.0];

According to my test output, your quicksort sorts these coords in ids into this order (which shows the point number, and the dist() function result to cx/cy: 3:=1371.25 19:=2934.25 14:=2491.25 18:=4817.25 9:=8800.25 13:=10083.25 20:=13781.25 1:=14636.25 5:=20320.25 12:=18354.25 16:=25665.25 4:=27884.25 8:=41184.25 10:=42381.25 15:=45016.25 2:=51891.25 7:=57151.25 11:=65757.25 0:=82251.25 17:=82441.25

[mourner]

Can you set this up as a minimal test case on the JS side so that I could verify?

[Av3r3tt]

I'm not very versed in JavaScript, so I just copied and pasted the code into this website (https://www.w3schools.com/js/tryit.asp?filename=tryjs_syntax_variables) and tried it out live. Using the script this way:

<script>
var result;

var coords = [12.0, 100.0, 181.0, 13.0, 72.0, 7.0, 274.0, 23.0, 133.0, 13.0, 153.0, 66.0, 248.0, 71.0, 56.0, 59.0, 97.0, 98.0, 207.0, 86.0, 496.0, 9.0, 546.0, 1.0, 427.0, 23.0, 384.0, 7.0, 332.0, 20.0, 507.0, 62.0, 450.0, 94.0, 578.0, 102.0, 364.0, 61.0, 317.0, 103.0, 400.0, 106.0];

    var minX = Infinity;
    var minY = Infinity;
    var maxX = -Infinity;
    var maxY = -Infinity;

    var n = coords.length >> 1;
    var ids = new Uint32Array(n);

    for (var i = 0; i < n; i++) {
        var x = coords[2 * i];
        var y = coords[2 * i + 1];
        if (x < minX) minX = x;
        if (y < minY) minY = y;
        if (x > maxX) maxX = x;
        if (y > maxY) maxY = y;
        ids[i] = i;
    }

    var cx = (minX + maxX) / 2;
    var cy = (minY + maxY) / 2;

    var minDist = Infinity;
    var i0, i1, i2;

    // pick a seed point close to the centroid
    for (i = 0; i < n; i++) {
        var d = dist(cx, cy, coords[2 * i], coords[2 * i + 1]);
        if (d < minDist) {
            i0 = i;
            minDist = d;
        }
    }

    minDist = Infinity;

    // find the point closest to the seed
    for (i = 0; i < n; i++) {
        if (i === i0) continue;
        d = dist(coords[2 * i0], coords[2 * i0 + 1], coords[2 * i], coords[2 * i + 1]);
        if (d < minDist && d > 0) {
            i1 = i;
            minDist = d;
        }
    }

    var minRadius = Infinity;

    // find the third point which forms the smallest circumcircle with the first two
    for (i = 0; i < n; i++) {
        if (i === i0 || i === i1) continue;

        var r = circumradius(
            coords[2 * i0], coords[2 * i0 + 1],
            coords[2 * i1], coords[2 * i1 + 1],
            coords[2 * i], coords[2 * i + 1]);

        if (r < minRadius) {
            i2 = i;
            minRadius = r;
        }
    }

    if (minRadius === Infinity) {
        throw new Error('No Delaunay triangulation exists for this input.');
    }

    // swap the order of the seed points for counter-clockwise orientation
    if (area(coords[2 * i0], coords[2 * i0 + 1],
        coords[2 * i1], coords[2 * i1 + 1],
        coords[2 * i2], coords[2 * i2 + 1]) < 0) {

        var tmp = i1;
        i1 = i2;
        i2 = tmp;
    }

    var i0x = coords[2 * i0];
    var i0y = coords[2 * i0 + 1];
    var i1x = coords[2 * i1];
    var i1y = coords[2 * i1 + 1];
    var i2x = coords[2 * i2];
    var i2y = coords[2 * i2 + 1];

    var center = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y);
    this._cx = center.x;
    this._cy = center.y;

    // sort the points by distance from the seed triangle circumcenter
    quicksort(ids, coords, 0, ids.length - 1, center.x, center.y);

function dist(ax, ay, bx, by) {
    var dx = ax - bx;
    var dy = ay - by;
    return dx * dx + dy * dy;
}

function area(px, py, qx, qy, rx, ry) {
    return (qy - py) * (rx - qx) - (qx - px) * (ry - qy);
}

function inCircle(ax, ay, bx, by, cx, cy, px, py) {
    ax -= px;
    ay -= py;
    bx -= px;
    by -= py;
    cx -= px;
    cy -= py;

    var ap = ax * ax + ay * ay;
    var bp = bx * bx + by * by;
    var cp = cx * cx + cy * cy;

    return ax * (by * cp - bp * cy) -
           ay * (bx * cp - bp * cx) +
           ap * (bx * cy - by * cx) < 0;
}

function circumradius(ax, ay, bx, by, cx, cy) {
    bx -= ax;
    by -= ay;
    cx -= ax;
    cy -= ay;

    var bl = bx * bx + by * by;
    var cl = cx * cx + cy * cy;

    if (bl === 0 || cl === 0) return Infinity;

    var d = bx * cy - by * cx;
    if (d === 0) return Infinity;

    var x = (cy * bl - by * cl) * 0.5 / d;
    var y = (bx * cl - cx * bl) * 0.5 / d;

    return x * x + y * y;
}

function circumcenter(ax, ay, bx, by, cx, cy) {
    bx -= ax;
    by -= ay;
    cx -= ax;
    cy -= ay;

    var bl = bx * bx + by * by;
    var cl = cx * cx + cy * cy;

    var d = bx * cy - by * cx;

    var x = (cy * bl - by * cl) * 0.5 / d;
    var y = (bx * cl - cx * bl) * 0.5 / d;

    return {
        x: ax + x,
        y: ay + y
    };
}

// create a new node in a doubly linked list
function insertNode(coords, i, prev) {
    var node = {
        i: i,
        x: coords[2 * i],
        y: coords[2 * i + 1],
        t: 0,
        prev: null,
        next: null,
        removed: false
    };

    if (!prev) {
        node.prev = node;
        node.next = node;

    } else {
        node.next = prev.next;
        node.prev = prev;
        prev.next.prev = node;
        prev.next = node;
    }
    return node;
}

function removeNode(node) {
    node.prev.next = node.next;
    node.next.prev = node.prev;
    node.removed = true;
    return node.prev;
}

function quicksort(ids, coords, left, right, cx, cy) {
    var i, j, temp;

    if (right - left <= 20) {
        for (i = left + 1; i <= right; i++) {
            temp = ids[i];
            j = i - 1;
            while (j >= left && compare(coords, ids[j], temp, cx, cy) > 0) ids[j + 1] = ids[j--];
            ids[j + 1] = temp;
        }
    } else {
        var median = (left + right) >> 1;
        i = left + 1;
        j = right;
        swap(ids, median, i);
        if (compare(coords, ids[left], ids[right], cx, cy) > 0) swap(ids, left, right);
        if (compare(coords, ids[i], ids[right], cx, cy) > 0) swap(ids, i, right);
        if (compare(coords, ids[left], ids[i], cx, cy) > 0) swap(ids, left, i);

        temp = ids[i];
        while (true) {
            do i++; while (compare(coords, ids[i], temp, cx, cy) < 0);
            do j--; while (compare(coords, ids[j], temp, cx, cy) > 0);
            if (j < i) break;
            swap(ids, i, j);
        }
        ids[left + 1] = ids[j];
        ids[j] = temp;

        if (right - i + 1 >= j - left) {
            quicksort(ids, coords, i, right, cx, cy);
            quicksort(ids, coords, left, j - 1, cx, cy);
        } else {
            quicksort(ids, coords, left, j - 1, cx, cy);
            quicksort(ids, coords, i, right, cx, cy);
        }
    }
}

function compare(coords, i, j, cx, cy) {
    var d1 = dist(coords[2 * i], coords[2 * i + 1], cx, cy);
    var d2 = dist(coords[2 * j], coords[2 * j + 1], cx, cy);
    return (d1 - d2) || (coords[2 * i] - coords[2 * j]) || (coords[2 * i + 1] - coords[2 * j + 1]);
}

function swap(arr, i, j) {
    var tmp = arr[i];
    arr[i] = arr[j];
    arr[j] = tmp;
}

n = ids.length;
result=cx + ' ' +cy + "<BR>"
for (i=0;i<n;i++) {
    j = ids[i];
    result = result + j + ':=' + dist(coords[2 * j], coords[2 * j + 1], cx, cy)+'<BR>';
}
result = result + ids;

document.getElementById("demo").innerHTML = result;
</script>

[mourner]

You calculate it for the wrong cx and cy. These are for the center of the bbox to find a good starting seed. But quicksort sorts points around center which is calculated later (see this._cx and this._cy)

[mourner]

The Delaunay looks correct for your input coords, and it's always the same (no element of randomness in the algorithm):

[Av3r3tt]

Thanks Vlad - I just discovered this too,...but my own (Swift) code still keeps running in an eternal loop. I guess the error is still somewhere else in my code - sorry for wasting your time!

[redblobgames]

Confirmed. I see the same thing (test page) — 

[MaxGraey]

@Av3r3tt By the way, for large data and 32-bit index types current solution for finding middle can cause to overflow:

var median = (left + right) >> 1;

Better use this:

var mid = left + (right - left) / 2;

Just keep in mind

[Av3r3tt]

Thanks Amit, thanks Max. As I said I'm trying to to refactor this all into Swift because of this algorithm's pure speed. Alas, a quarter of the time (even with small number of points) it seems to fail and gets stuck in an eternal loop after sorting. And three quarters of the time, it succeeds and produces the exact same version of triangulation my other (slower) algorithm does. It's weird... but clearly the error is mine (and is still there in my implementation). Thanks again everybody!

[mourner]

@MaxGraey to make it not safe to use a bitwise shift for indices, we'd have to run Delaunator on an array of at least 1 billion points — JS engines wouldn't even allow you to allocate that kind of array. So I think we're safe here.

[MaxGraey]

For JS yes, because op left + right actually add two doubles which will overflow only if sum >= 2 ^ 53 - 1. But for swift and int32 type this could be mater

[Av3r3tt]

Thanks all for your support. Realizing that my error was elsewhere, it took me a bit but I found it... now for refactoring Voronoi! :-D

[luizv]

@Av3r3tt, did you succeed with your port for Swift?

[Av3r3tt]

@Av3r3tt, did you succeed with your port for Swift?

luizv, yes I did -- but I had to switch to other stuff for several months now, but my swift code works -- but it is not up to date with what happened with delanaunator since, ie not as speed optimized as the current JS version. And I haven't tested it as extensively yet. I was working on a project that created blue noise first, and then did triangulating and voronoi-ing from there, so it might well be that various cases would not delaunay well under my version.