← Back to branch summary

~ubuntu-branches/ubuntu/lucid/python-scipy/lucid

« back to all changes in this revision

Viewing changes to Lib/sandbox/image/transforms.py

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/image/transforms.py
 
1
import numpy as sb
 
2
pi = sb.pi
 
3
 
 
4
# fast discrete cosine transforms of real sequences (using the fft)
 
5
#  These implement the DCT-II and inverse DCT-II (DCT-III)
 
6
#  described at http://en.wikipedia.org/wiki/Discrete_cosine_transform
 
7
 
 
8
def dct(x,axis=-1):
 
9
    """Discrete cosine transform based on the FFT.
 
10
 
 
11
    For even-length signals it uses an N-point FFT
 
12
    For odd-length signals it uses a 2N-point FFT.
 
13
    """
 
14
    n = len(x.shape)
 
15
    N = x.shape[axis]
 
16
    even = (N%2 == 0)
 
17
    slices = [None]*4
 
18
    for k in range(4):
 
19
        slices[k] = []
 
20
        for j in range(n):
 
21
            slices[k].append(slice(None))
 
22
    if even:
 
23
        xtilde = 0.0*x
 
24
        slices[0][axis] = slice(None,N/2)
 
25
        slices[1][axis] = slice(None,None,2)
 
26
        slices[2][axis] = slice(N/2,None)
 
27
        slices[3][axis] = slice(N,None,-2)
 
28
    else:
 
29
        newshape = list(x.shape)
 
30
        newshape[axis] = 2*N
 
31
        xtilde = sb.empty(newshape,sb.float64)
 
32
        slices[0][axis] = slice(None,N)
 
33
        slices[2][axis] = slice(N,None)
 
34
        slices[3][axis] = slice(None,None,-1)
 
35
 
 
36
    for k in range(4):
 
37
        slices[k] = tuple(slices[k])
 
38
    xtilde[slices[0]] = x[slices[1]]
 
39
    xtilde[slices[2]] = x[slices[3]]
 
40
    Xt = sb.fft(xtilde,axis=axis)
 
41
    pk = sb.exp(-1j*pi*sb.arange(N)/(2*N))
 
42
    newshape = sb.ones(n)
 
43
    newshape[axis] = N
 
44
    pk.shape = newshape
 
45
 
 
46
    if not even:
 
47
        pk /= 2;
 
48
        Xt = Xt[slices[0]]
 
49
 
 
50
    return sb.real(Xt*pk)
 
51
 
 
52
 
 
53
def idct(v,axis=-1):
 
54
    n = len(v.shape)
 
55
    N = v.shape[axis]
 
56
    even = (N%2 == 0)
 
57
    slices = [None]*4
 
58
    for k in range(4):
 
59
        slices[k] = []
 
60
        for j in range(n):
 
61
            slices[k].append(slice(None))
 
62
    k = sb.arange(N)
 
63
    if even:
 
64
        ak = sb.r_[1.0,[2]*(N-1)]*sb.exp(1j*pi*k/(2*N))
 
65
        newshape = sb.ones(n)
 
66
        newshape[axis] = N
 
67
        ak.shape = newshape
 
68
        xhat = sb.real(sb.ifft(v*ak,axis=axis))
 
69
        x = 0.0*v
 
70
        slices[0][axis] = slice(None,None,2)
 
71
        slices[1][axis] = slice(None,N/2)
 
72
        slices[2][axis] = slice(N,None,-2)
 
73
        slices[3][axis] = slice(N/2,None)
 
74
        for k in range(4):
 
75
            slices[k] = tuple(slices[k])
 
76
        x[slices[0]] = xhat[slices[1]]
 
77
        x[slices[2]] = xhat[slices[3]]
 
78
        return x
 
79
    else:
 
80
        ak = 2*sb.exp(1j*pi*k/(2*N))
 
81
        newshape = sb.ones(n)
 
82
        newshape[axis] = N
 
83
        ak.shape = newshape
 
84
        newshape = list(v.shape)
 
85
        newshape[axis] = 2*N
 
86
        Y = zeros(newshape,sb.Complex)
 
87
        #Y[:N] = ak*v
 
88
        #Y[(N+1):] = conj(Y[N:0:-1])
 
89
        slices[0][axis] = slice(None,N)
 
90
        slices[1][axis] = slice(None,None)
 
91
        slices[2][axis] = slice(N+1,None)
 
92
        slices[3][axis] = slice((N-1),0,-1)
 
93
        Y[slices[0]] = ak*v
 
94
        Y[slices[2]] = conj(Y[slices[3]])
 
95
        x = sb.real(sb.ifft(Y,axis=axis))[slices[0]]
 
96
        return x
 
97
 
 
98
def dct2(x,axes=(-1,-2)):
 
99
    return dct(dct(x,axis=axes[0]),axis=axes[1])
 
100
 
 
101
def idct2(v,axes=(-1,-2)):
 
102
    return idct(idct(v,axis=axes[0]),axis=axes[1])
 
103
 
 
104
def dctn(x,axes=None):
 
105
    if axes is None:
 
106
        axes = sb.arange(len(x.shape))
 
107
    res = x
 
108
    for k in axes:
 
109
        res = dct(res,axis=k)
 
110
    return res
 
111
 
 
112
def idctn(v,axes=None):
 
113
    if axes is None:
 
114
        axes = sb.arange(len(v.shape))
 
115
    res = v
 
116
    for k in axes:
 
117
        res = idct(res,axis=k)
 
118
    return res
 
119
 
 
120
 
 
121
def makeC(N):
 
122
    n,l = ogrid[:N,:N]
 
123
    C = sb.cos(pi*(2*n+1)*l/(2*N))
 
124
    return C
 
125
 
 
126
def dct2raw(x):
 
127
    M,N = x.shape
 
128
    CM = makeC(M)
 
129
    CN = makeC(N)
 
130
    return dot(transpose(CM),dot(x,CN))
 
131
 
 
132
def idct2raw(v):
 
133
    M,N = v.shape
 
134
    iCM = linalg.inv(makeC(M))
 
135
    iCN = linalg.inv(makeC(N))
 
136
    return dot(transpose(iCM),dot(v,iCN))
 
137
 
 
138
def makeS(N):
 
139
    n,k = ogrid[:N,:N]
 
140
    C = sin(pi*(k+1)*(n+1)/(N+1))
 
141
    return C
 
142
 
 
143
# DST-I
 
144
def dst(x,axis=-1):
 
145
    """Discrete Sine Transform (DST-I)
 
146
 
 
147
    Implemented using 2(N+1)-point FFT
 
148
    xsym = r_[0,x,0,-x[::-1]]
 
149
    DST = (-imag(fft(xsym))/2)[1:(N+1)]
 
150
 
 
151
    adjusted to work over an arbitrary axis for entire n-dim array
 
152
    """
 
153
    n = len(x.shape)
 
154
    N = x.shape[axis]
 
155
    slices = [None]*3
 
156
    for k in range(3):
 
157
        slices[k] = []
 
158
        for j in range(n):
 
159
            slices[k].append(slice(None))
 
160
    newshape = list(x.shape)
 
161
    newshape[axis] = 2*(N+1)
 
162
    xtilde = sb.zeros(newshape,sb.float64)
 
163
    slices[0][axis] = slice(1,N+1)
 
164
    slices[1][axis] = slice(N+2,None)
 
165
    slices[2][axis] = slice(None,None,-1)
 
166
    for k in range(3):
 
167
        slices[k] = tuple(slices[k])
 
168
    xtilde[slices[0]] = x
 
169
    xtilde[slices[1]] = -x[slices[2]]
 
170
    Xt = sb.fft(xtilde,axis=axis)
 
171
    return (-sb.imag(Xt)/2)[slices[0]]
 
172
 
 
173
def idst(v,axis=-1):
 
174
    n = len(v.shape)
 
175
    N = v.shape[axis]
 
176
    slices = [None]*3
 
177
    for k in range(3):
 
178
        slices[k] = []
 
179
        for j in range(n):
 
180
            slices[k].append(slice(None))
 
181
    newshape = list(v.shape)
 
182
    newshape[axis] = 2*(N+1)
 
183
    Xt = sb.zeros(newshape,sb.Complex)
 
184
    slices[0][axis] = slice(1,N+1)
 
185
    slices[1][axis] = slice(N+2,None)
 
186
    slices[2][axis] = slice(None,None,-1)
 
187
    val = 2j*v
 
188
    for k in range(3):
 
189
        slices[k] = tuple(slices[k])
 
190
    Xt[slices[0]] = -val
 
191
    Xt[slices[1]] = val[slices[2]]
 
192
    xhat = sb.real(sb.ifft(Xt,axis=axis))
 
193
    return xhat[slices[0]]
 
194
 
 
195
def dst2(x,axes=(-1,-2)):
 
196
    return dst(dst(x,axis=axes[0]),axis=axes[1])
 
197
 
 
198
def idst2(v,axes=(-1,-2)):
 
199
    return idst(idst(v,axis=axes[0]),axis=axes[1])
 
200
 
 
201
def dstn(x,axes=None):
 
202
    if axes is None:
 
203
        axes = sb.arange(len(x.shape))
 
204
    res = x
 
205
    for k in axes:
 
206
        res = dst(res,axis=k)
 
207
    return res
 
208
 
 
209
def idstn(v,axes=None):
 
210
    if axes is None:
 
211
        axes = sb.arange(len(v.shape))
 
212
    res = v
 
213
    for k in axes:
 
214
        res = idst(res,axis=k)
 
215
    return res
 
216
 
 
217
def digitrevorder(x,base):
 
218
    x = sb.asarray(x)
 
219
    rem = N = len(x)
 
220
    L = 0
 
221
    while 1:
 
222
        if rem < base:
 
223
            break
 
224
        intd = rem // base
 
225
        if base*intd != rem:
 
226
            raise ValueError, "Length of data must be power of base."
 
227
        rem = intd
 
228
        L += 1
 
229
    vec = r_[[base**n for n in range(L)]]
 
230
    newx = x[sb.newaxis,:]*vec[:,sb.newaxis]
 
231
    # compute digits
 
232
    for k in range(L-1,-1,-1):
 
233
        newx[k] = x // vec[k]
 
234
        x = x - newx[k]*vec[k]
 
235
    # reverse digits
 
236
    newx = newx[::-1,:]
 
237
    x = 0*x
 
238
    # construct new result from reversed digts
 
239
    for k in range(L):
 
240
        x += newx[k]*vec[k]
 
241
    return x
 
242
 
 
243
 
 
244
def bitrevorder(x):
 
245
    return digitrevorder(x,2)
 
246
 
 
247
 
 
248
# needs to be fixed
 
249
def wht(data):
 
250
    """Walsh-Hadamaard Transform (sequency ordered)
 
251
 
 
252
    adapted from MATLAB algorithm published on the web by
 
253
    Author: Gylson Thomas
 
254
    e-mail: gylson_thomas@yahoo.com
 
255
    Asst. Professor, Electrical and Electronics Engineering Dept.
 
256
    MES College of Engineering Kuttippuram,
 
257
    Kerala, India, February 2005.
 
258
    copyright 2005.
 
259
    Reference: N.Ahmed, K.R. Rao, "Orthogonal Transformations for
 
260
    Digital Signal Processing" Spring Verlag, New York 1975. page-111.
 
261
    """
 
262
    N = len(data)
 
263
    L = sb.log2(N);
 
264
    if ((L-sb.floor(L)) > 0.0):
 
265
        raise ValueError, "Length must be power of 2"
 
266
    x=bitrevorder(data);
 
267
 
 
268
    k1=N; k2=1; k3=N/2;
 
269
    for i1 in range(1,L+1):  #Iteration stage
 
270
        L1=1;
 
271
        for i2 in range(1,k2+1):
 
272
            for i3 in range(1,k3+1):
 
273
                i=i3+L1-1; j=i+k3;
 
274
                temp1= x[i-1]; temp2 = x[j-1];
 
275
                if (i2 % 2) == 0:
 
276
                    x[i-1] = temp1 - temp2;
 
277
                    x[j-1] = temp1 + temp2;
 
278
                    x[i-1] = temp1 + temp2;
 
279
                    x[j-1] = temp1 - temp2;
 
280
                L1=L1+k1;
 
281
            k1 = k1/2;  k2 = k2*2;  k3 = k3/2;
 
282
    x = x*1.0/N; # Delete this line for inverse wht
 
283
    return x