%----------------------------------------------- Implemented by Kwok-Shing Leung (ksleung@cs.cmu.edu), 1995 This is a O(n lg n) work O(lg^2 n) time convex hull algorithm based on Overmar's merging technique. -------------------------------------------------- % % ----------------------------------- % % overmar_cross_product(v0, v1): % % % % Cross product between v0 and v1. % % ----------------------------------- % function overmar_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 }; % --------------------------------------------------------------------- % % overmar_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 overmar_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 = overmar_cross_product(v_front, v_21); cross_back_21 = overmar_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); % --------------------------------------------------------------------- % % overmar_intersection(LL0, LL1, LR0, LR1): % % % % Determine the x coordinate of the intersecting point of the two % % tangents. % % --------------------------------------------------------------------- % function overmar_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; % --------------------------------------------------------------------- % % overmar_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 overmar_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); % --------------------------------------------------------------------- % % overmar_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 overmar_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); % --------------------------------------------------------------------- % % overmar_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 overmar_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]; % --------------------------------------------------------------------- % % overmar_binary(H1, H2, L1, R1, L2, R2, line, off): % % % % Given two hulls and their leftmost and rightmost indices, % % overmar_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 overmar_position() to find to cutting configuration, % % and use the information to determine the next candidate sets. % % --------------------------------------------------------------------- % function overmar_binary(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 overmar_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 = overmar_right(H1[M1_front], H1M1, H1[M1_back]); H1M1_back = overmar_left(H1[M1_back], H1M1, H1[M1_front]); H2M2_front = overmar_left(H2[M2_front], H2M2, H2[M2_back]); H2M2_back = overmar_right(H2[M2_back], H2M2, H2[M2_front]); ans = 3 * overmar_position(H1M1, H2M2, H1M1_front, H1M1_back) + overmar_position(H1M1, H2M2, H2M2_front, H2M2_back); in if ans == 0 then (M1, M2) else if ans == 4 then (if (overmar_intersection(H1[M1], H1[M1_front], H2[M2], H2[M2_front]) < line) then overmar_binary(H1, H2, M1, R1, L2, R2, line, off) else overmar_binary(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 overmar_binary(H1, H2, L1_new, R1_new, L2_new, R2_new, line, off); % --------------------------------------------------------------------- % % overmar_help_gather(points, starting, ending): % % % % Gather all the points from starting to ending in a clockwise fashion. % % --------------------------------------------------------------------- % function overmar_help_gather(points, start, end) = if (start <= end) then subseq(points, start, end + 1) else let size = #points; in subseq(points, start, size) ++ subseq(points, 0, end + 1); % --------------------------------------------------------------------- % % overmar_help(points): % % % % overmar_help(points) computes the convex hull on a set of x-sorted % % list of points. It first divides all the points into two halves of % % equal sizes, and then call overmar_help() on each on them. Then, it % % calls overmar_binary() to find the upper and lower tangents, and the % % results are sewed together. % % --------------------------------------------------------------------- % function overmar_help(points) = if (#points < 7) then points->quick_convex_hull(points) else let left, right = let p = bottop(points); in p[0], p[1]; (left_x, left_y), (right_x, right_y) = unzip(left), unzip(right); left_hull, right_hull = overmar_help(left), overmar_help(right); left_count, right_count = #left_hull, #right_hull; (left_hull_x, left_hull_y) = unzip(left_hull); (right_hull_x, right_hull_y) = unzip(right_hull); L1, R1 = min_index(left_hull_x), max_index(left_hull_x); L2, R2 = min_index(right_hull_x), max_index(right_hull_x); center_line = (left_x[#left_x - 1] + right_x[0]) / 2.0; (M1_top, M2_top) = overmar_binary(left_hull, right_hull, L1, R1, L2, R2, center_line, 0); inv_left_hull = { (x, -y) : (x, y) in left_hull } -> { left_count - 1 - i: i in [0:left_count] }; inv_right_hull = { (x, -y) : (x, y) in right_hull } -> { right_count - 1 - i: i in [0:right_count] }; inv_L1, inv_R1 = left_count - 1 - L1, left_count - 1 - R1; inv_L2, inv_R2 = right_count - 1 - L2, right_count - 1 - R2; inv_center_line = center_line; (inv_M1_top, inv_M2_top) = overmar_binary(inv_left_hull, inv_right_hull, inv_L1, inv_R1, inv_L2, inv_R2, inv_center_line, 0); (M1_bot, M2_bot) = left_count - 1 - inv_M1_top, right_count - 1 - inv_M2_top; in overmar_help_gather(left_hull, M1_bot, M1_top) ++ overmar_help_gather(right_hull, M2_top, M2_bot); % ===================================================================== % % overmar_convex_hull(points): % % % % overmar_convex_hull() takes a point set, sorts the points, and passes % % the point set to overmar_help() to find the convex hull using % % Overmars' algorithm. % % ===================================================================== % function overmar_convex_hull(points) = let points_x, points_y = unzip(points); rank_x = rank(points_x); in overmar_help(permute(points, rank_x)); % -------------------------------- % % overmar_random_points(n): % % % % Generate a set of random points. % % -------------------------------- % function overmar_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} ; % ------------------------------------------------------- % % overmar_draw_dots(points, window, r): % % % % Draw a point (as a rectangle for better visual effect). % % ------------------------------------------------------- % function overmar_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 overmar_draw_dots(subseq(points, 1, #points), window, r); % ===================================================================== % % overmar_demo(n): % % % % overmar_demo() randomly generates a set of points in the Euclidean % % plane, finds the convex hull, adn then displays all the points and % % the convex hull. % % ===================================================================== % function overmar_demo(n) = let points = overmar_random_points(n); 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 = overmar_draw_dots(points, window, 0.01); hull = overmar_convex_hull(points); ignore = draw_lines(hull ++ take(hull, 1), 1, window); in close_window(window);