Feature mapviz (#61)

* Fix framerate issue for gifs.

* Fix aeson version on lts-11

* Improve transformation performance.

* Vastly improve performance for SVG serialization.

* Improve geo projection performance.

* Add reference to cartogram program.

* Add rasterSized function.

* Better clippaths.

* Use newer reanimate-svg package.

* Show time spent and time remaining when rendering.

* Support external images and multiple rasters.

* Improve performance of bounding box calculations.

* Evict macos cache.

Former-commit-id: 3eb65d85fc2bb7dd9d67486d2da34d024767717d
This commit is contained in:
David Himmelstrup 2020-03-15 19:35:07 +08:00 committed by GitHub
commit 6608fb7410
19 changed files with 780 additions and 301 deletions

View file

@ -20,7 +20,7 @@ jobs:
steps:
- task: Cache@2
inputs:
key: ${{ parameters.name }} | "${{ parameters.vmImage }}" | $(STACK_YAML) | stack-root | $(Agent.OS) | version1
key: ${{ parameters.name }} | "${{ parameters.vmImage }}" | $(STACK_YAML) | stack-root | $(Agent.OS) | version2
path: $(STACK_ROOT)
cacheHitVar: CACHE_RESTORED
displayName: Cache stack root

View file

@ -6,7 +6,8 @@
module Main (main) where
import Reanimate
import Reanimate.Povray (povraySlow)
import Reanimate.Povray (povraySlow')
import Reanimate.Raster
import Codec.Picture
import Codec.Picture.Types
@ -32,15 +33,14 @@ main = reanimate $ parA bg $ sceneAnimation $ do
t <- spriteT
dur <- spriteDuration
pure $
povraySlow [] $
mkImage screenWidth screenHeight $ povraySlow' [] $
script (svgAsPngFile (texture (t/dur))) transZ getX getZ
wait 2
tweenVar zPos 9 (\t v -> fromToS v 8 (t/9))
tweenVar xRot 9 (\t v -> fromToS v 360 $ curveS 2 (t/9))
tweenVar zRot 9 (\t v -> fromToS v 360 $ curveS 2 (t/9))
fork $ tweenVar zPos 9 $ \v -> fromToS v 8
fork $ tweenVar xRot 9 $ \v -> fromToS v 360 . curveS 2
fork $ tweenVar zRot 9 $ \v -> fromToS v 360 . curveS 2
wait 10
tweenVar zPos 2 (\t v -> fromToS v 0 $ curveS 3 (t/2))
wait 2
tweenVar zPos 2 $ \v -> fromToS v 0 . curveS 3
where
bg = animate $ const $ mkBackgroundPixel $ PixelRGBA8 252 252 252 0xFF

View file

@ -145,6 +145,8 @@ reanimate animation = do
,"--target", target
,"+RTS", "-N", "-RTS"]
else do
raster <- selectRaster renderRaster
setRaster raster
setFPS fps
setWidth width
setHeight height
@ -153,9 +155,10 @@ reanimate animation = do
\ width: %d\n\
\ height: %d\n\
\ fmt: %s\n\
\ target: %s\n"
fps width height (showFormat fmt) target
raster <- selectRaster renderRaster
\ target: %s\n\
\ raster: %s\n"
fps width height (showFormat fmt) target (show raster)
render animation target raster fmt width height fps
selectRaster :: Raster -> IO Raster

View file

@ -1,4 +1,5 @@
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE MultiWayIf #-}
module Reanimate.GeoProjection
( Projection(..)
@ -6,6 +7,7 @@ module Reanimate.GeoProjection
, LonLat(..)
, project
, interpP
, interpBBP
, mergeP
, isValidP
, scaleP
@ -17,6 +19,7 @@ module Reanimate.GeoProjection
, mercatorP
, mollweideP
, hammerP
, cylindricalEqualAreaP
, lambertP
, bottomleyP
, sinusoidalP
@ -42,22 +45,33 @@ module Reanimate.GeoProjection
import Codec.Picture
import Codec.Picture.Types
import Control.Lens ((^.))
import Control.Lens ((^.))
import Control.Monad
import Control.Monad.ST
import Control.Monad.ST.Unsafe
import Data.Aeson
import Data.Foldable
import Data.Geospatial hiding (LonLat)
import Data.Geospatial hiding (LonLat)
import Data.Hashable
import Data.LinearRing
import qualified Data.LineString as Line
import qualified Data.LineString as Line
import Data.Vector.Storable (unsafeWith)
import qualified Data.Vector.Unboxed as V
import Debug.Trace
import Graphics.SvgTree (Tree (None))
import Linear (distance, lerp)
import Linear.V2 hiding (angle)
import Foreign
import GHC.Exts (Double (..), cosDouble#, sinDouble#,
(*##), (+##), (-##), (/##))
import Graphics.SvgTree (Tree (None))
import Linear (distance, lerp)
import Linear.V2 hiding (angle)
import Reanimate
import System.IO.Unsafe
{-# INLINE halfPi #-}
{-# INLINE sqrtPi #-}
{-# INLINE sqrt2 #-}
{-# INLINE epsilon #-}
{-# INLINE tau #-}
-- Constants
halfPi, sqrtPi, sqrt2, epsilon, tau :: Double
halfPi = pi/2
@ -70,136 +84,225 @@ toRads, cot :: Double -> Double
toRads dec = dec/180 * pi
cot = recip . tan
srcPixel :: Pixel pixel => Image pixel -> LonLat -> pixel
srcPixel src (LonLat lam phi) =
pixelAt src xPx yPx
-- pixelAt src xPx yPx
unsafePixelAt (imageData src) (pixelBaseIndex src xPx yPx)
where
!xPx = round $ ((lam+pi)/tau) * fromIntegral (imageWidth src-1)
!yPx = round $ (1-((phi+halfPi)/pi)) * fromIntegral (imageHeight src-1)
srcPixelFast :: Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
srcPixelFast src wMax hMax (LonLat lam phi) = unsafeIOToST $
unsafeWith (imageData src) $ \ptr -> do
let ptr' = plusPtr ptr idx
r <- peek ptr'
g <- peek $ plusPtr ptr' 1
b <- peek $ plusPtr ptr' 2
a <- peek $ plusPtr ptr' 3
return $ PixelRGBA8 r g b a
where
!idx = pixelBaseIndex src xPx yPx
!xPx = round $ ((lam+pi)/tau) * wMax
!yPx = round $ (1-((phi+halfPi)/pi)) * hMax
{- HLINT ignore -}
findValidCoord :: Image pixel -> Projection -> XYCoord -> XYCoord
findValidCoord !src !p (XYCoord x y) = foldr const (XYCoord x y)
[ XYCoord x' y'
| let !xi = round $ x * wMax
!yi = round $ y * hMax
, !ax <- [xi, xi-1, xi+1]
, ax >= 0
, ax < w
, let !x' = fromIntegral ax / wMax
, !ay <- [yi, yi-1, yi+1]
, ay >= 0
, ay < h
, let !y' = fromIntegral ay / hMax
, validLonLat $! projectionInverse p $! XYCoord x' y'
]
where
!w = imageWidth src
!h = imageHeight src
!wMax = fromIntegral (w-1)
!hMax = fromIntegral (h-1)
-- findValidCoord :: Int -> Int -> Double -> Double -> (XYCoord -> LonLat) -> XYCoord -> XYCoord
-- findValidCoord !w !h !wMax !hMax !p_inv (XYCoord x y) = foldr const (XYCoord x y)
-- [ XYCoord x' y'
-- | let !xi = round $ x * wMax
-- !yi = round $ y * hMax
-- , !ax <- [xi, xi-1, xi+1]
-- , ax >= 0
-- , ax < w
-- , let !x' = fromIntegral ax / wMax
-- , !ay <- [yi, yi-1, yi+1]
-- , ay >= 0
-- , ay < h
-- , let !y' = fromIntegral ay / hMax
-- , validLonLat $! p_inv $! XYCoord x' y'
-- ]
isInWorld :: Projection -> XYCoord -> Bool
isInWorld p coord =
odd $ length $ isInWorld' p coord
-- isInWorld :: Projection -> XYCoord -> Bool
-- isInWorld p coord =
-- odd $ length $ isInWorld' p coord
--
-- isInWorld' :: Projection -> XYCoord -> [(Double, Double)]
-- isInWorld' p (XYCoord x y) =
-- [ (x1, y1)
-- | (XYCoord x1 y1, XYCoord x2 y2) <- world
-- , (y1 > y) /= (y2 > y) -- y is between y1 and y2
-- , x < (x2 - x1) * (y - y1) / (y2 - y1) + x1 -- x is to the left of the line
-- ]
-- where
-- world = worldPolygon p
--
-- worldPolygon :: Projection -> [(XYCoord, XYCoord)]
-- worldPolygon p =
-- interp (-pi, -halfPi) (-pi, halfPi) ++
-- interp (-pi, halfPi) (pi, halfPi) ++
-- interp (pi, halfPi) (pi, -halfPi) ++
-- interp (pi, -halfPi) (-pi, -halfPi)
-- where
-- apply (lam, phi) = projectionForward p $ LonLat lam phi
-- steps = 100
-- interp (x1, y1) (x2,y2) =
-- [ ( apply (fromToS x1 x2 (n/steps), fromToS y1 y2 (n/steps))
-- , apply (fromToS x1 x2 ((n+1)/steps), fromToS y1 y2 ((n+1)/steps)))
-- | n <- [0..steps-1]]
isInWorld' :: Projection -> XYCoord -> [(Double, Double)]
isInWorld' p (XYCoord x y) =
[ (x1, y1)
| (XYCoord x1 y1, XYCoord x2 y2) <- world
, (y1 > y) /= (y2 > y) -- y is between y1 and y2
, x < (x2 - x1) * (y - y1) / (y2 - y1) + x1 -- x is to the left of the line
]
where
world = worldPolygon p
-- findNearestPixel :: MutableImage s PixelRGBA8 -> Int -> Int -> Int -> Int -> ST s PixelRGBA8
-- findNearestPixel src w h srcX srcY = worker $ take 20
-- [ (x, y)
-- | n <- [1..]
-- , x <- [srcX-n .. srcY+n]
-- , y <- if x == srcX-n || x == srcX+n then [srcY-n,srcY+n] else [srcY-n .. srcY+n]
-- , x >= 0
-- , y >= 0
-- , x < w
-- , y < h
-- ]
-- where
-- worker [] = pure $ PixelRGBA8 0xFF 0x00 0x00 0xFF
-- worker ((x,y):rest) = do
-- this <- readPixel src x y
-- if this == blank
-- then worker rest
-- else return this
-- blank = PixelRGBA8 0x00 0x00 0x00 0x00
worldPolygon :: Projection -> [(XYCoord, XYCoord)]
worldPolygon p =
interp (-pi, -halfPi) (-pi, halfPi) ++
interp (-pi, halfPi) (pi, halfPi) ++
interp (pi, halfPi) (pi, -halfPi) ++
interp (pi, -halfPi) (-pi, -halfPi)
where
apply (lam, phi) = projectionForward p $ LonLat lam phi
steps = 100
interp (x1, y1) (x2,y2) =
[ ( apply (fromToS x1 x2 (n/steps), fromToS y1 y2 (n/steps))
, apply (fromToS x1 x2 ((n+1)/steps), fromToS y1 y2 ((n+1)/steps)))
| n <- [0..steps-1]]
findNearestPixel :: MutableImage s PixelRGBA8 -> Int -> Int -> Int -> Int -> ST s PixelRGBA8
findNearestPixel src w h srcX srcY = worker $ take 20
[ (x, y)
| n <- [1..]
, x <- [srcX-n .. srcY+n]
, y <- if x == srcX-n || x == srcX+n then [srcY-n,srcY+n] else [srcY-n .. srcY+n]
, x >= 0
, y >= 0
, x < w
, y < h
]
where
worker [] = pure $ PixelRGBA8 0xFF 0x00 0x00 0xFF
worker ((x,y):rest) = do
this <- readPixel src x y
if this == blank
then worker rest
else return this
blank = PixelRGBA8 0x00 0x00 0x00 0x00
-- Original version: 134,925 pixels/second
-- Inlined projections: 136,332 pixels/second
-- TEST: no write pixels: 134,288 pixels/second !!!
-- Fast theta: 1,489,719 pixels/second
-- Cached theta: 3,622,254 pixels/second
-- to equirectangularP: 9,015,326 pixels/second
-- 9,680,217
-- 9,466,744
-- 14,830,330
-- to lambertP: 6,735,973
-- to mercatorP: 4,248,927
-- to mollweideP: 3,593,545
-- to hammerP: 3,237,699
-- to bottomleyP: 3,864,848
-- to sinusoidalP: 7,020,756
-- to wernerP: 4,433,295
-- to bonneP: 4,071,187
-- to augustP: 2,553,177
-- to collignonP: 5,486,849
-- to eckert1P: 7,428,849
-- to eckert3P: 6,666,936
-- to eckert5P: 6,106,492
-- to faheyP: 5,137,291
-- to foucautP: 3,983,151
-- to lagrangeP: 3,850,611
-- interpFastP :: Image PixelRGBA8 -> Projection -> Projection -> Double -> Image PixelRGBA8
-- interpFastP !src (Projection _ p1 p1_inv) (Projection _ p2 p2_inv) !t = runST $ do
-- unsafeIOToST $ putStrLn "Allocating new array"
-- !img <- newMutableImage w h
-- unsafeIOToST $ putStrLn "done"
-- start <- unsafeIOToST $ getCurrentTime
-- let factor = 2
-- total = w*factor * h*factor
-- let l1 =
-- loopTo (w*factor) $ \x -> do
-- loopTo (h*factor) $ \y -> do
-- let thisIndex = (x*h*factor+y)
-- when (thisIndex `mod` 1000000 == 0) $ unsafeIOToST $ do
-- now <- getCurrentTime
-- let diff = realToFrac (diffUTCTime now start):: Double
-- printf "%.2f pixels/second\n" (fromIntegral (total-thisIndex) / diff)
-- let !x1' = fromIntegral x / (wMax*fromIntegral factor)
-- !y1' = fromIntegral y / (hMax*fromIntegral factor)
-- !lonlat = p1_inv $! XYCoord x1' y1'
-- -- p = srcPixel src lonlat
-- -- unsafeIOToST (evaluate lonlat)
--
-- when (validLonLat lonlat) $ do
-- p <- srcPixelFast src wMax hMax lonlat
-- when (pixelOpacity p /= 0) $ do
-- let XYCoord !x1 !y1 = p1 lonlat
-- -- !coord = p2 lonlat
-- XYCoord !x2 !y2 = p2 lonlat -- findValidCoord w h wMax hMax p2_inv $ p2 lonlat
-- !x3 = round $ fromToS x1 x2 t * wMax
-- !y3 = round $ (1 - fromToS y1 y2 t) * hMax :: Int
-- -- unsafeIOToST (evaluate coord)
-- -- return ()
-- when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
-- writePixel img x3 y3 p
-- l1
-- unsafeFreezeImage img
-- where
-- loopTo m fn = go m
-- where go 0 = return ()
-- go n = fn (n-1) >> go (n-1)
-- !w = imageWidth src
-- !h = imageHeight src
-- !wMax = fromIntegral (w-1)
-- !hMax = fromIntegral (h-1)
interpP :: Image PixelRGBA8 -> Projection -> Projection -> Double -> Image PixelRGBA8
interpP !src !p1 !p2 !t = runST $ do
interpP src p1 _ 0 = project src p1
interpP src _ p2 1 = project src p2
interpP !src (Projection _label1 p1 p1_inv) !(Projection _label2 p2 p2_inv) !t = runST $ do
!img <- newMutableImage w h
let blank = PixelRGBA8 0x00 0x00 0x00 0x00
let isBlank pixel = pixel == blank
-- forM_ [0..w-1] $ \x ->
-- forM_ [0..h-1] $ \y -> do
-- let x1 = fromIntegral x / (wMax)
-- y1 = 1 - fromIntegral y / (hMax)
-- when (isInWorld (mergeP p1 p2 t) (XYCoord x1 y1)) $
-- writePixel img x y $ PixelRGBA8 0xFF 0x00 0x00 0xFF
let factor = 2
let l1 =
-- total = w*factor * h*factor
let l1 = do
-- start <- unsafeIOToST $ getCurrentTime
loopTo (w*factor) $ \x -> do
loopTo (h*factor) $ \y -> do
-- let thisIndex = (x*h*factor+y)
-- when (thisIndex `mod` 1000000 == 0) $ unsafeIOToST $ do
-- now <- getCurrentTime
-- let diff = realToFrac (diffUTCTime now start):: Double
-- printf "%.2f pixels/second: %s\n" (fromIntegral (total-thisIndex) / diff) label1
let !x1' = fromIntegral x / (wMax*fromIntegral factor)
!y1' = fromIntegral y / (hMax*fromIntegral factor)
!lonlat = projectionInverse p1 $! XYCoord x1' y1'
p = srcPixel src lonlat
when (validLonLat lonlat && pixelOpacity p /= 0) $ do
let XYCoord x1 y1 = projectionForward p1 lonlat
XYCoord x2 y2 = findValidCoord src p2 $ projectionForward p2 lonlat
!x3 = round $ fromToS x1 x2 t * wMax
!y3 = round $ (1 - fromToS y1 y2 t) * hMax
when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
writePixel img x3 y3 p
l2 =
!lonlat = p1_inv $! XYCoord x1' y1'
-- p = srcPixel src lonlat
when (validLonLat lonlat) $ do
p <- srcPixelFast src wMax hMax lonlat
when (pixelOpacity p /= 0) $ do
let XYCoord x1 y1 = p1 lonlat
-- XYCoord x2 y2 = findValidCoord w h wMax hMax p2_inv $ p2 lonlat
XYCoord x2 y2 = p2 lonlat
!x3 = round $ fromToS x1 x2 t * wMax
!y3 = round $ (1 - fromToS y1 y2 t) * hMax
when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
writePixel img x3 y3 p
l2 = do
-- start <- unsafeIOToST $ getCurrentTime
loopTo (w*factor) $ \x ->
loopTo (h*factor) $ \y -> do
-- let thisIndex = (x*h*factor+y)
-- when (thisIndex `mod` 1000000 == 0) $ unsafeIOToST $ do
-- now <- getCurrentTime
-- let diff = realToFrac (diffUTCTime now start):: Double
-- printf "%.2f pixels/second: %s\n" (fromIntegral (total-thisIndex) / diff) label2
let !x2' = fromIntegral x / (wMax*fromIntegral factor)
!y2' = fromIntegral y / (hMax*fromIntegral factor)
!lonlat = projectionInverse p2 (XYCoord x2' y2')
p = srcPixel src lonlat
when (validLonLat lonlat && pixelOpacity p /= 0) $ do
let XYCoord x2 y2 = projectionForward p2 lonlat
XYCoord x1 y1 = findValidCoord src p1 $ projectionForward p1 lonlat
!x3 = round $ fromToS x1 x2 t * wMax
!y3 = round $ (1 - fromToS y1 y2 t) * hMax
when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $ do
writePixel img x3 y3 p
!lonlat = p2_inv (XYCoord x2' y2')
-- p = srcPixel src lonlat
when (validLonLat lonlat) $ do
p <- srcPixelFast src wMax hMax lonlat
when (pixelOpacity p /= 0) $ do
let XYCoord x2 y2 = p2 lonlat
-- XYCoord x1 y1 = findValidCoord w h wMax hMax p1_inv $ p1 lonlat
XYCoord x1 y1 = p1 lonlat
!x3 = round $ fromToS x1 x2 t * wMax
!y3 = round $ (1 - fromToS y1 y2 t) * hMax
when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $ do
writePixel img x3 y3 p
if t < 0.5
then l1 >> l2
else l2 >> l1
when False $
forM_ [0..w-1] $ \x ->
forM_ [0..h-1] $ \y -> do
let x1 = fromIntegral x / (wMax)
y1 = 1 - fromIntegral y / (hMax)
this <- readPixel img x y
when (isBlank this) $
when (isInWorld (mergeP p1 p2 t) (XYCoord x1 y1)) $
writePixel img x y =<< findNearestPixel img w h x y
unsafeFreezeImage img
where
loopTo m fn = go m
@ -210,6 +313,76 @@ interpP !src !p1 !p2 !t = runST $ do
!wMax = fromIntegral (w-1)
!hMax = fromIntegral (h-1)
interpBBP :: Image PixelRGBA8 -> Projection -> Projection ->
(Double,Double,Double,Double) -> (Double,Double,Double,Double) -> Double -> Image PixelRGBA8
interpBBP !src (Projection _ p1 p1_inv) !(Projection _ p2 p2_inv) (fx,fy,fw,fh) (tx, ty, tw, th) !t = runST $ do
!img <- newMutableImage w h
let factor = 2
let l1 =
loopTo (w*factor) $ \x -> do
loopTo (h*factor) $ \y -> do
let !x1' = fromIntegral x / (wMax*fromIntegral factor)
!y1' = fromIntegral y / (hMax*fromIntegral factor)
when (x1' >= fx && x1' <= fx+fw && y1' >= fy && y1' <= fy+fh) $ do
let !lonlat = p1_inv $! XYCoord x1' y1'
-- p = srcPixel src lonlat
when (validLonLat lonlat) $ do
-- let LonLat lam phi = lonlat
-- !xPx = ((lam+pi)/tau)
-- !yPx = (((phi+halfPi)/pi))
-- when (xPx >= fx && xPx <= fx+fw && yPx >= fy && yPx <= fy+fh) $ do
p <- srcPixelFast src wMax hMax lonlat
when (pixelOpacity p /= 0) $ do
let XYCoord x1 y1 = p1 lonlat
XYCoord x2 y2 = p2 lonlat
-- XYCoord x2 y2 = findValidCoord w h wMax hMax p2_inv $ p2 lonlat
!x3 = round $ fromToS x1 x2 t * wMax
!y3 = round $ (1 - fromToS y1 y2 t) * hMax
when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
writePixel img x3 y3 p
l2 =
loopTo (w*factor) $ \x ->
loopTo (h*factor) $ \y -> do
let !x2' = fromIntegral x / (wMax*fromIntegral factor)
!y2' = fromIntegral y / (hMax*fromIntegral factor)
when (x2' >= tx && x2' <= tx+tw && y2' >= ty && y2' <= ty+th) $ do
let !lonlat = p2_inv (XYCoord x2' y2')
-- p = srcPixel src lonlat
when (validLonLat lonlat) $ do
p <- srcPixelFast src wMax hMax lonlat
when (pixelOpacity p /= 0) $ do
let XYCoord x2 y2 = p2 lonlat
-- XYCoord x1 y1 = findValidCoord w h wMax hMax p1_inv $ p1 lonlat
XYCoord x1 y1 = p1 lonlat
!x3 = round $ fromToS x1 x2 t * wMax
!y3 = round $ (1 - fromToS y1 y2 t) * hMax
when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $ do
writePixel img x3 y3 p
if t < 0.5
then l1 >> l2
else l2 >> l1
-- when False $
-- forM_ [0..w-1] $ \x ->
-- forM_ [0..h-1] $ \y -> do
-- let x1 = fromIntegral x / (wMax)
-- y1 = 1 - fromIntegral y / (hMax)
-- this <- readPixel img x y
-- when (isBlank this) $
-- when (isInWorld (mergeP p1 p2 t) (XYCoord x1 y1)) $
-- writePixel img x y =<< findNearestPixel img w h x y
unsafeFreezeImage img
where
loopTo m fn = go m
where go 0 = return ()
go n = fn (n-1) >> go (n-1)
!w = imageWidth src
!h = imageHeight src
!wMax = fromIntegral (w-1)
!hMax = fromIntegral (h-1)
eqLonLat :: LonLat -> LonLat -> Bool
eqLonLat (LonLat x1 y1) (LonLat x2 y2)
= eqDouble x1 x2 && eqDouble y1 y2
@ -225,14 +398,19 @@ data XYCoord = XYCoord !Double !Double -- 0 to 1
deriving (Read,Show,Eq,Ord)
data LonLat = LonLat !Double !Double -- -pi to +pi, -halfPi to +halfPi
deriving (Read,Show,Eq,Ord)
instance Hashable LonLat where
hashWithSalt s (LonLat a b) = hashWithSalt s (a,b)
data Projection = Projection
{ projectionForward :: !(LonLat -> XYCoord)
{ projectionLabel :: String
, projectionForward :: !(LonLat -> XYCoord)
--
-- (Double# -> Double# -> (# Double#, Double# #))
, projectionInverse :: !(XYCoord -> LonLat)
}
-- FIXME: Verify that 'src' has an aspect ratio of 2:1.
project :: Image PixelRGBA8 -> Projection -> Image PixelRGBA8
project src (Projection _ pInv) = generateImage fn w h
project src (Projection _label _ pInv) = generateImage fn w h
where
w = imageWidth src
h = imageHeight src
@ -253,7 +431,7 @@ _validXYCoord :: XYCoord -> Bool
_validXYCoord (XYCoord x y) = x >= 0 && x <= 1 && y >= 0 && y <= 1
isValidP :: Projection -> Bool
isValidP (Projection p pInv) = and
isValidP (Projection _label p pInv) = and
[ check x y
| x <- [0..w-1::Int]
, y <- [0..h-1::Int] ]
@ -269,7 +447,7 @@ isValidP (Projection p pInv) = and
|| trace (show (lonlat, lonlat2)) False
moveBottomP :: Double -> Projection -> Projection
moveBottomP offset (Projection p pInv) = Projection p' pInv'
moveBottomP offset (Projection label p pInv) = Projection label p' pInv'
where
p' (LonLat lon lat) =
case p (LonLat lon lat) of
@ -280,7 +458,7 @@ moveTopP :: Double -> Projection -> Projection
moveTopP offset = flipYAxisP . moveBottomP offset . flipYAxisP
flipYAxisP :: Projection -> Projection
flipYAxisP (Projection p pInv) = Projection p' pInv'
flipYAxisP (Projection label p pInv) = Projection label p' pInv'
where
p' (LonLat lam phi) =
let XYCoord x y = p (LonLat lam (negate phi))
@ -290,18 +468,17 @@ flipYAxisP (Projection p pInv) = Projection p' pInv'
in LonLat lam (negate phi)
scaleP :: Double -> Double -> Projection -> Projection
scaleP xScale yScale (Projection p pInv) = Projection forward inverse
scaleP xScale yScale (Projection label p pInv) = Projection label forward inverse
where
forward lonlat =
case p lonlat of
XYCoord x y -> XYCoord ((x-0.5)*xScale+0.5) ((y-0.5)*yScale+0.5)
inverse (XYCoord x y) =
let new = XYCoord ((x-0.5)/xScale+0.5) ((y-0.5)/yScale+0.5)
in pInv new
pInv $ XYCoord ((x-0.5)/xScale+0.5) ((y-0.5)/yScale+0.5)
mergeP :: Projection -> Projection -> Double -> Projection
mergeP p1 p2 t = Projection p pInv
mergeP p1 p2 t = Projection (projectionLabel p1 ++ "/" ++ projectionLabel p2) p pInv
where
p lonlat =
let XYCoord x1 y1 = projectionForward p1 lonlat
@ -317,7 +494,7 @@ mergeP p1 p2 t = Projection p pInv
-- | <<docs/gifs/doc_equirectangularP.gif>>
equirectangularP :: Projection
equirectangularP = Projection forward inverse
equirectangularP = Projection "equirectangular" forward inverse
where
forward (LonLat lam phi) = XYCoord ((lam+pi)/tau) ((phi+pi/2)/pi)
inverse (XYCoord x y) = LonLat xPi yPi
@ -327,7 +504,7 @@ equirectangularP = Projection forward inverse
-- | <<docs/gifs/doc_mercatorP.gif>>
mercatorP :: Projection
mercatorP = Projection forward inverse
mercatorP = Projection "mercator" forward inverse
where
forward (LonLat lam phi) =
XYCoord ((lam+pi)/tau)
@ -337,33 +514,65 @@ mercatorP = Projection forward inverse
xPi = fromToS (-pi) pi x
yPi = fromToS (-pi) pi y
thetas :: V.Vector Double
thetas = V.fromList $
map (find_theta_fast . fromIndex) [0 .. granularity]
granularity :: Int
granularity = 50000
toIndex :: Double -> Int
toIndex phi = round ((phi+halfPi)/pi * fromIntegral granularity)
fromIndex :: Int -> Double
fromIndex x = fromToS (-halfPi) halfPi (fromIntegral x / fromIntegral granularity)
granualize :: Double -> Double
granualize = fromIndex . toIndex
{-# INLINE mollweideP #-}
-- | <<docs/gifs/doc_mollweideP.gif>>
mollweideP :: Projection
mollweideP = Projection forward inverse
mollweideP = Projection "mollweide" forward inverse
where
forward (LonLat lam phi) =
forward (LonLat !lam !phi) =
XYCoord ((x+sqrt2*2)/(4*sqrt2)) ((y+sqrt2)/(2*sqrt2))
where
x = (2*sqrt2)/pi * lam * cos theta
y = sqrt2*sin theta
theta = find_theta 100
find_theta :: Int -> Double
find_theta 0 = phi
find_theta _ | abs phi == pi/2 = signum phi * pi/2
find_theta n =
let sub = find_theta (n-1)
in sub - (2*sub+sin (2*sub)-pi*sin phi)/(2+2*cos(2*sub))
theta = thetas V.! toIndex phi
-- find_theta :: Int -> Double
-- find_theta 0 = phi
-- find_theta _ | abs phi == pi/2 = signum phi * pi/2
-- find_theta n =
-- let !sub = find_theta (n-1)
-- in sub - (2*sub+sin (2*sub)-pi*sin phi)/(2+2*cos(2*sub))
inverse (XYCoord x' y') = LonLat lam phi
where
x = fromToS (-2*sqrt2) (2*sqrt2) x'
y = fromToS (-sqrt2) sqrt2 y'
theta = asin (y/sqrt2)
y = fromToS (-1) 1 y'
theta = granualize (asin y)
lam = pi*x/(2*sqrt2*cos theta)
phi = asin ((2*theta+sin(2*theta))/pi)
find_theta_fast :: Double -> Double
find_theta_fast !phi | abs phi == pi/2 = signum phi * halfPi
find_theta_fast !(D# phi) = go phi
where
!(D# pi#) = pi
go acc =
let c = cosDouble# (acc +## acc)
s = sinDouble# (acc +## acc)
next =
acc -##
(acc +## acc +## s -## pi# *## (sinDouble# phi))
/## (2.0## +## c +## c) in
if abs (D# (next -## acc)) < epsilon
then D# next
else go next
-- | <<docs/gifs/doc_hammerP.gif>>
hammerP :: Projection
hammerP = Projection forward inverse
hammerP = Projection "hammer" forward inverse
where
forward (LonLat lam phi) =
XYCoord ((x+sqrt2*2)/(4*sqrt2)) ((y+sqrt2)/(2*sqrt2))
@ -378,9 +587,20 @@ hammerP = Projection forward inverse
lam = 2 * atan2 (z*x) (2*(2*z**2-1))
phi = asin (z*y)
cylindricalEqualAreaP :: Double -> Projection
cylindricalEqualAreaP phi0 = Projection "lambert" forward inverse
where
cosPhi0 = cos phi0
forward (LonLat lam phi) =
XYCoord ((lam+pi)/tau) ((sin phi/cosPhi0+1/cosPhi0)/2/cosPhi0)
inverse (XYCoord x' y') = LonLat x (asin y / (asin cosPhi0 / halfPi))
where
x = fromToS (-pi) pi x'
y = fromToS (-1) 1 y' * cosPhi0
-- | <<docs/gifs/doc_lambertP.gif>>
lambertP :: Projection
lambertP = Projection forward inverse
lambertP = Projection "lambert" forward inverse
where
forward (LonLat lam phi) =
XYCoord ((lam+pi)/tau) ((sin phi+1)/2)
@ -391,7 +611,7 @@ lambertP = Projection forward inverse
-- | <<docs/gifs/doc_bottomleyP.gif>>
bottomleyP :: Double -> Projection
bottomleyP !phi_1 = Projection forward inverse
bottomleyP !phi_1 = Projection "bottomley" forward inverse
where
forward (LonLat lam phi) =
XYCoord ((x+pi)/tau) ((y+pi/2)/pi)
@ -413,7 +633,7 @@ bottomleyP !phi_1 = Projection forward inverse
-- | <<docs/gifs/doc_sinusoidalP.gif>>
sinusoidalP :: Projection
sinusoidalP = Projection forward inverse
sinusoidalP = Projection "sinusoidal" forward inverse
where
forward (LonLat lam phi) =
XYCoord ((x+pi)/tau) ((y+pi/2)/pi)
@ -427,7 +647,7 @@ sinusoidalP = Projection forward inverse
-- | <<docs/gifs/doc_wernerP.gif>>
wernerP :: Projection
wernerP = moveTopP 0.23 $ Projection forward inverse
wernerP = moveTopP 0.23 $ Projection "werner" forward inverse
where
forward (LonLat lam phi) =
XYCoord ((x+pi)/tau) ((y+pi/2)/pi)
@ -449,7 +669,8 @@ wernerP = moveTopP 0.23 $ Projection forward inverse
-- | <<docs/gifs/doc_bonneP.gif>>
bonneP :: Double -> Projection
bonneP 0 = sinusoidalP
bonneP phi_0 = moveTopP (-0.17*factor) $ scaleP 1 (fromToS 1 0.65 factor) $ Projection forward inverse
bonneP phi_0 = moveTopP (-0.17*factor) $ scaleP 1 (fromToS 1 0.65 factor) $
Projection "bonne" forward inverse
where
factor = sin phi_0 / sin (pi/4)
forward (LonLat lam phi ) = XYCoord ((x+pi)/tau) ((y+halfPi)/pi)
@ -471,7 +692,7 @@ bonneP phi_0 = moveTopP (-0.17*factor) $ scaleP 1 (fromToS 1 0.65 factor) $ Proj
-- | <<docs/gifs/doc_orthoP.gif>>
orthoP :: LonLat -> Projection
orthoP (LonLat lam_0 phi_0) = Projection forward inverse
orthoP (LonLat lam_0 phi_0) = Projection "ortho" forward inverse
where
forward (LonLat lam phi)
| cosc < 0 =
@ -503,7 +724,7 @@ orthoP (LonLat lam_0 phi_0) = Projection forward inverse
| otherwise = v
cassiniP :: Projection
cassiniP = Projection forward inverse
cassiniP = Projection "cassini" forward inverse
where
forward (LonLat lam phi) =
XYCoord ((asin (cos phi * sin lam)+halfPi)/pi) ((atan2 (tan phi) (cos lam)+pi)/tau)
@ -514,8 +735,9 @@ cassiniP = Projection forward inverse
lam = atan2 (tan x) (cos y)
phi = asin (sin y * cos x)
augustP :: Projection
augustP = scaleP 0.70 0.70 $ Projection forward inverse
augustP = scaleP 0.70 0.70 $ Projection "august" forward inverse
where
xHi = 16/3
xLo = -xHi
@ -550,7 +772,7 @@ augustP = scaleP 0.70 0.70 $ Projection forward inverse
-- | <<docs/gifs/doc_collignonP.gif>>
collignonP :: Projection
collignonP = Projection forward inverse
collignonP = Projection "collignon" forward inverse
where
yHi = sqrtPi
yLo = sqrtPi * (1 - sqrt2)
@ -571,7 +793,7 @@ collignonP = Projection forward inverse
-- | <<docs/gifs/doc_eckert1P.gif>>
eckert1P :: Projection
eckert1P = Projection forward inverse
eckert1P = Projection "eckert1" forward inverse
where
alpha = sqrt (8 / (3*pi))
yLo = -alpha * halfPi
@ -591,7 +813,7 @@ eckert1P = Projection forward inverse
-- | <<docs/gifs/doc_eckert3P.gif>>
eckert3P :: Projection
eckert3P = Projection forward inverse
eckert3P = Projection "eckert3" forward inverse
where
k = sqrt (pi * (4 + pi))
yLo = negate yHi
@ -611,7 +833,7 @@ eckert3P = Projection forward inverse
-- | <<docs/gifs/doc_eckert5P.gif>>
eckert5P :: Projection
eckert5P = Projection forward inverse
eckert5P = Projection "eckert5" forward inverse
where
k = sqrt (2 + pi)
yLo = negate yHi
@ -629,9 +851,10 @@ eckert5P = Projection forward inverse
lam = k * x / (1 + cos phi)
phi = y * k / 2
{-# INLINE faheyP #-}
-- | <<docs/gifs/doc_faheyP.gif>>
faheyP :: Projection
faheyP = Projection forward inverse
faheyP = Projection "fahey" forward inverse
where
faheyK = cos (toRads 35)
yLo = negate yHi
@ -651,8 +874,9 @@ faheyP = Projection forward inverse
lam = x / (faheyK * sqrt (1 - t*t))
phi = 2 * atan2 y (1 + faheyK)
{-# INLINE foucautP #-}
foucautP :: Projection
foucautP = Projection forward inverse
foucautP = Projection "foucaut" forward inverse
where
yLo = negate yHi
yHi = sqrtPi * tan (halfPi/2)
@ -673,8 +897,9 @@ foucautP = Projection forward inverse
phi = 2 * k
lam = x * sqrtPi / 2 / (cos phi * cosk * cosk)
{-# INLINE lagrangeP #-}
lagrangeP :: Projection
lagrangeP = Projection forward inverse
lagrangeP = Projection "lagrange" forward inverse
where
yLo = negate yHi
yHi = 2
@ -691,7 +916,9 @@ lagrangeP = Projection forward inverse
x = 2 * sin (lam*n) / c
y = (v - 1/v) /c
inverse (XYCoord x' y')
| abs (abs y'-1) < epsilon = LonLat 0 (signum y * halfPi)
| abs (y'-1) < epsilon
-- = LonLat 0 (signum y * halfPi)
= LonLat 100 100
| otherwise = LonLat lam phi
where
x = fromToS xLo xHi x' / 2

View file

@ -38,6 +38,10 @@ data CacheMap = CacheMap !(Map.Map DynamicName CacheMap) !(Map.Map DynamicName D
emptyCacheMap :: CacheMap
emptyCacheMap = CacheMap Map.empty Map.empty
-- sizeCacheMap :: CacheMap -> Int
-- sizeCacheMap (CacheMap sub vals) =
-- sum (map sizeCacheMap (Map.elems sub)) + Map.size vals
cacheMapLookup :: [DynamicName] -> CacheMap -> Maybe Dynamic
cacheMapLookup [] _ = Nothing
cacheMapLookup [k] (CacheMap _ vals) = Map.lookup k vals

View file

@ -1,17 +1,47 @@
module Reanimate.Parameters
( pFPS
( Raster(..)
, Width
, Height
, FPS
, pRaster
, pFPS
, pWidth
, pHeight
, pNoExternals
, pRootDirectory
, setRaster
, setFPS
, setWidth
, setHeight
, setNoExternals
, setRootDirectory
) where
import System.IO.Unsafe
import Data.IORef
import Reanimate.Render
type Width = Int
type Height = Int
type FPS = Int
data Raster
= RasterNone
| RasterAuto
| RasterInkscape
| RasterRSvg
| RasterConvert
deriving (Show)
{-# NOINLINE pRasterRef #-}
pRasterRef :: IORef Raster
pRasterRef = unsafePerformIO (newIORef RasterNone)
{-# NOINLINE pRaster #-}
pRaster :: Raster
pRaster = unsafePerformIO (readIORef pRasterRef)
setRaster :: Raster -> IO ()
setRaster = writeIORef pRasterRef
{-# NOINLINE pFPSRef #-}
pFPSRef :: IORef FPS
@ -58,3 +88,13 @@ pNoExternals = unsafePerformIO (readIORef pNoExternalsRef)
setNoExternals :: Bool -> IO ()
setNoExternals = writeIORef pNoExternalsRef
{-# NOINLINE pRootDirectoryRef #-}
pRootDirectoryRef :: IORef FilePath
pRootDirectoryRef = unsafePerformIO (newIORef (error "root directory not set"))
{-# NOINLINE pRootDirectory #-}
pRootDirectory :: FilePath
pRootDirectory = unsafePerformIO (readIORef pRootDirectoryRef)
setRootDirectory :: FilePath -> IO ()
setRootDirectory = writeIORef pRootDirectoryRef

View file

@ -1,5 +1,7 @@
module Reanimate.Raster
( embedImage
( mkImage
, cacheImage
, embedImage
, embedDynamicImage
, embedPng
, raster
@ -35,6 +37,26 @@ import System.IO
import System.IO.Temp
import System.IO.Unsafe
-- FIXME: Embed the image data as inline base64 iff no raster engine is specified.
mkImage :: Double -> Double -> FilePath -> SVG
mkImage width height path = unsafePerformIO $ do
exists <- doesFileExist target
unless exists $ copyFile path target
return $ flipYAxis $ ImageTree $ defaultSvg
& Svg.imageWidth .~ Svg.Num width
& Svg.imageHeight .~ Svg.Num height
& Svg.imageHref .~ ("file://"++target)
& Svg.imageCornerUpperLeft .~ (Svg.Num (-width/2), Svg.Num (-height/2))
& Svg.imageAspectRatio .~ Svg.PreserveAspectRatio False Svg.AlignNone Nothing
where
target = pRootDirectory </> show hashPath <.> takeExtension path
hashPath = hash path
cacheImage :: (PngSavable pixel, Hashable a) => a -> Image pixel -> FilePath
cacheImage key gen = unsafePerformIO $ cacheFile template $ \path ->
writePng path gen
where
template = show (hash key) <.> "png"
{-# INLINE embedImage #-}
embedImage :: PngSavable a => Image a -> Tree

View file

@ -10,16 +10,18 @@ module Reanimate.Render
import Control.Concurrent
import Control.Exception
import Control.Monad (forM_, void)
import Control.Monad (forM_, void, unless, forever)
import qualified Data.Text as T
import qualified Data.Text.IO as T
import Data.Time
import Graphics.SvgTree (Number (..))
import Numeric
import Reanimate.Animation
import Reanimate.Misc
import System.FilePath ((</>))
import System.FilePath (replaceExtension)
import Reanimate.Parameters
import System.Exit
import System.FilePath ((</>))
import System.FilePath (replaceExtension)
import System.IO
import Text.Printf (printf)
@ -71,21 +73,9 @@ filterFrameList seen nthFrame nFrames =
where
isSeen x = any (\y -> x `mod` y == 0) seen
data Raster
= RasterNone
| RasterAuto
| RasterInkscape
| RasterRSvg
| RasterConvert
deriving (Show)
data Format = RenderMp4 | RenderGif | RenderWebm
deriving (Show)
type Width = Int
type Height = Int
type FPS = Int
render :: Animation
-> FilePath
-> Raster
@ -131,18 +121,35 @@ render ani target raster format width height fps = do
generateFrames :: Raster -> Animation -> Width -> Height -> FPS -> (FilePath -> IO a) -> IO a
generateFrames raster ani width_ height_ rate action = withTempDir $ \tmp -> do
setRootDirectory tmp
done <- newMVar (0::Int)
let frameName nth = tmp </> printf nameTemplate nth
putStr $ "\r0/" ++ show frameCount
putStr $ "\rFrames rendered: 0/" ++ show frameCount ++ "\27[K\r"
hFlush stdout
handle h $ concurrentForM_ frames $ \n -> do
writeFile (frameName n) $ renderSvg width height $ nthFrame n
applyRaster raster (frameName n)
modifyMVar_ done $ \nDone -> do
putStr $ "\r" ++ show (nDone+1) ++ "/" ++ show frameCount
hFlush stdout
return (nDone+1)
putStrLn "\n"
start <- getCurrentTime
let statusPrinter = forever $ do
nDone <- readMVar done
now <- getCurrentTime
let spent = diffUTCTime now start
remaining = (spent / (fromIntegral nDone / fromIntegral frameCount)) - spent
putStr $ "\rFrames rendered: " ++ show nDone ++ "/" ++ show frameCount
putStr $ ", time spent: " ++ ppDiff spent
unless (nDone==0) $ do
putStr $ ", time remaining: " ++ ppDiff remaining
putStr $ ", total time: " ++ ppDiff (remaining+spent)
putStr $ "\27[K\r"
hFlush stdout
threadDelay 1000000
withBackgroundThread statusPrinter $ do
handle h $ concurrentForM_ frames $ \n -> do
writeFile (frameName n) $ renderSvg width height $ nthFrame n
applyRaster raster (frameName n)
modifyMVar_ done $ \nDone -> return (nDone+1)
now <- getCurrentTime
let spent = diffUTCTime now start
putStr $ "\rFrames rendered: " ++ show frameCount ++ "/" ++ show frameCount
putStr $ ", time spent: " ++ ppDiff spent
putStr $ "\27[K\n"
action (tmp </> rasterTemplate raster)
where
width = Just $ Px $ fromIntegral width_
@ -152,12 +159,25 @@ generateFrames raster ani width_ height_ rate action = withTempDir $ \tmp -> do
\Hit ctrl-c again to abort."
return ()
h other = throwIO other
frames = [0..frameCount-1]
-- frames = [0..frameCount-1]
frames = frameOrder rate frameCount
nthFrame nth = frameAt (recip (fromIntegral rate) * fromIntegral nth) ani
frameCount = round (duration ani * fromIntegral rate) :: Int
nameTemplate :: String
nameTemplate = "render-%05d.svg"
withBackgroundThread :: IO () -> IO a -> IO a
withBackgroundThread t = bracket (forkIO t) killThread . const
ppDiff :: NominalDiffTime -> String
ppDiff diff
| hours == 0 && mins == 0 = show secs ++ "s"
| hours == 0 = printf "%.2d:%.2d" mins secs
| otherwise = printf "%.2d:%.2d:%.2d" hours mins secs
where
(osecs, secs) = round diff `divMod` (60::Int)
(hours, mins) = osecs `divMod` 60
rasterTemplate :: Raster -> String
rasterTemplate RasterNone = "render-%05d.svg"
rasterTemplate _ = "render-%05d.png"

View file

@ -44,7 +44,8 @@ linePoints = worker zero
case x of
LineMove to -> worker to xs
-- LineDraw to -> from:to:worker to xs
-- FIXME: Use approximation from Geom2D.Bezier
LineBezier [p] ->
p : worker p xs
LineBezier ctrl -> -- approximation
[ last (partialBezierPoints (from:ctrl) 0 (recip chunks*i)) | i <- [0..chunks]] ++
worker (last ctrl) xs

View file

@ -0,0 +1,40 @@
#!/usr/bin/env stack
-- stack --resolver lts-13.14 runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RecordWildCards #-}
module EndScene (endScene) where
import Control.Lens ()
import Control.Monad
import qualified Data.ByteString as BS
import qualified Data.Map as Map
import Data.Monoid
import qualified Data.Text as T
import Codec.Picture
import Codec.Picture.Jpg
import Codec.Picture.Types
import Data.Maybe
import Data.Word
import Graphics.SvgTree hiding (Image, imageHeight, imageWidth)
import Graphics.SvgTree.Memo
import Numeric
import Reanimate
import Reanimate.Animation
import Reanimate.ColorMap
import Reanimate.ColorSpace
import Reanimate.Builtin.Images
import Reanimate.Constants
import Reanimate.Effect
import Reanimate.LaTeX
import Reanimate.Raster
import Reanimate.Scene
import Reanimate.Signal
import Reanimate.Svg
import System.IO.Unsafe
endScene :: Animation
endScene = mkAnimation 10 $ const $
mkGroup
[ mkBackground "black"
, scale 0.5 $ githubIcon ]

View file

@ -0,0 +1,54 @@
Globes give a good sense of the shape and size of countries but also limit your
view to less than half the world at any one time.
If we unwrap a globe and lay it flat, the result is a projection that neither
preserves the shapes or relative sizes of countries.
15s: unfolded
We can improve this by stretching the projection, giving true land sizes but
still distorting shapes.
25s: mercator
Mercator distorts everything and is, unfortunately, the projection used for
Google Maps.
30s: mollweide
Mollweide is a vast improvement over Mercator for most uses. Like Lambert,
it shows the true size of countries and continents at the cost of distorting
shapes. But this distortion is minimized by projecting on to an oval rather
than a rectangle.
50s: werner
Werner's projection isn't in much use anymore but can be found in certain
historical documents from the 15 hundreds. Also, you never know when you'll need
a heart-shaped map. Keep Werner in mind when wooing a carthographer.
65s: Collignon
Collignon might look strange but still manages to preserve relative sizes.
70s: Eckert
Eckert projections are some of my favorites. Each of Eckert's six projections
have a different trade-off between the accuracy of shapes and sizes but I find
all of them to be visually pleasing.
85s: Fahey
Some people, not to be confused with flat-earthers, know that the world is round
and therefore a round projection is needed.
90s: Foucaut
Other projections may appear of little practical use but they often have
interesting mathematical properties that can be important, especially when comparing
smaller map segments.
105: GitHub
This animation was created with Haskell and Blender. The full source code is
available on GitHub.

View file

@ -0,0 +1 @@
87949a7c6e3899b215eac29f902bbfb62d1c2203

View file

@ -0,0 +1 @@
5e56c0fadb910df0aef60061f69b5ba522abbd7a

View file

@ -0,0 +1 @@
757640d5d70bf61aca382de96203686240d2ad35

Binary file not shown.

After

Width:  |  Height:  |  Size: 55 KiB

View file

@ -0,0 +1 @@
df20addd14cbdc3289de40959a472ce7394c75e6

View file

@ -0,0 +1 @@
d705ac84752ae6af5963a76f65d99b2fb891a77e

View file

@ -1 +0,0 @@
088ba82a112f25607c756814500ad4246d7ced4c

View file

@ -1,21 +1,21 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE ApplicativeDo #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE BangPatterns #-}
module Main(main) where
import Codec.Picture
import Codec.Picture.Jpg
import Codec.Picture.Types
import Control.Lens ((^.), (%~), (&))
import Control.Exception
import Control.Lens ((%~), (&), (^.))
import Control.Monad
import Control.Monad.ST
import Data.Aeson
import qualified Data.ByteString as BS
import qualified Data.Sequence as Seq
import Data.Char
import Data.Foldable
import Data.Geospatial hiding (LonLat)
@ -24,18 +24,26 @@ import qualified Data.LineString as Line
import Data.Map (Map)
import qualified Data.Map as Map
import Data.Maybe
import qualified Data.Sequence as Seq
import Data.String.Here
import qualified Data.Text as T
import Debug.Trace
import Graphics.SvgTree (ElementRef (..), PathCommand (..),
Tree (None))
import Reanimate
import Reanimate.Animation
import Reanimate.Blender
import Reanimate.GeoProjection
import Reanimate.Scene
import Reanimate.Memo
import Reanimate.Raster
import Reanimate.Scene
import Reanimate.Transition
import Reanimate.Builtin.Flip
import System.IO.Unsafe
import EndScene
quick = False
{-
Show the earth.
@ -45,57 +53,66 @@ import System.IO.Unsafe
earth :: FilePath
-- earth = "earth-low.jpg"
-- earth = "earth-mid.jpg"
earth = "earth-high.jpg"
-- earth = "earth-high.jpg"
-- earth = "earth-extreme.jpg"
earth = "earth-1440.png"
earthMax :: FilePath
earthMax = "earth-max.jpg"
main :: IO ()
main = reanimate mainScene
main = do
-- putStrLn "Loading earth"
-- _ <- evaluate equirectangular
-- print (imageWidth equirectangular, imageHeight equirectangular)
-- putStrLn "Interpolating"
-- evaluate $ interpFastP equirectangular equirectangularP augustP 0.1
-- return ()
-- reanimate testScene
reanimate $
parA (staticFrame 1 $ mkBackground "darkgrey") $
overlapT 2 (signalT (curveS 2) flipTransition)
mainScene
endScene
testScene :: Animation
testScene = sceneAnimation $ do
bg <- newSpriteSVG $ mkBackground "white"
spriteZ bg (-1)
-- play $ animate $ const $ scale (1/halfPi * (4/4.5)) $ scaleToSize screenWidth screenHeight $
-- embedImage $ project src (orthoP $ LonLat 0 0)
bend <- newVar 1
let LonLat lam phi = newzealandLonLat
rotX <- newVar 0
rotY <- newVar 0
draw <- newVar 0
_ <- newSprite $ do
getBend <- unVar bend
getRotX <- unVar rotX
getRotY <- unVar rotY
getDraw <- unVar draw
return $
scale (fromToS 1 ((9*pi)/16) getBend) $
mkGroup
[ blender (script earth getBend getRotX getRotY)
-- , lowerTransformations $ scale (1/halfPi * (4/4.5)) $
-- grid $ orthoP $ LonLat (getRotY) (getRotX)
newSpriteSVG $ mkBackground "white"
grid <- newSpriteSVG $ mkGroup
[ withStrokeWidth (defaultStrokeWidth*0.2) $ withStrokeColor "black" $ mkGroup
[ mkLine (0,screenHeight) (0,-screenHeight)
, mkLine (screenHeight/2,screenHeight) (screenHeight/2,-screenHeight)
, mkLine (-screenHeight/2,screenHeight) (-screenHeight/2,-screenHeight)
]
]
spriteZ grid 1
Globe{..} <- newGlobe
let pos = brazilLonLat
writeVar globePosition pos
tweenVar globePosition 1 $ \v -> fromToLonLat v (LonLat 0 0) . curveS 2
tweenVar globeBend 2 $ \v -> fromToS v 0 . curveS 2
wait 1
-- destroySprite globeSprite
-- newSpriteSVG $
-- mkGroup
-- [ scaleToSize screenWidth screenHeight $
-- embedImage $ project src (orthoP pos)
-- ]
-- wait 1
-- tweenVar bend 1 $ \v -> fromToS v 0 . curveS 2
fork $ tweenVar rotY 1 $ \v -> fromToS v (lam) . curveS 2
tweenVar rotX 1 $ \v -> fromToS v (phi) . curveS 2
-- wait 1
fork $ tweenVar rotY 2 $ \v -> fromToS v 0 . curveS 2
fork $ tweenVar rotX 2 $ \v -> fromToS v 0 . curveS 2
wait 0.3
tweenVar bend 2 $ \v -> fromToS v 0 . curveS 2
where
src = equirectangular
mainScene :: Animation
mainScene = seq equirectangular $ -- takeA 5 $ dropA 19 $
mainScene = seq equirectangular $ -- takeA 10 $ dropA 21 $
mapA (withStrokeColor "black") $ sceneAnimation $ do
bg <- newSpriteSVG $ mkBackground "white"
spriteZ bg (-1)
let offset = translate 0 (-screenHeight/2 * 0.25)
orthoScale = 0.35
largeScale = 0.50
let offset = translate 0 (-screenHeight/2 * 0.20)
orthoScale = 0.50
largeScale = 0.75
Globe{..} <- newGlobe
morph <- newVar 0
@ -104,18 +121,24 @@ mainScene = seq equirectangular $ -- takeA 5 $ dropA 19 $
spriteModify globeSprite $ do
s <- unVar mapScale
pure $ \(svg, z) -> (scale s svg, z)
pure $ \(svg, z) -> (scale orthoScale svg, z)
spriteMap globeSprite (offset)
-- destroySprite globeSprite
let addRegion x y s proj lonlat@(LonLat lam phi) = do
let addRegion x y s proj lonlat@(LonLat lam phi) label = do
let idName = filter isAlphaNum $ show lonlat
outline = lowerTransformations $ scale orthoScale $
proj (orthoP lonlat)
region1 <- newSpriteSVG outline
destroySprite region1
srcWidth = imageWidth equirectangularMax
srcHeight = imageHeight equirectangularMax
subImg = trace ("subImg for: " ++ T.unpack label) $
convertRGBA8 $ rasterSized srcWidth srcHeight $ mkGroup
[ mkGroup []
, mkClipPath idName $
clipSvg
, withClipPathRef (Ref idName) $
scaleToSize screenWidth screenHeight $
embedImage equirectangularMax]
clipSvg = removeGroups $
lowerTransformations $ proj equirectangularP
spriteZ region1 2
move <- newVar 0
region1Shadow <- newSprite $ do
~(from, to) <- unVar projs
@ -125,15 +148,21 @@ mainScene = seq equirectangular $ -- takeA 5 $ dropA 19 $
t <- unVar move
pure $
let srcWidth = imageWidth equirectangularExtreme
srcHeight = imageHeight equirectangularExtreme
!subImg = convertRGBA8 $ rasterSized srcWidth srcHeight $ mkGroup
[ mkGroup []
, mkClipPath idName $
clipSvg
, withClipPathRef (Ref idName) $
scaleToSize screenWidth screenHeight $
embedImage $ project equirectangularExtreme equirectangularP]
let
(bx,by,bw,bh) = boundingBox $ proj (from lonlat)
fx = (bx+screenWidth/2)/screenWidth
fy = (by+screenHeight/2)/screenHeight
fw = bw/screenWidth
fh = bh/screenHeight
(bx',by',bw',bh') = boundingBox $ proj to
fx' = (bx'+screenWidth/2)/screenWidth
fy' = (by'+screenHeight/2)/screenHeight
fw' = bw'/screenWidth
fh' = bh'/screenHeight
imgKey = (projectionLabel (from lonlat), projectionLabel to, lonlat, earthMax, m)
imgFile = cacheImage imgKey $
trace ("interp for:" ++ show imgKey) $
interpBBP subImg (from lonlat) to (fx,fy,fw,fh) (fx',fy',fw',fh') m
setPos =
translate (fromToS 0 x $ curveS 2 t)
(fromToS 0 y $ curveS 2 t) .
@ -141,36 +170,50 @@ mainScene = seq equirectangular $ -- takeA 5 $ dropA 19 $
scale (fromToS 1 s $ curveS 2 t)
posSvg =
lowerTransformations $ proj (mergeP (from lonlat) to m)
clipSvg = removeGroups $
lowerTransformations $ proj equirectangularP in
finalSvg =
lowerTransformations $ proj to
toCenter = if projectionLabel (from lonlat) == "ortho"
then centerWithDelta m finalSvg
else centerWithDelta 1 posSvg
-- targetXY
-- translate ((-x-w/2)*d) ((-y-h/2)*d) t
in
mkGroup
[ mkGroup []
, setPos $ scale relScale $ centerWithDelta 1 posSvg $
scaleToSize screenWidth screenHeight $
embedImage $ interpP subImg (from lonlat) to m
, lowerTransformations $ setPos $ scale orthoScale $ center $
proj (orthoP lonlat)
, setPos $ scale relScale $ toCenter $
mkImage screenWidth screenHeight imgFile
, withStrokeWidth (fromToS 0 (defaultStrokeWidth*0.3) t) $
lowerTransformations $ setPos $ scale orthoScale $
(proj (orthoP lonlat))
]
-- destroySprite region1Shadow
fork $ tweenVar move 1 $ \v -> fromToS v 1 . curveS 2
fork $ play $ (staticFrame 2 $
withStrokeWidth 0 $
translate 0 (-screenHeight*0.40) $
center $ latex label)
# applyE (overBeginning 0.2 fadeInE)
# applyE (overEnding 0.2 fadeOutE)
tweenVar globePosition 2 $ \v -> fromToLonLat v usaLonLat . curveS 2
fork $ addRegion (-5.5) 4 3 america usaLonLat
fork $ addRegion (-5.5) 3.5 2 america usaLonLat "USA"
tweenVar globePosition 2 $ \v -> fromToLonLat v brazilLonLat . curveS 2
fork $ addRegion (-6) 0 3 brazil brazilLonLat
fork $ addRegion (-6) (-0.5) 2 brazil brazilLonLat "Brazil"
tweenVar globePosition 2 $ \v -> fromToLonLat v ukLonLat . curveS 2
fork $ addRegion (-1) 4 5 uk ukLonLat
fork $ addRegion (-1) 4 4 uk ukLonLat "UK, scaled 200\\%"
tweenVar globePosition 2 $ \v -> fromToLonLat v germanyLonLat . curveS 2
fork $ addRegion 1 4 5 germany germanyLonLat
fork $ addRegion 1 4 4 germany germanyLonLat "Germany, scaled 200\\%"
tweenVar globePosition 2 $ \v -> fromToLonLat v ausLonLat . curveS 2
fork $ addRegion 6 4 3 australia ausLonLat
fork $ addRegion 6 3.5 2 australia ausLonLat "Australia"
tweenVar globePosition 2 $ \v -> fromToLonLat v newzealandLonLat . curveS 2
fork $ addRegion 6 0 4 newzealand newzealandLonLat
fork $ addRegion 6 (-0.5) 3 newzealand newzealandLonLat "New Zealand, scaled 150\\%"
fork $ tweenVar globePosition 3 $ \v -> fromToLonLat v (LonLat 0 0) . curveS 2
wait 1
@ -184,45 +227,65 @@ mainScene = seq equirectangular $ -- takeA 5 $ dropA 19 $
mapS <- newSprite $ do
~(from, to) <- unVar projs
m <- unVar morph
relScale <- unVar mapScale
pure $ lowerTransformations $ scale relScale $ mkGroup
[ mkGroup []
, scaleToSize screenWidth screenHeight $
embedImage $ interpP src (from (LonLat 0 0)) to m
, grid $ mergeP (from (LonLat 0 0)) to m
]
-- relScale <- unVar mapScale
t <- spriteT
pure $
let imgKey = (projectionLabel (from (LonLat 0 0)), projectionLabel to, earth, m
,"global"::String)
imgFile = cacheImage imgKey $
trace ("global map for: " ++ show imgKey) $
interpP src (from (LonLat 0 0)) to m
in lowerTransformations $ scale orthoScale $ mkGroup
[ mkGroup []
, if quick then None else
mkImage screenWidth screenHeight imgFile
, withStrokeWidth (fromToS 0 (defaultStrokeWidth*0.2) $ min t 1) $
grid $ mergeP (from (LonLat 0 0)) to m
]
spriteMap mapS offset
wait 1
play $ (staticFrame (projMorphT+projWaitT) $
withStrokeWidth 0 $
translate 0 (-screenHeight*0.43) $
center $ latex "Equirectangular")
# applyE (overBeginning 0.2 fadeInE)
# applyE (overEnding 0.2 fadeOutE)
let push proj = do
let push proj label = do
fork $ play $ (staticFrame (projMorphT+projWaitT) $
withStrokeWidth 0 $
translate 0 (-screenHeight*0.43) $
center $ latex label)
# applyE (overBeginning 0.2 fadeInE)
# applyE (overEnding 0.2 fadeOutE)
(_, prev) <- readVar projs
writeVar projs (const $ prev, proj)
writeVar morph 0
tweenVar morph 1 $ \v -> fromToS v 1 . curveS 2
wait 1
tweenVar morph projMorphT $ \v -> fromToS v 1 . curveS 2
wait projWaitT
push lambertP
push mercatorP
push mollweideP
push hammerP
push (bottomleyP $ 30/180*pi)
push sinusoidalP
push wernerP
push (bonneP $ 45/180*pi)
push augustP
push collignonP
push eckert1P
push eckert3P
push eckert5P
push faheyP
push foucautP
push lagrangeP
push lambertP "Lambert"
push mercatorP "Mercator"
push mollweideP "Mollweide"
push hammerP "Hammer"
push (bottomleyP $ 30/180*pi) "Bottomley"
push sinusoidalP "Sinusoidal"
push wernerP "Werner"
push (bonneP $ 45/180*pi) "Bonne"
push augustP "August"
push collignonP "Collignon"
push eckert1P "Eckert 1"
push eckert3P "Eckert 3"
push eckert5P "Eckert 5"
push faheyP "Fahey"
push foucautP "Foucaut"
push lagrangeP "Lagrange"
wait 5
where
src = equirectangular
waitT = 2
morphT = 2
projMorphT = 2
projWaitT = 3
centerDelta :: Double -> Tree -> Tree
centerDelta d t = translate ((-x-w/2)*d) ((-y-h/2)*d) t
@ -256,14 +319,14 @@ renderLabel label =
equirectangular :: Image PixelRGBA8
equirectangular = unsafePerformIO $ do
dat <- BS.readFile earth
case decodeJpeg dat of
case decodeImage dat of
Left err -> error err
Right img -> return $ convertRGBA8 img
equirectangularExtreme :: Image PixelRGBA8
equirectangularExtreme = unsafePerformIO $ do
dat <- BS.readFile "earth-extreme.jpg"
case decodeJpeg dat of
equirectangularMax :: Image PixelRGBA8
equirectangularMax = unsafePerformIO $ do
dat <- BS.readFile earthMax
case decodeImage dat of
Left err -> error err
Right img -> return $ convertRGBA8 img
@ -297,7 +360,6 @@ fetchCountry p checker =
(screenHeight)
$
translate (-1/2) (-1/2) $
withStrokeWidth strokeWidth $
withFillOpacity 0 $
mkGroup
@ -324,7 +386,7 @@ filterCountries = do
geo' = geo & geofeatures %~ Seq.filter fn
fn feat =
case Map.lookup "NAME" (feat ^. properties) of
Nothing -> False
Nothing -> False
Just name -> name `elem` goodNames
encodeFile "countries-limited.json" geo'
where
@ -389,13 +451,12 @@ grid p =
(screenHeight)
$
translate (-1/2) (-1/2) $
withStrokeWidth strokeWidth $
withFillOpacity 0 $
-- withFillColor "black" $
mkGroup
[ mkGroup []
, withStrokeColor "black" $
, withStrokeColor "grey" $
applyProjection p $
svgPointsToRadians $
pathify $ landGeo annotate
@ -455,7 +516,8 @@ newGlobe = do
~(LonLat lam phi) <- unVar pos
pure $
scale (fromToS 1 ((9*pi)/16) getBend) $
blender $ script earth getBend phi lam
mkImage screenWidth screenHeight $
blender' $ script earthMax getBend phi lam
return $ Globe globe pos bend
script :: FilePath -> Double -> Double -> Double -> T.Text
@ -482,7 +544,7 @@ cam.select_set(True)
focus_target.select_set(True)
bpy.ops.object.parent_set()
focus_target.rotation_euler = (${negate rotX}, 0, 0)
focus_target.rotation_euler = (${negate rotX}, ${rotY}, 0)
origin = bpy.data.objects['Cube']
@ -494,7 +556,7 @@ x = ${bend}
bpy.ops.mesh.primitive_plane_add()
plane = bpy.context.object
plane.scale = (16/2,${fromToS (9/2) 4 bend},1)
bpy.ops.object.shade_smooth()
#bpy.ops.object.shade_smooth()
bpy.context.object.active_material = bpy.data.materials['Material']
mat = bpy.context.object.active_material
@ -536,7 +598,7 @@ plane.select_set(True);
bpy.ops.object.origin_clear()
bpy.ops.object.origin_set(type='GEOMETRY_ORIGIN')
plane.rotation_euler = (0, ${negate rotY}, 0)
#plane.rotation_euler = (${negate rotX}, ${negate rotY}, 0)
scn = bpy.context.scene
@ -546,6 +608,8 @@ scn = bpy.context.scene
scn.view_settings.view_transform = 'Standard'
scn.render.film_transparent = True
scn.render.resolution_x = ${pWidth} #3200
scn.render.resolution_y = ${pHeight} #1800
bpy.ops.render.render( write_still=True )
|]