from tqdm.notebook import tqdm, trange
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)
print(self.opt.figpath)
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):
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
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, do_color=True, videoname=None, dpi=50):
#hist = self.do_raytracing(z)
binsx, binsy = self.opt.nx//self.opt.bin_dens, self.opt.ny//self.opt.bin_dens
if videoname is None:
videoname=f'{self.opt.figpath}/{self.opt.tag}.{self.opt.vext}'
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 tqdm(range(self.opt.nframe)):
for ii_wavelength, i_wavelength in enumerate(range(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
# print(i_wavelength, N_wavelengths, modulation)
xv, yv = self.transform(z[:, :, i_frame], modulation=modulation)
hist_, edge_x, edge_y = np.histogram2d(xv.ravel(), yv.ravel(),
bins=[binsx, binsy],
range=[[0, 1], [0, self.ratio]],
density=True)
spec = np.zeros((self.N_wavelengths))
spec[i_wavelength] = 1
rgb = self.cs_srgb.spec_to_rgb(spec)
rgb *= self.spectrum_sky[i_wavelength]
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_, edge_x, edge_y = np.histogram2d(xv.ravel(), yv.ravel(),
bins=[binsx, binsy],
range=[[0, 1], [0, self.ratio]],
density=True)
#hist /= hist.max()
fnames = []
for i_frame in range(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.vext == 'gif':
return make_gif(videoname, fnames, fps=self.opt.fps)
else:
return make_mp4(videoname, fnames, fps=self.opt.fps)
def show(self, videoname, width=1024):
if self.opt.vext == 'gif':
from IPython.display import Image, display
return display(Image(url=videoname, width=width))
else:
#import moviepy.editor as mpy
#return mpy.ipython_display(videoname, width=width)
# https://github.com/NeuralEnsemble/MotionClouds/blob/master/MotionClouds/MotionClouds.py#L858
from IPython.display import HTML, display
opts = ' loop="1" autoplay="1" controls '
html = HTML(f'<video {opts} src="{videoname}" type="video/{self.opt.vext}" width={width}\>')
html.reload()
return display(html)
