从阴影恢复形状技术是计算机视觉中活跃的一个分支^[1]。这种算法通常采用建立图像亮度与反射图表面之间的几何函数的方法来实现3维重建的过程。传统方法大多都是将朗伯体模型作为反射模型，再通过变分法找到病态方程组中系数的最接近的值，从而得到图像的高度信息^[2]。但是SFS要求光源方向已知，或者加入其他的约束条件来估计光源的初始值。另外，朗伯体模型是一种理想的反射模型，然而现实中的情况大多数都不满足朗伯体模型，物体表面并非都是均匀变化的。此时对于光源参数的假设理论上是不太准确的，因此传统算法对于物体表面起伏较大，变化剧烈的情况，就不适用了。

近年来，神经网络广泛应用于各个领域^[3]。在1988年，Lehky等人^[4]首次在SFS问题中引入了神经网络，通过对多幅图像采用多层向前网络进行学习，来获取未知表面物体的形状；马颂德等人^[12]于1999年改进了Lehky的方法，利用神经网络对光源参数进行估计，只用单幅图像即可恢复未知光源参数的表面形状；Wei等人^[5]使用多层神经网络对已知光源参数的SFS问题进行了研究；俞鸿波等人^[8]利用向前神经网络建立了反射模型，对未知光源参数的SFS问题进行了讨论，但是这种方法隐含层不易确定，容易陷入局部极小值，所以还需要附加其他条件来避免网络出现局部极值行为。

RBF(radial basis function)神经网络由于其可动态决定网络结构和隐层单元的数据中心以及扩展常数、学习速度快、具有更强的生命力等优点，正在更多的领域内替代传统的神经网络。因此，本文引入RBF神经网络来建立反射模型，利用RBF的局部映射和函数逼近能力来处理SFS问题。它解决了传统算法在进行计算时，必须已知光源参数的限制。通过网络训练过程中的权值代替物体的光源信息，建立了径向基函数神经网络的反射模型。采用该模型，对表面高度$z$进行变分迭代，就可以得到最终的3维表面高度信息。这种模型可以在训练过程中自动的简化和调整权值，因此结果更接近实际的光源参数，恢复的3维效果更准确。另外，这种神经网络在自我学习和调整的过程中，对于一些在获取图像中不可避免引入的噪声也有一定的消除作用。

朗伯体(Lambertian)表面反射各个方向的光线，即所有的入射光都被均匀的散射出去^[6]。朗伯体表面的亮度与入射光的能量成正比。在入射光强度不改变的情况下，表面点的亮度不改变。朗伯体模型必须满足以下3个条件：

1) 光源为无限远处的点光源，或者均匀照明的平行光；

2) 成像几何关系为正交投影；

3) 物体表面为理想散射表面，从所有观察方向看，它都具有同样的亮度，并且完全反射所有的入射光。

物体成像后会丢失其3维信息，而SFS过程就是物体成像的一个逆过程，从2维的图像中去恢复物体的3维信息。因此物体的表面亮度与光源、物体的材质、摄像机和观察者的位置、成像几何关系等因素有着密切关系。设表面法线是$\mathit{\boldsymbol{n = }}\left( {p, q, -1} \right)$，$R\left( {p, q} \right)$表示反射模型，光源方向为${\mathit{\boldsymbol{n}}_s} = \left( {{p_s},{q_s}, - 1} \right)$，而$E\left( {x, y} \right)$表示图像的亮度，则反射等式为

$ R\left( {p, q} \right) = E\left( {x, y} \right) $

(1)

结合朗伯体余弦定律，则朗伯体反射模型可以表示为

$ R\left( {p, q} \right) = \frac{{{E_0}}}{\pi }\cos {\theta _{\rm{s}}} $

(2)

式中，${\theta _{\rm{s}}}$是表面法线与光源方向矢量夹角。

表面法线与光源方向的矢量夹角的余弦为两个矢量的点积除以向量的模的乘积，因此有

$ \begin{array}{l} \;\;\;\cos {\theta _{\rm{s}}} = \frac{{{\mathit{\boldsymbol{n}}_{\rm{s}}} \cdot \mathit{\boldsymbol{n}}}}{{\left\| {{\mathit{\boldsymbol{n}}_{\rm{s}}} \cdot \mathit{\boldsymbol{n}}} \right\|}} = \\ \frac{{p{p_{\rm{s}}} + q{q_{\rm{s}}} + 1}}{{\sqrt {{p^2} + {q^2} + 1} \sqrt {p_{\rm{s}}^2 + q_{\rm{s}}^2 + 1} }} \end{array} $

(3)

将式(3) 代入式(2) 并经过归一化可得图像某一点($x$，$y$)处的亮度与对应场景中点的表面梯度($p$，$q$)之间的对应关系

$ R\left( {p, q} \right) = \frac{{1 + p{p_{\rm{s}}} + q{q_{\rm{s}}}}}{{\sqrt {{p^2} + {q^2} + 1} \sqrt {p_{\rm{s}}^2 + q_{\rm{s}}^2 + 1} }} $

(4)

式中，$p = \frac{{\partial z}}{{\partial x}}, q = \frac{{\partial z}}{{\partial y}}$。

此时，在已知表面梯度和图像亮度的条件下，利用表面梯度与高度$z$的关系，即可通过迭代求出表面高度$z$($x$，$y$)^[7-9]。

RBF神经网络具有优秀的局部映射和函数逼近能力，是一种3层向前网络^[9]。

它具有3层结构：输入层、隐含层、输出层。其隐含层只有一层，输入层到隐含层的变换是非线性的，而从隐含层到输出层的变换则是线性的。网络的输出为隐含层单元各节点的输出的线性加权。

图 1所示为一个具有$n$个输入，隐含层具有$m$个节点以及单输出的的RBF网络。

图 1中，输入层为$\left( {{x_1}, {x_2}, {x_3} \cdots, {x_n}} \right)$，隐含层为$\left( {{\varphi _1}\left( x \right), {\varphi _2}\left( x \right), {\varphi _3}\left( x \right) \cdots {\varphi _m}\left( x \right)} \right)$，输出层为$f\left( x \right)$，隐含层与输出层之间的连接权重为$\left( {{\omega _1}, {\omega _2}, {\omega _3} \cdots {\omega _m}} \right)$。隐含层节点的作用函数称为径向基函数，通常采用非线性函数$\varphi \left( \cdot \right)$作为径向基函数^[10-11]。该函数有很多种，一般情况下采用Gauss函数，因为Gauss函数具有从中心到两边都为单调递减函数的特性，而且可以近似的去逼近任意连续函数。因此其局部映射效果更好，更适合模仿生物特征的神经网络。Gauss函数的表达式为

$ \varphi \left( x \right) = \exp \left( {-\frac{{{{\left\| {x-c} \right\|}^2}}}{{2{\sigma ^2}}}} \right) $

(5)

式中，$x$为输入层的样本，${\sigma ^2}$表示径向基函数的方差，$c$表示网络的中心，$\left\| \cdot \right\|$表示欧氏距离。

输入层将数据输入隐含层的每个节点，然后隐含层节点通过求线性加权的和就可得到该网络的输出。因此，该神经网络通过映射$f$得到的输出为

$ f\left( x \right) = {\omega _0}{\rm{ + }}\sum\limits_{i = 1}^m {{w_i}{\varphi _i}\left( x \right)} $

(6)

式中，${\omega _0}$为偏差(bias)，${\omega _i}$为隐含层到输出层之间的权值。

通过RBF神经网络的训练过程，得到隐含层与输出层节点之间的权值，进而使输出可以逼近给定的期望。从输入层到隐含层的非线性变化，可以得到隐含层到输出层的线性映射。

因为朗伯体反射模型必须满足3个假设条件，并且需要已知光源参数，最终才能建立起一组图像亮度与3维高度之间对应关系的非线性方程组。但是现实中物体往往不能满足朗伯体模型，实际情况往往更复杂，而且光源信息在大多数情况下都是未知的。因此，直接采用朗伯体模型，对于物体的3维重建效果产生了不良影响，造成结果的不准确^[12-13]。考虑到反射模型的本质是反映物体的表面方向与对应的图像灰度之间的非线性关系，为了得到更准确的结果，避免理想的简化模型带来的误差，采用RBF神经网络来建立二者之间的映射关系。

建立一个如图 2所示的RBF反射模型。在图 2中，映射函数为$f\left( \cdot \right)$，输入层为${\mathit{\boldsymbol{r}}_{i, j}}{\rm{ = }}\left( {{p_{i, j}}, {q_{i, j}}, 1} \right)$，隐含层共有$m$个节点，径向基函数$\varphi \left( \cdot \right)$选为Gauss函数，输入层到隐含层之间的连接权值为$\mathit{\boldsymbol{a = }}{\left[{{a_1}, {a_2}, {a_3}, \cdots {a_m}} \right]^{\rm{T}}}$，隐含层到输出层的权值矩阵为$\mathit{\boldsymbol{\omega = }}{\left[{{\omega _1}, {\omega _2}, {\omega _3}, \cdots, {\omega _m}, {\omega _0}} \right]^{\rm{T}}}$。

则图 2中的反射模型可以表示为

$ \begin{array}{l} {R_{i, j}}\left[{p, q} \right] = f\left( {\mathit{\boldsymbol{a}}, \mathit{\boldsymbol{\omega }}, {r_{i, j}}} \right) = \\ \;\;\;\;{\mathit{\boldsymbol{\omega }}_0} + \sum\limits_{m = 1}^m {{\omega _m}\phi \left( {{a_m}{r_{i, j}}} \right)} \end{array} $

(7)

在求解SFS问题时，大概分为最小化法、演化法、局部法和线性法等几大类经典算法。而最小化法得到的结果最稳定，效果最好。因此，文对Horn提出的最小化算法通过RBF神经网络做出改进。

引入RBF网络的Horn^[14]公式改写为

$ {E_{\rm{T}}} = \int {\int {\left[\begin{array}{l} \;\;\;{\left[{E-R\left( \mathit{\boldsymbol{r}} \right)} \right]^2} + \\ \lambda \left( {{p_{xx}}, {p_{yy}}, {q_{xx}}, {q_{yy}}} \right) \end{array} \right]} } {\rm{d}}x{\rm{d}}y $

(8)

经过有限差分和多网格对非线性项进行线性化可得参数的离散计算公式为

$ \begin{array}{l} \;\;\;{p_{i, j}}\left( {n + 1} \right) = {{\bar p}_{i, j}}\left( n \right) + \\ \frac{\varepsilon^2}{{4\lambda }}\left( {E\left( {i, j} \right)-R\left( {{\mathit{\boldsymbol{r}}_{i, j}}} \right)} \right){R_p}\left( {{\mathit{\boldsymbol{r}}_{i, j}}} \right) \end{array} $

(9)

$ \;\;\;{q_{i, j}}\left( {n + 1} \right) = {{\bar q}_{i, j}}\left( n \right) + \\ \frac{{{\varepsilon ^2}}}{{4\lambda }}\left( {E\left( {i, j} \right)-R\left( {{\mathit{\boldsymbol{r}}_{i, j}}} \right)} \right){R_q}\left( {{\mathit{\boldsymbol{r}}_{i, j}}} \right) $

(10)

$ {z_{i, j}}\left( {n + 1} \right) = {\bar z_{i, j}}-\frac{{{\varepsilon ^2}}}{4}\left( {{f_{i, j}} + {g_{i, j}}} \right) $

(11)

式中，${\bar p_{i, j}}$，${\bar q_{i, j}}$和${\bar z_{i, j}}$表示的是隐含层中相邻的4个神经元的平均值，${f_{i, j}}$和${g_{i, j}}$为$p$和$q$在局部的平均值。

因此，式(11) 可改写为离散形式

$ \begin{array}{l} \;\;\;\;{E_{\rm{T}}} = \sum\limits_{i, j \in \mathit{\boldsymbol{D}}} {{{\left[{E\left( {i, j} \right)-R\left( {{\mathit{\boldsymbol{r}}_{i, j}}} \right)} \right]}^2}} + \\ \lambda \sum {\left( {{\nabla _x}p_{i, j}^2 + {\nabla _y}p_{i, j}^2 + {\nabla _x}q_{i, j}^2 + {\nabla _y}q_{i, j}^2} \right)} \end{array} $

(12)

式中，${\nabla _x}{p_{i, j}}$表示的是对于$p$的交叉偏导，即${\nabla _x}{p_{i, j}} = \left( {{p_{i + 1, j}}-{p_{i-1, j}}} \right)/2$，其他3个参数与之意义类似。

此时，已经将SFS问题转换为求解价值函数${E_{\rm{T}}}$的最优解问题。为了避免网络体积过大、陷入局部最优解等问题，在中心${c_k}$的选择上采用正交最小二乘法，权重的计算方法采用最小均方算法(LMS)，网络的学习算法采用模拟退火算法。同时，为解决传统的RBF神经网络固定的学习率造成的收敛性慢的问题，加入最优学习率算法来根据迭代过程，在训练过程中动态地调整学习效率，加速了网络的收敛速度^[15]。

网络中的相关参数迭代计算公式为

$ \begin{array}{l} \omega \left( {n + 1} \right) = {\left[{\sum\limits_{i, j \in \mathit{\boldsymbol{D}}} {\left( {{\rho _{i, j}}\rho _{i, j}^{\rm{T}}} \right)} } \right]^{ - 1}}\\ \;\;\;\left[{\sum\limits_{i, j \in \mathit{\boldsymbol{D}}} {\left( {{E_{i, j}}\rho _{i, j}^T} \right)} } \right] + {\eta ^*}{\nabla _\omega }{E_{\rm{T}}} \end{array} $

(13)

$ \begin{array}{l} \;\;{c_k}\left( {n + 1} \right) = {c_k}\left( n \right) + \\ {\eta _c}\left[{-\frac{{\partial {E_{\rm{T}}}}}{{\partial {c_k}}} + \nabla {E_{pen}}\left( {{c_k}} \right)} \right] \end{array} $

(14)

$ \begin{array}{l} \;\;{\sigma _k}\left( {n + 1} \right) = {\sigma _k}\left( n \right) + \\ {\eta _c}\left[{-\frac{{\partial {E_{\rm{T}}}}}{{\partial {\sigma _k}}} + \nabla {E_{pen}}\left( {{\sigma _k}} \right)} \right] \end{array} $

(15)

$ {\eta ^*} = \frac{{{\mathit{\boldsymbol{e}}^{\rm{T}}}\rho {\rho ^{\rm{T}}}\mathit{\boldsymbol{e}}}}{{{\mathit{\boldsymbol{e}}^{\rm{T}}}\rho {\rho ^{\rm{T}}}{\rho ^{\rm{T}}}\rho \mathit{\boldsymbol{e}}}} $

(16)

$ \mathit{\boldsymbol{e}} = \mathit{\boldsymbol{\hat y}}-\mathit{\boldsymbol{y}} $

(17)

$ \begin{array}{l} \;\;\;\;\;{\rho _{i, j}} = \left( {1, \phi \left( {\left\| {{a_{i, j}}-{c_1}} \right\|} \right)} \right., \\ \phi \left( {\left\| {{a_{i, j}}-{c_2}} \right\|} \right), \; \cdots, \phi \left( {\left\| {{a_{i, j}}-{c_m}} \right\|} \right) \end{array} $

(18)

式中，$\mathit{\boldsymbol{D}}$为图像像素点的集合，${\eta _c}$和${\eta _\sigma }$为RBF神经网络的中心和半径的学习效率，${\eta ^*}$为学习效率的修正因子，$\hat y$为网络计算的输出值，$y$为实际输出值，$\nabla {E_{pen}}$为最优化过程中的处罚因子。

实验中，最终学习过程的停止标志设为$E\left( {i, j} \right)$与$R\left( {{\mathit{\boldsymbol{r}}_{i, j}}} \right)$之间的误差小于设定值或者达到最大迭代次数。

针对SFS算法处理的经典合成图像以及实际图像做了3维重建工作。图 3为两幅合成图像和两幅实际图像：图 4—图 8分别为传统的SFS算法恢复的3维图像与本文改进算法所恢复的3维图像以及文献[2, 6]中的恢复算法处理后的3维图像效果对比。图 9为地形图放大后截取右上角部分区域的对比图。

从几幅图像的实验结果来看，传统的SFS算法对于表面光滑的图像，其处理后恢复的小效果相对于表面不太光滑的效果更好。传统算法对图 4的结果中图 4(a)显示的花瓶3维效果图最好，图 5(a)是对莫扎特的3维重建，在鼻子的部分出现了严重缺陷，头部的轮廓效果也不太好；图 6(a)是对古文物的3维重建图，恢复了文物的基本形状，但是文物花纹的处理就出现了误差，较为模糊；图 7 (a)为对实际图像地形图的3维重建，尖峰点特别明显，只能分辨出图像中颜色有明显差异的部分，而对于山区的形状信息，没有重构出来。因此，传统的SFS算法对于合成图像效果尚可，但对实际图像就产生了严重的错误。

本文对于文献[2]和文献[6]中的算法进行了实验，并与传统算法和本文改进后的算法进行了对比。从结果来看，本文算法对于每幅图像的处理结果都得到了优于其他3种方法的性能。图 4中花瓶的3维重建图，图 4 (b)较为平滑，图 4 (c)在花瓶的颈部出现了一些扭曲错误，图 4 (d)的视觉效果最好，整体表面光滑，视觉效果最佳；图 5对于莫扎特图像的恢复，图 5 (b)出现了过于平滑的问题，面部特征几乎不能识别，图 5 (c)基本恢复了面部特征，但是细节处无法清晰的辨认，图 5 (d)不但恢复了面部特征，而且对于细节处的处理更好，五官更加细腻；图 6中，图 6 (b)对于文物的花纹出现了较多的尖峰点，图 6 (c)的结果3维效果很差，但是形状恢复比较好，图 6 (d)兼顾了3维视觉效果和文物的真实形状；图 7(b)依然出现了过于平滑的问题，其效果不佳，图 7 (c)出现了严重错误，基本不能识别图像的形状信息，图 7 (d)中山脉的基本形状清晰可见，下方的河流也基本可以识别。图 8为地形图右上角的部分区域经过放大后的对比，可以看到，图 8 (d)中包含着更多的形状信息，其效果更好。

本文对于文献[2]和文献[6]中的算法进行了实验，并与传统算法和本文改进后的算法进行了对比。从结果来看，本文算法对于每幅图像的处理结果都得到了优于其他3种方法的性能。图 4中花瓶的3维重建图，图 4 (b)较为平滑，图 4 (c)在花瓶的颈部出现了一些扭曲错误，图 4 (d)的视觉效果最好，整体表面光滑，视觉效果最佳；图 5对于莫扎特图像的恢复，图 5 (b)出现了过于平滑的问题，面部特征几乎不能识别，图 5 (c)基本恢复了面部特征，但是细节处无法清晰的辨认，图 5 (d)不但恢复了面部特征，而且对于细节处的处理更好，五官更加细腻；图 6中，图 6 (b)对于文物的花纹出现了较多的尖峰点，图 6 (c)的结果3维效果很差，但是形状恢复比较好，图 6 (d)兼顾了3维视觉效果和文物的真实形状；图 7(b)依然出现了过于平滑的问题，其效果不佳，图 7 (c)出现了严重错误，基本不能识别图像的形状信息，图 7 (d)中山脉的基本形状清晰可见，下方的河流也基本可以识别。图 8为地形图右上角的部分区域经过放大后的对比，可以看到，图 8 (d)中包含着更多的形状信息，其效果更好。

综合以上实验结果及分析可以发现，传统算法对于表面光滑的物体，其3维恢复效果较好，而对于表面变化剧烈的物体，就会产生严重错误。并且，传统算法适用于合成图像，对于实际图像，其误差不可控。而文献[2]的方法对于图像3维恢复产生了过于平滑的问题，对于表面变化较大的物体，过于平滑导致物体的形状信息失真严重，并且该方法对于实际图像的效果很差。文献[6]的方法在进行3维重建时，对于3维形状信息恢复较好，但在视觉效果方面差强人意，实验中的花瓶图像颈部产生了较大的扭曲，莫扎特脸部细节丢失，古文物视觉效果差，实际地形图出现严重错误，形状信息基本丢失。而本文改进后的算法对于表面平滑和变化剧烈的物体同样适用，都表现了优异的性能。不但在形状的恢复上细节处理更加到位，对于3维视觉效果方面也明显优于其他3种算法。

本文所改进的算法在两种图像的视觉效果都展现了优于传统算法的性能，下面在恢复的图像高度误差方面给出了更为定量的对比。选取表面变化最大的地形图作为定量对比的实验数据来源。随机提取了该地形图像中128个像素点的真实高度值，然后提取传统算法和本文改进算法处理后的图像中对应的128个像素点的灰度值。为更加易于对比高度的误差，对真实的高度值和处理后的灰度值采用min-max函数进行了归一化处理，将高度值归一化到0和1之间，然后求两个高度值的绝对误差，并做出其折线图。图 9显示了其误差范围。

从图 9中可以看出，传统的SFS算法误差范围较大，最大值达到0.987 5，而且误差不稳定，这源于其对光源的估计值不准确，而本文的算法误差范围小，只有极个别点达到0.3，最大误差比传统算法提升了60 %以上。算法的误差范围波动较为稳定，因为改进后的算法经过训练过程，对于光源的估计更准确，因此在误差和稳定性方面的提升显而易见。

由于本文算法中添加了调整学习率的最优学习算法，下面给出固定学习率和自动调整学习率下的算法运行时间对比。对于一幅128×128像素的图像，在Matlab2014a的软件平台下，其运行时间如表 1所示。

从表 1中的时间对比可以看出，由于加入了神经网络，所以相对于传统的SFS算法，无论是固定学习率还是自调整学习率，所花费的时间都要比传统方法多。但是，自调整的学习率大大缩短了训练过程，耗时不到固定学习率的一半，提高了RBF神经网络的收敛速度。

传统的SFS算法采用的理想反射模型—朗伯体模型在实际中难以满足需要，对光源估计不准确造成恢复的3维图像误差较大，尤其对于物体表面高度变化剧烈的图像误差更大，3维效果很差。本文引入RBF神经网络模型，针对该问题，建立了基于RBF网络的反射模型。然后基于该模型，对传统的SFS算法进行了改进。改进后的算法克服了传统方法对于光源估计不准确的缺点，而且对于合成图像和实际图像都展现了优秀的3维恢复能力，经过改进的算法在物体表面高度值变化剧烈的情况下仍然适用。实验结果表明，本文算法性能优越，而且效率高，耗时少，其稳定性得到了很大的提升。

本文改进的算法对于合成图像和实际图像都表现了优秀的性能，适用于多数非高精度场合的物体表面3维重建工作，不需要初始光源方向即可进行物体的表面3维恢复。该算法仍然存在某些网络初值需要依靠经验选取的问题，网络的训练时间和训练过程收到初值的影响，下一步可以对初值的选取进行建模，以获取更加准确的结果。