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