packages feed

mandulia-0.4: src/AmmannA3.hs

{-
Mandulia -- Mandelbrot/Julia explorer
Copyright (C) 2010  Claude Heiland-Allen <claudiusmaximus@goto10.org>

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
-}

module AmmannA3 (AmmannA3(), ammannA3, Tile(..), Tile'(..), tiles, zoom, zoomTo, stepIn) where

import Data.Maybe (fromMaybe)
import Data.Tree

import Bounds
import Utils
import Vector

data Tile = A | B | C
  deriving (Show, Read, Eq, Ord, Enum, Bounded)

bounds0 :: Tile -> Bounds
bounds0 x = bounds0' !! fromEnum x

bounds0' :: [Bounds]
bounds0' = [ sbound 1 phi' -- A
           , sbound 2 1    -- B
           , sbound 1 1    -- C
           ]

sbound :: R -> R -> Bounds
sbound w h =
  bounds [ V x y 1 | x <- [negate w, w], y <- [negate h, h] ]

transforms :: [( Tile, [( Tile, (M, Integer) )] )]
transforms =
  [ (A, [ (B, (t ( 0) (p * n/2 - f/2)                 0                             , 1))
        , (A, (t (-1) (f/2 - p * e/2)                 0                             , 2)) ])
  , (B, [ (A, (t ( 1) (p * e/2               - n/2)   (p * f/2 - k/2 - m/2)         , 3))
        , (C, (t ( 2) (p * (a + x/2 + v/2)   - n/2)   (p * (s/2 + u/2) - k/2 - m/2) , 4))
        , (A, (t (-1) (p * (a + v - c + e/2) - n/2)   (k/2 + m/2         - p * f/2) , 5))
        , (A, (t ( 0) (n/2               - p * f/2)   (p * e/2         - k/2 - m/2) , 6)) ])
  , (C, [ (C, (t ( 2) (p * (x/2 + v/2) - x/2 - v/2)   (p * (s/2 + u/2) - s/2 - u/2) , 7))
        , (A, (t (-1) (p * (v - c + e/2) - x/2 - v/2) (s/2 + u/2 - p * f/2)         , 8))
        , (A, (t ( 0) (x/2 + v/2 - p * f/2)           (p * e/2        - s/2 - u/2)  , 9)) ])
  ]
  where
    t da dx dy = translate (dx*2) (dy*2) ^^*^^ rotate (da * pi / 2) ^^*^^ scale p p
    p = phi'
    a = p * p
    c = p * p * p
    e = p
    f = 1
    k = p
    m = p * p
    n = (1 - p * p * p) / p
    s = p
    u = p * p
    v = p * p * p + p
    x = p * p * p * p

centerC :: V
centerC =
  let Just ts    = lookup C transforms
      Just (t,_) = lookup C ts
      ps = iterate (t ^^*^) (V 0 0 1)
  in  ps !! 256

inRadiusC :: R
inRadiusC =
  let cornerC = V (1/2 - phi') (1/2 - phi' * phi') 1
  in  cornerC ^|-|^ centerC

data Tile'' =
  Tile''
    { ttTile      :: !Tile
    , ttId        :: !Integer
    , ttTransform :: !M
    }

builder :: Tile'' -> (Tile'', [Tile''])
builder tm = tm `seq` (tm, map mkTile (mine transforms))
  where
    mine = concatMap snd . filter ((==) (ttTile tm) . fst)
    mkTile (x, (mm, j)) =
      Tile''
        { ttTile = x
        , ttId = 10 * ttId tm + j
        , ttTransform = ttTransform tm ^^*^^ mm
        }

data Tile' =
  Tile'
    { tTile   :: !Tile
    , tBounds :: !Bounds
    , tCenter :: !V
    , tDepth  :: !Int
    , tLevel  :: !Int
    , tId     :: !Integer
    }

tree :: R -> Tree Tile'
tree maxRadius =
  let s = maxRadius / inRadiusC
      V x y _ = centerC
      tr = scale s s ^^*^^ translate (-x) (-y)
      t0 = Tile''{ ttTile = C, ttId = 7, ttTransform = tr }
  in toTiles (Just (V 0 0 1)) C 0 (tree' t0)

tree' :: Tile'' -> Tree Tile''
tree' t = unfoldTree builder t

tB :: M -> Tile -> Bounds
tB m t = m `transform'` bounds0 t

toTiles :: Maybe V -> Tile -> Int -> Tree Tile'' -> Tree Tile'
toTiles v0 t0 level tr =
  let Tile''{ ttTile = t, ttId = n, ttTransform = m } = rootLabel tr
      ts  = subForest tr
      v1 = fromMaybe (centerPoint m) v0
      v2 = if t0 == C && t == C then v0 else Nothing
      b1 = tB m t
      nn = normalizeId n
      tile =
        Tile'
          { tTile = t
          , tBounds = b1
          , tCenter = v1
          , tDepth  = idToLevel' nn
          , tLevel  = level
          , tId = nn
          }
      level' = level + 1
      forest = level' `seq` map (toTiles v2 t level') ts
  in tile `seq` Node{ rootLabel = tile, subForest = forest }

centerPoint :: M -> V
centerPoint = (^^*^ centerC)

data LevelA3 =
  LevelA3
    { lInnerTiles :: Forest Tile'
    , lOuterTiles :: Forest Tile'
    , lBounds     :: Bounds
    }

data AmmannA3 =
  AmmannA3
    { aOuter  :: [LevelA3]
    , aFocus  ::  LevelA3
    , aBounds :: Bounds
    , aRadius :: R
    }

ammannA3 :: Bounds -> AmmannA3
ammannA3 box =
  let r = diagonal box / 2
      (is, os, _) = triPart box [tree r]
      l = LevelA3{ lInnerTiles = is, lOuterTiles = os, lBounds = box }
  in  AmmannA3{ aOuter = [l], aFocus = l, aBounds = box, aRadius = r }

triPart :: Bounds -> [Tree Tile'] -> ([Tree Tile'], [Tree Tile'], [Tree Tile'])
triPart box = foldr go ([],[],[])
  where
    go t (is, es, os)
      | b `insideOrEqual` box = (t:is, es, os)
      | b `outside`       box = (is, es, t:os)
      | otherwise             = (is, t:es, os)
      where b = tBounds . rootLabel $ t

zoomTo :: Bounds -> AmmannA3 -> Maybe AmmannA3
zoomTo box a3
  | box `insideOrEqual` region =
    (if factor >= phi'
      then Just
      else zoomTo box . (\a -> a{ aRadius = phi' * aRadius a }) . stepIn) $
        let focus = aFocus a3
            ots   = prune box (lOuterTiles focus)
            (its, ots', _) = triPart box (lInnerTiles focus)
        in  a3{ aFocus = focus{ lOuterTiles = ots' ++ ots
                              , lInnerTiles = its
                              , lBounds = box
                              } }
  | otherwise = zoomTo box =<< stepOut a3{ aRadius = phi  * aRadius a3 }
  where
    factor = radius / aRadius a3
    radius = diagonal box / 2
    region = lBounds . aFocus $ a3

zoom :: R -> AmmannA3 -> Maybe AmmannA3
zoom factor a3 = flip zoomTo a3 . expand factor . lBounds . aFocus $ a3

prune :: Bounds -> Forest Tile' -> Forest Tile'
prune box = filter (not . outside box . tBounds . rootLabel)

stepOut :: AmmannA3 -> Maybe AmmannA3
stepOut a3 =
  case aOuter a3 of
    []     -> Nothing
    os@[l] -> Just a3{ aOuter = os, aFocus = l, aRadius = aRadius a3 * phi' }
    (l:os) -> Just a3{ aOuter = os, aFocus = l }

stepIn  :: AmmannA3 -> AmmannA3
stepIn  a3 =
  let l0  =    aFocus a3
      os  = l0:aOuter a3
      box = lBounds l0
      its =                          concatMap subForest . lInnerTiles $ l0
      (its', ots, _) = triPart box . concatMap subForest . lOuterTiles $ l0
      l   = l0{ lInnerTiles = its' ++ its, lOuterTiles = ots }
  in a3{ aOuter = os, aFocus = l }

tiles :: Int -> AmmannA3 -> [Tile']
tiles lod = map rootLabel . (\l -> lOuterTiles l ++ lInnerTiles l) . aFocus . (!!lod) . iterate stepIn

normalizeId :: Integer -> Integer -- C=>C is transform 7
normalizeId n = let (d, m) = n `divMod` 10 in if m == 7 then normalizeId d else n

idToLevel' :: Integer -> Int -- n must be normalized
idToLevel' n = snd . head . dropWhile ((<n) . fst) $ tens

tens :: [(Integer, Int)]
tens = iterate (10*) 1 `zip` [0..]