% This example finds the convex hull in 2 dimensions for a set of points using the QuickHull algorith. The algorithm is described in the NESL language definition. This file also includes a function, hull_demo(n), that picks approximately n random points in a circle, plots the points, finds the hull, and plots the hull. % % Used to find the distance of a point (o) from a line (line). % function cross_product(o,line) = let (xo,yo) = o; ((x1,y1),(x2,y2)) = line; in (x1-xo)*(y2-yo) - (y1-yo)*(x2-xo); % Given two points on the convex hull (p1 and p2), hsplit finds all the points on the hull between p1 and p2 (clockwise), inclusive of p1 but not of p2. % function hsplit(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 [p1] ++ packed else let pm = points[max_index(cross)]; in flatten({hsplit(packed,p1,p2): p1 in [p1,pm]; p2 in [pm,p2]}); % Finds the points with minimum and maximum x coordinates, and then finds the upper and lower convex hull: the part clockwise from minx to maxx (upper) and clockwise from maxx to minx (lower). % function convex_hull(points) = let x = {x : (x,y) in points}; minx = points[min_index(x)]; maxx = points[max_index(x)]; in hsplit(points,minx,maxx) ++ hsplit(points,maxx,minx); % Generates approximately n random points in the unit circle. It works by generating 4n/pi points in a square and throwing out the points that are not within the circle. % function 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} ; % Demo of convex hull. N is the number approximate number of points. % function hull_demo(n) = let points = random_points(n); ignore = if (#points < 2) then print_string("Not enough points.") else f; box = bounding_box([(-1.0,-1.0),(1.0,1.0)]); window = w_make_window(((200,100),(600,600)), % Window position and size % "convex hull", % Window title % w_black, % Background color % display); % Default display % window,box = w_add_box(((50,50),(500,500)),% offset,size within window % ((-1.,-1.),(2.,2.)),% offset,size of virtual coords% "hull box", % box title % w_white, % box color % window); % window % ignore = w_draw_points(points,w_black,box); hull = convex_hull(points); ignore = w_draw_segments(zip(hull,rotate(hull,1)),2,w_black,box); ignore = w_add_text((250,20),"Click on window to remove",w_white,window); ignore = w_get_input(window); % wait for an event % in w_kill_window(window)$