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