-
Notifications
You must be signed in to change notification settings - Fork 0
/
2024-07-15_caustique-Granet.py
275 lines (238 loc) · 12.6 KB
/
2024-07-15_caustique-Granet.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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
import os
import numpy as np
# the higher PRECISION, the bigger the file
# PRECISION = 4 # debugging
# PRECISION = 10 # good quality
PRECISION = 12 # for production
PRECISION = 13 # for production
# https://docs.python.org/3/library/dataclasses.html?highlight=dataclass#module-dataclasses
from dataclasses import dataclass
@dataclass
class init:
figpath: str = '2024-07-15_caustique-Granet' # Folder to store images
phi: float = 1.61803 # beauty is gold
tag: str = 'caustique' # Tag
ext: str = 'png' # Extension for output
nx: int = 5*2**PRECISION # number of pixels (vertical)
ny: int = 5*2**PRECISION # number of pixels (horizontal)
nframe: int = 1 # number of frames
bin_dens: int = 1 # relative bin density
bin_spectrum: int = 1 # bin spacing in spectrum (lower is more CPU)
seed: int = 2024 # seed for RNG
H: float = 20.0 # depth of the pool
variation: float = .20 # variation of diffraction index: http://www.philiplaven.com/p20.html 1.40 at 400 nm and 1.37 at 700nm makes a 2% variation
scale: float = .50*2**PRECISION # period in pixels
B_sf: float = 0.75 # bandwidth in sf
V_Y: float = 0.3 # horizontal speed
V_X: float = 0.3 # vertical speed
B_V: float = 1.0 # bandwidth in speed
zmin: float = 0.2 # gradient of wave height
# theta: float = 2*np.pi*(2-1.61803) # angle with the horizontal
theta: float = np.pi/2 # angle with the horizontal
B_theta: float = 60*np.pi/180 # bandwidth in theta
min_lum: float = .1 # diffusion level for the rendering
gamma: float = 2.9 # Gamma exponant to convert luminosity to luminance
fps: float = 18 # frames per second
multispectral: bool = True # Compute caustics on the full spectrogram.
cache: bool = True # Cache intermediate output.
verbose: bool = True # Displays more verbose output.
do_display: bool = False # Displays images in notebook.
do_recompute: bool = True # Restart each computation
from lambda2color import Lambda2color, xyz_from_xy
# borrowed from https://github.com/gummiks/gummiks.github.io/blob/master/scripts/astro/planck.py
def planck(wav, T):
import scipy.constants as const
c = const.c # c = 3.0e+8
h = const.h # h = 6.626e-34
k = const.k # k = 1.38e-23
a = 2.0*h*c**2
b = h*c/(wav*k*T)
intensity = a / ( (wav**5) * (np.exp(b) - 1.0) )
return intensity
def scattering(wav, a=0.005, p=1.3, b=0.45):
"""
b is proportionate to the column density of aerosols
along the path of sunlight, from outside the atmosphere
to the point of observation
see https://laurentperrinet.github.io/sciblog/posts/2020-07-04-colors-of-the-sky.html for more details
"""
# converting wav in µm:
intensity = np.exp(-a/((wav/1e-6)**4)) # Rayleigh extinction by nitrogen
intensity *= (wav/1e-6)**-4
intensity *= np.exp(-b/((wav/1e-6)**p)) # Aerosols
return intensity
import matplotlib
import matplotlib.pyplot as plt
from tqdm import tqdm, trange
import shutil
import hashlib
import MotionClouds as mc
class Caustique:
def __init__(self, opt):
"""
Image coordinates follow 'ij' indexing, that is,
* their origin at the top left,
* the X axis is vertical and goes "down",
* the Y axis is horizontal and goes "right".
"""
self.mc = mc
self.ratio = opt.ny/opt.nx # ratio between height and width (>1 for portrait, <1 for landscape)
X = np.linspace(0, 1, opt.nx, endpoint=False) # vertical
Y = np.linspace(0, self.ratio, opt.ny, endpoint=False) # horizontal
self.xv, self.yv = np.meshgrid(X, Y, indexing='ij')
self.opt = opt
# https://stackoverflow.com/questions/16878315/what-is-the-right-way-to-treat-python-argparse-namespace-as-a-dictionary
self.d = vars(opt)
os.makedirs(self.opt.figpath, exist_ok=True)
self.cachepath = os.path.join('/tmp', self.opt.figpath)
if opt.verbose: print(f'{self.cachepath=}')
os.makedirs(self.cachepath, exist_ok=True)
# a standard white:
# illuminant_D65 = xyz_from_xy(0.3127, 0.3291),
illuminant_sun = xyz_from_xy(0.325998, 0.335354)
# color conversion class
self.cs_srgb = Lambda2color(red=xyz_from_xy(0.64, 0.33),
green=xyz_from_xy(0.30, 0.60),
blue=xyz_from_xy(0.15, 0.06),
white=illuminant_sun)
self.wavelengths = self.cs_srgb.cmf[:, 0]*1e-9
self.N_wavelengths = len(self.wavelengths)
# multiply by the spectrum of the sky
intensity5800 = planck(self.wavelengths, 5800.)
scatter = scattering(self.wavelengths)
self.spectrum_sky = intensity5800 * scatter
self.spectrum_sky /= self.spectrum_sky.max()
def wave(self):
filename = f'{self.cachepath}/{self.opt.tag}_wave.npy'
if os.path.isfile(filename) and not(self.opt.do_recompute):
z = np.load(filename)
else:
# A simplistic model of a wave using https://github.com/NeuralEnsemble/MotionClouds
fx, fy, ft = mc.get_grids(self.opt.nx, self.opt.ny, self.opt.nframe)
env = mc.envelope_gabor(fx, fy, ft, V_X=self.opt.V_Y, V_Y=self.opt.V_X, B_V=self.opt.B_V,
sf_0=1./self.opt.scale, B_sf=self.opt.B_sf/self.opt.scale,
theta=self.opt.theta, B_theta=self.opt.B_theta)
z = mc.rectif(mc.random_cloud(env, seed=self.opt.seed))
if self.opt.cache: np.save(filename, z)
return z * np.linspace(self.opt.zmin, 1., self.opt.nx)[:, None, None]
def transform(self, z_, modulation=1.):
xv, yv = self.xv.copy(), self.yv.copy()
dzdx = z_ - np.roll(z_, 1, axis=0)
dzdy = z_ - np.roll(z_, 1, axis=1)
xv = xv + modulation * self.opt.H * dzdx
yv = yv + modulation * self.opt.H * dzdy
xv = np.mod(xv, 1)
yv = np.mod(yv, self.ratio)
return xv, yv
def plot(self, z, image=None, do_color=True, dpi=50):
"""
dpi: output resolution - sets the figure size as we ensure there is a one to one correspondance between pixels in the data and the output image
"""
# output filename
md5 = hashlib.sha224((self.opt.figpath + self.opt.tag).encode()).hexdigest()[:8] # an unique identifier for future tagging
output_filename = f'{self.opt.figpath}/{self.opt.tag}_{md5}.{self.opt.ext}'
if os.path.isfile(output_filename) and not(self.opt.do_recompute):
return output_filename
else:
# 1/ do the raytracing of image through z:
binsx, binsy = self.opt.nx//self.opt.bin_dens, self.opt.ny//self.opt.bin_dens
# a fixed image in degree of contrast (from 0=black to 1=white)
if image is None: image = np.ones((self.opt.nx, self.opt.ny))
#hist = self.do_raytracing(z)
# binsx, binsy = self.opt.nx//self.opt.bin_dens, self.opt.ny//self.opt.bin_dens
subplotpars = matplotlib.figure.SubplotParams(left=0., right=1., bottom=0., top=1., wspace=0., hspace=0.,)
if self.opt.multispectral:
#image_rgb = self.cs_srgb.spec_to_rgb(hist)
image_rgb = np.zeros((self.opt.nx//self.opt.bin_dens, self.opt.ny//self.opt.bin_dens, 3, self.opt.nframe))
for i_frame in trange(self.opt.nframe):
for i_wavelength in trange(self.opt.bin_spectrum//2, self.N_wavelengths, self.opt.bin_spectrum):
modulation = 1. + self.opt.variation/2 - self.opt.variation*i_wavelength/self.N_wavelengths
xv, yv = self.transform(z[:, :, i_frame], modulation=modulation)
hist_, _, _ = np.histogram2d(xv.ravel(), yv.ravel(),
bins=[binsx, binsy],
weights=image.ravel(),
range=[[0, 1], [0, self.ratio]],
density=True)
# we convert the spectrum into a color
spec = np.zeros((self.N_wavelengths))
spec[i_wavelength] = 1
rgb = self.cs_srgb.spec_to_rgb(spec)
rgb *= self.spectrum_sky[i_wavelength]
# we add the color to the image
image_rgb[:, :, :, i_frame] += hist_[:, :, None] * rgb[None, None, :]
image_rgb -= image_rgb.min()
image_rgb /= image_rgb.max()
else:
hist = np.zeros((binsx, binsy, self.opt.nframe))
for i_frame in trange(self.opt.nframe):
xv, yv = self.transform(z[:, :, i_frame])
hist_, _, _ = np.histogram2d(xv.ravel(), yv.ravel(),
bins=[binsx, binsy],
range=[[0, 1], [0, self.ratio]],
density=True)
#hist /= hist.max()
# 2/ transform light into image:
fnames = []
for i_frame in trange(self.opt.nframe):
fig, ax = plt.subplots(figsize=(self.opt.ny/self.opt.bin_dens/dpi, self.opt.nx/self.opt.bin_dens/dpi), subplotpars=subplotpars)
if self.opt.multispectral:
ax.imshow(image_rgb[:, :, :, i_frame] ** (1/self.opt.gamma), vmin=0, vmax=1)
else:
if do_color:
bluesky = np.array([0.268375, 0.283377]) # xyz
sun = np.array([0.325998, 0.335354]) # xyz
# ax.pcolormesh(edge_y, edge_x, hist[:, :, i_frame], vmin=0, vmax=1, cmap=plt.cm.Blues_r)
# https://en.wikipedia.org/wiki/CIE_1931_color_space#Mixing_colors_specified_with_the_CIE_xy_chromaticity_diagram
L1 = 1 - hist[:, :, i_frame]
L2 = hist[:, :, i_frame]
image_denom = L1 / bluesky[1] + L2 / sun[1]
image_x = (L1 * bluesky[0] / bluesky[1] + L2 * sun[0] / sun[1]) / image_denom
image_y = (L1 + L2) / image_denom
image_xyz = np.dstack((image_x, image_y, 1 - image_x - image_y))
image_rgb = self.cs_srgb.xyz_to_rgb(image_xyz)
image_L = self.opt.min_lum + (1-self.opt.min_lum)* L2 ** .61803
ax.imshow(image_L[:, :, None]*image_rgb, vmin=0, vmax=1)
else:
ax.imshow(1-image_L, vmin=0, vmax=1)
fname = f'{self.cachepath}/{self.opt.tag}_frame_{i_frame:04d}.png'
fig.savefig(fname, dpi=dpi)
fnames.append(fname)
plt.close('all')
if self.opt.nframe==1:
shutil.copyfile(fname, output_filename)
return output_filename
else:
if self.opt.ext == 'gif':
return make_gif(output_filename, fnames, fps=self.opt.fps)
else:
return make_mp4(output_filename, fnames, fps=self.opt.fps)
def show(self, output_filename, width=1024):
from IPython.display import HTML, Image, display
if self.opt.nframe==1:
display(Image(url=output_filename.replace(self.opt.ext, 'png'), width=width))
else:
if self.opt.ext == 'gif':
return display(Image(url=output_filename, width=width))
else:
#import moviepy.editor as mpy
#return mpy.ipython_display(output_filename, width=width)
# https://github.com/NeuralEnsemble/MotionClouds/blob/master/MotionClouds/MotionClouds.py#L858
opts = ' loop="1" autoplay="1" controls '
html = HTML(f'<video {opts} width="{width}"> <source src="{output_filename}" type="video/{self.opt.ext}" /> </video>')
html.reload()
return display(html)
def generate_image(nx, ny, periods=6.25, threshold=0.45, radius=.9):
X, Y = np.meshgrid(np.linspace(-1, 1, nx, endpoint=True), np.linspace(-1, 1, ny, endpoint=True))
image = (np.cos(2*np.pi*X*periods) > threshold)*1.
# image += (np.cos(2*np.pi*Y*periods) > threshold)*1.
# image = (image>=1) * 1.
image *= (X**2 < radius**2) * (Y**2 < radius**2) * 1.
return image
opt = init()
c = Caustique(opt)
c.opt.tag = f'{c.opt.figpath}-{PRECISION=}'
z = c.wave()
image = generate_image(nx=c.opt.nx, ny=c.opt.ny)
output_filename = c.plot(z, image)
print(output_filename)
print('Done !')