1. 摘要

自从 Horn 在 20 世纪 70 年代初开发出第一种从阴影中恢复形状 (SfS) 技术以来，出现了许多不同的方法。在本文中，复述了1994年蔡博士的论文 Shape from shading using a linear approximation，对表面法线（分量p和分量q）使用有限差分法进行离散近似，然后线性化深度反射函数 Z(x,y)。并将该算法在示例图像上进行了测试，成功重建三维(Tsai et al., 1994)，并对试验结果进行简单分析。

1. 简介

Shape-from-Shading （从阴影中恢复形状）是一种经典的三维重建技术。它利用单目图像（或多目图像）中的图像强度信息，结合光照方向与传感器（例如照相机）方向对目标表面三维信息进行估计，不需要立体像对。SfS技术主要的估计算法分为四种：最小化方法、传播方法、局部方法和线性方法(Zhang et al., 1998)。本篇内容参考著名的Tsai & Shah线性化方法，在朗伯表面实现三维重建。包括公式推导与python代码实现。

1. 前置知识

3.1 朗伯模型

朗伯模型是一种最经典也是最简单的反射模型。朗伯模型假设从粗糙表面反射的光的强度仅取决于表面法线和光源方向，如图Fig. 1。因此，朗伯表面在所有方向上都同样明亮。

Fig. 1. 朗伯模型示意图

反射是通过对表面的单位法向量 和一个归一化的光照方向矢量 进行点积来计算的, 如Eq.(1)所示。

(1)

3.2 复合函数求导

(2)

3.3 一元函数泰勒展开

(3)

3.4 雅克比迭代法

雅克比迭代法是一种求解线性方程组的方法，可以在不断迭代后准确的求解线性方程组的仅此解。在此就不过多赘述，但可以简单地理解为不断将估算值输入进线性方程组中，得到的解再次输入线性方程组中，每次得到的解将更接近线性方程组的精确解。

1. 假设

我们SfS算法本质上是一种对表面深度的估算，所以我们很难在不进行任何理想化考虑的情况下进行准确的三维重建。因此，为了减少未知数同时也为了简化计算，我们将进行如下假设：

1. 朗伯模型

1. 平行光源

1. 暂不考虑地形相互遮挡

1. 重构的形状产生与输入图像相同的亮度

也就是说，在朗伯表面的每一个微元，都由相同方位的光源进行照射，并且产生的阴影不会由更远处的地形产生遮挡。以上四条假设也是本篇所介绍的SfS方法的核心之一。

1. 使用线性近似实现Shape-from-Shading

5.1 参数设置

: 位于图像 处的像素强度值

: 朗伯反射强度值

: 表面法线在X方向的偏导数

: 表面法线在Y方向的偏导数

: 光照方位角

: 光照入射角

: 光照向量在X方向的分量

: 光照向量在Y方向的分量

5.2 Shape-from-Shading

5.2.1 计算表面梯度

我们使用最简单的梯度计算方法，也就是某像素在X方向的上的梯度由它和它左边的像素进行计算，在Y方向上的梯度由它和它上边的像素进行计算，如Eq.(4)和Eq.(5)。因此，用这种方法计算得到的梯度图在斜向45°上拥有较好地效果，也导致整体SfS在斜向45°方向上有较好地重建效果。

(4)

(5)

5.2.2 计算反射强度

根据假设4，我们可以直接令初始图像强度 = 反射强度， 也就是Eq.(6)

(6)

将计算表面梯度的公式Eq.(4)和Eq.(5)代入朗伯反射模型公式Eq.(1)中，得到公式Eq.(7)

(7)

令:

(8)

将Eq.(7)代入Eq.(8)，并进行泰勒一阶展开得到Eq.(9)

(9)

其中，是在第次迭代中恢复的三维表面。令并进行移项，将移动到等号左边，得到Eq.(10)

(10)

其中，

(11)

这里涉及到复合函数求导与分式求导，具体知识可参考同济大学《高等数学》。经过以上计算，我们得到了最终的参数关系：

最后，将初始三维设置为0，也就是，这样就是可以在一个水平的表面逐渐恢复三维信息了。

1. Python实现

"""References: "Ping-Sing, T., & Shah, M. (1994). Shape from shading using linear approximation.""Ruo, Z., Ping-Sing, T., Cryer, J. E., & Shah, M. (1999). Shape-from-shading: a survey"
"""
from PIL import Image
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D# Path of the image.
img_path = "mozart.jpg"# Read image as an array.
E = Image.open(img_path) # Define the light source orientation.
tilt  = ...
slant = ...tilt  = np.deg2rad(tilt)
slant = np.deg2rad(slant)# Define the surface normal.
p = np.zeros_like(E)
q = np.zeros_like(E)# Define the baseline surface.
Z = np.zeros_like(E)# Define the surface derivatives in x and y directions.
Z_x = np.zeros_like(E)
Z_y = np.zeros_like(E)# Define the maximum number of iterations.
maxIter = 200# Calculate the normalized illumination direction
ix = np.cos(tilt) * np.tan(slant)
iy = np.sin(tilt) * np.tan(slant)# Get the image size
(h, w) = np.shape(E)# Jacobi Iteration
for i in range(maxIter):# Using the illumination direction and the currently estimate# surface normals, compute the corresponding reflectance map.# refer to Eq.(7) .R = (np.cos(slant) + p * np.cos(tilt)*np.sin(slant)+ q * np.sin(tilt) * np.sin(slant)) / np.sqrt(1 + p**2 + q**2)# At each iteration, make sure that the reflectance map is positive at# each pixel, set negative values to 0.R = np.maximum(0, R)# Compute our function f which is the deviation of the computed# reflectance map from the original image. Refer to Eq.(8)f = E - R# Compute the derivative of f with respect to our surface Z, refer to Eq.(11)df_dZ =(p+q) * (ix*p + iy*q + 1) / (np.sqrt((1 + p**2 + q**2)**3)*np.sqrt(1 + ix**2 + iy**2))-(ix+iy)/(np.sqrt(1 + p**2 + q**2)*np.sqrt(1 + ix**2 + iy**2))# Update the surface, refer to (10), and please avoid dividing by zeroZ = Z - f/(df_dZ + np.finfo(np.float32).eps)# Compute the surface derivatives with respect to x and y.Z_x[1:h,:] = Z[0:h-1,:]Z_y[:,1:w] = Z[:,0:w-1]# Using the updated surface, compute new surface normals, refer to Eq.(4) and Eq.(5)p = Z - Z_xq = Z - Z_y# Show the result
fig  = plt.figure()
ax3  = plt.axes(projection = '3d')
xx   = np.arange(0, h, 1)
yy   = np.arange(0, w, 1)
X, Y = np.meshgrid(xx, yy)
ax3.plot_surface(X, Y, Z, rstride = 4, cstride = 4, cmap = 'Greys_r')
plt.show()

1. 结果与拓展讨论

我们使用的图像是经典的计算机视觉示例图像Mozart，如图Fig.2。

Fig. 2. 所使用的输入图像"mozart.jpg"

经过200次迭代，所建立的三维模型如图Fig.3所示。

Fig. 3. 在不同角度下观察三维重建结果

可以看到，人物的鼻子是凹陷的，明显有悖于我们的常识。其原因是因为原始图像中人物的鼻子两侧有阴影，类似的阴影也可以被认为是鼻子对光线造成的遮挡，违背了我们的假设4,，并且鼻子部分的像素值与脸颊部分的像素值相近，而我们的算法是完全依照像素值进行计算，无法分辨实际真实情况，这也使此算法和朗伯模型的弊端。为此我们可以在设置时，也就是代码line33，将设置为一个已有的DEM数据。这种方法已经应用在相当多的领域，当然，在我们本次的实验中，此方法也很难纠正上述问题。

参考文献

Ping-Sing, T., & Shah, M. (1994). Shape from shading using linear approximation. Image and vision computing, 12(8), 487-498. https://doi.org/10.1016/0262-8856(94)90002-7

Ruo, Z., Ping-Sing, T., Cryer, J. E., & Shah, M. (1999). Shape-from-shading: a survey. IEEE transactions on pattern analysis and machine intelligence, 21(8), 690-706. https://doi.org/10.1109/34.784284