path: root/hulls/graham.hs

---------------------------------------------------------------------------------------------
--
-- 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;"]