web/www/study/study.py

#!/usr/bin/env python

import math, gd, random, sys, os

# Select which chapters to run
def chapter(n):
    if len(sys.argv) == 1:
        return True
    return str(n) in sys.argv

##############################################################################

# Tiny image class to make examples short and readable
class Image(gd.image):
    gd.gdMaxColors = 256 * 256 * 256
    def __init__(self, *args):
        if args[0].__class__ == str:
            print "[LOAD] %s" % (args[0],)
        gd.image.__init__(self, *args)
    def save(self, name):
        print "[PNG] %s" % (name,)
        self.writePng(name)
    def getGray(self, x, y):
        p = self.getPixel((x, y))
        c = self.colorComponents(p)[0] / 255.0
        return c
    def getRgb(self, x, y):
        p = self.getPixel((x, y))
        rgb = self.colorComponents(p)
        return [rgb[0] / 255.0, rgb[1] / 255.0, rgb[2] / 255.0]
    def setGray(self, x, y, t):
        p = (int)(t * 255.999)
        c = self.colorResolve((p, p, p))
        self.setPixel((x, y), c)
    def setRgb(self, x, y, r, g, b):
        r = (int)(r * 255.999)
        g = (int)(g * 255.999)
        b = (int)(b * 255.999)
        c = self.colorResolve((r, g, b))
        self.setPixel((x, y), c)
    def getRegion(self, x, y, w, h):
        dest = Image((w, h), True)
        self.copyTo(dest, (-x, -y))
        return dest
    def getZoom(self, z):
        (w, h) = self.size()
        dest = Image((w * z, h * z), True)
        for y in range(h):
            for x in range(w):
                rgb = self.getRgb(x, y)
                for j in range(z):
                    for i in range(z):
                        dest.setRgb(x * z + i, y * z + j, *rgb)
        return dest

# Manipulate gamma values
class Gamma:
    def CtoI(x):
        if x < 0:
            return - math.pow(-x, 2.2)
        return math.pow(x, 2.2)
    def ItoC(x):
        if x < 0:
            return - math.pow(-x, 1 / 2.2)
        return math.pow(x, 1 / 2.2)
    CtoI = staticmethod(CtoI)
    ItoC = staticmethod(ItoC)
    def CtoI3(x):
        return [Gamma.CtoI(x[0]), Gamma.CtoI(x[1]), Gamma.CtoI(x[2])]
    def ItoC3(x):
        return [Gamma.ItoC(x[0]), Gamma.ItoC(x[1]), Gamma.ItoC(x[2])]
    CtoI3 = staticmethod(CtoI3)
    ItoC3 = staticmethod(ItoC3)
    def Cto2(x):
        if x < Gamma.CtoI(0.50):
            return 0.
        return 1.
    def Cto3(x):
        if x < Gamma.CtoI(0.25):
            return 0.
        elif x < Gamma.CtoI(0.75):
            return Gamma.CtoI(0.5)
        return 1.
    def Cto4(x):
        if x < Gamma.CtoI(0.17):
            return 0.
        elif x < Gamma.CtoI(0.50):
            return Gamma.CtoI(0.3333)
        elif x < Gamma.CtoI(0.83):
            return Gamma.CtoI(0.6666)
        return 1.
    Cto2 = staticmethod(Cto2)
    Cto3 = staticmethod(Cto3)
    Cto4 = staticmethod(Cto4)

# Create matrices
def Matrix(w, h, val = 0):
    return [[val] * w for n in range(h)]

# Iterate in 2D space
def rangexy(w, h):
    for y in range(h):
        for x in range(w):
            yield (x, y)

##############################################################################

# Create SVG files from matrix data
class Svg:
    _data = ''
    _w, _h = 0, 0
    def _colorise(val):
        return '#fff'
    def _reduce(this, val):
        if type(val) == float:
            for x in range(1, 1000):
                if abs(val * x - round(val * x)) < 0.001:
                    return (int)(round(val * x)), x
        return val, 1
    def __init__(self, mat, colorise = _colorise):
        # Check whether it is an error diffusion matrix
        ed = False
        for l in mat:
            for x in l:
                if type(x) == str or math.floor(x) != x:
                    ed = True
        # Generate SVG file
        (w, h) = (len(mat[0]), len(mat))
        s = \
         '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n' \
         '<svg\n' \
         ' xmlns="http://www.w3.org/2000/svg"\n' \
         ' xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"\n' \
         ' xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"\n' \
         ' width="%u"\n' \
         ' height="%u">\n' \
         '  <sodipodi:namedview inkscape:document-units="px" />\n' \
         '  <g>\n' % (64 * w + 2, 64 * h + 2)
        line = \
         '  <path sodipodi:nodetypes="cc" d="M %u,%u L %u,%u" ' \
               'style="stroke:#000;stroke-width:1" />\n'
        box = \
         '  <path sodipodi:nodetypes="cccc" ' \
               'd="M %u,%u L %u,%u L %u,%u L %u,%u z" ' \
               'style="fill:%s;stroke:#000;stroke-width:2" />\n'
        text = \
         '  <text style="font-size:%upx;text-align:center;text-anchor:middle;' \
               'font-family:Bitstream Vera Serif" x="%u" y="%u"%s>%s</text>\n'
        for x, y in rangexy(w, h):
            val = mat[y][x]
            if not ed and val == -1:
                continue
            if ed and val == 0:
                continue
            # Put box
            (ix, iy) = (64. * x + 1, 64. * y + 1)
            if ed and val == -1:
                val = ''
                c = '#fbb'
            else:
                c = colorise(val)
            s += box % (ix, iy, ix + 64, iy, ix + 64, iy + 64, ix, \
                        iy + 64, c)
            # Put value
            a, b = self._reduce(val)
            extra = ''
            if b == 1:
                (tx, ty) = (ix + 32, iy + 44)
                n = len(str(val))
                if n > 3 and type(val) == str:
                    extra = ' transform="scale(0.7,0.7)"'
                    tx = tx / 0.7
                    ty = ty / 0.7
                elif n > 3:
                    extra = ' transform="scale(0.7,1.)"'
                    tx = tx / 0.7
                elif n > 2:
                    extra = ' transform="scale(0.8,1.)"'
                    tx = tx / 0.8
                s += text % (32, tx, ty, extra, val)
            else:
                s += line % (ix + 8, iy + 32, ix + 56, iy + 32)
                (tx, ty) = (ix + 32, iy + 26)
                s += text % (24, tx, ty, extra, a)
                (tx, ty) = (ix + 32, iy + 56)
                s += text % (24, tx, ty, extra, b)
        s += \
         '  </g>\n' \
         '</svg>\n'
        self._w = 64 * w + 2
        self._h = 64 * h + 2
        self._data = s
    def save(self, name, size):
        svgname = name + ".tmp.svg"
        f = open(svgname, 'w')
        f.write(self._data)
        f.close()
        f = os.popen("inkscape %s -a %u:%u:%u:%u -w%u -h%u -e %s >/dev/null 2>&1" % (svgname, 0, 0, self._w, self._h, size * self._w / 64., size * self._h / 64., name))
        print "[SVG] %s" % (name,)
        f.close()
        os.unlink(svgname)

##############################################################################
print "Initialisation"

# Load the original Lena image
lena512 = Image("lena512.png")
(w, h) = lena512.size()

# Image 1: greyscale conversion
# Read the compression FAQ [55] for the rationale behind using the green
# channel (http://www.faqs.org/faqs/compression-faq/part1/section-30.html)
if chapter(0):
    (w, h) = lena512.size()
    lena512bw = Image((w, h))
    for x, y in rangexy(w, h):
        rgb = lena512.getRgb(x, y)
        c = rgb[1]
        lena512bw.setGray(x, y, c)
    lena512bw.save("lena512bw.png")
else:
    lena512bw = Image("lena512bw.png")

def gammascale(src, scale):
    (w, h) = src.size()
    count = src.colorsTotal()
    dest = Image((w / scale, h / scale), count == 0 or count > 256)
    for x, y in rangexy(w / scale, h / scale):
        r = g = b = 0.
        for i, j in rangexy(scale, scale):
            rgb = src.getRgb(x * scale + i, y * scale + j)
            r += Gamma.CtoI(rgb[0])
            g += Gamma.CtoI(rgb[1])
            b += Gamma.CtoI(rgb[2])
        r = Gamma.ItoC(r / (scale * scale))
        g = Gamma.ItoC(g / (scale * scale))
        b = Gamma.ItoC(b / (scale * scale))
        dest.setRgb(x, y, r, g, b)
    return dest

# Image 2: 50% greyscale
# Image 3: 50% scaling
if chapter(0):
    lena256bw = gammascale(lena512bw, 2)
    lena256bw.save("lena256bw.png")
    lena256 = gammascale(lena512, 2)
    lena256.save("lena256.png")
else:
    lena256bw = Image("lena256bw.png")
    lena256 = Image("lena256.png")

# Create a 32x256 greyscale gradient
if chapter(0):
    grad256bw = Image((32, 256))
    for x, y in rangexy(32, 256):
        grad256bw.setGray(x, 255 - y, y / 255.)
    grad256bw.save("gradient256bw.png")
else:
    grad256bw = Image("gradient256bw.png")

# Create a 64x256 colour gradient
if chapter(0):
    grad256 = Image((64, 256), True)
    for x, y in rangexy(64, 256):
        grad256.setRgb(x, y, x / 63., (255. - y) / 255, x / 63.)
    grad256.save("gradient256.png")
else:
    grad256 = Image("gradient256.png")

##############################################################################
if chapter(1):
    print "Chapter 1. Colour quantisation"

# Output 1.1.1: 50% threshold
# Output 1.1.2: 40% threshold
# Output 1.1.3: 60% threshold
def test11x(src, threshold):
    (w, h) = src.size()
    dest = Image((w, h))
    for x, y in rangexy(w, h):
        c = src.getGray(x, y) > threshold
        dest.setGray(x, y, c)
    return dest

if chapter(1):
    test11x(grad256bw, 0.5).save("out/grad1-1-1.png")
    test11x(lena256bw, 0.5).save("out/lena1-1-1.png")
    test11x(grad256bw, 0.4).save("out/grad1-1-2.png")
    test11x(lena256bw, 0.4).save("out/lena1-1-2.png")
    test11x(grad256bw, 0.6).save("out/grad1-1-3.png")
    test11x(lena256bw, 0.6).save("out/lena1-1-3.png")

# Output 1.2.1: 3-colour threshold
# Output 1.2.2: 5-colour threshold
def test12x(src, colors):
    (w, h) = src.size()
    dest = Image((w, h))
    q = colors - 1
    p = -.00001 + colors
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        c = math.floor(c * p) / q
        dest.setGray(x, y, c)
    return dest

if chapter(1):
    test12x(grad256bw, 3).save("out/grad1-2-1.png")
    test12x(lena256bw, 3).save("out/lena1-2-1.png")
    test12x(grad256bw, 5).save("out/grad1-2-2.png")
    test12x(lena256bw, 5).save("out/lena1-2-2.png")

# Output 1.2.3: 3-colour threshold, minimal error
# Output 1.2.4: 5-colour threshold, minimal error
def test12y(src, colors):
    (w, h) = src.size()
    dest = Image((w, h))
    q = colors - 1
    p = -.00001 + colors - 1
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        c = math.floor((c + 0.5 / p) * p) / q
        dest.setGray(x, y, c)
    return dest

if chapter(1):
    test12y(grad256bw, 3).save("out/grad1-2-3.png")
    test12y(lena256bw, 3).save("out/lena1-2-3.png")
    test12y(grad256bw, 5).save("out/grad1-2-4.png")
    test12y(lena256bw, 5).save("out/lena1-2-4.png")

# Output 1.3.1: 2-colour threshold, dynamic thresholding
# Output 1.3.2: 5-colour threshold, dynamic thresholding
def test13x(src, n, fuzzy = None):
    random.seed(0)
    (w, h) = src.size()
    dest = Image((w, h))
    # Compute histogram
    histo = [0] * 256
    for x, y in rangexy(w, h):
        histo[(int)(src.getGray(x, y) * 255.9999)] += 1
    thresholds = [1. * (1. + i) / n for i in range(n - 1)]
    values = [i / (n - 1.) for i in range(n)]
    # Parse histogram
    total = 0
    t = 0
    for i in range(256):
        total += histo[i]
        if total > thresholds[t] * w * h:
            thresholds[t] = i / 255.0
            t += 1
            if t + 1 > n - 1:
                break
    # Compute image
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        for (i, t) in enumerate(thresholds):
            if fuzzy:
                t += fuzzy()
            if c < t:
                dest.setGray(x, y, values[i])
                break
        else:
            dest.setGray(x, y, values[n - 1])
    return dest

if chapter(1):
    test13x(grad256bw, 2).save("out/grad1-3-1.png")
    test13x(lena256bw, 2).save("out/lena1-3-1.png")
    test13x(grad256bw, 5).save("out/grad1-3-2.png")
    test13x(lena256bw, 5).save("out/lena1-3-2.png")

# Output 1.4.1: uniform random dithering
def test141(src):
    random.seed(0)
    (w, h) = src.size()
    dest = Image((w, h))
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        d = c > random.random()
        dest.setGray(x, y, d)
    return dest

if chapter(1):
    test141(grad256bw).save("out/grad1-4-1.png")
    test141(lena256bw).save("out/lena1-4-1.png")

# Output 1.4.2: gaussian random dithering
def test142(src):
    random.seed(0)
    (w, h) = src.size()
    dest = Image((w, h))
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        d = c > random.gauss(0.5, 0.15)
        dest.setGray(x, y, d)
    return dest

if chapter(1):
    test142(grad256bw).save("out/grad1-4-2.png")
    test142(lena256bw).save("out/lena1-4-2.png")

# Output 1.4.3: 3-colour threshold, dynamic thresholding, gaussian perturbation
if chapter(1):
    fuzzy = lambda : random.gauss(0., 0.08)
    test13x(grad256bw, 4, fuzzy).save("out/grad1-4-3.png")
    test13x(lena256bw, 4, fuzzy).save("out/lena1-4-3.png")

##############################################################################
if chapter(2):
    print "Chapter 2. Halftoning patterns"

# Pattern 2.1.1: a 50% halftone pattern with various block sizes
# Pattern 2.1.2: 25% and 75% halftone patterns with various block sizes
if chapter(2):
    dest = Image((320, 80))
    for x in range(320):
        d = 8 >> (x / 80)
        for y in range(80):
            c = (x / d + y / d) & 1
            dest.setGray(x, y, c)
    dest.save("out/pat2-1-1.png")

    dest = Image((320, 80))
    for x in range(320):
        d = 8 >> (x / 80)
        for y in range(40):
            c = ((x / d + y / d) & 1) or (y / d & 1)
            dest.setGray(x, y, c)
        for y in range(40, 80):
            c = ((x / d + y / d) & 1) and (y / d & 1)
            dest.setGray(x, y, c)
    dest.save("out/pat2-1-2.png")

# Output 2.1.1: 20/40/60/80% threshold with 25/50/75% patterns inbetween:
def test211(src):
    (w, h) = src.size()
    dest = Image((w, h))
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        if c < 0.2:
            c = 0.
        elif c < 0.4:
            c = ((x + y) & 1) and (y & 1)
        elif c < 0.6:
            c = (x + y) & 1
        elif c < 0.8:
            c = ((x + y) & 1) or (y & 1)
        else:
            c = 1.
        dest.setGray(x, y, c)
    return dest

if chapter(2):
    test211(grad256bw).save("out/grad2-1-1.png")
    test211(lena256bw).save("out/lena2-1-1.png")

# Pattern 2.2.1: vertical, mixed and horizontal black-white halftones
# Pattern 2.2.2: two different 25% patterns
if chapter(2):
    dest = Image((240, 80))
    for y in range(80):
        for x in range(80):
            c = x & 1
            dest.setGray(x, y, c)
        for x in range(80, 160):
            c = (x / d + y / d) & 1
            dest.setGray(x, y, c)
        for x in range(160, 240):
            c = y & 1
            dest.setGray(x, y, c)
    dest.save("out/pat2-2-1.png")

    dest = Image((320, 80))
    for y in range(80):
        for x in range(80):
            c = (x / 2 & 1) and (y / 2 & 1)
            dest.setGray(x, y, c)
        for x in range(80, 160):
            c = (x & 1) and (y & 1)
            dest.setGray(x, y, c)
        for x in range(160, 240):
            c = (x & 1) and ((y + x / 2) & 1)
            dest.setGray(x, y, c)
        for x in range(240, 320):
            c = (x / 2 & 1) and ((y / 2 + x / 4) & 1)
            dest.setGray(x, y, c)
    dest.save("out/pat2-2-2.png")

# Output 2.3.0: 2x2 Bayer dithering
# Output 2.3.1: 4x4 Bayer dithering
# Output 2.3.1b: 8x8 Bayer dithering
# Output 2.3.2: 4x4 cluster dot
# Output 2.3.2b: 8x8 cluster dot
# Output 2.3.3: 5x3 line dithering
def ordereddither(src, mat):
    (w, h) = src.size()
    dest = Image((w, h))
    dx = len(mat[0])
    dy = len(mat)
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        threshold = (1. + mat[y % dy][x % dx]) / (dx * dy + 1)
        c = c > threshold
        dest.setGray(x, y, c)
    return dest

def makebayer(rank, mat = False):
    if not mat:
        mat = Matrix(1, 1)
    if not rank:
        return mat
    n = len(mat)
    newmat = Matrix(n * 2, n * 2)
    for i, j in rangexy(n, n):
        x = mat[j][i]
        newmat[j * 2][i * 2] = x
        newmat[j * 2][i * 2 + 1] = x + n * n * 3
        newmat[j * 2 + 1][i * 2] = x + n * n * 2
        newmat[j * 2 + 1][i * 2 + 1] = x + n * n
    return makebayer(rank - 1, newmat)

DITHER_BAYER22 = makebayer(1)
DITHER_BAYER44 = makebayer(2)
DITHER_BAYER88 = makebayer(3)

DITHER_CLUSTER44 = \
    [[ 12,  5,  6, 13],
     [  4,  0,  1,  7],
     [ 11,  3,  2,  8],
     [ 15, 10,  9, 14]]
DITHER_CLUSTER88 = \
    [[ 24, 10, 12, 26, 35, 47, 49, 37],
     [  8,  0,  2, 14, 45, 59, 61, 51],
     [ 22,  6,  4, 16, 43, 57, 63, 53],
     [ 30, 20, 18, 28, 33, 41, 55, 39],
     [ 34, 46, 48, 36, 25, 11, 13, 27],
     [ 44, 58, 60, 50,  9,  1,  3, 15],
     [ 42, 56, 62, 52, 23,  7,  5, 17],
     [ 32, 40, 54, 38, 31, 21, 19, 29]]
DITHER_LINE53 = \
    [[  9,  3,  0,  6, 12],
     [ 10,  4,  1,  7, 13],
     [ 11,  5,  2,  8, 14]]

def bayercol(val): return ['#fdf', '#dfd', '#ffd', '#dff'][val % 4]

if chapter(2):
    ordereddither(grad256bw, DITHER_BAYER22).save("out/grad2-3-0.png")
    ordereddither(lena256bw, DITHER_BAYER22).save("out/lena2-3-0.png")

    Svg(DITHER_BAYER44, bayercol).save("out/fig2-3-2.png", 40)
    ordereddither(grad256bw, DITHER_BAYER44).save("out/grad2-3-1.png")
    ordereddither(lena256bw, DITHER_BAYER44).save("out/lena2-3-1.png")

    Svg(DITHER_BAYER88, bayercol).save("out/fig2-3-2b.png", 30)
    ordereddither(grad256bw, DITHER_BAYER88).save("out/grad2-3-1b.png")
    ordereddither(lena256bw, DITHER_BAYER88).save("out/lena2-3-1b.png")

    Svg(DITHER_CLUSTER44).save("out/fig2-3-3.png", 40)
    ordereddither(grad256bw, DITHER_CLUSTER44).save("out/grad2-3-2.png")
    ordereddither(lena256bw, DITHER_CLUSTER44).save("out/lena2-3-2.png")

    Svg(DITHER_CLUSTER88).save("out/fig2-3-3b.png", 30)
    ordereddither(grad256bw, DITHER_CLUSTER88).save("out/grad2-3-2b.png")
    ordereddither(lena256bw, DITHER_CLUSTER88).save("out/lena2-3-2b.png")

    Svg(DITHER_LINE53).save("out/fig2-3-4.png", 40)
    ordereddither(grad256bw, DITHER_LINE53).save("out/grad2-3-3.png")
    ordereddither(lena256bw, DITHER_LINE53).save("out/lena2-3-3.png")

# Output 2.4.1: 4x4 Bayer dithering with gaussian perturbation
def test241(src, mat):
    random.seed(0)
    (w, h) = src.size()
    dest = Image((w, h))
    dx = len(mat[0])
    dy = len(mat)
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        threshold = (1. + mat[y % dy][x % dx]) / (dx * dy + 1)
        threshold += random.gauss(0, 0.08)
        c = c > threshold
        dest.setGray(x, y, c)
    return dest

if chapter(2):
    test241(grad256bw, DITHER_BAYER88).save("out/grad2-4-1.png")
    test241(lena256bw, DITHER_BAYER88).save("out/lena2-4-1.png")

# Output 2.4.2: random dither matrice selection
def test242(src, mlist):
    random.seed(0)
    (w, h) = src.size()
    dest = Image((w, h))
    dx = len(mlist[0][0])
    dy = len(mlist[0])
    for x, y in rangexy(w / dx, h / dy):
        mat = random.choice(mlist)
        for i, j in rangexy(dx, dy):
            c = src.getGray(x * dx + i, y * dy + j)
            threshold = (1. + mat[j][i]) / (dx * dy + 1)
            d = c > threshold
            dest.setGray(x * dx + i, y * dy + j, d)
    return dest

if chapter(2):
    m1 = [[1, 4, 7],
          [6, 0, 2],
          [3, 8, 5]]
    m2 = [[4, 6, 3],
          [8, 1, 5],
          [0, 3, 7]]
    m3 = [[5, 0, 3],
          [2, 8, 6],
          [7, 4, 1]]
    m4 = [[8, 2, 5],
          [6, 4, 0],
          [1, 7, 3]]
    m5 = [[2, 5, 8],
          [0, 7, 3],
          [4, 1, 6]]
    m6 = [[7, 4, 1],
          [3, 6, 8],
          [2, 0, 5]]
    mlist = [m1, m2, m3, m4, m5, m6]
    test242(grad256bw, mlist).save("out/grad2-4-2.png")
    test242(lena256bw, mlist).save("out/lena2-4-2.png")

# Output 2.5.1: cross pattern
# Output 2.5.2: hex pattern
# Output 2.5.3: square pattern
def test25x(src, mat, vec):
    # 1. count positive cells
    n = 0
    for line in mat:
        for x in line:
            if x >= 0:
                n += 1
    # 2. create list of vectors
    l = []
    x = y = 0
    while (x, y) not in l:
        l.append((x, y))
        (x, y) = ((x + vec[0][0]) % n, (y + vec[0][1]) % n)
        if (x, y) in l:
            (x, y) = ((x + vec[1][0]) % n, (y + vec[1][1]) % n)
    # 3. create big matrix
    m = Matrix(n, n)
    for v in l:
        for x, y in rangexy(len(mat[0]), len(mat)):
            if mat[y][x] < 0:
                continue
            m[(v[1] + y + n) % n][(v[0] + x + n) % n] = mat[y][x]
    # 4. dither image
    (w, h) = src.size()
    dest = Image((w, h))
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        threshold = (1. + m[y % n][x % n]) / (1. + n)
        d = c > threshold
        dest.setGray(x, y, d)
    return dest

if chapter(2):
    mat = [[-1,  4, -1],
           [ 3,  0,  1],
           [-1,  2, -1]]
    vec = [(2, -1), (1, 2)]
    test25x(grad256bw, mat, vec).save("out/grad2-5-1.png")
    test25x(lena256bw, mat, vec).save("out/lena2-5-1.png")
    mat = [[-1,  4,  2, -1],
           [ 6,  0,  5,  3],
           [-1,  7,  1, -1]]
    vec = [(2, -2), (3, 1)]
    test25x(grad256bw, mat, vec).save("out/grad2-5-2.png")
    test25x(lena256bw, mat, vec).save("out/lena2-5-2.png")
    mat = [[-1, -1, -1,  7, -1],
           [-1,  2,  6,  9,  8],
           [ 3,  0,  1,  5, -1],
           [-1,  4, -1, -1, -1]]
    vec = [(2, 4), (3, 1)]
    test25x(grad256bw, mat, vec).save("out/grad2-5-3.png")
    test25x(lena256bw, mat, vec).save("out/lena2-5-3.png")
    mat = [[-1, -1,  2, -1],
           [ 0,  1,  3,  8],
           [ 5,  4,  7,  6],
           [-1,  9, -1, -1]]
    vec = [(2, 2), (0, 5)]
    test25x(grad256bw, mat, vec).save("out/grad2-5-4.png")
    test25x(lena256bw, mat, vec).save("out/lena2-5-4.png")

# Output 2.6.1: 4-wise cross pattern
# Output 2.6.2: 3-wise hex pattern
# Output 2.6.3: 4-wise square pattern
if chapter(2):
    mat = [[-1, -1, -1, 18, -1, -1],
           [-1, 16, 14,  2,  6, -1],
           [12,  0,  4, 10, 17, -1],
           [-1,  8, 19, 13,  1,  5],
           [-1, 15,  3,  7,  9, -1],
           [-1, -1, 11, -1, -1, -1]]
    vec = [(4, -2), (2, 4)]
    test25x(grad256bw, mat, vec).save("out/grad2-6-1.png")
    test25x(lena256bw, mat, vec).save("out/lena2-6-1.png")

    mat = [[-1, 12,  6, -1, -1, -1, -1],
           [18,  0, 15, 9,  14,  8, -1],
           [-1, 21,  3, 20,  2, 17, 11],
           [-1, -1, 13,  7, 23,  5, -1],
           [-1, 19,  1, 16, 10, -1, -1],
           [-1, -1, 22,  4, -1, -1, -1]]
    vec = [(5, -1), (-1, 5)]
    test25x(grad256bw, mat, vec).save("out/grad2-6-2.png")
    test25x(lena256bw, mat, vec).save("out/lena2-6-2.png")

    mat = [[-1, -1, -1, -1, 30, -1, -1, -1, -1],
           [-1, -1, 10, 26, 38, 34, -1, 29, -1],
           [-1, 14,  2,  6, 22,  9, 25, 37, 33],
           [-1, -1, 18, 28, 13,  1,  5, 21, -1],
           [-1,  8, 24, 36, 32, 17, 31, -1, -1],
           [12,  0,  4, 20, 11, 27, 39, 35, -1],
           [-1, 16, -1, 15,  3,  7, 23, -1, -1],
           [-1, -1, -1, -1, 19, -1, -1, -1, -1]]
    vec = [(6, 2), (-2, 6)]
    test25x(grad256bw, mat, vec).save("out/grad2-6-3.png")
    test25x(lena256bw, mat, vec).save("out/lena2-6-3.png")

    mat = [[-1, -1,  3, 35, -1, -1,  1, 33, -1, -1, -1, -1, -1],
           [-1, -1, -1, 19, 51, -1, -1, 17, 49, -1, -1, -1, -1],
           [ 7, 39, 11, 43,  5, 37,  9, 41, -1, -1, -1, -1, -1],
           [-1, 23, 55, 27, 59, 21, 53, 25, 57, -1, -1, -1, -1],
           [15, 47, -1, -1, 13, 45,  2, 34, -1, -1,  0, 32, -1],
           [-1, 31, 63, -1, -1, 29, 61, 18, 50, -1, -1, 16, 48],
           [-1, -1, -1, -1,  6, 38, 10, 42,  4, 36,  8, 40, -1],
           [-1, -1, -1, -1, -1, 22, 54, 26, 58, 20, 52, 24, 56],
           [-1, -1, -1, -1, 14, 46, -1, -1, 12, 44, -1, -1, -1],
           [-1, -1, -1, -1, -1, 30, 62, -1, -1, 28, 60, -1, -1]]
    vec = [(8, 0), (0, 8)]
    def colorise(val):
        if val == 0: return '#fff'
        if (val % 64) == 32: return '#ddf'
        if (val % 32) == 16: return '#fdd'
        if (val % 16) == 8: return '#dff'
        if (val % 8) == 4: return '#ffd'
        if (val % 4) == 2: return '#fdf'
        return '#dfd'
    Svg(mat, colorise).save("out/fig2-6-4.png", 25)
    test25x(grad256bw, mat, vec).save("out/grad2-6-4.png")
    test25x(lena256bw, mat, vec).save("out/lena2-6-4.png")

    mat = [[ -1, -1, -1, -1, -1, -1,  8, -1],
           [ -1, -1,  6, -1,  2,  5, 11, 26],
           [  0,  3,  9, 24, 17, 14, 23, 20],
           [ 15, 12, 21, 18,  7, 29, -1, -1],
           [ -1, 27,  1,  4, 10, 25, -1, -1],
           [ -1, -1, 16, 13, 22, 19, -1, -1],
           [ -1, -1, -1, 28, -1, -1, -1, -1]]
    vec = [(6, 1), (0, 5)]
    def colorise(val): return ['#dff', '#ffd', '#fdf'][val % 3]
    #Svg(mat, colorise).save("out/fig2-6-5.png", 30)
    test25x(grad256bw, mat, vec).save("out/grad2-6-5.png")
    test25x(lena256bw, mat, vec).save("out/lena2-6-5.png")

    mat = [[ -1, -1, -1, -1, -1, -1, 24, -1, -1, -1, -1, -1, -1, -1],
           [ -1, -1, 18, -1,  6, 15, 33, 78, -1, -1, -1, -1, 26, -1],
           [  0,  9, 27, 72, 51, 42, 69, 60, 20, -1,  8, 17, 35, 80],
           [ 45, 36, 63, 54, 21, 87,  2, 11, 29, 74, 53, 44, 71, 62],
           [ -1, 81,  3, 12, 30, 75, 47, 38, 65, 56, 23, 89, -1, -1],
           [ -1, -1, 48, 39, 66, 57, 25, 83,  5, 14, 32, 77, -1, -1],
           [ -1, -1, 19, 84,  7, 16, 34, 79, 50, 41, 68, 59, -1, -1],
           [  1, 10, 28, 73, 52, 43, 70, 61, -1, 86, -1, -1, -1, -1],
           [ 46, 37, 64, 55, 22, 88, -1, -1, -1, -1, -1, -1, -1, -1],
           [ -1, 82,  4, 13, 31, 76, -1, -1, -1, -1, -1, -1, -1, -1],
           [ -1, -1, 49, 40, 67, 58, -1, -1, -1, -1, -1, -1, -1, -1],
           [ -1, -1, -1, 85, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]]
    test25x(grad256bw, mat, vec).save("out/grad2-6-6.png")
    test25x(lena256bw, mat, vec).save("out/lena2-6-6.png")

# Output 2.7.1: void and cluster matrix generation
def makegauss(n):
    c = (-1. + n) / 2
    mat = Matrix(n, n)
    for x, y in rangexy(n, n):
        mat[y][x] = math.exp(- ((c - x) * (c - x) + (c - y) * (c - y)) / (0.05 * n * n))
    return mat

def countones(mat):
    total = 0
    for l in mat:
        for x in l:
            if x:
                total += 1
    return total

GAUSS77 = makegauss(7)
GAUSS99 = makegauss(9)

def getminmax(mat, c):
    min = 9999.
    max = 0.
    h = len(mat)
    w = len(mat[0])
    for x, y in rangexy(w, h):
        if mat[y][x] != c:
            continue
        total = 0.
        for i, j in rangexy(7, 7):
            total += mat[(y + j - 3 + h) % h][(x + i - 3 + w) % w] \
                      * GAUSS77[j][i]
        if total > max:
            (max, max_x, max_y) = (total, x, y)
        if total < min:
            (min, min_x, min_y) = (total, x, y)
    return (min_x, min_y, max_x, max_y)

def makeuniform(n):
    random.seed(0)
    mat = Matrix(n, n)
    for t in range(n * n / 10):
        x = (int)(random.random() * n)
        y = (int)(random.random() * n)
        mat[y][x] = 1
    while True:
        (dummy0, dummy1, x, y) = getminmax(mat, 1.)
        mat[y][x] = 0.
        (x2, y2, dummy0, dummy1) = getminmax(mat, 0.)
        mat[y2][x2] = 1.
        if x2 == x and y2 == y:
            break
    return mat

def makevoidandcluster(n):
    vnc = Matrix(n, n)
    # Step 1: step down to zero
    mat = makeuniform(n)
    rank = countones(mat)
    while rank > 0:
        rank -= 1
        (dummy0, dummy1, x, y) = getminmax(mat, 1.)
        mat[y][x] = 0.
        vnc[y][x] = rank
    # Step 2: step up to n * n
    mat = makeuniform(n)
    rank = countones(mat)
    while rank < n * n:
        (x, y, dummy0, dummy1) = getminmax(mat, 0.)
        mat[y][x] = 1.
        vnc[y][x] = rank
        rank += 1
    return vnc

def vnccol(n):
    return lambda val : ['#fff', '#ddd'][val < n * n / 10]

if chapter(2):
    tmp = makevoidandcluster(14)
    Svg(tmp, vnccol(14)).save("out/fig2-7-1.png", 25)
    ordereddither(grad256bw, tmp).save("out/grad2-7-1.png")
    ordereddither(lena256bw, tmp).save("out/lena2-7-1.png")
    tmp = makevoidandcluster(25)
    Svg(tmp, vnccol(25)).save("out/fig2-7-2.png", 20)
    ordereddither(grad256bw, tmp).save("out/grad2-7-2.png")
    ordereddither(lena256bw, tmp).save("out/lena2-7-2.png")

##############################################################################
if chapter(3):
    print "Chapter 3. Error diffusion"

# Output 3.0.1: naive error diffusion
# Output 3.1.1: standard Floyd-Steinberg
# Output 3.1.2: serpentine Floyd-Steinberg
# FIXME: serpentine only works if rows == offset * 2 + 1
# Output 3.2.1: Fan (modified Floyd-Steinberg)
# Output 3.2.1b: Shiau-Fan 1
# Output 3.2.1c: Shiau-Fan 2
# Output 3.2.2: Jarvis, Judice and Ninke
# Output 3.2.3: Stucki
# Output 3.2.4: Burkes
# Output 3.2.5: Sierra
# Output 3.2.6: Two-line Sierra
# Output 3.2.7: Sierra's Filter Lite
# Output 3.2.8: Atkinson
def test3xx(src, mat, serpentine):
    (w, h) = src.size()
    dest = Image((w, h))
    lines = len(mat)
    rows = len(mat[0])
    offset = mat[0].index(-1)
    ey = Matrix(w + rows - 1, lines, 0.)
    for y in range(h):
        ex = [0.] * (rows - offset)
        if serpentine and y & 1:
            xrange = range(w - 1, -1, -1)
        else:
            xrange = range(w)
        for x in xrange:
            # Set pixel
            c = src.getGray(x, y) + ex[0] + ey[0][x + offset]
            d = c > 0.5
            dest.setGray(x, y, d)
            error = c - d
            # Propagate first line of error
            for dx in range(rows - offset - 2):
                ex[dx] = ex[dx + 1] + error * mat[0][offset + 1 + dx]
            ex[rows - offset - 2] = error * mat[0][rows - 1]
            # Propagate next lines
            if serpentine and y & 1:
                for dy in range(1, lines):
                    for dx in range(rows):
                        ey[dy][x + dx] += error * mat[dy][rows - 1 - dx]
            else:
                for dy in range(1, lines):
                    for dx in range(rows):
                        ey[dy][x + dx] += error * mat[dy][dx]
        for dy in range(lines - 1):
            ey[dy] = ey[dy + 1]
        ey[lines - 1] = [0.] * (w + rows - 1)
    return dest

ERROR_NAIVE = \
    [[ -1, 1]]
ERROR_FSTEIN = \
    [[    0.,    -1, 7./16],
     [ 3./16, 5./16, 1./16]]
ERROR_JAJUNI = \
    [[    0.,    0.,    -1, 7./48, 5./48],
     [ 3./48, 5./48, 7./48, 5./48, 3./48],
     [ 1./48, 3./48, 5./48, 3./48, 1./48]]
ERROR_FAN = \
    [[    0.,    0.,    -1, 7./16],
     [ 1./16, 3./16, 5./16,    0.]]
ERROR_SHIAUFAN = \
    [[    0.,    0.,    -1, 8./16],
     [ 2./16, 2./16, 4./16,    0.]]
ERROR_SHIAUFAN2 = \
    [[    0.,    0.,    0.,    -1, 8./16],
     [ 1./16, 1./16, 2./16, 4./16,    0.]]
ERROR_STUCKI = \
    [[    0.,    0.,    -1, 8./42, 4./42],
     [ 2./42, 4./42, 8./42, 4./42, 2./42],
     [ 1./42, 2./42, 4./42, 2./42, 1./42]]
ERROR_BURKES = \
    [[    0.,    0.,    -1, 8./32, 4./32],
     [ 2./32, 4./32, 8./32, 4./32, 2./32]]
ERROR_SIERRA = \
    [[    0.,    0.,    -1, 5./32, 3./32],
     [ 2./32, 4./32, 5./32, 4./32, 2./32],
     [    0., 2./32, 3./32, 2./32,    0.]]
ERROR_SIERRA2 = \
    [[    0.,    0.,    -1, 4./16, 3./16],
     [ 1./16, 2./16, 3./16, 2./16, 1./16]]
ERROR_FILTERLITE = \
    [[   0.,   -1, 2./4],
     [ 1./4, 1./4,   0.]]
ERROR_ATKINSON = \
    [[   0.,   -1, 1./8, 1./8],
     [ 1./8, 1./8, 1./8,   0.],
     [   0., 1./8,   0.,   0.]]
## This is Stevenson-Arce in hex lattice
#ERROR_STAR = \
#    [[      0.,      0.,      0.,      -1,      0.,  32./200,      0.],
#     [ 12./200,      0., 26./200,      0., 30./200,       0., 16./200],
#     [      0., 12./200,      0., 26./200,      0.,  12./200,      0.],
#     [  5./200,      0., 12./200,      0., 12./200,       0.,  5./200]]
## This is an attempt at implementing Stevenson-Arce in square lattice
#ERROR_STAR = \
#    [[      0.,      0.,      -1, 32./200,      0.],
#     [  6./200, 19./200, 28./200, 23./200,  8./200],
#     [      0., 12./200, 26./200, 12./200,      0.],
#     [  2./200,  9./200, 12./200,  9./200,  2./200]]

if chapter(3):
    test3xx(grad256bw, ERROR_NAIVE, False).save("out/grad3-0-1.png")
    test3xx(lena256bw, ERROR_NAIVE, False).save("out/lena3-0-1.png")

    Svg(ERROR_FSTEIN).save("out/fig3-1-1.png", 40)
    test3xx(grad256bw, ERROR_FSTEIN, False).save("out/grad3-1-1.png")
    test3xx(lena256bw, ERROR_FSTEIN, False).save("out/lena3-1-1.png")
    test3xx(grad256bw, ERROR_FSTEIN, True).save("out/grad3-1-2.png")
    test3xx(lena256bw, ERROR_FSTEIN, True).save("out/lena3-1-2.png")

    Svg(ERROR_JAJUNI).save("out/fig3-1-3.png", 40)
    test3xx(grad256bw, ERROR_JAJUNI, False).save("out/grad3-1-3.png")
    test3xx(lena256bw, ERROR_JAJUNI, False).save("out/lena3-1-3.png")

    Svg(ERROR_FAN).save("out/fig3-2-1.png", 40)
    test3xx(grad256bw, ERROR_FAN, False).save("out/grad3-2-1.png")
    test3xx(lena256bw, ERROR_FAN, False).save("out/lena3-2-1.png")
    Svg(ERROR_SHIAUFAN).save("out/fig3-2-1b.png", 40)
    test3xx(grad256bw, ERROR_SHIAUFAN, False).save("out/grad3-2-1b.png")
    test3xx(lena256bw, ERROR_SHIAUFAN, False).save("out/lena3-2-1b.png")
    Svg(ERROR_SHIAUFAN2).save("out/fig3-2-1c.png", 40)
    test3xx(grad256bw, ERROR_SHIAUFAN2, False).save("out/grad3-2-1c.png")
    test3xx(lena256bw, ERROR_SHIAUFAN2, False).save("out/lena3-2-1c.png")

    Svg(ERROR_STUCKI).save("out/fig3-2-3.png", 40)
    test3xx(grad256bw, ERROR_STUCKI, False).save("out/grad3-2-3.png")
    test3xx(lena256bw, ERROR_STUCKI, False).save("out/lena3-2-3.png")

    Svg(ERROR_BURKES).save("out/fig3-2-4.png", 40)
    test3xx(grad256bw, ERROR_BURKES, False).save("out/grad3-2-4.png")
    test3xx(lena256bw, ERROR_BURKES, False).save("out/lena3-2-4.png")

    Svg(ERROR_SIERRA).save("out/fig3-2-5.png", 40)
    test3xx(grad256bw, ERROR_SIERRA, False).save("out/grad3-2-5.png")
    test3xx(lena256bw, ERROR_SIERRA, False).save("out/lena3-2-5.png")

    Svg(ERROR_SIERRA2).save("out/fig3-2-6.png", 40)
    test3xx(grad256bw, ERROR_SIERRA2, False).save("out/grad3-2-6.png")
    test3xx(lena256bw, ERROR_SIERRA2, False).save("out/lena3-2-6.png")

    Svg(ERROR_FILTERLITE).save("out/fig3-2-7.png", 40)
    test3xx(grad256bw, ERROR_FILTERLITE, False).save("out/grad3-2-7.png")
    test3xx(lena256bw, ERROR_FILTERLITE, False).save("out/lena3-2-7.png")

    Svg(ERROR_ATKINSON).save("out/fig3-2-8.png", 40)
    test3xx(grad256bw, ERROR_ATKINSON, False).save("out/grad3-2-8.png")
    test3xx(lena256bw, ERROR_ATKINSON, False).save("out/lena3-2-8.png")

    #test3xx(grad256bw, ERROR_STAR, False).save("out/grad3-2-9.png")
    #test3xx(lena256bw, ERROR_STAR, False).save("out/lena3-2-9.png")

    #test3xx(grad256bw, ERROR_STAR, False).save("out/grad3-2-9.png")
    #test3xx(lena256bw, ERROR_STAR, False).save("out/lena3-2-9.png")

# Output 3.3.1: Floyd-Steinberg on grey 90%
# Output 3.3.2: serpentine Floyd-Steinberg on grey 90%
if chapter(3):
    tmp = Image((128, 128))
    for x, y in rangexy(128, 128):
        tmp.setGray(x, y, 0.90)
    test3xx(tmp, ERROR_FSTEIN, True).getZoom(2).save("out/lena3-3-2.png")
    test3xx(tmp, ERROR_FSTEIN, False).getZoom(2).save("out/lena3-3-1.png")

# Output 3.3.3: Riemersma dither on a Hilbert curve
# Output 3.3.4: Riemersma dither on a Hilbert 2 curve
def hilbert(x, y, n):
    d1 = [0, 4, 3, 2, 1]
    d2 = [0, 3, 4, 1, 2]
    m = n/2

    if x == n - 1 and y == 0: return 0
    if x == 0 and y == m - 1: return 1
    if x == m - 1 and y == m: return 4
    if x == n - 1 and y == m: return 2

    if y >= m: return hilbert(x % m, y - m, m)
    if x < m: return d1[hilbert(y, x, m)]
    else: return d2[hilbert(m - 1 - y, n - 1 - x, m)]

def hilbert2(x, y, n):
    d1 = [0, 2, 1, 3, 4]
    d2 = [0, 1, 2, 4, 3]
    d3 = [0, 2, 1, 4, 3]
    m = n/3

    if x == n - 1 and y == n - 1: return 0
    if x == m - 1 and y == m - 1: return 4
    if x == 2 * m - 1 and y == 0: return 4
    if x == n - 1 and y == m - 1: return 1
    if x == 0 and y == 2 * m - 1: return 1
    if x == m and y == m: return 3
    if x == m * 2 and y == m * 2 - 1: return 3
    if x == m - 1 and y == n - 1: return 4
    if x == 2 * m - 1 and y == 2 * m: return 4

    if (x < m or x >= m * 2) and (y < m or y >= m * 2):
        return hilbert2(x % m, y % m, m)
    if x >= m and x < m * 2 and y >= m and y < m * 2:
        return d3[hilbert2(2 * m - 1 - x, 2 * m - 1 - y, m)]
    if x >= m and x < m * 2:
        return d1[hilbert2(x - m, m - 1 - (y % m), m)]
    else: # if y >= m and y < m * 2
        return d2[hilbert2(m - 1 - (x % m), y - m, m)]

def riemersma(src, iterator, order):
    (w, h) = src.size()
    dest = Image((w, h))
    q = 16
    r = 16
    size = 1
    while size < w or size < h: size *= order
    coord = [(0, 0), (0, 1), (0, -1), (-1, 0), (1, 0)]
    err = [0] * q
    list = [math.exp(math.log(r) * i / (q - 1)) / r for i in range(q)]
    (x, y) = (0, 0)
    out = False
    while True:
        if not out:
            a = c = src.getGray(x, y)
            for i in range(q):
                a += err[i] * list[i]
            d = a > 0.5
            dest.setGray(x, y, d)
            for i in range(q - 1):
                err[i] = err[i + 1]
            err[q - 1] = c - d
        t = iterator(x, y, size)
        if t == 0:
            break
        dx, dy = coord[t]
        x += dx
        y += dy
        if out:
            if x < w and y < h:
                err = [0] * q
                out = False
            continue
        # Did we fall off the screen?
        out = (x > w + q or y > h + q)
    return dest

if chapter(3):
    riemersma(grad256bw, hilbert, 2).save("out/grad3-3-3.png")
    riemersma(lena256bw, hilbert, 2).save("out/lena3-3-3.png")
    riemersma(grad256bw, hilbert2, 3).save("out/grad3-3-4.png")
    riemersma(lena256bw, hilbert2, 3).save("out/lena3-3-4.png")

# Output 3.3.5: spatial Hilbert dither on a Hilbert curve
# Output 3.3.6: spatial Hilbert dither on a Hilbert 2 curve
def spatialhilbert(src, iterator, order):
    (w, h) = src.size()
    dest = Image((w, h))
    q = 16
    size = 1
    while size < w or size < h: size *= order
    coord = [(0, 0), (0, 1), (0, -1), (-1, 0), (1, 0)]
    err = [0] * q
    dx = [0] * q
    dy = [0] * q
    dist = [0] * q
    (x, y) = (0, 0)
    out = False
    while True:
        if not out:
            c = src.getGray(x, y) + err[0]
            d = c > 0.5
            dest.setGray(x, y, d)
        t = iterator(x, y, size)
        if t == 0:
            break
        dx[0], dy[0] = coord[t]
        if not out:
            error = c - d
            errdiv = 0.
            for i in range(q - 1):
                t = coord[iterator(x + dx[i], y + dy[i], size)]
                dx[i + 1] = dx[i] + t[0]
                dy[i + 1] = dy[i] + t[1]
            for i in range(q):
                dist[i] = dx[i] * dx[i] + dy[i] * dy[i]
                errdiv += 1. / dist[i]
            error /= errdiv
            for i in range(q - 1):
                err[i] = err[i + 1] + error / dist[i]
            err[q - 1] = error / dist[q - 1]
        x += dx[0]
        y += dy[0]
        if out:
            if x < w and y < h:
                err = [0] * q
                out = False
            continue
        # Did we fall off the screen?
        out = (x > w + q or y > h + q)
    return dest

if chapter(3):
    spatialhilbert(grad256bw, hilbert, 2).save("out/grad3-3-5.png")
    spatialhilbert(lena256bw, hilbert, 2).save("out/lena3-3-5.png")
    spatialhilbert(grad256bw, hilbert2, 3).save("out/grad3-3-6.png")
    spatialhilbert(lena256bw, hilbert2, 3).save("out/lena3-3-6.png")

# Output 3.3.7: Knuth's dot diffusion
# Output 3.3.8: Knuth's dot diffusion, sharpen = 0.9
# Output 3.3.9: Knuth's dot diffusion, sharpen = 0.9, zeta = 0.2
# Output 3.3.10: serpentine Floyd-Steinberg, sharpen = 0.9
def sharpen(src, sharpening):
    (w, h) = src.size()
    dest = Image((w, h))
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        total = 0.
        for j in range(-1, 2):
            for i in range(-1, 2):
                total += src.getGray(x + i, y + j)
        total /= 9.
        d = (c - sharpening * total) / (1.0 - sharpening)
        if d < 0.:
            d = 0.
        elif d > 1.:
            d = 1.
        dest.setGray(x, y, d)
    return dest

def test337(src, mat, propagate, zeta):
    (w, h) = src.size()
    dest = Image((w, h))
    dx = len(mat[0])
    dy = len(mat)
    # 0. analyse diffusion matrix to speed up things later
    diff = []
    cx, cy = -1, -1
    for x, y in rangexy(len(propagate[0]), len(propagate)):
        if propagate[y][x] == -1:
            cx, cy = x, y
    for x, y in rangexy(len(propagate[0]), len(propagate)):
        diff.append((x - cx, y - cy, propagate[y][x]))
    # 1. analyse order matrix to get equivalence classes
    nclasses = 0
    for l in mat:
        for v in l:
            if v + 1 > nclasses:
                nclasses = v + 1
    classes = [[] for n in range(nclasses)]
    for x, y in rangexy(w, h):
        classes[mat[y % dy][x % dx]].append((x, y))
    # 2. copy image
    img = [[src.getGray(x, y) for x in range(w)] for y in range(h)]
    aa = Matrix(w, h, 1.)
    # 3. parse all classes
    for l in classes:
        for x, y in l:
            c = img[y][x]
            if aa[y][x] == 1.:
                error = c + 4. * zeta
            else:
                error = c - zeta
                if x > 0 and aa[y][x-1] == 1.: error += zeta
                if y > 0 and aa[y-1][x] == 1.: error += zeta
                if x < w-1 and aa[y][x+1] == 1.: error += zeta
                if y < h-1 and aa[y+1][x] == 1.: error += zeta
            if c + error < 1.:
                aa[y][x] = 0.
                if x > 0 and aa[y][x-1] == 1.: aa[y][x-1] = .5
                if y > 0 and aa[y-1][x] == 1.: aa[y-1][x] = .5
                if x < w-1 and aa[y][x+1] == 1.: aa[y][x+1] = .5
                if y < h-1 and aa[y+1][x] == 1.: aa[y+1][x] = .5
            else:
                error = c - 1.
            # Propagate first line of error
            total = 0
            err = []
            for (i, j, e) in diff:
                if x + i < 0 or x + i >= w \
                    or y + j < 0 or y + j >= h:
                    continue
                n = mat[(y + j) % dy][(x + i) % dx] - mat[y % dy][x % dx]
                if n <= 0:
                    continue
                err.append((i, j, e))
                total += e
            for (i, j, e) in err:
                img[y + j][x + i] += error * e / total
    # 4. copy image, replacing grey with white
    for x, y in rangexy(w, h):
        dest.setGray(x, y, aa[y][x] > 0.)
    return dest

ERROR_DOT = \
    [[1./12, 1./6, 1./12],
     [ 1./6,   -1,  1./6],
     [1./12, 1./6, 1./12]]

if chapter(3):
    Svg(ERROR_DOT).save("out/fig3-3-7b.png", 40)
    test337(grad256bw, DITHER_CLUSTER88, ERROR_DOT, 0.).save("out/grad3-3-7.png")
    test337(lena256bw, DITHER_CLUSTER88, ERROR_DOT, 0.).save("out/lena3-3-7.png")
    tmp = sharpen(grad256bw, 0.9)
    test337(tmp, DITHER_CLUSTER88, ERROR_DOT, 0.).save("out/grad3-3-8.png")
    test337(tmp, DITHER_CLUSTER88, ERROR_DOT, 0.2).save("out/grad3-3-9.png")
    test3xx(tmp, ERROR_FSTEIN, True).save("out/grad3-3-10.png")
    tmp = sharpen(lena256bw, 0.9)
    test337(tmp, DITHER_CLUSTER88, ERROR_DOT, 0.).save("out/lena3-3-8.png")
    test337(tmp, DITHER_CLUSTER88, ERROR_DOT, 0.2).save("out/lena3-3-9.png")
    test3xx(tmp, ERROR_FSTEIN, True).save("out/lena3-3-10.png")

# Output 3.3.11: omni-directional error diffusion
def test3311(src):
    w, h = src.size()
    g = { 0: 0, 1: 1, 2: 1 }
    g[-1] = g[2] - g[1]
    g[-2] = g[1] - g[0]
    n = 2
    while g[n] < w and g[n] < h:
        n += 1
        g[n] = g[n - 1] + g[n - 3]
        g[-n] = g[-n + 3] - g[-n + 2]
    u = g[n - 2]
    v = g[n - 1]
    a = [[(i * u + j * v) % g[n] for i in range(g[n])] for j in range(g[n])]
    return a

ERROR_OMNI = \
    [[0.1, 0.2, 0.1],
     [0.1,  -1, 0.1],
     [0.1, 0.2, 0.1]]

if chapter(3):
    Svg(ERROR_OMNI).save("out/fig3-3-11.png", 40)

    mat = test3311(grad256bw)
    test337(grad256bw, mat, ERROR_OMNI, 0.).save("out/grad3-3-11.png")

    mat = test3311(lena256bw)
    tmp = [[str(mat[y][x]) for x in range(16)] for y in range(12)]
    for x in range(16): tmp[11][x] = "..."
    for y in range(12): tmp[y][15] = "..."
    Svg(tmp).save("out/fig3-3-11b.png", 20)
    test337(lena256bw, mat, ERROR_OMNI, 0.).save("out/lena3-3-11.png")

# Output 3.4.1: Ostromoukhov's variable error diffusion
def test341(src, serpentine):
    m = [[13, 0, 5], [13, 0, 5], [21, 0, 10], [7, 0, 4],
         [8, 0, 5], [47, 3, 28], [23, 3, 13], [15, 3, 8],
         [22, 6, 11], [43, 15, 20], [7, 3, 3], [501, 224, 211],
         [249, 116, 103], [165, 80, 67], [123, 62, 49], [489, 256, 191],
         [81, 44, 31], [483, 272, 181], [60, 35, 22], [53, 32, 19],
         [237, 148, 83], [471, 304, 161], [3, 2, 1], [481, 314, 185],
         [354, 226, 155], [1389, 866, 685], [227, 138, 125], [267, 158, 163],
         [327, 188, 220], [61, 34, 45], [627, 338, 505], [1227, 638, 1075],
         [20, 10, 19], [1937, 1000, 1767], [977, 520, 855], [657, 360, 551],
         [71, 40, 57], [2005, 1160, 1539], [337, 200, 247], [2039, 1240, 1425],
         [257, 160, 171], [691, 440, 437], [1045, 680, 627], [301, 200, 171],
         [177, 120, 95], [2141, 1480, 1083], [1079, 760, 513], [725, 520, 323],
         [137, 100, 57], [2209, 1640, 855], [53, 40, 19], [2243, 1720, 741],
         [565, 440, 171], [759, 600, 209], [1147, 920, 285], [2311, 1880, 513],
         [97, 80, 19], [335, 280, 57], [1181, 1000, 171], [793, 680, 95],
         [599, 520, 57], [2413, 2120, 171], [405, 360, 19], [2447, 2200, 57],
         [11, 10, 0], [158, 151, 3], [178, 179, 7], [1030, 1091, 63],
         [248, 277, 21], [318, 375, 35], [458, 571, 63], [878, 1159, 147],
         [5, 7, 1], [172, 181, 37], [97, 76, 22], [72, 41, 17],
         [119, 47, 29], [4, 1, 1], [4, 1, 1], [4, 1, 1],
         [4, 1, 1], [4, 1, 1], [4, 1, 1], [4, 1, 1],
         [4, 1, 1], [4, 1, 1], [65, 18, 17], [95, 29, 26],
         [185, 62, 53], [30, 11, 9], [35, 14, 11], [85, 37, 28],
         [55, 26, 19], [80, 41, 29], [155, 86, 59], [5, 3, 2],
         [5, 3, 2], [5, 3, 2], [5, 3, 2], [5, 3, 2],
         [5, 3, 2], [5, 3, 2], [5, 3, 2], [5, 3, 2],
         [5, 3, 2], [5, 3, 2], [5, 3, 2], [5, 3, 2],
         [305, 176, 119], [155, 86, 59], [105, 56, 39], [80, 41, 29],
         [65, 32, 23], [55, 26, 19], [335, 152, 113], [85, 37, 28],
         [115, 48, 37], [35, 14, 11], [355, 136, 109], [30, 11, 9],
         [365, 128, 107], [185, 62, 53], [25, 8, 7], [95, 29, 26],
         [385, 112, 103], [65, 18, 17], [395, 104, 101], [4, 1, 1]]
    (w, h) = src.size()
    dest = Image((w, h))
    ey = [0.] * (w + 2)
    for y in range(h):
        ex = 0
        newey = [0.] * (w + 2)
        if serpentine and y & 1:
            xrange = range(w - 1, -1, -1)
        else:
            xrange = range(w)
        for x in xrange:
            # Set pixel
            c = src.getGray(x, y) + ex + ey[x + 1]
            d = c > 0.5
            dest.setGray(x, y, d)
            error = c - d
            i = (int)(c * 255.9999)
            if i > 127:
                i = 255 - i
            (d1, d2, d3) = m[i]
            t = d1 + d2 + d3
            # Propagate error
            ex = error * d1 / t
            if serpentine and y & 1:
                newey[x + 2] += error * d3 / t
                newey[x + 1] += error * d2 / t
            else:
                newey[x] += error * d2 / t
                newey[x + 1] += error * d3 / t
        ey = newey
    return dest

if chapter(3):
    mat = [[0, -1,
            '<tspan style="font-style:italic">d1(i)</tspan>'],
           ['<tspan style="font-style:italic">d2(i)</tspan>',
            '<tspan style="font-style:italic">d3(i)</tspan>', 0]]
    Svg(mat).save("out/fig3-4-1.png", 40)
    test341(grad256bw, True).save("out/grad3-4-1.png")
    test341(lena256bw, True).save("out/lena3-4-1.png")

def test351(src, mat, tiles, diff, serpentine, glob):
    random.seed(0)
    (w, h) = src.size()
    dest = Image((w, h))
    ntiles = len(tiles)
    ty = len(tiles[0])
    tx = len(tiles[0][0])
    cur = Matrix(tx, ty, 0.)
    w, h = w / tx, h / ty
    lines = len(mat)
    rows = len(mat[0])
    offset = mat[0].index(-1)
    ey = Matrix(w + rows - 1, lines, 0.)
    for y in range(h):
        ex = [0.] * (rows - offset)
        if serpentine and y & 1:
            xrange = range(w - 1, -1, -1)
        else:
            xrange = range(w)
        for x in xrange:
            # Get block value
            for i, j in rangexy(tx, ty):
                e = ex[0] + ey[0][x + offset]
                cur[j][i] = src.getGray(x * tx + i, y * ty + j) + diff[j][i] * e
            # Compute closest block and associated error
            dist = tx * ty
            for n in range(ntiles):
                e = 0.
                d1 = 0.
                d2 = random.random() / 1000.
                for i, j in rangexy(tx, ty):
                    e += cur[j][i] - tiles[n][j][i]
                    d1 += diff[j][i] * (cur[j][i] - tiles[n][j][i])
                    d2 += diff[j][i] * abs(cur[j][i] - tiles[n][j][i])
                if glob:
                    d = abs(d1) + d2 / 1000.
                else:
                    d = abs(d1) / (tx * ty) + d2
                if d < dist:
                    dist = d
                    error = e
                    best = n
            # Set pixel
            for i, j in rangexy(tx, ty):
                dest.setGray(x * tx + i, y * ty + j, tiles[best][j][i])
            # Propagate first line of error
            for dx in range(rows - offset - 2):
                ex[dx] = ex[dx + 1] + error * mat[0][offset + 1 + dx]
            ex[rows - offset - 2] = error * mat[0][rows - 1]
            # Propagate next lines
            if serpentine and y & 1:
                for dy in range(1, lines):
                    for dx in range(rows):
                        ey[dy][x + dx] += error * mat[dy][rows - 1 - dx]
            else:
                for dy in range(1, lines):
                    for dx in range(rows):
                        ey[dy][x + dx] += error * mat[dy][dx]
        for dy in range(lines - 1):
            ey[dy] = ey[dy + 1]
        ey[lines - 1] = [0.] * (w + rows - 1)
    return dest

LINES22 = \
    [[[0., 0.], [0., 0.]],
     [[0., 1.], [0., 1.]],
     [[1., 0.], [1., 0.]],
     [[1., 1.], [0., 0.]],
     [[0., 0.], [1., 1.]],
     [[1., 1.], [1., 1.]]]

SQUARES33 = \
    [[[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]],
     [[1., 0., 0.], [0., 0., 0.], [0., 0., 0.]],
     [[0., 1., 0.], [0., 0., 0.], [0., 0., 0.]],
     [[0., 0., 1.], [0., 0., 0.], [0., 0., 0.]],
     [[0., 0., 0.], [1., 0., 0.], [0., 0., 0.]],
     [[0., 0., 0.], [0., 1., 0.], [0., 0., 0.]],
     [[0., 0., 0.], [0., 0., 1.], [0., 0., 0.]],
     [[0., 0., 0.], [0., 0., 0.], [1., 0., 0.]],
     [[0., 0., 0.], [0., 0., 0.], [0., 1., 0.]],
     [[0., 0., 0.], [0., 0., 0.], [0., 0., 1.]],
     [[1., 1., 0.], [1., 1., 0.], [0., 0., 0.]],
     [[0., 1., 1.], [0., 1., 1.], [0., 0., 0.]],
     [[0., 0., 0.], [1., 1., 0.], [1., 1., 0.]],
     [[0., 0., 0.], [0., 1., 1.], [0., 1., 1.]],
     [[1., 1., 1.], [1., 1., 1.], [1., 1., 1.]]]

TILES22 = []
for n in range(1 << 4):
    mat = Matrix(2, 2)
    for x, y in rangexy(2, 2):
        if n & (1 << (y * 2 + x)):
            mat[y][x] = 1.
    TILES22.append(mat)

DIFF_EVEN22 = \
    [[1./4, 1./4],
     [1./4, 1./4]]

DIFF_EVEN33 = \
    [[1./9, 1./9, 1./9],
     [1./9, 1./9, 1./9],
     [1./9, 1./9, 1./9]]

if chapter(3):
    def colorise(val):
        if type(val) == str:
            return '#fbb'
        return '#fff'

    mat = [[.25, .25],
           [.25, .25]]
    Svg(mat).save("out/fig3-5-2b.png", 40)
    mat = [[0, 0, 'a', 'b', 7./64, 7./64],
           [0, 0, 'c', 'd', 7./64, 7./64],
           [3./64, 3./64, 5./64, 5./64, 1./64, 1./64],
           [3./64, 3./64, 5./64, 5./64, 1./64, 1./64]]
    Svg(mat, colorise).save("out/fig3-5-2.png", 40)

    test351(grad256bw, [[-1, 0]],
            LINES22, DIFF_EVEN22, True, True).save("out/grad3-5-1.png")
    test351(lena256bw, [[-1, 0]],
            LINES22, DIFF_EVEN22, True, True).save("out/lena3-5-1.png")
    test351(grad256bw, ERROR_FSTEIN,
            LINES22, DIFF_EVEN22, True, True).save("out/grad3-5-2.png")
    test351(lena256bw, ERROR_FSTEIN,
            LINES22, DIFF_EVEN22, True, True).save("out/lena3-5-2.png")

    test351(grad256bw, ERROR_FSTEIN,
            SQUARES33, DIFF_EVEN33, True, True).save("out/grad3-5-3.png")
    test351(lena256bw, ERROR_FSTEIN,
            SQUARES33, DIFF_EVEN33, True, True).save("out/lena3-5-3.png")

    test351(grad256bw, ERROR_FSTEIN,
            TILES22, DIFF_EVEN22, True, True).save("out/grad3-5-4.png")
    test351(lena256bw, ERROR_FSTEIN,
            TILES22, DIFF_EVEN22, True, True).save("out/lena3-5-4.png")

DIFF_WEIGHTED22 = \
    [[51./128, 33./128],
     [25./128, 19./128]]

if chapter(3):
    mat = [[6./16, 5./16],
           [3./16, 2./16]]
    Svg(mat).save("out/fig3-5-5b.png", 40)
    mat = [[0, 0, 'a', 'b', 6*7./256, 5*7./256],
           [0, 0, 'c', 'd', 3*7./256, 2*7./256],
           [6*3./256, 5*3./256, 6*5./256, 5*5./256, 6*1./256, 5*1./256],
           [3*3./256, 2*3./256, 3*5./256, 2*5./256, 3*1./256, 2*1./256]]
    Svg(mat, colorise).save("out/fig3-5-5.png", 40)

    test351(grad256bw, ERROR_FSTEIN,
            TILES22, DIFF_WEIGHTED22, True, False).save("out/grad3-5-5.png")
    test351(lena256bw, ERROR_FSTEIN,
            TILES22, DIFF_WEIGHTED22, True, False).save("out/lena3-5-5.png")

# Output 3.6.1: sub-block error diffusion
def subblock(src, tiles, propagate, diff, gamma):
    (w, h) = src.size()
    # Gamma correction
    if gamma:
        ctoi = Gamma.CtoI
        itoc = Gamma.ItoC
    else:
        ctoi = itoc = lambda x : x
    # Propagating the error to a temporary buffer is becoming more and
    # more complicated. We decide to use an intermediate matrix instead.
    tmp = Matrix(w, h, 0.)
    for x, y in rangexy(w, h):
        tmp[y][x] = ctoi(src.getGray(x, y))
    dest = Image((w, h))
    # Analyse tile list
    ntiles = len(tiles)
    ty = len(tiles[0])
    tx = len(tiles[0][0])
    cur = Matrix(tx, ty, 0.)
    w, h = w / tx, h / ty
    # Analyse error propagate list
    for x, y in rangexy(w, h):
        # Get block value
        for i, j in rangexy(tx, ty):
            cur[j][i] = itoc(tmp[y * ty + j][x * tx + i])
        # Select closest block
        dist = tx * ty
        for n in range(ntiles):
            d = 0.
            e = 0.
            for i, j in rangexy(tx, ty):
                d += cur[j][i] - tiles[n][j][i]
                e += diff[j][i] * abs(cur[j][i] - tiles[n][j][i])
            if abs(d) / (tx * ty) + e < dist:
                dist = abs(d) / (tx * ty) + e
                best = n
        # Set pixel
        for i, j in rangexy(tx, ty):
            dest.setGray(x * tx + i, y * ty + j, tiles[best][j][i])
        # Propagate error
        for i, j in rangexy(tx, ty):
            e = ctoi(cur[j][i]) - ctoi(tiles[best][j][i])
            m = propagate[j][i]
            for px, py in rangexy(len(m[0]), len(m)):
                if m[py][px] == 0:
                    continue
                if m[py][px] == -1:
                    cx, cy = px, py
                    continue
                tmpx = x * tx + i + px - cx
                tmpy = y * ty + j + py - cy
                if tmpx > w * tx - 1 or tmpy > h * ty - 1:
                    continue
                tmp[tmpy][tmpx] += m[py][px] * e
    return dest

ERROR_SUBFS22 = \
    [[[[0, -1, 0, 8./64],
       [0, 0, 0, 10./64],
       [7./64, 22./64, 15./64, 2./64]],
      [[0, 0, -1, 20./64],
       [0, 0, 0, 14./64],
       [2./64, 11./64, 15./64, 2./64]]],
     [[[0, 0, 0, 0./64],
       [0, -1, 0, 6./64],
       [12./64, 32./64, 13./64, 1./64]],
      [[0, 0, 0, 0./64],
       [0, 0, -1, 20./64],
       [0./64, 12./64, 28./64, 4./64]]]]

if chapter(3):
    subblock(grad256bw, TILES22,
             ERROR_SUBFS22, DIFF_WEIGHTED22, False).save("out/grad3-6-1.png")
    subblock(lena256bw, TILES22,
             ERROR_SUBFS22, DIFF_WEIGHTED22, False).save("out/lena3-6-1.png")

    subblock(grad256bw, LINES22,
             ERROR_SUBFS22, DIFF_WEIGHTED22, False).save("out/grad3-6-2.png")
    subblock(lena256bw, LINES22,
             ERROR_SUBFS22, DIFF_WEIGHTED22, False).save("out/lena3-6-2.png")

    def colorise(val):
        if val == '':
            return '#ccc'
        if type(val) == str:
            return '#fbb'
        return '#fff'

    # XXX: hack, we modify ERROR_SUBFS22 because it's so much more convenient
    for x, y in rangexy(2, 2):
        for j in range(2):
            for i in range(1, 3):
                ERROR_SUBFS22[y][x][j][i] = ''
        ERROR_SUBFS22[y][x][y][1 + x] = chr(ord('a') + 2 * y + x)
    Svg(ERROR_SUBFS22[0][0], colorise).save("out/fig3-6-1a.png", 40)
    Svg(ERROR_SUBFS22[0][1], colorise).save("out/fig3-6-1b.png", 40)
    Svg(ERROR_SUBFS22[1][0], colorise).save("out/fig3-6-1c.png", 40)
    Svg(ERROR_SUBFS22[1][1], colorise).save("out/fig3-6-1d.png", 40)
    for x, y in rangexy(2, 2):
        for j in range(2):
            for i in range(1, 3):
                ERROR_SUBFS22[y][x][j][i] = 0
        ERROR_SUBFS22[y][x][y][1 + x] = -1

ERROR_SUBFS33 = \
    [[[[     0,     -1,      0,      0,  2./64],
       [     0,      0,      0,      0,  5./64],
       [     0,      0,      0,      0,  6./64],
       [ 5./64, 17./64, 17./64,  9./64,  1./64]],
      [[     0,      0,     -1,      0,  6./64],
       [     0,      0,      0,      0,  9./64],
       [     0,      0,      0,      0,  8./64],
       [ 2./64, 11./64, 16./64, 11./64,  1./64]],
      [[     0,      0,      0,     -1, 20./64],
       [     0,      0,      0,      0, 14./64],
       [     0,      0,      0,      0,  8./64],
       [     0,  3./64,  9./64,  9./64,  1./64]]],
     [[[     0,      0,      0,      0,      0],
       [     0,     -1,      0,      0,  2./64],
       [     0,      0,      0,      0,  5./64],
       [ 7./64, 23./64, 18./64,  8./64,  1./64]],
      [[     0,      0,      0,      0,      0],
       [     0,      0,     -1,      0,  6./64],
       [     0,      0,      0,      0,  9./64],
       [ 2./64, 12./64, 21./64, 13./64,  1./64]],
      [[     0,      0,      0,      0,      0],
       [     0,      0,      0,     -1, 20./64],
       [     0,      0,      0,      0, 14./64],
       [     0,  2./64, 11./64, 15./64,  2./64]]],
     [[[     0,      0,      0,      0,      0],
       [     0,      0,      0,      0,      0],
       [     0,     -1,      0,      0,  2./64],
       [12./64, 32./64, 14./64,  4./64,      0]],
      [[     0,      0,      0,      0,      0],
       [     0,      0,      0,      0,      0],
       [     0,      0,     -1,      0,  6./64],
       [     0, 12./64, 32./64, 13./64,  1./64]],
      [[     0,      0,      0,      0,      0],
       [     0,      0,      0,      0,      0],
       [     0,      0,      0,     -1, 20./64],
       [     0,      0, 12./64, 28./64,  4./64]]]]

TILES33 = []
for n in range(1 << 9):
    mat = Matrix(3, 3)
    for x, y in rangexy(3, 3):
        if n & (1 << (y * 3 + x)):
            mat[y][x] = 1.
    TILES33.append(mat)

DIFF_WEIGHTED33 = \
    [[15./64, 10./64,  6./64],
     [10./64,  6./64,  4./64],
     [ 6./64,  4./64,  3./64]]

if chapter(3):
    subblock(grad256bw, TILES33,
             ERROR_SUBFS33, DIFF_WEIGHTED33, False).save("out/grad3-6-3.png")
    subblock(lena256bw, TILES33,
             ERROR_SUBFS33, DIFF_WEIGHTED33, False).save("out/lena3-6-3.png")
    subblock(grad256bw, SQUARES33,
             ERROR_SUBFS33, DIFF_WEIGHTED33, False).save("out/grad3-6-4.png")
    subblock(lena256bw, SQUARES33,
             ERROR_SUBFS33, DIFF_WEIGHTED33, False).save("out/lena3-6-4.png")

    # XXX: hack, we modify ERROR_SUBFS33 because it's so much more convenient
    for x, y in rangexy(3, 3):
        for j in range(3):
            for i in range(1, 4):
                ERROR_SUBFS33[y][x][j][i] = ''
        ERROR_SUBFS33[y][x][y][1 + x] = chr(ord('a') + 3 * y + x)
    Svg(ERROR_SUBFS33[0][0], colorise).save("out/fig3-6-3a.png", 30)
    Svg(ERROR_SUBFS33[0][1], colorise).save("out/fig3-6-3b.png", 30)
    Svg(ERROR_SUBFS33[0][2], colorise).save("out/fig3-6-3c.png", 30)
    Svg(ERROR_SUBFS33[1][0], colorise).save("out/fig3-6-3d.png", 30)
    Svg(ERROR_SUBFS33[1][1], colorise).save("out/fig3-6-3e.png", 30)
    Svg(ERROR_SUBFS33[1][2], colorise).save("out/fig3-6-3f.png", 30)
    Svg(ERROR_SUBFS33[2][0], colorise).save("out/fig3-6-3g.png", 30)
    Svg(ERROR_SUBFS33[2][1], colorise).save("out/fig3-6-3h.png", 30)
    Svg(ERROR_SUBFS33[2][2], colorise).save("out/fig3-6-3i.png", 30)
    for x, y in rangexy(3, 3):
        for j in range(3):
            for i in range(1, 4):
                ERROR_SUBFS33[y][x][j][i] = 0
        ERROR_SUBFS33[y][x][y][1 + x] = -1

##############################################################################
if chapter(4):
    print "Chapter 4. Model-based dithering"

# Output 4.1.1: gaussian HVS applied to 8x8 Bayer dither, sigma = 1
# Output 4.1.2: gaussian HVS applied to 8x8 Bayer dither, sigma = 2
def gaussian(n, sigma):
    m = Matrix(n, n, 0.)
    t = 0.
    for x, y in rangexy(n, n):
        i = x - (float)(n - 1.) / 2.
        j = y - (float)(n - 1.) / 2.
        v = math.pow(math.e, - (i * i + j * j) / (2. * sigma * sigma))
        m[y][x] = v
        t += v
    for x, y in rangexy(n, n):
        m[y][x] /= t
    return m

def convolution(src, m):
    (w, h) = src.size()
    dest = Image((w, h))
    dy = len(m)
    dx = len(m[0])
    srcmat = [[src.getGray(x, y) for x in range(w)] for y in range(h)]
    for x, y in rangexy(w, h):
        c = t = 0.
        for i, j in rangexy(dx, dy):
            u = i - (dx - 1) / 2
            v = j - (dy - 1) / 2
            if x + u >= w or y + v >= h or x + u < 0 or y + v < 0:
                continue
            c += srcmat[y + v][x + u] * m[j][i]
            t += m[j][i]
        dest.setGray(x, y, c / t)
    return dest

if chapter(4):
    tmp = Image("out/grad2-3-1b.png")
    convolution(tmp, gaussian(11, 1.)).save("out/grad4-1-1.png")
    convolution(tmp, gaussian(11, 2.)).save("out/grad4-1-2.png")
    tmp = Image("out/lena2-3-1b.png")
    convolution(tmp, gaussian(11, 1.)).save("out/lena4-1-1.png")
    convolution(tmp, gaussian(11, 2.)).save("out/lena4-1-2.png")

# Output 4.2.1: direct binary search, iteration 0
# Output 4.2.2: direct binary search, iteration 1
# Output 4.2.3: direct binary search, iteration 2
# Output 4.2.4: direct binary search, iteration 5
def test42x(src):
    random.seed(0)
    (w, h) = src.size()
    dest = Image((w, h))
    for x, y in rangexy(w, h):
        c = src.getGray(x, y) + random.random() - 0.5
        d = c > 0.5
        dest.setGray(x, y, d)
    return dest

def test42y(src, dest, hvs):
    threshold = 0.4
    kernel = Matrix(8, 8, 0.) # have two borders of zeroes
    for i, j in rangexy(6, 6):
         kernel[j][i] = hvs(i * i + j * j)
    (w, h) = src.size()
    # Build fast pixel lookup tables
    srcmat = Matrix(w, h, 0.)
    destmat = Matrix(w, h, 0.)
    for x, y in rangexy(w, h):
        srcmat[y][x] = src.getGray(x, y)
        destmat[y][x] = dest.getGray(x, y)
    # Build human perception model for both source and destination
    srchvs = Matrix(w, h, 0.)
    desthvs = Matrix(w, h, 0.)
    for x, y in rangexy(w, h):
        srcp = destp = 0.
        for j in range(-5, 6):
            if y + j < 0 or y + j >= h:
                continue
            for i in range(-3, 4):
                if x + i < 0 or x + i >= w:
                    continue
                m = kernel[abs(j)][abs(i)]
                srcp += m * srcmat[y + j][x + i]
                destp += m * destmat[y + j][x + i]
        srchvs[y][x] = srcp
        desthvs[y][x] = destp
    swaps = toggles = 0
    for x, y in rangexy(w, h):
        d = destmat[y][x]
        best = 0.
        # Compute the effect of a toggle
        e = 0.
        for j in range(-5, 6):
            if y + j < 0 or y + j >= h:
                continue
            for i in range(-5, 6):
                if x + i < 0 or x + i >= w:
                    continue
                m = kernel[abs(j)][abs(i)]
                p = srchvs[y + j][x + i]
                q1 = desthvs[y + j][x + i]
                q2 = q1 - m * d + m * (1. - d)
                e += abs(q1 - p) - abs(q2 - p)
        if e > best:
            best = e
            op = False
        # Compute the effect of swaps
        for dx, dy in [(0, 1), (0, -1), (-1, 0), (1, 0)]:
            if y + dy < 0 or y + dy >= h or x + dx < 0 or x + dx >= w:
                continue
            d2 = destmat[y + dy][x + dx]
            if d2 == d:
                continue
            e = 0.
            for j in range(-6, 7):
                for i in range(-6, 7):
                    if y + j < 0 or y + j >= h or x + i < 0 or x + i >= w:
                        continue
                    ma = kernel[abs(j)][abs(i)]
                    mb = kernel[abs(j - dy)][abs(i - dx)]
                    p = srchvs[y + j][x + i]
                    q1 = desthvs[y + j][x + i]
                    q2 = q1 - ma * d + ma * d2 - mb * d2 + mb * d
                    e += abs(q1 - p) - abs(q2 - p)
            if e > best:
                best = e
                op = (dx, dy)
        # Apply the change if interesting
        if best <= 0.:
            continue
        if op:
            dx, dy = op
            d2 = destmat[y + dy][x + dx]
            destmat[y + dy][x + dx] = d
        else:
            d2 = 1. - d
        destmat[y][x] = d2
        for j in range(-5, 6):
            for i in range(-5, 6):
                m = kernel[abs(j)][abs(i)]
                if y + j >= 0 and y + j < h and x + i >= 0 and x + i < w:
                    desthvs[y + j][x + i] -= m * d
                    desthvs[y + j][x + i] += m * d2
                if op and y + dy + j >= 0 and y + dy + j < h \
                   and x + dx + i >= 0 and x + dx + i < w:
                    desthvs[y + dy + j][x + dx + i] -= m * d2
                    desthvs[y + dy + j][x + dx + i] += m * d
    for x, y in rangexy(w, h):
        dest.setGray(x, y, destmat[y][x])
    return dest

if chapter(4):
    hvs = lambda x: math.pow(math.e, - math.sqrt(x))

    tmp = test42x(grad256bw)
    tmp.save("out/grad4-2-1.png")
    tmp = test42y(grad256bw, tmp, hvs)
    tmp.save("out/grad4-2-2.png")
    tmp = test42y(grad256bw, tmp, hvs)
    tmp.save("out/grad4-2-3.png")
    for n in range(3): tmp = test42y(grad256bw, tmp, hvs)
    tmp.save("out/grad4-2-4.png")

    tmp = test42x(lena256bw)
    tmp.save("out/lena4-2-1.png")
    tmp = test42y(lena256bw, tmp, hvs)
    tmp.save("out/lena4-2-2.png")
    tmp = test42y(lena256bw, tmp, hvs)
    tmp.save("out/lena4-2-3.png")
    for n in range(3): tmp = test42y(lena256bw, tmp, hvs)
    tmp.save("out/lena4-2-4.png")

if chapter(4):
    hvs = lambda x: math.pow(math.e, -x / 2.)
    tmp = test42x(grad256bw)
    for n in range(5): tmp = test42y(grad256bw, tmp, hvs)
    tmp.save("out/grad4-2-5.png")
    tmp = test42x(lena256bw)
    for n in range(5): tmp = test42y(lena256bw, tmp, hvs)
    tmp.save("out/lena4-2-5.png")

if chapter(4):
    hvs = lambda x: math.pow(math.e, -x / 4.5)
    tmp = test42x(grad256bw)
    for n in range(5): tmp = test42y(grad256bw, tmp, hvs)
    tmp.save("out/grad4-2-6.png")
    tmp = test42x(lena256bw)
    for n in range(5): tmp = test42y(lena256bw, tmp, hvs)
    tmp.save("out/lena4-2-6.png")

if chapter(4):
    hvs = lambda x: math.pow(math.e, -x / 8.)
    tmp = test42x(grad256bw)
    for n in range(5): tmp = test42y(grad256bw, tmp, hvs)
    tmp.save("out/grad4-2-7.png")
    tmp = test42x(lena256bw)
    for n in range(5): tmp = test42y(lena256bw, tmp, hvs)
    tmp.save("out/lena4-2-7.png")

if chapter(4):
    hvs = lambda x: math.pow(math.e, -x / 8.) + 2. * math.pow(math.e, -x / 1.5)
    tmp = test42x(grad256bw)
    for n in range(5): tmp = test42y(grad256bw, tmp, hvs)
    tmp.save("out/grad4-2-8.png")
    tmp = test42x(lena256bw)
    for n in range(5): tmp = test42y(lena256bw, tmp, hvs)
    tmp.save("out/lena4-2-8.png")

# Chapter 4.3: RMSE computations
def rmse(gray, dither):
    (w, h) = gray.size()
    error = 0.
    for y in range(5, h - 5):
        for x in range(5, w - 5):
            c = gray.getGray(x, y)
            d = dither.getGray(x, y)
            error += (c - d) * (c - d)
    return error / ((w - 10) * (h - 10))

if chapter(4):
    sigmas = [1., 1.5, 2.]
    tmp = [convolution(lena256bw, gaussian(11, s)) for s in sigmas]
    for name in ["1-1-1", "1-1-2", "1-1-3", "1-3-1", "1-4-1", "1-4-2",
                 "2-1-1", "2-3-0", "2-3-1", "2-3-1b", "2-3-2", "2-3-2b",
                 "2-3-3", "2-4-1", "2-4-2", "2-5-1", "2-5-2", "2-5-3",
                 "2-5-4", "2-6-1", "2-6-2", "2-6-3", "2-6-4", "2-6-5",
                 "2-6-6", "2-7-1", "2-7-2",
                 "3-0-1", "3-1-1", "3-1-2", "3-1-3", "3-2-1", "3-2-1b",
                 "3-2-1c", "3-2-3", "3-2-4", "3-2-5", "3-2-6", "3-2-7",
                 "3-2-8", "3-3-3", "3-3-4", "3-3-5", "3-3-6", "3-3-7",
                 "3-3-8", "3-3-9", "3-3-10", "3-3-11", "3-4-1", "3-5-1",
                 "3-5-2", "3-5-3", "3-5-4", "3-5-5", "3-6-1", "3-6-2",
                 "3-6-3", "3-6-4",
                 "4-2-2", "4-2-3", "4-2-4", "4-2-5", "4-2-6", "4-2-7",
                 "4-2-8"]:
        img = Image("out/lena" + name + ".png")
        for i in range(3):
            tmp2 = convolution(img, gaussian(11, sigmas[i]))
            err = 100. * rmse(tmp[i], tmp2)
            print "[ERR] %.5f" % (err,)

##############################################################################
if chapter(5):
    print "Chapter 5. Greyscale dithering"

# Output 5.0.1: 4x4 Bayer dithering, 3 colours
def test501(src, mat):
    (w, h) = src.size()
    dest = Image((w, h))
    dx = len(mat[0])
    dy = len(mat)
    for x, y in rangexy(w, h):
        c = src.getGray(x, y)
        threshold = (1. + mat[y % dy][x % dx]) / (dx * dy + 1)
        if c < 0.5:
            c = 0.5 * (c > threshold / 2)
        else:
            c = 0.5 + 0.5 * (c > 0.5 + threshold / 2)
        dest.setGray(x, y, c)
    return dest

if chapter(5):
    test501(grad256bw, DITHER_BAYER88).save("out/grad5-0-1.png")
    test501(lena256bw, DITHER_BAYER88).save("out/lena5-0-1.png")

# Output 5.0.2: standard Floyd-Steinberg, 3 colours
def test502(src, mat, serpentine):
    (w, h) = src.size()
    dest = Image((w, h))
    lines = len(mat)
    rows = len(mat[0])
    offset = mat[0].index(-1)
    ey = Matrix(w + rows - 1, lines, 0.)
    for y in range(h):
        ex = [0.] * (rows - offset)
        if serpentine and y & 1:
            xrange = range(w - 1, -1, -1)
        else:
            xrange = range(w)
        for x in xrange:
            # Set pixel
            c = src.getGray(x, y) + ex[0] + ey[0][x + offset]
            d = 0.5 * (c > 0.25) + 0.5 * (c > 0.75)
            dest.setGray(x, y, d)
            error = c - d
            # Propagate first line of error
            for dx in range(rows - offset - 2):
                ex[dx] = ex[dx + 1] + error * mat[0][offset + 1 + dx]
            ex[rows - offset - 2] = error * mat[0][rows - 1]
            # Propagate next lines
            if serpentine and y & 1:
                for dy in range(1, lines):
                    for dx in range(rows):
                        ey[dy][x + dx] += error * mat[dy][rows - 1 - dx]
            else:
                for dy in range(1, lines):
                    for dx in range(rows):
                        ey[dy][x + dx] += error * mat[dy][dx]
        for dy in range(lines - 1):
            ey[dy] = ey[dy + 1]
        ey[lines - 1] = [0.] * (w + rows - 1)
    return dest

if chapter(5):
    test502(grad256bw, ERROR_FSTEIN, True).save("out/grad5-0-2.png")
    test502(lena256bw, ERROR_FSTEIN, True).save("out/lena5-0-2.png")

# Pattern 5.1.1: gamma-corrected 50% grey, black-white halftone, 50% grey
if chapter(5):
    dest = Image((240, 80))
    for y in range(80):
        for x in range(80):
            dest.setGray(x, y, Gamma.ItoC(0.5))
        for x in range(80, 160):
            c = (x + y) & 1
            dest.setGray(x, y, c)
        for x in range(160, 240):
            dest.setGray(x, y, 0.5)
    dest.save("out/pat5-1-1.png")

# Output 5.2.1: gamma-corrected 2-colour Floyd-Steinberg
# Output 5.2.2: gamma-corrected 3-colour Floyd-Steinberg
# Output 5.2.3: gamma-corrected 4-colour Floyd-Steinberg
def test52x(src, mat, serpentine, threshold):
    (w, h) = src.size()
    dest = Image((w, h))
    lines = len(mat)
    rows = len(mat[0])
    offset = mat[0].index(-1)
    ey = Matrix(w + rows - 1, lines, 0.)
    for y in range(h):
        ex = [0.] * (rows - offset)
        if serpentine and y & 1:
            xrange = range(w - 1, -1, -1)
        else:
            xrange = range(w)
        for x in xrange:
            # Set pixel
            c = Gamma.CtoI(src.getGray(x, y)) + ex[0] + ey[0][x + offset]
            d = threshold(c)
            dest.setGray(x, y, Gamma.ItoC(d))
            error = c - d
            # Propagate first line of error
            for dx in range(rows - offset - 2):
                ex[dx] = ex[dx + 1] + error * mat[0][offset + 1 + dx]
            ex[rows - offset - 2] = error * mat[0][rows - 1]
            # Propagate next lines
            if serpentine and y & 1:
                for dy in range(1, lines):
                    for dx in range(rows):
                        ey[dy][x + dx] += error * mat[dy][rows - 1 - dx]
            else:
                for dy in range(1, lines):
                    for dx in range(rows):
                        ey[dy][x + dx] += error * mat[dy][dx]
        for dy in range(lines - 1):
            ey[dy] = ey[dy + 1]
        ey[lines - 1] = [0.] * (w + rows - 1)
    return dest

if chapter(5):
    test52x(grad256bw, ERROR_FSTEIN, True, Gamma.Cto2).save("out/grad5-2-1.png")
    test52x(lena256bw, ERROR_FSTEIN, True, Gamma.Cto2).save("out/lena5-2-1.png")
    test52x(grad256bw, ERROR_FSTEIN, True, Gamma.Cto3).save("out/grad5-2-2.png")
    test52x(lena256bw, ERROR_FSTEIN, True, Gamma.Cto3).save("out/lena5-2-2.png")
    test52x(grad256bw, ERROR_FSTEIN, True, Gamma.Cto4).save("out/grad5-2-3.png")
    test52x(lena256bw, ERROR_FSTEIN, True, Gamma.Cto4).save("out/lena5-2-3.png")

# Output 5.3.1: full 4-colour block error diffusion
GREY22 = []
for n in range(4*4*4*4):
    vals = [0., 0.333, 0.666, 1.]
    a, b, c, d = n & 3, (n >> 2) & 3, (n >> 4) & 3, (n >> 6) & 3
    GREY22.append([[vals[a], vals[b]], [vals[c], vals[d]]])

if chapter(5):
    subblock(grad256bw, GREY22,
             ERROR_SUBFS22, DIFF_WEIGHTED22, True).save("out/grad5-3-1.png")
    subblock(lena256bw, GREY22,
             ERROR_SUBFS22, DIFF_WEIGHTED22, True).save("out/lena5-3-1.png")

# Output 5.3.2: 4-colour block error diffusion with only line tiles
GREYLINES22 = []
for n in range(4*4*4*4):
    vals = [0., 0.333, 0.666, 1.]
    a, b, c, d = n & 3, (n >> 2) & 3, (n >> 4) & 3, (n >> 6) & 3
    if (a != b or c != d) and (a != c or b != d):
        continue
    GREYLINES22.append([[vals[a], vals[b]], [vals[c], vals[d]]])

if chapter(5):
    subblock(grad256bw, GREYLINES22,
             ERROR_SUBFS22, DIFF_WEIGHTED22, True).save("out/grad5-3-2.png")
    subblock(lena256bw, GREYLINES22,
             ERROR_SUBFS22, DIFF_WEIGHTED22, True).save("out/lena5-3-2.png")

##############################################################################
if chapter(6):
    print "Chapter 6. Colour dithering"

# Pattern 6.1.1: 8-colour palette
if chapter(6):
    dest = Image((512, 64))
    for x in range(512):
        d = x / 64
        r = (d & 2) >> 1
        g = (d & 4) >> 2
        b = d & 1
        for y in range(64):
            dest.setRgb(x, y, r, g, b)
    dest.save("out/pat6-1-1.png")

# Figure 6.1.1: 128x128 Lena
# Figure 6.1.2a: red channel
# Figure 6.1.2b: green channel
# Figure 6.1.2c: blue channel
# Figure 6.1.3a: red channel promoted to greyscale
# Figure 6.1.3b: green channel promoted to greyscale
# Figure 6.1.3c: blue channel promoted to greyscale
# Figure 6.1.4a: dithered red channel
# Figure 6.1.4b: dithered green channel
# Figure 6.1.4c: dithered blue channel
# Figure 6.1.5: combined dithered channels
if chapter(6):
    tmp = gammascale(lena256, 2)
    tmp.save("out/fig6-1-1.png")
    (w, h) = tmp.size()
    dst = [Image((w, h), True) for i in range(3)]
    for x, y in rangexy(w, h):
        rgb = tmp.getRgb(x, y)
        dst[0].setRgb(x, y, rgb[0], 0, 0)
        dst[1].setRgb(x, y, 0, rgb[1], 0)
        dst[2].setRgb(x, y, 0, 0, rgb[2])
    dst[0].save("out/fig6-1-2a.png")
    dst[1].save("out/fig6-1-2b.png")
    dst[2].save("out/fig6-1-2c.png")
    for x, y in rangexy(w, h):
        for i in range(3):
            rgb = dst[i].getRgb(x, y)
            dst[i].setRgb(x, y, rgb[i], rgb[i], rgb[i])
    dst[0].save("out/fig6-1-3a.png")
    dst[1].save("out/fig6-1-3b.png")
    dst[2].save("out/fig6-1-3c.png")
    for i in range(3):
        dst[i] = test52x(dst[i], ERROR_FSTEIN, True, Gamma.Cto2)
    dst[0].save("out/fig6-1-4a.png")
    dst[1].save("out/fig6-1-4b.png")
    dst[2].save("out/fig6-1-4c.png")
    for x, y in rangexy(w, h):
        for i in range(3):
            rgb = [0., 0., 0.]
            rgb[i] = (dst[i].getRgb(x, y))[i]
            dst[i].setRgb(x, y, *rgb)
    dst[0].save("out/fig6-1-5a.png")
    dst[1].save("out/fig6-1-5b.png")
    dst[2].save("out/fig6-1-5c.png")
    for x, y in rangexy(w, h):
        rgb = [0., 0., 0.]
        for i in range(3):
            rgb[i] = (dst[i].getRgb(x, y))[i]
        tmp.setRgb(x, y, *rgb)
    tmp.save("out/fig6-1-6.png")

# Output 6.1.1: 8-colour Floyd-Steinberg RGB dithering
# Output 6.1.2: 8-colour gamma-corrected Floyd-Steinberg RGB dithering
def test61x(src, mat, func):
    (w, h) = src.size()
    dest = Image((w, h))
    tmp = [Image((w, h)), Image((w, h)), Image((w, h))]
    for x, y in rangexy(w, h):
        rgb = src.getRgb(x, y)
        for i in range(3):
            tmp[i].setGray(x, y, rgb[i])
    for i in range(3):
        tmp[i] = func(tmp[i], mat, True, Gamma.Cto2)
    for x, y in rangexy(w, h):
        (r, g, b) = [tmp[i].getGray(x, y) for i in range(3)]
        dest.setRgb(x, y, r, g, b)
    return dest

def test61y(src, mat, serpentine, threshold):
    return test3xx(src, mat, serpentine)

if chapter(6):
    test61x(grad256, ERROR_FSTEIN, test61y).save("out/grad6-1-1.png")
    test61x(lena256, ERROR_FSTEIN, test61y).save("out/lena6-1-1.png")
    test61x(grad256, ERROR_FSTEIN, test52x).save("out/grad6-1-2.png")
    test61x(lena256, ERROR_FSTEIN, test52x).save("out/lena6-1-2.png")

# Pattern 6.2.1: different colours give the same result
if chapter(6):
    dest = Image((320, 160))
    for x in range(80):
        for y in range(80):
            r = DITHER_BAYER44[(y / 8) % 4][(x / 8) % 4] > 7
            g = DITHER_BAYER44[(y / 8) % 4][(x / 8) % 4] > 13
            b = DITHER_BAYER44[(y / 8) % 4][(x / 8) % 4] > 13
            dest.setRgb(x, y, b, g, r)
        for y in range(80, 160):
            r = DITHER_BAYER44[y % 4][x % 4] > 7
            g = DITHER_BAYER44[y % 4][x % 4] > 13
            b = DITHER_BAYER44[y % 4][x % 4] > 13
            dest.setRgb(x, y, b, g, r)
    for x in range(80, 160):
        for y in range(80):
            r = DITHER_BAYER44[(y / 8) % 4][(x / 8 - 1) % 4] > 7
            g = DITHER_BAYER44[(y / 8) % 4][(x / 8) % 4] > 13
            b = DITHER_BAYER44[(y / 8) % 4][(x / 8) % 4] > 13
            dest.setRgb(x, y, b, g, r)
        for y in range(80, 160):
            r = DITHER_BAYER44[y % 4][(x - 1) % 4] > 7
            g = DITHER_BAYER44[y % 4][x % 4] > 13
            b = DITHER_BAYER44[y % 4][x % 4] > 13
            dest.setRgb(x, y, b, g, r)
    for x in range(160, 240):
        for y in range(80):
            r = DITHER_BAYER44[(y / 8 + 1) % 4][(x / 8 + 1) % 4] > 7
            g = DITHER_BAYER44[(y / 8) % 4][(x / 8) % 4] > 13
            b = DITHER_BAYER44[(y / 8 + 1) % 4][(x / 8) % 4] > 13
            dest.setRgb(x, y, b, g, r)
        for y in range(80, 160):
            r = DITHER_BAYER44[(y + 1) % 4][(x + 1) % 4] > 7
            g = DITHER_BAYER44[y % 4][x % 4] > 13
            b = DITHER_BAYER44[(y + 1) % 4][x % 4] > 13
            dest.setRgb(x, y, b, g, r)
    for x in range(240, 320):
        for y in range(80):
            r = DITHER_BAYER44[(y / 8 + 1) % 4][(x / 8) % 4] > 7
            g = DITHER_BAYER44[(y / 8) % 4][(x / 8) % 4] > 13
            b = DITHER_BAYER44[(y / 8) % 4][(x / 8 + 2) % 4] > 13
            dest.setRgb(x, y, b, g, r)
        for y in range(80, 160):
            r = DITHER_BAYER44[(y + 1) % 4][x % 4] > 7
            g = DITHER_BAYER44[y % 4][x % 4] > 13
            b = DITHER_BAYER44[y % 4][(x + 2) % 4] > 13
            dest.setRgb(x, y, b, g, r)
    dest.save("out/pat6-2-1.png")

# Pattern 6.2.2: 16-colour palette
if chapter(6):
    dest = Image((512, 128))
    for x, y in rangexy(64, 64):
        dest.setGray(x, y, 0)
        dest.setGray(448 + x, y, 0.7)
        dest.setGray(x, 64 + y, 0.3)
        dest.setGray(448 + x, 64 + y, 1.)
    for x in range(64, 448):
        d = x / 64
        r = (d & 2) >> 1
        g = (d & 4) >> 2
        b = d & 1
        for y in range(64, 128):
            dest.setRgb(x, y, r, g, b)
        r /= 2.
        g /= 2.
        b /= 2.
        for y in range(64):
            dest.setRgb(x, y, r, g, b)
    dest.save("out/pat6-2-2.png")

# Output 6.2.1: gamma-corrected Floyd-Steinberg with ANSI palette (sigma-abs)
# Output 6.2.2: gamma-corrected Floyd-Steinberg with ANSI palette (euclidian)
def test62x(src, mat, cpal, distance, serpentine):
    (w, h) = src.size()
    ipal = [[Gamma.CtoI(c[i]) for i in range(3)] for c in cpal]
    dest = Image((w, h))
    lines = len(mat)
    rows = len(mat[0])
    offset = mat[0].index(-1)
    ey = [[[0.] * 3 for n in range(w + rows - 1)] for n in range(lines)]
    for y in range(h):
        ex = [[0.] * 3 for n in range(rows - offset)]
        if serpentine and y & 1:
            xrange = range(w - 1, -1, -1)
        else:
            xrange = range(w)
        for x in xrange:
            # Get pixel, add error, set pixel
            crgb = src.getRgb(x, y)
            irgb = [Gamma.CtoI(crgb[i]) + ex[0][i] + ey[0][x + offset][i] \
                        for i in range(3)]
            crgb = [Gamma.ItoC(irgb[i]) for i in range(3)]
            max = 999999.
            for n in range(len(cpal)):
                d = distance(crgb, cpal[n])
                if d < max:
                    max = d
                    cbest = cpal[n]
                    ibest = ipal[n]
            dest.setRgb(x, y, *cbest)
            # Compute error and propagate it
            for i in range(3):
                # Propagate first line of error
                error = irgb[i] - ibest[i]
                for dx in range(rows - offset - 2):
                    ex[dx][i] = ex[dx + 1][i] + error * mat[0][offset + 1 + dx]
                ex[rows - offset - 2][i] = error * mat[0][rows - 1]
                # Propagate next lines
                if serpentine and y & 1:
                    for dy in range(1, lines):
                        for dx in range(rows):
                            ey[dy][x + dx][i] += error * mat[dy][rows - 1 - dx]
                else:
                    for dy in range(1, lines):
                        for dx in range(rows):
                            ey[dy][x + dx][i] += error * mat[dy][dx]
        for dy in range(lines - 1):
            ey[dy] = ey[dy + 1]
        ey[lines - 1] = [[0.] * 3 for n in range(w + rows - 1)]
    return dest

RGB_PALETTE = \
    [[0, 0, 0],
     [0, 0, 1],
     [0, 1, 0],
     [0, 1, 1],
     [1, 0, 0],
     [1, 0, 1],
     [1, 1, 0],
     [1, 1, 1]]

ANSI_PALETTE = \
    [[0, 0, 0],
     [0, 0, 0.5],
     [0, 0.5, 0],
     [0, 0.5, 0.5],
     [0.5, 0, 0],
     [0.5, 0, 0.5],
     [0.5, 0.5, 0],
     [0.7, 0.7, 0.7],
     [0.3, 0.3, 0.3],
     [0, 0, 1],
     [0, 1, 0],
     [0, 1, 1],
     [1, 0, 0],
     [1, 0, 1],
     [1, 1, 0],
     [1, 1, 1]]

def distmax(u, v):
    r, g, b = [abs(u[i] - v[i]) for i in range(3)]
    return r + g + b

def disteuclidian(u, v):
    r, g, b = [u[i] - v[i] for i in range(3)]
    return r*r + g*g + b*b

if chapter(6):
    test62x(grad256, ERROR_FSTEIN,
            ANSI_PALETTE, distmax, True).save("out/grad6-2-1.png")
    test62x(lena256, ERROR_FSTEIN,
            ANSI_PALETTE, distmax, True).save("out/lena6-2-1.png")
    test62x(grad256, ERROR_FSTEIN,
            ANSI_PALETTE, disteuclidian, True).save("out/grad6-2-2.png")
    test62x(lena256, ERROR_FSTEIN,
            ANSI_PALETTE, disteuclidian, True).save("out/lena6-2-2.png")

def rgb2hsv(r, g, b):
    m = (float)(min(r, g, b))
    v = (float)(max(r, g, b))
    if v == m or v == 0:
        return 0., 0., v
    s = (v - m) / v
    if r == v:
        h = (g - b) / (v - m)
        if h < 0.:
            h += 6
    elif g == v:
        h = 2. + (b - r) / (v - m)
    elif b == v:
        h = 4. + (r - g) / (v - m)
    return math.pi * h / 3, s, v

def disthsv(u, v):
    u = rgb2hsv(*u)
    v = rgb2hsv(*v)
    d1 = u[2] - v[2]
    d2 = u[2] * u[1] * math.cos(u[0]) - v[2] * v[1] * math.cos(v[0])
    d3 = u[2] * u[1] * math.sin(u[0]) - v[2] * v[1] * math.sin(v[0])
    return d1*d1 + d2*d2 + d3*d3

def disthsv3(u, v):
    u = rgb2hsv(*u)
    v = rgb2hsv(*v)
    d1 = u[2] - v[2]
    d2 = u[2] * u[1] * math.cos(u[0]) - v[2] * v[1] * math.cos(v[0])
    d3 = u[2] * u[1] * math.sin(u[0]) - v[2] * v[1] * math.sin(v[0])
    return 9*d1*d1 + d2*d2 + d3*d3

if chapter(6):
    test62x(grad256, ERROR_FSTEIN,
            ANSI_PALETTE, disthsv, True).save("out/grad6-2-3.png")
    test62x(lena256, ERROR_FSTEIN,
            ANSI_PALETTE, disthsv, True).save("out/lena6-2-3.png")
    test62x(grad256, ERROR_FSTEIN,
            ANSI_PALETTE, disthsv3, True).save("out/grad6-2-4.png")
    test62x(lena256, ERROR_FSTEIN,
            ANSI_PALETTE, disthsv3, True).save("out/lena6-2-4.png")

# Output 6.4.1: colour sub-block error diffusion
def colorsubblock(src, tiles, propagate, diff):
    (w, h) = src.size()
    # Propagating the error to a temporary buffer is becoming more and
    # more complicated. We decide to use an intermediate matrix instead.
    tmp = Matrix(w, h, None)
    for x, y in rangexy(w, h):
        tmp[y][x] = Gamma.CtoI3(src.getRgb(x, y))
    dest = Image((w, h))
    # Analyse tile list
    ntiles = len(tiles)
    ty = len(tiles[0])
    tx = len(tiles[0][0])
    cur = Matrix(tx, ty, None)
    w, h = w / tx, h / ty
    # Analyse error propagate list
    for x, y in rangexy(w, h):
        # Get block value
        for i, j in rangexy(tx, ty):
            cur[j][i] = Gamma.ItoC3(tmp[y * ty + j][x * tx + i])
        # Select closest block
        dist = tx * ty
        for n in range(ntiles):
            d = [0., 0., 0.]
            e = 0.
            for i, j in rangexy(tx, ty):
                tmpe = 0.
                for k in range(3):
                    delta = cur[j][i][k] - tiles[n][j][i][k]
                    d[k] += delta
                    tmpe += delta * delta
                e += diff[j][i] * math.sqrt(tmpe)
            # Without / 3. ugly colour bleeding artifacts appear
            absd = (abs(d[0]) + abs(d[1]) + abs(d[2])) / 3.
            if absd / (tx * ty) + e < dist:
                dist = absd / (tx * ty) + e
                best = n
        # Set pixel
        for i, j in rangexy(tx, ty):
            dest.setRgb(x * tx + i, y * ty + j, *(tiles[best][j][i]))
        # Propagate error
        for i, j in rangexy(tx, ty):
            curp = Gamma.CtoI3(cur[j][i])
            bestp = Gamma.CtoI3(tiles[best][j][i])
            m = propagate[j][i]
            for k in range(3):
                e = curp[k] - bestp[k]
                for px, py in rangexy(len(m[0]), len(m)):
                    if m[py][px] == 0:
                        continue
                    if m[py][px] == -1:
                        cx, cy = px, py
                        continue
                    tmpx = x * tx + i + px - cx
                    tmpy = y * ty + j + py - cy
                    if tmpx > w * tx - 1 or tmpy > h * ty - 1:
                        continue
                    tmp[tmpy][tmpx][k] += m[py][px] * e
    return dest

RGBLINES22 = []
for n in range(8*8*8*8):
    a, b, c, d = n & 7, (n >> 3) & 7, (n >> 6) & 7, (n >> 9) & 7
    if (a != b or c != d) and (a != c or b != d):
        continue
    RGBLINES22.append([[RGB_PALETTE[a], RGB_PALETTE[b]],
                       [RGB_PALETTE[c], RGB_PALETTE[d]]])

if chapter(6):
    colorsubblock(grad256, RGBLINES22,
                  ERROR_SUBFS22, DIFF_WEIGHTED22).save("out/grad6-4-1.png")
    colorsubblock(lena256, RGBLINES22,
                  ERROR_SUBFS22, DIFF_WEIGHTED22).save("out/lena6-4-1.png")

ANSILINES22 = []
for n in range(16*16*16*16):
    a, b, c, d = n & 15, (n >> 4) & 15, (n >> 8) & 15, (n >> 12) & 15
    if (a != b or c != d) and (a != c or b != d):
        continue
    ANSILINES22.append([[ANSI_PALETTE[a], ANSI_PALETTE[b]],
                        [ANSI_PALETTE[c], ANSI_PALETTE[d]]])

if chapter(6):
    colorsubblock(grad256, ANSILINES22,
                  ERROR_SUBFS22, DIFF_WEIGHTED22).save("out/grad6-4-2.png")
    colorsubblock(lena256, ANSILINES22,
                  ERROR_SUBFS22, DIFF_WEIGHTED22).save("out/lena6-4-2.png")

##############################################################################
if chapter(7):
    print "Chapter 7. Photographic mosaics"

# Output 7.0.1: create a mosaic from Lena
def mosaic_split(src, tnw, tnh):
    random.seed(0)
    thumbs = []
    (w, h) = src.size()
    sw = w / tnw
    sh = h / tnh
    for x, y in rangexy(sw, sh):
        thumbs.append(src.getRegion(x * tnw, y * tnh, tnw, tnh))
    random.shuffle(thumbs)
    return thumbs

def mosaic_analyse(tnlist, dx, dy):
    coeffs = []
    for (n, img) in enumerate(tnlist):
        tmp = [[[0] * 3 for x in range(dx)] for y in range(dy)]
        (w, h) = img.size()
        for x, y in rangexy(w, h):
            my = y * dy / h
            mx = x * dx / w
            (r, g, b) = img.getRgb(x, y)
            tmp[my][mx][0] += Gamma.CtoI(r) / (w / dx * h / dy)
            tmp[my][mx][1] += Gamma.CtoI(g) / (w / dx * h / dy)
            tmp[my][mx][2] += Gamma.CtoI(b) / (w / dx * h / dy)
        coeffs.append(tmp)
    return coeffs

def test701(tnlist, cols):
    (tnw, tnh) = tnlist[0].size()
    dw = cols
    dh = (len(tnlist) + cols - 1) / cols
    dest = Image((dw * tnw + 8 * (dw + 1), dh * tnh + 8 * (dh + 1)), True)
    for (n, img) in enumerate(tnlist):
        di = 8 + (n % dw) * (tnw + 8)
        dj = 8 + (n / dw) * (tnh + 8)
        img.copyTo(dest, (di, dj))
    return dest

if chapter(7):
    tnlist = mosaic_split(lena256, 32, 32)
    test701(tnlist, 10).save("out/lena7-0-1.png")

# Output 7.1.1: extract 1 colour feature from mosaic tiles
# Output 7.1.2: crop Lena
# Output 7.1.3: generate a mosaic from the 1-feature database
# Output 7.1.4: extract 4 colour features from mosaic tiles
# Output 7.1.5: generate a mosaic from the 4-feature database
def test71x(coeffs, cols, tnw, tnh):
    dx = len(coeffs[0][0])
    dy = len(coeffs[0])
    dw = cols
    dh = (len(coeffs) + cols - 1) / cols
    dest = Image((dw * tnw + 8 * (dw + 1), dh * tnh + 8 * (dh + 1)), True)
    for (n, tab) in enumerate(coeffs):
        di = 8 + (n % dw) * (tnw + 8)
        dj = 8 + (n / dw) * (tnh + 8)
        for x, y in rangexy(tnw, tnh):
            (r, g, b) = tab[y * dy / tnh][x * dx / tnw]
            dest.setRgb(di + x, dj + y, Gamma.ItoC(r), Gamma.ItoC(g), Gamma.ItoC(b))
    return dest

def test71y(src, sqw, sqh, tnlist, coeffs):
    (w, h) = src.size()
    (tnw, tnh) = tnlist[0].size()
    dx = len(coeffs[0][0])
    dy = len(coeffs[0])
    nx = w / sqw
    ny = h / sqh
    dest = Image((nx * tnw, ny * tnh), True)
    for X, Y in rangexy(nx, ny):
        # 1. create statistics about the current square
        cur = [[[0] * 3 for x in range(dx)] for y in range(dy)]
        for x, y in rangexy(sqw, sqh):
            my = y * dy / sqh
            mx = x * dx / sqw
            (r, g, b) = src.getRgb(X * sqw + x, Y * sqh + y)
            cur[my][mx][0] += Gamma.CtoI(r) / (sqw / dx * sqh / dy)
            cur[my][mx][1] += Gamma.CtoI(g) / (sqw / dx * sqh / dy)
            cur[my][mx][2] += Gamma.CtoI(b) / (sqw / dx * sqh / dy)
        # 2. find the best mosaic part
        best = -1
        dist = 5.
        for (n, tmp) in enumerate(coeffs):
            d = 0
            for i, j in rangexy(dx, dy):
                for c in range(3):
                    t = cur[j][i][c] - tmp[j][i][c]
                    d += t * t
            if d < dist:
                dist = d
                best = n
        # 3. blit mosaic chunk
        tnlist[best].copyTo(dest, (X * tnw, Y * tnh))
    return dest

if chapter(7):
    coeffs1x1 = mosaic_analyse(tnlist, 1, 1)
    test71x(coeffs1x1, 10, 8, 8).save("out/lena7-1-1.png")
    out712 = lena256.getRegion(100, 90, 80, 80)
    out712.save("out/lena7-1-2.png")
    test71y(out712, 6, 6, tnlist, coeffs1x1).save("out/lena7-1-3.png")

    coeffs2x2 = mosaic_analyse(tnlist, 2, 2)
    test71x(coeffs2x2, 10, 16, 16).save("out/lena7-1-4.png")
    test71y(out712, 6, 6, tnlist, coeffs2x2).save("out/lena7-1-5.png")

##############################################################################
print "Finished"

# Place temporary cruft below this
sys.exit(0)

def simusubblock(src):
    (w, h) = src.size()
    dest = Image((w, h))
    tmp = Matrix(w, h, 0.)
    for x, y in rangexy(w, h):
        tmp[y][x] = src.getRgb(x, y)
    dest = Image((w, h))
    # Analyse tile list
    ntiles = len(tiles)
    ty = len(tiles[0])
    tx = len(tiles[0][0])
    cur = Matrix(tx, ty, None)
    w, h = w / tx, h / ty
    # Analyse error propagate list
    for x, y in rangexy(w, h):
        # Get block value
        for i, j in rangexy(tx, ty):
            cur[j][i] = Gamma.ItoC3(tmp[y * ty + j][x * tx + i])
        # Select closest block
        dist = tx * ty
        for n in range(ntiles):
            d = [0., 0., 0.]
            e = 0.
            for i, j in rangexy(tx, ty):
                tmpe = 0.
                for k in range(3):
                    delta = cur[j][i][k] - tiles[n][j][i][k]
                    d[k] += delta
                    tmpe += delta * delta
                e += diff[j][i] * math.sqrt(tmpe)
            # Without / 3. ugly colour bleeding artifacts appear
            absd = (abs(d[0]) + abs(d[1]) + abs(d[2])) / 3.
            if absd / (tx * ty) + e < dist:
                dist = absd / (tx * ty) + e
                best = n
        # Set pixel
        for i, j in rangexy(tx, ty):
            dest.setRgb(x * tx + i, y * ty + j, *(tiles[best][j][i]))
        # Propagate error
        for i, j in rangexy(tx, ty):
            curp = Gamma.CtoI3(cur[j][i])
            bestp = Gamma.CtoI3(tiles[best][j][i])
            m = propagate[j][i]
            for k in range(3):
                e = curp[k] - bestp[k]
                for px, py in rangexy(len(m[0]), len(m)):
                    if m[py][px] == 0:
                        continue
                    if m[py][px] == -1:
                        cx, cy = px, py
                        continue
                    tmpx = x * tx + i + px - cx
                    tmpy = y * ty + j + py - cy
                    if tmpx > w * tx - 1 or tmpy > h * ty - 1:
                        continue
                    tmp[tmpy][tmpx][k] += m[py][px] * e
    return dest


# XXX: test -- ranked dither -- it SUCKS
def test26x(src, mat):
    (w, h) = src.size()
    dest = Image((w, h))
    dx = len(mat[0])
    dy = len(mat)
    for x, y in rangexy(w / dx, h / dy):
        # Step 1: get the pixels and count groups
        groups = {}
        for i, j in rangexy(dx, dy):
            p = src.getGray(x * dx + i, y * dy + j)
            if groups.has_key(p):
                groups[p].append((i, j))
            else:
                groups[p] = [(i, j)]
        # Step 2: create the ranked dither
        ranked = Matrix(dx, dy)
        for p, g in groups.items():
            n = (int)(round(p * len(g)))
            if not n:
                continue
            v = [(mat[j][i], (i, j)) for (i, j) in g]
            v.sort()
            v = v[0 : n - 1]
            for (k, (i, j)) in v:
                ranked[j][i] = 1
        # Step 3: blit the ranked dither
        for j in range(dy):
            for i in range(dx):
                dest.setGray(x * dx + i, y * dy + j, ranked[j][i])
    return dest

if chapter(2):
    #test26x(lena256bw, DITHER_BAYER88).save("out/lena2-6-1.png")
    #test26x(grad256bw, DITHER_BAYER88).save("out/grad2-6-1.png")
    test26x(lena256bw, DITHER_CLUSTER88).save("out/lena2-6-1.png")
    test26x(grad256bw, DITHER_CLUSTER88).save("out/grad2-6-1.png")

#####################
#CIE-L*a*b* transformation -- euclidian distance doesn't seem to work great
def rgb2lab(r, g, b):
    Xw50 = 0.9642;
    Yw50 = 1.0000;
    Zw50 = 0.8249;
    # RGB to sRGB
    r = Gamma.CtoI(r)
    g = Gamma.CtoI(g)
    b = Gamma.CtoI(b)
    # sRGB to XYZ(D65)
    x65 =  0.4124*r + 0.3576*g + 0.1805*b
    y65 =  0.2126*r + 0.7152*g + 0.0722*b
    z65 =  0.0193*r + 0.1192*g + 0.9505*b
    # XYZ(D65) to XYZ(D50)
    x50 =  1.0282015*x65 + 0.0500707*y65 - 0.0579688*z65
    y50 =  0.0197032*x65 + 0.9871848*y65 - 0.0054285*z65
    z50 = -0.0002329*x65 + 0.0006862*y65 + 0.7573070*z65
    # XYZ(D50) to Lab(D50)
    if x50 / Xw50 > 0.008856:
        xx50 = math.pow(x50 / Xw50, 1. / 3.)
    else:
        xx50 = 7.78 * (x50 / Xw50) + 16. / 116.
    if y50 / Yw50 > 0.008856:
        yy50 = math.pow(y50 / Yw50, 1. / 3.)
    else:
        yy50 = 7.78 * (y50 / Yw50) + 16. / 116.
    if z50 / Zw50 > 0.008856:
        zz50 = math.pow(z50 / Zw50, 1. / 3.)
    else:
        zz50 = 7.78 * (z50 / Zw50) + 16. / 116.
    lo = 116. * yy50 - 16.
    ao = 500. * (xx50 - yy50)
    bo = 200. * (yy50 - zz50)
    return lo, ao, bo

def distlab(u, v):
    u = rgb2lab(*u)
    v = rgb2lab(*v)
    l, a, b = u[0] - v[0], u[1] - v[1], u[2] - v[2]
    return l*l + a*a + b*b