Title: Hiding Images in Plain Sight: The Physics Of Magic Windows (2021/08/18)

Author: Matt Ferraro

Link: https://mattferraro.dev/posts/caustics-engineering

3b1b颁奖典礼.jpg（03:16开始）

12:40

今年夏天一些优秀的数学知识布道者

11.7万

171

以下为译文（前排提醒：公式较多）

Hiding Images in Plain Sight: The Physics Of Magic Windows

August 18, 2021 by Matt Ferraro

我最近做了违背直觉的物体，它是一个正方形的亚克力板，两面都很光滑，并且完全透明。

就像一个小窗户。

但它有一个神奇的特性：如果你用手电筒照它，它就会形成一个图像：

如果你把它拿到阳光下，它就会产生这种 3D 全息图：

由于b站不支持视频插入，原视频转成了gif，有画质损失

本文阐述了制作这个东西的数学原理，你也可以自己动手做！

首先，这是怎么做到的？

在谈论全息图之前，我们先关注一下二维图像。

这种物理现象被称作「焦散」（Caustics）。

焦散是指当光线穿过一个透明物体时看到的明亮的光斑。所有没能直接穿过物体的光子，投影形成了阴影。而这些光子总归要去往一个地方，它们就形成了这些焦散光斑图案。

焦散最有趣的方面是，即便是表面平整度最微小的变化，都能反映在最终的焦散图像中。即使是水池表面最温和的波浪，也会形成强大的透镜，在地板上投射出强烈的焦散现象。

我的亚克力方块之所以能产生图像，就是因为其表面的凹凸度分布恰到好处，使折射的光线能够生成特定的焦散图像。

为了能获得一点直观感受，我们首先考虑一个传统的凸透镜：

这种透镜形成了最简单的焦散。如果所有的入射光线都是来自一个无穷远处的点光源，例如太阳，这个透镜就会把所有的入射光线会聚到一个点。这个透镜的焦散图像就是周围漆黑，中心特别明亮的一个点。

放大透镜的一小部分，我们注意到一些特性：

1. 透镜的厚度对出射光线角度没有直接影响。我们可以加厚透镜左侧，这没有任何影响。第一个折射，从空气到玻璃（左边界），可以完全忽略

2. 入射光线与玻璃-空气边界（右边界）之间的角度对折射光线角度有很大影响

3. 两条光线是会聚还是发散是由透镜在玻璃与空气交界处的弯曲程度所控制的

换句话说，玻璃的厚度  并不重要，重要的是厚度的变化率 . 通过斯涅尔定律（折射定律），我们可以得到光线的出射角度。

光线会聚的地方，图像比光源更亮；光线发散的地方，图像比光源更暗。因此，图像的亮度（该点光线落下的地方）与  有关。

我的亚克力板的厚度在整个  平面上是连续变化的，因此我将其写做 ，将其看做一个高度图（heightmap）。

通过控制  和 ，我们可以将所有出射光线引导到图像的正确位置，同时产生正确的亮度，使其可识别。通过一些简化的假设，我们可以保证得到的高度图是平滑和连续的。

下图所示的「Magic Window」规格为 10cm × 10cm，表面凹凸的变化约为 2.0mm.

看到表面高度的轻微变化扭曲了地板上的直线吗？我们的「Magic Window」工作原理就像其他透镜一样，是通过弯曲光线实现的。

目录

（待补充）

抽象化问题

我们想找到一个高度图  对应焦散图像的亮度分布与某个输入图像相同。为了实现这一点，我们可以想象在亚克力透镜表面上有一个类似像素的单元格。透镜上的每个“像素”都对应着图像中的一个像素。图像像素和它们对应的透镜空间的“像素”被共同标记为共享的 (i,j) 坐标。

记住，  是标记像素列和行的整数，而  和  是以米或英寸为单位的实数。

解决步骤

第一步：对透镜上的单元格进行变形，使其变大或变小，从而使透镜单元  面积与图像单元  亮度成正比。由此产生的透镜网格不再是方块——必须引入大量的扭曲和倾斜以保持连续性。这一步是迄今为止最难的部分，必须迭代求解。

第二步：对每个单元  ，我们需要找到从透镜单元到图像单元的角度，并用折射定律找到所需的表面法线。这一步是纯粹的几何学。

第三步：整合所有表面法线，找到一个连续的高度图  。这也需要迭代求解，但是如果我们在第一步应用了某些限制条件，这一步会变得简单又快速。

对单元格进行变形

对于一个  像素的图像，透镜网格需要  个点。这是因为透镜网格的每个单元由四个点定义。从技术上讲，我们应该采用另一个坐标系来标记透镜网格的「点」，因为它们与透镜网格中的「单元」不同，但我认为只用  更简单。我们可以说，对于网格单元 ，左上方的点被定义为网格点 。

这一处理给我们留下了最底部的一行和最右边的一列多余网格点，不过这在实际处理的时候都不重要了。

透镜网格中的每个点 (i,j) 都有一个 (x,y) 坐标。一个点的 (i,j) 坐标永远不变，但 (x,y) 坐标会随着单元格变形而改变。

计算损失

已知所有透镜网格点的 (x,y) 坐标，通过简单的几何学可以计算出每个透镜网格单元的面积。当然一开始每个单元的面积都是一样的，但当我们开始变形后，情况就会改变。

我们想要的条件是：每个透镜网格单元 (i,j) 都有一个面积，它与图像像素的亮度 b(u,v) 成正比。

面积和亮度单位并不相同，所以用整个窗口的面积来规范化单元的面积，用图像的总亮度来规范化像素的亮度，这样就能用一个无单位的“百分比”来衡量每个单元。

直观地说，这意味着：

如果一个像素对整个图像亮度有 x% 的贡献，那么相应的窗口单元应该占整个窗口面积的 x%.

方程 (1.0) 是我们的目标，但是只有对窗口网格进行变形，它才能实现。要做到这一点，我们需要计算一个损失函数（loss function），告诉我们错过目标的程度。例如：

代码：

# Julia 风格伪代码
img = read_image("cat.png") #读取图像
brightness = convert_to_grayscale(img) #转灰度图
total_brightness = sum(brightness) 
brightness = brightness ./ total_brightness #单元亮度占比

w = .1 #米
h = .1 #米
area_total = w * h
loss = compute_pixel_area(grid) ./ area_total - brightness

上图右，我对损失函数进行了着色，红色表示网格单元需要扩大的区域，蓝色表示网格单元需要缩小的区域。

这张图就是损失函数 L ，我经常参考它。

迭代减小损失

可以认为损失图像是一个标量场 。标量场的梯度会产生一个矢量场，写做 。我们可以照着这个梯度场的方向移动每个网格点。太小的单元会扩大，太大的单元会缩小。最终减小损失值（loss），我们就能够创造出我们的Kitty🐱图像！

首先要做的是计算  并查看向量场：

我超。。。

 向量场表现太差了。。。有噪音、不连续，而且在很多地方等于零。

几乎在任何地方，相邻的点都需要在截然不同的方向上移动。这就导致改善一个单元的损失必然会恶化相邻单元的损失，这意味着在实践中，这种方法永远无法收敛。这是一个死胡同。

相反，让我们用计算流体力学（Computational Fluid Dynamics，CFD）来做个比喻。我们需要根据一个亮度函数来伸缩单元，这类似于可压缩气流建模，其中每个单元的压力都定义为一个压力函数。

如果二维网格中的每个单元都有一些初始压力，系统会怎么随时间而变化？高压区域扩大，低压区域收缩，中等压力的区域在某种整体的拉锯战中被推来推去。很明显，我们的问题是类似的。

那么，这个问题在 CFD 模拟中是怎样解决的呢？一个常用方法是定义一个称为  的速度势（Velocity Potential）。速度势是一个标量场，单位为 m^2/s . 乍一看可能很难理解，但是之所以说  很方便，是因为它的空间导数是以 m/s 来衡量的。换句话说，速度势  的梯度给出了一个矢量，其单位就是速度：

上图是一个  的例子。它只是一个标量场，因此最好看成一个高度图。

上图是同样  的梯度图。这些向量是指向上坡（uphill）的速度向量。在计算流体力学中，这些向量将表明流体如何从高压区域流向低压区域。

注意这个向量场表现得多好！在整个区域内变化温和平滑，相邻区域变化相似，并且所有的箭头都没有穿过边界。

在本文这个案例中，我们没有流体压力，但是有光的压力。图像中太亮的区域有很高的的光压，这在我们的损失函数 L 中被量化了。

如果我们能用 L 找到一个描述光压分布的  ，那么接下来需要做的就是计算出 . 然后用计算出的  对所有透镜网格点进行变形，从而减少损失。

那么怎么找到合适的  呢？我们已经知道的每个单元的属性就是它的损失（loss），它反映了单元需要增长或收缩的程度。

这个属性，即单元在随速度场运动时，随着时间推移膨胀或收缩的大小，被称作该场的散度（divergence）。

散度写作  ，所以在我们的例子中，我们要找到一个速度场  的散度等于损失，即满足：

不幸的是，并没有 “逆散度” 算子，因此我们不能直接倒转这个方程来求出  。不过我们能够将方程 (1.2) 代入 (1.3) 中，得到：

我们可以将其理解为：「速度势场  梯度的散度等于损失 loss」

这个等式在物理和数学的许多分支中经常出现，它通常被简写为：

你可能会认识到这就是 泊松方程！

这是一个好消息，因为泊松方程非常容易求解！如果你还不是很熟悉，你就可以把这一步想象成大矩阵求逆，或者对ODE（常微分方程）数值积分，或者实数开方。这是一项复杂而繁琐的工作，用纸笔来算会很痛苦，不过计算机正是用来处理这些的。

现在我们已经把这个问题写成了泊松方程，问题基本上就已经解决了。因为我们可以使用任何现成的求解器，插入已知的 L(x,y) ，使用诺依曼边界条件（Neumann boundary condition），然后砰的一声，就像变魔术一样求出了 .

（译者注：诺依曼边界条件就是第二类边界条件：给出未知函数在边界外法线的方向导数）

隐约可以看到一只小猫

你知道为什么猫🐱在这个  的 3D 渲染图中这么清晰吗？这个图中是什么控制了每个像素的亮度？（作者留的一个小思考题）

我们把求得的  插入方程（1.2）得到  ，看看这个向量场长什么样子：

有点失望。。。看起来并不像一只小猫咪🐱

从技术上讲，如果我们想减少损失 L ，我们就要向  的负方向上移动。这是  的向量场图：

不过好消息是，这个向量场十分平滑，并且表现良好。我们只需要沿着这个向量场移动网格点，就可以得到我们需要的东西。

如果你眯起眼睛从远处看，你几乎可以看到明亮的背景将扩大，而猫的黑毛将缩小。（想象成「光压」的概念）。

我们将所有的透镜网格点向着  的方向移动一小段。在对网格进行小幅度变形后，我们重新计算损失函数 L ，计算出新的  和新的  ，并再向前走一小步。

# Julia 风格伪代码
image = read_image("cat.png") #读取图像
gray = convert_to_grayscale(image) #转成灰度图
grid = create_initial_grid(gray.size + 1) #初始化矩阵

L = compute_loss(gray, grid) #计算损失函数

# 迭代计算直到 loss < 0.01
while max(L) > 0.01
    ϕ = poisson_solver(L, "neumann", 0)
    v = compute_gradient(ϕ)
    grid = step_grid(grid, -v)
    L = compute_loss(gray, grid)
end

经过三四次迭代后，损失变得非常小，我们就得到了变形后的网格！

看看这只猫🐱的下巴是如何膨胀起来的，但是鼻子和额头却缩小了。由于白色背景的膨胀，她的左耳明显被拉长了。她的瞳孔从椭圆形变成了尖尖的纺锤形。

请注意右边的图片只是Fusion360默认的网格渲染界面截图，它的线框是默认打开的：

Fusion360截图

某些区域颜色较深是因为这些区域的网格更密集。让我们放大来看看：

看看这多么的细节！我们甚至捕捉到了她眼睛里的明亮反光。进一步放大瞳孔：

我们可以看到网格单元的精细结构。（我们这次研究只考虑作为网格点的四边形，你在图中看到的三角形只是为了方便软件处理而从四边形简单分割而成的）

总结一下流程图：

（应该都能看懂，就没翻译了）

遵循这些步骤就能成功对网格点进行变形。现在，我们要做一些几何学的事情了！

折射定律和法线向量

折射定律告诉我们光线从一种介质到另一种介质是如何弯曲的。

其中 n1=1.49 是亚克力（acrylic，也做有机玻璃，丙烯酸塑料）的折射率，而 n2=1 是空气的折射率。如果已知 θ2，那么根据折射定律就能求出 θ1.

折射定律并不是什么物理学的公理，它是费马最短时间原理的直接结论，也是射线光学和波动光线的桥梁。但这是另一个话题了。

在我们的例子中，每个透镜单元 (i,j) 都已经移动到了 (x,y) 的位置，它需要把光线送到 (u,v) 的图像平面，两者相距某个距离 d 。

我们首先定义一个三维法向量 ，在每个点都垂直于我们的高度图。

法向量总是垂直于它们所处的平面。它们通常用方向而不是长度来定义，所以我们可以自由缩放大小。人们通常将法向量长度定为1，即单位法向量。

但是我们可以通过将 z 坐标定为 -1 来实现归一化：

这样就只用考虑 x 和 y，我们认识到：

这是一个在计算机图形学应用中经常使用的性质，同时也在涉及 DEM（数字高程模型）的地理空间应用中经常出现。（作者有关于DEM的进一步介绍文章）

利用折射定律、小角度近似和大量繁琐的几何学计算，我们可以找到法向量  的 x 分量和 y 分量：

这个推导没什么意思，我在这里就跳过了，感兴趣可以自己验证。

寻找高度图

目前为止，我们已经有了变形的网格单元，以及所有的表面法线。我们要做的就是找到一个高度图 h(x,y) ，使其满足所需的表面法线要求。

不幸的是，这并不是一个一般情况下能解决的问题。

我们可以尝试手动积分法线，从一个角开始，沿着网格往下走，但这种方法通常得不到一个物理上可实现的物体。

如果从左到右的法线积分将表面拉高，但是从上到下的法线积分却将表面压低，那么就没办法得到一个坚实的、连续的表面。

一个更好的方法是回到方程（2.2），重新写一遍：

并取两边的散度：

你还认识这个方程的形式吗？简写并交换两边：

我们得到了泊松方程的又一个例子！在上一节我们得到了 ，计算一个已知向量场的散度（divergence）是很容易的：

用代码实现是这样的：

δx = (Nx[i+1, j] - Nx[i, j])
δy = (Ny[i, j+1] - Ny[i, j])
divergence[i, j] = δx + δy

接下来就只是把已知的  代入到具有诺依曼边界条件的泊松求解器中，算得 h(x,y) 就可以使用了！

不过，还有一点需要改进。修改每个点的高度，我们实际上已经改变了每个透镜点到图像的距离，所以透镜—图像的距离不再是一个常数 d ，它实际上是一个函数 D(x,y) 。已知高度图了我们很容易计算出：

并不断重复这个过程，用 D(x,y) 而不是 d 来计算新的法线，然后又得到一个新的高度图。

不断循环这个过程并且确保收敛，实际上只需要两三次迭代就够了。

# Julia 风格伪代码
d = .2 #初始距离，单位：m
D = d .* array_of_ones(n, n) #距离函数 D(x,y)

# 迭代三次即可
for i in 1:3
    Nx, Ny = compute_normals(grid, D)
    divergence = compute_divergence(Nx, Ny)
    h = poisson_solver(divergence, "neumann", 0)
    D = copy(h)
end

所得到的高度图可以通过三角形网格分割，并闭合背面，来转换为一个实体对象。

注意到从正面看时，图像看起来像是镜像的。这是因为高度图形成的是后表面，前表面是一个平面。

高度的变化很微小，但足以完成工作：

制造

「Magic Window」的制造过程和任何其它 2.5D 物体的制造过程都相同。

我们将模型导入 Fusion360 或任何其他 CAM 软件（计算机辅助制造软件）。设置一个从左到右的粗加工刀具路径，以及一个从上到下的精加工刀具路径，就像你在大多数教程中看到的那样。

任何旧的数控雕刻机或者铣床（CNC router / mill）都可以。我去年设计并建造了自己的雕刻机。如果你也想做一台，建议从这里开始。（作者的另一篇博客：mattferraro.dev/posts/cnc-router）

我用的是一台四分之一英寸（6.35mm）直径的球鼻硬质合金钻头，同时用来粗加工和精加工，分别花了10分钟和90分钟。

雕刻后表面非常粗糙，并且是半透明的。我们需要用200、400、600、1000和1500目的砂纸进行湿磨，然后用软布和一些汽车抛光剂进行抛光。打磨和抛光对于一个 10cm×10cm 的「Magic Window」来说大约需要半个小时。

致谢

这篇文章所有的数学知识都来自于 Yue 等（2014）年发表的一篇现象级论文《Poisson-Based Continuous Surface Generation for Goal-Based Caustics》（原文见：http://nishitalab.org/user/egaku/tog14/yue-continuous-caustics-lens.pdf）

代码

源代码开源于此：https://github.com/MattFerraro/causticsEngineering

我是一个用 Julia 编程的新手，所以如果你有关于如何改进这段代码的建议，请联系我或直接Pull Request！

注意事项：我的代码有很多问题。我在几个地方混淆了 x 和 y。我插入了一些额外的负号，才使代码得以运行，但我不知道为什么。我的单位和符号在整个过程中是不一致的。原论文提出了一个更好的计算损失的方法，但我没有实现它，因为这个方法更简单。但我还是写出了我自己的网格工具和泊松求解器，因为我喜欢这种挑战 😆。

简而言之：对我来说，这个代码是一个有趣的 side project。如果你想通过这段代码来开展业务，你可能需要雇用一个知道如何用 Julia 进行专业编程的人。

许可证

MIT license 开源代码。无论是业余爱好、教育还是商业用途都欢迎自由使用。我只要求，如果你做出了什么东西，PLEASE SHOW ME!

该篇博文以 CC-BY 2.0 授权转载。

文章中的猫名叫 Mitski，她同意你使用她的照片作为图像处理论文的新标准参考图片。是时候让 Lenna 退休了。

One Last Thing

我知道你在想什么。那全息图呢！？

上面的数学推导意思是全息图总是会展现出来，还是这只猫的全息图只是一个不可思议的巧合？

好吧，你看，其实我已经找到了一个及其巧妙的证明，不幸的是，这个专栏字数限制，容纳不下这个证明 :)

联系作者

如果你用我的代码也做出来了一个「Magic Window」，或者在制造过程中有任何的困难，都欢迎联系我！

Twitter：@mferraro89

Email：mattferraro.dev@gmail.com

译者表示真希望自己能做一个出来，把模型发给 tb 代加工说不定可以？哈哈哈，希望有手工大佬实现出来，然后发个视频，分分钟上热门好吧🤣，我给你三连！（搜了一下还真没类似题材）

接下来尝试翻译一下3b1b其它获奖作品，欢迎关注~~~