Implementing the Closest Pair Algorithm in OpenCL and WebCL

I was looking for a fun geometry problem to solve and came across the closest pair of points problem.

The problem statement is:
Given a set of points, in 2D, compute the closest pair of points.

The brute force algorithm takes O(n^2) time. There’s a better solution described at the Wikipedia page Closest pair of points problem, which takes O(n \cdot \textrm{log}(n)) time.
As far as I know, there aren’t any posted solutions using OpenCL to compute the closest pair. So I’ve implemented one, it’s posted at closest-pair.

There are a few challenges in adapting the algorithm for OpenCL. Namely, we can’t use recursion so we must convert the recursive algorithm to a procedural one. In this case it’s not to complicated because the structure of the solution is easily converted from the top-down recursive algorithm to a bottom-up parallel algorithm. There are some tricky issues when the number of points are not a power of 2. I’ve commented the code for those cases.

The second challenge is using WebCL. WebCL has an additional restriction that you can’t pass structures between Javascript and OpenCL kernels. Because of this I had to use dumb arrays of simple uint’s instead of using arrays of uint2, uint3, and uint4. This made the code more verbose. To help reduce the verbosity I added macros in the OpenCL ckStrips kernel. The WebCL version is posted at

I hope the solution is useful to someone. Enjoy reading the code, it requires careful thought.
Check it out at closest-pair or view the WebCL one online at closest_pair_ocl.html.

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.

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];

Solving Tangrams Using JTS

The project, 2dfit, solves Tangram puzzles using the Java Topology Suite (JTS). The algorithm implementation is based on what I outlined in item 3 of the post ““.

The implementation difficulties are from using floating point arithmetic, which is not robust for geometric operations. The JTS library attempts to minimize this by a coordinate snapping technique. But for the operations used in solving Tangrams the provided snapping was not sufficient.

There’s an option in JTS to specify the snapping tolerance (it has a fairly small default). I added small wrapper functions for the two operations of Boolean intersection and Boolean difference. The wrapper functions apply successively larger snapping tolerances up to a factor of epsilon, where epsilon = 1e-5. The below code shows the wrapper functions (in the code, g1 is the Tangram and g2 is a puzzle piece).

    public static Geometry SemiRobustGeoOp(Geometry g1, Geometry g2, int op) throws Exception {
        double e1 = EPSILON/10;
        double snapTolerance = GeometrySnapper.computeOverlaySnapTolerance(g1, g2);
        while (snapTolerance < e1) {
            try {
                Geometry[] geos = GeometrySnapper.snap(g1, g2, snapTolerance);
                switch (op) {                    
                case DIF_OP:   // difference
                    return geos[0].difference(geos[1]);
                case UNION_OP: // union
                    return geos[0].union(geos[1]);
                    throw new Exception("unhandled semirobustgeoop: " + op);
            } catch (TopologyException e){
                snapTolerance *= 2;
        return null;

    public static boolean SemiRobustGeoPred(Geometry g1, Geometry g2, int pred) throws Exception {
        double e1 = EPSILON/10;
        double snapTolerance = GeometrySnapper.computeOverlaySnapTolerance(g1, g2);
        while (snapTolerance < e1) {
            try {
                Geometry[] geos = GeometrySnapper.snap(g1, g2, snapTolerance);
                switch (pred) {                    
                case COVER_PRED: // 
                    return geos[0].covers(geos[1]);
                    throw new Exception("unhandled semirobustgeopred: " + pred);
            } catch (TopologyException e){
                snapTolerance *= 2;
        return false;

Using the wrapper functions was key to a more robust implementation. The below figure shows a solved Tangram puzzle (from, in the figure the puzzle pieces are labeled l1.dat, l2.dat,…, l7.dat (I was lazy in naming the files). It’s the result of running FitTest_ToSingleLargeTriange() in and plotting the result using gnuplot.

I used the symmetry of each puzzle piece and a heuristic for choosing which puzzle piece to fit in reducing the number of permutations used for solving a Tangram.

For the puzzle piece symmetry, I used that the square is completely symmetric so only its first line segment needs to be used when fitting it. The triangle pieces are only partially symmetric so two of their three line segments need to be tried.

The below figure illustrates the symmetry of the square:

I used two heuristics for which pieces to try fitting, first try larger pieces before smaller ones and two skip a piece if another identical piece has already failed to be fitted (there are two identical small triangles and similarly two identical large triangles).

With the above two optimizations it takes ~1min to solve a Tangram. Without the optimizations the algorithm did not complete for the Tangrams I tried.

Quick Heuristic For Tangram Polygon Intersection

For the combinatorial approach to solving Tangram puzzles one frequent operation is to test if a puzzle piece, p, is contained in the silhouette s. Before this test is made the puzzle piece, p, is translated to a vertex, v, of s and one of the edges of p is aligned with an edge, e, coming out of vertex v.

For clarity I’ve included a figure below illustrating this.

puzzle piece and silhouette

In the above figure it’s obvious that the puzzle piece, p, does not fit in the silhouette. One way to see this computationally is to look at the edge, e, and notice that at vertex, v_1, the next edge of the silhouette, s, makes a left turn with respect to the edge e. This means that if the length of e is less than the length of all the edges in p then there is no possible way to fit p in s. That is as long as we have the condition that one of the edges of p is aligned with the edge e of s.

The above statement provides a quick way to skip vertices of the silhouette s when trying to fit the puzzle piece p in the silhouette s. Of course before skipping the vertex the same test needs to be done on the other edge coming out of v. Also we note that the above test only works if the next edge coming out of v_1 makes a left turn with respect to the edge e, otherwise the full test for polygon intersection has to be run. Also if one of the edges of the puzzle piece p has length less than or equal to the length of e then the full test for polygon intersection has to be run.

In summary this is a very specific heuristic that in general isn’t very useful but if you ever find yourself needing it, it can be just what the doctor ordered.

A Heuristic Solution to the Tangram Puzzle

In this post I present a short summary of the paper “A Heuristic Solution to the Tangram Puzzle” by E.S. Deutsch and K. C. Hayes Jr. I’ve uploaded a scanned copy of the paper to scribd here. Uploading the scanned copy probably violates some copyright law so if they or anyone else contacts me I’ll remove the document. Ironically it’s taken me over two-years since my post here to finally read the paper.

In their paper they talk about two possible approaches to solving tangrams. One is the combinatorial approach that I’ve taken and the other is the heuristic method that they describe in their paper. The heuristic method that they use is a little clever, what they do is take the outline of the tangram puzzle and then try to break it out into sub-puzzles. The following figure illustrates where two puzzles have both been separated out into sub-puzzles. In both cases one of the sub-puzzles exactly matches a tangram piece thus allowing the algorithm to solve that sub-puzzle.

Puzzle Illustration

The dashed lines in the figure show where the algorithm is splitting the sub-puzzles out from the original puzzle. First note that the dashed lines are called “extension lines” in the paper. They are generated to allow the extraction of the sub-puzzles. To create the extension lines the algorithm goes along the outline of the puzzle and at convex corners it generates an “extension line” that allows the possible separation of the puzzle into sub-puzzles. The following figure illustrates this. In the paper they outline various ways for eliminating extension lines that are unlikely to be of use.

Extension Lines

Once the above extension lines are created there is a set of heuristics for splitting out sub-puzzles, the heuristics are used in the following order.

1. direct-match rule

The algorithm attempts to locate puzzle pieces fully described by edges, rather than by extension lines or by combinations of edges and extension lines.

2. 2 1/2 – 3 1/2 edge-match rule
The direct-match rule while being very reliable in finding matches cannot be used in most cases this means that a more relaxed rule is needed. The 2 1/2 – 3 1/2 edge-match rule allows the location of puzzle pieces where most of the periphery is described by edges and a small portion is described by extension lines. Specifically the “2 1/2 – 3 1/2 edge-match rules requires that in the case of triangular shapes two complete sides must be defined by edges and the remaining side can be defined by a combination of collinear edges and extension lines. Moreover, the combination must include at least one portion of edge. For four sided puzzle pieces the rule requires that the additional side also be fully described by an edge.” [1]

There are an additional 8 extraction rules that are used but the above two give you their general flavor.

The organization of the extraction rules is to run them recursively and at each step try to apply the most specific one possible first. If the algorithm gets to a step where there are no possible extractions it backtracks and tries another extraction rule. This process is not guaranteed to find a solution and in fact the paper gives an example of a puzzle on page 239 for which no solution is found.

The paper is a niece example of where a little clever thinking and some very specific heuristics can solve tangram puzzles.


[1] Deutsch, E. S., and K. C. Hayes Jr. “A Heuristic Solution to the Tangram Puzzle”, Machine Intelligence vol. 7, 1972, p. 205–240.

Convex Hulls

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

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.