Large-Scale Bounded Distortion Mappings 编写思路
Large-Scale Bounded Distortion Mappings 软件编写思路
目录
https://services.math.duke.edu/~shaharko//LargeScaleBD.html
这是Lipman组2015年的一篇SIggraph Asia文章
一、原文优化问题以及方法
1、paper中要处理的问题
\[\min_x \|Tx-Tx_0\|_2 \\
s.t. Ax=b\\
Tx \in D
\]
\(x\) 表示一个映射，一般是分片仿射变换。\(T\) 是微分算子。比如在一个三角形上，某一点\(v_i\) ,\(x_i(v_i)=A_iv_i+\delta_i \) ，那么\(T(x_i)=A_i\) 。然后第二项是线性约束，不解释了。第三项是给定的扭曲的界限，就是对Jacobi矩阵的奇异值的一些要求，文章中给的是Conformal distortion,也就是给定一个常数\(K\) ，然后\(A_i\) 的奇异值\(\sigma_{i,j}\) （从大到小排序）满足\(\sigma_{i,1}<K\sigma_{i,d}\) ，\(d\) 表示维度。
2、paper中解决此问题的方法!
paper中主要把原来的非线性不等式越是转换成线性等式约束。
\[\min_x \|Tx-Tx_0\|_2 \\
s.t. Ax=b\\
Tx \in D
\]
随着\(x\) 的变化，\(Tx\)形成一条线性结构\(R(T)\) 原问题就变成找到满足于等式约束的，并且在\(D\)与\(R(x)\)相交部分中与\(Tx_0\)最接近的一点。这样的话，可以由\(Tx_0\)找到\(D\)上的投影点（最近点），然后找到相切的超平面\(H_0\)然后再找超平面与\(R(T)\)交点，如此迭代几次，很快便可以收敛到\(Tx_*\)。
二、程序及实现细节
通过以上思路，可以把原问题转化成迭代计算线性结构的优化：
\[\min_x \|Tx-\Pi_D(Tx_0)\|_2 \\
s.t. Ax=b\\
Tx \in H
\]
1、给定初始值（参数化）
这里仅仅用三角网格参数化到二维之后，然后通过此方法降低扭曲的例子。
参数化之后的网格
2、初始映射的Jacobi矩阵向区域\(D\)中投影
矩阵\(C\)的投影计算如下：
\[\Pi_{D^K}(C)=U diag(\tau_1,\dots,\tau_d)V^T \\
\min_{\tau_i} \sum_{i=1}^d (\sigma_i-\tau_i)^2\\
s.t. \tau_1<K\tau_d
\]
其中\(K\)为共形扭曲的界限，\(\sigma_i\)为从大到小的符号奇异值。
二维的计算程序（使用Eigen库）：
Jacobi2 BoundedDistortion::Project_OneTri(Jacobi2 Z, double Dk, double &D_max)
{
Matrix2d U, V, W; Vector2d A; W.fill(0);
JacobiSVD<MatrixXd> svd(Z, ComputeThinU | ComputeThinV);
W.fill(0);
U = svd.matrixU();
V = svd.matrixV();
A = svd.singularValues();
if (D_max < (A(0) / A(1)))
{
D_max = (A(0) / A(1));
}
if (A[0] <= Dk*A[1]) return Z;
else
{
W(1, 1) = (Dk*A(0) + A(1)) / (Dk*Dk + 1);
W(0, 0) = Dk*W(1, 1);
Z = U*W*(V.transpose());
}
return Z;
}
3、等式约束的最小二乘问题等价于求解线性方程组
KKT线性方程
\[\begin{vmatrix}
T^TT & A^T & T^T\mathbf{n}_0\\
A & 0 & 0 \\
\mathbf{n}_0^TT & 0 & 0 \\
\end{vmatrix}
\begin{vmatrix}
x\\ \lambda\\ \mu
\end{vmatrix}
=
\begin{vmatrix}
T^TTx_0\\
b\\
\mathbf{n}_0^T \Pi_D(z_0)
\end{vmatrix}
\]
重新分块
\[\begin{vmatrix}
M & \eta_0\\
\eta_0^T & 0 \\
\end{vmatrix}
\begin{vmatrix}
x\\ \lambda\\ \mu
\end{vmatrix}
=
\begin{vmatrix}
c_0\\
d_0\\
\end{vmatrix}
\tag{1}
\]
可以把(1)分解可得分解成两个方程组：
\[M
\begin{vmatrix}
x\\ \lambda
\end{vmatrix}
+ \eta_0^T\mu
=
c_0
\tag{2.1}
\]
\[\eta_0^T
\begin{vmatrix}
x\\ \lambda\\
\end{vmatrix}
=
d_0
\tag{2.2}
\]
由(2.1)解出\(\begin{vmatrix} x\\ \lambda\end{vmatrix}\)代入(2.2)即可解出\(\mu\)，然后再代入(2.1)求解出\(x\)。
void BoundedDistortion::GetMatrix()
{
//主要计算和分解M和T矩阵
SparseMatrix<double >T(nf * 2 * 2, nv * 2), TT(nv * 2, nv * 2), M(nv * 2 + nc * 2, nv * 2 + nc * 2);
vector<Jacobi2>coef = para_->GetJacobiMatrix();
T_trip.clear();
int nf = para_->parameter_mesh.n_faces();
for (int i = 0; i < nf; i++)
{
Get_TMatrix_OneTri(T_trip, i,coef[i]);
}
T.setFromTriplets(T_trip.begin(), T_trip.end());//T
TT = T.transpose()*T;
M_trip = SparseToTriplet(TT);
M_trip = Get_linear_constrain(M_trip, nv * 2, 0);//加上线性约束
M.setFromTriplets(M_trip.begin(), M_trip.end());
solver_M.compute(M);
}vector<Triplet<double>> BoundedDistortion::Get_KKT_Matrix(vector<Triplet<double>> triple)
{//计算M的LU分解
int nf = map_mesh.n_faces();
int nv = map_mesh.n_vertices();
int nc = para_->DeliveryFixed().size();//约束点个数
SparseMatrix<double >T0(nf * 2 * 2, nv * 2), T(nv * 2, nv* 2), M(nf * 2 * 2 + nc * 2, nf * 2 * 2 + nc * 2);
vector<Trip> T_trip;
T_trip = Get_diff_operator(T_trip);
T0.setFromTriplets(T_trip.begin(), T_trip.end());//T
T = T0.transpose()*T0;
T_trip = SparseToTriplet(T);
auto CC=Get_linear_constrain(T_trip, nv * 2,0);
M.setFromTriplets(T_trip.begin(), T_trip.end());
solver_M.compute(M);
return T_trip;
}
vector<Jacobi2> ParaMeterization::GetJacobiMatrix()
{
coordinates_matrix.clear();
vector<Jacobi2> coefficient;
int nt = ptr_mesh->n_faces();
jacobi_matrix.clear();
for (int i=0;i<nt;i++)
{
OpenMesh::FaceHandle face = ptr_mesh->face_handle(i);
vector<Mesh::Normal >x, u; Mesh::Normal x_normal,u_normal;
x_normal = ptr_mesh->normal(face);
for (auto it = ptr_mesh->fv_ccwbegin(face); it != ptr_mesh->fv_ccwend(face);it++)
{
auto vertex = ptr_mesh->vertex_handle(it->idx());
x.push_back( ptr_mesh->point(vertex));
}
face = parameter_mesh.face_handle(i);
u_normal = parameter_mesh.normal(face);
for (auto it = parameter_mesh.fv_ccwbegin(face); it != parameter_mesh.fv_ccwend(face); it++)
{
auto vertex = ptr_mesh->vertex_handle(it->idx());
u.push_back(parameter_mesh.point(vertex));
}
Jacobi2 Ax,Au,J;
Coordinats3_2 A;
Ax(0, 0) = (x[1] - x[0]).length();
Ax(1, 0) = 0;
Ax(0, 1) = (x[2] - x[0]) | (x[1] - x[0]) / (x[1] - x[0]).length();
auto e = x_normal % (x[1] - x[0]);
Ax(1, 1) = (x[2] - x[0]) | e / (e).length();
Au(0, 0) = (u[1] - u[0])[0];
Au(1, 0) = (u[1] - u[0])[1];
Au(0, 1) = (u[2] - u[0])[0];
Au(1, 1) = (u[2] - u[0])[1];
//计算每个三角形的局部坐标系
Vector3d e1, e2;
for (int k = 0; k < 3;k++)
{
e1[k] = (x[0] - x[1])[k];
e2[k] = e[k];
}
A.col(0) = e1/e1.norm();
A.col(1) = e2/e2.norm();//每个面的局部坐标系的单位基底
coordinates_matrix.push_back(A);
coefficient.push_back(Ax.inverse());
J =Au *coefficient.back();
jacobi_matrix.push_back(J);
}
return coefficient;
}
vector<Triplet<double>> BoundedDistortion::SparseToTriplet(SparseMatrix<double> A)
{
vector<Triplet<double>> A_Triplet;
for (int k = 0; k < A.outerSize(); ++k)
{
for (SparseMatrix<double>::InnerIterator it(A, k); it; ++it)
{
A_Triplet.push_back(Triplet<double>(it.row(), it.col(), it.value()));
}
}
return A_Triplet;
}
vector<Triplet<double>> BoundedDistortion::Get_linear_constrain(vector<Triplet<double>> &triple, int n, int m)//A,from the nth
{
vector<int >linearconstrain_idx = para_->DeliveryFixed();
for (int i = 0; i < linearconstrain_idx.size(); i++)
{
//A
triple.push_back(Triplet<double>(n + i * 2, m + 2 * linearconstrain_idx[i], 1));
triple.push_back(Triplet<double>(n + i * 2 + 1, m + 2 * linearconstrain_idx[i] + 1, 1));
//AT
triple.push_back(Triplet<double>(m + 2 * linearconstrain_idx[i], n + i * 2, 1));
triple.push_back(Triplet<double>(m + 2 * linearconstrain_idx[i] + 1, n + i * 2 + 1, 1));
}
return triple;
}
以上函数执行可以得到M的分解以及矩阵T了，下面开始迭代求解
void BoundedDistortion::GetBoundedDistortionMap()
{
Initial();
VectorXd Xn(2 * nv), Xn1(2 * nv), X_temp(2 * nv);
for(int i = 0; i < nv; i++)
{
Xn(2 * i) = para_->parameter_mesh.point(para_->parameter_mesh.vertex_handle(i))[0];
Xn(2 * i + 1) = para_->parameter_mesh.point(para_->parameter_mesh.vertex_handle(i))[1];
}
GetDistortion(Xn);
for(int i = 0; i <100;i++)
{
VectorXd PIDTXn = Project(Xn);//投影系数矩阵
Xn1 = SolveSystem(Xn, PIDTXn);
Xn = Xn1;
if (MAX_distortion<=receive_distortion*1.03) break;
}
for (int i = 0; i < nv; i++)
{
para_->parameter_mesh.set_point(para_->parameter_mesh.vertex_handle(i), OpenMesh::Vec3d(Xn(2 * i), Xn(2 * i + 1), 0.0));
}
GetDistortion(Xn1);
*para_->ptr_mesh = para_->parameter_mesh;
}
VectorXd BoundedDistortion::SolveSystem(VectorXd Xn, VectorXd PIDTXn)
{
vector<int > contrains = para_->DeliveryFixed();//约束点个数
VectorXd b1(2 * nc), b0(2 * nv), n0(nf * 4), z0(nf * 4); double b2;
for (int i = 0; i < nc; i++)
{
b1(2 * i) = para_->parameter_mesh.point(para_->parameter_mesh.vertex_handle(contrains[i]))[0];
b1(2 * i + 1) = para_->parameter_mesh.point(para_->parameter_mesh.vertex_handle(contrains[i]))[1];
}
SparseMatrix<double >T(nf * 4, nv * 2);
T.setFromTriplets(T_trip.begin(), T_trip.end());
b0 = T.transpose()*T*Xn;
n0 = T*Xn-PIDTXn;
b2 = n0.dot(PIDTXn);
VectorXd c0(2 * nv + 2 * nc), eta(2 * nv + 2 * nc);
c0.segment(0, 2 * nv) = b0;
c0.segment(2 * nv, 2 * nc) = b1;
eta.fill(0);
eta.segment(0, 2 * nv) = T.transpose()*n0;
double mu = eta.dot(solver_M.solve(c0)) - b2;
mu = mu / (eta.dot(solver_M.solve(eta)));
VectorXd Mb = c0 - mu*eta;
VectorXd x_labda = solver_M.solve(Mb);
VectorXd x_ = x_labda.segment(0, 2 * nv);
return x_;
}
void BoundedDistortion::Initial()
{
nf = para_->parameter_mesh.n_faces();
nv = para_->parameter_mesh.n_vertices();
nc = para_->DeliveryFixed().size();//约束点个数
receive_distortion = 11;
K.resize(nf, receive_distortion);
GetMatrix();
}
结果展示
原映射：
优化后的映射：
