Embedded Deformation for shape manipulation
发表期刊:siggraph 2007
1. 概述
本论文通过建立“deformation graph”的变形模型,并在计算变形模型在操作过程中仿射变换的最优值,来实现三维模型的直接操控变形,支持细节保留、粒子系统变形及模型局部变化。
图 1. 操控恐龙闭嘴
图 1 上面的网格图称为deformation graph,下面是操控恐龙闭嘴的过程,其中模型绿色区域是固定点区域,黄色代表可以操作的点,移动黄色的点使恐龙闭嘴。
2. 算法原理
2.1 Deformation Graph
获得deformation graph的过程,就是通过在原来的三维模型上均匀采样顶点,然后根据设定的采样密度来删除一定数量的顶点,直到到达密度要求,并使这些顶点构成网格。Nodes数量在200-300个之间较为合适。
Deformation Graph是一种新型的可变形模型,可作用于所有欧式三维空间
R
3
\boldsymbol{R}^3
R3的点node,node的位置可表示为
g
j
∈
R
3
,
j
∈
1
⋯
m
g_j \in R^3, j \in 1\cdots m
gj∈R3,j∈1⋯m ,当nodes位置确定后,通过在影响同一vertex的两个nodes之间建立edge来确定连通性。Node
j
j
j的邻点
N
(
j
)
N(j)
N(j)包含了所有与
j
j
j有edge相连的nodes。
空间中的变形由一系列仿射变换所定义,仿射变换包含了旋转项和位移项,定义为
3
×
3
3 \times 3
3×3的旋转矩阵
R
j
R_{j}
Rj和
3
×
1
3 \times 1
3×1的位移向量
t
j
t_{j}
tj。模型上每个点
p
p
p的位置都受到与这个nodes有边相连的顶点的影响,变形后的位置
p
~
\tilde{p}
p~由公式 1 得到:
p
~
=
R
j
(
p
−
g
j
)
+
g
j
+
t
j
⃗
(1)
\tilde{p} = R_{j}(p-g_{j})+g_{j}+\vec{t_{j}} \tag1
p~=Rj(p−gj)+gj+tj(1)
Deformation graph也可以对用vertex定义的模型进行变形,用每个vertex位置加权和来计算变形后新的顶点位置
v
~
i
\tilde{v}_i
v~i:
v
~
j
=
∑
j
=
1
m
w
j
(
v
j
)
[
R
j
(
v
i
−
g
j
)
+
g
j
+
t
j
]
(2)
\tilde v_j=\sum^m_{j=1}w_j(v_j)[R_j(v_i-g_j)+g_j+t_j] \tag2
v~j=j=1∑mwj(vj)[Rj(vi−gj)+gj+tj](2)
w
j
(
v
i
)
=
(
1
−
∥
v
i
−
g
j
∥
/
d
m
a
x
)
2
(3)
w_j(v_i)= (1-\Vert v_i-g_j\Vert / d_{max})^2 \tag3
wj(vi)=(1−∥vi−gj∥/dmax)2(3)
权重 w j w_j wj根据vertex空间位置而改变,用于限制变形图中vertex对最近 k 个nodes 的影响,越远的点施加的影响就越小。 d m a x d_{max} dmax表示 v i v_i vi到第 k+1 近的点的距离,这样可以控制只受最近 k 个点影响。例如, v i v_i vi的邻点与 v i v_i vi的距离分别为 1, 3, 5, 7, 9, 11, 13…。当我们取 k = 4时, d m a x d_{max} dmax = 9,此时我们只计算前四个点对 v i v_i vi的影响,过远的点我们就不考虑了。
2.2 Optimizer
确定deformation graph后,即可通过选择vertex,并移动它们来操纵变形。通过最小化能量方程来求解每个顶点的旋转矩阵和位移向量,能量方程包含旋转、正则化和约束三个误差项。
2.2.1 Rotation
为了保证变形后的细节,仿射变换的
R
R
R,应该是旋转矩阵,也就是满足旋转群
S
O
(
3
)
SO(3)
SO(3)条件:矩阵的三个列必须都是单位长度,且两两正交,公式表示为
Rot
(
R
)
=
(
c
1
⋅
c
2
)
2
+
(
c
1
⋅
c
3
)
2
+
(
c
2
⋅
c
3
)
2
+
(
c
1
⋅
c
1
−
1
)
2
+
(
c
2
⋅
c
2
−
1
)
2
+
(
c
3
⋅
c
3
−
1
)
2
(4)
\begin{aligned} \operatorname{Rot}(\mathbf{R})=&\left(\mathbf{c}_{1} \cdot \mathbf{c}_{2}\right)^{2}+\left(\mathbf{c}_{1} \cdot \mathbf{c}_{3}\right)^{2}+\left(\mathbf{c}_{2} \cdot \mathbf{c}_{3}\right)^{2}+\\ &\left(\mathbf{c}_{1} \cdot \mathbf{c}_{1}-1\right)^{2}+\left(\mathbf{c}_{2} \cdot \mathbf{c}_{2}-1\right)^{2}+\left(\mathbf{c}_{3} \cdot \mathbf{c}_{3}-1\right)^{2} \end{aligned} \tag4
Rot(R)=(c1⋅c2)2+(c1⋅c3)2+(c2⋅c3)2+(c1⋅c1−1)2+(c2⋅c2−1)2+(c3⋅c3−1)2(4)
公式 4 的
c
1
,
c
2
,
c
3
c_1, c_2, c_3
c1,c2,c3分别是
R
R
R的三个列向量,
E
r
o
t
E_{rot}
Erot是变形图中所有变化的旋转误差的总和:
E
r
o
t
=
∑
j
=
1
m
R
o
t
(
R
j
)
(5)
E_{rot}=\sum^m_{j=1}Rot(R_j) \tag5
Erot=j=1∑mRot(Rj)(5)
2.2.2 Regularization
每个仿射变换实际上代表了以某个node为中心的局部变形,由于相邻的变换会存在相互交叠的影响,我们必须保证所计算的变化彼此一致,因此加入正则化项。
E
r
e
g
=
∑
j
=
1
m
∑
k
∈
N
(
j
)
α
j
k
∥
R
j
(
g
k
−
g
j
)
+
g
j
+
t
j
−
(
g
k
+
t
k
)
∥
2
2
(6)
E_{reg}=\sum^m_{j=1} \sum_{k \in N(j)} \alpha_{jk} \Vert R_j(g_k-g_j)+g_j+t_j-(g_k+t_k)\Vert_2^2 \tag6
Ereg=j=1∑mk∈N(j)∑αjk∥Rj(gk−gj)+gj+tj−(gk+tk)∥22(6)
其中
α
j
k
=
1.0
\alpha_{jk}=1.0
αjk=1.0。上式说明,如果nodes
j
j
j和
k
k
k相邻,那么它们将影响统一公共区域。从
j
j
j的仿射变换所推断出的
k
k
k的位置(
R
j
(
g
k
−
g
j
)
+
g
j
+
t
j
R_j(g_k-g_j)+g_j+t_j
Rj(gk−gj)+gj+tj),应该与
k
k
k由自身变换得到的位置(
g
k
+
t
k
g_k+t_k
gk+tk)是相一致的,然后计算它们的误差,得到
E
r
e
g
E_{reg}
Ereg。
2.2.3 Constraint
约束分两种:
a) handle constraint,操作点约束,用户选择一部分点,操作这部分点进行变化;
b) fixed constraint,固定点约束,用户选择一部分点,固定这些点不发生变化。
操作点约束会影响计算,因此求约束误差项:
E
c
o
n
=
∑
l
=
1
p
∥
v
~
i
n
d
e
x
(
l
)
−
q
l
∥
2
2
(7)
E_{con} = \sum_{l=1}^p\Vert \tilde v_{index(l)} - q_l \Vert_2^2 \tag7
Econ=l=1∑p∥v~index(l)−ql∥22(7)
其中
v
~
i
n
d
e
x
(
l
)
\tilde v_{index(l)}
v~index(l)是由公式 2 计算得到的的点位置,
q
l
q_l
ql是用户变形得到的位置。
固定点约束不会影响计算,反而可以减少未知数的数量来加快计算速度。
图 2. 三个误差项不同组合作用在变形图上的效果。每个节点的四边形说明相应的仿射变换引起的变形。
图 2 在变形图将
g
4
g_4
g4移动到目标位置,添加不同约束项会影响相邻点的变化。只添加约束项,则只有
g
4
g_4
g4到达目标位置,邻点不会发生改变;使用约束项和正则化项的组合,会产生不自然的剪切变形;使用了旋转项、约束项和正则化项的能量方程,能够获得保留细节的变形。
2.2.4 Numerics
将旋转项、约束项和正则项求一个加权和,求出使能量方程最小的
R
R
R和
t
t
t的值,即得到对应每个顶点所需要做的变换。
min
R
1
,
t
1
…
R
m
,
t
m
w
rot
E
rot
+
w
reg
E
reg
+
w
con
E
con
(8)
\min _{\mathbf{R}_{1}, \mathbf{t}_{1} \ldots \mathbf{R}_{m}, \mathbf{t}_{m}} w_{\text {rot }} \mathrm{E}_{\text {rot }}+w_{\text {reg }} \mathrm{E}_{\text {reg }}+w_{\text {con }} \mathrm{E}_{\text {con }} \tag8
R1,t1…Rm,tmminwrot Erot +wreg Ereg +wcon Econ (8)
s
.
t
.
R
q
=
I
,
t
q
=
0
,
∀
q
∈
f
i
x
e
d
i
d
s
s.t.\ \mathbf{R}_{q}=\mathbf{I}, \mathbf{t}_{q}=\mathbf{0}, \forall q \in fixed\ ids
s.t. Rq=I,tq=0,∀q∈fixed ids
其中
w
r
o
t
=
1
,
w
r
e
g
=
10
,
w
c
o
n
=
100
w_{rot}=1,w_{reg}=10, w_{con}=100
wrot=1,wreg=10,wcon=100
上式表明要求出所有nodes的旋转矩阵和位移向量,对于固定点,旋转矩阵为单位矩阵,位移向量为0。
原文中公式 8使用高斯-牛顿法求解的,但针对非线性最小二乘问题,使用Google开发的非线性优化的库ceres能够更快地求解,而且精度也高。
3. 算法实现
3.1 实验环境
电脑:Ubuntu 20.04
编程环境:c++
使用库:libigl 用于读写文件和对vertices均匀采样得到nodes;用ceres库的莱文伯格-马夸特方法迭代求解公式 8,在参考文献1中有详细介绍如何使用这个库。
3.2 核心代码(可能参考意义不大,因为原cpp文件很多内容)
算法流程:
- 对三维模型进行均匀采样;
/*
* V: vertex of the surface
* F: faces of the surface
* N: nodes of the deformation graph
* E: edges of the deformation graph
*/
EmbeddedDeformation::EmbeddedDeformation(Eigen::MatrixXd V, Eigen::MatrixXi F)
{
nodes_connectivity = 4; // num of k-nearest nodes = 4
w_rot = 1;
w_reg = 10;
w_con = 100;
// define nodes as subset
Eigen::MatrixXd N;
N = igl::random_points_on_mesh(250, V,F); // 250 为采样后的nodes数量
// define edges
Eigen::MatrixXi E(N.rows()*(nodes_connectivity_+1), 2);
std::vector<int> closest_points;
int counter = 0;
for (int i = 0; i < N.rows(); ++i) {
closest_points = search_object.return_k_closest_points( indexes_of_deformation_graph_in_V_[i], nodes_connectivity_+2 );
for (int j=0; j< closest_points.size(); j++) {
E(counter, 0) = i;
E(counter, 1) = closest_points[j];
counter ++;
}
}
// define deformation graph
deformation_graph_ptr_ = new libgraphcpp::Graph(N, E);
}
- 构建deformation graph变形图,即计算边和点的连通性,得到邻接矩阵;
Graph::Graph(Eigen::MatrixXd nodes, Eigen::MatrixXi edges)
{
num_nodes_ = nodes_.rows();
num_edges_ = edges_.rows();
std::vector< std::vector<int> > adjacency_list_; // contains for each nodes, its nodes neighbors
std::vector< std::vector<int> > adjacency_edge_list_; // contains for each nodes, its edges neighbors
adjacency_list_.resize(num_nodes_);
adjacency_edge_list_.resize(num_nodes_);
for (int i=0; i<num_edges_; i++)
{
adjacency_list_[edges_(i, 0)].push_back(edges_(i, 1));
adjacency_list_[edges_(i, 1)].push_back(edges_(i, 0));
adjacency_edge_list_[edges_(i, 0)].push_back(i);
adjacency_edge_list_[edges_(i, 1)].push_back(i);
}
};
- 选择操控点并移动,并用ceres库优化求解旋转矩阵 R R R和位移向量 t t t,最后回代公式 2 得到每个vertex变形后的坐标。
// initialize the parameters that needs to be optimized (nodes rotations and translations)
std::vector< std::vector <double> > params;
params.resize(deformation_graph_ptr_->num_nodes());
for (int i = 0; i < params.size(); ++i)
{
params[i].resize(12);
// 初始化参数 R = params[i][0-8] 和 t = params[i][9-11]
params[i][0 ] = 1;
params[i][1 ] = 0;
params[i][2 ] = 0;
params[i][3 ] = 0;
params[i][4 ] = 1;
params[i][5 ] = 0;
params[i][6 ] = 0;
params[i][7 ] = 0;
params[i][8 ] = 1;
params[i][9 ] = 0;
params[i][10] = 0;
params[i][11] = 0;
}
// run levenberg marquardt optimization
ceres::Problem problem;
std::cout << "progress : add the residuals for Rotation CostFunction ..." << std::endl;
for (int i = 0; i < deformation_graph_ptr_->num_nodes(); ++i)
{
RotCostFunction* cost_function = new RotCostFunction(w_rot_); // equation (2)
problem.AddResidualBlock(cost_function, NULL, ¶ms[i][0]);
}
std::cout << "progress : add the residuals for Regularization CostFunction ..." << std::endl;
for (int i = 0; i < deformation_graph_ptr_->num_nodes(); ++i)
{
for (int j = 0; j < deformation_graph_ptr_->get_adjacency_list(i).size(); ++j)
{
int k = deformation_graph_ptr_->get_adjacency_list(i, j);
Eigen::Vector3d g_j = deformation_graph_ptr_->get_node(i);
Eigen::Vector3d g_k = deformation_graph_ptr_->get_node(k);
RegCostFunction* cost_function = new RegCostFunction(w_reg_, g_j, g_k);
// add residual block
problem.AddResidualBlock(cost_function, NULL, ¶ms[i][0], ¶ms[k][0]);
}
}
std::cout << "progress : add the residuals for Constraint CostFunction ..." << std::endl;
for (int i = 0; i < sources.rows(); ++i)
{
// stack the parameters
std::vector<double *> parameter_blocks;
std::vector<Eigen::Vector3d> vector_g;
// define cost function and associated parameters
for (int j = 0; j < nodes_connectivity_; ++j)
parameter_blocks.push_back(¶ms[ sources_nodes_neighbours[i][j] ][0]);
for (int j = 0; j < sources_nodes_neighbours[i].size(); ++j)
vector_g.push_back( deformation_graph_ptr_->get_node(sources_nodes_neighbours[i][j]) );
// create the cost function
ConCostFunction* cost_function = new ConCostFunction(w_con_, vector_g, sources.row(i), targets.row(i));
// add residual block
problem.AddResidualBlock(cost_function, NULL, parameter_blocks);
}
// Run the solver
ceres::Solver::Options options;
options.check_gradients = false;
options.gradient_check_relative_precision = 0.01;
options.minimizer_progress_to_stdout = verbose_; // 打印输出,相当于cout
options.num_threads = 1;
options.max_num_iterations = 100;
options.function_tolerance = 1e-10;
options.parameter_tolerance = 1e-10;
ceres::Solver::Summary summary; // 打印不断更新的优化参数的信息
Solve(options, &problem, &summary);
// Extract the optimized parameters
rotation_matrices_.resize(deformation_graph_ptr_->num_nodes());
translations_.resize(deformation_graph_ptr_->num_nodes());
for (int i = 0; i < deformation_graph_ptr_->num_nodes(); ++i)
{
rotation_matrices_[i] = Eigen::Map<Eigen::Matrix3d>( &(params[i][0]) );
translations_[i] = Eigen::Map<Eigen::Vector3d> ( &(params[i][0]) + 9);
}
// redefine the deformed mesh
Eigen::MatrixXd New_V = V_;
std::vector<double> w_j;
Eigen::Vector3d new_point;
for (int i = 0; i < V_.rows(); ++i)
{
// get closest nodes
std::vector< int > neighbours_nodes;
neighbours_nodes = deformation_graph_greedy_search->return_k_closest_points(i, nodes_connectivity_+1);
// equation (4)
for (int j = 0; j < nodes_connectivity_; ++j)
w_j.push_back( pow(1 - (V_.row(i) - deformation_graph_ptr_->get_node( neighbours_nodes[j] ).transpose() ).squaredNorm()
/ (V_.row(i) - deformation_graph_ptr_->get_node( neighbours_nodes.back() ).transpose() ).squaredNorm(), 2) );
// 权重归一化
double sum = 0;
for (int j = 0; j < nodes_connectivity_; ++j)
sum += w_j[j];
for (int j = 0; j < nodes_connectivity_; ++j)
w_j[j] /= sum;
// equation (2)
for (int j = 0; j < nodes_connectivity_; ++j)
new_vertex += w_j[j] * (rotation_matrices_[ neighbours_nodes[j] ] * (V_.row(i).transpose() - deformation_graph_ptr_->get_node( neighbours_nodes[j] ) )
+ deformation_graph_ptr_->get_node( neighbours_nodes[j] ) + translations_[ neighbours_nodes[j] ]);
New_V.row(i) = new_vertex;
}
}
3.3 算法结果
(a) 原始模型 (b) 兔子右边的耳朵向下弯的模型
图 3 算法结果
4. 心得体会
- 不要畏惧一门语言。我之前都是用python写的程序,因为我觉得c++要写变量的数据类型很麻烦,并且有很多指针,所以一直没用c++。但是在图形学上,由于进行大量矩阵运算,python的运行速度很慢,而c++运行速度比python快,而且c++有很多图形学优秀的库,比如libigl, Eigen, ceres。所以需要多练习来掌握一门语言。