{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Caustics ([wikipedia](https://en.wikipedia.org/wiki/Caustic_%28optics%29) are luminous patterns which are resulting from the superposition of smoothly deviated light rays. It is for instance the heart-shaped pattern in your cup of coffee which is formed as the rays of the sun are reflected on the cup's surface. It is also the wiggly pattern of light curves that you will see on the floor of a pool as the sun's light is *refracted* at the surface of the water. Here, we simulate that particular physical phenomenon. Simply because they are mesmerizingly beautiful, but also as it is of interest in visual neuroscience. Indeed, it speaks to how images are formed (more on this later), hence how the brain may understand images. \n", "\n", "In this post, I will develop a simple formalism to generate such patterns, with the paradoxical result that it is *very* simple to code yet generates patterns with great complexity, such as: \n", "
\n", "
\n", "
\n", "\n", "This is joint work with artist [Etienne Rey](https://laurentperrinet.github.io/authors/etienne-rey/), in which I especially follow the ideas put forward in the series [Turbulence](http://ondesparalleles.org/projets/turbulences/).\n", "\n", "\"DOI\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "Let's begin with an import of the libraries that we will need: [numpy](https://numpy.org/) and [matplotlib](https://matplotlib.org/) (essentials!) but also [MotionClouds](https://github.com/NeuralEnsemble/MotionClouds) to generate a simple waveform (more on this later in this post). These are all installable using [pip](https://pip.pypa.io/en/stable/) (see the [requirements.txt](https://github.com/NaturalPatterns/2020_caustiques/blob/master/requirements.txt) file)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", "import matplotlib\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import MotionClouds as mc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's define an unique identifier" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Tagging our simulations with tag=2020-06-19_caustique\n" ] } ], "source": [ "import datetime\n", "date = datetime.datetime.now().date().isoformat() # https://en.wikipedia.org/wiki/ISO_8601\n", "tag = f'{date}_caustique'\n", "tag = '2020-06-19_caustique' # overwriting tag\n", "print(f'Tagging our simulations with tag={tag}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And a `Namespace` object to store all parameters:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def init(args=[], ds=1):\n", " import argparse\n", "\n", " parser = argparse.ArgumentParser()\n", " parser.add_argument(\"--tag\", type=str, default='caustique', help=\"Tag\")\n", " parser.add_argument(\"--figpath\", type=str, default='../files/2020-06-19_caustique', help=\"Folder to store images\")\n", " parser.add_argument(\"--nx\", type=int, default=5*2**8, help=\"number of pixels (vertical)\")\n", " parser.add_argument(\"--ny\", type=int, default=8*2**8, help=\"number of pixels (horizontal)\")\n", " parser.add_argument(\"--bin_dens\", type=int, default=4, help=\"relative bin density\")\n", " parser.add_argument(\"--nframe\", type=int, default=2**7, help=\"number of frames\")\n", " parser.add_argument(\"--seed\", type=int, default=42, help=\"seed for RNG\")\n", " parser.add_argument(\"--H\", type=float, default=10., help=\"depth of the pool\")\n", " parser.add_argument(\"--sf_0\", type=float, default=0.004, help=\"sf\")\n", " parser.add_argument(\"--B_sf\", type=float, default=0.002, help=\"bandwidth in sf\")\n", " parser.add_argument(\"--V_Y\", type=float, default=0.3, help=\"horizontal speed\")\n", " parser.add_argument(\"--V_X\", type=float, default=0.3, help=\"vertical speed\")\n", " parser.add_argument(\"--B_V\", type=float, default=4.0, help=\"bandwidth in speed\")\n", " parser.add_argument(\"--theta\", type=float, default=2*np.pi*(2-1.61803), help=\"angle with the horizontal\")\n", " parser.add_argument(\"--B_theta\", type=float, default=np.pi/8, help=\"bandwidth in theta\")\n", " parser.add_argument(\"--min_lum\", type=float, default=.2, help=\"diffusion level for the rendering\")\n", " parser.add_argument(\"--fps\", type=float, default=18, help=\"bandwidth in theta\")\n", " parser.add_argument(\"--cache\", type=bool, default=False, help=\"Cache intermediate output.\")\n", " parser.add_argument(\"--verbose\", type=bool, default=False, help=\"Displays more verbose output.\")\n", "\n", " opt = parser.parse_args(args=args)\n", "\n", " if opt.verbose:\n", " print(opt)\n", " return opt\n", "\n", "opt = init()\n", "opt.tag = tag" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Namespace(B_V=4.0, B_sf=0.002, B_theta=0.39269908169872414, H=10.0, V_X=0.3, V_Y=0.3, bin_dens=4, cache=False, figpath='../files/2020-06-19_caustique', fps=18, min_lum=0.2, nframe=128, nx=1280, ny=2048, seed=42, sf_0=0.004, tag='2020-06-19_caustique', theta=2.399988291783386, verbose=False)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This object is always handy when we will try to explore some parameters.\n", "\n", "A last utility function is a way to convert the resulting numpy array to a set of images and then to an animated GIF using [imageio](https://imageio.github.io/) and [pygifsicle](https://github.com/LucaCappelletti94/pygifsicle):" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def make_gif(gifname, fnames, fps):\n", " import imageio\n", "\n", " with imageio.get_writer(gifname, mode='I', fps=fps) as writer:\n", " for fname in fnames:\n", " writer.append_data(imageio.imread(fname))\n", "\n", " from pygifsicle import optimize\n", " optimize(str(gifname))\n", " return gifname" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from IPython.display import Image, display\n", "width = 1024" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# What is a caustic?\n", "\n", "A caustic usually refers to the image which is formed by the convergence or divergence of rays on a surface. Mathematically, this may be understood as the *enveloppe* of a set of rays, and intuitively, it is best understood by looking at the best known example: When the sun hits a glass or cup and you see heart-shaped pattern (actually more a kidney than a heart, hence the name \"nephroid\"):\n", "\n", "
\n", "
\n", "
\n", "\n", "(image by [User:Paul venter](https://commons.wikimedia.org/wiki/User:Paul_venter)).\n", "\n", "Their study is essential in optics, as these are the bases for lenses and therefore of central importance in our image-centered society :-)\n", "\n", "They are also of interest simply because they create mesmerizing patterns in nature, for instance the pattern of wiggly lines on swimming pool bottoms, or on the sea floor at shallow depths. \n", "\n", "The formation of that specific pattern is very well illustrated by this animation (credit [Jacopo Bertolotti](https://twitter.com/j_bertolotti), check-out that twitter [thread](https://twitter.com/j_bertolotti/status/1229743952594640897)) where light comes from the bottom, travels through air, hits the surface and gets deviated:\n", "\n", "
\n", "
\n", "
\n", "\n", "What happens? Light travels fast (*very* fast) but it will travel slower in water than in the air. A consequence of this is that a ray changes slightly its direction as it hits the water: the light ray is *refracted*. That is why a pen which is half immerged looks \"as if\" it is broken. In fact, it's only its image that gets distorted:\n", " \n", "
\n", "
\"Pen
\n", "
\n", "\n", "The laws governing this deviations are well known as [Snell's law](https://en.wikipedia.org/wiki/Snell%27s_law) and state that\n", "\n", "$$\n", "\\frac{\\sin\\theta_2}{\\sin\\theta_1} = \\frac{n_1}{n_2}\n", "$$\n", "Where $n_1\\approx 1$ and $n_2\\approx1.33$ are the refractive indices of air and water, respectively. (Note that, as a consequence, the light in the water is always slightly closer to the normal of the surface.)\n", "\n", "
\n", "
\n", "\"Snells \n", "
\n", "
\n", " \n", "Going back to the swimming pool, considering that the sun is emitting an uniform pattern of parallel rays from an incident angle $\\theta_I$, we may apply this law to all of these light rays. The important ingredient is that the surface of the water (given as the surface height $z(x, y)$) changes smoothly and will deviate the light rays non-uniformly, thus creating the patterns. Writing down the equation, the displacement resulting from the angle of refraction is given by:\n", "\n", "$$\n", "\\Delta x(x, y) = H \\cdot \\arcsin [ \\frac 1 n \\sin( \\theta_I - \\arctan(\\frac{\\partial}{\\partial x} z(x, y)) + \\arctan(\\frac{\\partial}{\\partial x} z(x, y) ]\n", "$$\n", "\n", "Where $H$ is the pool's depth, $n$ is the refractive index of water (that of air being set to $1$) and $\\theta_I$ the incident angle.\n", "The same formula applies in the other spatial dimension $y$ by computing the angle of the surface along this dimension as $\\arctan(\\frac{\\partial}{\\partial y} z(x, y))$. \n", "\n", "In the limit of a zenithal light source (that is, the sun at the vertical of the pool) and in the limit of small angles, the equations simplifies to:\n", "$$\n", "\\Delta x(x, y) \\approx H' \\cdot \\frac{\\partial z(x, y)}{\\partial x}\n", "$$\n", "and\n", "$$\n", "\\Delta y(x, y) \\approx H' \\cdot \\frac{\\partial z(x, y)}{\\partial y}\n", "$$\n", "\n", "with $ H' = \\frac{n-1}{n} H$. What is interesting here is that we will transform an incoming \"image\" (here the uniform light source from the sun) by a set of small and local deformations of the support of the image. The equations are ultra-simple, but the consequences will paradoxically look complex.\n", "Thinking about images as densities which we may discretize on a rectangular grid (in the `__init__` function), it is very easy to generate a wave pattern ass a smooth surface in space and time (in the `wave` function), simulate the transformation on that surface (in the `transform` function) and compute the resulting image as an histogram (in the `plot` function):" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "class Caustique:\n", " def __init__(self, opt):\n", " \"\"\"\n", " Image coordinates follow 'ij' indexing, that is,\n", " * their origin at the top left,\n", " * the X axis is vertical and goes \"down\",\n", " * the Y axis is horizontal and goes \"right\".\n", "\n", " \"\"\"\n", " self.ratio = opt.ny/opt.nx # ratio between height and width (>1 for portrait, <1 for landscape)\n", " X = np.linspace(0, 1, opt.nx, endpoint=False) # vertical\n", " Y = np.linspace(0, self.ratio, opt.ny, endpoint=False) # horizontal\n", " self.xv, self.yv = np.meshgrid(X, Y, indexing='ij')\n", " self.opt = opt\n", " # https://stackoverflow.com/questions/16878315/what-is-the-right-way-to-treat-python-argparse-namespace-as-a-dictionary\n", " self.d = vars(opt)\n", " os.makedirs(self.opt.figpath, exist_ok=True)\n", " if self.opt.cache: \n", " self.cachepath = os.path.join('/tmp', self.opt.figpath)\n", " os.makedirs(self.cachepath, exist_ok=True)\n", "\n", " # a standard white:\n", " illuminant_D65 = xyz_from_xy(0.3127, 0.3291), \n", " illuminant_sun = xyz_from_xy(0.325998, 0.335354)\n", " # color conversion class\n", " self.cs_srgb = ColourSystem(red=xyz_from_xy(0.64, 0.33),\n", " green=xyz_from_xy(0.30, 0.60),\n", " blue=xyz_from_xy(0.15, 0.06),\n", " white=illuminant_sun) \n", " def wave(self):\n", " filename = f'{self.cachepath}/{self.opt.tag}_wave.npy'\n", " if self.opt.cache and os.path.isfile(filename):\n", " z = np.load(filename)\n", " else:\n", " # A simplistic model of a wave using https://github.com/NeuralEnsemble/MotionClouds\n", " import MotionClouds as mc\n", " fx, fy, ft = mc.get_grids(self.opt.nx, self.opt.ny, self.opt.nframe)\n", " 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,\n", " sf_0=self.opt.sf_0, B_sf=self.opt.B_sf,\n", " theta=self.opt.theta, B_theta=self.opt.B_theta)\n", " z = mc.rectif(mc.random_cloud(env, seed=self.opt.seed))\n", " if self.opt.cache: np.save(filename, z)\n", " return z\n", "\n", " def transform(self, z_):\n", " xv, yv = self.xv.copy(), self.yv.copy()\n", "\n", " dzdx = z_ - np.roll(z_, 1, axis=0)\n", " dzdy = z_ - np.roll(z_, 1, axis=1)\n", " xv = xv + self.opt.H * dzdx\n", " yv = yv + self.opt.H * dzdy\n", "\n", " xv = np.mod(xv, 1)\n", " yv = np.mod(yv, self.ratio)\n", "\n", " return xv, yv\n", "\n", " def plot(self, z, gifname=None, dpi=150):\n", " if gifname is None:\n", " gifname=f'{self.opt.figpath}/{self.opt.tag}.gif'\n", "\n", " filename = f'{self.cachepath}/{self.opt.tag}_hist.npy'\n", " binsx, binsy = self.opt.nx//self.opt.bin_dens, self.opt.ny//self.opt.bin_dens\n", " if self.opt.cache and os.path.isfile(filename):\n", " hist = np.load(filename)\n", " else:\n", " hist = np.zeros((binsx, binsy, self.opt.nframe))\n", " for i_frame in range(self.opt.nframe):\n", " xv, yv = self.transform(z[:, :, i_frame])\n", " hist[:, :, i_frame], edge_x, edge_y = np.histogram2d(xv.ravel(), yv.ravel(),\n", " bins=[binsx, binsy],\n", " range=[[0, 1], [0, self.ratio]],\n", " density=True)\n", " hist /= hist.max()\n", " if self.opt.cache: np.save(filename, hist)\n", " \n", " subplotpars = matplotlib.figure.SubplotParams(left=0., right=1., bottom=0., top=1., wspace=0., hspace=0.,)\n", " fnames = []\n", " for i_frame in range(self.opt.nframe):\n", " fig, ax = plt.subplots(figsize=(binsy/dpi, binsx/dpi), subplotpars=subplotpars)\n", " bluesky = np.array([0.268375, 0.283377]) # xyz\n", " sun = np.array([0.325998, 0.335354]) # xyz\n", " # ax.pcolormesh(edge_y, edge_x, hist[:, :, i_frame], vmin=0, vmax=1, cmap=plt.cm.Blues_r)\n", " # https://en.wikipedia.org/wiki/CIE_1931_color_space#Mixing_colors_specified_with_the_CIE_xy_chromaticity_diagram\n", " L1 = 1 - hist[:, :, i_frame]\n", " L2 = hist[:, :, i_frame]\n", " image_denom = L1 / bluesky[1] + L2 / sun[1] \n", " image_x = (L1 * bluesky[0] / bluesky[1] + L2 * sun[0] / sun[1]) / image_denom\n", " image_y = (L1 + L2) / image_denom \n", " image_xyz = np.dstack((image_x, image_y, 1 - image_x - image_y))\n", " image_rgb = self.cs_srgb.xyz_to_rgb(image_xyz)\n", " image_L = self.opt.min_lum + (1-self.opt.min_lum)* L2 ** .61803\n", " \n", " ax.imshow(image_L[:, :, None]*image_rgb, vmin=0, vmax=1)\n", "\n", " fname = f'{self.cachepath}/{self.opt.tag}_frame_{i_frame}.png'\n", " fig.savefig(fname, dpi=dpi)\n", " fnames.append(fname)\n", " plt.close('all')\n", "\n", " return make_gif(gifname, fnames, fps=self.opt.fps)\n", "\n", " \n", "# the rest is borrowed from https://github.com/gummiks/gummiks.github.io/blob/master/scripts/astro/planck.py\n", "# as synthesized in https://laurentperrinet.github.io/sciblog/posts/2020-07-04-colors-of-the-sky.html\n", "\n", "def planck(wav, T):\n", " import scipy.constants as const\n", " c = const.c # c = 3.0e+8\n", " h = const.h # h = 6.626e-34\n", " k = const.k # k = 1.38e-23\n", " a = 2.0*h*c**2\n", " b = h*c/(wav*k*T)\n", " intensity = a / ( (wav**5) * (np.exp(b) - 1.0) )\n", " return intensity\n", "\n", "def scattering(wav, a=0.005, p=1.3, b=0.45):\n", " \"\"\"\n", " b is proportionate to the column density of aerosols \n", " along the path of sunlight, from outside the atmosphere \n", " to the point of observation\n", " \n", " N_O3 is the ozone column density along the path of sunlight, \n", " sigma_O3 is the wavelength dependent ozone absorption cross-section.\n", " \n", " \"\"\"\n", " # converting wav in µm:\n", " intensity = np.exp(-a/((wav/1e-6)**4)) # Rayleigh extinction by nitrogen\n", " intensity *= (wav/1e-6)**-4\n", " intensity *= np.exp(-b/((wav/1e-6)**p)) # Aerosols\n", " return intensity\n", "\n", "def xyz_from_xy(x, y):\n", " \"\"\"Return the vector (x, y, 1-x-y).\"\"\"\n", " return np.array((x, y, 1-x-y))\n", "\n", "class ColourSystem:\n", " \"\"\"A class representing a colour system.\n", "\n", " A colour system defined by the CIE x, y and z=1-x-y coordinates of\n", " its three primary illuminants and its \"white point\".\n", "\n", " TODO: Implement gamma correction\n", "\n", " \"\"\"\n", "\n", " # CMF is the CIE colour matching function for 380 - 780 nm in 5 nm intervals\n", "\n", " def __init__(self, red, green, blue, white):\n", " \"\"\"Initialise the ColourSystem object.\n", "\n", " Pass vectors (ie NumPy arrays of shape (3,)) for each of the\n", " red, green, blue chromaticities and the white illuminant\n", " defining the colour system.\n", "\n", " \"\"\"\n", "\n", " # Chromaticities\n", " self.red, self.green, self.blue = red, green, blue\n", " self.white = white\n", " # The chromaticity matrix (rgb -> xyz) and its inverse\n", " self.M = np.vstack((self.red, self.green, self.blue)).T \n", " self.MI = np.linalg.inv(self.M)\n", " # White scaling array\n", " self.wscale = self.MI.dot(self.white)\n", " # xyz -> rgb transformation matrix\n", " self.T = self.MI / self.wscale[:, np.newaxis]\n", "\n", " CMF_str = \"\"\"380 0.0014 0.0000 0.0065\n", " 385 0.0022 0.0001 0.0105\n", " 390 0.0042 0.0001 0.0201\n", " 395 0.0076 0.0002 0.0362\n", " 400 0.0143 0.0004 0.0679\n", " 405 0.0232 0.0006 0.1102\n", " 410 0.0435 0.0012 0.2074\n", " 415 0.0776 0.0022 0.3713\n", " 420 0.1344 0.0040 0.6456\n", " 425 0.2148 0.0073 1.0391\n", " 430 0.2839 0.0116 1.3856\n", " 435 0.3285 0.0168 1.6230\n", " 440 0.3483 0.0230 1.7471\n", " 445 0.3481 0.0298 1.7826\n", " 450 0.3362 0.0380 1.7721\n", " 455 0.3187 0.0480 1.7441\n", " 460 0.2908 0.0600 1.6692\n", " 465 0.2511 0.0739 1.5281\n", " 470 0.1954 0.0910 1.2876\n", " 475 0.1421 0.1126 1.0419\n", " 480 0.0956 0.1390 0.8130\n", " 485 0.0580 0.1693 0.6162\n", " 490 0.0320 0.2080 0.4652\n", " 495 0.0147 0.2586 0.3533\n", " 500 0.0049 0.3230 0.2720\n", " 505 0.0024 0.4073 0.2123\n", " 510 0.0093 0.5030 0.1582\n", " 515 0.0291 0.6082 0.1117\n", " 520 0.0633 0.7100 0.0782\n", " 525 0.1096 0.7932 0.0573\n", " 530 0.1655 0.8620 0.0422\n", " 535 0.2257 0.9149 0.0298\n", " 540 0.2904 0.9540 0.0203\n", " 545 0.3597 0.9803 0.0134\n", " 550 0.4334 0.9950 0.0087\n", " 555 0.5121 1.0000 0.0057\n", " 560 0.5945 0.9950 0.0039\n", " 565 0.6784 0.9786 0.0027\n", " 570 0.7621 0.9520 0.0021\n", " 575 0.8425 0.9154 0.0018\n", " 580 0.9163 0.8700 0.0017\n", " 585 0.9786 0.8163 0.0014\n", " 590 1.0263 0.7570 0.0011\n", " 595 1.0567 0.6949 0.0010\n", " 600 1.0622 0.6310 0.0008\n", " 605 1.0456 0.5668 0.0006\n", " 610 1.0026 0.5030 0.0003\n", " 615 0.9384 0.4412 0.0002\n", " 620 0.8544 0.3810 0.0002\n", " 625 0.7514 0.3210 0.0001\n", " 630 0.6424 0.2650 0.0000\n", " 635 0.5419 0.2170 0.0000\n", " 640 0.4479 0.1750 0.0000\n", " 645 0.3608 0.1382 0.0000\n", " 650 0.2835 0.1070 0.0000\n", " 655 0.2187 0.0816 0.0000\n", " 660 0.1649 0.0610 0.0000\n", " 665 0.1212 0.0446 0.0000\n", " 670 0.0874 0.0320 0.0000\n", " 675 0.0636 0.0232 0.0000\n", " 680 0.0468 0.0170 0.0000\n", " 685 0.0329 0.0119 0.0000\n", " 690 0.0227 0.0082 0.0000\n", " 695 0.0158 0.0057 0.0000\n", " 700 0.0114 0.0041 0.0000\n", " 705 0.0081 0.0029 0.0000\n", " 710 0.0058 0.0021 0.0000\n", " 715 0.0041 0.0015 0.0000\n", " 720 0.0029 0.0010 0.0000\n", " 725 0.0020 0.0007 0.0000\n", " 730 0.0014 0.0005 0.0000\n", " 735 0.0010 0.0004 0.0000\n", " 740 0.0007 0.0002 0.0000\n", " 745 0.0005 0.0002 0.0000\n", " 750 0.0003 0.0001 0.0000\n", " 755 0.0002 0.0001 0.0000\n", " 760 0.0002 0.0001 0.0000\n", " 765 0.0001 0.0000 0.0000\n", " 770 0.0001 0.0000 0.0000\n", " 775 0.0001 0.0000 0.0000\n", " 780 0.0000 0.0000 0.0000\"\"\"\n", " CMF = np.zeros((len(CMF_str.split('\\n')), 4))\n", " for i, line in enumerate(CMF_str.split('\\n')): \n", " CMF[i, :] = np.fromstring(line, sep=' ')\n", "\n", " self.cmf = CMF\n", " \n", " self.wavelengths = np.linspace(200e-9, 2000e-9, 200) \n", " self.intensity5800 = planck(self.wavelengths, 5800.)\n", " self.intensity5800 /= self.intensity5800.max()\n", " self.intensitysky = self.intensity5800 * scattering(self.wavelengths)\n", "\n", "\n", " def xyz_to_rgb(self, xyz):\n", " \"\"\"Transform from xyz to rgb representation of colour.\n", "\n", " The output rgb components are normalized on their maximum\n", " value. If xyz is out the rgb gamut, it is desaturated until it\n", " comes into gamut.\n", "\n", " By default, fractional rgb components are returned.\n", "\n", " \"\"\"\n", "\n", " # rgb = self.T.dot(xyz)\n", " rgb = np.tensordot(xyz, self.T.T, axes=1)\n", " if np.any(rgb < 0):\n", " # We're not in the RGB gamut: approximate by desaturating\n", " w = - np.min(rgb)\n", " rgb += w\n", " if not np.all(rgb==0):\n", " # Normalize the rgb vector\n", " rgb /= np.max(rgb)\n", "\n", " return rgb\n", "\n", " def spec_to_xyz(self, spec):\n", " \"\"\"Convert a spectrum to an xyz point.\n", "\n", " The spectrum must be on the same grid of points as the colour-matching\n", " function, self.cmf: 380-780 nm in 5 nm steps.\n", "\n", " \"\"\"\n", "\n", " XYZ = np.sum(spec[:, np.newaxis] * self.cmf, axis=0)\n", " den = np.sum(XYZ)\n", " if den == 0.:\n", " return XYZ\n", " return XYZ / den\n", "\n", " def spec_to_rgb(self, spec, out_fmt=None):\n", " \"\"\"Convert a spectrum to an rgb value.\"\"\"\n", "\n", " xyz = self.spec_to_xyz(spec)\n", " return self.xyz_to_rgb(xyz, out_fmt)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# a first instance\n", "\n", "Let's generate a first instance of this caustic with the default parameters. Let's first have a look at the (simplistic) waveform we are using:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "
\n", "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "if not os.path.isfile(f'{opt.figpath}/{opt.tag}{mc.vext}'):\n", " print('Doing', f'{opt.figpath}/{opt.tag}.{mc.vext}')\n", " c = Caustique(opt)\n", " z = c.wave()\n", " mc.anim_save(z.swapaxes(0, 1), f'{opt.figpath}/{opt.tag}')\n", "mc.in_show_video(f'{opt.tag}', figpath=opt.figpath)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It has many desirable properties: it is smooth in time and space, it is stochastic and periodic (no border effect). As such we may apply our simple transform and as a result we get:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "gifname = f'{opt.figpath}/{opt.tag}.gif'\n", "if not os.path.isfile(gifname):\n", " c = Caustique(opt)\n", " z = c.wave()\n", " gifname = c.plot(z)\n", "display(Image(url=gifname, width=width))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This image ressembles the wiggling lines that we expected to find, great! Some properties:\n", "\n", "* it is periodic in space and time by construction which may look unrealistic for a pool, but looks better as you do not have border effects...\n", "* there are different classes of visual objects: veil-like smooth curves, traveling waves with a difference in their group/phase speed, nephroid-like intersection points, smooth gradients, ... which mainly depend on the surface's ondulation,\n", "* when looking at this video long enough (I recommend opening this animated GIF in a separate tab and to zoom on it fullscreen), you may perceive some three-dimensionnal percepts, especially when some features change their speed..\n", "\n", "## exploring parameters\n", "Is there more to see in our system?\n", "The first thing to do is to explore is the effect of the pool's depth:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "N_scan = 9\n", "base = 2" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H = 5.000\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "H = 5.946\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "H = 7.071\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "H = 8.409\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "H = 10.000\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "H = 11.892\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "H = 14.142\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "H = 16.818\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "H = 20.000\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "opt = init()\n", "opt.tag = tag\n", "\n", "c = Caustique(opt)\n", "z = None\n", "for H_ in c.opt.H*np.logspace(-1, 1, N_scan, base=base):\n", " print(f'H = {H_:.3f}')\n", " c.opt.H = H_\n", " gifname = f'{opt.figpath}/{opt.tag}_H_{H_:.3f}.gif'\n", " if not os.path.isfile(gifname):\n", " if z is None:\n", " z = c.wave()\n", " url=c.plot(z, gifname=gifname)\n", " display(Image(url=gifname, width=width))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This shows that for this waveform, there is an optimal depth to generate more concentrated caustics. This corresponds to the depth where the curvature of the waveform concentrates light as a lens. Note that at great distortions some ripples appear which are due to our sampling of space.\n", "\n", "We may now explore the parameters of this waveform to see their effect on the resulting caustic:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "======sf_0======\n", "sf_0=sf_0(default)*0.500=2.000E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "sf_0=sf_0(default)*0.595=2.378E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "sf_0=sf_0(default)*0.707=2.828E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "sf_0=sf_0(default)*0.841=3.364E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "sf_0=sf_0(default)*1.000=4.000E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "sf_0=sf_0(default)*1.189=4.757E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "sf_0=sf_0(default)*1.414=5.657E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "sf_0=sf_0(default)*1.682=6.727E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "sf_0=sf_0(default)*2.000=8.000E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "======B_sf======\n", "B_sf=B_sf(default)*0.500=1.000E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_sf=B_sf(default)*0.595=1.189E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_sf=B_sf(default)*0.707=1.414E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_sf=B_sf(default)*0.841=1.682E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_sf=B_sf(default)*1.000=2.000E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_sf=B_sf(default)*1.189=2.378E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_sf=B_sf(default)*1.414=2.828E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_sf=B_sf(default)*1.682=3.364E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_sf=B_sf(default)*2.000=4.000E-03\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "======theta======\n", "theta=theta(default)*0.500=1.200E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "theta=theta(default)*0.595=1.427E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "theta=theta(default)*0.707=1.697E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "theta=theta(default)*0.841=2.018E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "theta=theta(default)*1.000=2.400E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "theta=theta(default)*1.189=2.854E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "theta=theta(default)*1.414=3.394E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "theta=theta(default)*1.682=4.036E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "theta=theta(default)*2.000=4.800E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "======B_theta======\n", "B_theta=B_theta(default)*0.500=1.963E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_theta=B_theta(default)*0.595=2.335E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_theta=B_theta(default)*0.707=2.777E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_theta=B_theta(default)*0.841=3.302E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_theta=B_theta(default)*1.000=3.927E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_theta=B_theta(default)*1.189=4.670E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_theta=B_theta(default)*1.414=5.554E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_theta=B_theta(default)*1.682=6.604E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_theta=B_theta(default)*2.000=7.854E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "======V_X======\n", "V_X=V_X(default)*0.500=1.500E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_X=V_X(default)*0.595=1.784E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_X=V_X(default)*0.707=2.121E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_X=V_X(default)*0.841=2.523E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_X=V_X(default)*1.000=3.000E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_X=V_X(default)*1.189=3.568E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_X=V_X(default)*1.414=4.243E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_X=V_X(default)*1.682=5.045E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_X=V_X(default)*2.000=6.000E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "======V_Y======\n", "V_Y=V_Y(default)*0.500=1.500E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_Y=V_Y(default)*0.595=1.784E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_Y=V_Y(default)*0.707=2.121E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_Y=V_Y(default)*0.841=2.523E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_Y=V_Y(default)*1.000=3.000E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_Y=V_Y(default)*1.189=3.568E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_Y=V_Y(default)*1.414=4.243E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_Y=V_Y(default)*1.682=5.045E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "V_Y=V_Y(default)*2.000=6.000E-01\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "======B_V======\n", "B_V=B_V(default)*0.500=2.000E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_V=B_V(default)*0.595=2.378E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_V=B_V(default)*0.707=2.828E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_V=B_V(default)*0.841=3.364E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_V=B_V(default)*1.000=4.000E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_V=B_V(default)*1.189=4.757E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_V=B_V(default)*1.414=5.657E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_V=B_V(default)*1.682=6.727E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "B_V=B_V(default)*2.000=8.000E+00\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "opt = init()\n", "\n", "for variable in ['sf_0', 'B_sf', 'theta', 'B_theta', 'V_X', 'V_Y', 'B_V']:\n", " print(f'======{variable}======')\n", " for modul in np.logspace(-1, 1, N_scan, base=base):\n", " opt = init()\n", " opt.tag = tag\n", "\n", " c = Caustique(opt)\n", " c.d[variable] *= modul\n", "\n", " print(f'{variable}={variable}(default)*{modul:.3f}={c.d[variable]:.3E}')\n", " gifname = f'{opt.figpath}/{opt.tag}_{variable}_modul_{modul:.3f}.gif'\n", " if not os.path.isfile(gifname):\n", " print('Doing ', gifname)\n", " z = c.wave()\n", " mcname = f'{opt.figpath}/{opt.tag}_{variable}_modul_{modul:.3f}'\n", " if not os.path.isfile(f'{mcname}{mc.vext}'): \n", " print('Doing ', f'{mcname}{mc.vext}')\n", " mc.anim_save(z.swapaxes(0, 1), f'{mcname}')\n", " url=c.plot(z, gifname=gifname)\n", " display(Image(url=gifname, width=width))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Future developments should:\n", "\n", "* check if the approximation of small angles qualitatively changes the resulting images (should not, what maters is the transport of the density, not it's precision),\n", "* use a more realistic waveform such as in https://laurentperrinet.github.io/sciblog/posts/2016-04-24_a-wave-going-backwards.html\n", "* use different indices for different wavelengths (+add gamma correction for rendering),\n", "* make it live (using pyglet for instance)\n", "\n", "This development happens in https://github.com/NaturalPatterns/2020_caustiques - I am happy to welcome any people that would be interested in exploring this.\n", "\n", "This shows that for this waveform, there is an optimal depth to generate more concentrated caustics. This corresponds to the depth for which the curvature of the waveform concentrates light as a lens. \n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 4 }