CN102819866B - 一种非朗伯表面快速三维重构方法 - Google Patents

Abstract

本发明属于图像理解与计算机视觉技术领域，特别涉及一种非朗伯表面快速三维重构方法。提供的技术方案是：一种非朗伯表面快速三维重构方法，依次包括下述步骤：一．构建适于非朗伯表面且接近实际成像条件下的成像模型；二．由上述成像模型和摄像机获取的灰度图像，建立非朗伯表面下SFS问题的图像辐照度方程；三．利用粘性解理论计算上述图像辐照度方程的粘性解；四．根据图像辐照度方程解的结果，重构出非朗伯表面的三维形状信息。本发明的优点是：直接对非朗伯表面下建立的图像辐照度方程进行求解，同时利用弱解理论——粘性解理论计算图像辐照度方程的粘性解，使非朗伯表面的SFS问题良性化；具有快速、准确的优点。

Description

一种非朗伯表面快速三维重构方法

技术领域

本发明属于图像理解与计算机视觉技术领域，主要涉及到由单幅图像明暗变化实现三维重构的技术，特别涉及到一种非朗伯表面快速三维重构方法。

背景技术

为了满足工农业产品检测、文档图像复原、医学图像分析、生物特征识别以及地形、地貌测量等领域中三维建模的需求，有效重构物体表面的三维形状信息，近年来已经成为图像理解与计算机视觉技术领域中一个非常重要的研究方向。然而，如何快速、准确重构非朗伯表面的三维形状，至今仍未得到较好解决。

目前，三维形状信息的获取手段主要有两大类：接触式测量和非接触式测量。尽管接触式测量方法已经比较成熟，而且还可以获得很高精度的数据，但其本身也存在较多不足：①接触式测量使得测头容易磨损及划伤被测表面，测量过程需要人工干预，因而测量成本较高，自动化程度较低；②接触式测量采用逐点扫描的方法，且扫描速度受到机械运动的限制，故测量时间较长、测量效率较低等，不太适于快速测量的场合；③对大型物件，尤其是软质材料测量效果不好，对测头不能触及的表面无法测量，而且对测量环境要求较高，故其应用范围往往受到限制。总之，接触式测量难以满足快速、有效的测量需求。

在申请号为03153504.6的中国发明专利中，申请人提到了非接触式测量的几种方法：光学传感器法、激光扫描法、立体视觉法、投影栅相位法等。随着计算机视觉和图像技术的发展，利用图像数据获取目标三维形状信息的非接触式测量技术已在社会生产的各个方面显示出越来越重要的地位和作用。由单幅图像明暗变化恢复形状(简称SFS)技术是实现物体表面三维重构的关键方法之一，其目标是利用单幅图像的明暗变化恢复物体表面各点的法向向量，进而可以获得其相对高度。

许多研究者试图找出最有效的实现SFS技术的方法，并取得了一定的进展。目前应用最广泛的主要有四种方法：演化和偏微分方程方法、最优化方法、局部化方法以及线性化方法等。演化和偏微分方程方法是从图像中一组已知高度值的边界点或奇异点出发，逐步确定出所有图像像素点对应的物体表面的形状信息，它是对SFS问题建立的一阶非线性偏微分方程本身直接求解。最优化方法就是将SFS问题中的图像辐照度方程表示为能量函数形式，然后通过附加约束条件将其转化为求解泛函极值的问题。局部化方法的思路就是在物体表面的每一点处，假定表面形状为平面、球面或抛物面等，并将表面反射模型与假设的局部形状结合起来，构成关于局部形状参数的线性偏微分方程组，然后利用已知的边界条件来求解该方程组得到物体表面的三维形状。线性化方法就是通过将反射图线性化，使非线性问题转化为线性问题，进而求解得到图像辐照度方程的解。

在上述这些实现SFS技术的方法中，各种方法存在一定的优点，也不可避免地存在不足之处。如线性化方法速度较快，但重构精度不高，求解过程缺乏可信度；局部化方法使用范围受到限制，且精度不高；演化和偏微分方程方法虽然精度较高，但需要考虑边界条件；最优化方法不需要考虑边界条件，但需要引入各种约束条件等等。

同时，上述方法通常只是针对朗伯表面设计的，如果应用到非朗伯表面，往往误差较大。

另外，SFS方法中建立的图像辐照度方程通常是Hamilton-Jacobi类型的一阶非线性偏微分方程，一般情况下不存在经典意义下的光滑解，通常认为是病态问题。为了使SFS问题良性化，考虑图像辐照度方程在弱解意义下的特性。Hamilton-Jacobi方程粘性解理论是一种非常好的弱解理论，可定量、直观地描述一阶或二阶偏微分方程的弱解。法国学者Prados提出了使用粘性解理论对朗伯表面下建立的图像辐照度方程解的特性进行了分析。但遗憾的是，对于非朗伯表面建立的图像辐照度方程的研究，尚未发现国内外有相关文献报导。

发明内容

本发明的目的在于针对现有技术的不足，提出一种非朗伯表面快速三维重构方法。本发明所述的方法不需要特定的边界条件，而且使SFS问题良性化，具有快速、准确的优点。

为了实现上述目的，本发明采用的技术方案如下：

一种非朗伯表面快速三维重构方法，依次包括下述步骤：

一.构建适于非朗伯表面且接近实际成像条件下的成像模型：从三个阶段入手，光源入射过程、非朗伯表面对入射光的反射过程、摄像机对反射光的记录过程，涉及到光源的分布、非朗伯表面的反射特性以及摄像机的投影方式的建模，其中：所述光源为近点光源，而且光源辐照度随着物体与光源之间距离的增大而存在衰减；非朗伯表面的反射特性使用Oren-Nayer反射模型来近似；摄像机的投影方式采用透视投影方式；

二.由上述成像模型和摄像机获取的灰度图像，建立非朗伯表面下SFS问题的图像辐照度方程：该步骤里依次包括下述步骤，

(一)建立三维笛卡尔直角坐标系，设成像平面为x＝(x₁，x₂)位于Z＝-f处，f＞0为摄像机的焦距，摄像机的光轴与Z轴重合，光源位于投影中心，即光心O处,

非朗伯表面S可以表示为Z＝Z(X)，X＝(X₁，X₂)，像平面上一点(x，-f)是物体表面S上点P＝(X₁，X₂，Z)在成像平面的投影点，r为光心O到P点的距离，根据透视投影成像原理有

$\frac{X_1}{x_1}=\frac{X_2}{x_2}=\frac{Z}{-f}=\frac{r}{{({\left|\right|x\left|\right|}^2+f^2)}^{1/2}}---\left(1\right)$

令r＝u(x)≥0，于是物体表面S可以用函数表示

$S=\{S\left(x\right)=\frac{u\left(x\right)}{{({\left|\right|x\left|\right|}^2+f^2)}^{1/2}}\left(\begin{array}{c}x\\ -f\end{array}\right),x=\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\∈\stackrel{\‾}{\mathrm{\Ω}}\}---\left(2\right)$

其中：Ω是定义在实数集上的一个开集，代表图像的区域,

(二)计算光源入射到S(x)点的单位方向向量为

$L\left(x\right)=\frac{1}{{({\left|\right|x\left|\right|}^2+f^2)}^{1/2}}\left(\begin{array}{c}-x\\ f\end{array}\right)---\left(3\right)$

(三)利用Oren-Nayer反射模型计算非朗伯表面辐射亮度L_s

$L_s({\mathrm{\θ}}_i,{\mathrm{\φ}}_i;{\mathrm{\θ}}_r,{\mathrm{\φ}}_r;\mathrm{\σ})=\frac{I_0}{r^2}\frac{\mathrm{\ρ}}{\mathrm{\π}}\cos {\mathrm{\θ}}_i\×(A+\mathrm{M\; m\; ax}[0,\cos ({\mathrm{\φ}}_r-{\mathrm{\φ}}_i)\left]\sin \mathrm{\α}\tan \mathrm{\β}\right)---\left(4\right)$

其中：θ_i，φ_i和θ_r，φ_r是非朗伯表面某一点的局部坐标系，分别为光源方向向量L和摄像机方向向量V的天顶角和方位角；σ为Gauss分布函数的标准差，代表了表面的粗糙程度；ρ为漫反射率；I₀为点光源的辐射强度；

$A=1-0.5\frac{{\mathrm{\σ}}^2}{{\mathrm{\σ}}^2+0.33};B=0.45\frac{{\mathrm{\σ}}^2}{{\mathrm{\σ}}^2+0.09};$ α＝max[θ_i，θ_r]；β＝min[θ_i，θ_r]，

由于光源位于投影中心，故θ_i＝θ_r＝α＝β，φ_i＝φ_r，式(4)简化为

$I\left(x\right)=\frac{1}{r^2}(A\cos {\mathrm{\θ}}_i+Bs{\mathrm{in}}^2{\mathrm{\θ}}_i)---\left(5\right)$

其中：I(x)为摄像机获取的灰度图像；r＝fu(x)，

(四)计算摄像机遵循透视投影方式时，非朗伯表面点S(x)处的法向向量n(x)为

$n\left(x\right)=\left(\begin{array}{c}f\▿u\left(x\right)-\frac{\mathrm{fu}\left(x\right)}{{\left|\right|x\left|\right|}^2+f^2}x\\ x\·\▿u\left(x\right)+\frac{\mathrm{fu}\left(x\right)}{{\left|\right|x\left|\right|}^2+f^2}f\end{array}\right)---\left(6\right)$

(五)由于θ_i为n(x)与L(x)之间的夹角，且令v(x)＝ln u(x)，故有

${\mathrm{\θ}}_i=\arccos [\frac{n\left(x\right)}{\left|\right|n\left(x\right)\left|\right|}\·L\left(x\right)]=\arccos \left[\frac{Q\left(x\right)}{{(f^2{\left|\right|\▿v\left(x\right)\left|\right|}^2+{(x\·\▿\left(x\right))}^2+Q{\left(x\right)}^2)}^{1/2}}\right]---\left(7\right)$

Q(x)＝f/(||x||²+f²)^1/2.

其中：

将式(7)代入式(5)，得到非朗伯表面下SFS问题的图像辐照度方程

$I\left(x\right)=\frac{1}{f^2{e}^{2v}}[A\frac{Q\left(x\right)}{F(x,\▿v)}+B\frac{F{(x,\▿v)}^2-Q{\left(x\right)}^2}{F{(x,\▿v)}^2}]---\left(8\right)$

其中：I(x)为摄像机获取的灰度图像；

三.利用粘性解理论计算上述图像辐照度方程的粘性解：该步骤依次包括下述步骤

(一)将图像辐照度方程(8)看作是关于的一元二次方程

$[I\left(x\right)f^2{e}^{2v}-B]F{(x,\▿v)}^2-\mathrm{AQ}\left(x\right)F(x,\▿v)+\mathrm{BQ}{\left(x\right)}^2=0---\left(9\right)$

(二)求解一元二次方程(9)得

$F(x,\▿v)=\frac{A+\sqrt{A^2-4[I\left(x\right)f^2{e}^{2v}-B]B}}{2[I\left(x\right)f^2{e}^{2v}-B]}Q\left(x\right)---\left(10\right)$

(三)于是，图像辐照度方程(8)等价为

$-{(A-\sqrt{A^2-4[I\left(x\right)f^2{e}^{2v}-B]B})}^{-1}+{\left[2\mathrm{BQ}\left(x\right)\right]}^{-1}\sqrt{f^2{\left|\right|\▿v\left(x\right)\left|\right|}^2+Q{\left(x\right)}^2}=0---\left(11\right)$

(四)利用粘性解理论计算上述图像辐照度方程的解,

显然，上式是一个Hamilton-Jacobi类型的一阶非线性PDE，将简记为p，可以得到相应的Hamilton函数

$H(x,v,p)=-{(A-\sqrt{A^2-4[I\left(x\right)f^2{e}^{2v}-B]B})}^{-1}+{\left[2\mathrm{BQ}\left(x\right)\right]}^{-1}\sqrt{f^2{\left|\right|p\left|\right|}^2+{(x\·p)}^2+Q{\left(x\right)}^2}---\left(12\right)$

(五)应用最优控制理论，式(12)可以转化为如下的控制形式

$H(x,v,p)=-{(A-\sqrt{A^2-4[I\left(x\right)f^2{e}^{2v}-B]B})}^{-1}+\underset{a\∈{\stackrel{\‾}{B}}_2\left(\mathrm{0,1}\right)}{\sup }\{-f_c(x,a)\·p-l_c(x,a)\}---\left(13\right)$

其中：f_c(x，a)＝-[2BQ(x)]^-1R^TDR a；l_c(x，a)＝-1/(2B)(1-||a||²)^1/2；是定义在上的单位圆面；D和R满足

$D=\left(\begin{array}{cc}{({\left|\right|x\left|\right|}^2+f^2)}^{1/2}& 0\\ 0& f\end{array}\right),R=\left\{\begin{array}{cc}\left(\begin{array}{cc}{x}_1/\left|\right|x\left|\right|& x_2/\left|\right|x\left|\right|\\ x_2/\left|\right|x\left|\right|& -x_1/\left|\right|x\left|\right|\end{array}\right)& \left|\right|x\left|\right|\≠0\\ \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)& \left|\right|x\left|\right|=0\end{array}\right.$

(六)图像区域Ω＝(1，m)×(1，n)的均匀离散网格点：x_i，j＝(ih₁，jh₂)，i＝1，2，…，m，j＝1，2，…，n，其中：h＝(h₁，h₂)定义了数值算法中离散网格的大小,

定义对其在时域应用前向Euler公式得到图像辐照度方程求解算法

$V_{i,j}^{n+1}=v_{i,j}^n-\mathrm{\ΔtH}(x_{i,j},V_{i,j},p)---\left(14\right)$

(七)使用如下方法逼近H(x_i，j，V_i，j，p)

$\begin{array}{c}H(x_{i,j},V_{i,j},p)\≈-{(A-\sqrt{A^2-4[I\left(x_{i,j}\right)f^2{e}^{{2V}_{i,j}}-B]B})}^{-1}+\underset{a\∈{\stackrel{\‾}{B}}_2\left(\mathrm{0,1}\right)}{\sup }\left\{\max \right[-f_1(x_{i,j},a),0]p_1^{-}\\ +\min [{-f}_1(x_{i,j},a),0]p_1^{+}+\max [{-f}_2(x_{i,j},a),0]p_2^{-}+\min [{-f}_2(x_{i,j},a),0]p_2^{+}-l_c(x_{i,j},a)\}\end{array}---\left(15\right)$

其中：f₁(x_i，j，a)、f₂(x_i，j，a)为f_c(x，a)的分量；分别为p的后向、前向差分,

(八)对于逼近方法(15)，计算以下最优问题：

$M=\underset{a\∈{\stackrel{\‾}{B}}_2\left(\mathrm{0,1}\right)}{\sup }\{-\left({-\left[2\mathrm{BQ}\left(x_{i,j}\right)\right]}^{-1}{R}^T\mathrm{DR\; a}\right)\·p+1/2B\sqrt{1-{\left|\right|a\left|\right|}^2}\}---\left(16\right)$

(九)求解算法(14)转化为

$V_{i,j}^{n+1}-\mathrm{\Δt}{(A-\sqrt{A^2-4[I\left(x_{i,j}\right)f^2{e}^{{2V}_{i,j}^{n+1}}-B]B})}^{-1}=V_{i,j}^n-\mathrm{\ΔtM}---\left(17\right)$

使用Newton法来进行求解，得到图像辐照度方程的粘性解,

(十)使用迭代Fast Marching对上述求解算法(17)进行加速收敛,

四.根据图像辐照度方程解的结果，重构出非朗伯表面的三维形状信息。

与现有技术相比，本发明的优点是：

1.本发明可实现由单幅图像明暗变化重构非朗伯表面的三维形状，避免了多幅图像之间的匹配问题；

2.本发明考虑光源辐照度随着非朗伯表面与光源之间距离的增大而存在衰减的影响，同时使用Oren-Nayer反射模型描述非朗伯表面的反射特性，采用透视投影方式近似摄像机的投影方式，避免了由于成像模型不准确导致的重构精度不高的问题，因此可达到准确重构的需求；

3.本发明直接对非朗伯表面下建立的图像辐照度方程进行求解，避免了需要引入约束条件以及由于数值算法不合适导致的重构精度不高的问题；

4.本发明利用弱解理论——粘性解理论计算图像辐照度方程的粘性解，克服了传统SFS方法的病态性，从而使非朗伯表面的SFS问题良性化，因此进一步保障了准确重构的需求；

5.本发明不需要特定的边界条件，这是目前其他方法所无法实现的，从而使得本发明具有更广阔的适用范围；

6.本发明使用迭代Fast Marching对求解数值算法进行加速收敛，因而CPU运行时间较少，使本发明具有快速的特点。

附图说明：

图1是本发明的流程图；

图2是本发明建立的三维笛卡尔直角坐标系；

图3是本发明中非朗伯表面某一点的局部坐标系；

图4是非朗伯半球表面图像；

图5是非朗伯花瓶表面图像；

图6是本发明重构的非朗伯半球表面的三维形状；

图7是本发明重构的非朗伯花瓶表面的三维形状；

图8是非朗伯半球表面的真实高度值；

图9是非朗伯花瓶表面的真实高度值；

图10是本发明重构的非朗伯半球表面的高度误差；

图11是本发明重构的非朗伯花瓶表面的高度误差。

具体实施方式：

下面结合附图和实施例，对本发明做进一步详细说明，本实施例以本发明的技术方案为前提进行实施，给出了详细的实施方式，但本发明的保护范围不限于下述的实施例。

参照附图1，本发明所述非朗伯表面快速三维重构方法包括如下步骤：

一.构建适于非朗伯表面且接近实际成像条件下的成像模型

构建成像模型的方法为：从三个阶段入手：光源入射过程、非朗伯表面对入射光的反射过程、摄像机对反射光的记录过程，涉及到光源的分布、非朗伯表面的反射特性以及摄像机的投影方式的建模。本方法中，光源为近点光源，而且光源辐照度随着物体与光源之间距离的增大而存在衰减；非朗伯表面的反射特性使用Oren-Nayer反射模型来近似；摄像机的投影方式采用透视投影方式。

二.由上述成像模型和摄像机获取的灰度图像，建立非朗伯表面下SFS问题的图像辐照度方程，具体步骤是：

(一)建立如附图2所示的三维笛卡尔直角坐标系，设成像平面为x＝(x₁，x₂)位于Z＝-f处，f＞0为摄像机的焦距，摄像机的光轴与Z轴重合，光源位于投影中心，即光心O处。

非朗伯表面S可以表示为Z＝Z(X)，X＝(X₁，X₂)，像平面上一点(x，-f)是物体表面S上点P＝(X₁，X₂，Z)在成像平面的投影点，r为光心O到P点的距离，根据透视投影成像原理有

$\frac{X_1}{x_1}=\frac{X_2}{x_2}=\frac{Z}{-f}=\frac{r}{{({\left|\right|x\left|\right|}^2+f^2)}^{1/2}}---\left(1\right)$

令r＝u(x)≥0，于是物体表面S可以用函数表示

$S=\{S\left(x\right)=\frac{u\left(x\right)}{{({\left|\right|x\left|\right|}^2+f^2)}^{1/2}}\left(\begin{array}{c}x\\ -f\end{array}\right),x=\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\∈\stackrel{\‾}{\mathrm{\Ω}}\}---\left(2\right)$

其中：Ω是定义在实数集上的一个开集，代表图像的区域。

(二)计算光源入射到S(x)点的单位方向向量为

$L\left(x\right)=\frac{1}{{({\left|\right|x\left|\right|}^2+f^2)}^{1/2}}\left(\begin{array}{c}-x\\ f\end{array}\right)---\left(3\right)$

(三)利用Oren-Nayer反射模型计算非朗伯表面辐射亮度L_s

$L_s({\mathrm{\θ}}_i,{\mathrm{\φ}}_i;{\mathrm{\θ}}_r,{\mathrm{\φ}}_r;\mathrm{\σ})=\frac{I_0}{r^2}\frac{\mathrm{\ρ}}{\mathrm{\π}}\cos {\mathrm{\θ}}_i\×(A+\mathrm{M\; m\; ax}[0,\cos ({\mathrm{\φ}}_r-{\mathrm{\φ}}_i)\left]\sin \mathrm{\α}\tan \mathrm{\β}\right)---\left(4\right)$

其中：θ_i，φ_i和θ_r，φ_r是非朗伯表面某一点的局部坐标系，分别为光源方向向量L和摄像机方向向量V的天顶角和方位角；σ为Gauss分布函数的标准差，代表了表面的粗糙程度；ρ为漫反射率；I₀为点光源的辐射强度；

$A=1-0.5\frac{{\mathrm{\σ}}^2}{{\mathrm{\σ}}^2+0.33};B=0.45\frac{{\mathrm{\σ}}^2}{{\mathrm{\σ}}^2+0.09};$ α＝max[θ_i，θ_r]；β＝min[θ_i，θ_r]。

本发明中由于光源位于投影中心，故θ_i＝θ_r＝α＝β，φ_i＝φ_r，式(4)简化为

$I\left(x\right)=\frac{1}{r^2}(A\cos {\mathrm{\θ}}_i+Bs{\mathrm{in}}^2{\mathrm{\θ}}_i)---\left(5\right)$

其中：I(x)为摄像机获取的灰度图像；r＝fu(x)。

(四)计算摄像机遵循透视投影方式时，非朗伯表面点S(x)处的法向向量n(x)为

$n\left(x\right)=\left(\begin{array}{c}f\▿u\left(x\right)-\frac{\mathrm{fu}\left(x\right)}{{\left|\right|x\left|\right|}^2+f^2}x\\ x\·\▿u\left(x\right)+\frac{\mathrm{fu}\left(x\right)}{{\left|\right|x\left|\right|}^2+f^2}f\end{array}\right)---\left(6\right)$

(五)由于θ_i为n(x)与L(x)之间的夹角，且令v(x)＝ln u(x)，故有

${\mathrm{\θ}}_i=\arccos [\frac{n\left(x\right)}{\left|\right|n\left(x\right)\left|\right|}\·L\left(x\right)]=\arccos \left[\frac{Q\left(x\right)}{{(f^2{\left|\right|\▿v\left(x\right)\left|\right|}^2+{(x\·\▿\left(x\right))}^2+Q{\left(x\right)}^2)}^{1/2}}\right]---\left(7\right)$

其中：Q(x)＝f/(||x||²+f²)^1/2.。将式(7)代入式(5)，得到非朗伯表面下SFS问题的图像辐照度方程

$I\left(x\right)=\frac{1}{f^2{e}^{2v}}[A\frac{Q\left(x\right)}{F(x,\▿v)}+B\frac{F{(x,\▿v)}^2-Q{\left(x\right)}^2}{F{(x,\▿v)}^2}]---\left(8\right)$

其中：I(x)为摄像机获取的灰度图像；

$F(x,\▿v)={(f^2{\left|\right|\▿v\left(x\right)\left|\right|}^2+{(x\·\▿v\left(x\right))}^2+Q{\left(x\right)}^2)}^{1/2}.$

三.利用粘性解理论计算上述图像辐照度方程的粘性解：该步骤依次包括下述步骤

(一)将图像辐照度方程(8)看作是关于的一元二次方程

$[I\left(x\right)f^2{e}^{2v}-B]F{(x,\▿v)}^2-\mathrm{AQ}\left(x\right)F(x,\▿v)+\mathrm{BQ}{\left(x\right)}^2=0---\left(9\right)$

(二)求解一元二次方程(9)得

$F(x,\▿v)=\frac{A+\sqrt{A^2-4[I\left(x\right)f^2{e}^{2v}-B]B}}{2[I\left(x\right)f^2{e}^{2v}-B]}Q\left(x\right)---\left(10\right)$

(三)于是，图像辐照度方程(8)等价为

$-{(A-\sqrt{A^2-4[I\left(x\right)f^2{e}^{2v}-B]B})}^{-1}+{\left[2\mathrm{BQ}\left(x\right)\right]}^{-1}\sqrt{f^2{\left|\right|\▿v\left(x\right)\left|\right|}^2+Q{\left(x\right)}^2}=0---\left(11\right)$

(四)一般情况下，方程(11)不存在经典解，可以利用粘性解理论计算上述图像辐照度方程的解。显然，上式是一个Hamilton-Jacobi类型的一阶非线性PDE，将简记为p，可以得到相应的Hamilton函数

$H(x,v,p)=-{(A-\sqrt{A^2-4[I\left(x\right)f^2{e}^{2v}-B]B})}^{-1}+{\left[2\mathrm{BQ}\left(x\right)\right]}^{-1}\sqrt{f^2{\left|\right|p\left|\right|}^2+{(x\·p)}^2+Q{\left(x\right)}^2}---\left(12\right)$

(五)应用最优控制理论，式(12)可以转化为如下的控制形式

$H(x,v,p)=-{(A-\sqrt{A^2-4[I\left(x\right)f^2{e}^{2v}-B]B})}^{-1}+\underset{a\∈{\stackrel{\‾}{B}}_2\left(\mathrm{0,1}\right)}{\sup }\{-f_c(x,a)\·p-l_c(x,a)\}---\left(13\right)$

其中：f_c(x，a)＝-[2BQ(x)]^-1R^TDR a；l_c(x，a)＝-1/(2B)(1-||a||²)^1/2；是定义在上的单位圆面；D和R满足

$D=\left(\begin{array}{cc}{({\left|\right|x\left|\right|}^2+f^2)}^{1/2}& 0\\ 0& f\end{array}\right),R=\left\{\begin{array}{cc}\left(\begin{array}{cc}{x}_1/\left|\right|x\left|\right|& x_2/\left|\right|x\left|\right|\\ x_2/\left|\right|x\left|\right|& -x_1/\left|\right|x\left|\right|\end{array}\right)& \left|\right|x\left|\right|\≠0\\ \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)& \left|\right|x\left|\right|=0\end{array}\right.$

(六)图像区域Ω＝(1，m)×(1，n)的均匀离散网格点：x_i，j＝(ih₁，jh₂)，i＝1，2，…，m，j＝1，2，…，n，其中：h＝(h₁，h₂)定义了数值算法中离散网格的大小。定义对其在时域应用前向Euler公式得到图像辐照度方程求解算法

$V_{i,j}^{n+1}=v_{i,j}^n-\mathrm{\ΔtH}(x_{i,j},V_{i,j},p)---\left(14\right)$

(七)本发明中使用如下方法逼近H(x_i，j，V_i，j，p)

$\begin{array}{c}H(x_{i,j},V_{i,j},p)\≈-{(A-\sqrt{A^2-4[I\left(x_{i,j}\right)f^2{e}^{{2V}_{i,j}}-B]B})}^{-1}+\underset{a\∈{\stackrel{\‾}{B}}_2\left(\mathrm{0,1}\right)}{\sup }\left\{\max \right[-f_1(x_{i,j},a),0]p_1^{-}\\ +\min [{-f}_1(x_{i,j},a),0]p_1^{+}+\max [{-f}_2(x_{i,j},a),0]p_2^{-}+\min [{-f}_2(x_{i,j},a),0]p_2^{+}-l_c(x_{i,j},a)\}\end{array}---\left(15\right)$

其中：f₁(x_i，j，a)、f₂(x_i，j，a)为f_c(x，a)的分量；分别为p的后向、前向差分。

(八)对于逼近方法(15)，计算以下最优问题：

$M=\underset{a\∈{\stackrel{\‾}{B}}_2\left(\mathrm{0,1}\right)}{\sup }\{-\left({-\left[2\mathrm{BQ}\left(x_{i,j}\right)\right]}^{-1}{R}^T\mathrm{DR\; a}\right)\·p+1/2B\sqrt{1-{\left|\right|a\left|\right|}^2}\}---\left(16\right)$

(九)求解算法(14)转化为

$V_{i,j}^{n+1}-\mathrm{\Δt}{(A-\sqrt{A^2-4[I\left(x_{i,j}\right)f^2{e}^{{2V}_{i,j}^{n+1}}-B]B})}^{-1}=V_{i,j}^n-\mathrm{\ΔtM}---\left(17\right)$

使用Newton法来进行求解，得到图像辐照度方程的粘性解。

(十)使用迭代Fast Marching对上述求解算法(17)进行加速收敛。

四.根据图像辐照度方程解的结果，重构出非朗伯表面的三维形状信息。

本发明效果可以通过以下实验进一步证实：

一.实验条件：

实验所使用的输入图像如图4和图5所示，图4是非朗伯半球表面图像，图5是非朗伯花瓶表面图像。实验中的方法使用Matlab结合C语言mex函数实现，其中输入输出部分使用Matlab实现，图像辐照度方程求解算法使用C语言实现。

二.实验结果：

为了直观说明本发明的性能，实验中首先用本发明提出的方法对图4和图5进行了三维重构，算法分别迭代8步、10步后收敛，非朗伯半球表面的重构结果如图6所示，非朗伯花瓶表面的重构结果如图7所示。为了便于说明重构误差，实验中给出了非朗伯半球、花瓶表面的真实高度值，如图8和图9所示；并给出了本发明重构结果与真实值之间的高度误差图，如图10和图11所示，可以看出，误差绝大部分集中在物体表面的边界之处，其他地方的误差相对较小。

同时，为了定量地分析本发明的性能，从重构时间、重构精度等角度出发，计录了CPU运行时间并计算了高度平均绝对误差和高度均方根误差等指标，具体见表1。

表1 本发明重构结果的性能评价

输入图像	CPU运行时间/秒	高度平均绝对误差/像素	高度均方根误差/像素
				图4	0.25	0.4162	0.5337
图5	0.28	0.3812	0.4835

从图6、7、10、11及表1可以看出，本发明能有效地重构出非朗伯半球、花瓶表面的三维形状，而且CPU运行时间较少，高度平均绝对误差和高度均方根误差均较低，达到了快速、准确的重构需求。

Claims (1)

1.一种非朗伯表面快速三维重构方法，依次包括下述步骤：

一.构建适于非朗伯表面且接近实际成像条件下的成像模型：从三个阶段入手，光源入射过程、非朗伯表面对入射光的反射过程、摄像机对反射光的记录过程，涉及到光源的分布、非朗伯表面的反射特性以及摄像机的投影方式的建模，其中：所述光源为近点光源，而且光源辐照度随着物体与光源之间距离的增大而存在衰减；非朗伯表面的反射特性使用Oren-Nayer反射模型来近似；摄像机的投影方式采用透视投影方式；

二.由上述成像模型和摄像机获取的灰度图像，建立非朗伯表面下SFS问题的图像辐照度方程：该步骤里依次包括下述步骤，

(一)建立三维笛卡尔直角坐标系，设成像平面为x＝(x₁，x₂)位于Z＝-f处，f＞0为摄像机的焦距，摄像机的光轴与Z轴重合，光源位于投影中心，即光心O处，

非朗伯表面S可以表示为Z＝Z(X)，X＝(X₁，X₂)，像平面上一点(x，-f).是物体表面S上点P＝(X₁，X₂，Z)在成像平面的投影点，r为光心O到P点的距离，根据透视投影成像原理有

(1)

令r＝u(x)≥0，于是物体表面S可以用函数表示

(2)

其中：Ω：是定义在实数集上的一个开集，代表图像的区域，

(二)计算光源入射到S(x).点的单位方向向量为

(3)

(三)利用Oren-Nayer反射模型计算非朗伯表面辐射亮度L₃

(4)

其中：θ_i，φ_i和θ_r，φ_r是非朗伯表面某一点的局部坐标系，分别为光源方向向量L和摄像机方向向量V的天顶角和方位角；σ为Gauss分布函数的标准差，代表了表面的粗糙程度；ρ′为漫反射率；I₀为点光源的辐射强度；

α＝max[θ_i，θ_r]；β＝min[θ_i，θ_r]，

由于光源位于投影中心，故θ_i＝θ_r＝α＝β，φ_i＝φ_r，式(4)简化为

(5)

其中：I(x)为摄像机获取的灰度图像；r＝fu(x)，

(四)计算摄像机遵循透视投影方式时，非朗伯表面点S(x)处的法向向量n(x)为

(6)

(五)由于θ_i为n(x)与L(x)之间的夹角，且令v(x)＝lnu(x)，故有

(7)

Q(x)＝f/(||x||²+f²)^1/2. 

其中：

将式(7)代入式(5)，得到非朗伯表面下SFS问题的图像辐照度方程

(8)

其中：I(x)为摄像机获取的灰度图像；

三.利用粘性解理论计算上述图像辐照度方程的粘性解：该步骤依次包括下述步骤

(一)将图像辐照度方程(8)看作是关于的一元二次方程

(9)

(二)求解一元二次方程(9)得

(10)

(三)于是，图像辐照度方程(8)等价为

(11)

(四)利用粘性解理论计算上述图像辐照度方程的解，

显然，上式是一个Hamilton-Jacobi类型的一阶非线性PDE，将简记为p，可以得到相应的Hamilton函数

(12)

(五)应用最优控制理论，式(12)可以转化为如下的控制形式

(13)

其中：f_c(x，a)＝-[2BQ（x)]^-1R^TDRa；l_c(x，a)＝-1/(2B)(1-||a||²)^1/2；是定义在上的单位圆面；D和R满足

(六)图像区域Ω＝(1，m)×(1，n)的均匀离散网格点：x_i，j＝(ih₁，jh₂)，i＝1，2，…，m，j＝1，2，…，n，其中：h＝(h₁，h₂)定义了数值算法中离散网格的大小，

定义对其在时域应用前向Euler公式得到图像辐照度方程求解算法

(14)

(七)使用如下方法逼近H(x_i，j，V_i，j，p)

(15)

其中：f₁(x_i，j，a)、f₂(x_i，j，a)为f₀(x，a)的分量；分别为p的后向、前向差分，

(八)对于逼近方法(15)，计算以下最优问题：

(九)求解算法(14)转化为

(17)

使用Newton法来进行求解，得到图像辐照度方程的粘性解，

(十)使用迭代Fast Marching对上述求解算法(17)进行加速收敛，

四.根据图像辐照度方程解的结果，重构出非朗伯表面的三维形状信息。