diff options
Diffstat (limited to 'hulls/graham.hs')
| -rw-r--r-- | hulls/graham.hs | 268 |
1 files changed, 268 insertions, 0 deletions
diff --git a/hulls/graham.hs b/hulls/graham.hs new file mode 100644 index 0000000..6fc1bdc --- /dev/null +++ b/hulls/graham.hs @@ -0,0 +1,268 @@ +--------------------------------------------------------------------------------------------- +-- +-- 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 <m.i@gmx.at> +-- 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 $ "<pre>" + 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 $ "</pre>" + + 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++"</"++tag++">" + 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;"] |