# HG changeset patch
# User pmeerw@pan
# Date 1210107688 -7200
# Node ID 63af49cca5d27e099caee4b11ab2feeb90f16c76

initial import

diff -r 000000000000 -r 63af49cca5d2 .hgignore
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Tue May 06 23:01:28 2008 +0200
@@ -0,0 +1,23 @@
+syntax: glob
+
+*.elc
+*.orig
+*.rej
+*~
+*.mergebackup
+*.o
+*.so
+*.pyc
+*.swp
+*.prof
+tests/.coverage*
+tests/annotated
+tests/*.err
+build
+MANIFEST
+.DS_Store
+tags
+cscope.*
+
+syntax: regexp
+^\.pc/
diff -r 000000000000 -r 63af49cca5d2 README
diff -r 000000000000 -r 63af49cca5d2 pydct/__init__.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pydct/__init__.py	Tue May 06 23:01:28 2008 +0200
@@ -0,0 +1,6 @@
+from dct import fdct, idct, zigzag
+from t4x4 import fdct4x4, idct4x4
+from t8x8 import fdct8x8, idct8x8, fdct8x8x8, idct8x8x8, do_f8x8, do_i8x8, build_vector_from_8x8, set_8x8_from_vector
+from tNxM import fdct2d, idct2d
+from tNxN import fdctNxN, idctNxN
+from perceptual import watson_slack
diff -r 000000000000 -r 63af49cca5d2 pydct/dct.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pydct/dct.py	Tue May 06 23:01:28 2008 +0200
@@ -0,0 +1,54 @@
+import numpy, math
+
+__dctN = {}
+
+def __init_dctN(n):
+    global __dctN
+    if not __dctN.has_key(n): 
+        c = numpy.empty((n, n), numpy.float)
+        c[0,:] = 1.0 / math.sqrt(n)
+        __cosf = lambda i, j: numpy.cos(math.pi * (2*j + 1) * (i + 1) / (2.0 * n))
+        c[1:,:] = math.sqrt(2.0 / n) * numpy.fromfunction(__cosf, (n-1, n))
+        __dctN[n] = c
+    return __dctN[n]
+
+def fdct(v):
+    '''
+    Forward DCT on vector.
+    '''
+    c = __init_dctN(len(v))
+    return numpy.dot(c, v)
+
+def idct(v):
+    '''
+    Inverse DCT on vector.
+    '''
+    c = __init_dctN(len(v))
+    return numpy.dot(numpy.transpose(c), v)
+
+def zigzag(x):
+    """Generates zig-zag scan sequence for two dimensional array."""
+    j = 0
+    i = 0
+    while True:
+
+        if i == x.shape[0]: j += 2; i = x.shape[0]-1
+        else: j = 0
+        while i >= 0 and j < x.shape[1]:
+            # run up
+            yield i, j
+            i -= 1
+            j += 1
+
+        if j > x.shape[1]: break
+
+        if j == x.shape[1]: j = x.shape[1]-1; i += 2
+        else: i = 0
+        while j >= 0 and i < x.shape[0]:
+            # run down
+            yield i, j
+            i += 1
+            j -= 1
+
+        if i > x.shape[0]: break
+
diff -r 000000000000 -r 63af49cca5d2 pydct/perceptual.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pydct/perceptual.py	Tue May 06 23:01:28 2008 +0200
@@ -0,0 +1,39 @@
+import numpy
+
+JPEG_luma_Q = numpy.array(
+    [[16,  11,  10,  16,  24,  40,  51,  61],
+    [12,  12,  14,  19,  26,  58,  60,  55],
+    [14,  13,  16,  24,  40,  57,  69,  56],
+    [14,  17,  22,  29,  51,  87,  80,  62],
+    [18,  22,  37,  56,  68, 109, 103,  77],
+    [24,  35,  55,  64,  81, 104, 113,  92],
+    [49,  64,  78,  87, 103, 121, 120, 101],
+    [72,  92,  95,  98, 112, 100, 103,  99]])
+
+JPEG_chroma_Q = numpy.array(
+    [[17,  18,  24,  47,  99,  99,  99,  99],
+    [18,  21,  26,  66,  99,  99,  99,  99],
+    [24,  26,  56,  99,  99,  99,  99,  99],
+    [47,  66,  99,  99,  99,  99,  99,  99],
+    [99,  99,  99,  99,  99,  99,  99,  99],
+    [99,  99,  99,  99,  99,  99,  99,  99],
+    [99,  99,  99,  99,  99,  99,  99,  99],
+    [99,  99,  99,  99,  99,  99,  99,  99]])
+
+def watson_slack(b, c00):
+    # table of sensitivity values
+    t = numpy.array( \
+        [1.404,  1.011,  1.169,  1.664,  2.408,  3.433,  4.796,  6.563,\
+         1.011,  1.452,  1.323,  1.529,  2.006,  2.716,  3.679,  4.939,\
+         1.169,  1.323,  2.241,  2.594,  2.988,  3.649,  4.604,  5.883,\
+         1.664,  1.529,  2.594,  3.773,  4.559,  5.305,  6.281,  7.600,\
+         2.408,  2.006,  2.988,  4.559,  6.152,  7.463,  8.713, 10.175,\
+         3.433,  2.716,  3.649,  5.305,  7.463,  9.625, 11.588, 13.519,\
+         4.796,  3.679,  4.604,  6.281,  8.713, 11.588, 14.500, 17.294,\
+         6.563,  4.939,  5.883,  7.600, 10.175, 13.519, 17.294, 21.156])
+
+    tl = t * math.pow(b.ravel()[0] / c00, 0.649)
+  
+    tl[1:] = numpy.maximum(tl[1:], (numpy.abs(b.ravel()[1:])**0.7) * (tl[1:]**0.3))
+    return tl.reshape((8, 8))
+
diff -r 000000000000 -r 63af49cca5d2 pydct/t4x4.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pydct/t4x4.py	Tue May 06 23:01:28 2008 +0200
@@ -0,0 +1,39 @@
+import numpy
+
+def fdct4x4(b):
+    '''
+    Compute the approximate 4x4 DCT coefficients of an array 
+    as defined by H.264/AVC.
+    '''
+    c, e = \
+        numpy.array(
+            [[ 1,  1,  1,  1], 
+             [ 2,  1, -1, -2], 
+             [ 1, -1, -1,  1], 
+             [ 1, -2,  2, -1]], numpy.float), \
+        numpy.array(
+            [[0.25, 0.15811388, 0.25, 0.15811388], 
+             [0.15811388, 0.1, 0.15811388, 0.1],
+             [0.25, 0.15811388, 0.25, 0.15811388], 
+             [0.15811388, 0.1, 0.15811388, 0.1]], numpy.float)
+        
+    return numpy.dot(c, numpy.dot(b, numpy.transpose(c))) * e
+
+def idct4x4(b):
+    '''
+    Compute the inverse 4x4 DCT of the array.
+    '''
+    c, e = \
+        numpy.array(
+            [[1, 1, 1, 0.5], 
+             [1, 0.5, -1, -1], 
+             [1, -0.5, -1, 1], 
+             [1, -1, 1, -0.5]], numpy.float), \
+        numpy.array(
+            [[0.25, 0.31622777, 0.25, 0.31622777], 
+             [0.31622777, 0.4, 0.31622777, 0.4], 
+             [0.25, 0.31622777, 0.25, 0.31622777],
+             [0.31622777, 0.4, 0.31622777, 0.4]], numpy.float)
+
+    return numpy.dot(c, numpy.dot(b*e, numpy.transpose(c)))
+
diff -r 000000000000 -r 63af49cca5d2 pydct/t8x8.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pydct/t8x8.py	Tue May 06 23:01:28 2008 +0200
@@ -0,0 +1,75 @@
+import tNxN, numpy
+
+def fdct8x8(b):
+    '''
+    Forward 8x8 DCT.
+    '''
+    return fdctNxN(b)    
+    
+def idct8x8(b):
+    '''
+    Inverse 8x8 DCT.
+    '''
+    return idctNxN(b)
+
+def fdct8x8x8(b):
+    '''
+    Forward 8x8x8 DCT.
+    '''
+    d = numpy.empty((8, 8, 8), numpy.float)
+    for i in xrange(8):
+        d[i,:,:] = fdct8x8(b[i,:,:])    
+    for r in xrange(8):
+        for c in xrange(8):
+            d[:,r,c] = fdct(d[:,r,c])            
+    return d
+    
+def idct8x8x8(b):
+    '''
+    Inverse 8x8x8 DCT.
+    '''
+    d = numpy.empty((8, 8, 8), numpy.float)
+    for r in xrange(8):
+        for c in xrange(8):
+            d[:,r,c] = idct(b[:,r,c])            
+    for i in xrange(8):
+        d[i,:,:] = idct8x8(d[i,:,:])    
+    return d
+
+def do_f8x8(s):
+    for r in xrange(0, s.shape[0], 8):
+        for c in xrange(0, s.shape[1], 8): 
+            s[r:r+8,c:c+8] = fdct8x8(s[r:r+8,c:c+8])
+    return s
+
+def _get_zz():
+    zz = []
+    for i, j in wm.zigzag(numpy.empty((8,8))):
+        zz.append(i*8 + j)
+    return zz
+
+def build_vector_from_8x8(s):
+    v = numpy.zeros((64, (s.shape[0]//8)*(s.shape[1]//8)), numpy.float)
+    i = 0
+    zz = _get_zz()
+    for r in xrange(0, s.shape[0], 8):
+        for c in xrange(0, s.shape[1], 8): 
+            v[:,i] = s[r:r+8,c:c+8].flat[zz]
+            i += 1
+    return v
+
+def do_i8x8(s):
+    for r in xrange(0, s.shape[0], 8):
+        for c in xrange(0, s.shape[1], 8): 
+            s[r:r+8,c:c+8] = idct8x8(s[r:r+8,c:c+8])
+    return s
+
+def set_8x8_from_vector(s, v):
+    i = 0
+    zz = _get_zz()
+    for r in xrange(0, s.shape[0], 8):
+        for c in xrange(0, s.shape[1], 8): 
+            s[r:r+8,c:c+8].flat[zz] = v[:,i]
+            i += 1
+    return s            
+
diff -r 000000000000 -r 63af49cca5d2 pydct/tNxM.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pydct/tNxM.py	Tue May 06 23:01:28 2008 +0200
@@ -0,0 +1,56 @@
+import numpy, math
+
+__dctNxM = {}
+
+def __init_dctNxM(b):
+    global __dctNxM
+    if not __dctNxM.has_key(b.shape): 
+        c = numpy.sqrt([2.0 / b.shape[0], 2.0 / b.shape[1]])
+        __cosf0 = lambda j, i: numpy.cos((math.pi * (2*i + 1) * j) / (2.0 * b.shape[0]))
+        y = c[0] * numpy.fromfunction(__cosf0, (b.shape[0], b.shape[0]))
+        __cosf1 = lambda j, i: numpy.cos((math.pi * (2*i + 1) * j) / (2.0 * b.shape[1]))
+        x = c[1] * numpy.fromfunction(__cosf1, (b.shape[1], b.shape[1]))
+        __dctNxM[b.shape] = (y, x), c / math.sqrt(2.0)
+    return __dctNxM[b.shape]
+    
+def fdct2d(b):
+    '''
+    Compute forward 2D DCT on the array of arbitrary size.
+    Note: slow.
+    '''
+    global __dctNxM
+    ct, c = __init_dctNxM(b)
+    
+    d = numpy.empty(b.shape, numpy.float)
+    d[0][0] = c[0] * c[1] * numpy.sum(numpy.sum(b))
+
+    for i in xrange(1, b.shape[1]):
+        d[0][i] = c[0] * numpy.sum(ct[1][i] * numpy.sum(b, 0))
+
+    for j in xrange(1, b.shape[0]):
+        d[j][0] = c[1] * numpy.sum(ct[0][j] * numpy.sum(b, 1))
+
+    bt = numpy.transpose(b)
+    for j in xrange(1, b.shape[0]):
+        d[j,1:] = numpy.sum(ct[1][1:] * numpy.sum(bt * ct[0][j], 1), 1)
+    
+    return d
+
+def idct2d(b):
+    '''
+    Compute inverse 2D DCT on array of arbitrary size.
+    '''
+    global __dctNxM
+    ct, c = __init_dctNxM(b)
+
+    p0 = c[0] * c[1] * b[0][0]
+    bt = numpy.transpose(b[1:,1:])
+    ctt = ct[1][1:,:].transpose()
+    d = numpy.empty(b.shape, numpy.float)
+    for y in xrange(b.shape[0]):
+        d[y,:] = p0 + \
+            c[0] * numpy.sum(b[0,1:] * ctt, 1) + \
+            c[1] * numpy.sum(b[1:,0] * ct[0][1:,y]) + \
+            numpy.sum(numpy.sum(bt * ct[0][1:,y],1) * ctt,1)
+    return d
+
diff -r 000000000000 -r 63af49cca5d2 pydct/tNxN.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pydct/tNxN.py	Tue May 06 23:01:28 2008 +0200
@@ -0,0 +1,15 @@
+import dct, numpy
+
+def fdctNxN(b):
+    '''
+    Forward 2D DCT on NxN array.
+    '''
+    c = __init_dctN(b.shape[0])
+    return numpy.dot(c, numpy.dot(b, numpy.transpose(c)))
+
+def idctNxN(b):
+    '''
+    Inverse 2D DCT on NxN array.
+    '''
+    c = __init_dctN(b.shape[0])
+    return numpy.dot(numpy.transpose(c), numpy.dot(b, c))
diff -r 000000000000 -r 63af49cca5d2 setup.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/setup.py	Tue May 06 23:01:28 2008 +0200
@@ -0,0 +1,10 @@
+from distutils.core import setup
+
+setup (name = 'pydct',
+       version = '1.0',
+       author = 'Peter Meerwald',
+       author_email = 'pmeerw@pmeerw.net',
+       description = 'pydct: Discrete Cosine Transform routines for Python',
+       packages=['pydct'],
+       license = 'BSD'
+)