login

<     >

2023-10-06 15:22:28 (UTC-03:00)

Marcel Rodrigues <marcelgmr@gmail.com>

change directory structure

diff --git a/aa.lua b/aa.lua
deleted file mode 100644
index 1063328..0000000
--- a/aa.lua
+++ /dev/null
@@ -1,103 +0,0 @@
-local bit = require "bit"
-
-local surf = require "surf"
-
-local lshift, rshift =  bit.lshift,  bit.rshift
-
--- ByteMaps (and GIFs) have a maximum of 256 colors
--- we're anti-aliasing by mixing four pixels in one (k=4)
--- the maximum n where C(n+k-1, k) < 256 is 7
--- so we only need factorials up to 7+4-1=10
-local fact = {1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800}
-fact[0] = 1
-
--- binomial coefficient
--- returns zero if n < k
-local function C(n, k)
-    return n >= k and fact[n] / fact[k] / fact[n-k] or 0
-end
-
--- index of a combination with repetition
--- n must be ordered
-local function I(n, base)
-    local a, b, c, d = unpack(n)
-    local N = base + 3
-    return C(N, 4) - C(N-a-1, 4) - C(N-b-2, 3) - C(N-c-3, 2) - C(N-d-4, 1) - 1
-end
-
--- helper to generate combinations with repetition in lexicographical order
-local function next_mix(n, base)
-    for i = 4, 1, -1 do
-        if n[i] < base-1 then
-            n[i] = n[i] + 1
-            for j = i+1, 4 do
-                n[j] = n[j-1]
-            end
-            break
-        end
-    end
-end
-
-local function test_combinations()
-    local base = 7
-    local n = {0, 0, 0, 0}
-    local i = 0
-    while i < C(base+3, 4) do
-        assert(I(n, base) == i)
-        i = i + 1
-        next_mix(n, base)
-    end
-    print(i.." tests passed")
-end
-
-local function mixed_rgb(colors, n)
-    local c1 = colors[n[1]+1]
-    local c2 = colors[n[2]+1]
-    local c3 = colors[n[3]+1]
-    local c4 = colors[n[4]+1]
-    -- compute the average RGB of the four colors
-    local r = rshift(c1[1] + c2[1] + c3[1] + c4[1], 2)
-    local g = rshift(c1[2] + c2[2] + c3[2] + c4[2], 2)
-    local b = rshift(c1[3] + c2[3] + c3[3] + c4[3], 2)
-    return {r, g, b}
-end
-
-local function get_mixed_colors(colors)
-    local base = #colors
-    assert(base >= 2 and base <= 7)
-    local n = {0, 0, 0, 0}
-    local palette = {mixed_rgb(colors, n)}
-    for i = 2, C(base+3, 4) do
-        next_mix(n, base)
-        table.insert(palette, mixed_rgb(colors, n))
-    end
-    return palette
-end
-
-local function antialias(bytemap, base)
-    assert(base >= 2 and base <= 7)
-    -- output has half dimensions
-    -- discard last row/column for odd dimensions
-    local w = rshift(bytemap.w, 1)
-    local h = rshift(bytemap.h, 1)
-    local outmap = surf.new_bytemap(w, h)
-    local n = {} -- indices of colors to be mixed
-    local dx, dy -- coordinates in input bytemap (d is for "double")
-    for y = 0, h-1 do
-        dy = lshift(y, 1)
-        for x = 0, w-1 do
-            dx = lshift(x, 1)
-            n[1] = bytemap:pget(dx  , dy  )
-            n[2] = bytemap:pget(dx+1, dy  )
-            n[3] = bytemap:pget(dx  , dy+1)
-            n[4] = bytemap:pget(dx+1, dy+1)
-            assert(math.max(unpack(n)) < base)
-            table.sort(n)
-            local v = I(n, base)
-            outmap:pset(x, y, v)
-        end
-    end
-    return outmap
-end
-
-return {get_mixed_colors=get_mixed_colors, antialias=antialias}

diff --git a/aff.lua b/aff.lua
deleted file mode 100644
index 5ed98f6..0000000
--- a/aff.lua
+++ /dev/null
@@ -1,68 +0,0 @@
---[[
-Affine transformation matrix:
-  | a c e |
-  | b d f |
-  | 0 0 1 |
-]]
-
-local Affine = {}
-Affine.__index = Affine
-
-function Affine:reset()
-    self.a = 1; self.c = 0; self.e = 0
-    self.b = 0; self.d = 1; self.f = 0
-end
-
-function Affine:add_custom(a, b, c, d, e, f)
-    local na = a*self.a + c*self.b
-    local nb = b*self.a + d*self.b
-    local nc = a*self.c + c*self.d
-    local nd = b*self.c + d*self.d
-    local ne = a*self.e + c*self.f + e
-    local nf = b*self.e + d*self.f + f
-    self.a = na; self.c = nc; self.e = ne
-    self.b = nb; self.d = nd; self.f = nf
-end
-
-function Affine:add_squeeze(k)
-    self:add_custom(k, 0, 0, 1/k, 0, 0)
-end
-
-function Affine:add_scale(x, y)
-    y = y or x
-    self:add_custom(x, 0, 0, y, 0, 0)
-end
-
-function Affine:add_hshear(h)
-    self:add_custom(1, 0, h, 1, 0, 0)
-end
-
-function Affine:add_vshear(v)
-    self:add_custom(1, v, 0, 1, 0, 0)
-end
-
-function Affine:add_rotate(a)
-    local c, s = math.cos(a), math.sin(a)
-    self:add_custom(c, -s, s, c, 0, 0)
-end
-
-function Affine:add_translate(x, y)
-    self:add_custom(1, 0, 0, 1, x, y)
-end
-
-function Affine:apply(points)
-    local x, y
-    for i, p in ipairs(points) do
-        x = self.a*p[1] + self.c*p[2] + self.e
-        y = self.b*p[1] + self.d*p[2] + self.f
-        p[1], p[2] = x, y
-    end
-end
-
-local function new_affine()
-    local self = setmetatable({}, Affine)
-    self:reset()
-    return self
-end
-
-return {new_affine=new_affine}

diff --git a/anim.lua b/anim.lua
deleted file mode 100644
index 0e827ae..0000000
--- a/anim.lua
+++ /dev/null
@@ -1,10 +0,0 @@
-local anim = {}
-anim.surf = require "surf"
-anim.poly = require "poly"
-anim.aff  = require "aff"
-anim.map  = require "map"
-anim.aa   = require "aa"
-anim.ttf  = require "ttf"
-anim.shp  = require "shp"
-anim.gif  = require "gif"
-return anim

diff --git a/bio.lua b/bio.lua
deleted file mode 100644
index 4652f12..0000000
--- a/bio.lua
+++ /dev/null
@@ -1,174 +0,0 @@
--- binary input/output
-
-local ffi = require "ffi"
-local bit = require "bit"
-
--- ZigZag encoding to map signed to unsigned integers
-local function s2u(n)
-    return bit.bxor(bit.lshift(n, 1), bit.arshift(n, 31))
-end
-
--- ZigZag decoding to map unsigned back to signed integers
-local function u2s(n)
-    return bit.bxor(bit.rshift(n, 1), -bit.band(n, 1))
-end
-
-local function read_byte(fp)
-    return fp:read(1):byte(1)
-end
-
-local function read_leu16(fp)
-    return read_byte(fp) + bit.lshift(read_byte(fp), 8)
-end
-
-local function read_leu32(fp)
-    return read_leu16(fp) + bit.lshift(read_leu16(fp), 16)
-end
-
-local function read_beu16(fp)
-    return bit.lshift(read_byte(fp), 8) + read_byte(fp)
-end
-
-local function read_beu32(fp)
-    return bit.lshift(read_beu16(fp), 16) + read_beu16(fp)
-end
-
-local function signed(u, nbits)
-    local p = 2^(nbits-1)
-    if u >= p then u = u - p - p end
-    return u
-end
-
-local function read_lei16(fp)
-    return signed(read_leu16(fp), 16)
-end
-
-local function read_lei32(fp)
-    return signed(read_leu32(fp), 32)
-end
-
-local function read_bei16(fp)
-    return signed(read_beu16(fp), 16)
-end
-
-local function read_bei32(fp)
-    return signed(read_beu32(fp), 32)
-end
-
-local function read_led64(fp)
-    return ffi.cast("double *", ffi.new("char[8]", fp:read(8)))[0]
-end
-
-local function write_byte(fp, n)
-    fp:write(string.char(n))
-end
-
-local function write_beu16(fp, n)
-    write_byte(fp, bit.rshift(n, 8))
-    write_byte(fp, bit.band(n, 0xFF))
-end
-
-local function write_beu32(fp, n)
-    write_beu16(fp, bit.rshift(n, 16))
-    write_beu16(fp, bit.band(n, 0xFFFF))
-end
-
-local function write_lei16(fp, n)
-    write_byte(fp, bit.band(n, 0xFF))
-    write_byte(fp, bit.rshift(n, 8))
-end
-
-local function write_bei16(fp, n)
-    if n < 0 then
-        n = n + 0x10000
-    end
-    write_beu16(fp, n)
-end
-
-local RiceR = {}
-RiceR.__index = RiceR
-
-function RiceR:get_bit()
-    if self.n == 0 then
-        self.b = read_byte(self.fp)
-        self.n = 8
-    end
-    self.n = self.n - 1
-    return bit.band(bit.rshift(self.b, self.n), 1)
-end
-
-function RiceR:get_unsigned()
-    local q = 0
-    while self:get_bit() == 1 do
-        q = q + 1
-    end
-    local r = 0
-    for i = 1, self.k do
-        r = bit.bor(bit.lshift(r, 1), self:get_bit())
-    end
-    return bit.bor(bit.lshift(q, self.k), r)
-end
-
-function RiceR:get_signed()
-    return u2s(self:get_unsigned())
-end
-
-local function rice_r(fp, k)
-    local self = setmetatable({}, RiceR)
-    self.fp = fp    -- already opened file, read mode
-    self.k = k or 0 -- rice parameter
-    self.b = 0      -- value of last byte read
-    self.n = 0      -- number of bits available in self.b
-    return self
-end
-
-local RiceW = {}
-RiceW.__index = RiceW
-
-function RiceW:put_bit(b)
-    self.n = self.n - 1
-    self.b = bit.bor(self.b, bit.lshift(b, self.n))
-    if self.n == 0 then
-        self:flush()
-    end
-end
-
-function RiceW:flush()
-    if self.n ~= 8 then
-        write_byte(self.fp, self.b)
-        self.b = 0
-        self.n = 8
-    end
-end
-
-function RiceW:put_unsigned(n)
-    for i = 1, bit.rshift(n, self.k) do
-        self:put_bit(1)
-    end
-    self:put_bit(0)
-    for i = self.k-1, 0, -1 do
-        self:put_bit(bit.band(bit.rshift(n, i), 1))
-    end
-end
-
-function RiceW:put_signed(n)
-    self:put_unsigned(s2u(n))
-end
-
-local function rice_w(fp, k)
-    local self = setmetatable({}, RiceW)
-    self.fp = fp    -- already opened file, write mode
-    self.k = k or 0 -- rice parameter
-    self.b = 0      -- value of next byte to write
-    self.n = 8      -- number of bits available in self.b
-    return self
-end
-
-return {
-    read_byte=read_byte, read_leu16=read_leu16, read_leu32=read_leu32,
-    read_beu16=read_beu16, read_beu32=read_beu32, read_lei16=read_lei16,
-    read_lei32=read_lei32, read_bei16=read_bei16, read_bei32=read_bei32,
-    read_led64=read_led64, write_byte=write_byte, write_beu16=write_beu16,
-    write_beu32=write_beu32, write_lei16=write_lei16, write_bei16=write_bei16,
-    rice_r=rice_r, rice_w=rice_w
-}

diff --git a/gif.lua b/gif.lua
deleted file mode 100644
index 1409d4e..0000000
--- a/gif.lua
+++ /dev/null
@@ -1,232 +0,0 @@
-local ffi = require "ffi"
-local bit = require "bit"
-
-local bio = require "bio"
-local surf = require "surf"
-local lzw = require "lzw"
-
-local bnot = bit.bnot
-local bor, band = bit.bor, bit.band
-local lshift, rshift =  bit.lshift,  bit.rshift
-
--- == Encoder ==
-
-local write_num = bio.write_lei16
-
-local function write_nums(f, ...)
-    local nums = {...}
-    for i, n in pairs(nums) do
-        write_num(f, n)
-    end
-end
-
-local function get_depth(n)
-    local m, e = math.frexp(n-1)
-    return math.max(e, 1)
-end
-
-local GIFout = {}
-GIFout.__index = GIFout
-
-local function new_gif(f, w, h, colors)
-    if type(f) == "string" then f = io.open(f, "wb") end
-    local depth = get_depth(#colors)
-    local self = setmetatable({f=f, w=w, h=h, d=depth, gct=colors}, GIFout)
-    if depth == 1 then
-        self.back = surf.new_bitmap(w, h)
-    else
-        self.back = surf.new_bytemap(w, h)
-    end
-    f:write("GIF89a")
-    write_nums(f, w, h)
-    f:write(string.char(0xF0 + depth - 1, 0, 0)) -- FDSZ, BGINDEX, ASPECT
-    local i = 1
-    while i <= #colors do
-        f:write(string.char(unpack(colors[i])))
-        i = i + 1
-    end
-    while i <= 2^depth do             -- GCT size must be a power of two
-        f:write(string.char(0, 0, 0)) -- fill unused colors as black
-        i = i + 1
-    end
-    self.n = 0 -- # of frames added
-    return self
-end
-
-function GIFout:set_loop(n)
-    n = n or 0
-    self.f:write("!")
-    self.f:write(string.char(0xFF, 0x0B))
-    self.f:write("NETSCAPE2.0")
-    self.f:write(string.char(0x03, 0x01))
-    write_num(self.f, n)
-    self.f:write(string.char(0))
-end
-
-function GIFout:set_delay(d)
-    self.f:write("!")
-    self.f:write(string.char(0xF9, 0x04, 0x04))
-    write_num(self.f, d)
-    self.f:write(string.char(0, 0))
-end
-
-function GIFout:get_bbox(frame)
-    local w, h = self.w, self.h
-    if self.n == 0 then return 0, 0, w, h end
-    local xmin, ymin = w, h
-    local xmax, ymax = 0, 0
-    local back = self.back
-    for y = 0, h-1 do
-        for x = 0, w-1 do
-            if frame:pget(x, y) ~= back:pget(x, y) then
-                if x < xmin then xmin = x end
-                if y < ymin then ymin = y end
-                if x > xmax then xmax = x end
-                if y > ymax then ymax = y end
-            end
-        end
-    end
-    if xmin == w or ymin == h then return 0, 0, 1, 1 end
-    return xmin, ymin, xmax-xmin+1, ymax-ymin+1
-end
-
-function GIFout:put_image(frame, x, y, w, h)
-    self.f:write(",")
-    write_nums(self.f, x, y, w, h)
-    self.f:write(string.char(0))
-    lzw.encode(self.f, self.d, frame, x, y, w, h) -- IP (Appendix F)
-end
-
-function GIFout:add_frame(frame, delay)
-    if delay then self:set_delay(delay) end
-    local x, y, w, h = self:get_bbox(frame)
-    self:put_image(frame, x, y, w, h)
-    self.n = self.n + 1
-    self.back:blit(x, y, frame, x, y, w, h)
-end
-
-function GIFout:close()
-    self.f:write(";")
-    self.f:close()
-end
-
--- == Decoder ==
-
-local read_num = bio.read_lei16
-
-local GIFin = {}
-GIFin.__index = GIFin
-
-local function open_gif(f)
-    if type(f) == "string" then f = io.open(f, "rb") end
-    assert(f:read(6) == "GIF89a", "invalid signature")
-    local w, h = read_num(f), read_num(f)
-    local fdsz = bio.read_byte(f)
-    local has_gct = rshift(fdsz, 7) == 1
-    local d = band(fdsz, 7) + 1
-    local bg = bio.read_byte(f)
-    f:seek("cur", 1) -- ASPECT:u8
-    local self = setmetatable({f=f, w=w, h=h, d=d, bg=bg}, GIFin)
-    self.gct = has_gct and self:read_color_table(d) or {}
-    if d == 1 then
-        self.surf = surf.new_bitmap(w, h)
-    else
-        self.surf = surf.new_bytemap(w, h)
-    end
-    return self
-end
-
-function GIFin:read_color_table(d)
-    local ct = {}
-    local ncolors = lshift(1, d)
-    for i = 1, ncolors do
-        local r = bio.read_byte(self.f)
-        local g = bio.read_byte(self.f)
-        local b = bio.read_byte(self.f)
-        table.insert(ct, {r, g, b})
-    end
-    return ct
-end
-
-function GIFin:discard_sub_blocks()
-    repeat
-        local size = bio.read_byte(self.f)
-        self.f:seek("cur", size)
-    until size == 0
-end
-
-function GIFin:read_graphic_control_ext()
-    assert(bio.read_byte(self.f) == 0x04, "invalid GCE block size")
-    local rdit = bio.read_byte(self.f)
-    self.disposal = band(rshift(rdit, 2), 3)
-    self.input = band(rshift(rdit, 1), 1) == 1
-    local transp = band(rdit, 1) == 1
-    self.delay = read_num(self.f)
-    local tindex = bio.read_byte(self.f)
-    self.tindex = transp and tindex or -1
-    self.f:seek("cur", 1) -- end-of-block
-end
-
-function GIFin:read_application_ext()
-    assert(bio.read_byte(self.f) == 0x0B, "invalid APP block size")
-    local app_ip = self.f:read(8)
-    local app_auth_code = self.f:read(3)
-    if app_ip == "NETSCAPE" then
-        self.f:seek("cur", 2) -- always 0x03, 0x01
-        self.loop = read_num(self.f)
-        self.f:seek("cur", 1) -- end-of-block
-    else
-        self:discard_sub_blocks()
-    end
-end
-
-function GIFin:read_ext()
-    local label = bio.read_byte(self.f)
-    if label == 0xF9 then
-        self:read_graphic_control_ext()
-    elseif label == 0xFF then
-        self:read_application_ext()
-    else
-        error(("unknown extension: %02X"):format(label))
-    end
-end
-
-function GIFin:read_image()
-    local x, y = read_num(self.f), read_num(self.f)
-    local w, h = read_num(self.f), read_num(self.f)
-    local fisrz = bio.read_byte(self.f)
-    local has_lct = band(rshift(fisrz, 7), 1) == 1
-    --~ assert(not has_lct, "unsupported GIF feature: Local Color Table")
-    local interlace = band(rshift(fisrz, 6), 1) == 1
-    assert(not interlace, "unsupported GIF feature: Interlaced Frame")
-    local d = band(fisrz, 7) + 1
-    local lct = has_lct and self:read_color_table(d) or {}
-    d = has_lct and d or self.d
-    lzw.decode(self.f, d, self.surf, x, y, w, h) -- IP (Appendix F)
-end
-
-function GIFin:get_frame()
-    local sep = self.f:read(1)
-    while sep ~= "," do
-        if sep == ";" then
-            return nil
-        elseif sep == "!" then
-            self:read_ext()
-        else
-            return nil
-        end
-        sep = self.f:read(1)
-    end
-    self:read_image()
-    return self.surf
-end
-
-function GIFin:frames()
-    return function() return self:get_frame() end
-end
-
-function GIFin:close()
-    self.f:close()
-end
-
-return {new_gif=new_gif, open_gif=open_gif}

diff --git a/lib/anim.lua b/lib/anim.lua
new file mode 100644
index 0000000..fbb6a1c
--- /dev/null
+++ b/lib/anim.lua
@@ -0,0 +1,10 @@
+local anim = {}
+anim.surf = require "anim.surf"
+anim.poly = require "anim.poly"
+anim.aff  = require "anim.aff"
+anim.map  = require "anim.map"
+anim.aa   = require "anim.aa"
+anim.ttf  = require "anim.ttf"
+anim.shp  = require "anim.shp"
+anim.gif  = require "anim.gif"
+return anim

diff --git a/lib/anim/aa.lua b/lib/anim/aa.lua
new file mode 100644
index 0000000..7ec70ea
--- /dev/null
+++ b/lib/anim/aa.lua
@@ -0,0 +1,103 @@
+local bit = require "bit"
+
+local surf = require "anim.surf"
+
+local lshift, rshift =  bit.lshift,  bit.rshift
+
+-- ByteMaps (and GIFs) have a maximum of 256 colors
+-- we're anti-aliasing by mixing four pixels in one (k=4)
+-- the maximum n where C(n+k-1, k) < 256 is 7
+-- so we only need factorials up to 7+4-1=10
+local fact = {1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800}
+fact[0] = 1
+
+-- binomial coefficient
+-- returns zero if n < k
+local function C(n, k)
+    return n >= k and fact[n] / fact[k] / fact[n-k] or 0
+end
+
+-- index of a combination with repetition
+-- n must be ordered
+local function I(n, base)
+    local a, b, c, d = unpack(n)
+    local N = base + 3
+    return C(N, 4) - C(N-a-1, 4) - C(N-b-2, 3) - C(N-c-3, 2) - C(N-d-4, 1) - 1
+end
+
+-- helper to generate combinations with repetition in lexicographical order
+local function next_mix(n, base)
+    for i = 4, 1, -1 do
+        if n[i] < base-1 then
+            n[i] = n[i] + 1
+            for j = i+1, 4 do
+                n[j] = n[j-1]
+            end
+            break
+        end
+    end
+end
+
+local function test_combinations()
+    local base = 7
+    local n = {0, 0, 0, 0}
+    local i = 0
+    while i < C(base+3, 4) do
+        assert(I(n, base) == i)
+        i = i + 1
+        next_mix(n, base)
+    end
+    print(i.." tests passed")
+end
+
+local function mixed_rgb(colors, n)
+    local c1 = colors[n[1]+1]
+    local c2 = colors[n[2]+1]
+    local c3 = colors[n[3]+1]
+    local c4 = colors[n[4]+1]
+    -- compute the average RGB of the four colors
+    local r = rshift(c1[1] + c2[1] + c3[1] + c4[1], 2)
+    local g = rshift(c1[2] + c2[2] + c3[2] + c4[2], 2)
+    local b = rshift(c1[3] + c2[3] + c3[3] + c4[3], 2)
+    return {r, g, b}
+end
+
+local function get_mixed_colors(colors)
+    local base = #colors
+    assert(base >= 2 and base <= 7)
+    local n = {0, 0, 0, 0}
+    local palette = {mixed_rgb(colors, n)}
+    for i = 2, C(base+3, 4) do
+        next_mix(n, base)
+        table.insert(palette, mixed_rgb(colors, n))
+    end
+    return palette
+end
+
+local function antialias(bytemap, base)
+    assert(base >= 2 and base <= 7)
+    -- output has half dimensions
+    -- discard last row/column for odd dimensions
+    local w = rshift(bytemap.w, 1)
+    local h = rshift(bytemap.h, 1)
+    local outmap = surf.new_bytemap(w, h)
+    local n = {} -- indices of colors to be mixed
+    local dx, dy -- coordinates in input bytemap (d is for "double")
+    for y = 0, h-1 do
+        dy = lshift(y, 1)
+        for x = 0, w-1 do
+            dx = lshift(x, 1)
+            n[1] = bytemap:pget(dx  , dy  )
+            n[2] = bytemap:pget(dx+1, dy  )
+            n[3] = bytemap:pget(dx  , dy+1)
+            n[4] = bytemap:pget(dx+1, dy+1)
+            assert(math.max(unpack(n)) < base)
+            table.sort(n)
+            local v = I(n, base)
+            outmap:pset(x, y, v)
+        end
+    end
+    return outmap
+end
+
+return {get_mixed_colors=get_mixed_colors, antialias=antialias}

diff --git a/lib/anim/aff.lua b/lib/anim/aff.lua
new file mode 100644
index 0000000..5ed98f6
--- /dev/null
+++ b/lib/anim/aff.lua
@@ -0,0 +1,68 @@
+--[[
+Affine transformation matrix:
+  | a c e |
+  | b d f |
+  | 0 0 1 |
+]]
+
+local Affine = {}
+Affine.__index = Affine
+
+function Affine:reset()
+    self.a = 1; self.c = 0; self.e = 0
+    self.b = 0; self.d = 1; self.f = 0
+end
+
+function Affine:add_custom(a, b, c, d, e, f)
+    local na = a*self.a + c*self.b
+    local nb = b*self.a + d*self.b
+    local nc = a*self.c + c*self.d
+    local nd = b*self.c + d*self.d
+    local ne = a*self.e + c*self.f + e
+    local nf = b*self.e + d*self.f + f
+    self.a = na; self.c = nc; self.e = ne
+    self.b = nb; self.d = nd; self.f = nf
+end
+
+function Affine:add_squeeze(k)
+    self:add_custom(k, 0, 0, 1/k, 0, 0)
+end
+
+function Affine:add_scale(x, y)
+    y = y or x
+    self:add_custom(x, 0, 0, y, 0, 0)
+end
+
+function Affine:add_hshear(h)
+    self:add_custom(1, 0, h, 1, 0, 0)
+end
+
+function Affine:add_vshear(v)
+    self:add_custom(1, v, 0, 1, 0, 0)
+end
+
+function Affine:add_rotate(a)
+    local c, s = math.cos(a), math.sin(a)
+    self:add_custom(c, -s, s, c, 0, 0)
+end
+
+function Affine:add_translate(x, y)
+    self:add_custom(1, 0, 0, 1, x, y)
+end
+
+function Affine:apply(points)
+    local x, y
+    for i, p in ipairs(points) do
+        x = self.a*p[1] + self.c*p[2] + self.e
+        y = self.b*p[1] + self.d*p[2] + self.f
+        p[1], p[2] = x, y
+    end
+end
+
+local function new_affine()
+    local self = setmetatable({}, Affine)
+    self:reset()
+    return self
+end
+
+return {new_affine=new_affine}

diff --git a/lib/anim/bio.lua b/lib/anim/bio.lua
new file mode 100644
index 0000000..4652f12
--- /dev/null
+++ b/lib/anim/bio.lua
@@ -0,0 +1,174 @@
+-- binary input/output
+
+local ffi = require "ffi"
+local bit = require "bit"
+
+-- ZigZag encoding to map signed to unsigned integers
+local function s2u(n)
+    return bit.bxor(bit.lshift(n, 1), bit.arshift(n, 31))
+end
+
+-- ZigZag decoding to map unsigned back to signed integers
+local function u2s(n)
+    return bit.bxor(bit.rshift(n, 1), -bit.band(n, 1))
+end
+
+local function read_byte(fp)
+    return fp:read(1):byte(1)
+end
+
+local function read_leu16(fp)
+    return read_byte(fp) + bit.lshift(read_byte(fp), 8)
+end
+
+local function read_leu32(fp)
+    return read_leu16(fp) + bit.lshift(read_leu16(fp), 16)
+end
+
+local function read_beu16(fp)
+    return bit.lshift(read_byte(fp), 8) + read_byte(fp)
+end
+
+local function read_beu32(fp)
+    return bit.lshift(read_beu16(fp), 16) + read_beu16(fp)
+end
+
+local function signed(u, nbits)
+    local p = 2^(nbits-1)
+    if u >= p then u = u - p - p end
+    return u
+end
+
+local function read_lei16(fp)
+    return signed(read_leu16(fp), 16)
+end
+
+local function read_lei32(fp)
+    return signed(read_leu32(fp), 32)
+end
+
+local function read_bei16(fp)
+    return signed(read_beu16(fp), 16)
+end
+
+local function read_bei32(fp)
+    return signed(read_beu32(fp), 32)
+end
+
+local function read_led64(fp)
+    return ffi.cast("double *", ffi.new("char[8]", fp:read(8)))[0]
+end
+
+local function write_byte(fp, n)
+    fp:write(string.char(n))
+end
+
+local function write_beu16(fp, n)
+    write_byte(fp, bit.rshift(n, 8))
+    write_byte(fp, bit.band(n, 0xFF))
+end
+
+local function write_beu32(fp, n)
+    write_beu16(fp, bit.rshift(n, 16))
+    write_beu16(fp, bit.band(n, 0xFFFF))
+end
+
+local function write_lei16(fp, n)
+    write_byte(fp, bit.band(n, 0xFF))
+    write_byte(fp, bit.rshift(n, 8))
+end
+
+local function write_bei16(fp, n)
+    if n < 0 then
+        n = n + 0x10000
+    end
+    write_beu16(fp, n)
+end
+
+local RiceR = {}
+RiceR.__index = RiceR
+
+function RiceR:get_bit()
+    if self.n == 0 then
+        self.b = read_byte(self.fp)
+        self.n = 8
+    end
+    self.n = self.n - 1
+    return bit.band(bit.rshift(self.b, self.n), 1)
+end
+
+function RiceR:get_unsigned()
+    local q = 0
+    while self:get_bit() == 1 do
+        q = q + 1
+    end
+    local r = 0
+    for i = 1, self.k do
+        r = bit.bor(bit.lshift(r, 1), self:get_bit())
+    end
+    return bit.bor(bit.lshift(q, self.k), r)
+end
+
+function RiceR:get_signed()
+    return u2s(self:get_unsigned())
+end
+
+local function rice_r(fp, k)
+    local self = setmetatable({}, RiceR)
+    self.fp = fp    -- already opened file, read mode
+    self.k = k or 0 -- rice parameter
+    self.b = 0      -- value of last byte read
+    self.n = 0      -- number of bits available in self.b
+    return self
+end
+
+local RiceW = {}
+RiceW.__index = RiceW
+
+function RiceW:put_bit(b)
+    self.n = self.n - 1
+    self.b = bit.bor(self.b, bit.lshift(b, self.n))
+    if self.n == 0 then
+        self:flush()
+    end
+end
+
+function RiceW:flush()
+    if self.n ~= 8 then
+        write_byte(self.fp, self.b)
+        self.b = 0
+        self.n = 8
+    end
+end
+
+function RiceW:put_unsigned(n)
+    for i = 1, bit.rshift(n, self.k) do
+        self:put_bit(1)
+    end
+    self:put_bit(0)
+    for i = self.k-1, 0, -1 do
+        self:put_bit(bit.band(bit.rshift(n, i), 1))
+    end
+end
+
+function RiceW:put_signed(n)
+    self:put_unsigned(s2u(n))
+end
+
+local function rice_w(fp, k)
+    local self = setmetatable({}, RiceW)
+    self.fp = fp    -- already opened file, write mode
+    self.k = k or 0 -- rice parameter
+    self.b = 0      -- value of next byte to write
+    self.n = 8      -- number of bits available in self.b
+    return self
+end
+
+return {
+    read_byte=read_byte, read_leu16=read_leu16, read_leu32=read_leu32,
+    read_beu16=read_beu16, read_beu32=read_beu32, read_lei16=read_lei16,
+    read_lei32=read_lei32, read_bei16=read_bei16, read_bei32=read_bei32,
+    read_led64=read_led64, write_byte=write_byte, write_beu16=write_beu16,
+    write_beu32=write_beu32, write_lei16=write_lei16, write_bei16=write_bei16,
+    rice_r=rice_r, rice_w=rice_w
+}

diff --git a/lib/anim/gif.lua b/lib/anim/gif.lua
new file mode 100644
index 0000000..9a8b22a
--- /dev/null
+++ b/lib/anim/gif.lua
@@ -0,0 +1,232 @@
+local ffi = require "ffi"
+local bit = require "bit"
+
+local bio = require "anim.bio"
+local surf = require "anim.surf"
+local lzw = require "anim.lzw"
+
+local bnot = bit.bnot
+local bor, band = bit.bor, bit.band
+local lshift, rshift =  bit.lshift,  bit.rshift
+
+-- == Encoder ==
+
+local write_num = bio.write_lei16
+
+local function write_nums(f, ...)
+    local nums = {...}
+    for i, n in pairs(nums) do
+        write_num(f, n)
+    end
+end
+
+local function get_depth(n)
+    local m, e = math.frexp(n-1)
+    return math.max(e, 1)
+end
+
+local GIFout = {}
+GIFout.__index = GIFout
+
+local function new_gif(f, w, h, colors)
+    if type(f) == "string" then f = io.open(f, "wb") end
+    local depth = get_depth(#colors)
+    local self = setmetatable({f=f, w=w, h=h, d=depth, gct=colors}, GIFout)
+    if depth == 1 then
+        self.back = surf.new_bitmap(w, h)
+    else
+        self.back = surf.new_bytemap(w, h)
+    end
+    f:write("GIF89a")
+    write_nums(f, w, h)
+    f:write(string.char(0xF0 + depth - 1, 0, 0)) -- FDSZ, BGINDEX, ASPECT
+    local i = 1
+    while i <= #colors do
+        f:write(string.char(unpack(colors[i])))
+        i = i + 1
+    end
+    while i <= 2^depth do             -- GCT size must be a power of two
+        f:write(string.char(0, 0, 0)) -- fill unused colors as black
+        i = i + 1
+    end
+    self.n = 0 -- # of frames added
+    return self
+end
+
+function GIFout:set_loop(n)
+    n = n or 0
+    self.f:write("!")
+    self.f:write(string.char(0xFF, 0x0B))
+    self.f:write("NETSCAPE2.0")
+    self.f:write(string.char(0x03, 0x01))
+    write_num(self.f, n)
+    self.f:write(string.char(0))
+end
+
+function GIFout:set_delay(d)
+    self.f:write("!")
+    self.f:write(string.char(0xF9, 0x04, 0x04))
+    write_num(self.f, d)
+    self.f:write(string.char(0, 0))
+end
+
+function GIFout:get_bbox(frame)
+    local w, h = self.w, self.h
+    if self.n == 0 then return 0, 0, w, h end
+    local xmin, ymin = w, h
+    local xmax, ymax = 0, 0
+    local back = self.back
+    for y = 0, h-1 do
+        for x = 0, w-1 do
+            if frame:pget(x, y) ~= back:pget(x, y) then
+                if x < xmin then xmin = x end
+                if y < ymin then ymin = y end
+                if x > xmax then xmax = x end
+                if y > ymax then ymax = y end
+            end
+        end
+    end
+    if xmin == w or ymin == h then return 0, 0, 1, 1 end
+    return xmin, ymin, xmax-xmin+1, ymax-ymin+1
+end
+
+function GIFout:put_image(frame, x, y, w, h)
+    self.f:write(",")
+    write_nums(self.f, x, y, w, h)
+    self.f:write(string.char(0))
+    lzw.encode(self.f, self.d, frame, x, y, w, h) -- IP (Appendix F)
+end
+
+function GIFout:add_frame(frame, delay)
+    if delay then self:set_delay(delay) end
+    local x, y, w, h = self:get_bbox(frame)
+    self:put_image(frame, x, y, w, h)
+    self.n = self.n + 1
+    self.back:blit(x, y, frame, x, y, w, h)
+end
+
+function GIFout:close()
+    self.f:write(";")
+    self.f:close()
+end
+
+-- == Decoder ==
+
+local read_num = bio.read_lei16
+
+local GIFin = {}
+GIFin.__index = GIFin
+
+local function open_gif(f)
+    if type(f) == "string" then f = io.open(f, "rb") end
+    assert(f:read(6) == "GIF89a", "invalid signature")
+    local w, h = read_num(f), read_num(f)
+    local fdsz = bio.read_byte(f)
+    local has_gct = rshift(fdsz, 7) == 1
+    local d = band(fdsz, 7) + 1
+    local bg = bio.read_byte(f)
+    f:seek("cur", 1) -- ASPECT:u8
+    local self = setmetatable({f=f, w=w, h=h, d=d, bg=bg}, GIFin)
+    self.gct = has_gct and self:read_color_table(d) or {}
+    if d == 1 then
+        self.surf = surf.new_bitmap(w, h)
+    else
+        self.surf = surf.new_bytemap(w, h)
+    end
+    return self
+end
+
+function GIFin:read_color_table(d)
+    local ct = {}
+    local ncolors = lshift(1, d)
+    for i = 1, ncolors do
+        local r = bio.read_byte(self.f)
+        local g = bio.read_byte(self.f)
+        local b = bio.read_byte(self.f)
+        table.insert(ct, {r, g, b})
+    end
+    return ct
+end
+
+function GIFin:discard_sub_blocks()
+    repeat
+        local size = bio.read_byte(self.f)
+        self.f:seek("cur", size)
+    until size == 0
+end
+
+function GIFin:read_graphic_control_ext()
+    assert(bio.read_byte(self.f) == 0x04, "invalid GCE block size")
+    local rdit = bio.read_byte(self.f)
+    self.disposal = band(rshift(rdit, 2), 3)
+    self.input = band(rshift(rdit, 1), 1) == 1
+    local transp = band(rdit, 1) == 1
+    self.delay = read_num(self.f)
+    local tindex = bio.read_byte(self.f)
+    self.tindex = transp and tindex or -1
+    self.f:seek("cur", 1) -- end-of-block
+end
+
+function GIFin:read_application_ext()
+    assert(bio.read_byte(self.f) == 0x0B, "invalid APP block size")
+    local app_ip = self.f:read(8)
+    local app_auth_code = self.f:read(3)
+    if app_ip == "NETSCAPE" then
+        self.f:seek("cur", 2) -- always 0x03, 0x01
+        self.loop = read_num(self.f)
+        self.f:seek("cur", 1) -- end-of-block
+    else
+        self:discard_sub_blocks()
+    end
+end
+
+function GIFin:read_ext()
+    local label = bio.read_byte(self.f)
+    if label == 0xF9 then
+        self:read_graphic_control_ext()
+    elseif label == 0xFF then
+        self:read_application_ext()
+    else
+        error(("unknown extension: %02X"):format(label))
+    end
+end
+
+function GIFin:read_image()
+    local x, y = read_num(self.f), read_num(self.f)
+    local w, h = read_num(self.f), read_num(self.f)
+    local fisrz = bio.read_byte(self.f)
+    local has_lct = band(rshift(fisrz, 7), 1) == 1
+    --~ assert(not has_lct, "unsupported GIF feature: Local Color Table")
+    local interlace = band(rshift(fisrz, 6), 1) == 1
+    assert(not interlace, "unsupported GIF feature: Interlaced Frame")
+    local d = band(fisrz, 7) + 1
+    local lct = has_lct and self:read_color_table(d) or {}
+    d = has_lct and d or self.d
+    lzw.decode(self.f, d, self.surf, x, y, w, h) -- IP (Appendix F)
+end
+
+function GIFin:get_frame()
+    local sep = self.f:read(1)
+    while sep ~= "," do
+        if sep == ";" then
+            return nil
+        elseif sep == "!" then
+            self:read_ext()
+        else
+            return nil
+        end
+        sep = self.f:read(1)
+    end
+    self:read_image()
+    return self.surf
+end
+
+function GIFin:frames()
+    return function() return self:get_frame() end
+end
+
+function GIFin:close()
+    self.f:close()
+end
+
+return {new_gif=new_gif, open_gif=open_gif}

diff --git a/lib/anim/lzw.lua b/lib/anim/lzw.lua
new file mode 100644
index 0000000..377d23a
--- /dev/null
+++ b/lib/anim/lzw.lua
@@ -0,0 +1,215 @@
+local ffi = require "ffi"
+local bit = require "bit"
+
+local bnot = bit.bnot
+local bor, band = bit.bor, bit.band
+local lshift, rshift =  bit.lshift,  bit.rshift
+
+-- == Encoder ==
+
+local BUFout = {}
+BUFout.__index = BUFout
+
+local function new_buffer_out(f)
+    local self = setmetatable({f=f, offset=0, partial=0}, BUFout)
+    self.buf = ffi.new("char[255]")
+    return self
+end
+
+function BUFout:put_key(key, size)
+    local offset, partial = self.offset, self.partial
+    local f, buf = self.f, self.buf
+    local byte_offset, bit_offset = math.floor(offset / 8), offset % 8
+    partial = bor(partial, lshift(key, bit_offset))
+    local bits_to_write = bit_offset + size
+    while bits_to_write >= 8 do
+        buf[byte_offset] = band(partial, 0xFF)
+        byte_offset = byte_offset + 1
+        if byte_offset == 0xFF then -- flush
+            f:write(string.char(0xFF))
+            f:write(ffi.string(buf, 0xFF))
+            byte_offset = 0
+        end
+        partial = rshift(partial, 8)
+        bits_to_write = bits_to_write - 8
+    end
+    self.offset = (offset + size) % (0xFF * 8)
+    self.partial = partial
+end
+
+function BUFout:end_key()
+    local offset, partial = self.offset, self.partial
+    local f, buf = self.f, self.buf
+    local byte_offset = math.floor(offset / 8)
+    if offset % 8 ~= 0 then
+        buf[byte_offset] = band(partial, 0xFF)
+        byte_offset = byte_offset + 1
+    end
+    if byte_offset > 0 then
+        f:write(string.char(byte_offset))
+        f:write(ffi.string(buf, byte_offset))
+    end
+    f:write(string.char(0))
+    self.offset, self.partial = 0, 0
+end
+
+
+local function new_trie(degree)
+    local children = {}
+    for key = 0, degree-1 do
+        children[key] = {key=key, children={}}
+    end
+    return {children=children}
+end
+
+local function encode(f, d, s, x, y, w, h)
+    local buf = new_buffer_out(f)
+    local code_size = math.max(d, 2)
+    f:write(string.char(code_size))
+    local degree = 2 ^ code_size
+    local root = new_trie(degree)
+    local clear, stop = degree, degree + 1
+    local nkeys = degree + 2 -- skip clear code and stop code
+    local node = root
+    local key_size = code_size + 1
+    buf:put_key(clear, key_size)
+    for j = y, y+h-1 do
+        for i = x, x+w-1 do
+            local index = band(s:pget(i, j), degree-1)
+            local child = node.children[index]
+            if child ~= nil then
+                node = child
+            else
+                buf:put_key(node.key, key_size)
+                if nkeys < 0x1000 then
+                    if nkeys == 2 ^ key_size then
+                        key_size = key_size + 1
+                    end
+                    node.children[index] = {key=nkeys, children={}}
+                    nkeys = nkeys + 1
+                else
+                    buf:put_key(clear, key_size)
+                    root = new_trie(degree)
+                    node = root
+                    nkeys = degree + 2
+                    key_size = code_size + 1
+                end
+                node = root.children[index]
+            end
+        end
+    end
+    buf:put_key(node.key, key_size)
+    buf:put_key(stop, key_size)
+    buf:end_key()
+end
+
+-- == Decoder ==
+
+local BUFin = {}
+BUFin.__index = BUFin
+
+local function new_buffer_in(f)
+    local self = setmetatable({}, BUFin)
+    self.f = f      -- already opened file, read mode
+    self.s = 0      -- number of bytes available in block
+    self.b = 0      -- value of last byte read
+    self.n = 0      -- number of bits available in self.b
+    return self
+end
+
+function BUFin:get_key(size)
+    local key = 0
+    for i = 1, size do
+        if self.s == 0 then
+            self.s = self.f:read(1):byte(1)
+            assert(self.s > 0, "unexpected end-of-block")
+        end
+        if self.n == 0 then
+            self.b = self.f:read(1):byte(1)
+            self.n = 8
+            self.s = self.s - 1
+        end
+        key = bor(key, lshift(band(rshift(self.b, 8-self.n), 1), i-1))
+        self.n = self.n - 1
+    end
+    return key
+end
+
+local CodeTable = {}
+CodeTable.__index = CodeTable
+
+local function new_code_table(key_size)
+    local self = setmetatable({}, CodeTable)
+    self.len = lshift(1, key_size)
+    self.tab = {}
+    for key = 0, self.len+1 do
+        self.tab[key] = {length=1, prefix=0xFFF, suffix=key}
+    end
+    self.len = self.len + 2 -- clear & stop
+    return self
+end
+
+function CodeTable:add_entry(length, prefix, suffix)
+    self.tab[self.len] = {length=length, prefix=prefix, suffix=suffix}
+    self.len = self.len + 1
+    if band(self.len, self.len-1) == 0 then
+        return 1
+    end
+    return 0
+end
+
+local function decode(f, d, s, fx, fy, w, h)
+    local key_size = f:read(1):byte(1)
+    assert(key_size == math.max(d, 2), "invalid code size")
+    local buf = new_buffer_in(f)
+    local clear = lshift(key_size, 1)
+    local stop = clear + 1
+    key_size = key_size + 1
+    local init_key_size = key_size
+    local key = buf:get_key(key_size)
+    assert(key == clear, "expected clear code, got "..key)
+    local code_table, table_is_full, entry, str_len, ret
+    local frm_off = 0 -- pixels read
+    local frm_size = w * h
+    while frm_off < frm_size do
+        if key == clear then
+            key_size = init_key_size
+            code_table = new_code_table(key_size-1)
+            table_is_full = false
+        elseif not table_is_full then
+            ret = code_table:add_entry(str_len+1, key, entry.suffix)
+            if code_table.len == 0x1000 then
+                ret = 0
+                table_is_full = true
+            end
+        end
+        key = buf:get_key(key_size)
+        if key ~= clear then
+            if key == stop or key == 0x1000 then break end
+            if ret == 1 then key_size = key_size + 1 end
+            entry = code_table.tab[key]
+            str_len = entry.length
+            for i = 1, str_len do
+                local p = frm_off + entry.length - 1
+                local x = p % w
+                local y = math.floor(p / w)
+                s:pset(fx+x, fy+y, entry.suffix)
+                if entry.prefix == 0xFFF then
+                    break
+                else
+                    entry = code_table.tab[entry.prefix]
+                end
+            end
+            frm_off = frm_off + str_len
+            if key < code_table.len-1 and not table_is_full then
+                code_table.tab[code_table.len-1].suffix = entry.suffix
+            end
+        end
+    end
+    while buf.s > 0 do
+        f:seek("cur", buf.s)
+        buf.s = f:read(1):byte(1)
+    end
+end
+
+return {encode=encode, decode=decode}

diff --git a/lib/anim/map.lua b/lib/anim/map.lua
new file mode 100644
index 0000000..0fc6930
--- /dev/null
+++ b/lib/anim/map.lua
@@ -0,0 +1,310 @@
+local util = require "anim.util"
+local bio = require "anim.bio"
+
+local abs = math.abs
+local rad, deg = math.rad, math.deg
+local cos, sin, tan = math.cos, math.sin, math.tan
+local acos, asin, atan, atan2 = math.acos, math.asin, math.atan, math.atan2
+local sqrt = math.sqrt
+local pi, huge = math.pi, math.huge
+
+-- == Utilities ==
+
+-- region is a list of polygons in geographic coordinates.
+
+local function distance(lon1, lat1, lon2, lat2, r)
+    r = r or 6378137
+    local dlat = rad(lat2 - lat1)
+    local dlon = rad(lon2 - lon1)
+    lat1, lat2 = rad(lat1), rad(lat2)
+    local a1, a2, a, c
+    a1 = sin(dlat/2) * sin(dlat/2)
+    a2 = sin(dlon/2) * sin(dlon/2) * cos(lat1) * cos(lat2)
+    a = a1 + a2
+    c = 2 * atan2(sqrt(a), sqrt(1-a))
+    return r * c
+end
+
+local function bbox(polys)
+    local x0, y0, x1, y1 = huge, huge, -huge, -huge
+    for poly in util.func_iter(polys) do
+        for point in util.func_iter(poly) do
+            local x, y = unpack(point)
+            x0 = x < x0 and x or x0
+            y0 = y < y0 and y or y0
+            x1 = x > x1 and x or x1
+            y1 = y > y1 and y or y1
+        end
+    end
+    return {x0=x0, y0=y0, x1=x1, y1=y1}
+end
+
+local function centroid(region)
+    local epsilon = 1e-10
+    local bb = bbox(region)
+    local x0, y0, x1, y1 = bb.x0, bb.y0, bb.x1, bb.y1
+    local lon0 = (x0 + x1) / 2
+    local lat0 = (y0 + y1) / 2
+    local lon1, lat1
+    while true do
+        local prj = Proj("AzimuthalEqualArea", {lon0, lat0})
+        local cw = {}
+        for i, polygon in ipairs(region) do
+            local xys = {}
+            for j, point in ipairs(polygon) do
+                xys[j] = {prj:map(unpack(point))}
+            end
+            if xys[#xys][0] ~= xys[1][0] or xys[#xys][1] ~= xys[1][1] then
+                xys[#xys+1] = xys[1]
+            end
+            -- http://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
+            local cx, cy, sa = 0, 0, 0
+            for j = 1, #xys-1 do
+                local x0, y0 = unpack(xys[j])
+                local x1, y1 = unpack(xys[j+1])
+                local f = x0 * y1 - x1 * y0
+                cx = cx + (x0 + x1) * f
+                cy = cy + (y0 + y1) * f
+                sa = sa + f
+            end
+            cx = cx / (3 * sa)
+            cy = cy / (3 * sa)
+            cw[#cw+1] = {cx, cy, sa}
+        end
+        local cx, cy, sw = 0, 0, 0
+        for i = 1, #cw do
+            local x, y, w = unpack(cw[i])
+            cx = cx + x * w
+            cy = cy + y * w
+            sw = sw + w
+        end
+        cx = cx / sw
+        cy = cy / sw
+        lon1, lat1 = prj:inv(cx, cy)
+        if abs(lon1-lon0) <= epsilon and abs(lat1-lat0) <= epsilon then
+            break
+        end
+        lon0, lat0 = lon1, lat1
+    end
+    return lon1, lat1
+end
+
+-- == Projections ==
+
+-- Lambert Azimuthal Equal-Area Projection for the Spherical Earth.
+
+local AzimuthalEqualArea = {}
+AzimuthalEqualArea.__index = AzimuthalEqualArea
+
+function AzimuthalEqualArea:map(lon, lat)
+    lon, lat = rad(lon), rad(lat)
+    lon = lon - self.lon
+    local k, x, y
+    k = sqrt(2 / (1 + sin(self.lat) * sin(lat) + cos(self.lat) * cos(lat) * cos(lon)))
+    x = self.r * k * cos(lat) * sin(lon)
+    y = self.r * k * (cos(self.lat) * sin(lat) - sin(self.lat) * cos(lat) * cos(lon))
+    return x, y
+end
+
+function AzimuthalEqualArea:inv(x, y)
+    local p = sqrt(x*x + y*y)
+    local c = 2 * asin(p / (2 * self.r))
+    local lon, lat
+    -- FIXME: In the formulas below, should it be atan or atan2?
+    if self.lat == pi / 2 then
+        -- North Polar Aspect.
+        lon = self.lon + atan(x/(-y))
+    elseif self.lat == -pi / 2 then
+        -- South Polar Aspect.
+        lon = self.lon + atan(x/y)
+    else
+        -- Any other Oblique Aspect.
+        local den = p * cos(self.lat) * cos(c) - y * sin(self.lat) * sin(c)
+        lon = self.lon + atan(x * sin(c) / den)
+    end
+    lat = asin(cos(c) * sin(self.lat) + y * sin(c) * cos(self.lat) / p)
+    lon, lat = deg(lon), deg(lat)
+    return lon, lat
+end
+
+function AzimuthalEqualArea:model()
+    return {type="sphere", r=self.r}
+end
+
+local projs = {
+    AzimuthalEqualArea=AzimuthalEqualArea
+}
+
+-- Generic projection interface.
+
+local function Proj(name, origin, radius)
+    local proj = {}
+    proj.name = name
+    proj.lon, proj.lat = unpack(origin or {0, 0})
+    proj.lon, proj.lat = rad(proj.lon), rad(proj.lat)
+    proj.r = radius or 6378137
+    return setmetatable(proj, projs[name])
+end
+
+-- == Frames ==
+
+--[[
+ A frame stores information that specify the relation between the Earth's
+geometry (geodesy) and a particular map geometry (usually in 2D). Each map has a
+frame over which multiple layers of entities are drawn. There are two parameters
+that define a unique frame:
+* a projection;
+* a bounding box on the projected (plane) space.
+
+ The projection parameter must be fully specified, including the center of the
+projection, the orientation of the projection and the Earth model used (sphere
+or ellipsoid) along with its parameters.
+
+ The bounding box is specified as two points, determining minimal an maximal
+coordinates. The coordinate system for the bounding box is the projected one,
+but without scale, i.e. with meters as unit.
+]]
+
+local sep = ":"
+
+local Frame = {}
+Frame.__index = Frame
+
+function Frame:fit(x, y)
+    x = (x - self.bbox.x0) * self.s
+    y = self.h - (y - self.bbox.y0) * self.s
+    return x, y
+end
+
+function Frame:map(lon, lat)
+    local x, y = self.proj:map(lon, lat)
+    return self:fit(x, y)
+end
+
+function Frame:set_height(h)
+    local mw = self.bbox.x1 - self.bbox.x0
+    local mh = self.bbox.y1 - self.bbox.y0
+    self.h = h
+    self.s = h / mh
+    self.w = math.floor(mw * self.s + 0.5)
+end
+
+function Frame:add_margin(m)
+    local f = (self.h + m) / self.h
+    self.s = self.s / f
+    local mw = self.bbox.x1 - self.bbox.x0
+    local mh = self.bbox.y1 - self.bbox.y0
+    local cx = (self.bbox.x0 + self.bbox.x1) / 2
+    local cy = (self.bbox.y0 + self.bbox.y1) / 2
+    self.bbox.x0 = cx - mw/2 * f
+    self.bbox.x1 = cx + mw/2 * f
+    self.bbox.y0 = cy - mh/2 * f
+    self.bbox.y1 = cy + mh/2 * f
+end
+
+function Frame:fitted(polys)
+    self.touched = false
+    return function()
+        local points = polys()
+        if points then
+            return function()
+                local point = points()
+                if point then
+                    local x, y = unpack(point)
+                    x, y = self:fit(x, y)
+                    if x >= 0 and x < self.w and y >= 0 and y < self.h then
+                        self.touched = true
+                    end
+                    return {x, y}
+                end
+            end
+        end
+    end
+end
+
+function Frame:mapped(polys)
+    self.touched = false
+    return function()
+        local points = polys()
+        if points then
+            return function()
+                local point = points()
+                if point then
+                    local lon, lat = unpack(point)
+                    local x, y = self:map(lon, lat)
+                    if x >= 0 and x < self.w and y >= 0 and y < self.h then
+                        self.touched = true
+                    end
+                    return {x, y}
+                end
+            end
+        end
+    end
+end
+
+function Frame:save(fname)
+    local frm = io.open(fname, "w")
+    local model = self.proj:model()
+    frm:write("type", sep, model.type, "\n")
+    if model.type == "ellipsoid" then
+        frm:write("a", sep, model.a, "\n")
+        frm:write("b", sep, model.b, "\n")
+        frm:write("e", sep, model.e, "\n")
+        frm:write("f", sep, model.f, "\n")
+    elseif model.type == "sphere" then
+        frm:write("r", sep, model.r, "\n")
+    end
+    frm:write("proj", sep, self.proj.name, "\n")
+    frm:write("lon", sep, deg(self.proj.lon), "\n")
+    frm:write("lat", sep, deg(self.proj.lat), "\n")
+    frm:write("x0", sep, self.bbox.x0, "\n")
+    frm:write("y0", sep, self.bbox.y0, "\n")
+    frm:write("x1", sep, self.bbox.x1, "\n")
+    frm:write("y1", sep, self.bbox.y1, "\n")
+    frm:close()
+end
+
+local function new_frame(proj, bbox)
+    local self = setmetatable({}, Frame)
+    self.proj = proj
+    self.bbox = bbox
+    self.w, self.h, self.s = bbox.x1-bbox.x0, bbox.y1-bbox.y0, 1
+    return self
+end
+
+local function load_frame(fname)
+    local self = setmetatable({}, Frame)
+    local frm = io.open(fname, "r")
+    local function get(field)
+        local line = frm:read()
+        local got = line:sub(1, #field)
+        assert(got == field, "expected field "..field.." but got "..got)
+        return line:sub(#field+#sep+1)
+    end
+    local model = {}
+    model.type = get "type"
+    if model.type == "ellipsoid" then
+        model.a = tonumber(get "a")
+        model.b = tonumber(get "b")
+        model.e = tonumber(get "e")
+        model.f = tonumber(get "f")
+    elseif model.type == "sphere" then
+        model.r = tonumber(get "r")
+    end
+    local proj_name = get "proj"
+    local proj_lon = tonumber(get "lon")
+    local proj_lat = tonumber(get "lat")
+    proj = Proj(proj_name, {proj_lon, proj_lat}, model.r)
+    local bbox = {}
+    bbox.x0 = tonumber(get "x0")
+    bbox.y0 = tonumber(get "y0")
+    bbox.x1 = tonumber(get "x1")
+    bbox.y1 = tonumber(get "y1")
+    frm:close()
+    return new_frame(proj, bbox)
+end
+
+return {
+    distance=distance, bbox=bbox, centroid=centroid, Proj=Proj,
+    new_frame=new_frame, load_frame=load_frame
+}

diff --git a/lib/anim/poly.lua b/lib/anim/poly.lua
new file mode 100644
index 0000000..6dbaa55
--- /dev/null
+++ b/lib/anim/poly.lua
@@ -0,0 +1,197 @@
+local util = require "anim.util"
+
+local pi = math.pi
+local sqrt2 = math.sqrt(2)
+
+local function dashed(points, pattern)
+    local x0, y0, x1, y1
+    local cx, cy
+    local i, j
+    local d, h
+    local polylines = {}
+    local polyline
+    local draw = true
+    x0, y0 = unpack(points[1])
+    polyline = {{x0, y0}}
+    i = 2
+    x1, y1 = unpack(points[i])
+    j = 1
+    d = pattern[j]
+    while true do
+        h = util.hypot(x1-x0, y1-y0)
+        if d < h then
+            cx = x0 + (x1-x0)*d/h
+            cy = y0 + (y1-y0)*d/h
+            if draw then
+                table.insert(polyline, {cx, cy})
+                table.insert(polylines, polyline)
+            else
+                polyline = {{cx, cy}}
+            end
+            x0, y0 = cx, cy
+            draw = not draw
+            if j < #pattern then j = j + 1 else j = 1 end
+            d = pattern[j]
+        else
+            if draw then
+                table.insert(polyline, {x1, y1})
+            end
+            d = d - h
+            if i < #points then i = i + 1 else break end
+            x0, y0 = x1, y1
+            x1, y1 = unpack(points[i])
+        end
+    end
+    if draw then
+        table.insert(polyline, points[i])
+        table.insert(polylines, polyline)
+    end
+    return polylines
+end
+
+-- convert bezier curve {{ax, ay}, {bx, by}, {cx, cy}} to polyline
+local function bezier(curve)
+    local a, b, c, h
+    local ax, ay, bx, by, cx, cy
+    local dx, dy, ex, ey, fx, fy
+    local points = {}
+    local stack = {curve}
+    while #stack > 0 do
+        a, b, c = unpack(table.remove(stack))
+        ax, ay = unpack(a)
+        bx, by = unpack(b)
+        cx, cy = unpack(c)
+        h = math.abs((ax-cx)*(by-ay)-(ax-bx)*(cy-ay))/util.hypot(cx-ax, cy-ay)
+        if h > 1 then -- split curve
+            dx, dy = (ax+bx)/2, (ay+by)/2
+            fx, fy = (bx+cx)/2, (by+cy)/2
+            ex, ey = (dx+fx)/2, (dy+fy)/2
+            table.insert(stack, {{ex, ey}, {fx, fy}, {cx, cy}})
+            table.insert(stack, {{ax, ay}, {dx, dy}, {ex, ey}})
+        else -- add point to polyline
+            table.insert(points, {cx, cy})
+        end
+    end
+    return points
+end
+
+-- convert a sequence of linked Bézier curves to a polyline
+local function unfold(control_points)
+    local s, a, b
+    local px, py, qx, qy, rx, ry
+    s = control_points[1]
+    b = control_points[2]
+    local points = {s}
+    local sub
+    if b[3] then
+        table.insert(points, b)
+        s = b
+    end
+    for i = 3, #control_points do
+        a = b
+        b = control_points[i]
+        if a[3] then
+            if b[3] then
+                table.insert(points, b)
+                s = b
+            end
+        else
+            if b[3] then
+                px, py, _ = unpack(s)
+                qx, qy, _ = unpack(a)
+                rx, ry, _ = unpack(b)
+                sub = bezier({{px, py}, {qx, qy}, {rx, ry}})
+                for i, p in ipairs(sub) do
+                    table.insert(points, p)
+                end
+                s = b
+            else
+                px, py, _ = unpack(s)
+                qx, qy, _ = unpack(a)
+                rx, ry, _ = unpack(b)
+                rx, ry = (qx+rx)/2, (qy+ry)/2
+                sub = bezier({{px, py}, {qx, qy}, {rx, ry}})
+                for i, p in ipairs(sub) do
+                    table.insert(points, p)
+                end
+                s = {rx, ry, true}
+            end
+        end
+    end
+    return points
+end
+
+-- bezigon that approximates a circle
+local function bcircle(x, y, r)
+    local a = r * (sqrt2-1)
+    return {
+        {x  , y+r, true},
+        {x+a, y+r, false},
+        {x+r, y+a, false},
+        {x+r, y-a, false},
+        {x+a, y-r, false},
+        {x-a, y-r, false},
+        {x-r, y-a, false},
+        {x-r, y+a, false},
+        {x-a, y+r, false},
+        {x  , y+r, true}
+    }
+end
+
+-- regular polygon
+local function ngon(x, y, r, n, mina, maxa)
+    mina = mina or 0
+    maxa = maxa or 2 * pi
+    local a = 2 * pi / n -- angle between points
+    local pgon = {}
+    local px, py
+    local pa = mina
+    while pa < maxa do
+        px = x + r * math.cos(pa)
+        py = y + r * math.sin(pa)
+        table.insert(pgon, {px, py})
+        pa = pa + a
+    end
+    px = x + r * math.cos(maxa)
+    py = y + r * math.sin(maxa)
+    table.insert(pgon, {px, py})
+    return pgon
+end
+
+-- approximate an arc between angles mina and maxa
+local function parc(x, y, r, mina, maxa)
+    local h = 0.5 -- maximum radius-apothem allowed
+    local n = math.ceil(pi / math.acos(1 - h/r)) -- # of sides
+    return ngon(x, y, r, n, mina, maxa)
+end
+
+-- regular polygon that approximates a circle
+local function pcircle(x, y, r)
+    return parc(x, y, r, 0, 2 * pi)
+end
+
+local function arrow_head(x0, y0, x1, y1, w, h)
+    local dx, dy = x1-x0, y1-y0
+    local a = math.atan2(dy, dx)    -- line angle
+    local b = math.atan2(-dx, dy)   -- perpendicular angle
+    local mx, my = math.cos(a), math.sin(a)
+    local nx, ny = math.cos(b), math.sin(b)
+    local bx, by = x1 - mx * h, y1 - my * h     -- back of arrow
+    local lx, ly = bx - nx * w/2, by - ny * w/2 -- left point
+    local rx, ry = bx + nx * w/2, by + ny * w/2 -- right point
+    bx, by = bx + mx * h/2, by + my * h/2       -- off-curve point
+    local control_points = {
+      {lx, ly, true},
+      {x1, y1, true},
+      {rx, ry, true},
+      {bx, by, false},
+      {lx, ly, true},
+    }
+    return unfold(control_points)
+    --~ return control_points
+end
+
+return {
+  dashed=dashed, unfold=unfold, bcircle=bcircle, ngon=ngon,
+  parc=parc, pcircle=pcircle, arrow_head=arrow_head
+}

diff --git a/lib/anim/shp.lua b/lib/anim/shp.lua
new file mode 100644
index 0000000..3092305
--- /dev/null
+++ b/lib/anim/shp.lua
@@ -0,0 +1,375 @@
+local bit = require "bit"
+
+local util = require "anim.util"
+local bio = require "anim.bio"
+
+local SF = {}
+SF.__index = SF
+
+function SF:read_dbf()
+    local fp = io.open(self.path .. ".dbf", "rb")
+    local version = bit.band(bio.read_byte(fp), 0x07)
+    if version ~= 3 then
+        error("only DBF version 5 is supported")
+    end
+    self.year = 1900 + bio.read_byte(fp)
+    self.month = bio.read_byte(fp)
+    self.day = bio.read_byte(fp)
+    self.nrecs = bio.read_leu32(fp)
+    fp:seek("cur", 24)
+    local reclen = 0
+    local fields = {}
+    local byte = bio.read_byte(fp)
+    while byte ~= 0x0D do
+        local field_name = util.rtrim(string.char(byte) .. fp:read(10), "\000")
+        local field_type = fp:read(1)
+        fp:seek("cur", 4)
+        local field_length = bio.read_byte(fp)
+        reclen = reclen + field_length
+        local field_dec_count = bio.read_byte(fp)
+        local field = {
+          name=field_name,
+          type=field_type,
+          length=field_length,
+          dec_count=field_dec_count
+        }
+        table.insert(fields, field)
+        fp:seek("cur", 14)
+        byte = bio.read_byte(fp)
+    end
+    local tab = {}
+    for i = 1, self.nrecs do
+        if fp:read(1) == "*" then
+            fp:seek("cur", reclen)
+        else
+            local row = {}
+            for j = 1, #fields do
+                local field = fields[j]
+                local cell = util.rtrim(fp:read(field.length), " ")
+                if field.type == "F" or field.type == "N" then
+                    cell = tonumber(cell) or 0
+                end
+                table.insert(row, cell)
+            end
+            table.insert(tab, row)
+        end
+    end
+    self.fields = fields
+    self.tab = tab
+    io.close(fp)
+end
+
+function SF:search(field_name, value)
+    local col
+    for i = 1, #self.fields do
+        if self.fields[i].name == field_name then
+            col = i
+            break
+        end
+    end
+    if col ~= nil then
+        for i = 1, self.nrecs do
+            if self.tab[i][col] == value then
+                return i
+            end
+        end
+    end
+end
+
+function SF:print_summary(n)
+    n = n or 5
+    local sep = ":"
+    for i = 1, #self.fields do
+        local field = self.fields[i]
+        local row = {field.name}
+        for j = 1, n do
+            table.insert(row, self.tab[j][i])
+        end
+        print(table.concat(row, sep))
+    end
+    print("records".. sep .. #self.tab)
+    print("shape".. sep .. self.header.shape)
+end
+
+function SF:tab2csv(sep, fp)
+    -- TODO: refactor to use table.concat(str_list, sep)
+    sep = sep or ":"
+    fp = fp or io.output()
+    for i = 1, #self.fields do
+        if i > 1 then
+            fp:write(sep)
+        end
+        fp:write(self.fields[i].name)
+    end
+    fp:write("\n")
+    for i = 1, #self.tab do
+        for j = 1, #self.fields do
+            if j > 1 then
+                fp:write(sep)
+            end
+            local cell = self.tab[i][j]
+            fp:write(cell)
+        end
+        fp:write("\n")
+    end
+end
+
+local shape_name = {
+    [ 0] = "null",          [31] = "multipatch",
+    [ 1] = "point",   [ 3] = "polyline",   [ 5] = "polygon",   [ 8] = "multipoint",
+    [11] = "point-z", [13] = "polyline-z", [15] = "polygon-z", [18] = "multipoint-z",
+    [21] = "point-m", [23] = "polyline-m", [25] = "polygon-m", [28] = "multipoint-m",
+}
+
+local function read_bbox(fp)
+    local bbox = {}
+    bbox.xmin = bio.read_led64(fp)
+    bbox.ymin = bio.read_led64(fp)
+    bbox.xmax = bio.read_led64(fp)
+    bbox.ymax = bio.read_led64(fp)
+    return bbox
+end
+
+-- SHX and SHP share the same header structure!
+local function read_sf_header(fp)
+    local file_code = bio.read_bei32(fp)
+    if file_code ~= 9994 then
+        error("does not seem to be a shapefile (invalid FileCode)")
+    end
+    fp:seek("cur", 20) -- unused, always all-zero
+    local file_len = bio.read_bei32(fp)
+    local version = bio.read_lei32(fp)
+    if version ~= 1000 then
+        error("invalid shapefile version")
+    end
+    local header = {}
+    header.shape = shape_name[bio.read_lei32(fp)]
+    header.bbox = read_bbox(fp)
+    header.zmin = bio.read_led64(fp)
+    header.zmax = bio.read_led64(fp)
+    header.mmin = bio.read_led64(fp)
+    header.mmax = bio.read_led64(fp)
+    return header
+end
+
+function SF:read_shx()
+    local fp = io.open(self.path .. ".shx", "rb")
+    self.header = read_sf_header(fp)
+    self.index = {}
+    for i = 1, self.nrecs do
+        local offset = bio.read_bei32(fp)
+        local length = bio.read_bei32(fp)
+        table.insert(self.index, {offset=offset, length=length})
+    end
+    io.close(fp)
+end
+
+function SF:read_record_header(index)
+    local fp = self.fp
+    fp:seek("set", self.index[index].offset * 2)
+    local num = bio.read_bei32(fp)
+    local len = bio.read_bei32(fp)
+    local shape = shape_name[bio.read_lei32(fp)]
+    return num, len, shape
+end
+
+function SF:read_bbox(index)
+    local shape_name = self.header.shape
+    assert(
+        util.startswith(shape_name, "polyline") or
+        util.startswith(shape_name, "polygon") or
+        util.startswith(shape_name, "multipoint"),
+        ("%s shape doesn't have bbox"):format(shape_name)
+    )
+    local num, len, shape = self:read_record_header(index)
+    if shape == "null" then
+        return nil
+    end
+    return read_bbox(self.fp)
+end
+
+function SF:get_point(index, proj)
+    assert(util.startswith(self.header.shape, "point"), "type mismatch")
+    local fp = self.fp
+    local num, len, shape = self:read_record_header(index)
+    if shape == "null" then
+        return nil
+    end
+    local x = bio.read_led64(fp)
+    local y = bio.read_led64(fp)
+    if proj ~= nil then
+        x, y = proj:map(x, y)
+    end
+    return {x, y}
+end
+
+function SF:get_polys(index, proj)
+    assert(util.startswith(self.header.shape, "poly"), "type mismatch")
+    local fp = self.fp
+    local num, len, shape = self:read_record_header(index)
+    if shape == "null" then
+        return nil
+    end
+    local bbox = read_bbox(fp)
+    local nparts = bio.read_lei32(fp)
+    --~ print(nparts)
+    local npoints = bio.read_lei32(fp)
+    local total = bio.read_lei32(fp) -- first index is always zero?
+    local lengths = {}
+    for i = 2, nparts do
+        local index = bio.read_lei32(fp)
+        table.insert(lengths, index - total)
+        total = index
+    end
+    table.insert(lengths, npoints - total)
+    local i = 0
+    return bbox, lengths, function()
+        i = i + 1
+        if i <= #lengths then
+            local j = 0
+            local length = lengths[i]
+            return function()
+                j = j + 1
+                if j <= length then
+                    local x = bio.read_led64(fp)
+                    local y = bio.read_led64(fp)
+                    if proj ~= nil then
+                        x, y = proj:map(x, y)
+                    end
+                    return {x, y}
+                end
+            end
+        end
+    end
+end
+
+function SF:close()
+    io.close(self.fp)
+end
+
+function SF:save_cache(fname, k, proj, scale, filter)
+    local indices = {}
+    for i = 1, #self.tab do
+        local row = self.tab[i]
+        local rec = {}
+        for j = 1, #self.fields do
+            rec[self.fields[j].name] = row[j]
+        end
+        local key = filter(rec)
+        if key ~= nil then
+            table.insert(indices, {i, key:sub(1, 16)})
+        end
+    end
+    local cache = io.open(fname, "w")
+    bio.write_beu32(cache, util.round(1000 / scale))
+    bio.write_byte(cache, k)
+    for i = 1, #indices do
+        local index, key = unpack(indices[i])
+        cache:write(key, string.rep("\0", 20 - #key))
+    end
+    cache:write(string.rep("\0", 20)) -- end list of entries
+    local function geo2grid(p)
+        local x, y = unpack(p)
+        x, y = proj:map(x, y)
+        x, y = util.round(x * scale), util.round(y * scale)
+        return x, y
+    end
+    for i = 1, #indices do
+        local index, key = unpack(indices[i])
+        local offset = cache:seek()
+        cache:seek("set", i * 20 + 1)
+        bio.write_beu32(cache, offset)
+        cache:seek("set", offset)
+        local bb, lens, polys = self:get_polys(index)
+        bio.write_beu16(cache, #lens)
+        for poly in polys do
+            local px, py = geo2grid(poly())
+            bio.write_bei16(cache, px)
+            bio.write_bei16(cache, py)
+            local qx, qy = geo2grid(poly())
+            local rice = bio.rice_w(cache, k)
+            for point in poly do
+                local rx, ry = geo2grid(point)
+                if (px-rx)*(qy-py) - (px-qx)*(ry-py) ~= 0 then
+                    local dx, dy = qx-px, qy-py
+                    rice:put_signed(dx)
+                    rice:put_signed(dy)
+                    px, py = qx, qy
+                end
+                qx, qy = rx, ry
+            end
+            rice:put_signed(0)
+            rice:put_signed(0)
+            rice:flush()
+        end
+    end
+    cache:close()
+end
+
+local function open_shapefile(path)
+    local self = setmetatable({path=path}, SF)
+    self:read_dbf()
+    self:read_shx()
+    self.fp = io.open(path .. ".shp", "rb")
+    return self
+end
+
+local Cache = {}
+Cache.__index = Cache
+
+function Cache:keys()
+    local cache = self.fp
+    local offset = 5
+    return function()
+        cache:seek("set", offset)
+        offset = offset + 20
+        local ckey = cache:read(16)
+        if ckey:byte() ~= 0 then
+            return ckey:sub(1, ckey:find("\0")-1)
+        end
+    end
+end
+
+function Cache:get_polys(key)
+    local cache = self.fp
+    cache:seek("set", 5)
+    local ckey = cache:read(16)
+    local offset = -1
+    while ckey:byte() ~= 0 do
+        if ckey:sub(1, ckey:find("\0")-1) == key then
+            offset = bio.read_beu32(cache)
+            break
+        end
+        cache:seek("cur", 4)
+        ckey = cache:read(16)
+    end
+    assert(offset > 0, ("key '%s' not found in cache"):format(key))
+    cache:seek("set", offset)
+    local npolys = bio.read_beu16(cache)
+    return function()
+        if npolys > 0 then
+            npolys = npolys - 1
+            local ox, oy = bio.read_bei16(cache), bio.read_bei16(cache)
+            local rice = bio.rice_r(cache, self.k)
+            local x, y = 0, 0
+            return function()
+                if x ~= ox or y ~= oy then
+                    x, y = ox, oy
+                    local dx, dy = rice:get_signed(), rice:get_signed()
+                    ox, oy = ox+dx, oy+dy
+                    return {x * self.s, y * self.s}
+                end
+            end
+        end
+    end
+end
+
+local function load_cache(fname)
+    local self = setmetatable({}, Cache)
+    self.fp = io.open(fname, "r")
+    self.s = bio.read_beu32(self.fp) / 1000
+    self.k = bio.read_byte(self.fp)
+    return self
+end
+
+return {open_shapefile=open_shapefile, load_cache=load_cache}

diff --git a/lib/anim/surf.lua b/lib/anim/surf.lua
new file mode 100644
index 0000000..d49a5d2
--- /dev/null
+++ b/lib/anim/surf.lua
@@ -0,0 +1,266 @@
+local ffi = require "ffi"
+local bit = require "bit"
+
+local util = require "anim.util"
+
+local bnot = bit.bnot
+local bor, band = bit.bor, bit.band
+local lshift, rshift =  bit.lshift,  bit.rshift
+
+local Surf = {}
+
+function Surf:inside(x, y)
+    return x >= 0 and x < self.w and y >= 0 and y < self.h
+end
+
+function Surf:hline(x, y, w, v)
+    for i = x, x+w-1 do
+        self:pset(i, y, v)
+    end
+end
+
+function Surf:vline(x, y, h, v)
+    for i = y, y+h-1 do
+        self:pset(x, i, v)
+    end
+end
+
+function Surf:disk_int(cx, cy, r, v)
+    if r == 0 then
+        self:pset(cx, cy, v)
+        return
+    end
+    local x, y, d = r, 0, 1-r
+    while x >= y do
+        self:hline(cx-x, cy+y, 2*x, v)
+        self:hline(cx-y, cy+x, 2*y, v)
+        self:hline(cx-x, cy-y, 2*x, v)
+        self:hline(cx-y, cy-x, 2*y, v)
+        y = y + 1
+        if d <= 0 then
+            d = d + 2*y + 1
+        else
+            x = x - 1
+            d = d + 2*(y-x) + 1
+        end
+    end
+end
+
+function Surf:disk(cx, cy, r, v)
+    self:disk_int(util.round(cx), util.round(cy), r, v)
+end
+
+function Surf:line(x0, y0, x1, y1, v, r)
+    r = r or 0
+    local dx, dy = x1-x0, y1-y0
+    local n = math.max(math.abs(dx), math.abs(dy))
+    local sx, sy = dx/n, dy/n
+    local x, y = x0, y0
+    self:disk_int(math.floor(x), math.floor(y), r, v)
+    for i = 1, n do
+        x = x + sx
+        y = y + sy
+        self:disk_int(math.floor(x), math.floor(y), r, v)
+    end
+end
+
+function Surf:polyline(points, v, r)
+    points = util.func_iter(points)
+    local x0, y0, x1, y1
+    x0, y0 = unpack(points())
+    for point in points do
+        x1, y1 = unpack(point)
+        self:line(x0, y0, x1, y1, v, r)
+        x0, y0 = x1, y1
+    end
+end
+
+function Surf:polylines(polylines, v, r)
+    polylines = util.func_iter(polylines)
+    for points in polylines do
+        self:polyline(points, v, r)
+    end
+end
+
+local function cross_comp(a, b)
+    return a[1] < b[1]
+end
+
+function Surf:scan(points)
+    points = util.func_iter(points)
+    if not self.scans then
+        self.scans = {}
+        for i = 0, self.h-1 do
+            self.scans[i] = {}
+        end
+    end
+    local x0, y0, x1, y1
+    local ax, ay, bx, by -- same line as above, but enforce ay < by
+    local sign
+    x0, y0 = unpack(points())
+    for point in points do
+        x1, y1 = unpack(point)
+        if y1 ~= y0 then
+            if y0 < y1 then
+                ax, ay, bx, by = x0, y0, x1, y1
+                sign =  1
+            else
+                ax, ay, bx, by = x1, y1, x0, y0
+                sign = -1
+            end
+            local slope = (bx-ax) / (by-ay)
+            ay, by = util.round(ay), util.round(by)
+            while ay < by do
+                if ay >= 0 and ay < self.h then
+                    table.insert(self.scans[ay], {ax, sign})
+                end
+                ax = ax + slope
+                ay = ay + 1
+            end
+        end
+        x0, y0 = x1, y1
+    end
+end
+
+function Surf:fill_scans(v)
+    local scan, wind
+    local x, sign
+    local ax, bx
+    for i = 0, self.h-1 do
+        scan = self.scans[i]
+        table.sort(scan, cross_comp)
+        wind = 0
+        for j, cross in ipairs(scan) do
+            x, sign = unpack(cross)
+            if wind == 0 then
+                ax = math.floor(x)
+            end
+            wind = wind + sign
+            if wind == 0 then
+                bx = math.ceil(x)
+                self:hline(ax, i, bx-ax, v)
+            end
+        end
+    end
+    self.scans = nil
+end
+
+function Surf:polygon(points, v)
+    self:scan(points)
+    self:fill_scans(v)
+end
+
+function Surf:polygons(polygons, v)
+    polygons = util.func_iter(polygons)
+    for points in polygons do
+        self:scan(points)
+    end
+    self:fill_scans(v)
+end
+
+function Surf:blit(x, y, surf, sx, sy, w, h)
+    for j = 0, h-1 do
+        for i = 0, w-1 do
+            self:pset(x+i, y+j, surf:pget(sx+i, sy+j))
+        end
+    end
+end
+
+function Surf:compose(layers)
+    for y = 0, self.h-1 do
+        for x = 0, self.w-1 do
+            for _, layer in ipairs(layers) do
+                local p = layer.map[layer.surf:pget(x, y)+1]
+                if p >= 0 then
+                    self:pset(x, y, p)
+                end
+            end
+        end
+    end
+end
+
+local BitMap = {}
+
+function BitMap:fill(v)
+    ffi.fill(self.p, self.t * self.h, v)
+end
+
+function BitMap:_index_shift_mask(x, y)
+    local index = y * self.t + math.floor(x / 8)
+    local shift = (7-x) % 8
+    local mask = lshift(1, shift)
+    return index, shift, mask
+end
+
+function BitMap:pget(x, y)
+    if not self:inside(x, y) then return 0 end
+    local index, shift, mask = self:_index_shift_mask(x, y)
+    return rshift(band(self.p[index], mask), shift)
+end
+
+function BitMap:pset(x, y, v)
+    if not self:inside(x, y) then return end
+    local index, shift, mask = self:_index_shift_mask(x, y)
+    local byte = self.p[index]
+    if v > 0 then
+        byte = bor(byte, mask)
+    else
+        byte = band(byte, bnot(mask))
+    end
+    self.p[index] = byte
+end
+
+function BitMap:save_pbm(fname)
+    local pbm = io.open(fname, "wb")
+    pbm:write("P4\n", self.w, " ", self.h, "\n")
+    local row = self.p + 0
+    for y = 0, self.h-1 do
+        pbm:write(ffi.string(row, self.t))
+        row = row + self.t
+    end
+    pbm:close()
+end
+
+setmetatable(BitMap, {__index=Surf})
+
+local function new_bitmap(w, h)
+    local self = setmetatable({w=w, h=h}, {__index=BitMap})
+    self.t = math.ceil(w / 8) -- stride, i.e., bytes/row
+    self.p = ffi.new("unsigned char[?]", self.t * self.h)
+    return self
+end
+
+local ByteMap = {}
+
+function ByteMap:fill(v)
+    ffi.fill(self.p, self.w * self.h, v)
+end
+
+function ByteMap:pget(x, y)
+    if not self:inside(x, y) then return 0 end
+    return self.p[y*self.w+x]
+end
+
+function ByteMap:pset(x, y, v)
+    if not self:inside(x, y) then return end
+    self.p[y*self.w+x] = v
+end
+
+function ByteMap:save_ppm(fname, colors)
+    local ppm = io.open(fname, "wb")
+    ppm:write("P6\n", self.w, " ", self.h, "\n", 255, "\n")
+    for i = 0, self.w*self.h-1 do
+        ppm:write(string.char(unpack(colors[self.p[i]+1])))
+    end
+    ppm:close()
+end
+
+setmetatable(ByteMap, {__index=Surf})
+
+local function new_bytemap(w, h)
+    local self = setmetatable({w=w, h=h}, {__index=ByteMap})
+    self.p = ffi.new("unsigned char[?]", self.w * self.h)
+    return self
+end
+
+return {new_bitmap=new_bitmap, new_bytemap=new_bytemap}

diff --git a/lib/anim/ttf.lua b/lib/anim/ttf.lua
new file mode 100644
index 0000000..6aca3b5
--- /dev/null
+++ b/lib/anim/ttf.lua
@@ -0,0 +1,624 @@
+local bit = require "bit"
+
+local aff = require "anim.aff"
+local poly = require "anim.poly"
+
+local bnot = bit.bnot
+local bor, band = bit.bor, bit.band
+local lshift, rshift =  bit.lshift,  bit.rshift
+
+local function utf8to32(utf8str)
+    assert(type(utf8str) == "string")
+    local res, seq, val = {}, 0, nil
+    for i = 1, #utf8str do
+        local c = string.byte(utf8str, i)
+        if seq == 0 then
+            table.insert(res, val)
+            seq = c < 0x80 and 1 or c < 0xE0 and 2 or c < 0xF0 and 3 or
+                  c < 0xF8 and 4 or --c < 0xFC and 5 or c < 0xFE and 6 or
+                  error("invalid UTF-8 character sequence")
+            val = band(c, 2^(8-seq) - 1)
+        else
+            val = bor(lshift(val, 6), band(c, 0x3F))
+        end
+        seq = seq - 1
+    end
+    table.insert(res, val)
+    return res
+end
+
+-- Note: TrueType uses big endian for everything.
+
+local function uint(s)
+    local x = 0
+    for i = 1, #s do
+        x = x * 256 + s:byte(i)
+    end
+    return x
+end
+
+local function int(s)
+    local x = uint(s)
+    local p = 2^(#s*8)
+    if x >= p/2 then x = x - p end
+    return x
+end
+
+local function str(x)
+    local s = ""
+    while x > 0 do
+        s = string.char(x % 256) .. s
+        x = math.floor(x / 256)
+    end
+    return s
+end
+
+local Face = {}
+Face.__index = Face
+
+function Face:log(s)
+    if self.dbg then
+        io.stderr:write(s .. "\n")
+    end
+end
+
+function Face:str(n)    return self.fp:read(n) end
+function Face:uint8()   return uint(self.fp:read(1)) end
+function Face:uint16()  return uint(self.fp:read(2)) end
+function Face:uint32()  return uint(self.fp:read(4)) end
+function Face:uint64()  return uint(self.fp:read(8)) end
+function Face:int8()    return  int(self.fp:read(1)) end
+function Face:int16()   return  int(self.fp:read(2)) end
+function Face:int32()   return  int(self.fp:read(4)) end
+function Face:int64()   return  int(self.fp:read(8)) end
+
+function Face:goto(tag, offset)
+    self.fp:seek("set", self.dir[tag].offset + (offset or 0))
+end
+
+function Face:getpos()
+    return self.fp:seek()
+end
+
+function Face:setpos(pos)
+    self.fp:seek("set", pos)
+end
+
+function Face:offset()
+    local scaler_type = self:uint32()
+    assert(scaler_type == 0x74727565 or scaler_type == 0x00010000,
+        ("invalid scaler type for TrueType: 0x%08X"):format(scaler_type))
+    local num_tables = self:uint16()
+    --  The entries for search_range, entry_selector and range_shift are used to
+    -- facilitate quick binary searches of the table directory. Unless a font 
+    -- has a large number of tables, a sequential search will be fast enough.
+    local search_range = self:uint16()
+    local entry_selector = self:uint16()
+    local range_shift = self:uint16()
+    self.num_tables = num_tables
+end
+
+function Face:directory()
+    local dir = {}
+    for i = 1, self.num_tables do
+        local tag = self:str(4)
+        local checksum = self:uint32()
+        local offset = self:uint32()
+        local length = self:uint32()
+        -- TODO: verify checksums
+        dir[tag] = {checksum=checksum, offset=offset, length=length}
+    end
+    self.dir = dir
+end
+
+function Face:head()
+    self:goto("head")
+    local version = self:uint32()
+    assert(version == 0x00010000, ("invalid version: 0x%08X"):format(version))
+    local revision = self:uint32()
+    local checksum_adj = self:uint32()
+    local magic = self:uint32()
+    assert(magic == 0x5F0F3CF5, ("invalid magic: 0x%08X"):format(magic))
+    local flags = self:uint16()
+    self.units_per_em = self:uint16()
+    local created = self:int64()
+    local modified = self:int64()
+    local xmin = self:int16()
+    local ymin = self:int16()
+    local xmax = self:int16()
+    local ymax = self:int16()
+    local mac_style = self:uint16()
+    local lowest_rec_ppem = self:uint16()
+    local direction_hint = self:int16()
+    self.index_to_loc_fmt = self:int16()
+    local glyph_data_fmt = self:int16()
+end
+
+function Face:maxp()
+    self:goto("maxp")
+    local version = self:uint32()
+    if version == 0x00005000 then
+        self.num_glyphs = self:uint16()
+    elseif version == 0x00010000 then
+        self.num_glyphs = self:uint16()
+        local max_points = self:uint16()
+        local max_contours = self:uint16()
+        local max_composite_points = self:uint16()
+        local max_composite_contours = self:uint16()
+        local max_zones = self:uint16()
+        local max_twilight_points = self:uint16()
+        local max_storage = self:uint16()
+        local max_function_defs = self:uint16()
+        local max_instruction_defs = self:uint16()
+        local max_stack_elements = self:uint16()
+        local max_size_of_instructions = self:uint16()
+        local max_component_elements = self:uint16()
+        local max_component_depth = self:uint16()
+    else
+        error(("invalid maxp version: 0x%08X"):format(version))
+    end
+end
+
+function Face:cmap()
+    self:goto("cmap")
+    local version = self:uint16()
+    assert(version == 0, ("invalid cmap version: %d"):format(version))
+    local num_subtables = self:uint16()
+    local encoding, suboffset
+    local ok = false
+    for i = 1, num_subtables do
+        local platform_id = self:uint16()
+        encoding = self:uint16()
+        suboffset = self:uint32()
+        if platform_id == 0 then        -- platform == Unicode
+            ok = true
+            break
+        elseif platform_id == 3 then    -- platform == Microsoft
+            if encoding == 10 or encoding == 1 then
+                ok = true
+                break
+            end
+        end
+    end
+    if not ok then
+        error(("could not find Unicode cmap in %d subtables"):format(num_subtables))
+    end
+    self:goto("cmap", suboffset)
+    self:subcmap()
+end
+
+function Face:subcmap()
+    local format = self:uint16()
+    local segs = {}
+    local gia = {}
+    if format == 4 then
+        local length = self:uint16()
+        local language = self:uint16()
+        assert(language == 0, ("invalid subcmap language: %d"):format(language))
+        local seg_count = self:uint16() / 2
+        local search_range = self:uint16()
+        local entry_selector = self:uint16()
+        local range_shift = self:uint16()
+        for i = 1, seg_count do
+            segs[i] = {end_code=self:uint16()}
+        end
+        local last = segs[seg_count].end_code
+        assert(last == 0xFFFF, ("invalid subcmap last end code: %d"):format(last))
+        local pad = self:uint16()
+        assert(pad == 0, ("invalid subcmap reserved pad: %d"):format(pad))
+        for i = 1, seg_count do
+            segs[i].start_code = self:uint16()
+        end
+        for i = 1, seg_count do
+            segs[i].id_delta = self:uint16()
+        end
+        for i = 1, seg_count do
+            segs[i].id_range_offset = self:uint16()
+        end
+        local gia_len = (length - (16+8*seg_count)) / 2
+        for i = 1, gia_len do
+            gia[i] = self:uint16()
+        end
+    -- TODO: support other formats, specially 6 and 12
+    else
+        error(("unsupported subcmap format: %d"):format(format))
+    end
+    self.segs, self.gia = segs, gia
+end
+
+function Face:os2()
+    self:goto("OS/2")
+    local version = self:uint16()
+    assert(version >= 2, ("invalid OS/2 version: %u"):format(version))
+    self.avg_char_width = self:int16()
+    self.weight_class = self:uint16()
+    self.width_class = self:uint16()
+    self.fp:seek("cur", 78)
+    self.x_height = self:int16()
+    self.cap_height = self:int16()
+    self.default_char = self:uint16()
+    self.break_char = self:uint16()
+end
+
+function Face:hhea()
+    self:goto("hhea")
+    local versionH = self:uint16()
+    local versionL = self:uint16()
+    assert(versionH == 1 and versionL == 0,
+        ("invalid hhea version: %d.%d"):format(versionH, versionL))
+    self.ascent  = self:int16()
+    self.descent = self:int16()
+    local line_gap = self:int16()
+    local advance_width_max = self:uint16()
+    local min_left_side_bearing  = self:int16()
+    local min_right_side_bearing = self:int16()
+    local x_max_extent = self:int16()
+    local caret_slope_rise = self:int16()
+    local caret_slope_run  = self:int16()
+    local caret_offset     = self:int16()
+    for i = 1, 4 do
+        local reserved = self:uint16()
+        assert(reserved == 0, "nonzero reserved field in hhea")
+    end
+    local metric_data_format = self:int16()
+    assert(metric_data_format == 0,
+        ("invalid metric data format: %d"):format(metric_data_format))
+    self.nlong_hor_metrics = self:uint16()
+end
+
+function Face:hmetrics(id)
+    local n = self.nlong_hor_metrics -- for readability of expressions below
+    local advance, bearing
+    if id < n then
+        self:goto("hmtx", 4*id)
+        advance = self:uint16()
+        bearing = self:int16()
+    else
+        self:goto("hmtx", 4*(n-1))
+        advance = self:uint16()
+        self:goto("hmtx", 4*n+2*(id-n))
+        bearing = self:int16()
+    end
+    return advance, bearing
+end
+
+function Face:kern()
+    self:goto("kern")
+    local version = self:uint16()
+    assert(version == 0, ("invalid kern table version: %d"):format(version))
+    local ntables = self:uint16()
+    for i = 1, ntables do
+        local version = self:uint16()
+        local length = self:uint16()
+        local format = self:uint8() -- usually 0
+        local coverage = self:uint8()
+        local horizontal   = band(coverage, 2^0) > 0 -- usually true
+        local minimum      = band(coverage, 2^1) > 0 -- usually false (kerning)
+        local cross_stream = band(coverage, 2^2) > 0 -- usually false (regular)
+        local override     = band(coverage, 2^3) > 0 -- usually false (accumulate)
+        assert(band(coverage, 0xF0) == 0, "invalid coverage bits set")
+        if format == 0 then
+            self.num_kernings = self:uint16()
+            local search_range = self:uint16()
+            local entry_selector = self:uint16()
+            local range_shift = self:uint16()
+            local kerning = {}
+            for j = 1, self.num_kernings do
+                -- glyph indices (left * 2^16 + right)
+                local key  = self:uint32()
+                -- kerning value
+                local value = self:int16()
+                kerning[key] = value
+            end
+            self.kerning = kerning
+        else
+            self:log(("unsupported kerning table format: %d"):format(format))
+        end
+    end
+end
+
+function Face:get_kerning(left_id, right_id)
+    return self.kerning[left_id * 2^16 + right_id] or 0
+end
+
+-- Convert a character code to its glyph id.
+function Face:char_index(code)
+    local i = 1
+    while code > self.segs[i].end_code do i = i + 1 end
+    if self.segs[i].start_code > code then return 0 end
+    local iro = self.segs[i].id_range_offset
+    if iro == 0 then
+        return (code + self.segs[i].id_delta) % 0x10000
+    else
+        local idx = iro + 2 * (code - self.segs[i].start_code)
+        idx = idx - (#self.segs - i + 1) * 2
+        local id = self.gia[idx/2+1]
+        if id > 0 then
+            id = (id + self.segs[i].id_delta) % 0x10000
+        end
+        return id
+    end
+end
+
+function Face:glyph(id)
+    local suboffset
+    if self.index_to_loc_fmt == 0 then      -- short offsets
+        self:goto("loca", 2*id)
+        suboffset = self:uint16() * 2
+    else                                    -- long offsets
+        self:goto("loca", 4*id)
+        suboffset = self:uint32()
+    end
+    self:goto("glyf", suboffset)
+    local num_contours = self:int16()
+    local xmin = self:int16()
+    local ymin = self:int16()
+    local xmax = self:int16()
+    local ymax = self:int16()
+    local points, end_points = {}, {}
+    if num_contours > 0 then        -- simple glyph
+        for i = 1, num_contours do
+            end_points[i] = self:uint16() + 1
+        end
+        local num_points = end_points[#end_points]
+        local instruction_length = self:uint16()
+        local instructions = self:str(instruction_length)
+        local i = 0
+        while i < num_points do
+            i = i + 1
+            local flags = self:uint8()
+            assert(flags < 64, "point flag with higher bits set")
+            local point = {
+                on_curve    = band(flags, 2^0) > 0,
+                x_short     = band(flags, 2^1) > 0,
+                y_short     = band(flags, 2^2) > 0,
+                repeated    = band(flags, 2^3) > 0,
+                x_sign_same = band(flags, 2^4) > 0,
+                y_sign_same = band(flags, 2^5) > 0
+            }
+            points[i] = point
+            if point.repeated then
+                local repeats = self:uint8()
+                for j = 1, repeats do
+                    i = i + 1
+                    points[i] = {
+                        on_curve    = point.on_curve,
+                        x_short     = point.x_short,
+                        y_short     = point.y_short,
+                        x_sign_same = point.x_sign_same,
+                        y_sign_same = point.y_sign_same
+                    }
+                end
+            end
+        end
+        local last_x, last_y = 0, 0
+        for i = 1, #points do
+            if points[i].x_short then
+                local x = self:uint8()
+                if not points[i].x_sign_same then x = -x end
+                points[i].x = last_x + x
+            else
+                if not points[i].x_sign_same then
+                    points[i].x = last_x + self:int16()
+                else
+                    points[i].x = last_x
+                end
+            end
+            last_x = points[i].x
+        end
+        for i = 1, #points do
+            if points[i].y_short then
+                local y = self:uint8()
+                if not points[i].y_sign_same then y = -y end
+                points[i].y = last_y + y
+            else
+                if not points[i].y_sign_same then
+                    points[i].y = last_y + self:int16()
+                else
+                    points[i].y = last_y
+                end
+            end
+            last_y = points[i].y
+        end
+    elseif num_contours < 0 then    -- compound glyph
+        local more = true
+        while more do
+            local flags = self:uint16()
+            local args_are_words    = band(flags, 2^0x0) > 0
+            local args_are_xy       = band(flags, 2^0x1) > 0
+            local round_xy_to_grid  = band(flags, 2^0x2) > 0
+            local regular_scale     = band(flags, 2^0x3) > 0
+            local obsolete          = band(flags, 2^0x4) > 0
+            more                    = band(flags, 2^0x5) > 0
+            local irregular_scale   = band(flags, 2^0x6) > 0
+            local two_by_two        = band(flags, 2^0x7) > 0
+            local instructions      = band(flags, 2^0x8) > 0
+            local use_my_metrics    = band(flags, 2^0x9) > 0
+            local overlap           = band(flags, 2^0xA) > 0
+            local scaled_offset     = band(flags, 2^0xB) > 0
+            local unscaled_offset   = band(flags, 2^0xC) > 0
+            if obsolete then
+                self:log("warning: glyph component using obsolete flag")
+            end
+            local gid = self:uint16()
+            local pos = self:getpos()
+            local sub_points, sub_end_points = self:glyph(gid)
+            self:setpos(pos)
+            local e, f
+            if args_are_xy then
+                if args_are_words then
+                    e = self:int16()
+                    f = self:int16()
+                else
+                    e = self:int8()
+                    f = self:int8()
+                end
+                if round_xy_to_grid then
+                    self:log("warning: ignoring request to round component offset")
+                end
+            else
+                local i, j
+                if args_are_words then
+                    i = self:uint16()
+                    j = self:uint16()
+                else
+                    i = self:uint8()
+                    j = self:uint8()
+                end
+                e = points[i].x - sub_points[j].x
+                f = points[i].y - sub_points[j].y
+            end
+            local a, b, c, d = 1, 0, 0, 1
+            if regular_scale then
+                self:log("regular scale")
+                a = self:int16() / 0x4000
+                d = a
+            elseif irregular_scale then
+                self:log("irregular scale")
+                a = self:int16() / 0x4000
+                d = self:int16() / 0x4000
+            elseif two_by_two then
+                self:log("2x2 transformation")
+                a = self:int16() / 0x4000
+                b = self:int16() / 0x4000
+                c = self:int16() / 0x4000
+                d = self:int16() / 0x4000
+            end
+            local m = math.max(math.abs(a), math.abs(b))
+            local n = math.max(math.abs(c), math.abs(d))
+            if math.abs(math.abs(a)-math.abs(c)) <= 0x21/0x10000 then
+                m = m * 2
+            end
+            if math.abs(math.abs(c)-math.abs(d)) <= 0x21/0x10000 then
+                n = n * 2
+            end
+            for i, p in ipairs(sub_points) do
+                points[#points+1] = {
+                    x=m*(a*p.x/m + c*p.y/m + e),
+                    y=n*(b*p.x/n + d*p.y/n + f),
+                    on_curve=p.on_curve
+                }
+            end
+            local offset = end_points[#end_points] or 0
+            for i, e in ipairs(sub_end_points) do
+                end_points[#end_points+1] = offset + e
+            end
+        end
+        -- TODO: read instructions for composite character
+    end
+    return points, end_points
+end
+
+function Face:pack_outline(points, end_points)
+    local outline = {}
+    local h = self.cap_height
+    local p, q
+    local j = 1
+    for i = 1, #end_points do
+        local contour = {}
+        while j <= end_points[i] do
+            p = points[j]
+            q = {p.x, h-p.y, p.on_curve}
+            table.insert(contour, q)
+            j = j + 1
+        end
+        table.insert(outline, contour)
+    end
+    return outline
+end
+
+function Face:string(s, pt, x, y, anchor, a)
+    anchor = anchor or "tl"
+    a = a or 0
+    local codes = utf8to32(s)
+    local cur_x = 0
+    local contours = {}
+    local outline
+    local li, ri
+    local advance, bearing
+    for i, code in ipairs(codes) do
+        ri = self:char_index(code)
+        if i > 1 and self.num_kernings > 0 then
+            cur_x = cur_x + self:get_kerning(li, ri)
+        end
+        if code ~= 32 then
+            outline = self:pack_outline(self:glyph(ri))
+            for j, contour in ipairs(outline) do
+                for k, point in ipairs(contour) do
+                    point[1] = point[1] + cur_x
+                end
+                table.insert(contours, contour)
+            end
+        end
+        advance, bearing = self:hmetrics(ri)
+        cur_x = cur_x + advance
+        li = ri
+    end
+    local ax, ay -- anchor position
+    local av, ah = anchor:sub(1, 1), anchor:sub(2, 2)
+    if av == "t" then
+        ay = 0
+    elseif av == "m" then
+        ay = self.cap_height / 2
+    elseif av == "b" then
+        ay = self.cap_height
+    end
+    if ah == "l" then
+        ax = 0
+    elseif ah == "c" then
+        ax = cur_x / 2
+    elseif ah == "r" then
+        ax = cur_x
+    end
+    if ax == nil or ay == nil then
+        error("invalid anchor: "..anchor)
+    end
+    local scl = pt * self.resolution / (72 * self.units_per_em)
+    local t = aff.new_affine()
+    t:add_translate(-ax, -ay)
+    t:add_scale(scl)
+    t:add_rotate(a)
+    t:add_translate(x, y)
+    local polygons = {}
+    local polygon
+    for i, contour in ipairs(contours) do
+        if #contour > 2 then
+            t:apply(contour)
+            polygon = poly.unfold(contour)
+            table.insert(polygon, polygon[1]) -- close
+            table.insert(polygons, polygon)
+        end
+    end
+    return polygons
+end
+
+local function load_face(f, dbg)
+    if type(f) == "string" then f = io.open(f, "rb") end
+    local self = setmetatable({fp=f}, Face)
+    self.dbg = dbg or false
+    self:offset()
+    self:directory()
+    self:head()
+    self:maxp()
+    self:cmap()
+    if self.dir["hhea"] then
+        self:hhea()
+    else
+        self:log("no horizontal metrics (hhea+hmtx)")
+    end
+    if self.dir["kern"] then
+        self:kern()
+    else
+        self.num_kernings = 0
+        self:log("no kerning table (kern)")
+    end
+    if self.dir["OS/2"] then
+        self:os2()
+    else
+        self:log("no x-height and Cap-Height (OS/2)")
+    end
+    self.resolution = 300 -- dpi
+    return self
+end
+
+return {load_face=load_face}

diff --git a/lib/anim/util.lua b/lib/anim/util.lua
new file mode 100644
index 0000000..63f141f
--- /dev/null
+++ b/lib/anim/util.lua
@@ -0,0 +1,45 @@
+local ffi = require "ffi"
+
+ffi.cdef[[
+double hypot(double x, double y);
+double copysign(double x, double y);
+]]
+local hypot = ffi.C.hypot
+local copysign = ffi.C.copysign
+
+local function round(x)
+    local i, f = math.modf(x + copysign(0.5, x))
+    return i
+end
+
+local function rtrim(s, c)
+    local i = #s
+    while s:sub(i, i) == c do
+        i = i - 1
+    end
+    return s:sub(1, i)
+end
+
+local function startswith(s1, s2)
+    return s1:sub(1, #s2) == s2
+end
+
+local function func_iter(seq)
+    if type(seq) == "table" then
+        local i = 0
+        return function()
+            i = i + 1
+            if i <= #seq then
+                return seq[i]
+            end
+        end
+    else
+        return seq
+    end
+end
+
+return {
+    hypot=hypot, copysign=copysign, round=round,
+    rtrim=rtrim, startswith=startswith,
+    func_iter=func_iter
+}

diff --git a/lzw.lua b/lzw.lua
deleted file mode 100644
index 377d23a..0000000
--- a/lzw.lua
+++ /dev/null
@@ -1,215 +0,0 @@
-local ffi = require "ffi"
-local bit = require "bit"
-
-local bnot = bit.bnot
-local bor, band = bit.bor, bit.band
-local lshift, rshift =  bit.lshift,  bit.rshift
-
--- == Encoder ==
-
-local BUFout = {}
-BUFout.__index = BUFout
-
-local function new_buffer_out(f)
-    local self = setmetatable({f=f, offset=0, partial=0}, BUFout)
-    self.buf = ffi.new("char[255]")
-    return self
-end
-
-function BUFout:put_key(key, size)
-    local offset, partial = self.offset, self.partial
-    local f, buf = self.f, self.buf
-    local byte_offset, bit_offset = math.floor(offset / 8), offset % 8
-    partial = bor(partial, lshift(key, bit_offset))
-    local bits_to_write = bit_offset + size
-    while bits_to_write >= 8 do
-        buf[byte_offset] = band(partial, 0xFF)
-        byte_offset = byte_offset + 1
-        if byte_offset == 0xFF then -- flush
-            f:write(string.char(0xFF))
-            f:write(ffi.string(buf, 0xFF))
-            byte_offset = 0
-        end
-        partial = rshift(partial, 8)
-        bits_to_write = bits_to_write - 8
-    end
-    self.offset = (offset + size) % (0xFF * 8)
-    self.partial = partial
-end
-
-function BUFout:end_key()
-    local offset, partial = self.offset, self.partial
-    local f, buf = self.f, self.buf
-    local byte_offset = math.floor(offset / 8)
-    if offset % 8 ~= 0 then
-        buf[byte_offset] = band(partial, 0xFF)
-        byte_offset = byte_offset + 1
-    end
-    if byte_offset > 0 then
-        f:write(string.char(byte_offset))
-        f:write(ffi.string(buf, byte_offset))
-    end
-    f:write(string.char(0))
-    self.offset, self.partial = 0, 0
-end
-
-
-local function new_trie(degree)
-    local children = {}
-    for key = 0, degree-1 do
-        children[key] = {key=key, children={}}
-    end
-    return {children=children}
-end
-
-local function encode(f, d, s, x, y, w, h)
-    local buf = new_buffer_out(f)
-    local code_size = math.max(d, 2)
-    f:write(string.char(code_size))
-    local degree = 2 ^ code_size
-    local root = new_trie(degree)
-    local clear, stop = degree, degree + 1
-    local nkeys = degree + 2 -- skip clear code and stop code
-    local node = root
-    local key_size = code_size + 1
-    buf:put_key(clear, key_size)
-    for j = y, y+h-1 do
-        for i = x, x+w-1 do
-            local index = band(s:pget(i, j), degree-1)
-            local child = node.children[index]
-            if child ~= nil then
-                node = child
-            else
-                buf:put_key(node.key, key_size)
-                if nkeys < 0x1000 then
-                    if nkeys == 2 ^ key_size then
-                        key_size = key_size + 1
-                    end
-                    node.children[index] = {key=nkeys, children={}}
-                    nkeys = nkeys + 1
-                else
-                    buf:put_key(clear, key_size)
-                    root = new_trie(degree)
-                    node = root
-                    nkeys = degree + 2
-                    key_size = code_size + 1
-                end
-                node = root.children[index]
-            end
-        end
-    end
-    buf:put_key(node.key, key_size)
-    buf:put_key(stop, key_size)
-    buf:end_key()
-end
-
--- == Decoder ==
-
-local BUFin = {}
-BUFin.__index = BUFin
-
-local function new_buffer_in(f)
-    local self = setmetatable({}, BUFin)
-    self.f = f      -- already opened file, read mode
-    self.s = 0      -- number of bytes available in block
-    self.b = 0      -- value of last byte read
-    self.n = 0      -- number of bits available in self.b
-    return self
-end
-
-function BUFin:get_key(size)
-    local key = 0
-    for i = 1, size do
-        if self.s == 0 then
-            self.s = self.f:read(1):byte(1)
-            assert(self.s > 0, "unexpected end-of-block")
-        end
-        if self.n == 0 then
-            self.b = self.f:read(1):byte(1)
-            self.n = 8
-            self.s = self.s - 1
-        end
-        key = bor(key, lshift(band(rshift(self.b, 8-self.n), 1), i-1))
-        self.n = self.n - 1
-    end
-    return key
-end
-
-local CodeTable = {}
-CodeTable.__index = CodeTable
-
-local function new_code_table(key_size)
-    local self = setmetatable({}, CodeTable)
-    self.len = lshift(1, key_size)
-    self.tab = {}
-    for key = 0, self.len+1 do
-        self.tab[key] = {length=1, prefix=0xFFF, suffix=key}
-    end
-    self.len = self.len + 2 -- clear & stop
-    return self
-end
-
-function CodeTable:add_entry(length, prefix, suffix)
-    self.tab[self.len] = {length=length, prefix=prefix, suffix=suffix}
-    self.len = self.len + 1
-    if band(self.len, self.len-1) == 0 then
-        return 1
-    end
-    return 0
-end
-
-local function decode(f, d, s, fx, fy, w, h)
-    local key_size = f:read(1):byte(1)
-    assert(key_size == math.max(d, 2), "invalid code size")
-    local buf = new_buffer_in(f)
-    local clear = lshift(key_size, 1)
-    local stop = clear + 1
-    key_size = key_size + 1
-    local init_key_size = key_size
-    local key = buf:get_key(key_size)
-    assert(key == clear, "expected clear code, got "..key)
-    local code_table, table_is_full, entry, str_len, ret
-    local frm_off = 0 -- pixels read
-    local frm_size = w * h
-    while frm_off < frm_size do
-        if key == clear then
-            key_size = init_key_size
-            code_table = new_code_table(key_size-1)
-            table_is_full = false
-        elseif not table_is_full then
-            ret = code_table:add_entry(str_len+1, key, entry.suffix)
-            if code_table.len == 0x1000 then
-                ret = 0
-                table_is_full = true
-            end
-        end
-        key = buf:get_key(key_size)
-        if key ~= clear then
-            if key == stop or key == 0x1000 then break end
-            if ret == 1 then key_size = key_size + 1 end
-            entry = code_table.tab[key]
-            str_len = entry.length
-            for i = 1, str_len do
-                local p = frm_off + entry.length - 1
-                local x = p % w
-                local y = math.floor(p / w)
-                s:pset(fx+x, fy+y, entry.suffix)
-                if entry.prefix == 0xFFF then
-                    break
-                else
-                    entry = code_table.tab[entry.prefix]
-                end
-            end
-            frm_off = frm_off + str_len
-            if key < code_table.len-1 and not table_is_full then
-                code_table.tab[code_table.len-1].suffix = entry.suffix
-            end
-        end
-    end
-    while buf.s > 0 do
-        f:seek("cur", buf.s)
-        buf.s = f:read(1):byte(1)
-    end
-end
-
-return {encode=encode, decode=decode}

diff --git a/map.lua b/map.lua
deleted file mode 100644
index 9976c4b..0000000
--- a/map.lua
+++ /dev/null
@@ -1,310 +0,0 @@
-local util = require "util"
-local bio = require "bio"
-
-local abs = math.abs
-local rad, deg = math.rad, math.deg
-local cos, sin, tan = math.cos, math.sin, math.tan
-local acos, asin, atan, atan2 = math.acos, math.asin, math.atan, math.atan2
-local sqrt = math.sqrt
-local pi, huge = math.pi, math.huge
-
--- == Utilities ==
-
--- region is a list of polygons in geographic coordinates.
-
-local function distance(lon1, lat1, lon2, lat2, r)
-    r = r or 6378137
-    local dlat = rad(lat2 - lat1)
-    local dlon = rad(lon2 - lon1)
-    lat1, lat2 = rad(lat1), rad(lat2)
-    local a1, a2, a, c
-    a1 = sin(dlat/2) * sin(dlat/2)
-    a2 = sin(dlon/2) * sin(dlon/2) * cos(lat1) * cos(lat2)
-    a = a1 + a2
-    c = 2 * atan2(sqrt(a), sqrt(1-a))
-    return r * c
-end
-
-local function bbox(polys)
-    local x0, y0, x1, y1 = huge, huge, -huge, -huge
-    for poly in util.func_iter(polys) do
-        for point in util.func_iter(poly) do
-            local x, y = unpack(point)
-            x0 = x < x0 and x or x0
-            y0 = y < y0 and y or y0
-            x1 = x > x1 and x or x1
-            y1 = y > y1 and y or y1
-        end
-    end
-    return {x0=x0, y0=y0, x1=x1, y1=y1}
-end
-
-local function centroid(region)
-    local epsilon = 1e-10
-    local bb = bbox(region)
-    local x0, y0, x1, y1 = bb.x0, bb.y0, bb.x1, bb.y1
-    local lon0 = (x0 + x1) / 2
-    local lat0 = (y0 + y1) / 2
-    local lon1, lat1
-    while true do
-        local prj = Proj("AzimuthalEqualArea", {lon0, lat0})
-        local cw = {}
-        for i, polygon in ipairs(region) do
-            local xys = {}
-            for j, point in ipairs(polygon) do
-                xys[j] = {prj:map(unpack(point))}
-            end
-            if xys[#xys][0] ~= xys[1][0] or xys[#xys][1] ~= xys[1][1] then
-                xys[#xys+1] = xys[1]
-            end
-            -- http://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
-            local cx, cy, sa = 0, 0, 0
-            for j = 1, #xys-1 do
-                local x0, y0 = unpack(xys[j])
-                local x1, y1 = unpack(xys[j+1])
-                local f = x0 * y1 - x1 * y0
-                cx = cx + (x0 + x1) * f
-                cy = cy + (y0 + y1) * f
-                sa = sa + f
-            end
-            cx = cx / (3 * sa)
-            cy = cy / (3 * sa)
-            cw[#cw+1] = {cx, cy, sa}
-        end
-        local cx, cy, sw = 0, 0, 0
-        for i = 1, #cw do
-            local x, y, w = unpack(cw[i])
-            cx = cx + x * w
-            cy = cy + y * w
-            sw = sw + w
-        end
-        cx = cx / sw
-        cy = cy / sw
-        lon1, lat1 = prj:inv(cx, cy)
-        if abs(lon1-lon0) <= epsilon and abs(lat1-lat0) <= epsilon then
-            break
-        end
-        lon0, lat0 = lon1, lat1
-    end
-    return lon1, lat1
-end
-
--- == Projections ==
-
--- Lambert Azimuthal Equal-Area Projection for the Spherical Earth.
-
-local AzimuthalEqualArea = {}
-AzimuthalEqualArea.__index = AzimuthalEqualArea
-
-function AzimuthalEqualArea:map(lon, lat)
-    lon, lat = rad(lon), rad(lat)
-    lon = lon - self.lon
-    local k, x, y
-    k = sqrt(2 / (1 + sin(self.lat) * sin(lat) + cos(self.lat) * cos(lat) * cos(lon)))
-    x = self.r * k * cos(lat) * sin(lon)
-    y = self.r * k * (cos(self.lat) * sin(lat) - sin(self.lat) * cos(lat) * cos(lon))
-    return x, y
-end
-
-function AzimuthalEqualArea:inv(x, y)
-    local p = sqrt(x*x + y*y)
-    local c = 2 * asin(p / (2 * self.r))
-    local lon, lat
-    -- FIXME: In the formulas below, should it be atan or atan2?
-    if self.lat == pi / 2 then
-        -- North Polar Aspect.
-        lon = self.lon + atan(x/(-y))
-    elseif self.lat == -pi / 2 then
-        -- South Polar Aspect.
-        lon = self.lon + atan(x/y)
-    else
-        -- Any other Oblique Aspect.
-        local den = p * cos(self.lat) * cos(c) - y * sin(self.lat) * sin(c)
-        lon = self.lon + atan(x * sin(c) / den)
-    end
-    lat = asin(cos(c) * sin(self.lat) + y * sin(c) * cos(self.lat) / p)
-    lon, lat = deg(lon), deg(lat)
-    return lon, lat
-end
-
-function AzimuthalEqualArea:model()
-    return {type="sphere", r=self.r}
-end
-
-local projs = {
-    AzimuthalEqualArea=AzimuthalEqualArea
-}
-
--- Generic projection interface.
-
-local function Proj(name, origin, radius)
-    local proj = {}
-    proj.name = name
-    proj.lon, proj.lat = unpack(origin or {0, 0})
-    proj.lon, proj.lat = rad(proj.lon), rad(proj.lat)
-    proj.r = radius or 6378137
-    return setmetatable(proj, projs[name])
-end
-
--- == Frames ==
-
---[[
- A frame stores information that specify the relation between the Earth's
-geometry (geodesy) and a particular map geometry (usually in 2D). Each map has a
-frame over which multiple layers of entities are drawn. There are two parameters
-that define a unique frame:
-* a projection;
-* a bounding box on the projected (plane) space.
-
- The projection parameter must be fully specified, including the center of the
-projection, the orientation of the projection and the Earth model used (sphere
-or ellipsoid) along with its parameters.
-
- The bounding box is specified as two points, determining minimal an maximal
-coordinates. The coordinate system for the bounding box is the projected one,
-but without scale, i.e. with meters as unit.
-]]
-
-local sep = ":"
-
-local Frame = {}
-Frame.__index = Frame
-
-function Frame:fit(x, y)
-    x = (x - self.bbox.x0) * self.s
-    y = self.h - (y - self.bbox.y0) * self.s
-    return x, y
-end
-
-function Frame:map(lon, lat)
-    local x, y = self.proj:map(lon, lat)
-    return self:fit(x, y)
-end
-
-function Frame:set_height(h)
-    local mw = self.bbox.x1 - self.bbox.x0
-    local mh = self.bbox.y1 - self.bbox.y0
-    self.h = h
-    self.s = h / mh
-    self.w = math.floor(mw * self.s + 0.5)
-end
-
-function Frame:add_margin(m)
-    local f = (self.h + m) / self.h
-    self.s = self.s / f
-    local mw = self.bbox.x1 - self.bbox.x0
-    local mh = self.bbox.y1 - self.bbox.y0
-    local cx = (self.bbox.x0 + self.bbox.x1) / 2
-    local cy = (self.bbox.y0 + self.bbox.y1) / 2
-    self.bbox.x0 = cx - mw/2 * f
-    self.bbox.x1 = cx + mw/2 * f
-    self.bbox.y0 = cy - mh/2 * f
-    self.bbox.y1 = cy + mh/2 * f
-end
-
-function Frame:fitted(polys)
-    self.touched = false
-    return function()
-        local points = polys()
-        if points then
-            return function()
-                local point = points()
-                if point then
-                    local x, y = unpack(point)
-                    x, y = self:fit(x, y)
-                    if x >= 0 and x < self.w and y >= 0 and y < self.h then
-                        self.touched = true
-                    end
-                    return {x, y}
-                end
-            end
-        end
-    end
-end
-
-function Frame:mapped(polys)
-    self.touched = false
-    return function()
-        local points = polys()
-        if points then
-            return function()
-                local point = points()
-                if point then
-                    local lon, lat = unpack(point)
-                    local x, y = self:map(lon, lat)
-                    if x >= 0 and x < self.w and y >= 0 and y < self.h then
-                        self.touched = true
-                    end
-                    return {x, y}
-                end
-            end
-        end
-    end
-end
-
-function Frame:save(fname)
-    local frm = io.open(fname, "w")
-    local model = self.proj:model()
-    frm:write("type", sep, model.type, "\n")
-    if model.type == "ellipsoid" then
-        frm:write("a", sep, model.a, "\n")
-        frm:write("b", sep, model.b, "\n")
-        frm:write("e", sep, model.e, "\n")
-        frm:write("f", sep, model.f, "\n")
-    elseif model.type == "sphere" then
-        frm:write("r", sep, model.r, "\n")
-    end
-    frm:write("proj", sep, self.proj.name, "\n")
-    frm:write("lon", sep, deg(self.proj.lon), "\n")
-    frm:write("lat", sep, deg(self.proj.lat), "\n")
-    frm:write("x0", sep, self.bbox.x0, "\n")
-    frm:write("y0", sep, self.bbox.y0, "\n")
-    frm:write("x1", sep, self.bbox.x1, "\n")
-    frm:write("y1", sep, self.bbox.y1, "\n")
-    frm:close()
-end
-
-local function new_frame(proj, bbox)
-    local self = setmetatable({}, Frame)
-    self.proj = proj
-    self.bbox = bbox
-    self.w, self.h, self.s = bbox.x1-bbox.x0, bbox.y1-bbox.y0, 1
-    return self
-end
-
-local function load_frame(fname)
-    local self = setmetatable({}, Frame)
-    local frm = io.open(fname, "r")
-    local function get(field)
-        local line = frm:read()
-        local got = line:sub(1, #field)
-        assert(got == field, "expected field "..field.." but got "..got)
-        return line:sub(#field+#sep+1)
-    end
-    local model = {}
-    model.type = get "type"
-    if model.type == "ellipsoid" then
-        model.a = tonumber(get "a")
-        model.b = tonumber(get "b")
-        model.e = tonumber(get "e")
-        model.f = tonumber(get "f")
-    elseif model.type == "sphere" then
-        model.r = tonumber(get "r")
-    end
-    local proj_name = get "proj"
-    local proj_lon = tonumber(get "lon")
-    local proj_lat = tonumber(get "lat")
-    proj = Proj(proj_name, {proj_lon, proj_lat}, model.r)
-    local bbox = {}
-    bbox.x0 = tonumber(get "x0")
-    bbox.y0 = tonumber(get "y0")
-    bbox.x1 = tonumber(get "x1")
-    bbox.y1 = tonumber(get "y1")
-    frm:close()
-    return new_frame(proj, bbox)
-end
-
-return {
-    distance=distance, bbox=bbox, centroid=centroid, Proj=Proj,
-    new_frame=new_frame, load_frame=load_frame
-}

diff --git a/poly.lua b/poly.lua
deleted file mode 100644
index 01d132d..0000000
--- a/poly.lua
+++ /dev/null
@@ -1,197 +0,0 @@
-local util = require "util"
-
-local pi = math.pi
-local sqrt2 = math.sqrt(2)
-
-local function dashed(points, pattern)
-    local x0, y0, x1, y1
-    local cx, cy
-    local i, j
-    local d, h
-    local polylines = {}
-    local polyline
-    local draw = true
-    x0, y0 = unpack(points[1])
-    polyline = {{x0, y0}}
-    i = 2
-    x1, y1 = unpack(points[i])
-    j = 1
-    d = pattern[j]
-    while true do
-        h = util.hypot(x1-x0, y1-y0)
-        if d < h then
-            cx = x0 + (x1-x0)*d/h
-            cy = y0 + (y1-y0)*d/h
-            if draw then
-                table.insert(polyline, {cx, cy})
-                table.insert(polylines, polyline)
-            else
-                polyline = {{cx, cy}}
-            end
-            x0, y0 = cx, cy
-            draw = not draw
-            if j < #pattern then j = j + 1 else j = 1 end
-            d = pattern[j]
-        else
-            if draw then
-                table.insert(polyline, {x1, y1})
-            end
-            d = d - h
-            if i < #points then i = i + 1 else break end
-            x0, y0 = x1, y1
-            x1, y1 = unpack(points[i])
-        end
-    end
-    if draw then
-        table.insert(polyline, points[i])
-        table.insert(polylines, polyline)
-    end
-    return polylines
-end
-
--- convert bezier curve {{ax, ay}, {bx, by}, {cx, cy}} to polyline
-local function bezier(curve)
-    local a, b, c, h
-    local ax, ay, bx, by, cx, cy
-    local dx, dy, ex, ey, fx, fy
-    local points = {}
-    local stack = {curve}
-    while #stack > 0 do
-        a, b, c = unpack(table.remove(stack))
-        ax, ay = unpack(a)
-        bx, by = unpack(b)
-        cx, cy = unpack(c)
-        h = math.abs((ax-cx)*(by-ay)-(ax-bx)*(cy-ay))/util.hypot(cx-ax, cy-ay)
-        if h > 1 then -- split curve
-            dx, dy = (ax+bx)/2, (ay+by)/2
-            fx, fy = (bx+cx)/2, (by+cy)/2
-            ex, ey = (dx+fx)/2, (dy+fy)/2
-            table.insert(stack, {{ex, ey}, {fx, fy}, {cx, cy}})
-            table.insert(stack, {{ax, ay}, {dx, dy}, {ex, ey}})
-        else -- add point to polyline
-            table.insert(points, {cx, cy})
-        end
-    end
-    return points
-end
-
--- convert a sequence of linked Bézier curves to a polyline
-local function unfold(control_points)
-    local s, a, b
-    local px, py, qx, qy, rx, ry
-    s = control_points[1]
-    b = control_points[2]
-    local points = {s}
-    local sub
-    if b[3] then
-        table.insert(points, b)
-        s = b
-    end
-    for i = 3, #control_points do
-        a = b
-        b = control_points[i]
-        if a[3] then
-            if b[3] then
-                table.insert(points, b)
-                s = b
-            end
-        else
-            if b[3] then
-                px, py, _ = unpack(s)
-                qx, qy, _ = unpack(a)
-                rx, ry, _ = unpack(b)
-                sub = bezier({{px, py}, {qx, qy}, {rx, ry}})
-                for i, p in ipairs(sub) do
-                    table.insert(points, p)
-                end
-                s = b
-            else
-                px, py, _ = unpack(s)
-                qx, qy, _ = unpack(a)
-                rx, ry, _ = unpack(b)
-                rx, ry = (qx+rx)/2, (qy+ry)/2
-                sub = bezier({{px, py}, {qx, qy}, {rx, ry}})
-                for i, p in ipairs(sub) do
-                    table.insert(points, p)
-                end
-                s = {rx, ry, true}
-            end
-        end
-    end
-    return points
-end
-
--- bezigon that approximates a circle
-local function bcircle(x, y, r)
-    local a = r * (sqrt2-1)
-    return {
-        {x  , y+r, true},
-        {x+a, y+r, false},
-        {x+r, y+a, false},
-        {x+r, y-a, false},
-        {x+a, y-r, false},
-        {x-a, y-r, false},
-        {x-r, y-a, false},
-        {x-r, y+a, false},
-        {x-a, y+r, false},
-        {x  , y+r, true}
-    }
-end
-
--- regular polygon
-local function ngon(x, y, r, n, mina, maxa)
-    mina = mina or 0
-    maxa = maxa or 2 * pi
-    local a = 2 * pi / n -- angle between points
-    local pgon = {}
-    local px, py
-    local pa = mina
-    while pa < maxa do
-        px = x + r * math.cos(pa)
-        py = y + r * math.sin(pa)
-        table.insert(pgon, {px, py})
-        pa = pa + a
-    end
-    px = x + r * math.cos(maxa)
-    py = y + r * math.sin(maxa)
-    table.insert(pgon, {px, py})
-    return pgon
-end
-
--- approximate an arc between angles mina and maxa
-local function parc(x, y, r, mina, maxa)
-    local h = 0.5 -- maximum radius-apothem allowed
-    local n = math.ceil(pi / math.acos(1 - h/r)) -- # of sides
-    return ngon(x, y, r, n, mina, maxa)
-end
-
--- regular polygon that approximates a circle
-local function pcircle(x, y, r)
-    return parc(x, y, r, 0, 2 * pi)
-end
-
-local function arrow_head(x0, y0, x1, y1, w, h)
-    local dx, dy = x1-x0, y1-y0
-    local a = math.atan2(dy, dx)    -- line angle
-    local b = math.atan2(-dx, dy)   -- perpendicular angle
-    local mx, my = math.cos(a), math.sin(a)
-    local nx, ny = math.cos(b), math.sin(b)
-    local bx, by = x1 - mx * h, y1 - my * h     -- back of arrow
-    local lx, ly = bx - nx * w/2, by - ny * w/2 -- left point
-    local rx, ry = bx + nx * w/2, by + ny * w/2 -- right point
-    bx, by = bx + mx * h/2, by + my * h/2       -- off-curve point
-    local control_points = {
-      {lx, ly, true},
-      {x1, y1, true},
-      {rx, ry, true},
-      {bx, by, false},
-      {lx, ly, true},
-    }
-    return unfold(control_points)
-    --~ return control_points
-end
-
-return {
-  dashed=dashed, unfold=unfold, bcircle=bcircle, ngon=ngon,
-  parc=parc, pcircle=pcircle, arrow_head=arrow_head
-}

diff --git a/shp.lua b/shp.lua
deleted file mode 100644
index 55a37f5..0000000
--- a/shp.lua
+++ /dev/null
@@ -1,375 +0,0 @@
-local bit = require "bit"
-
-local util = require "util"
-local bio = require "bio"
-
-local SF = {}
-SF.__index = SF
-
-function SF:read_dbf()
-    local fp = io.open(self.path .. ".dbf", "rb")
-    local version = bit.band(bio.read_byte(fp), 0x07)
-    if version ~= 3 then
-        error("only DBF version 5 is supported")
-    end
-    self.year = 1900 + bio.read_byte(fp)
-    self.month = bio.read_byte(fp)
-    self.day = bio.read_byte(fp)
-    self.nrecs = bio.read_leu32(fp)
-    fp:seek("cur", 24)
-    local reclen = 0
-    local fields = {}
-    local byte = bio.read_byte(fp)
-    while byte ~= 0x0D do
-        local field_name = util.rtrim(string.char(byte) .. fp:read(10), "\000")
-        local field_type = fp:read(1)
-        fp:seek("cur", 4)
-        local field_length = bio.read_byte(fp)
-        reclen = reclen + field_length
-        local field_dec_count = bio.read_byte(fp)
-        local field = {
-          name=field_name,
-          type=field_type,
-          length=field_length,
-          dec_count=field_dec_count
-        }
-        table.insert(fields, field)
-        fp:seek("cur", 14)
-        byte = bio.read_byte(fp)
-    end
-    local tab = {}
-    for i = 1, self.nrecs do
-        if fp:read(1) == "*" then
-            fp:seek("cur", reclen)
-        else
-            local row = {}
-            for j = 1, #fields do
-                local field = fields[j]
-                local cell = util.rtrim(fp:read(field.length), " ")
-                if field.type == "F" or field.type == "N" then
-                    cell = tonumber(cell) or 0
-                end
-                table.insert(row, cell)
-            end
-            table.insert(tab, row)
-        end
-    end
-    self.fields = fields
-    self.tab = tab
-    io.close(fp)
-end
-
-function SF:search(field_name, value)
-    local col
-    for i = 1, #self.fields do
-        if self.fields[i].name == field_name then
-            col = i
-            break
-        end
-    end
-    if col ~= nil then
-        for i = 1, self.nrecs do
-            if self.tab[i][col] == value then
-                return i
-            end
-        end
-    end
-end
-
-function SF:print_summary(n)
-    n = n or 5
-    local sep = ":"
-    for i = 1, #self.fields do
-        local field = self.fields[i]
-        local row = {field.name}
-        for j = 1, n do
-            table.insert(row, self.tab[j][i])
-        end
-        print(table.concat(row, sep))
-    end
-    print("records".. sep .. #self.tab)
-    print("shape".. sep .. self.header.shape)
-end
-
-function SF:tab2csv(sep, fp)
-    -- TODO: refactor to use table.concat(str_list, sep)
-    sep = sep or ":"
-    fp = fp or io.output()
-    for i = 1, #self.fields do
-        if i > 1 then
-            fp:write(sep)
-        end
-        fp:write(self.fields[i].name)
-    end
-    fp:write("\n")
-    for i = 1, #self.tab do
-        for j = 1, #self.fields do
-            if j > 1 then
-                fp:write(sep)
-            end
-            local cell = self.tab[i][j]
-            fp:write(cell)
-        end
-        fp:write("\n")
-    end
-end
-
-local shape_name = {
-    [ 0] = "null",          [31] = "multipatch",
-    [ 1] = "point",   [ 3] = "polyline",   [ 5] = "polygon",   [ 8] = "multipoint",
-    [11] = "point-z", [13] = "polyline-z", [15] = "polygon-z", [18] = "multipoint-z",
-    [21] = "point-m", [23] = "polyline-m", [25] = "polygon-m", [28] = "multipoint-m",
-}
-
-local function read_bbox(fp)
-    local bbox = {}
-    bbox.xmin = bio.read_led64(fp)
-    bbox.ymin = bio.read_led64(fp)
-    bbox.xmax = bio.read_led64(fp)
-    bbox.ymax = bio.read_led64(fp)
-    return bbox
-end
-
--- SHX and SHP share the same header structure!
-local function read_sf_header(fp)
-    local file_code = bio.read_bei32(fp)
-    if file_code ~= 9994 then
-        error("does not seem to be a shapefile (invalid FileCode)")
-    end
-    fp:seek("cur", 20) -- unused, always all-zero
-    local file_len = bio.read_bei32(fp)
-    local version = bio.read_lei32(fp)
-    if version ~= 1000 then
-        error("invalid shapefile version")
-    end
-    local header = {}
-    header.shape = shape_name[bio.read_lei32(fp)]
-    header.bbox = read_bbox(fp)
-    header.zmin = bio.read_led64(fp)
-    header.zmax = bio.read_led64(fp)
-    header.mmin = bio.read_led64(fp)
-    header.mmax = bio.read_led64(fp)
-    return header
-end
-
-function SF:read_shx()
-    local fp = io.open(self.path .. ".shx", "rb")
-    self.header = read_sf_header(fp)
-    self.index = {}
-    for i = 1, self.nrecs do
-        local offset = bio.read_bei32(fp)
-        local length = bio.read_bei32(fp)
-        table.insert(self.index, {offset=offset, length=length})
-    end
-    io.close(fp)
-end
-
-function SF:read_record_header(index)
-    local fp = self.fp
-    fp:seek("set", self.index[index].offset * 2)
-    local num = bio.read_bei32(fp)
-    local len = bio.read_bei32(fp)
-    local shape = shape_name[bio.read_lei32(fp)]
-    return num, len, shape
-end
-
-function SF:read_bbox(index)
-    local shape_name = self.header.shape
-    assert(
-        util.startswith(shape_name, "polyline") or
-        util.startswith(shape_name, "polygon") or
-        util.startswith(shape_name, "multipoint"),
-        ("%s shape doesn't have bbox"):format(shape_name)
-    )
-    local num, len, shape = self:read_record_header(index)
-    if shape == "null" then
-        return nil
-    end
-    return read_bbox(self.fp)
-end
-
-function SF:get_point(index, proj)
-    assert(util.startswith(self.header.shape, "point"), "type mismatch")
-    local fp = self.fp
-    local num, len, shape = self:read_record_header(index)
-    if shape == "null" then
-        return nil
-    end
-    local x = bio.read_led64(fp)
-    local y = bio.read_led64(fp)
-    if proj ~= nil then
-        x, y = proj:map(x, y)
-    end
-    return {x, y}
-end
-
-function SF:get_polys(index, proj)
-    assert(util.startswith(self.header.shape, "poly"), "type mismatch")
-    local fp = self.fp
-    local num, len, shape = self:read_record_header(index)
-    if shape == "null" then
-        return nil
-    end
-    local bbox = read_bbox(fp)
-    local nparts = bio.read_lei32(fp)
-    --~ print(nparts)
-    local npoints = bio.read_lei32(fp)
-    local total = bio.read_lei32(fp) -- first index is always zero?
-    local lengths = {}
-    for i = 2, nparts do
-        local index = bio.read_lei32(fp)
-        table.insert(lengths, index - total)
-        total = index
-    end
-    table.insert(lengths, npoints - total)
-    local i = 0
-    return bbox, lengths, function()
-        i = i + 1
-        if i <= #lengths then
-            local j = 0
-            local length = lengths[i]
-            return function()
-                j = j + 1
-                if j <= length then
-                    local x = bio.read_led64(fp)
-                    local y = bio.read_led64(fp)
-                    if proj ~= nil then
-                        x, y = proj:map(x, y)
-                    end
-                    return {x, y}
-                end
-            end
-        end
-    end
-end
-
-function SF:close()
-    io.close(self.fp)
-end
-
-function SF:save_cache(fname, k, proj, scale, filter)
-    local indices = {}
-    for i = 1, #self.tab do
-        local row = self.tab[i]
-        local rec = {}
-        for j = 1, #self.fields do
-            rec[self.fields[j].name] = row[j]
-        end
-        local key = filter(rec)
-        if key ~= nil then
-            table.insert(indices, {i, key:sub(1, 16)})
-        end
-    end
-    local cache = io.open(fname, "w")
-    bio.write_beu32(cache, util.round(1000 / scale))
-    bio.write_byte(cache, k)
-    for i = 1, #indices do
-        local index, key = unpack(indices[i])
-        cache:write(key, string.rep("\0", 20 - #key))
-    end
-    cache:write(string.rep("\0", 20)) -- end list of entries
-    local function geo2grid(p)
-        local x, y = unpack(p)
-        x, y = proj:map(x, y)
-        x, y = util.round(x * scale), util.round(y * scale)
-        return x, y
-    end
-    for i = 1, #indices do
-        local index, key = unpack(indices[i])
-        local offset = cache:seek()
-        cache:seek("set", i * 20 + 1)
-        bio.write_beu32(cache, offset)
-        cache:seek("set", offset)
-        local bb, lens, polys = self:get_polys(index)
-        bio.write_beu16(cache, #lens)
-        for poly in polys do
-            local px, py = geo2grid(poly())
-            bio.write_bei16(cache, px)
-            bio.write_bei16(cache, py)
-            local qx, qy = geo2grid(poly())
-            local rice = bio.rice_w(cache, k)
-            for point in poly do
-                local rx, ry = geo2grid(point)
-                if (px-rx)*(qy-py) - (px-qx)*(ry-py) ~= 0 then
-                    local dx, dy = qx-px, qy-py
-                    rice:put_signed(dx)
-                    rice:put_signed(dy)
-                    px, py = qx, qy
-                end
-                qx, qy = rx, ry
-            end
-            rice:put_signed(0)
-            rice:put_signed(0)
-            rice:flush()
-        end
-    end
-    cache:close()
-end
-
-local function open_shapefile(path)
-    local self = setmetatable({path=path}, SF)
-    self:read_dbf()
-    self:read_shx()
-    self.fp = io.open(path .. ".shp", "rb")
-    return self
-end
-
-local Cache = {}
-Cache.__index = Cache
-
-function Cache:keys()
-    local cache = self.fp
-    local offset = 5
-    return function()
-        cache:seek("set", offset)
-        offset = offset + 20
-        local ckey = cache:read(16)
-        if ckey:byte() ~= 0 then
-            return ckey:sub(1, ckey:find("\0")-1)
-        end
-    end
-end
-
-function Cache:get_polys(key)
-    local cache = self.fp
-    cache:seek("set", 5)
-    local ckey = cache:read(16)
-    local offset = -1
-    while ckey:byte() ~= 0 do
-        if ckey:sub(1, ckey:find("\0")-1) == key then
-            offset = bio.read_beu32(cache)
-            break
-        end
-        cache:seek("cur", 4)
-        ckey = cache:read(16)
-    end
-    assert(offset > 0, ("key '%s' not found in cache"):format(key))
-    cache:seek("set", offset)
-    local npolys = bio.read_beu16(cache)
-    return function()
-        if npolys > 0 then
-            npolys = npolys - 1
-            local ox, oy = bio.read_bei16(cache), bio.read_bei16(cache)
-            local rice = bio.rice_r(cache, self.k)
-            local x, y = 0, 0
-            return function()
-                if x ~= ox or y ~= oy then
-                    x, y = ox, oy
-                    local dx, dy = rice:get_signed(), rice:get_signed()
-                    ox, oy = ox+dx, oy+dy
-                    return {x * self.s, y * self.s}
-                end
-            end
-        end
-    end
-end
-
-local function load_cache(fname)
-    local self = setmetatable({}, Cache)
-    self.fp = io.open(fname, "r")
-    self.s = bio.read_beu32(self.fp) / 1000
-    self.k = bio.read_byte(self.fp)
-    return self
-end
-
-return {open_shapefile=open_shapefile, load_cache=load_cache}

diff --git a/surf.lua b/surf.lua
deleted file mode 100644
index a175e6b..0000000
--- a/surf.lua
+++ /dev/null
@@ -1,266 +0,0 @@
-local ffi = require "ffi"
-local bit = require "bit"
-
-local util = require "util"
-
-local bnot = bit.bnot
-local bor, band = bit.bor, bit.band
-local lshift, rshift =  bit.lshift,  bit.rshift
-
-local Surf = {}
-
-function Surf:inside(x, y)
-    return x >= 0 and x < self.w and y >= 0 and y < self.h
-end
-
-function Surf:hline(x, y, w, v)
-    for i = x, x+w-1 do
-        self:pset(i, y, v)
-    end
-end
-
-function Surf:vline(x, y, h, v)
-    for i = y, y+h-1 do
-        self:pset(x, i, v)
-    end
-end
-
-function Surf:disk_int(cx, cy, r, v)
-    if r == 0 then
-        self:pset(cx, cy, v)
-        return
-    end
-    local x, y, d = r, 0, 1-r
-    while x >= y do
-        self:hline(cx-x, cy+y, 2*x, v)
-        self:hline(cx-y, cy+x, 2*y, v)
-        self:hline(cx-x, cy-y, 2*x, v)
-        self:hline(cx-y, cy-x, 2*y, v)
-        y = y + 1
-        if d <= 0 then
-            d = d + 2*y + 1
-        else
-            x = x - 1
-            d = d + 2*(y-x) + 1
-        end
-    end
-end
-
-function Surf:disk(cx, cy, r, v)
-    self:disk_int(util.round(cx), util.round(cy), r, v)
-end
-
-function Surf:line(x0, y0, x1, y1, v, r)
-    r = r or 0
-    local dx, dy = x1-x0, y1-y0
-    local n = math.max(math.abs(dx), math.abs(dy))
-    local sx, sy = dx/n, dy/n
-    local x, y = x0, y0
-    self:disk_int(math.floor(x), math.floor(y), r, v)
-    for i = 1, n do
-        x = x + sx
-        y = y + sy
-        self:disk_int(math.floor(x), math.floor(y), r, v)
-    end
-end
-
-function Surf:polyline(points, v, r)
-    points = util.func_iter(points)
-    local x0, y0, x1, y1
-    x0, y0 = unpack(points())
-    for point in points do
-        x1, y1 = unpack(point)
-        self:line(x0, y0, x1, y1, v, r)
-        x0, y0 = x1, y1
-    end
-end
-
-function Surf:polylines(polylines, v, r)
-    polylines = util.func_iter(polylines)
-    for points in polylines do
-        self:polyline(points, v, r)
-    end
-end
-
-local function cross_comp(a, b)
-    return a[1] < b[1]
-end
-
-function Surf:scan(points)
-    points = util.func_iter(points)
-    if not self.scans then
-        self.scans = {}
-        for i = 0, self.h-1 do
-            self.scans[i] = {}
-        end
-    end
-    local x0, y0, x1, y1
-    local ax, ay, bx, by -- same line as above, but enforce ay < by
-    local sign
-    x0, y0 = unpack(points())
-    for point in points do
-        x1, y1 = unpack(point)
-        if y1 ~= y0 then
-            if y0 < y1 then
-                ax, ay, bx, by = x0, y0, x1, y1
-                sign =  1
-            else
-                ax, ay, bx, by = x1, y1, x0, y0
-                sign = -1
-            end
-            local slope = (bx-ax) / (by-ay)
-            ay, by = util.round(ay), util.round(by)
-            while ay < by do
-                if ay >= 0 and ay < self.h then
-                    table.insert(self.scans[ay], {ax, sign})
-                end
-                ax = ax + slope
-                ay = ay + 1
-            end
-        end
-        x0, y0 = x1, y1
-    end
-end
-
-function Surf:fill_scans(v)
-    local scan, wind
-    local x, sign
-    local ax, bx
-    for i = 0, self.h-1 do
-        scan = self.scans[i]
-        table.sort(scan, cross_comp)
-        wind = 0
-        for j, cross in ipairs(scan) do
-            x, sign = unpack(cross)
-            if wind == 0 then
-                ax = math.floor(x)
-            end
-            wind = wind + sign
-            if wind == 0 then
-                bx = math.ceil(x)
-                self:hline(ax, i, bx-ax, v)
-            end
-        end
-    end
-    self.scans = nil
-end
-
-function Surf:polygon(points, v)
-    self:scan(points)
-    self:fill_scans(v)
-end
-
-function Surf:polygons(polygons, v)
-    polygons = util.func_iter(polygons)
-    for points in polygons do
-        self:scan(points)
-    end
-    self:fill_scans(v)
-end
-
-function Surf:blit(x, y, surf, sx, sy, w, h)
-    for j = 0, h-1 do
-        for i = 0, w-1 do
-            self:pset(x+i, y+j, surf:pget(sx+i, sy+j))
-        end
-    end
-end
-
-function Surf:compose(layers)
-    for y = 0, self.h-1 do
-        for x = 0, self.w-1 do
-            for _, layer in ipairs(layers) do
-                local p = layer.map[layer.surf:pget(x, y)+1]
-                if p >= 0 then
-                    self:pset(x, y, p)
-                end
-            end
-        end
-    end
-end
-
-local BitMap = {}
-
-function BitMap:fill(v)
-    ffi.fill(self.p, self.t * self.h, v)
-end
-
-function BitMap:_index_shift_mask(x, y)
-    local index = y * self.t + math.floor(x / 8)
-    local shift = (7-x) % 8
-    local mask = lshift(1, shift)
-    return index, shift, mask
-end
-
-function BitMap:pget(x, y)
-    if not self:inside(x, y) then return 0 end
-    local index, shift, mask = self:_index_shift_mask(x, y)
-    return rshift(band(self.p[index], mask), shift)
-end
-
-function BitMap:pset(x, y, v)
-    if not self:inside(x, y) then return end
-    local index, shift, mask = self:_index_shift_mask(x, y)
-    local byte = self.p[index]
-    if v > 0 then
-        byte = bor(byte, mask)
-    else
-        byte = band(byte, bnot(mask))
-    end
-    self.p[index] = byte
-end
-
-function BitMap:save_pbm(fname)
-    local pbm = io.open(fname, "wb")
-    pbm:write("P4\n", self.w, " ", self.h, "\n")
-    local row = self.p + 0
-    for y = 0, self.h-1 do
-        pbm:write(ffi.string(row, self.t))
-        row = row + self.t
-    end
-    pbm:close()
-end
-
-setmetatable(BitMap, {__index=Surf})
-
-local function new_bitmap(w, h)
-    local self = setmetatable({w=w, h=h}, {__index=BitMap})
-    self.t = math.ceil(w / 8) -- stride, i.e., bytes/row
-    self.p = ffi.new("unsigned char[?]", self.t * self.h)
-    return self
-end
-
-local ByteMap = {}
-
-function ByteMap:fill(v)
-    ffi.fill(self.p, self.w * self.h, v)
-end
-
-function ByteMap:pget(x, y)
-    if not self:inside(x, y) then return 0 end
-    return self.p[y*self.w+x]
-end
-
-function ByteMap:pset(x, y, v)
-    if not self:inside(x, y) then return end
-    self.p[y*self.w+x] = v
-end
-
-function ByteMap:save_ppm(fname, colors)
-    local ppm = io.open(fname, "wb")
-    ppm:write("P6\n", self.w, " ", self.h, "\n", 255, "\n")
-    for i = 0, self.w*self.h-1 do
-        ppm:write(string.char(unpack(colors[self.p[i]+1])))
-    end
-    ppm:close()
-end
-
-setmetatable(ByteMap, {__index=Surf})
-
-local function new_bytemap(w, h)
-    local self = setmetatable({w=w, h=h}, {__index=ByteMap})
-    self.p = ffi.new("unsigned char[?]", self.w * self.h)
-    return self
-end
-
-return {new_bitmap=new_bitmap, new_bytemap=new_bytemap}

diff --git a/ttf.lua b/ttf.lua
deleted file mode 100644
index 69bad1f..0000000
--- a/ttf.lua
+++ /dev/null
@@ -1,624 +0,0 @@
-local bit = require "bit"
-
-local aff = require "aff"
-local poly = require "poly"
-
-local bnot = bit.bnot
-local bor, band = bit.bor, bit.band
-local lshift, rshift =  bit.lshift,  bit.rshift
-
-local function utf8to32(utf8str)
-    assert(type(utf8str) == "string")
-    local res, seq, val = {}, 0, nil
-    for i = 1, #utf8str do
-        local c = string.byte(utf8str, i)
-        if seq == 0 then
-            table.insert(res, val)
-            seq = c < 0x80 and 1 or c < 0xE0 and 2 or c < 0xF0 and 3 or
-                  c < 0xF8 and 4 or --c < 0xFC and 5 or c < 0xFE and 6 or
-                  error("invalid UTF-8 character sequence")
-            val = band(c, 2^(8-seq) - 1)
-        else
-            val = bor(lshift(val, 6), band(c, 0x3F))
-        end
-        seq = seq - 1
-    end
-    table.insert(res, val)
-    return res
-end
-
--- Note: TrueType uses big endian for everything.
-
-local function uint(s)
-    local x = 0
-    for i = 1, #s do
-        x = x * 256 + s:byte(i)
-    end
-    return x
-end
-
-local function int(s)
-    local x = uint(s)
-    local p = 2^(#s*8)
-    if x >= p/2 then x = x - p end
-    return x
-end
-
-local function str(x)
-    local s = ""
-    while x > 0 do
-        s = string.char(x % 256) .. s
-        x = math.floor(x / 256)
-    end
-    return s
-end
-
-local Face = {}
-Face.__index = Face
-
-function Face:log(s)
-    if self.dbg then
-        io.stderr:write(s .. "\n")
-    end
-end
-
-function Face:str(n)    return self.fp:read(n) end
-function Face:uint8()   return uint(self.fp:read(1)) end
-function Face:uint16()  return uint(self.fp:read(2)) end
-function Face:uint32()  return uint(self.fp:read(4)) end
-function Face:uint64()  return uint(self.fp:read(8)) end
-function Face:int8()    return  int(self.fp:read(1)) end
-function Face:int16()   return  int(self.fp:read(2)) end
-function Face:int32()   return  int(self.fp:read(4)) end
-function Face:int64()   return  int(self.fp:read(8)) end
-
-function Face:goto(tag, offset)
-    self.fp:seek("set", self.dir[tag].offset + (offset or 0))
-end
-
-function Face:getpos()
-    return self.fp:seek()
-end
-
-function Face:setpos(pos)
-    self.fp:seek("set", pos)
-end
-
-function Face:offset()
-    local scaler_type = self:uint32()
-    assert(scaler_type == 0x74727565 or scaler_type == 0x00010000,
-        ("invalid scaler type for TrueType: 0x%08X"):format(scaler_type))
-    local num_tables = self:uint16()
-    --  The entries for search_range, entry_selector and range_shift are used to
-    -- facilitate quick binary searches of the table directory. Unless a font 
-    -- has a large number of tables, a sequential search will be fast enough.
-    local search_range = self:uint16()
-    local entry_selector = self:uint16()
-    local range_shift = self:uint16()
-    self.num_tables = num_tables
-end
-
-function Face:directory()
-    local dir = {}
-    for i = 1, self.num_tables do
-        local tag = self:str(4)
-        local checksum = self:uint32()
-        local offset = self:uint32()
-        local length = self:uint32()
-        -- TODO: verify checksums
-        dir[tag] = {checksum=checksum, offset=offset, length=length}
-    end
-    self.dir = dir
-end
-
-function Face:head()
-    self:goto("head")
-    local version = self:uint32()
-    assert(version == 0x00010000, ("invalid version: 0x%08X"):format(version))
-    local revision = self:uint32()
-    local checksum_adj = self:uint32()
-    local magic = self:uint32()
-    assert(magic == 0x5F0F3CF5, ("invalid magic: 0x%08X"):format(magic))
-    local flags = self:uint16()
-    self.units_per_em = self:uint16()
-    local created = self:int64()
-    local modified = self:int64()
-    local xmin = self:int16()
-    local ymin = self:int16()
-    local xmax = self:int16()
-    local ymax = self:int16()
-    local mac_style = self:uint16()
-    local lowest_rec_ppem = self:uint16()
-    local direction_hint = self:int16()
-    self.index_to_loc_fmt = self:int16()
-    local glyph_data_fmt = self:int16()
-end
-
-function Face:maxp()
-    self:goto("maxp")
-    local version = self:uint32()
-    if version == 0x00005000 then
-        self.num_glyphs = self:uint16()
-    elseif version == 0x00010000 then
-        self.num_glyphs = self:uint16()
-        local max_points = self:uint16()
-        local max_contours = self:uint16()
-        local max_composite_points = self:uint16()
-        local max_composite_contours = self:uint16()
-        local max_zones = self:uint16()
-        local max_twilight_points = self:uint16()
-        local max_storage = self:uint16()
-        local max_function_defs = self:uint16()
-        local max_instruction_defs = self:uint16()
-        local max_stack_elements = self:uint16()
-        local max_size_of_instructions = self:uint16()
-        local max_component_elements = self:uint16()
-        local max_component_depth = self:uint16()
-    else
-        error(("invalid maxp version: 0x%08X"):format(version))
-    end
-end
-
-function Face:cmap()
-    self:goto("cmap")
-    local version = self:uint16()
-    assert(version == 0, ("invalid cmap version: %d"):format(version))
-    local num_subtables = self:uint16()
-    local encoding, suboffset
-    local ok = false
-    for i = 1, num_subtables do
-        local platform_id = self:uint16()
-        encoding = self:uint16()
-        suboffset = self:uint32()
-        if platform_id == 0 then        -- platform == Unicode
-            ok = true
-            break
-        elseif platform_id == 3 then    -- platform == Microsoft
-            if encoding == 10 or encoding == 1 then
-                ok = true
-                break
-            end
-        end
-    end
-    if not ok then
-        error(("could not find Unicode cmap in %d subtables"):format(num_subtables))
-    end
-    self:goto("cmap", suboffset)
-    self:subcmap()
-end
-
-function Face:subcmap()
-    local format = self:uint16()
-    local segs = {}
-    local gia = {}
-    if format == 4 then
-        local length = self:uint16()
-        local language = self:uint16()
-        assert(language == 0, ("invalid subcmap language: %d"):format(language))
-        local seg_count = self:uint16() / 2
-        local search_range = self:uint16()
-        local entry_selector = self:uint16()
-        local range_shift = self:uint16()
-        for i = 1, seg_count do
-            segs[i] = {end_code=self:uint16()}
-        end
-        local last = segs[seg_count].end_code
-        assert(last == 0xFFFF, ("invalid subcmap last end code: %d"):format(last))
-        local pad = self:uint16()
-        assert(pad == 0, ("invalid subcmap reserved pad: %d"):format(pad))
-        for i = 1, seg_count do
-            segs[i].start_code = self:uint16()
-        end
-        for i = 1, seg_count do
-            segs[i].id_delta = self:uint16()
-        end
-        for i = 1, seg_count do
-            segs[i].id_range_offset = self:uint16()
-        end
-        local gia_len = (length - (16+8*seg_count)) / 2
-        for i = 1, gia_len do
-            gia[i] = self:uint16()
-        end
-    -- TODO: support other formats, specially 6 and 12
-    else
-        error(("unsupported subcmap format: %d"):format(format))
-    end
-    self.segs, self.gia = segs, gia
-end
-
-function Face:os2()
-    self:goto("OS/2")
-    local version = self:uint16()
-    assert(version >= 2, ("invalid OS/2 version: %u"):format(version))
-    self.avg_char_width = self:int16()
-    self.weight_class = self:uint16()
-    self.width_class = self:uint16()
-    self.fp:seek("cur", 78)
-    self.x_height = self:int16()
-    self.cap_height = self:int16()
-    self.default_char = self:uint16()
-    self.break_char = self:uint16()
-end
-
-function Face:hhea()
-    self:goto("hhea")
-    local versionH = self:uint16()
-    local versionL = self:uint16()
-    assert(versionH == 1 and versionL == 0,
-        ("invalid hhea version: %d.%d"):format(versionH, versionL))
-    self.ascent  = self:int16()
-    self.descent = self:int16()
-    local line_gap = self:int16()
-    local advance_width_max = self:uint16()
-    local min_left_side_bearing  = self:int16()
-    local min_right_side_bearing = self:int16()
-    local x_max_extent = self:int16()
-    local caret_slope_rise = self:int16()
-    local caret_slope_run  = self:int16()
-    local caret_offset     = self:int16()
-    for i = 1, 4 do
-        local reserved = self:uint16()
-        assert(reserved == 0, "nonzero reserved field in hhea")
-    end
-    local metric_data_format = self:int16()
-    assert(metric_data_format == 0,
-        ("invalid metric data format: %d"):format(metric_data_format))
-    self.nlong_hor_metrics = self:uint16()
-end
-
-function Face:hmetrics(id)
-    local n = self.nlong_hor_metrics -- for readability of expressions below
-    local advance, bearing
-    if id < n then
-        self:goto("hmtx", 4*id)
-        advance = self:uint16()
-        bearing = self:int16()
-    else
-        self:goto("hmtx", 4*(n-1))
-        advance = self:uint16()
-        self:goto("hmtx", 4*n+2*(id-n))
-        bearing = self:int16()
-    end
-    return advance, bearing
-end
-
-function Face:kern()
-    self:goto("kern")
-    local version = self:uint16()
-    assert(version == 0, ("invalid kern table version: %d"):format(version))
-    local ntables = self:uint16()
-    for i = 1, ntables do
-        local version = self:uint16()
-        local length = self:uint16()
-        local format = self:uint8() -- usually 0
-        local coverage = self:uint8()
-        local horizontal   = band(coverage, 2^0) > 0 -- usually true
-        local minimum      = band(coverage, 2^1) > 0 -- usually false (kerning)
-        local cross_stream = band(coverage, 2^2) > 0 -- usually false (regular)
-        local override     = band(coverage, 2^3) > 0 -- usually false (accumulate)
-        assert(band(coverage, 0xF0) == 0, "invalid coverage bits set")
-        if format == 0 then
-            self.num_kernings = self:uint16()
-            local search_range = self:uint16()
-            local entry_selector = self:uint16()
-            local range_shift = self:uint16()
-            local kerning = {}
-            for j = 1, self.num_kernings do
-                -- glyph indices (left * 2^16 + right)
-                local key  = self:uint32()
-                -- kerning value
-                local value = self:int16()
-                kerning[key] = value
-            end
-            self.kerning = kerning
-        else
-            self:log(("unsupported kerning table format: %d"):format(format))
-        end
-    end
-end
-
-function Face:get_kerning(left_id, right_id)
-    return self.kerning[left_id * 2^16 + right_id] or 0
-end
-
--- Convert a character code to its glyph id.
-function Face:char_index(code)
-    local i = 1
-    while code > self.segs[i].end_code do i = i + 1 end
-    if self.segs[i].start_code > code then return 0 end
-    local iro = self.segs[i].id_range_offset
-    if iro == 0 then
-        return (code + self.segs[i].id_delta) % 0x10000
-    else
-        local idx = iro + 2 * (code - self.segs[i].start_code)
-        idx = idx - (#self.segs - i + 1) * 2
-        local id = self.gia[idx/2+1]
-        if id > 0 then
-            id = (id + self.segs[i].id_delta) % 0x10000
-        end
-        return id
-    end
-end
-
-function Face:glyph(id)
-    local suboffset
-    if self.index_to_loc_fmt == 0 then      -- short offsets
-        self:goto("loca", 2*id)
-        suboffset = self:uint16() * 2
-    else                                    -- long offsets
-        self:goto("loca", 4*id)
-        suboffset = self:uint32()
-    end
-    self:goto("glyf", suboffset)
-    local num_contours = self:int16()
-    local xmin = self:int16()
-    local ymin = self:int16()
-    local xmax = self:int16()
-    local ymax = self:int16()
-    local points, end_points = {}, {}
-    if num_contours > 0 then        -- simple glyph
-        for i = 1, num_contours do
-            end_points[i] = self:uint16() + 1
-        end
-        local num_points = end_points[#end_points]
-        local instruction_length = self:uint16()
-        local instructions = self:str(instruction_length)
-        local i = 0
-        while i < num_points do
-            i = i + 1
-            local flags = self:uint8()
-            assert(flags < 64, "point flag with higher bits set")
-            local point = {
-                on_curve    = band(flags, 2^0) > 0,
-                x_short     = band(flags, 2^1) > 0,
-                y_short     = band(flags, 2^2) > 0,
-                repeated    = band(flags, 2^3) > 0,
-                x_sign_same = band(flags, 2^4) > 0,
-                y_sign_same = band(flags, 2^5) > 0
-            }
-            points[i] = point
-            if point.repeated then
-                local repeats = self:uint8()
-                for j = 1, repeats do
-                    i = i + 1
-                    points[i] = {
-                        on_curve    = point.on_curve,
-                        x_short     = point.x_short,
-                        y_short     = point.y_short,
-                        x_sign_same = point.x_sign_same,
-                        y_sign_same = point.y_sign_same
-                    }
-                end
-            end
-        end
-        local last_x, last_y = 0, 0
-        for i = 1, #points do
-            if points[i].x_short then
-                local x = self:uint8()
-                if not points[i].x_sign_same then x = -x end
-                points[i].x = last_x + x
-            else
-                if not points[i].x_sign_same then
-                    points[i].x = last_x + self:int16()
-                else
-                    points[i].x = last_x
-                end
-            end
-            last_x = points[i].x
-        end
-        for i = 1, #points do
-            if points[i].y_short then
-                local y = self:uint8()
-                if not points[i].y_sign_same then y = -y end
-                points[i].y = last_y + y
-            else
-                if not points[i].y_sign_same then
-                    points[i].y = last_y + self:int16()
-                else
-                    points[i].y = last_y
-                end
-            end
-            last_y = points[i].y
-        end
-    elseif num_contours < 0 then    -- compound glyph
-        local more = true
-        while more do
-            local flags = self:uint16()
-            local args_are_words    = band(flags, 2^0x0) > 0
-            local args_are_xy       = band(flags, 2^0x1) > 0
-            local round_xy_to_grid  = band(flags, 2^0x2) > 0
-            local regular_scale     = band(flags, 2^0x3) > 0
-            local obsolete          = band(flags, 2^0x4) > 0
-            more                    = band(flags, 2^0x5) > 0
-            local irregular_scale   = band(flags, 2^0x6) > 0
-            local two_by_two        = band(flags, 2^0x7) > 0
-            local instructions      = band(flags, 2^0x8) > 0
-            local use_my_metrics    = band(flags, 2^0x9) > 0
-            local overlap           = band(flags, 2^0xA) > 0
-            local scaled_offset     = band(flags, 2^0xB) > 0
-            local unscaled_offset   = band(flags, 2^0xC) > 0
-            if obsolete then
-                self:log("warning: glyph component using obsolete flag")
-            end
-            local gid = self:uint16()
-            local pos = self:getpos()
-            local sub_points, sub_end_points = self:glyph(gid)
-            self:setpos(pos)
-            local e, f
-            if args_are_xy then
-                if args_are_words then
-                    e = self:int16()
-                    f = self:int16()
-                else
-                    e = self:int8()
-                    f = self:int8()
-                end
-                if round_xy_to_grid then
-                    self:log("warning: ignoring request to round component offset")
-                end
-            else
-                local i, j
-                if args_are_words then
-                    i = self:uint16()
-                    j = self:uint16()
-                else
-                    i = self:uint8()
-                    j = self:uint8()
-                end
-                e = points[i].x - sub_points[j].x
-                f = points[i].y - sub_points[j].y
-            end
-            local a, b, c, d = 1, 0, 0, 1
-            if regular_scale then
-                self:log("regular scale")
-                a = self:int16() / 0x4000
-                d = a
-            elseif irregular_scale then
-                self:log("irregular scale")
-                a = self:int16() / 0x4000
-                d = self:int16() / 0x4000
-            elseif two_by_two then
-                self:log("2x2 transformation")
-                a = self:int16() / 0x4000
-                b = self:int16() / 0x4000
-                c = self:int16() / 0x4000
-                d = self:int16() / 0x4000
-            end
-            local m = math.max(math.abs(a), math.abs(b))
-            local n = math.max(math.abs(c), math.abs(d))
-            if math.abs(math.abs(a)-math.abs(c)) <= 0x21/0x10000 then
-                m = m * 2
-            end
-            if math.abs(math.abs(c)-math.abs(d)) <= 0x21/0x10000 then
-                n = n * 2
-            end
-            for i, p in ipairs(sub_points) do
-                points[#points+1] = {
-                    x=m*(a*p.x/m + c*p.y/m + e),
-                    y=n*(b*p.x/n + d*p.y/n + f),
-                    on_curve=p.on_curve
-                }
-            end
-            local offset = end_points[#end_points] or 0
-            for i, e in ipairs(sub_end_points) do
-                end_points[#end_points+1] = offset + e
-            end
-        end
-        -- TODO: read instructions for composite character
-    end
-    return points, end_points
-end
-
-function Face:pack_outline(points, end_points)
-    local outline = {}
-    local h = self.cap_height
-    local p, q
-    local j = 1
-    for i = 1, #end_points do
-        local contour = {}
-        while j <= end_points[i] do
-            p = points[j]
-            q = {p.x, h-p.y, p.on_curve}
-            table.insert(contour, q)
-            j = j + 1
-        end
-        table.insert(outline, contour)
-    end
-    return outline
-end
-
-function Face:string(s, pt, x, y, anchor, a)
-    anchor = anchor or "tl"
-    a = a or 0
-    local codes = utf8to32(s)
-    local cur_x = 0
-    local contours = {}
-    local outline
-    local li, ri
-    local advance, bearing
-    for i, code in ipairs(codes) do
-        ri = self:char_index(code)
-        if i > 1 and self.num_kernings > 0 then
-            cur_x = cur_x + self:get_kerning(li, ri)
-        end
-        if code ~= 32 then
-            outline = self:pack_outline(self:glyph(ri))
-            for j, contour in ipairs(outline) do
-                for k, point in ipairs(contour) do
-                    point[1] = point[1] + cur_x
-                end
-                table.insert(contours, contour)
-            end
-        end
-        advance, bearing = self:hmetrics(ri)
-        cur_x = cur_x + advance
-        li = ri
-    end
-    local ax, ay -- anchor position
-    local av, ah = anchor:sub(1, 1), anchor:sub(2, 2)
-    if av == "t" then
-        ay = 0
-    elseif av == "m" then
-        ay = self.cap_height / 2
-    elseif av == "b" then
-        ay = self.cap_height
-    end
-    if ah == "l" then
-        ax = 0
-    elseif ah == "c" then
-        ax = cur_x / 2
-    elseif ah == "r" then
-        ax = cur_x
-    end
-    if ax == nil or ay == nil then
-        error("invalid anchor: "..anchor)
-    end
-    local scl = pt * self.resolution / (72 * self.units_per_em)
-    local t = aff.new_affine()
-    t:add_translate(-ax, -ay)
-    t:add_scale(scl)
-    t:add_rotate(a)
-    t:add_translate(x, y)
-    local polygons = {}
-    local polygon
-    for i, contour in ipairs(contours) do
-        if #contour > 2 then
-            t:apply(contour)
-            polygon = poly.unfold(contour)
-            table.insert(polygon, polygon[1]) -- close
-            table.insert(polygons, polygon)
-        end
-    end
-    return polygons
-end
-
-local function load_face(f, dbg)
-    if type(f) == "string" then f = io.open(f, "rb") end
-    local self = setmetatable({fp=f}, Face)
-    self.dbg = dbg or false
-    self:offset()
-    self:directory()
-    self:head()
-    self:maxp()
-    self:cmap()
-    if self.dir["hhea"] then
-        self:hhea()
-    else
-        self:log("no horizontal metrics (hhea+hmtx)")
-    end
-    if self.dir["kern"] then
-        self:kern()
-    else
-        self.num_kernings = 0
-        self:log("no kerning table (kern)")
-    end
-    if self.dir["OS/2"] then
-        self:os2()
-    else
-        self:log("no x-height and Cap-Height (OS/2)")
-    end
-    self.resolution = 300 -- dpi
-    return self
-end
-
-return {load_face=load_face}

diff --git a/util.lua b/util.lua
deleted file mode 100644
index 63f141f..0000000
--- a/util.lua
+++ /dev/null
@@ -1,45 +0,0 @@
-local ffi = require "ffi"
-
-ffi.cdef[[
-double hypot(double x, double y);
-double copysign(double x, double y);
-]]
-local hypot = ffi.C.hypot
-local copysign = ffi.C.copysign
-
-local function round(x)
-    local i, f = math.modf(x + copysign(0.5, x))
-    return i
-end
-
-local function rtrim(s, c)
-    local i = #s
-    while s:sub(i, i) == c do
-        i = i - 1
-    end
-    return s:sub(1, i)
-end
-
-local function startswith(s1, s2)
-    return s1:sub(1, #s2) == s2
-end
-
-local function func_iter(seq)
-    if type(seq) == "table" then
-        local i = 0
-        return function()
-            i = i + 1
-            if i <= #seq then
-                return seq[i]
-            end
-        end
-    else
-        return seq
-    end
-end
-
-return {
-    hypot=hypot, copysign=copysign, round=round,
-    rtrim=rtrim, startswith=startswith,
-    func_iter=func_iter
-}