# Computing the Smallest Enclosing Disk

I recently read in chapter 4 of Computational Geometry by de Berg et al. the problem of computing the smallest enclosing disk for a set of points.

I’ve shamelessly stolen the algorithm from there and done a simple conversion to Javascript.

The code is under the canvas-geolib GitHub repository in the geometry.js file, there’s also an example @ enclosingdisk.html. The example initially consists of three points and the smallest disk enclosing them. Click anywhere on the canvas and a new point will be added and the disk redrawn.

I don’t want to go over the general algorithm but I do want explain computing the unique disk with three points given for its boundary. In geometry.js the function “enclosingDisk3Points” takes three points and returns the unique disk that has those points on its boundary.

The below figure shows the two defining characteristics of the disk, which are (1) the disk center is centered at (x,y) and (2) the distance from the center to all three points (p0, p1, & p3) is equal i.e. the distance from the center to all three points is d.

From this computing (x,y) and d is simple. For simplicity we assume $p_0 = (0, 0)$ and $p_1 = (0, {p_y})$, also we let $p_2 = ({q}_x, {q}_y)$.

The set of equations to solve is:

$d^2 = x^2 + y^2$, $d^2 = {(x - 0)}^2 + {(y - p_y)}^2 = x^2 + y^2 - {p_y}^2 - 2\cdot p_y y$, and $d^2 = {(x - q_x)}^2 + {(y- q_y)}^2 = x^2 + y^2 + {q_x}^2 + {q_y}^2 - 2 q_x x - 2 q_y y.$

Solving for y, we have $y = p_y/2$, using this we can solve for x, which yields $x = \frac{{\|q\|}^2 - q_y p_y}{2 q_x}.$ Finally we also have $d = \sqrt{x^2 + y^2}$, which finishes our computation.

The below Javascript code implements the above computation and also adds the preprocessing steps of (1) translating the origin to $p_0$ and (2) rotating the coordinate system so that $p_1$ is of the form $p_1 = (0, {p_1}_y).$

// return the unique disk with p1, p2, and p3 as boundary points.
function enclosingDisk3Points(_p1, _p2, _p3){

var p1 = [_p1[0], _p1[1]];
var p2 = [_p2[0], _p2[1]];
var p3 = [_p3[0], _p3[1]];
if (dist(p1, p3) > dist(p1, p2)){
var p = p2;
p2 = p3;
p3 = p;
}

var p = p1;
// make p1 the origin
p2[0] = p2[0] - p1[0];
p2[1] = p2[1] - p1[1];
p3[0] = p3[0] - p1[0];
p3[1] = p3[1] - p1[1];

// apply rotation matrix to make p2.x = 0
// the rotation matrix is
// | p2[1]/dist(p2), -1 * p2[0]/dist(p2) |
// | p2[0]/dist(p2), p2[1]/dist(p2) |
//
var original_p2 = [p2[0], p2[1]];
p2[0] = 0;
p2[1] = d(original_p2);

// apply rotation matrix to p3
var original_p3 = [p3[0], p3[1]]
p3[0] = original_p2[1]/d(original_p2) * original_p3[0] - original_p2[0]/d(original_p2) * original_p3[1]
p3[1] = original_p2[0]/d(original_p2) * original_p3[0] + original_p2[1]/d(original_p2) * original_p3[1]

// the unique disk with the points p1, p2, and p3 as boundary points is
// defined by the equation y = p2.y/2 & x = (d(p3)^2 + p3.y * p2.y)/(2 * p3.x)
var y = p2[1]/2.0;
var x = (d(p3) * d(p3) - p3[1] * p2[1])/(2 * p3[0]);

// apply inverse of rotation matrix
var rotated_x = original_p2[1]/d(original_p2) * x + original_p2[0]/d(original_p2) * y
var rotated_y = -1 * original_p2[0]/d(original_p2) * x + original_p2[1]/d(original_p2) * y;

// translate back
rotated_x = rotated_x + p1[0];
rotated_y = rotated_y + p1[1];

var radius = d([rotated_x - p1[0], rotated_y - p1[1]]);