%----------------------------------------------- Implemented by Kwok-Shing Leung (ksleung@cs.cmu.edu), 1995 This file contains an implementation of an algorithm described in: M. J. Atallah and M. T. Goodrich. Parallel algorithms for some functions of two convex polygons. Algorithmica, 3(4):535-548, 1988. This is a O(n lg n) work O(lg n) time convex hull algorithm based on Overmar's merging technique. -------------------------------------------------- % % -------------------------------- % % fast_cross_product(v0, v1): % % % % Cross product between v0 and v1. % % -------------------------------- % function fast_cross_product(v0, v1) = let (v0_x, v0_y) = v0; (v1_x, v1_y) = v1; in v0_x * v1_y - v0_y * v1_x; % ---------------------------- % % quick_cross_product(v0, v1): % % % % Cross product for QuickHull. % % ---------------------------- % function quick_cross_product(o, line) = let (xo, yo, zo) = o; (x1, y1, z1), (x2, y2, z2) = line; in (x1 - xo) * (y2 - yo) - (y1 - yo) * (x2 - xo); % ---------------------------------------------------------- % % quick_hsplit(points, p1, p2): % % % % Find the upper hull bounded between p1 and p2 recursively. % % ---------------------------------------------------------- % function quick_hsplit(points, p1, p2) = let cross = { quick_cross_product(points, (p1, p2)) : points }; packed = { points : points; cross | plusp(cross) }; in if (#packed < 2) then [p1] ++ packed else let pm = points[max_index(cross)]; in flatten( { quick_hsplit(packed, p1, p2) : p1 in [p1, pm]; p2 in [pm, p2]} ); % --------------------------------------------------------------------- % % quick_convex_hull(points): % % % % Find the upper tangent of the convex hull for points when at most % % four points are involved. % % --------------------------------------------------------------------- % function quick_convex_hull(points) = let new_points = { (x, y, z) : (x, y) in points; z in [0:#points] }; x = { x : (x, y) in points }; minx = new_points[min_index(x)]; maxx = new_points[max_index(x)]; result = quick_hsplit(new_points, minx, maxx) ++ quick_hsplit(new_points, maxx, minx); in { z : (x, y, z) in result }; % --------------------------------------------------------------------- % % fast_position(M1, M2, M1_front, M1_back): % % % % Given M1, M2, M1_front, and M1_back, determine in O(1) time how the % % line from M1 to M2 cuts the convex hull. The line can (1) touch % % the convex hull, cut the convex hull at the back, and or cut the % % convex hull at the front. This is assuming that the points are on % % the right upper hull. % % % % back_21: +ve +ve -ve % % % % front_21: +ve -ve +ve % % % % config: 0 1 2 % % % % / \ % % shape: -+- -+- -+- % % / \ / \ % % --------------------------------------------------------------------- % function fast_position(M1, M2, M1_front, M1_back) = let (M1_x, M1_y) = M1; (M2_x, M2_y) = M2; (M1_back_x, M1_back_y) = M1_back; (M1_front_x, M1_front_y) = M1_front; v_front = (M1_front_x - M1_x, M1_front_y - M1_y); v_back = (M1_back_x - M1_x, M1_back_y - M1_y); v_21 = (M2_x - M1_x, M2_y - M1_y); cross_front_21 = fast_cross_product(v_front, v_21); cross_back_21 = fast_cross_product(v_back, v_21); in if (cross_back_21 < 0.0) then 2 else (if (cross_front_21 < 0.0) then 1 else 0); % --------------------------------------------------------------------- % % fast_intersection(LL0, LL1, LR0, LR1): % % % % Determine the x coordinate of the intersecting point of the two % % tangents. % % --------------------------------------------------------------------- % function fast_intersection(LL0, LL1, LR0, LR1) = let (LL0_x, LL0_y) = LL0; (LL1_x, LL1_y) = LL1; (LR0_x, LR0_y) = LR0; (LR1_x, LR1_y) = LR1; nominator = (LR1_x - LR0_x) * (LL0_x * LL1_y - LL1_x * LL0_y) - (LL1_x - LL0_x) * (LR0_x * LR1_y - LR1_x * LR0_y); denominator = (LR1_x - LR0_x) * (LL1_y - LL0_y) - (LL1_x - LL0_x) * (LR1_y - LR0_y); in nominator / denominator; % --------------------------------------------------------------------- % % fast_binary_two_two(H1, H2, L1, R1, L2, R2): % % % % Find the upper tangent when the number of candidate points on each % % side is at most two. This avoids the non-convergent median halfing % % problem (that is, reducing (L, R) to (L, M) or (M, R) may not change % % the size of the candidate set, when there are two or less candidates. % % --------------------------------------------------------------------- % function fast_binary_two_two(H1, H2, L1, R1, L2, R2) = let size_left = if (L1 == R1) then 1 else 2; result = quick_convex_hull([H1[L1], H1[R1], H2[L2], H2[R2]]); result_warp = subseq(result, 1, #result) ++ [result[0]]; answer = { result, result_warp : result; result_warp | result < size_left and result_warp >= size_left }; in let L, R = answer[0]; array = [L1, R1, L2, R2]; in array[L], array[R]; % --------------------------------------------------------------------- % % fast_left(P, P_ref, P_opposite): % % % % If the "left" point of P is to the right, replace it with an % % artificial point to the left of P so as to make the algorithm work. % % --------------------------------------------------------------------- % function fast_left(P, P_ref, P_opposite) = let P_x, P_y = P; P_ref_x, P_ref_y = P_ref; P_opp_x, P_opp_y = P_opposite; in if (P_x < P_ref_x) then P else (2.0 * P_ref_x - P_opp_x, 2.0 * P_ref_y - P_opp_y - 1.0); % --------------------------------------------------------------------- % % fast_right(P, P_ref, P_opposite): % % % % If the "right" point of P is to the left, replace it with an % % artificial point to the right of P so as to make the algorithm work. % % --------------------------------------------------------------------- % function fast_right(P, P_ref, P_opposite) = let P_x, P_y = P; P_ref_x, P_ref_y = P_ref; P_opp_x, P_opp_y = P_opposite; in if (P_x > P_ref_x) then P else (2.0 * P_ref_x - P_opp_x, 2.0 * P_ref_y - P_opp_y - 1.0); % --------------------------------------------------------------------- % % fast_binary_help(H1, H2, L1, R1, L2, R2, line, off): % % % % Given two hulls and their leftmost and rightmost indices, % % fast_binary_help compute the median nodes and then performs Overmars' % % binary search procedure to find the upper tangent. First of all, the % % median is found. off is used to fix the problem that happens when % % the candidate set has two or less elements. For example, L = 0, and % % R = 1, then M = (0+1)/2 = 0. If Overmars' algorithm decides that the % % next candidate set is [M,R], it means [0,1] which is the same as the % % one before. Once the medians are found, it calls fast_position() to % % find to cutting configuration, and use the information to determine % % the next candidate sets. % % --------------------------------------------------------------------- % function fast_binary_help(H1, H2, L1, R1, L2, R2, line, off) = let D1 = if (L1 > R1) then (R1 + #H1 - L1) else R1 - L1; D2 = if (L2 > R2) then (R2 + #H2 - L2) else R2 - L2; off = 1 - off; in if (D1 < 2 and D2 < 2) then fast_binary_two_two(H1, H2, L1, R1, L2, R2) else let (S1, S2) = (#H1, #H2); M1 = if (L1 > R1) then mod((L1 + S1 + R1 + off) / 2, S1) else (L1 + R1 + off) / 2; M2 = if (L2 > R2) then mod((L2 + S2 + R2 + off) / 2, S2) else (L2 + R2 + off) / 2; (M1_front, M1_back) = mod(M1 + 1, S1), mod(M1 + S1 - 1, S1); (M2_back, M2_front) = mod(M2 + 1, S2), mod(M2 + S2 - 1, S2); H1M1, H2M2 = H1[M1], H2[M2]; H1M1_front = fast_right(H1[M1_front], H1M1, H1[M1_back]); H1M1_back = fast_left(H1[M1_back], H1M1, H1[M1_front]); H2M2_front = fast_left(H2[M2_front], H2M2, H2[M2_back]); H2M2_back = fast_right(H2[M2_back], H2M2, H2[M2_front]); ans = 3 * fast_position(H1M1, H2M2, H1M1_front, H1M1_back) + fast_position(H1M1, H2M2, H2M2_front, H2M2_back); in if ans == 0 then (M1, M2) else if ans == 4 then (if (fast_intersection(H1M1, H1M1_front, H2M2, H2M2_front) < line) then fast_binary_help(H1, H2, M1, R1, L2, R2, line, off) else fast_binary_help(H1, H2, L1, R1, L2, M2, line, off)) else let L1_new = [-1, L1, L1, M1, -1, L1, L1, L1, L1][ans]; R1_new = [-1, M1, M1, R1, -1, R1, M1, M1, M1][ans]; L2_new = [-1, L2, M2, M2, -1, M2, M2, L2, M2][ans]; R2_new = [-1, M2, R2, R2, -1, R2, R2, R2, R2][ans]; in fast_binary_help(H1, H2, L1_new, R1_new, L2_new, R2_new, line, off); % -------------------------------------------------- % % fast_angle(p1, p2): % % % % Find the (tangent of the) angle between p1 and p2. % % -------------------------------------------------- % function fast_angle(p1, p2) = let (x1, y1), (x2, y2) = p1, p2; in (y2 - y1) / (x2 - x1); % ------------------------------------------------------------------ % % fast_binary(H1, H2, L1, R1, L2, R2, line, off): % % % % Find the upper tangent of H1 and H2 by calling fast_binary_help(). % % ------------------------------------------------------------------ % function fast_binary(H1, H2, L1, R1, L2, R2, line, off) = let M1, M2 = fast_binary_help(H1, H2, L1, R1, L2, R2, line, off); in M1, M2, fast_angle(H1[M1], H2[M2]); % --------------------------------------------------------------------- % % fast_help_gather(points, starting, ending): % % % % Gather all the points from starting to ending in a clockwise fashion. % % --------------------------------------------------------------------- % function fast_help_gather(points, starting, ending) = let size = #points; in if (starting <= ending) then subseq(points, starting, ending + 1) else subseq(points, starting, size) ++ subseq(points, 0, ending + 1); % -------------------------------------- % % fast_sort(index_angle): % % % % Find the index with the maximum angle. % % -------------------------------------- % function fast_sort(index_angle) = if (#index_angle == 0) then (-1, 0.0) else let index, angle = unzip(index_angle); best_index = max_index(angle); in index_angle[best_index]; % --------------------------------------------------------------------- % % fast_valid(index_angle1, index_angle2, R, S): % % % % Check whether the upper hull is in the final hull by checking whether % % whether the sum of the tangents of the angles are less than zero. % % --------------------------------------------------------------------- % function fast_valid(index_angle1, index_angle2, R, S) = let index1, angle1 = index_angle1; index2, angle2 = index_angle2; in angle1 + angle2 < 0.0; % --------------------------------------------------------------------- % % fast_help_help(hull, L, R, S, line): % % % % Find the upper final hull, given a list of hulls, their Ls and Rs, % % their sizes, and the separating line. First of all, all tangents % % between each pair of hulls are found (in result). Then, each hull % % collects all its tangents to the left (collect_H1) and all the % % tangents to the right (collect_H2). The best left and best right % % tangents are then found (in H1_result and H2_result). Depending on % % their angles, some hulls may or may not be on the final hull, and % % this is done by calling fast_valid(), and the result is in valid. % % Finally, all the valid hulls are sewed together and returned. % % --------------------------------------------------------------------- % function fast_help_help(hull, L, R, S, line) = let size = #hull; I_H1 = flatten({ dist(i, size) : i in [0:size] }); I_H2 = flatten(let H2 = [0:size]; in dist(H2, size)); collect_H1_I = { [i:i * size:size] : i in [0:size] }; collect_H2_I = { [i * size + i + 1:i * size + size] : i in [0:size] }; flag = { I_H1 < I_H2 : I_H1; I_H2 }; result = { if (not(flag)) then (-1, -1, 0.0) else fast_binary(H1, H2, L1, R1, L2, R2, line, 0) : H1 in hull -> I_H1; H2 in hull -> I_H2; L1 in L -> I_H1; R1 in R -> I_H1; L2 in L -> I_H2; R2 in R -> I_H2; line in line -> I_H1; flag }; collect_H1 = { { p2, -angle : (p1, p2, angle) in result -> collect_H1_I } : collect_H1_I }; collect_H2 = { { p1, angle : (p1, p2, angle) in result -> collect_H2_I } : collect_H2_I }; H1_result = { fast_sort(collect_H1) : collect_H1 }; H2_result = { fast_sort(collect_H2) : collect_H2 }; valid = { fast_valid(H1_result, H2_result, R, S) : H1_result; H2_result; R; S } <- [(0, f), (size - 1, f)]; list = { hull, i, j : hull; (i, theta1) in H1_result; (j, theta2) in H2_result; valid | valid }; top_L = let (x, a) = H2_result[0]; in x; top_R = let (x, a) = H1_result[size - 1]; in x; in top_L, top_R, flatten({ fast_help_gather(list) : list }); % --------------------------------------------------------------------- % % fast_help(points): % % % % fast_help(points) computes the convex hull on a set of x-sorted list % % of points. It first divides all the points into sqrt(n) pieces of % % equal sizes, and then call fast_help() on each on them. Then % % fast_help_help() is called to find the upper and lower hulls of the % % final hull, and the result is sewed together. % % --------------------------------------------------------------------- % function fast_help(points) = if (#points == 1) then points else let total_size = #points; size_set = let size = floor(sqrt(float(total_size))); in dist(size, total_size / size) ++ let remainder = mod(total_size, size); in if (remainder == 0) then [] int else [remainder]; size = #size_set; point_set = { p : p in partition(points, size_set) }; first_x = { let (x, y) = point_set[0]; in x : point_set }; last_x = { let (x, y) = point_set[#point_set - 1]; in x : point_set }; line = { (first_x + last_x) / 2.0 : first_x in subseq(first_x, 1, #first_x); last_x in subseq(last_x, 0, #last_x - 1) } ++ [0.0]; hull = { fast_help(point_set) : point_set }; hull_x = { { x : (x, y) in hull } : hull }; L, R = { min_index(hull_x) : hull_x }, { max_index(hull_x) : hull_x }; S = { #hull_x : hull_x }; inv_hull = { { (x, -y) : (x, y) in hull } -> { S - 1 - i : i in [0:S] } : hull; S }; inv_L, inv_R = { S - 1 - L : S; L }, { S - 1 - R : S; R }; top_L, top_R, top_list = fast_help_help(hull, L, R, S, line); inv_bot_L, inv_bot_R, inv_bot_list = fast_help_help(inv_hull, inv_L, inv_R, S, line); bot_L, bot_R, bot_list = S[0] - 1 - inv_bot_L, S[size - 1] - 1 - inv_bot_R, { (x, -y) : (x, y) in reverse(inv_bot_list) }; in fast_help_gather(hull[0], bot_L, top_L) ++ top_list ++ fast_help_gather(hull[size - 1], top_R, bot_R) ++ bot_list; % ===================================================================== % % fast_convex_hull(points): % % % % fast_convex_hull() takes a point set, sorts the points, and passes % % the point set to fast_help() to find the convex hull using the fast % % variant of the Overmars' algorithm. % % ===================================================================== % function fast_convex_hull(points) = let points_x, points_y = unzip(points); rank_x = rank(points_x); in fast_help(permute(points, rank_x)); % -------------------------------- % % fast_random_points(n): % % % % Generate a set of random points. % % -------------------------------- % function fast_random_points(n) = let m = round(float(n) * 4.0 / pi); points = {rand(i) - 1.0, rand(i) - 1.0: i in dist(2.0,m)}; in {x,y in points | x^2 + y^2 < 1.0} ; % ------------------------------------------------------- % % fast_draw_dots(points, window, r): % % % % Draw a point (as a rectangle for better visual effect). % % ------------------------------------------------------- % function fast_draw_dots(points, window, r) = if #points == 0 then 0 else let (x, y) = points[0]; ignore = draw_lines([(x-r, y-r), (x-r, y+r), (x+r, y+r), (x+r, y-r), (x-r,y-r)], 1, window); in fast_draw_dots(subseq(points, 1, #points), window, r); % --------------------------------------------------------------------- % % fast_demo_help(points): % % % % fast_demo() takes the points, calls convex hull, and then display the % % points and the hull. % % --------------------------------------------------------------------- % function fast_demo_help(points) = let ignore = if (#points < 2) then error("Not enough points.") else f; box = bounding_box([(-1.0, -1.0), (1.0, 1.0)]); window = make_window(((200, 100), (600, 600)), box, "convex hull", display); ignore = fast_draw_dots(points, window, 0.01); hull = fast_convex_hull(points); ignore = draw_lines(hull ++ take(hull, 1), 1, window); in close_window(window); % ===================================================================== % % fast_demo(n): % % % % fast_demo() randomly generates a set of points in the Euclidean % % plane, and passes the point set to fast_convex_hull() to find the % % convex hull. The points as well as the convex hull are then % % displayed. % % ===================================================================== % function fast_demo(n) = fast_demo_help(fast_random_points(n));