← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

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