% ************************* CONVEX HULL -- Used as a subroutine for Delaunay ************************** % % Used to find the distance of a point (o) from a line (line). % function cross_product(o,line) = let (i,xo,yo) = o; ((j,x1,y1),(k,x2,y2)) = line; in (x1-xo)*(y2-yo) - (y1-yo)*(x2-xo) \$ % Given two points on the convex hull (p1 and p2), hull finds all the points on the hull between p1 and p2 (clockwise), inclusive of p1 but not of p2. % function hull(points,p1,p2) = let cross = {cross_product(p,(p1,p2)): p in points}; packed = {p in points; c in cross | plusp(c)}; in if (#packed < 2) then [first(p1)] ++ {i: i,x,y in packed} else let pm = points[max_index(cross)]; in flatten({hull(packed,p1,p2): p1 in [p1,pm]; p2 in [pm,p2]}) \$ function convex_hull(points) = let x = {x : (i,x,y) in points}; minx = points[min_index(x)]; maxx = points[max_index(x)]; in hull(points,minx,maxx) ++ hull(points,maxx,minx) \$ function lower_hull(points) = let x = {x : (i,x,y) in points}; minx = points[min_index(x)]; maxx = points[max_index(x)]; in hull(points,maxx,minx) ++ [first(minx)] \$ % ************************* N^2 DELAUNAY ROUTINE --- used when groups are small ************************** % FUNCTION perp_line_from_point(x0,y0) = let sq = x0^2+y0^2 in x0/sq,y0/sq \$ FUNCTION delaunay_neighbors((i,(x0,y0)),pts) = let % Generate half planes from points % npts = {j,perp_line_from_point(x-x0,y-y0) : j,(x,y) in pts}; % Use convex hull to find intersection % hull = convex_hull([(i,(0.,0.))]++npts); % remove self from the hull, if it is in there % in {i,j: j in hull | j /= i} \$ FUNCTION nested_get(a,i) = let lens = {#i:i}; vals = a->flatten(i) in partition(vals,lens) \$ % Given a set of points and a set of edges (pairs of indices to the points), this returns for each point its delaunay edges sorted clockwise. It assumes that all Delaunay edges are included in the input edgelist but that they don't have to appear in both directions % FUNCTION delaunay_from_edgelist(points,edges) = let % orient the edges and remove duplicates % edges = remove_duplicates({max(i,j),min(i,j): (i,j) in edges}); % put in the back edges % edges = edges ++ {j,i : (i,j) in edges}; % create an adjacency list for each node % adj_lists = {e : i,e in int_collect(edges)}; % tag the points with indices % pts = zip([0:#points],points); % for each point subselect the delaunay edges and sort clockwise % adj_lists = {delaunay_neighbors(pt,npts): pt in pts; npts in nested_get(pts,adj_lists)}; in adj_lists \$ function slow_delaunay(pts) = if #pts == 0 then [] (int,int) else let rest = drop(pts,1); in delaunay_neighbors(pts[0],rest)++slow_delaunay(rest) \$ % ************************* THE FULL DELAUNAY ROUTINE ************************** % block_size = 64; % This is a divide-and-conquer algorithm for breaking up a set of points into groups such that points on the edges of groups are Delaunay edges % function delaunay_divide(points,previous_n) = if (#points <= block_size) or #points == previous_n % terminate if either points is smaller than block_size or if no progress was made in the previous step % then [points] else let n = #points; % flip x and y coordinates -- this makes it so that we alternate between cutting in x and in y % points = {i,y,x : i,x,y in points}; % find x median % med = median({x : i,x,y in points}); (i,xm,ym) = {i,x,y in points | x == med}[0]; % project points onto a parabola around median point % proj = {j,y,(x-xm)^2+(y-ym)^2: i,x,y in points; j in index(n)}; % find the lower hull of this parabola and mark these points % lower_hull_indices = lower_hull(proj); hull_flags = dist(f,n)<-{i,t: i in lower_hull_indices}; % divide points into two sets based on median and such that the hull points belong to both sets % down_points = {i,x,y in points; fl in hull_flags | x < med or fl}; up_points = {i,x,y in points; fl in hull_flags | x >= med or fl}; % Recurse % in flatten({delaunay_divide(x,n) : x in [down_points,up_points]}) \$ function delaunay(pts) = let % Tag the points with an index % ipts = {i,p : i in index(#pts) ; p in pts}; % Break into components using divide and conquer until point sets are of size block_size % point_groups = delaunay_divide(ipts,#ipts+1); % Now find delaunay edges within each block % all_edges = flatten({slow_delaunay(group) : group in point_groups}); % Finally put all blocks together and remove redundant edges % edges = delaunay_from_edgelist(pts,all_edges); in edges \$ % ************************* TEST ROUTINES ************************** % function random_points_in_a_square(n) = {rand(j),rand(j) : j in dist(1.0,n)} \$ function test_delaunay(n) = let pts = random_points_in_a_square(n); edges = flatten(delaunay(pts)); in #edges \$ function test_delaunay_and_plot_results(n) = let pts = random_points_in_a_square(n); edges = flatten(delaunay(pts)); window = w_make_window(((200,100),(600,600)), % Window position and size % "delaunay graph", % Window title % w_black, % Background color % display); % Default display % bound_box = w_bounding_box(pts); window,box = w_add_box(((50,50),(500,500)),% offset,size within window % bound_box, % offset,size of virtual coords% "graph box", % box title % w_white, % box color % window); % window % e1,e2 = unzip(edges); lines = zip(pts->e1,pts->e2); _ = w_draw_segments(lines,1,w_red,box); _ = w_draw_big_points(pts,2,w_blue,box); _ = w_add_text((250,20),"Click on window to remove",w_white,window); _ = w_get_input(window); % wait for an event % _ = w_kill_window(window); in #edges \$