Calculate and render light polygons
[L-seed.git] / src / Lseed / Geometry.hs
1 {-# LANGUAGE ScopedTypeVariables #-}
2 module Lseed.Geometry where
3
4 import Lseed.Data
5 import Lseed.Constants
6 import Lseed.Geometry.Generator
7 import Data.List
8 import Data.Maybe
9
10 import Debug.Trace
11
12 type Point = (Double, Double)
13 type Line  = (Point, Point)
14
15 lightFalloff = 0.7
16
17 eps = 1e-9
18
19 lineLength ((x1,y1),(x2,y2)) = sqrt (x1*x2 + y1*y2)
20
21 -- | from http://www.pdas.com/lineint.htm
22 crossPoint :: Line -> Line -> Maybe Point
23 crossPoint ((x1,y1),(x2,y2)) ((x3,y3),(x4,y4)) =
24         let a1 = y2-y1
25             b1 = x1-x2
26             c1 = x2*y1 - x1*y2  -- { a1*x + b1*y + c1 = 0 is line 1 }
27             a2 = y4-y3
28             b2 = x3-x4
29             c2 = x4*y3 - x3*y4  -- { a2*x + b2*y + c2 = 0 is line 2 }
30             denom = a1*b2 - a2*b1
31         in if abs denom > eps
32            then let x = (b1*c2 - b2*c1)/denom
33                     y = (a2*c1 - a1*c2)/denom
34                 in if  x1 <= x && x <= x2 &&
35                        y1 <= y && y <= y2 &&
36                        x3 <= x && x <= x4 &&
37                        y3 <= y && y <= y4
38                    then Just (x,y)
39                    else Nothing
40            else Nothing
41
42
43 plantedToLines :: Planted -> [Line]
44 plantedToLines planted = runGeometryGenerator (plantPosition planted, 0) 0 $
45                 plantToGeometry (phenotype planted)
46
47 plantToGeometry :: Plant -> GeometryGenerator ()
48 plantToGeometry Bud = return ()
49 plantToGeometry (Stipe p) = addLine ((0,0),(0,stipeLength)) >>
50                             translated (0,stipeLength) (plantToGeometry p)
51 plantToGeometry (Fork p1 p2) = rotated (-pi/4) (plantToGeometry p1) >>
52                                rotated ( pi/4) (plantToGeometry p2)
53
54 -- | Lines are annotated with its plant, identified by the position
55 gardenToLines :: Garden -> [(Line, Double)]
56 gardenToLines = concatMap (\planted -> map (\line -> (line, plantPosition planted)) (plantedToLines planted))
57
58 -- | Add lightning from a given angle
59 lightenLines :: (Ord a, Show a) => Double -> [(Line, a)] -> [(Line, a, Double)]
60 lightenLines angle lines = let (lighted,_) = allKindsOfStuffWithAngle angle lines
61                            in lighted
62
63 lightPolygons :: (Ord a, Show a) => Double -> [(Line, a)] -> [(Point,Point,Point,Point,Double)]
64 lightPolygons angle lines = let (_,polygons) = allKindsOfStuffWithAngle angle lines
65                             in polygons
66
67 allKindsOfStuffWithAngle :: forall a. (Ord a, Show a) => Double -> [(Line, a)] ->
68                             (  [(Line, a, Double)]
69                             , [(Point,Point,Point,Point,Double)] )
70 allKindsOfStuffWithAngle angle lines = (lighted, polygons)
71   where projectLine :: Line -> (Double, Double)
72         projectLine (p1, p2) = (projectPoint p1, projectPoint p2)
73         projectPoint :: Point -> Double
74         projectPoint (x,y) = x + y * 1 / tan (pi-angle)
75         
76         -- False means Beginning of Line
77         sweepPoints :: [(Double, Bool, (Line, a))]
78         sweepPoints = sort $ concatMap (\l@((p1,p2),i) -> 
79                         if projectPoint p1 == projectPoint p2
80                         then []
81                         else if projectPoint p1 < projectPoint p2
82                              then [(projectPoint p1,False,l), (projectPoint p2,True,l)]
83                              else [(projectPoint p2,False,l), (projectPoint p1,True,l)]
84                 ) lines
85
86         -- Find all crossing points
87         crossings :: [Double]
88         crossings = case mapAccumL step [] sweepPoints of
89                         ([],crosses) -> nub (sort (concat crosses))
90                         _            -> error "Lines left open after sweep"
91           where -- accumulator is open lines, return is list of cross points
92                 step :: [Line] -> (Double, Bool, (Line, a)) -> ([Line], [Double])
93                 step [] (_, True, _)      = error $ "Line ends with no lines open"
94                 -- Beginning of a new line, mark it as open, and mark it as a cross-point
95                 step ol (x, False, (l,_)) = (l:ol, [x]) 
96                 -- End of a line. Calculate crosses with all open line, and remove it from the
97                 -- list of open lines
98                 step ol (x, True, (l,_)) = 
99                         let ol' = filter (/= l) ol
100                             crosses = map projectPoint $ mapMaybe (crossPoint l) ol'
101                         in (ol', x:crosses)
102
103         -- Cross points inverval
104         intervals = zip crossings (tail crossings)
105
106         unlighted = map (\(l,i) -> (l,i,0)) lines
107         
108         unprojectPoint x (p1@(x1,y1),p2@(x2,y2)) = 
109                 let t = (x - projectPoint p1) /
110                         (projectPoint p2 - projectPoint p1)
111                 in (x1 + t * (x2-x1), y1 + t * (y2-y1))
112
113         lineAtRay x l = let (x1',x2') = projectLine l
114                       in x1' <= x && x <= x2' || x2' <= x && x <= x1'
115
116         aboveFirst x l1 l2 =
117                 let (_,y1) = unprojectPoint x l1
118                     (_,y2) = unprojectPoint x l2
119                 in y2 `compare` y1
120
121         lighted :: [(Line, a, Double)]
122         lighted = foldl go unlighted intervals
123           where go llines (x1,x2) = curlines' ++ otherlines
124                   where -- Calculation based on the ray at the mid point
125                         mid = (x1 + x2) / 2
126                         -- Light intensity
127                         width = (x2 - x1) * sin angle
128                         (curlines, otherlines) = partition (\(l,_,_) -> lineAtRay mid l)
129                                                            llines
130                         sorted = sortBy (\(l1,_,_) (l2,_,_) -> aboveFirst mid l1 l2)
131                                         curlines
132                         curlines' = snd $ mapAccumL shine 1 sorted
133                         shine intensity (l,i,amount) = ( intensity * lightFalloff
134                                                        , (l,i,amount + intensity * width))
135
136         polygons = concatMap go intervals
137           where go (x1,x2) = if null sorted then [nothingPoly] else lightedPolys
138                   where mid = (x1 + x2) / 2
139                         curlines = filter (lineAtRay mid) (map fst lines)
140                         sorted = sortBy (aboveFirst mid) curlines
141                         ceiling = ((0,10),(1,10))
142                         floor = ((0,0),(1,0))
143                         nothingPoly = let p1 = unprojectPoint x1 ceiling
144                                           p2 = unprojectPoint x1 floor
145                                           p3 = unprojectPoint x2 floor
146                                           p4 = unprojectPoint x2 ceiling
147                                       in (p1,p2,p3,p4,1)
148                         firstPoly = let p1 = unprojectPoint x1 ceiling
149                                         p2 = unprojectPoint x1 (head sorted)
150                                         p3 = unprojectPoint x2 (head sorted)
151                                         p4 = unprojectPoint x2 ceiling
152                                     in (p1,p2,p3,p4)
153                         lastPoly =  let p1 = unprojectPoint x1 (last sorted)
154                                         p2 = unprojectPoint x1 floor
155                                         p3 = unprojectPoint x2 floor
156                                         p4 = unprojectPoint x2 (last sorted)
157                                     in (p1,p2,p3,p4)
158                         polys = zipWith (\l1 l2 ->
159                                          let p1 = unprojectPoint x1 l1
160                                              p2 = unprojectPoint x1 l2
161                                              p3 = unprojectPoint x2 l2
162                                              p4 = unprojectPoint x2 l1
163                                          in (p1,p2,p3,p4)) sorted (tail sorted)
164                         polys' = [firstPoly] ++ polys ++ [lastPoly]
165                         lightedPolys = snd $ mapAccumL shine (1*lightFalloff) polys'
166                         shine intensity (p1,p2,p3,p4) = ( intensity * lightFalloff
167                                                         , (p1,p2,p3,p4,intensity))
168
169
170                       
171 {-
172 -- Yay, this is a sweep-line-algorithm from Kognitive Systeme, whe would have guessed
173 mergeLines lines = catMaybes $ snd $ mapAccumL step (Nothing, 0) points
174  where points = sort $ concatMap (
175                 \(p1,p2) -> if p1 < p2 then [(p1,False),(p2,True)] else [(p2,False),(p1,True)]
176                 ) lines
177        step (Nothing, 0) (p, False) = ((Just p, 1),    Nothing)
178        step (Nothing, 0) (p, True)  = error $ "End before start, point " ++ show p
179        step (Just p1, 1) (p2,True)  = ((Nothing, 0),   Just (p1,p2))
180        step (Just p1, n) (_, False) = ((Just p1, n+1), Nothing)
181        step (Just p1, n) (_, True)  = ((Just p1, n-1), Nothing)
182
183
184 -}