0引言随着计算机硬件性能的不断提升以及虚拟植物建模、景观设计在影视和游戏等领域研究的不断深入，植物的形态建模成为可能。叶片作为植物的重要器官之一，具有复杂的生理特性和形态特征，在虚拟场景中很难真实表现(陆声链等，2009)。基于图像的重建方法由于操作简单且恢复效果较好成为热点，而从单幅图像恢复目标区域的3维形状，由于图像所含信息有限，是3维重建中很有挑战的工作之一(孙宇阳，2011)。从图像中获得3维形状的方法大多是基于双目视觉原理，在同一场景中从不同视角获取两幅或多幅图像，利用图像间的视差，得到目标的3维形状。但从单幅图像推导3维形状，可利用的信息十分有限，只能依据图像的阴影、纹理和梯度等几何信息来估计表面深度(赵子阳等，2019)。Horn(1970)首先提出了从阴影估计物体形状，并定义了亮度图与表面形状之间映射关系的反射函数，使根据图像亮度恢复物体的深度成为可能。此后的相关研究多是在该理论模型的基础上添加约束条件使问题更容易求解。在Horn(1970)模型的基础上，Zheng和Chellappa(1991)引入梯度约束，既能保证表面的连续性，又保留了表面可能存在的尖锐特征，抗噪性较强，易于稳定收敛(廖熠和赵荣椿，2001；须明等，2004)。Tsai和Shah(1994)提出线性化反射函数，将非线性问题转换为线性问题，提高了计算效率但易受图像噪声的影响(Zhang等，1994；Zhang等，1999；胡志勇等，2007；Durou等，2008)。钟方洁等人(2012)利用图像中的明暗信息结合纹理信息恢复3维形状。Barron和Malik(2015)提出SIRFS(shape, illumination, and reflectance from shading)方法，利用训练集对反射率先验、形状先验和光照先验进行学习，然后利用单幅图像恢复3维形状。Quéau等人(2017)构造了自然光照下的变分框架，提出基于偏微分方程的变分求解方法。张淑军和刘文晓(2018)通过对目标提取轮廓及骨架，寻找最优“轮廓—骨架”点对，使点对自动生长，网格化得到3维模型。对于叶片的建模，王立臣等人(2010)提出利用Delaunay三角剖分的优化算法将叶片网格化，然后以叶脉作为变形控制骨架，模拟叶片的弯曲变形效果。李书钦等人(2016)通过田间实测数据获取叶片形态特征参数，确定NURBS(non-uniform rational B-splines)曲面控制点构建叶片3维表面。汪丽萍等人(2017)选取轮廓点创建样条曲线，将样条曲线的交点作为控制点构造叶片表面。在图像数据量有限，且不引入外部数据的情况下，大多数基于单幅图像恢复3维形状的方法对表面细节的恢复结果较好，但对整体宏观几何形状的恢复效果不突出(苗旺，2020)。为此，本文在从明暗恢复形状(shape from shading, SFS)方法的基础上，通过对图像进行预处理，为3维形状的推算增加先验信息，并将预处理中检测到的叶片骨架作为先验条件，增强目标物体的表面细节。同时，为了改进物体宏观几何形状的恢复效果，根据图像亮度信息的统计分布和叶片中轴(中线附近的主叶脉)的距离场约束推导宏观几何形状，最后合并多尺度3维形状信息生成目标表面。1研究基础由于成像过程将3维物体以2维图像的方式保存，丢失了大量3维信息，所以从单幅图像生成目标3维表面作为成像的逆过程，初始信息有限，问题求解存在不适定性(黄军等，2019)，且因为植物叶片表层并不完全遵循朗伯反射定理，需要从图像中提取SFS恢复3维形状需要的明暗信息和其他先验条件。相机捕获的图像与场景光照、物体材质和表面形状等因素有很大关系，因此，对图像进行本征图像分解(郭银斌，2014；Barron和Malik，2015)操作，得到反射率图和亮度图。反射率图反映物体颜色、材质等与周围环境无关的特性，亮度图反映光照环境与物体结构作用后的明暗信息。为了避免镜面反射造成的高光，对分解得到的亮度图利用双边滤波器消除高光(卢桂荣等，2014)。亮度图反映了光源与物体表面作用后的结果，SFS是在已知光源的条件下构造2维到3维之间的映射。对光源和反射率的估计主要分为全局方法和局部方法两种，区别在于全局方法假设反射率是常数，局部方法假定各点的反射率不同。Zheng和Chellappa(1991)假设表面由球形面片拼合，局部平坦，表面法向量的偏角和倾角互相独立且分布均匀，然后根据图像强度的统计矩，求出表面倾角和偏角的概率密度分布，求解出光源方向，最后利用图像亮度和亮度平方均值求解表面反射率。SFS从2维图像推导3维形状的原理是当入射光线照射到物体表面，表面会吸收或反射部分光，被表面反射的光进入相机镜头，在底片上成像，因此，图像亮度与光源、物体表面特性和表面形状相关，由于亮度图能获取的只有像素点的位置和亮度值，根据该反射原理，SFS的反射函数模型为 1 $\begin{gathered}E(x, y) \approx R(p, q)= \\\rho \frac{1+p \times p_{s}+q \times q_{s}}{\sqrt{p^{2}+q^{2}+1} \times \sqrt{p_{s}^{2}+q_{s}^{2}+1}}\end{gathered}$ 式中，$ \rho $是表面反射率，反映物体本身的特性，只与物体材质相关，$ \left(p_{s}, q_{s}, 1\right) $为光源方向，$ E\left({x, y} \right) $为图像亮度值，$ R\left({p, q} \right) $为反射函数，$ \left({p, q} \right) $为表面梯度，设$ Z\left({x, y} \right) $是物体表面深度，则 2 $p=\frac{\partial Z(x, y)}{\partial x}$ 3 $q=\frac{\partial Z(x, y)}{\partial y}$ 反射函数表示图像亮度与几何形状的关系，但对于每个像素点来说，有2个待求的未知量$ \left({p, q} \right) $，而只有一个已知的亮度值$ E\left({x, y} \right) $，所以要引入其他的图像特性作为约束条件限制解的范围。SFS算法主要分为演化方法、局部方法、最小化方法和线性化方法(廖熠和赵荣椿，2001)。演化方法是从已知形状或者可以求得形状的点出发，向周围扩展，演化出整个表面。局部方法将局部表面的形状假设与反射模型相结合，构造线性偏微分方程组进行求解。最小化方法是将反射模型和其他约束条件构造为能量函数的形式，转化为最优化问题，对该问题求最小值，引入其他约束条件是为了限制解的范围。线性化方法通过线性化反射函数，将非线性问题简化为线性问题，便于求解。演化方法中，图像的奇点易与图像噪声混淆；局部方法的假设对自然表面很难满足(廖熠和赵荣椿，2001)。SFS的线性化方法首先将函数$ f $关于深度$ Z $的线性逼近用泰勒级数展开，然后利用雅可比迭代求解，简化后可得 4 $\begin{gathered}0 \approx f(Z(x, y)) \approx f\left(Z^{n-1}(x, y)\right)+Z(x, y)- \\Z^{n-1}(x, y) \frac{\mathrm{d} f\left(Z^{n-1}(x, y)\right)}{\mathrm{d} Z(x, y)}\end{gathered}$ 式中，第$ n $次迭代的深度求解为 5 $\left.Z^{n}(x, y)\right)=Z^{n-1}(x, y)-\frac{\mathrm{d} f\left(Z^{n-1}(x, y)\right)}{\mathrm{d} Z(x, y)}$ 假设所有点的初始深度值为0，然后迭代式(5)可得到最终的$ Z $。SFS的最小化方法是通过构造包含亮度约束和其他约束的能量函数并使其最小化，从而求解出深度值。为了减少实际图像中噪声对结果的影响，约束方程引入光强误差作为亮度约束，具体为 6 $E_{1}=\iint(I-R)^{2} \mathrm{~d} x \mathrm{d} y$ 式中，$ I $是测量强度，$ R $是估计的反射强度。引入其他约束条件，构造整体误差表达式，具体为 7 $e=E_{1}+\lambda \times E_{2}+\mu \times E_{3}$ 式中，$ E_2 $和$ E_3 $表示可选择的其他约束表达式，$ \lambda $和$ \mu $是对应的加权因子。SFS的线性化方法计算简单，但由于该方法将反射函数线性化，导致最终结果只是对真正解的一种近似，始终存在系统误差，且受图像质量影响较大；最小化方法的鲁棒性相对较高，在理论上可以无限逼近正解(廖熠和赵荣椿，2001)。两种方法对表面细节的恢复较好，但几乎没有宏观几何形状的起伏，对全局重建效果不佳(苗旺，2020)。2单幅叶片图像3维重建方法对图像进行预处理，以便从图像中提取更多有效信息，将3维形状恢复分解为表面细节和整体宏观几何形状恢复两个方面，算法步骤如图 1所示。 图1 单幅叶子图像恢复3维形状的实现步骤 Steps for realizing the restoration of three-dimensional shape from a single leaf imageFig 12.1图像预处理对输入图像进行预处理，提取目标主体，改善图像质量，为推导3维形状提供更多的线索。首先，使用Chan-Vese方法检测图像中目标主体并标记前景和背景，最大程度降低背景对恢复结果的影响，并找到目标主体的边界，增加约束条件，限制解的范围，有助于形状的恢复。其次，对图像进行本征图像分解，得到亮度图和反射率图，然后对亮度图利用双边滤波器消除高光，接下来对亮度图进行多尺度细节提升操作，选取不同尺度的高斯模糊来增强图像的细节特征。最后，将Canny算子边缘检测的结果视为目标主体的骨架。2.22维图像推导3维形状SFS恢复的3维形状有较清晰的表面细节，但是几乎没有明显的宏观几何形状表现。在图像预处理的基础上，利用图像骨架特征作为约束条件提升表面细节的恢复效果，同时，基于叶片在中轴处弯折的特性和图像亮度分布规律对宏观几何形状进行恢复。2.2.1表面细节恢复SFS的最小化方法鲁棒性高，但相较于线性化方法求解速度慢，线性化方法的求解速度快，但易受图像噪声的影响(须明等，2004)，因此，将线性化方法与最小化方法结合，在图像预处理结果的基础上，利用线性化方法获得最小化方法的初始值。在已知光源方向和表面反射率等信息的条件下，先利用线性化方法初步恢复出物体的深度值，以该值作为最小化方法的初始值，迭代求解出最优深度值。本文使用的最小化方法的能量函数包含亮度约束(重建图像与输入图像的总亮度误差)、光滑性约束(保证曲面的连续性)和强度梯度约束(保证曲面本身锐化的性质)。光滑性约束定义为 8 $E_{2}=\iint\left(p_{x}^{2}+p_{y}^{2}+q_{x}^{2}+q_{y}^{2}\right) \mathrm{d} x \mathrm{d} y$ 式中，$ p_{x}, p_{y}, q_{x}, q_{y} $分别为$ p, q $在$ x, y $方向上的偏导。因为图像是离散的像素点，离散化的表面方向$ \left({{Z_x}, {Z_y}} \right) $与连续的表面方向$ \left({p, q} \right) $之间存在一定误差，所以引入的梯度约束表示为 9 $E_{3}=\iint\left(Z_{x}-p\right)^{2}+\left(Z_{y}-q\right)^{2} \mathrm{~d} x \mathrm{d} y$ 叶脉作为植物叶片的重要特征对3维形状的推导有重要意义，分为主脉、侧脉和细脉。主脉是从叶根到叶尖的最宽脉络，也叫中轴，侧脉和细脉呈羽状分布。经边缘检测得出的叶脉骨架附近存在较大的亮度变化，亮度变化意味着几何形状的变化，因此，在最小化方法中加入叶脉骨架距离场偏置增强表面细节的恢复。距离场是各像素距离边缘检测得出的骨架的最小像素距离，单位像素距离指相邻两个像素之间的距离。距离场偏置既要突显出骨架处的特征，又要保证表面连续性。骨架距离场偏置示意图如图 2所示，其中五角星表示检测的骨架位置，正方形表示距离骨架较近的点，圆形表示距离较远的点，骨架及距离骨架较近的点形状变化明显，而距离较远的点形状几乎不变。距离场偏置定义为 10 $E_{\mathrm{DC}}=a \times \mathrm{e}^{\frac{b}{s}}+c$ 图2 骨架距离场偏置示意图 Skeleton distance field offset diagramFig 2式中，$ s $是当前像素到最近骨架的像素距离，$ a $控制骨架处的凸显程度，$ b $控制随着距离的增加，深度变化的趋势，$ c $为常数，本文中$ c=0 $。偏置后的结果如图 2中三角形位置所示。在第1次迭代时加入距离场偏置，凸显表面细节，之后的迭代过程中利用光滑约束和梯度约束来不断优化表面细节，图 3(a)是仅使用骨架距离场偏置的结果，表面连续，骨架特征较清晰。利用上述方法对不同高斯模糊尺度下细节特征增强的结果恢复出多尺度深度值，得到包含不同形状信息的细节层，最后将这些细节层融合得到最终的表面细节$ {Z_{\rm{D}}} $，结果如图 3(b)所示。 图3 表面细节恢复 Surface detail restorationFig 3 ((a)skeleton distance field; (b)final surface details) 2.2.2宏观几何形状恢复在宏观几何形状的推理中，需要考虑多种影响宏观几何形状的因素，实现叶面凹凸和弯曲的特征。本文利用亮度分布和到叶片中轴的距离恢复宏观几何形状。图 4(a)为单幅图像获得的亮度信息统计。利用差分法计算得到相邻亮度值对应的像素点数量差的最大值和最小值，将最大值和最小值对应的亮度值作为亮暗分区的阈值，计算方法为 11 $\phi_{+}=\arg \max \left\{f_{i}-f_{i+1}\right\},\ i \in[0,254]$ 12 $\phi_{-}=\arg \min \left\{f_{i}-f_{i+1}\right\},\ i \in[0,254]$ 图4 亮度统计规律分区结果 Division results of brightness statisticsFig 4 ((a)brightness statistics histogram; (b)bright and dark partition results) 式中，$ f_i $表示亮度$ i $的频数，图像中亮度值小于$ {\phi _-} $的为暗区，大于$ {\phi _+} $的为亮区，分区结果如图 4(b)所示。将亮暗分区不连通区域的质心作为控制点(见图 5(a))，利用三次贝塞尔曲面恢复出宏观几何形状$ {Z_{\rm{b}}} $(如图 5(b))，可以看出叶片中间的暗区有凹陷，叶尖的亮区部分有凸起，图像中的亮暗分区与3维形状的凹凸基本对应。 图5 控制点恢复的宏观几何形状 Macroscopic geometry of control point restorationFig 5 ((a)control points; (b)restored macro geometry) 大多数叶片关于中轴具有一定程度的弯折，且中轴处曲率变化最大。因此，基于叶片各点到中轴的距离，可以推测出叶片整体的弯折情况。采取交互方法在中轴处手动选取操作点，然后用曲线拟合得到叶片中轴。将叶片弯折的方向定义为正面，则距离叶片中轴越近的位置曲率变化较大，深度值越小。叶片中轴的距离场约束定义为 13 $Z_{\mathrm{mDC}}=\sqrt{d}$ 式中，$ d $是当前像素点到中线的像素距离。对叶子中轴距离场约束的宏观几何形态恢复结果如图 6所示，可以看出叶片关于中轴有一定的形态变化，较为符合叶子的真实形态。 图6 中轴距离场约束恢复的宏观几何形状 Macroscopic geometric shapes recovered from the constrained recovery of the central axis distance fieldFig 6 ((a)blade axis; (b)recovery result from view 1;(c)recovery result from view 2) 将计算得到的表面形状$ {Z_{\rm{b}}} $和$ {Z_{\rm{mDC}}} $分别代入反射函数(式(1))求得亮度图$ {\boldsymbol{R}_{\rm{b}}} $和$ {\boldsymbol{R}_{\rm{mDC}}} $，与原图像的亮度值进行相似度对比，并以对比结果作为其对最终宏观几何形状的权重$ {\lambda_{\rm{b}}} $和$ {\lambda_{\rm{mDC}}} $，最终的宏观几何形状计算公式为 14 $Z_{\mathrm{H}}=\lambda_{\mathrm{b}} \times Z_{\mathrm{b}}+\lambda_{\mathrm{mDC}} \times Z_{\mathrm{mDC}}$ 最后将计算所得的微观表面结果作为细节层，添加到宏观几何形状上，得到叶片的3维表面形状恢复结果$ Z $(如图 7所示)，即 15 $Z=Z_{\mathrm{b}}+Z_{\mathrm{H}}$ 图7 最终3维形状 Final three-dimensional shapeFig 7 ((a)view 1;(b)view 2) 3实验结果与讨论为了验证本文方法的有效性和准确性，设计了7组对比实验，展示了本文方法对叶子3维形状及其他物体表面细节的恢复效果。图 8为本文方法与Tsai线性化方法(Tsai和Shah，1994)、Zheng最小化方法(Zheng和Chellappa，1991)、SIRFS方法(Barron和Malik，2015)和变分方法(Quéau等，2017)对单幅植物叶片图像表面细节的恢复结果对比，从上至下依次为叶片1、叶片2、叶片3、叶片4和叶片5，图中展示了两个视角的恢复结果。图 9为本文方法与Tsai线性化方法和Zheng最小化方法对单幅硬币和恐龙图像表面细节的恢复结果对比。图 8和图 9对本文方法有优势的地方用红圈标示。Tsai线性化方法的恢复结果如图 8(b)所示，纹理比较清晰，但该方法易受图像噪声的影响，鲁棒性低。如图 9(b)中硬币的恢复结果，在数字“1”附近由于噪声的影响恢复效果较差。Zheng(Zheng和Chellappa, 1991)的结果如图 8(c)所示，叶片表面纹理不明显，相比于Tsai线性化方法的恢复结果具有一定差距，但鲁棒性较高，如图 9(c)硬币的恢复结果，数字“1”附近的噪声对恢复结果影响不大。上述方法对叶片表面细节的恢复有较好的表现，但没有考虑叶片的宏观几何形态(叶片的凹凸和弯折)。本文方法基于线性化方法提供初始深度值，利用最小化方法优化叶片表面细节，使用骨架距离场偏置增强最小化方法恢复表面细节，并基于亮度统计找出分布于不同亮度区域的控制点对叶片表面进行三次贝塞尔曲面拟合，利用叶片中轴的距离场条件约束表面形状，恢复叶片的宏观几何形态。本文方法对5种不同叶片的表面细节恢复结果和整体表面形态如图 8(g)所示，相较于前两种方法，本文方法的恢复结果表面细节更清晰，并具有更明显的宏观几何形状。SIRFS方法利用训练集得到关于反射率、形状和光照的先验，根据图像的先验信息约束来推断单个图像最有可能的3维形状。变分方法构造了自然光照下SFS通用的数值框架，增加了形状先验和非线性平滑项。通过观察实验结果，注意到SIRFS方法表面细节不太清晰，但有明显的几何形状变化(图 8(d))，而变分法恢复结果有一部分细节明显，但无明显的宏观形状变化(图 8(e))。使用本文的控制点恢复结果(图 8(f))与SIRFS和变分法结果进行比较，本文方法的结果更符合真实的叶片凹凸特性。 图8 不同方法对单幅植物叶片图像的表面细节恢复结果 Recovery detail restoration results of single plant leaf images by different methodsFig 8 ((a)original images; (b)Tsai linear approach(Tsai and Shah, 1994); (c)Zheng minimization approach(Zheng and Chellappa, 1991); (d)SIRFS(Barron and Malilk, 2015); (e)variational approach(Quéau et al., 2017); (f)our control point method; (g)ours) 图9 不同方法对单幅硬币和恐龙图像的表面细节恢复结果 Surface detail restoration results of single coin and dinosaur images by different methodsFig 9 ((a)original images; (b)Tsai linear approach(Tsai and Shah, 1994); (c)Zheng minimization approach(Zheng and Chellappa, 1991); (d)ours from view 1;(e)ours from view 2) 在仅有单幅图像，不引入外部数据的情况下，SFS算法只能得出相对的3维形状，所以其重建效果与立体视觉重建效果或其他数据进行定量对比是无意义的(苗旺，2020)。为了验证本文方法的有效性，用误差—信息熵之比来衡量恢复效果。误差是将所求3维形状根据成像原理利用反射函数计算得到的亮度结果与输入图像的亮度信息进行误差分析；信息熵用来描述信息的丰富程度，熵值越大表明恢复的深度之间差异越大，表面几何形态变化越大。误差越小，信息熵越大时，误差—信息熵之比越小，恢复效果越好。误差—信息熵之比计算为 16 $e r r=\frac{\sum\limits_{i=1, j=1}^{M, N}\left|R_{i, j}-I_{i, j}\right|}{M \times N \sum\limits_{Z=0}^{Z_{\max }}\left|P(Z) \log _{2} P(Z)\right|}$ 式中，$ R_{i, j} $表示像素$ \left({i, j} \right) $处估计的亮度值，$ I_{i, j} $为测量亮度值，$ M $和$ N $表示图像的像素大小$ P\left(Z \right) = {n_Z}/(M \times N), {n_Z} $表示对应深度值的数量。不同方法对叶子3维重建的误差—信息熵之比如表 1所示。表1 不同方法对叶子3维重建的误差—信息熵之比 方法 叶片1 叶片2 叶片3 叶片4 叶片5 Tsai线性化(Tsai和Shah, 1994) 0.282 1 0.143 7 0.217 0 0.111 6 0.249 0 Zheng最小化(Zheng和Chellappa, 1991) 0.226 8 0.262 4 0.383 5 0.160 3 0.356 2 SIRFS(Barron和Malik, 2015) 0.107 1 0.141 0 0.124 2 0.152 4 0.283 5 变分法(Quéau等, 2017) 0.285 4 0.136 0 0.117 8 0.064 9 0.245 9 本文控制点 0.055 5 0.178 0 0.064 5 0.015 0 0.175 5 本文 0.053 2 0.085 0 0.092 7 0.082 6 0.191 4 Error-Information entropy ratio of 3D reconstruction of leaves by different methodsTable 1 加粗字体表示各列最优结果。虽然本文方法是针对单幅叶片图像的3维形状恢复，但也适用于其他物体基于单幅图像的表面细节恢复。本文方法与Tsai线性化方法(Tsai和Shah, 1994)和Zheng最小化方法(Zheng和Chellappa, 1991)对硬币和恐龙的单幅图像进行3维重建的误差—信息熵之比如表 2所示。可以看出，对其他物体表面细节的恢复，本文方法相较对比方法仍具有一定优势。表2 不同方法对其他物体3维重建的误差—信息熵之比 方法 硬币 恐龙 Tsai线性化(Tsai和Shah, 1994) 0.028 9 0.146 6 Zheng最小化(Zheng和Chellappa, 1991) 0.330 1 0.335 2 本文 0.023 5 0.143 1 Error-Information entropy ratio of 3D reconstruction of other objects by different methodsTable 2 加粗字体表示各列最优结果。4结论针对单幅图像的3维形状恢复，SFS方法对表面细节恢复有比较好的表现，但是未考虑物体的宏观几何形状变化。本文在图像数据量有限，且不引入外部数据的情况下，针对单幅叶片图像进行3维形状恢复，利用SFS线性化方法的深度值结果作为最小化方法的初始深度值，同时为最小化方法引入距离场偏置，既保证了表面的连续性，又使得恢复的3维表面细节更加突出。此外，根据图像明暗信息的统计分布，将图像分割为不同亮度区域，并在不同亮度区域确定控制点，在此基础上，利用三次贝塞尔曲面对叶片表面进行拟合，恢复叶片的凹凸特征。为了恢复叶片整体的弯折特征，提出了关于叶片中轴的距离场约束，两者结合对叶片的整体形状具有较好的恢复效果。针对算法的可靠性设计了5组对比实验。结果表明，本文方法恢复的叶片表面细节清晰可靠，并且有较好的宏观几何形状。另外设计了2组实验表明本文方法同样适用于其他物体的表面细节恢复。但基于中轴的距离场约束是针对叶片的特性提出的，并不具备通用性。此外，若图像中存在阴影遮挡，则遮挡部分的细节难以恢复。未来将进一步优化本文方法，提高计算效率，同时基于本文方法的恢复结果，结合不同的环境模型，模拟叶片在不同条件下的形态变化，使其有更广泛的应用。