--------------------------------------------------------------------------------------------- -- -- GRAHAM SCAN AND ALPHA SHAPES IMPLEMENTATIONS -- -- This little haskell programm calulates the Convex Hull for a set of 2D points. -- It ships wit a main function that feeds the Graham Scan algorithm with some -- random points and generates a simple SVG of the input points and resulting envelope. -- A simple SVG encoder is included. -- -- Alogrithm used: https://en.wikipedia.org/wiki/Graham_scan -- -- CREDITS -- -- -- Michal Idziorek -- 16 December 2017 -- --------------------------------------------------------------------------------------------- -- Added triangulation -- may 4th/6th 2018/ -- https://en.wikipedia.org/wiki/Point_set_triangulation -- ------------------------------------------------------------ --todo: better algos/ datastructures/ benchmark/optimize cleanup/modularize --voronoi/alphashapes/real alpha shapes with bended edges.. --use some lattice, especially for findNonDel! import Data.List import Data.List.Split import Data.Maybe import System.Random import Debug.Trace import System.IO import System.Environment -- SETTINGS -- --point_count=50 rnd_min=0 rnd_max=200 svg_w=500 svg_h=500 --------------------------------------------------------------------------------------------- -- GRAHAM SCAN -- -- Three points are clockwise if ccw < 0. ccw (p1x,p1y) (p2x,p2y) (p3x,p3y) = (p2x - p1x)*(p3y - p1y) - (p2y - p1y)*(p3x - p1x) -- Calculate the slope defined by 2 points. (return Infinity, if points are identical). slope (ax,ay) (bx,by) | (ax,ay ) == (bx,by) = 1/0 -- Infinity | otherwise = (bx-ax)/(by-ay) -- Comparison function to sort points counterclockwise (given a reference point). slope_cmp a b c = compare (slope a c) (slope a b) -- Comparison function using the y and x coordinates for ordering. graham_cmp (ax,ay) (bx,by) | ay /= by = compare ay by | otherwise = compare ax bx -- Graham scan on prepared data. this will calculate the convex hull. graham_calc [] hs = hs graham_calc (x1:xs) hh | length(hh) < 2 = graham_calc xs (x1:hh) graham_calc xx@(x1:xs) hh@(h1:h2:hs) | ccw x1 h1 h2 < 0 = graham_calc xs (x1:hh) | otherwise = graham_calc xx (h2:hs) -- Find the starting point, sort all points counterclockwise and perform the graham scan. graham xs = graham_calc sortedPoints [] where minPoint = minimumBy graham_cmp xs sortedPoints = sortBy (slope_cmp minPoint) xs -- -- ---- draws SVG points and lines (in hardcoded sizes and colors) -- -- --svg_draw w h p l = xml_enc_shit "svg" [style,["width",show w],["height",show h]] body -- -- -- where style = ["style", -- -- -- "background-color:black;border:3px solid green;margin:2px;"] -- -- -- body = (unlines (map point_to_svg_1 p )) ++ -- -- -- (unlines (map (line_to_svg (255,0,0)) l)) -- -- -- -- -- ---- calculate convex hull and generate svg -- -- --svg_graham xs = svg_draw 600 600 xs (zip hull hull_open) -- -- -- where hull_open = graham xs -- -- -- hull = (last hull_open) : hull_open -- TRIANGULATION -- flipper c tri = case findNonDel tri of Nothing -> (tri,c) Just (t1,t2,f1,f2) -> flipper (c+1) $ f1 : f2 : (delete t1 $ delete t2 tri) --todo check angles findNonDel tri = head' nondel where nondel= [(a,b,fst (flipTri a b),snd (flipTri a b)) | a<-tri , b<-tri, a/=b, angleBig a b] angleBig t1 t2 = case commonEdge t1 t2 of Just com -> an t1 com + an t2 com > pi _ -> False an (a,b,c) edge = vangle pt0 pt1 pt2 where (pt1:pt2:[])= filter (`elem` [fst edge,snd edge]) [a,b,c] (pt0:[])= filter (not.(`elem` [fst edge,snd edge])) [a,b,c] vangle v0 v1 v2 =acos $ vdot a b / (vmag a * vmag b) where a=vsub v1 v0 b=vsub v2 v0 vsub (x1,y1) (x2,y2) = (x1-x2,y1-y2) vmag (x,y) = sqrt (x*x+y*y) vdot (ax,ay) (bx,by) = ax*bx+ay*by flipTri t1@(a,b,c) t2@(d,e,f) = (t1',t2') where (t1',t2')=((v3,v4,v1),(v3,v4,v2)) (v1,v2) =fromJust $ commonEdge t1 t2 v3=head $ filter (not.(`elem` [v1,v2])) [a,b,c] v4=head $ filter (not.(`elem` [v1,v2])) [d,e,f] commonEdge t1@(a,b,c) t2@(d,e,f) = head' [l1 | l1@(a1,a2)<-edges, l2@(b1,b2)<-delete l1 edges, (a1==b1&&a2==b2)||(a1==b2&&a2==b1)] where edges = [(a,b), (b,c), (c,a), (d,e), (e,f), (f,d)] -- http://www.geom.uiuc.edu/~samuelp/del_project.html -- http://web.colby.edu/thegeometricviewpoint/2014/11/23/algorithms-for-working-with-triangulated-surfaces/ -- DIVIDE and CONQUER -- helping head head' (x:xs) = Just x head' _ = Nothing -- trivial triangulation stuff -- triangulate xs =foldl (f p) [] (zip pts$ tail pts) where (p:pts)=xs f s a (x,y) = (s,x,y):a subtriangulate tri xs = foldl f tri xs where f a x = case findInside a x of Just t@(t1,t2,t3) -> (x,t1,t2):(x,t1,t3):(x,t2,t3):delete t a Nothing -> error ("fail on subtriang " ++ show x) undefined -- (x,(0,0),(0,0)): -- TRIANGLE INTERIOR TEST -- -- https://stackoverflow.com/questions/2049582/how-to-determine-if-a-point-is-in-a-2d-triangle findInside tri p = foldl f Nothing tri where f Nothing t = if isInside t p then Just t else Nothing f x _ = x sign (p1x,p1y) (p2x,p2y) (p3x,p3y)= (p1x - p3x) * (p2y - p3y) - (p2x - p3x) * (p1y - p3y) isInside (v1,v2,v3) pt = b1==b2 && b2==b3 where b1 = sign pt v1 v2 < 0 b2 = sign pt v2 v3 < 0 b3 = sign pt v3 v1 < 0 -- ALPHA COMPLEX -- -- http://codeforces.com/blog/entry/17313 -- https://en.wikipedia.org/wiki/Circumscribed_circle#Cartesian_coordinates -- http://www.ambrsoft.com/trigocalc/circle3d.htm -- center_from (bx,by) (cx,cy)=((cy*_B-by*_C)/(2*_D), (bx*_C-cx*_B)/(2*_D)) where _B=bx*bx+by*by _C=cx*cx+cy*cy _D=bx*cy-by*cx circle_from ((ax,ay),(bx,by),(cx,cy)) = (radius,(ax+x,ay+y)) where radius =sqrt((x-(bx-ax))**2+(y-(by-ay))**2) (x,y) = center_from (bx-ax,by-ay) (cx-ax, cy-ay) radius_from = fst.circle_from -- rturn circle(I+A, abs(I)); -- where xc = ((ax * ax + ay * ay) * (by - cy) + (bx * bx + by * by) * (cy - ay) + (cx * cx + cy * cy) * (ay - by)) / d -- yc = ((ax * ax + ay * ay) * (cx - bx) + (bx * bx + by * by) * (ax - cx) + (cx * cx + cy * cy) * (bx - ax)) / d -- RANDOMIZING -- randomPoints g cnt = take cnt (zip r10a r10b) where r5 = randomRs (rnd_min,rnd_max) g :: [Double] r10a =zipWith (+) r5 (drop cnt r5) r10b =zipWith (+) (drop (2*cnt) r5) (drop (3*cnt) r5) --------------------------------------------------------------------------------------------- -- MAIN -- -- let svgData = svgGenerate svgExample1 main = do g <- newStdGen args <- getArgs let point_count = read $ args!!0 let points2=randomPoints g 99999 dat<-readFile "data.txt" outp<-openFile "nice.html" WriteMode let points = [(10,10),(100,100),(100,10)] let points = randomPoints g point_count let rnds = chunksOf 2 $ (randomRs (-1,1) g :: [Double]) let points = map (\((a:b:_),(x:y:_))->(read x*10+a,read y*10+b)) $zip rnds $ map (splitOn ",") $ lines dat print points let cvx_hull=graham points let cvx_hull_tri=triangulate cvx_hull let points_inside=filter (not.(`elem` graham points)) points let full_tri=subtriangulate cvx_hull_tri points_inside let (delaunay,delaunay_count)=flipper 0 full_tri let delaunay_circles=map circle_from delaunay let alfa_tri alf=filter ((<1/alf).radius_from) delaunay let putOutput = hPutStr outp putOutput $ "
"
           putOutput $ "Input data set: " ++ show (length points) ++ " points.\n"
           putOutput $ "Calculated convex hull having: " ++ show (length cvx_hull) ++ " points.\n"
           putOutput $ "Convex hull traingulated to: " ++ show (length cvx_hull_tri) ++ " triangles.\n"
           putOutput $ "Fully traingulated to: " ++ show (length full_tri) ++ " triangles."
           putOutput $ " (expected count: " ++ show (length points * 2 - 2 - length cvx_hull) ++ " according to: 2*n-h-2)\n"
           putOutput $ "Converted to delaunay triangulation after: " ++ show delaunay_count ++ " flips.\n"
           putOutput $ "
" putOutput $ svg_draw_body svg_w svg_h $ (concat $ map (tri_to_svg (100,100,255)) cvx_hull_tri) ++ (concat $ map (point_to_svg (255,100,100)) points) putOutput $ svg_draw_body svg_w svg_h $ (concat $ map (tri_to_svg (100,100,255)) full_tri) ++ (concat $ map (point_to_svg (255,100,100)) points) putOutput $ svg_draw_body svg_w svg_h $ (concat $ map (tri_to_svg (100,100,255)) delaunay) ++ (concat $ map (point_to_svg (255,100,100)) points) putOutput $ svg_draw_body svg_w svg_h $ (concat $ map (tri_to_svg (100,100,255)) delaunay) ++ (concat $ map (point_to_svg (255,100,100)) points) ++ (concat $ map (circle_to_svg (255,100,100)) delaunay_circles) let getAlpha alpha= svg_draw_body svg_w svg_h $ (concat $ map (tri_to_svg (100,100,100)) delaunay) ++ (concat $ map (tri_fill_svg (100,100,0)) (alfa_tri alpha)) putOutput $ concat $ map getAlpha [0.1,0.2..1.7] -- SOME BASIC SVG ENCODING FUNCS -- xml_enc_shit tag attrs body = "<"++tag++" "++xml_attrs++">"++body++"" where xml_attrs = unlines $ map xml_attr attrs xml_attr (x:xs) = x++"=\""++(head xs)++"\" " line_to_svg (r,g,b) ((x1,y1),(x2,y2)) = xml_enc_shit "line" [["x1",show lx1],["y1",show ly1], ["x2",show lx2],["y2",show ly2], ["style", "stroke:rgb("++show r++","++show g++","++show b++");stroke-width:2"]] "" where lx1=x1 lx2=x2 ly1=y1 ly2=y2 tri_to_svg col (a,b,c) = line_to_svg col (a,b)++line_to_svg col (b,c)++line_to_svg col (c,a) tri_fill_svg col ((x1,y1),(x2,y2),(x3,y3)) = xml_enc_shit "polygon" [["points",pts],["style","fill:rgb(100,100,200)"]] "" where pts=show (x1)++","++show (y1)++" "++show (x2)++","++show y2++" "++show (x3)++","++show (y3) point_to_svg_1 (x,y) = xml_enc_shit "circle" [["cx",show cx],["cy",show cy],["r","5"], ["fill","rgb(30,150,"++(show (floor dist))++")"]] "" where cx=x cy=y dist= (sqrt ((x-5)*(x-5) + (y-5)*(y-5)))*255/8 point_to_svg (r,g,b) (x,y) = xml_enc_shit "circle" [["cx",show cx],["cy",show cy],["r","5"], ["fill","rgb("++col++")"]] "" where cx=x cy=y dist= (sqrt ((x-5)*(x-5) + (y-5)*(y-5)))*255/8 col=show r++","++show g++","++show b circle_to_svg (r,g,b) (radius,(x,y)) = xml_enc_shit "circle" [["cx",show cx],["cy",show cy],["r",show rad], ["stroke","rgb("++col++")"],["fill","none"]] "" where cx=x cy=y rad=radius col=show r++","++show g++","++show b -- draws SVG points and lines (in hardcoded sizes and colors) svg_draw_body w h body = xml_enc_shit "svg" [style,["width",show w],["height",show h]] body where style = ["style", "background-color:black;border:3px solid green;margin:2px;"]