Last week, I read Michael Berry’s paper, “Laplacian Magic Windows.” Over the years, I have read many interesting papers by this longtime Mathematica user, but this one stood out for its maximizing of the product of simplicity and unexpectedness. Michael discusses what he calls the magic window. For 70+ years, we have known about holograms, and now we know about magic windows. So what exactly is a magic window? Here is a sketch of the optics of one:

Parallel light falls onto a glass sheet that is planar on the one side and has some gentle surface variation on the other side (bumps in the above image are vastly overemphasized; the bumps of a real magic window would be minuscule). The light gets refracted by the magic window (the deviation angles of the refracted light rays are also overemphasized in the graphic) and falls onto a wall. Although the window bumpiness shows no recognizable shape or pattern, the light density variations on the wall show a clearly recognizable image. Starting with the image that one wants to see on the wall, one can always construct a window that shows the image one has selected. The variations in the thickness of the glass are assumed to be quite small, and the imaging plane is assumed to be not too far away so that the refracted light does not form caustics—as one sees them, for instance, at the bottom of a swimming pool in sunny conditions.

Now, how should the window surface look to generate any pre-selected image on the wall? It turns out that the image visible on the wall is the Laplacian of the window surface. Magic windows sound like magic, but they are just calculus (differentiation, to be precise) in action. Isn’t this a neat application of multivariate calculus? Schematically, these are the mathematical steps involved in a magic window.

Implementation-wise, the core steps are the following:

And while magic windows are a 2017 invention, their roots go back hundreds of years to so-called magic mirrors. Magic mirrors are the mirror equivalent of magic windows: they too can act as optical Laplace operators (see the following).

Expressed more mathematically: Let the height of the bumpy side of the glass surface be f(x,y). Then the intensity of the light brightness on the wall is approximately Δ_x,y f(x,y), where Δ is the Laplacian ∂²./∂x²+∂²./∂y². Michael calls such a window a “magic window.” It is magic because the glass surface height f(x,y) does not in any way resemble Δ _x,y f(x,y).

It sounds miraculous that a window can operate as a Laplace operator. So let’s do some numerical experiments to convince ourselves that this does really work. Let’s start with a goat that we want to use as the image to be modeled. We just import a cute-looking Oberhasli dairy goat from the internet.

The gray values of the pixels can be viewed as a function h: ℝ*ℝ->[0,1]. Interpolation allows us to use this function constructively.

Here is a 3D plot of the goat function ifGoat.

And we can solve the Poisson equation with the image as the right-hand side: Δ f = image using NDSolveValue.

We will use Dirichlet boundary conditions for now. (But the boundary conditions will not matter for the main argument.)

The Poisson equation solution is quite a smooth function; the inverse of the Laplace operator is a smoothing operation. No visual trace of the goat seems to be left.

Overall it is smooth, and it is also still smooth when zoomed in.

Even when repeatedly zoomed in.

The overall shape of the Poisson equation solution can be easily understood through the Green’s function of the Laplace operator.

We calculate and visualize the first few terms (individually) of the double sum from the Green’s function.

Taking 25² terms into account, we have the following approximations for the Poisson equation solution and its Laplacian. The overall shape is the same as the previous numerical solution of the Poisson equation.

For this small number of Fourier modes, the outline of goat is recognizable, but its details aren’t.

Applying the Laplace operator to the PDE solutions recovers (by construction) a version of the goat. Due to finite element discretization and numerical differentiation, the resulting goat is not quite the original one.

A faster and less discretization-dependent way to solve the Poisson equation uses the fast Fourier transform (FFT).

This solution recovers the goat more faithfully. Here is the recovered goat after interpolating the function values.

Taking into account that any physical realization of a magic window made from glass will unavoidably have imperfections, a natural question to ask is: What happens if one adds small perturbations to the solution of the Poisson equation?

The next input modifies each grid point randomly by a perturbation of relative size 10^–p. We see that for this goat, the relative precision of the surface has to be on the order of 10^-6 or better—a challenging but realizable mechanical accuracy.

To see how the goat emerges after differentiation (Δ = ∂² ./∂x²+∂² ./∂y²), here are the partial derivatives shown.

And because we have ∂² ./∂x²+∂² ./∂y² =(∂./∂x+ⅈ ∂./∂y)(∂./∂x–ⅈ ∂./∂y), we also look at the Wirtinger derivatives.

We could also just use a simple finite difference formula to get the goat back. This avoids any interpolation artifacts and absolutely faithfully reproduces the original goat.

The differentiation can even be carried out as an image processing operation.

So far, nothing really interesting. We integrated and differentiated a function. Let’s switch gears and consider the refraction of a set of parallel light rays on a glass sheet.

We consider the lower parts of the glass sheet planar and the upper part slightly wavy with an explicit description height = f(x,y). The index of refraction is n, and we follow the light rays (coming from below) up to the imaging plane at height = Z. Here is a small Manipulate that visualizes this situation for the example surface f(x,y) = 1 + ε (cos(α x) + cos(β y).

We do want the upper surface of the glass nearly planar, so we use the factor ε in the previous equation.

The reason we want the upper surface mostly planar is that we want to avoid rays that “cross” near the surface and form caustics. We want to be in a situation where the density of the rays is position dependent, but the rays do not yet cross. This restricts the values of n, Z and the height of the surface modulation.

Now let’s do the refraction experiment with the previous solution of the Laplace equation as the height of the upper glass surface. To make the surface variations small, we multiply that solution by 0.0001.

We use the median refractive index of glass, n = 1.53.

Instead of using lightRay, we will use a compiled version for faster numerical evaluation.

In absolute units, say the variations in glass height are at most 1 mm; we look at the refracted rays a few meters behind the glass window. We will use about 3.2 million lights rays (4² per pixel).

Displaying all endpoints of the rays gives a rather strong Moiré effect. But the goat is visible—a true refraction goat!

If we accumulate the number of points that arrive in a small neighborhood of the given points {X,Y} in the plane height=Z, the goat becomes much more visible. (This is what would happen if we would observe the brightness of the light that goes through the glass sheet and we assume that the intensities are additive.) To do the accumulation efficiently, we use the function Nearest.

Note that looking into the light that comes through the window would not show the goat because the light that would fall into the eye would mostly come from a small spatial region due to the mostly parallel light rays.

The appearance of the Laplacian of the surface of the glass sheet is not restricted to only parallel light. In the following, we use a point light source instead of parallel light. This means that the effect would also be visible by using artificial light sources, rather than sunlight with a magic window.

So why is the goat visible in the density of rays after refraction? At first, it seems quite surprising whether either a parallel or point source shines on the window.

On second thought, one remembers Maxwell’s geometric meaning of the Laplace operator:

… where indicates the average of f on a sphere centered at with radius ρ. Here is a quick check of the last identity for two and three dimensions.

At a given point in the imaging plane, we add up the light rays from different points of the glass surface. This means we carry out some kind of averaging operation.

So let’s go back to the general refraction formula and have a closer look. Again we assume that the upper surface is mostly flat and that the parameter ε is small. The position {X,Y} of the light ray in the imaging plane can be calculated in closed form as a function of the surface g(x,y), the starting coordinates of the light ray {x,y}, the index of refraction n and the distance of the imaging plane Z.

That is a relatively complicated-looking formula. For a nearly planar upper glass surface (small ε), we have the following approximate coordinates for the {X,Y} coordinates of the imaging plane where we observe the light rays in terms of the coordinate {x,y} of the glass surface.

This means in zeroth order we have {X,Y} ≈ {x,y}. And the deviation of the light ray position in the imaging plane is proportional (n–1)Z. (Higher-order corrections to {X,Y} ≈ {x,y} we could get from Newton iterations, but we do not need them here.)

The density of rays is the inverse of the Jacobian for going from {x,y} to {X,Y}. (Think on the change of variable formulas for 1:1 transforms for multivariate integration.)

Quantifying the size of the resulting expression shows that it is indeed a large expression. This is quite a complex formula. For a quadratic function in x and y, we can get some feeling for the density as a function of the physical parameters ε, n and Z as well as the parameters that describe the surface by varying them in an interactive demonstration. For large values of n, Z and ε, we see how caustics arise.

For nearly planar surfaces (first order in ε), the density is equal to the Laplacian of the surface heights (in x,y coordinates). This is the main “trick” in the construction of magic windows.

This explains why the goat appears as the intensity pattern of the light rays after refraction. This means glass sheets act effectively as a Laplace operator.

Using Newton’s root-finding method, we could calculate the intensity in X,Y coordinates, but the expression explains heuristically why refraction on a nearly planar surface behaves like an optical Laplace operator. For more details, see this article.

Now we could model a better picture of the light ray density by pre-generating a matrix of points in the imaging plane using, say, 10 million rays, and record where they fall within the imaging plane. This time we model the solution of the Poisson equation using ListDeconvolve.

The approximate solution of the Poisson equation is not quite as smooth as the global solutions, but the goat is nevertheless invisible.

We adjust the brightness/darkness through a power law (a crude approximation for a Weber–Fechner perception).

If the imaging plane is too far away, we do get caustics (that remind me of the famous cave paintings from Lascaux).

If the image plane is even further away, the goat slowly becomes unrecognizable.

Although not practically realizable, we also show what the goat would look like for negative Z; now it seems much more sheep-like.

Here is a small animation showing the shape of the goat as a function of the distance Z of the imaging plane from the upper surface.

Even if the image is just made from a few lines (rather than each pixel having a non-white or non-black value), the solution of the Poisson equation is a smooth function, and the right-hand side is not recognizable in a plot of the solution.

But after refraction on a glass sheet (or applying the Laplacian), we see Homer quite clearly.

Despite the very localized curve-like structures that make the Homer image, the resulting Poisson equation solution again looks quite smooth. Here is the solution textured with its second derivative (the purple line will be used in the next input).

The next graphic shows a cross-section of the Poisson equation solution together with its (scaled) first and second derivatives with respect to x along the purple line of the last graphic. The lines show up quite pronounced in the second derivatives.

Let’s repeat a modification of the previous experiment to see how precise the surface would have to be to show Homer. We add some random waves to the Homer solution.

Again we see that the surface would have to be correct at the (10^-6) level or better.

Or one can design a nearly planar window that will project one’s most favorite physics equation on the wall when the Sun is shining.

When looking at the window, one will not notice any formulas. But this time, the solution of the Poisson equation has more overall structures.

But the refracted light will make physics equations. The resulting window is perfect for the entrance of, say, physics department buildings.

Now that we’re at the end of this post, let us mention that one can also implement the Laplacian through a mirror, rather than a window. See Michael Berry’s paper from 2006, “Oriental Magic Mirrors and the Laplacian Image” (see this article as well). Modifying the above function for refracting a light ray to reflecting a light ray and assuming a mostly flat mirror surface, we see the Laplacian of the mirror surface in the reflected light intensity.

Making transparent materials and mirrors of arbitrary shape, now called free-form optics, is considered the next generation of modern optics and will have wide applications in science, technology, architecture and art (see here). I think that a few years from now, when the advertising industry recognizes their potential, we will see magic windows with their unexpected images behind them everywhere.

---

Download this post as a Computable Document Format (CDF) file. New to CDF? Get your copy for free with this one-time download.