Large-Scale Bounded Distortion Mappings 编写思路

最新推荐文章于 2022-05-30 16:01:02 发布

weixin_30369041

最新推荐文章于 2022-05-30 16:01:02 发布

阅读量129

点赞数

原文链接：http://www.cnblogs.com/zhanghedong/p/7187486.html

版权

Large-Scale Bounded Distortion Mappings 软件编写思路
- 一、原文优化问题以及方法
  - 1、paper中要处理的问题
  - 2、paper中解决此问题的方法!
- 二、程序及实现细节

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();
}

结果展示

原映射：

优化后的映射：

转载于:https://www.cnblogs.com/zhanghedong/p/7187486.html

weixin_30369041

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
Large-Scale Bounded Distortion Mappings 编写思路

目录 Large-Scale Bounded Distortion Mappings 软件编写思路一、原文优化问题以及方法 1、paper中要处理的问题 2、paper中解决此问题的方法! 二、程序及实现细节 1、给定初始值（参数化） ...
复制链接

扫一扫