% *********************************** CONVEX HULL -- The following three routines are for finding a convex hull *********************************** % % Find the cross product a point with a line % FUNCTION cross_product(o,line) = let (xo,yo) = o; ((x1,y1),(x2,y2)) = line; in (x1-xo)*(y2-yo) - (y1-yo)*(x2-xo) $ % Assuming that p1 and p2 are on the hull this routine returns the clockwise part of a convex hull from p1 to p2. It uses the quickhull algorithm % FUNCTION sub_hull(points,p1,p2) = let cross = {cross_product(p,(rest(p1),rest(p2))): (i,p) in points}; packed = {p in points; c in cross | plusp(c)}; in if (#packed < 2) then [first(p1)] ++ {i: i,pt in packed} else let pm = points[max_index(cross)]; in flatten({sub_hull(packed,p1,p2): p1 in [p1,pm]; p2 in [pm,p2]}) $ % Returns the convex hull of a set of points. In particular it returns a clockwise sorted sequence of the indices of the points on the hull using the quickhull algorithm% FUNCTION convex_hull(points) = let x = {x : i,(x,y) in points}; minx = points[min_index(x)]; maxx = points[max_index(x)]; % Find the bottom and top subhulls and append them % in sub_hull(points,minx,maxx) ++ sub_hull(points,maxx,minx) $ % *********************************** SOME UTILITIES FOR LINES *********************************** % % Given a point x0,y0 this figures out the (a,b) in the equation of (a x + b y = 1) for the line that runs through (x0,y0) and is perpendicular to the line from (0,0) to (x0,y0). % FUNCTION perp_line_from_point(x0,y0) = let sq = x0^2+y0^2 in x0/sq,y0/sq $ % Find the intersection of two lines described by the equations (a1 x + b1 y = 1) and (a2 x + b2 y = 1). % FUNCTION line_intersect((a1,b1),(a2,b2)) = let denom = b1*a2 - b2*a1; x = (b1 - b2)/denom; y = (a2 - a1)/denom; in x,y $ % *********************************** SUBSELECTING DELAUNAY EDGES *********************************** % % For a particular Delaunay point this routine subselects from a set of other points which ones are actually neighbors in the Delaunay diagram. It does this by (1) converting the points into lines that bisect the edge from the center point to the points (i.e. potential voronoi edges). These lines represent a set of half spaces. (2) Finding the intersection of the half spaces of these lines, these form the voronoi cell of the point. This is done using the standard dual of solving plane intersection with convex hull. This technique will return only the voronoi edges on the edge of the insersection. These in turn correspond to the delaunay edges. The inputs are (i,(x0,y0)) -- the index of the point along with its location pts --- the potential delaunay neighbors represented as a sequence of (index,location) pairs (as above) Ouput The indices of the Delaunay edges sorted clockwise % 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 {j in hull | j /= i} $ % A utility function % 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 $ % *********************************** RETURNING THE DELAUNAY TRIANGLES *********************************** % FUNCTION cell_triangles(i,neighbors,neighbor_neighbors) = {i,j,k : j in rotate(neighbors,-1) ; k in neighbors ; kn in neighbor_neighbors | -1 /= find(j,kn) and i < j and i < k} $ FUNCTION delaunay_triangles_from_edges(adj_list) = flatten({cell_triangles(i,cell,cc) : i in [0:#adj_list] ; cell in adj_list ; cc in nested_get(adj_list,adj_list)}) $ % *********************************** RECOGNIZING BOUNDARY EDGES *********************************** % % For a given point it checks to see which of the edges from the point to neighbors are on the boundary of the region (only have a triangle on one side). The routine assumes the edges are sorted around the point. For each point neighbor point e_i it checks if e_i s adjacency list contains e_{i-1} and e_{i+1}. n -- the neighbors nn -- the neighbors' neighbors The function returns the neighbors along a boundary, if any. % FUNCTION cell_boundary_edges(n,nn) = let nleft = rotate(n,1); c = {i in n ; l in nleft; nb in nn | not(any({l==nb:nb}))}; in c $ FUNCTION boundary_edges(points,edges) = let % return Delaunay edges in sorted order clockwise around each point % adj_list = delaunay_from_edgelist(points,edges); % for each point return its left boundary edge(s) if any % boundaries = {cell_boundary_edges(n,nn) : n in adj_list; nn in nested_get(adj_list,adj_list)}; edges = {{i,j : j in boundary}: boundary in boundaries; i in [0:#points]} % flatten into a set of directed edges % in flatten(edges) $ % *********************************** SORTING BOUNDARY EDGES INTO CYCLES *********************************** % % For a set of linked cycles represented by the pointers pt this finds the maximum val in each cycle and returns it to all elements in the cycle % FUNCTION cycle_max(pt,val) = let nval = {max(v,nv): v in val; nv in val->pt}; in if eql(nval,val) then val else cycle_max(pt->pt,nval) $ % Simple list ranking routine based on Wylies algorithm% FUNCTION list_rank_rec(pt,val) = let npt = pt->pt; in if eql(npt,pt) then val else list_rank_rec(npt,{v + nv: v in val; nv in val->pt}) $ FUNCTION list_rank(pt) = list_rank_rec(pt,{if i==pt then 0 else 1: pt; i in [0:#pt]}) $ % This routine returns a nested sequence each of which is a list of nodes along a boundary. The list is sorted by order along the boundary going to the left when looking from the inside of a triangulated region to the outside. The cycle is broken at the node with the highest node number. % FUNCTION list_boundaries(points,edges) = let node_nums, next = unzip(boundary_edges(points,edges)); % Link the boundary nodes into cyclic linked lists (cycles) % next = (dist(0,#points)<-zip(node_nums,[0:#node_nums]))->next; % Find the maximum index in each boundary cycle using pointer jumping % maxes = cycle_max(next,node_nums); % Break the cycles at the max index to form a terminated linked list % next = {if (m == i) then j else nxt : m in maxes ; i in node_nums ; j in [0:#maxes]; nxt in next}; % Get each element's rank within their list % ranks = list_rank(next); % Collect the lists into their respective cycles % cycles = collect(zip3(maxes,ranks,node_nums)); % Sort the cycles by the rank of the elements % sorted_cycles = {permute({y: x,y in vals},rank({x: x,y in vals})) : (z,vals) in cycles}; in sorted_cycles $ % *********************************** FINDING THE VORONOI EDGES FROM THE DELAUNAY NEIGHBORS *********************************** % FUNCTION voronoi_point(pt,line,nline,i,neighbors) = let flag = minusp(find(i,neighbors)); vpt = if flag then pt else line_intersect(line,nline); in flag,vpt $ FUNCTION voronoi_points(n,(x0,y0),pts,lines,neighbors,neighbor_neighbors) = let rlines = rotate(lines,n); rneighbors = rotate(neighbors,n); vpt = {voronoi_point(pt,line,rline,rn,nn): pt in pts; line in lines; rline in rlines; rn in rneighbors; nn in neighbor_neighbors}; vpt = {fl,(x0+x,y0+y): (fl,(x,y)) in vpt}; in vpt $ % p0 -- the location of the point of the voronoi cell pts -- the location of the neighbors neighbors -- the indices of the neighbors neighbor_neighbors -- the indices of the neighbors neigbors -- (this will be a nested sequence) % FUNCTION voronoi_edges_for_cell(p0,pts,neighbors,neighbor_neighbors) = let (x0,y0) = p0; rel_pts = {((x-x0)/2.,(y-y0)/2.) : (x,y) in pts}; lines = {perp_line_from_point(pt) : pt in rel_pts}; lpt = voronoi_points(1,p0,rel_pts,lines,neighbors,neighbor_neighbors); rpt = voronoi_points(-1,p0,rel_pts,lines,neighbors,neighbor_neighbors); in {(n,fl1 or fl2),(pt1,pt2): (fl1,pt1) in lpt; (fl2,pt2) in rpt ;n in neighbors} $ % The Delaunay edges neet to be in sorted order around the points % FUNCTION voronoi_edges_from_delaunay(points,delaunay_edges) = let % generate voronoi edges for each cell % vedges = {voronoi_edges_for_cell(p0,pts,n,nn) : p0 in points ; pts in nested_get(points,delaunay_edges) ; n in delaunay_edges ; nn in nested_get(delaunay_edges,delaunay_edges)}; in vedges $