Tag Archives: javascript

Computing the Smallest Enclosing Disk

Posted by

2

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.

center point

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]]);
    return [[rotated_x, rotated_y], radius];
    
}

Convex Hulls

Posted by

1

The algorithm I’m going to showcase is taken from the book “Computational Geometry algorithms and applications” by Berg, Cheong, Kreveld, and Overmars. They give an excellent analysis of the algorithm in the book so I’ll only go over the highlights and show my implementation in Javascript using the html5 canvas.

The problem statement is given a set of points, S = \{p_1, p_2, \ldots, p_n\} in R^2 compute the convex hull of S. The input of the algorithm is the set, S, of points and the output of the algorithm is the set of points that are the vertices of the convex hull of S.

Our algorithm will use a standard design technique to generate what’s called an incremental algorithm. That is we compute the solution for the first p_1, \ldots, p_i points then add the point p_{i+1} and compute the new solution using the previous solution.

Algorithm ConvexHull(P) (the pseudocode is borrowed from “Computational Geometry algorithms and applications”)
Input. A set P of points in the plane
Output. A list containing the vertices of \mathcal{CH}(P) in clockwise order.
1. Sort the points by x-coordinate, resulting in a sequence p_1,\ldots, p_n.
2. Put the points p_1 and p_2 in a list \mathcal{L}_\textrm{upper}, with p_1 as the first point.
3. for i \leftarrow 3 to n
4.        do Append p_i to \mathcal{L}_\textrm{upper}.
5.                        while \mathcal{L}_\textrm{upper} contains more than two points and the last three points in \mathcal{L}_\textrm{upper} do not make a right turn
6.                           do Delete the middle of the last three points from \mathcal{L}_\textrm{upper}.
7. Put the points p_n and p_{n-1} in a list \mathcal{L}_\textrm{lower}, with p_n as the first point.
8. for i \leftarrow n -2 downto 1
9.           do Append p_i to \mathcal{L}_\textrm{lower}.
10.                     while \mathcal{L}_\textrm{lower} contains more than 2 points and the last three points in \mathcal{L}_\textrm{lower} do not make a right turn
11.                     do Delete the middle of the last three points from \mathcal{L}_\textrm{lower}.
12. Remove the first and last point from \mathcal{L}_\textrm{lower}. to ovoid duplication of the points where the upper and lower hull meet.
13. Append \mathcal{L}_\textrm{lower} to \mathcal{L}_\textrm{upper}, and call the resulting list \mathcal{L}.
14. return \mathcal{L}.

You can view an implementation using Javascript and the html5 canvas element, here. To add a new point simply click the left mouse button and it’ll add a point where the mouse is.

Right and Left Turns

Posted by

0

In my previous post on line segment intersection I introduced the two dimensional cross product as v1 \times v2 = {v1}_x \cdot {v2}_y - {v2}_x \cdot {v1}_y. The cross product can also be used to determine if a set of three points p_1, p_2, \textrm{ and } p_3 make a right turn.

First note that if v_1 \times v_2 > 0 then the angle between v_1 and v_2 is strictly less than \pi. For example the two vectors (1, 0) \times (0, 1) = 1 \cdot 1 - 0 = 1 > 0. Similarly if v_1 \times v_2 < 0 then the angle between them is strictly greater than \pi. In the below image I've shown an example where the three points p_1, p_2, \textrm{ and } p_3 make a right turn.

Points p1, p2, p3 making a right turn

To mathematically determine if the points make a right turn we let v_1 = p_1 - p_2 \textrm{ and } v_2 = p_3 - p_2. Then taking the cross product of v_1 \textrm{ and } v_2, we have v_1 \times v_2 > 0 which implies that the angle between them is less than \pi that is they make a right turn.

This method for determining whether the points make a right or left turn is very useful for determining the convex hull of a set of points.

I’ve posted an example webpage, here, where using javascript the lines change color depending on whether the turn is to the right or the left. Please note that the page uses the html5 canvas element so it will not work in internet explorer 7, or 6.

Line Segment Intersection

Posted by

11

We want to determine if two line segments, s_1 and s_2, intersect and if they do intersect, the point of intersection. First we analyze the problem. Case five in the list below is difficult to analyze which is why we leave it to the end.

  1. The segments are not parallel (easy)
  2. The segments are parallel and do not intersect (easy)
  3. The segments are collinear and do not intersect (easy)
  4. The segments are collinear and do intersect (easy)
  5. The segments are nearly parallel or nearly collinear (hard)

In this post I handle cases one through four. I leave case five for later posts.

For the problem setup let segment 1 be given by

s_1 = p + t \cdot r where 0 <= t <= 1 and p, r are 2d vectors, similarly let segment 2 be given by

s_2 = q + u \cdot s where 0 <= u <= 1 and q, s are 2d vectors.

This representation of the line segments yields a very natural method for computing their intersection.

Computing cross products is the heart of the algorithm for determining intersections. The two dimensional cross product of p_1 \times p_2 is given by

\textrm{det} \left( \begin{array} {cc} x_1 & x_2 \\ y_1 & y_2\end{array} \right) = x_1y_2 - x_2y_1 = -p_2 \times p_1

If the cross product is zero then the two vectors are collinear that is pointing in either the same direction or in opposite directions.

The two lines intersect when p + t \cdot r = q + u \cdot s, by crossing both sides with s we get

(p + t \cdot r) \times s = (q + u \cdot s) \times s = q \times s + u \cdot s \times s = q \times s.

From this we solve for t obtaining

t = \frac{(q - p) \times s}{r \times s}.

We can similarly solve for u, obtaining

u = \frac{(q -p) \times r}{r \times s}.

If both t and u are between 0 and 1, then the two line segments intersect and the intersection point is given by p + t \cdot r, where the value of t is the one we solved for. If either t or u are not between 0 and 1 then the line segments do not intersect.

But suppose r \times s = 0, then we can’t solve for t and u. This is because r \times s = 0 means that r and s are parallel which means the two line segments are parallel. Please also note that r \times s = 0 implies that r and s are scalar multiplies of each other. If (p -q) \times r \not= 0 then the segments are parallel but not collinear hence they do not intersect. If (p -q) \times r = 0 then the segments are collinear and we can simply project both of them onto to x-axis and determine if their projections intersect.

For case 5 when r \times s is close to zero the analysis needs to be well thought out and I’ll leave it to future posts.

I’ve posted a sample implementation in Javascript using the html5 canvas at
http://cloud.github.com/downloads/bjwbell/canvas-geolib/main.html
(please note it won’t work in internet explorer since it doesn’t support the canvas element).

Credit Gareth Rees for his post at stackoverflow.com which is based on the method for 3D line intersection algorithm from the article “Intersection of two lines in three-space” by Ronald Graham, published in Graphics Gems, page 304.