#!/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 = \ '\n' \ '\n' \ ' \n' \ ' \n' % (64 * w + 2, 64 * h + 2) line = \ ' \n' box = \ ' \n' text = \ ' %s\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 += \ ' \n' \ '\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, 'd1(i)'], ['d2(i)', 'd3(i)', 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