CN101739714A - 一种从图像中重建枝干模型的方法 - Google Patents

Abstract

本发明涉及计算机视觉领域，具体涉及一种从图像中重建植物枝干模型的方法，所述方法包括步骤：获取若干幅图像以及所述图像间的相对关系；获取所述图像的alpha matte作为后续操作的依赖图像；手工指定其中两幅图像中枝干底部的位置和其中一幅图像中对应的决定半径的像素；依据上述图像中指定的枝干初始位置计算三维初始位置以及半径；以所得三维初始位置以及半径为基础逐步生长三维枝干模型。本发明方法简单且最大限度地利用了图像上的信息，能够生成与图像较为符合的枝干模型。

Description

一种从图像中重建枝干模型的方法

技术领域

本发明涉及计算机视觉和图形学领域，具体涉及一种从图像中重建枝干模型的方法。

背景技术

植物的真实感建模在场景绘制以及珍稀植物原生态保存等方面都具有极其重要的意义。

近十年内，基于真实测量信息的植物重建方法得到广泛关注，其根本原因在于这类方法能够以较强的真实感重现自然界中植物。然而由于植物类型繁多，结构复杂，目前的重建方法多局限于特定形态植物。其中以针对室外较大树木的枝干建模方法居多。Shlyakher.等人(Shlyakher I.，Rozenoer M.，Dorsey J.，and Teller S.Reconstructing 3Dtree models from instrumented photographs，IEEE Computer Graphicsand Applications，2001)用图像恢复的visual hull及其骨架来控制L-system相关参数，生长得到较符合实际的树木模型。Han等人(HanF.，and Zhu S.C.Bayesian reconstruction of 3D shapes and scenes from asingle image，Proceedings of IEEE Workshop on Higher-LevelKnowledge in 3D modeling and Motion Analysis，2003)利用Bayesian方法结合简单的空间分布先验和图像的硬性观测，能够从纯枝干的单幅图像中推测被遮挡部分并恢复完整的枝干模型。Teng等人(Teng C.H.，Chen Y.S.，and Hsu W.H.Constructing a 3D trunk model from twoimages，Graphical Models，2007)利用立体视觉方法重建主枝干。Neubert等人(Neubert B.，Franken T.，and Deussen O.Approximateimage-based tree modeling using particle flows，ACM Transactions onGraphics，2007)利用图像上方向场的约束，用粒子流模拟生成枝干模型。Tan等人(Tan P.，Zeng G.，Wang J.D.，Wang S.B.，and Quan L.Image-based tree modeling，ACM Transactions on Graphics，2007)采用Structure from motion(SFM)技术，从多幅图像中得到三维点云模型，以此为基础重建可见枝干，并合成不可见枝干。除了图像数据之外，激光数据也是重建枝干的重要凭据。例如，Xu等人(Xu H，GossettN，and Chen B.Knowledge and heuristic-based modeling oflaser-scanned trees.ACM Transactions on Graphics，2007)从激光点云数据中抽取主枝干，并依据主枝干合成细小的难以从激光数据中可靠抽取的较细枝干。然而，以上方法在方便性以及真实度上都欠缺，尤其是重建被叶片遮挡的枝干部分。

发明内容

本发明的目的是提供一种能够保证真实度的植物重建方法，获得与图像符合的三维枝干模型。

为了达到上述发明目的，本发明提供了一种从图像中重建植物枝干模型的方法，所述方法包括以下步骤：

S1：获取若干幅图像以及所述图像间的相对关系；

S2：获取所述图像的alpha matte作为后续操作的依赖图像；

S3：手工指定其中两幅图像中枝干底部的位置和其中一幅图像中对应的决定半径的像素；

S4：依据上述图像中指定的枝干初始位置计算三维初始位置以及半径；

S5：以所得三维初始位置以及半径为基础逐步生长三维枝干模型。

其中，所述步骤S5包括：

S5-1：定义当前生长阶段的节点状态

其中

是第i个枝干当前生长阶段节点的状态，P_it表示什么，D_it表示什么，r_it表示什么；

S5-2：利用下式求解当前生长阶段的所有枝干节点的状态，

P(k_t，X_t|Z_t)∝P(k_t，X_t|k_t-1，X_t-1)P(Z_t|k_t，X_kt)

其中，所述步骤S5-2估计当前生长阶段的所有枝干节点状态的方法包括：

S5-2-1：初始化当前生长t状态的状态(k_t，X_kt)；

S5-2-2：按照proposal函数提议一个新的状态(k′_t，X′_kt)；

S5-2-3：利用下式计算新状态的接受率

$\mathrm{\α}=\frac{P\left(X_t^{\′}\right|Z_t)Q(X_t^s;X_t^{\′})}{P\left(X_t^s\right|Z_t)Q(X_t^{\′};X_t^s)}$

S5-2-4：当接受率大于1时，直接接受；否则，以概率a接受，其中a∈[0，1]；

S5-2-5：按照上述步骤迭代，迭代次数s＝1...N-1，生成若干采样；

S5-2-6：求取采样的期望作为t阶段的状态。

其中，在步骤5-2求解每个生长时段状态之前，检测该时段增加或者减少的枝。

其中，在步骤5-2求解每个生长时段状态之前，检测该时段生长或者停长的枝。

其中，在步骤5-2求解每个生长时段状态之前，对于保持生长的枝进行更新。

与现有技术相比，本发明的技术方案构建树木图像的柔性似然函数，第一次结合图像的柔性似然函数与植物模型的先验知识重建带叶片植物的枝干模型，方法简单且最大限度地利用了图像上的信息，能够生成与图像较为符合的枝干模型。

附图说明

图1为本发明的从图像中重建枝干模型的方法流程图。

具体实施方式

以下实施例用于说明本发明，但不用来限制本发明的范围。

本发明的技术方案包括如下步骤：获取若干幅图像以及所述图像间的相对关系；获取所述图像的alpha matte作为后续操作的依赖图像；指定其中两幅图像中枝干底部的位置和其中一幅图像中对应的决定半径的像素；依据上述图像中指定的枝干初始位置计算三维初始位置以及半径；以所得三维初始位置以及半径为基础逐步生长三维枝干模型。

具体叙述如下：输入两幅图像，依据ImageModeller软件得到两幅图像的相对关系，即确定每幅图像对应的相机光心位置、焦距和主点，然后对图像做alpha matting得到alpha值作为象素值的图像。标定和alpha matting采用现有技术。通常，一幅图像由前景和背景组成，图像中不同象素处两者参与的比例不一样。比如，在只有背景的地方，前景所占比例为0％，在只有前景的地方，前景所占比例为100％，这个比例即alpha值。alpha值作为象素值的图像即alpha matte。而求取alpha matte的过程，称之为alpha matting。

Alpha matting技术采用CVPR06上由Levin等人发表的″AClosed Form Solution to Natural Image Matting″。然后，在两幅图像上给定初始的根部生长点，并在其中一幅上指定决定半径的位于根部生长点一侧的一个象素点。为方便选取不同图像上的对应生长点，可以在拍摄植物之前在其底部系上标志。根据初始标注的信息找到的三维枝干信息。接下来的基于贝叶斯方法从根部逐段生长枝干的过程是本发明的核心。

本发明的核心处理过程如下。定义当前生长阶段的节点状态

其中

是第i个枝干当前生长阶段节点的状态，包括该枝干在当前阶段的空间坐标，生长方向以及半径。本发明估计当前生长阶段的所有枝干节点最优状态的方法是最大化贝叶斯后验概率

P(X_t|Z_t)∝P(Z_t|X_t)P(X_t|X_t-1)。    (1)

在实际的生长过程中枝干数目是可变的，因此，加入指示集k_t，用来标识当前阶段枝干的编号集合。比如，如果当前枝干分别为：1、3号枝干，则k_t＝{1，3}。按照生长过程中出现的先后顺序编号，初始指定的枝干编号为1。所举例子中，指示集中包括1、3号枝干说明：生长过程中总共加入3个枝干，而当前生长阶段2号枝干已经停止生长，只有1、3号枝干活跃。在可变枝干情况下，我们需要最大化的贝叶斯后验概率为：

P(k_t，X_t|Z_t)∝P(k_t，X_t|k_t-1，X_t-1)P(Z_t|k_t，X_kt)。(2)

以下，对公式2中的逐项给定本发明的定义方法。

1.公式2中第一项为取决于前一生长阶段的先验模型，可以分解并表示为：

P(k_t，X_t|k_t-1，X_t-1)＝P(X_t|k_t，k_t-1，X_t-1)P(k_t|k_t-1，X_t-1)。(3)

公式3右边第一项P(X_t|k_t，k_t-1，X_t-1)又可以进一步分解为新加入的枝干和原有的枝干，即：

P(X_t|k_t，k_t-1，X_t-1)＝P(X_ke)P(X_ks|X_ks(t-1))。(4)

这里之所以能够分解是假定新加入的枝干和原有枝干是独立互不影响的。其中，P(X_ke)表示新加枝干的位置先验分布，P(X_ks|X_ks(t-1))则是原有枝干的先验生长概率模型。定义为：

$P\left(X_{\mathrm{ks}}\right|X_{\mathrm{ks}(t-1)})\∝\underset{i}{\mathrm{\Π}}P\left(X_{\mathrm{it}}\right|X_{i(t-1)})\underset{i,j\∈E}{\mathrm{\Π}}\mathrm{\Ψ}(X_{\mathrm{it}},X_{\mathrm{jt}}).---\left(5\right)$

公式5右侧第一项为第i个枝干的生长先验模型，定义如下：

$[{\stackrel{\→}{P}}_{\mathrm{it}},{\stackrel{\→}{D}}_{\mathrm{it}},r_{\mathrm{it}}]~[N({\stackrel{\→}{P}}_{i(t-1)}+{\stackrel{\→}{D}}_{\mathrm{it}}{r}_{\mathrm{it}},{\mathrm{\σ}}_p^2),{\stackrel{\→}{D}}_{i(t-1)},N(r_{i(t-1)},{\mathrm{\σ}}_r^2)].---\left(6\right)$

公式(5)右侧第二项Ψ(X_it，X_jt)＝exp(-g(X_it，X_jt))是X_it与X_jt之间的Gibbs分布形式的势函数，g(X_it，X_jt)是一个惩罚函数，在本发明中定义如下。首先，令g′＝(r_it+r_jt-||P_it-P_jt||)，表示X_it与X_jt之间重叠部分大小，那么

1)当g′≥0时，g(X_it，X_jt)＝N(r_it+r_jt-||P_it-P_jt||)/(r_it+r_jt)，其中，N是一个较大的数值，如100，使得当两个枝干完全重合时，即g(X_it，X_jt)＝N时，exp(-N)是一个接近0的概率值。

2)当g′＜0时，直接令g(X_it，X_jt)＝0，使得距离没有冲突的枝干都可以得到exp(0)＝1，表示其生长或被加入均不影响其他现有枝干。

公式(3)的第二项P(k_t|k_t-1，X_t-1)指定哪个枝干可能进入，离去或者继续保留的概率值。

公式(2)的第二项为似然概率模型，其定义方式为

$P\left(Z_t\right|k_t,X_{\mathrm{kt}})=\underset{i}{\mathrm{\Π}}C({\stackrel{\→}{P}}_{\mathrm{itr}},{\stackrel{\←}{P}}_{\mathrm{itl}})G\left({\stackrel{\→}{P}}_{\mathrm{itr}}\right)G\left({\stackrel{\→}{P}}_{\mathrm{itl}}\right),---\left(7\right)$

其中，

C(a，b)函数是描述a到b的连通性，R是旋转函数，G是梯度函数，返回alpha图像上对应位置处的梯度值。

本发明采用RJMCMC方法求解MAP(Maximum a posteriori)分布作为当前生长阶段的状态。问题形式类似于Khan等人的文章MCMC-based particle filtering for tracking a variable number ofinteracting targets。该文在处理多目标跟踪时同样采用RJMCMC方法。

MCMC即Markov Chain Monte Carlo，其工作方式为：在状态空间X中定义马尔可夫链，相对于其他求解概率分布函数的期望值的算法，如拒绝采样，以及重要性采样等，MCMC可以更有效地处理高维度问题。在我们的应用当中，枝干的数目是比较多的，因此，采用MCMC的方法较为合适。在本发明中，MCMC所定义的马尔可夫链的稳态分布即是公式(2)后验分布。我们的目标是求解该稳态分布的期望值。Metropolis-Hastings(MH)是较常用的一种模拟并求解马尔可夫链稳态分布期望值的方法。对于每个生长阶段，采用MH方法生成一系列采样，作为估计该阶段最优解的依据。

Metropolis-Hastings算法(可参考Hastings的Monte carlo samplingmethods using Markov chains and their applications以及Bishop撰写的Pattern Recognition and Machine Learning一书)如下：任给一个初始状态X_t ¹，然后一直迭代，迭代次数s＝1...N-1：

根据proposal函数：Q(X′_t；X_t ^s)提议一个新的状态X′_t。

计算新状态的接收率：

P(X′_t|Z_t)是新状态的后验概率值，P(X_t ^s|Z_t)是原状态的概率值，Q(X_t ^s；X′_t)是依据proposal函数从新状态回到原状态的概率值，Q(X′_t；X_t ^s)是依据proposal函数从原状态到新状态的概率值。

如果接受率大于1，则直接接受。否则，以概率a接受，a是[0，1]之间的随机数，也就是从[0，1]随机采样一个数，如果大于接受率α，则接受，否则，拒绝。

具体处理过程中，考虑到一开始的MH算法可能不满足稳态分布，舍弃MH初始阶段的B个采样点。另外，为了去掉MH采样时邻近采样点的相关性，采取间断取样的方法。假定采样间隔为M，需要的样本个数为N，则MH算法的总采样次数为B+MN。

从MH的算法步骤也可以看出Proposal函数对于MH算法性能起着决定性作用。Proposal函数应当满足易估计性(计算该函数值)和有效采样性(从该函数中采样)。实际上是把对复杂的稳态函数的采样转为对proposal函数的采样。对于连续的状态空间，常用的proposal函数即以当前状态为中心的高斯函数。而高斯函数的方差选择对于MH算法的速度影响甚大。方差小了，则新状态的接受率较高，但是则会导致MH算法以非常小的随机游走速度遍历状态空间，速度方面会收到影响。但是，如果方差较大，则拒绝率较高。因为我们的稳态分布函数较为复杂，方差大的proposal函数所建议的新状态很大程度上会落在稳态函数的概率值小的地方，导致直接被拒绝。因此，方差过大或过小都会导致MH算法的收敛速度降低。本发明proposal方式为，按照一致分布随机选取k_t集合中的枝干，对单独的枝干做单独的proposal，proposal函数的方差选择为该枝干半径的十分之一。

RJMCMC即Reversible-Jump Markov Chain Monte Carlo，是可变维度状态空间的MCMC。Jump意味着可以改变维度，Reversible意思指每个jump都对应的这种移动方式(可参考文献P.Green.Reversible jump Markov chain Monte Carlo computation and Bayesianmodel determination.Biometrika，82：711-732，1995.以及P.J.Green.Trans-dimensional Markov chain Monte Carlo.Chapter in HighlyStructured Stochastic Systems，2003)。在用RJMCMC求解过程中，任意选取一个初始状态，(k_t，X_kt)∈χ。然后随机选取一种移动方式m，包括增加状态维度(新加入枝干)，降低状态维度(部分枝干停止生长)，或者保持维度不变。在此过程中，新状态(k′_t，X′_kt)的提出采用的proposal函数为Q_m(k′_t，X′_t；k_t，X_t)。假定逆向proposal函数为Q_m′(k_t，X_t；k′_t，X′_t)(m′是对应的m的逆向运动)，则新状态的接受率为：

$\mathrm{\α}=\frac{P(k_t^{\′},X_t^{\′}|Z_t)p_{m^{\′}}{Q}_{m^{\′}}(k_t,X_t;k_t^{\′},X_t^{\′})}{P(k_t,X_t|Z_t)p_m{Q}_m(k_t^{\′},X_t^{\′};k_t,X_t)}.---\left(8\right)$

为了保证可逆性，相对应于增加新枝干的proposal应当有一个去除枝干的proposal。考虑到模拟枝干生长过程中，有新枝干加入，也有原有枝干停止生长。因此，本发明采用了两对jump的可逆操作，分别为：增、减和长、停，再加上对于保持继续生长的枝干有一个同维度空间内的更新操作，共5种操作。本发明分别将这5种操作的选择概率表示为P_A，P_D，P_G，P_S，P_U。那么决定这五个概率值时应当满足P_A+P_D+P_G+P_S+P_U＝1。

1)在求解每个生长时段状态前，检测在该时段的所有可能增加的枝。在这些枝基础上做增减的proposal(即可提高RJMCMC效率的数据驱动方式：Data Driven RJMCMC，可参考文献Z.W.Tu and S.C.Zhu.Image segmentation by data-driven Markov chain Monte Carlo)。

增操作如下：如果检测目标a不在标识集合k_t中，则建议增加它。具体方法就是以一致概率1/|k_d\k_t|随机采样一个候选目标a，k_d\k_t是检测到的但是没有在k_t集合中的。将a加入当前k_t集合中。如果k_d\k_t集合为空，则令P_A＝0。其接受率为：

$a_A=P\left(Z_t\right|X_a)=\frac{P(k_t^{\′},X_{\mathrm{kt}}^{\′}|Z_{t-1})}{P(k_t,X_{\mathrm{kt}}|Z_{t-1})}\frac{p_D}{P_A}\frac{|k_d\backslash k_t|}{|k_t^{\′}\∩k_d|}.---\left(9\right)$

其中，第一项为新加入枝干的似然概率，第二项分数形式是新状态的生长先验分布比上原状态的生长先验分布，定义为：

$P(k_t,X_{\mathrm{kt}}|Z_{t-1})=P\left(X_{\mathrm{ka}}\right)P\left(X_{\mathrm{kut}}\right|X_{\mathrm{ku}(t-1)})P\left(k_t\right|k_{t-1},X_{k(t-1)})\underset{i,j\∈E}{\mathrm{\Π}}\mathrm{\Ψ}(X_{\mathrm{it}},X_{\mathrm{jt}})$

上式第一项建模了哪个位置出现新枝干的可能性比较大，我们设置其为一致分布，此项不起作用。第二项是上一生长阶段就存在的枝干的生长先验分布，其定义方式见公式(6)。第三项建模了哪个枝干加入的可能性，也是一致分布，不起作用。第四项是公式(5)定义的枝干间的相互影响。P_A与P_D通常是相等的。第四项分数是为了平衡增减操作对应的目标个数。最终，如果接受率大于1，则直接接受，否则以概率a接受。如果被拒绝，则恢复原状态。

减操作目标是使得状态可逆，即加入的枝干有可能恢复原来的状态。其过程如下：以一致概率1/|k_t∩k_d|随机采样一个候选目标d，k_t∩k_d是即在检测集合又在k_t集合中的目标集。将d从当前k_t集合中删除。如果k_t∩k_d集合为空，则令P_D＝0。其接受率为：

$a_D=\frac{1}{P\left(Z_t\right|X_d)}\frac{P(k_t^{\′},X_{\mathrm{kt}}^{\′}|Z_{t-1})}{P(k_t,X_{\mathrm{kt}}|Z_{t-1})}\frac{p_D}{P_A}\frac{|k_t\∩k_d|}{|k_d\backslash k_t^{\′}|}.---\left(10\right)$

2)停操作是枝干停止状态的模拟。长操作是其可逆过程，使得被停止的枝干有可能恢复原有状态。

长操作如下：从集合k_t-1\k_t中一致采样即以概率1/|k_t-1\k_t|随机采样一个目标g，其中k_t-1\k_t表示在前一个生长阶段集合但不在当前生长阶段集合中的目标枝干。在增加该目标时，依据其在前一生长阶段的状态值决定的当前状态分布，即该目标枝的先验生长模型，P(X_gt|X_g(t-1))，采样得到其当前状态。如果k_t-1\k_t集合为空，则令P_G＝0。其接受率为：

停操作如下：以一致概率1/k_t\k_d|采样候选目标枝s，其中k_t\k_d表示在当前标识集合k_t但不在检测集合k_d中的目标集合。改变其状态，即从当前k_t集合中去除。如果k_t\k_d集合为空，则令P_S＝0。计算接受率：

$a_S=\frac{1}{P\left(Z_t\right|X_s)}\frac{P(k_t^{\′},X_{\mathrm{kt}}^{\′}|Z_{t-1})}{P(k_t,X_{\mathrm{kt}}|Z_{t-1})}\frac{p_G}{P_S}\frac{|k_t\backslash k_d|P\left(X_{\mathrm{st}}\right|X_{s(t-1)})}{|k_{t-1}\backslash k_t|}.---\left(12\right)$

3)同维度空间的更新操作如下：以一致概率，1/k_t，从k_t集合中选取一个目标枝u，其状态为X_ut，采样proposal函数Q(X′_ut，；X_ut)提议一个新的状态X′_ut。如果，k_t集合为空则令P_U＝0。Q(X′_ut；X_ut)定义如前述，是以X_ut为中心的，以u的半径十分之一为方差的高斯函数。计算接受率：

${\mathrm{\α}}_U=\frac{P\left(X_{\mathrm{ut}}^{\′}\right|Z_t)Q(X_{\mathrm{ut}};X_{\mathrm{ut}}^{\′})}{P\left(X_{\mathrm{ut}}\right|Z_t)Q(X_{\mathrm{ut}}^{\′};X_{\mathrm{ut}})}=\frac{P\left(Z_t\right|X_{\mathrm{ut}}^{\′})P\left(X_{\mathrm{ut}}^{\′}\right|Z_{t-1})Q(X_{\mathrm{ut}};X_{\mathrm{ut}}^{\′})}{P\left(Z_t\right|X_{\mathrm{ut}})P\left(X_{\mathrm{ut}}\right|Z_{t-1})Q(X_{\mathrm{ut}}^{\′};X_{\mathrm{ut}})}.---\left(13\right)$

其中，P(X′_ut|Z^t-1)是公式5定义的考虑到枝干间相互影响的预测先验。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由其权利要求限定。

Claims (8)

1.一种从图像中重建枝干模型的方法，其特征在于，所述方法包括以下步骤：

S1：获取若干幅图像以及所述图像间的相对关系；

S2：获取所述图像的alpha matte作为后续操作的依赖图像；

S3：指定其中两幅图像中枝干底部的位置和其中一幅图像中对应的决定半径的像素；

S4：依据上述图像中指定的枝干初始位置计算三维初始位置以及半径；

S5：以所得三维初始位置以及半径为基础建立三维枝干模型。

2.如权利要求1所述的从图像中重建枝干模型的方法，其特征在于，所述步骤S5包括：

S5-1：定义当前生长阶段的节点状态其中是第i个枝干当前生长阶段节点的状态，

即第i个枝干在t生长阶段的节点三维空间位置，是节点的生长方向，r_it为节点对应半径；

S5-2：利用下式求解当前生长阶段的所有枝干节点的状态，

P(k_t，X_t|Z_t)∝P(k_t，X_t|k_t-1，X_t-1)P(Z_t|k_t，X_kt)

其中，Z_t是图像上的观测，k_t用于指示当前活跃的枝干。

3.如权利要求2所述的从图像中重建枝干模型的方法，其特征在于，所述步骤S5-2中求解当前生长阶段的所有枝干节点状态的方法包括：

S5-2-1：初始化当前生长t阶段的状态(k_t，X_kt)；

S5-2-2：按照proposal函数提议一个新的状态(k′_t，X′_kt)；

S5-2-3：计算新状态的接受率；

S5-2-4：当接受率大于1时，接受该状态为采样；否则，以概率a接受，其中a∈[0，1]；

S5-2-5：按照上述步骤迭代，迭代次数s＝1...N-1，生成若干采样；

S5-2-6：求取采样的期望作为t阶段的状态。

4.如权利要求2所述的从图像中重建枝干模型的方法，其特征在于，在步骤5-2求解每个生长时段状态之前，检测该时段增加的枝。

5.如权利要求2所述的从图像中重建枝干模型的方法，其特征在于，在步骤5-2求解每个生长时段状态之前，检测该时段减少的枝。

6.如权利要求2所述的从图像中重建枝干模型的方法，其特征在于，在步骤5-2求解每个生长时段状态之前，检测该时段处于生长状态的枝。

7.如权利要求2所述的从图像中重建枝干模型的方法，其特征在于，在步骤5-2求解每个生长时段状态之前，检测该时段处于停长状态的枝。

8.如权利要求2所述的从图像中重建枝干模型的方法，其特征在于，在步骤5-2求解每个生长时段状态之前，对于保持生长的枝进行更新。

Images