docs: void run(boolean hullonly) {
      if(points.size() < 3) {
        throw new UnsupportedOperationException("There is no delaunay triangulation for less than three objects!");
      }
      int len = points.size() - 1;
      hull = new LinkedList<>();
      tris = hullonly ? null : new ArrayList<Triangle>(len);
  
      // 1. Seed point x_0
      final double[] seed;
      final int seedid = 0;
  
      final double[] sortd = new double[len];
      final int[] sorti = new int[len];
      Arrays.fill(sorti, -42); // To cause errors
      // TODO: remove duplicates.
  
      // 2. sort by squared Euclidean distance
      {
        Iterator<double[]> iter = points.iterator();
        seed = iter.next();
        for(int i = 0, j = 1; iter.hasNext(); j++, i++) {
          double dist = quadraticEuclidean(seed, iter.next());
          if(dist <= 0.) { // Duplicate.
            --len; // Decrease candidate set size
            --i; // Increase j, but not i.
            continue;
          }
          sortd[i] = dist;
          sorti[i] = j;
        }
        DoubleIntegerArrayQuickSort.sort(sortd, sorti, len);
      }
      // Detect some degenerate situations:
      if(len < 2) {
        hull.add(new IntIntPair(seedid, -1));
        if(len == 1) {
          hull.add(new IntIntPair(sorti[0], -1));
        }
        return;
      }
      assert (sortd[0] > 0);
      // final double[] seed2 = points.get(sort[0].second);
      final int seed2id = sorti[0];
  
      // 3. Find minimal triangle for these two points:
      Triangle besttri = findSmallest(seedid, seed2id, sortd, sorti, len);
      if(besttri == null) { // Degenerate
        hull.add(new IntIntPair(seedid, -1));
        hull.add(new IntIntPair(seed2id, -1));
        return;
      }
      // Note: sortd no longer accurate, recompute below!
      int start = 2; // First two points have already been processed.
  
      // 5. Make right-handed:
      besttri.makeClockwise(points);
      // Seed triangulation
      if(!hullonly) {
        tris.add(besttri);
      }
      // Seed convex hull (point, triangle)
      hull.add(new IntIntPair(besttri.a, 0));
      hull.add(new IntIntPair(besttri.b, 0));
      hull.add(new IntIntPair(besttri.c, 0));
  
      if(LOG.isDebuggingFinest()) {
        debugHull();
      }
  
      // 6. Resort from triangle circumcircle center
      double[] center = besttri.m;
      for(int i = start; i < len; i++) {
        sortd[i] = quadraticEuclidean(center, points.get(sorti[i]));
      }
      DoubleIntegerArrayQuickSort.sort(sortd, sorti, start, len);
  
      // Grow hull and triangles
      for(int i = start; i < len; i++) {
        final int pointId = sorti[i];
        final double[] newpoint = points.get(pointId);
  
        LinkedList<Triangle> newtris = hullonly ? null : new LinkedList<Triangle>();
        // We identify edges by their starting point. -1 is invalid.
        int hstart = -1, hend = -1;
        // Find first and last consecutive visible edge, backwards:
        {
          Iterator<IntIntPair> iter = hull.descendingIterator();
          IntIntPair next = hull.getFirst();
          double[] nextV = points.get(next.first);
          for(int pos = hull.size() - 1; iter.hasNext(); pos--) {
            IntIntPair prev = iter.next();
            double[] prevV = points.get(prev.first);
            // Not yet visible:
            if(hend < 0) {
              if(leftOf(prevV, nextV, newpoint)) {
                hstart = hend = pos;
                if(!hullonly) {
                  // Clockwise, A is new point!
                  Triangle tri = new Triangle(pointId, next.first, prev.first);
                  assert (tri.isClockwise(points));
                  assert (prev.second >= 0);
                  tri.updateCircumcircle(points);
                  tri.bc = prev.second;
                  newtris.addFirst(tri);
                }
              }
            }
            else {
              if(leftOf(prevV, nextV, newpoint)) {
                hstart = pos;
                // Add triad:
                if(!hullonly) {
                  // Clockwise, A is new point!
                  Triangle tri = new Triangle(pointId, next.first, prev.first);
                  assert (tri.isClockwise(points));
                  assert (prev.second >= 0);
                  tri.updateCircumcircle(points);
                  tri.bc = prev.second;
                  newtris.addFirst(tri);
                }
              }
              else {
                break;
              }
            }
            next = prev;
            nextV = prevV;
          }
        }
        // If the last edge was visible, we also need to scan forwards:
        if(hend == hull.size() - 1) {
          Iterator<IntIntPair> iter = hull.iterator();
          IntIntPair prev = iter.next();
          double[] prevV = points.get(prev.first);
          while(iter.hasNext()) {
            IntIntPair next = iter.next();
            double[] nextV = points.get(next.first);
            if(leftOf(prevV, nextV, newpoint)) {
              hend++;
              // Add triad:
              if(!hullonly) {
                // Clockwise, A is new point!
                Triangle tri = new Triangle(pointId, next.first, prev.first);
                assert (tri.isClockwise(points));
                assert (prev.second >= 0);
                tri.updateCircumcircle(points);
                tri.bc = prev.second;
                newtris.addLast(tri);
              }
            }
            else {
              break;
            }
            prev = next;
            prevV = nextV;
          }
        }
        assert (hstart >= 0 && hend >= hstart);
        // Note that hend can be larger than hull.size() now, interpret as
        // "hend % hull.size()"
        // Update hull, remove points
        final int firsttri, lasttri;
        if(hullonly) {
          firsttri = lasttri = -1;
        }
        else {
          final int tristart = tris.size();
          firsttri = tristart;
          lasttri = tristart + newtris.size() - 1;
        }
        final int hullsize = hull.size();
        if(LOG.isDebuggingFinest()) {
          LOG.debugFinest("Size: " + hullsize + " start: " + hstart + " end: " + hend);
        }
        if(hend < hullsize) {
          ListIterator<IntIntPair> iter = hull.listIterator();
          int p = 0;
          // Skip
          for(; p <= hstart; p++) {
            iter.next();
          }
          // Remove
          for(; p <= hend; p++) {
            iter.next();
            iter.remove();
          }
          // Insert, and update edge->triangle mapping
          iter.add(new IntIntPair(pointId, lasttri));
          iter.previous();
          if(!hullonly) {
            (iter.hasPrevious() ? iter.previous() : hull.getLast()).second = firsttri;
          }
        }
        else {
          ListIterator<IntIntPair> iter = hull.listIterator();
          // Remove end
          int p = hullsize;
          for(; p <= hend; p++) {
            iter.next();
            iter.remove();
          }
          // Insert
          iter.add(new IntIntPair(pointId, lasttri));
          // Wrap around
          p -= hullsize;
          IntIntPair pre = null;
          for(; p <= hstart; p++) {
            pre = iter.next();
          }
          assert (pre != null);
          pre.second = firsttri;
          // Remove remainder
          while(iter.hasNext()) {
            iter.next();
            iter.remove();
          }
        }
        if(LOG.isDebuggingFinest()) {
          debugHull();
        }
        if(!hullonly) {
          final int tristart = tris.size();
          // Connect triads (they are ordered)
          Iterator<Triangle> iter = newtris.iterator();
          for(int o = 0; iter.hasNext(); o++) {
            // This triangle has num tristart + o.
            Triangle cur = iter.next();
            cur.ca = o > 0 ? tristart + o - 1 : -1; // previously added triangle
            cur.ab = iter.hasNext() ? tristart + o + 1 : -1; // next triangle
            // cur.bc was set upon creation
            assert (cur.bc >= 0);
            Triangle other = tris.get(cur.bc);
            Orientation orient = cur.findOrientation(other);
            assert (orient != null) : "Inconsistent triangles: " + cur + " " + other;
            switch(orient){
            case ORIENT_BC_BA:
              assert (other.ab == -1) : "Inconsistent triangles: " + cur + " " + other;
              other.ab = tristart + o;
              break;
            case ORIENT_BC_CB:
              assert (other.bc == -1) : "Inconsistent triangles: " + cur + " " + other;
              other.bc = tristart + o;
              break;
            case ORIENT_BC_AC:
              assert (other.ca == -1) : "Inconsistent triangles: " + cur + " " + other;
              other.ca = tristart + o;
              break;
            default:
              assert (cur.isClockwise(points));
              assert (other.isClockwise(points));
              throw new RuntimeException("Inconsistent triangles: " + cur + " " + other + " size:" + tris.size());
            }
            tris.add(cur);
          }
          assert (tris.size() == lasttri + 1);
        }
      }
      // Now check for triangles that need flipping.
      if(!hullonly) {
        final int size = tris.size();
        long[] flippedA = BitsUtil.zero(size), flippedB = BitsUtil.zero(size);
        // Initial flip
        if(flipTriangles(flippedA) > 0) {
          for(int iterations = 1; iterations < 1000; iterations += 2) {
            if(LOG.isDebuggingFinest()) {
              debugHull();
            }
            if(flipTriangles(flippedA, flippedB) == 0) {
              break;
            }
            if(LOG.isDebuggingFinest()) {
              debugHull();
            }
            if(flipTriangles(flippedB, flippedA) == 0) {
              break;
            }
          }
        }
      }
    }