-
Notifications
You must be signed in to change notification settings - Fork 0
/
psClouds.py
119 lines (94 loc) · 2.32 KB
/
psClouds.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
import numpy as np
import testImage as t
from pylab import *
def zpad(image):
"""Stupid old version of numpy doesn't have pad module in it """
xsize,ysize=image.shape
padded = np.zeros( (xsize*2,ysize*2))
padded[xsize/2:xsize*1.5,ysize/2:ysize*1.5] = image
return padded
def myplot(ack):
figure()
subplot(2,2,1)
imshow(ack.image)
title('Image')
subplot(2,2,2)
imshow(ack.fimage.real, vmin=-1,vmax=1)
title('real FFT')
colorbar()
subplot(2,2,3)
imshow(log(ack.psd2d), vmin=-10,vmax=10,origin='lower')
title('log 2D Power Spectrum')
colorbar()
subplot(2,2,4)
good = np.where(ack.psd1d != 0)
semilogy(ack.rcenters[good],ack.psd1d[good])
xlabel('Spatial Frequency')
ylabel('1-D Power Spectrum')
ylim([ack.psd1d[ack.rcenters.size-20], np.max(ack.psd1d)])
oldCLouds = np.load('oldclouds.npy')
newClouds = np.load('newclouds.npy')
ack = t.TestImage()
cloudim = zpad(oldCLouds.copy())
ack.setImage(cloudim)
ack.calcFft()
ack.calcPsd2d()
ack.calcPsd1d()
figure()
subplot(2,2,1)
imshow(ack.image)
title('Image')
subplot(2,2,2)
imshow(ack.fimage.real, vmin=-1,vmax=1)
title('real FFT')
colorbar()
subplot(2,2,3)
imshow(log(ack.psd2d), vmin=-10,vmax=10,origin='lower')
title('log 2D Power Spectrum')
colorbar()
subplot(2,2,4)
semilogy(ack.psd1d)
xlabel('Spatial Frequency')
ylabel('1-D Power Spectrum')
savefig('oldclouds.png',type='png')
clf()
cloudim = zpad(newClouds.copy())
ack.setImage(cloudim)
ack.calcFft()
ack.calcPsd2d()
ack.calcPsd1d()
figure()
subplot(2,2,1)
imshow(ack.image)
title('Image')
subplot(2,2,2)
imshow(ack.fimage.real, vmin=-1,vmax=1)
title('real FFT')
colorbar()
subplot(2,2,3)
imshow(log(ack.psd2d), vmin=-10,vmax=10,origin='lower')
title('log 2D Power Spectrum')
colorbar()
subplot(2,2,4)
semilogy(ack.psd1d)
xlabel('Spatial Frequency')
ylabel('1-D Power Spectrum')
savefig('newclouds.png',type='png')
xx = ack.xx.copy()
yy = ack.yy.copy()
im = t.TestImage()
scales = np.array([5.,10.,500.])
for i in np.arange(np.size(scales)):
scale = scales[i]
imsin = np.sin(2.*np.pi/scale*xx)#*im.xx**2.
imss = imsin*np.sin(2.*np.pi/scale*yy)
im.setImage(zpad(imss))
im.calcFft()
im.calcPsd2d()
im.calcPsd1d()
myplot(im)
ylim([1e-5,1e7])
savefig('sin_ps%d.png'%scale, type='png')
clf()
#im.showFft()
#im.showPsd1d()