welzl algorithm

welzl.js

function shuffle(arr) {
  for (let i = arr.length - 1; i >= 0; i--) {
    let j = Math.floor(Math.random() * (i + 1));
    j = Math.max(Math.min(j, i), 0);
    const tmp  = arr[i];
    arr[i] = arr[j];
    arr[j] = tmp;
  }
  return arr;
}


export default function med (circles) {
  const data = shuffle(circles.slice());
  const n = data.length;
  let bounds = [];

  let i = 0, p, circle = null;

  while (i < n) {
    const p = data[i];
    if (circle && containsWeak(circle, p)) i++;
    else {
      bounds = extend(bounds, p);
      circle = enclose(bounds);
      i = 0;
    }
  }

  return circle;
}


function extend (S, p) {
  if (containsWeakAll(p, S)) return [p];

  // If we get here then B must have at least one element.
  for (let i = 0; i < S.length; ++i) {
    if (!contains(p, S[i])
        && containsWeakAll(from2(S[i], p), S)) {
      return [S[i], p];
    }
  }

  // If we get here then B must have at least two elements.
  for (let i = 0; i < S.length - 1; i++) {
    for (let j = i + 1; j < S.length; j++) {
      if (!contains(from2(S[i], S[j]), p) &&
          !contains(from2(S[i], p), S[j]) &&
          !contains(from2(S[j], p), S[i]) &&
          containsWeakAll(from3(S[i], S[j], p), S)) {
        return [S[i], S[j], p];
      }
    }
  }

  // If we get here then something is very wrong.
  throw new Error;
}

function contains (a, b) {
  var dr = a.r - b.r, dx = b.x - a.x, dy = b.y - a.y;
  return dr >= 0 && dr * dr > dx * dx + dy * dy;
}

function containsWeak (a, b) {
  var dr = a.r - b.r + 1e-6, dx = b.x - a.x, dy = b.y - a.y;
  return dr > 0 && dr * dr > dx * dx + dy * dy;
}

function containsWeakAll (a, S) {
  for (let i = 0; i < S.length; ++i) {
    if (!containsWeak(a, S[i])) return false;
  }
  return true;
}


function enclose(B) {
  switch (B.length) {
    case 1: return from1(B[0]);
    case 2: return from2(B[0], B[1]);
    case 3: return from3(B[0], B[1], B[2]);
  }
}

function from1(a) {
  return {
    x: a.x,
    y: a.y,
    r: a.r
  };
}

function from2(a, b) {
  const x1 = a.x, y1 = a.y, r1 = a.r;
  const x2 = b.x, y2 = b.y, r2 = b.r;
  const x21 = x2 - x1, y21 = y2 - y1, r21 = r2 - r1;
  const l = Math.sqrt(x21 * x21 + y21 * y21);

  return {
    x: (x1 + x2 + x21 / l * r21) / 2,
    y: (y1 + y2 + y21 / l * r21) / 2,
    r: (l + r1 + r2) / 2
  };
}


function from3 (a, b, c) {
  const x1 = a.x, y1 = a.y, r1 = a.r;
  const x2 = b.x, y2 = b.y, r2 = b.r;
  const x3 = c.x, y3 = c.y, r3 = c.r;
  const a2 = x1 - x2;
  const a3 = x1 - x3;
  const b2 = y1 - y2;
  const b3 = y1 - y3;
  const c2 = r2 - r1;
  const c3 = r3 - r1;

  const d1 = x1 * x1 + y1 * y1 - r1 * r1;
  const d2 = d1 - x2 * x2 - y2 * y2 + r2 * r2;
  const d3 = d1 - x3 * x3 - y3 * y3 + r3 * r3;

  const ab = a3 * b2 - a2 * b3;
  const xa = (b2 * d3 - b3 * d2) / (ab * 2) - x1;
  const xb = (b3 * c2 - b2 * c3) / ab;
  const ya = (a3 * d2 - a2 * d3) / (ab * 2) - y1;
  const yb = (a2 * c3 - a3 * c2) / ab;

  const A = xb * xb + yb * yb - 1;
  const B = 2 * (r1 + xa * xb + ya * yb);
  const C = xa * xa + ya * ya - r1 * r1;
  const r = -(A ? (B + Math.sqrt(B * B - 4 * A * C)) / (2 * A) : C / B);

  return {
    x: x1 + xa + xb * r,
    y: y1 + ya + yb * r,
    r: r
  };
}



function test (N = 10000) {
  console.time('c');
  const r = Math.sqrt(N);
  const c = new Array(N).fill(0).map(() => {
    return {
      x: (Math.random() - 0.5) * N,
      y: (Math.random() - 0.5) * N,
      r: Math.random() * r,
    };
  });

  const m = med(c);

  console.timeEnd('c');
  console.log(m);
}