Libigl学习笔记——二、第四章——形状变形

第四章

  1. 现代基于网格的形状变形方法满足用户在手柄(网格上的选定顶点或区域)处的变形约束,并将这些手柄变形平滑地传播到形状的其余部分,而不会删除或扭曲细节。
  2. Libigl 提供各种最先进的变形技术实现,从基于二次网格的能量最小化器到蒙皮方法,再到非线性弹性启发技术。

4.1 双谐波变形 Biharmonic Deformation

  1. 2000年至2010年的研究期间产生了一系列技术,这些技术将基于手柄的形状变形问题视为二次能量最小化问题或等效于线性偏微分方程的解。
  2. 这些技术有很多种,但原型子集是那些考虑双拉普拉斯方程解的子集,即双谐波函数。这种四阶偏微分方程在边界条件下提供了足够的灵活性,以确保手柄约束(在细化极限中) 的 C^1 连续性。

4.1.1 双谐波曲面 Biharmonic Surfaces

  1. 让我们首先通过考虑双谐波表面来开始讨论双谐波变形。

  2. 我们将随意地将双谐波曲面定义为其位置函数相对于某些初始参数化为双谐波的曲面:
    在这里插入图片描述

  3. 并受一些句柄约束的约束,概念化为“边界条件”:
    在这里插入图片描述

  • 其中 x′ ,表面上点的未知 3D 位置。因此,我们要求每个空间坐标函数的双拉普拉斯函数为零。
  1. 在libigl中,可以解决双谐波问题 igl::harmonic 和设置 k=2 (双谐波):
// U_bc contains deformation of boundary vertices b
igl::harmonic(V,F,b,U_bc,2,U);
  1. 这将生成一个平滑的表面,用于插值手柄约束,但表面上的所有原始细节都将被平滑掉。最明显的是,如果原始表面尚未是双谐波,则为所有手柄提供恒等变形(使它们保持在静止位置)将不会重现原始表面。相反,结果将是双谐波表面,该表面确实插入了这些手柄位置。
  2. 因此,我们可以得出结论,这不是一种直观的形状变形技术。

4.1.2 双谐波变形场 Biharmonic Deformation Fields

  1. 现在我们知道变形技术的一个有用属性是“静止姿势再现”:对手柄不施加变形应该不对形状施加变形。

  2. 为了通过构造来保证这一点,我们可以处理变形场(即位移), d 而不是直接使用位置 x 。然后变形的位置可以恢复为:
    在这里插入图片描述

  3. 插入手柄约束的变形场的平滑变形场 d 将施加平滑变形形状 x′ 。当然,我们考虑双谐波变形场:
    在这里插入图片描述

  4. 受相同的句柄约束,但根据其边界(句柄)处的隐含变形场重写:
    在这里插入图片描述

  5. 同样,我们可以使用 igl::harmonic ,但这次求解变形场 k=2 ,然后恢复变形位置:

// U_bc contains deformation of boundary vertices b
D_bc = U_bc - igl::slice(V,b,1);
igl::harmonic(V,F,b,D_bc,2,D);
U = V+D;
  1. 401案例完整展示:
#include <igl/colon.h>
#include <igl/harmonic.h>
#include <igl/readOBJ.h>
#include <igl/readDMAT.h>
#include <igl/opengl/glfw/Viewer.h>
#include <algorithm>
#include <iostream>

double bc_frac = 1.0;
double bc_dir = -0.03;
bool deformation_field = false;
Eigen::MatrixXd V,U,V_bc,U_bc;
Eigen::VectorXd Z;
Eigen::MatrixXi F;
Eigen::VectorXi b;

bool pre_draw(igl::opengl::glfw::Viewer & viewer)
{
  using namespace Eigen;
  // Determine boundary conditions
  if(viewer.core().is_animating)
  {
    bc_frac += bc_dir;
    bc_dir *= (bc_frac>=1.0 || bc_frac<=0.0?-1.0:1.0);
  }

  const MatrixXd U_bc_anim = V_bc+bc_frac*(U_bc-V_bc);
  if(deformation_field)
  {
    MatrixXd D;
    MatrixXd D_bc = U_bc_anim - V_bc;
    igl::harmonic(V,F,b,D_bc,2,D);
    U = V+D;
  }else
  {
    igl::harmonic(V,F,b,U_bc_anim,2.,U);
  }
  viewer.data().set_vertices(U);
  viewer.data().compute_normals();
  return false;
}

bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
{
  switch(key)
  {
    case ' ':
      viewer.core().is_animating = !viewer.core().is_animating;
      return true;
    case 'D':
    case 'd':
      deformation_field = !deformation_field;
      return true;
  }
  return false;
}

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  igl::readOBJ(TUTORIAL_SHARED_PATH "/decimated-max.obj",V,F);
  U=V;
  // S(i) = j: j<0 (vertex i not in handle), j >= 0 (vertex i in handle j)
  VectorXi S;
  igl::readDMAT(TUTORIAL_SHARED_PATH "/decimated-max-selection.dmat",S);
  igl::colon<int>(0,V.rows()-1,b);
  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(),
   [&S](int i)->bool{return S(i)>=0;})-b.data());

  // Boundary conditions directly on deformed positions
  U_bc.resize(b.size(),V.cols());
  V_bc.resize(b.size(),V.cols());
  for(int bi = 0;bi<b.size();bi++)
  {
    V_bc.row(bi) = V.row(b(bi));
    switch(S(b(bi)))
    {
      case 0:
        // Don't move handle 0
        U_bc.row(bi) = V.row(b(bi));
        break;
      case 1:
        // move handle 1 down
        U_bc.row(bi) = V.row(b(bi)) + RowVector3d(0,-50,0);
        break;
      case 2:
      default:
        // move other handles forward
        U_bc.row(bi) = V.row(b(bi)) + RowVector3d(0,0,-25);
        break;
    }
  }

  // Pseudo-color based on selection
  MatrixXd C(F.rows(),3);
  RowVector3d purple(80.0/255.0,64.0/255.0,255.0/255.0);
  RowVector3d gold(255.0/255.0,228.0/255.0,58.0/255.0);
  for(int f = 0;f<F.rows();f++)
  {
    if( S(F(f,0))>=0 && S(F(f,1))>=0 && S(F(f,2))>=0)
    {
      C.row(f) = purple;
    }else
    {
      C.row(f) = gold;
    }
  }

  // Plot the mesh with pseudocolors
  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(U, F);
  viewer.data().show_lines = false;
  viewer.data().set_colors(C);
  viewer.core().trackball_angle = Eigen::Quaternionf(sqrt(2.0),0,sqrt(2.0),0);
  viewer.core().trackball_angle.normalize();
  viewer.callback_pre_draw = &pre_draw;
  viewer.callback_key_down = &key_down;
  //viewer.core().is_animating = true;
  viewer.core().animation_max_fps = 30.;
  cout<<
    "Press [space] to toggle deformation."<<endl<<
    "Press 'd' to toggle between biharmonic surface or displacements."<<endl;
  viewer.launch();
}

在这里插入图片描述双谐波变形示例将雕像的头部变形为双谐波表面(顶部)并使用双谐波位移(底部)

4.1.3 与“微分坐标”的关系和拉普拉斯曲面编辑 Relationship To “differential Coordinates” And Laplacian Surface Editing

  1. 双谐波函数(无论是位置还是位移)是双拉普拉斯方程的解,也是“拉普拉斯能量”的最小化器。例如,对于位移 d ,能量读取
    在这里插入图片描述

  2. 我们定义 Δd 简单地应用拉普拉斯坐标。

  3. 通过拉普拉斯(-Beltrami)算子的线性,我们可以用原始位置和未知位置 x x′=x−d 来重新表达这种能量:
    在这里插入图片描述

  4. 在索尔金等人的早期工作中,数量 Δx′ 和 Δx 被称为“微分坐标”。 因此,它们的变形(没有线性旋转)等效于双谐波变形场。

4.2 多谐波变形 Polyharmonic Deformation

  1. 我们可以通过考虑拉普拉斯量的不同幂来推广双谐波变形,从而产生一系列形式的偏微分方程:
    在这里插入图片描述

  2. 与 k∈1,2,3,… .的选择 k 决定了手柄的连续性水平。特别是,在边界处暗示, k=1 k=2 暗示, k=3 暗示 C^ 1、C^ 0 、C^ 2 和一般 k 暗示 C^(k−1) 。

int k = 2;// or 1,3,4,...
igl::harmonic(V,F,b,bc,k,Z);
  1. 402案例完整展示:
#include <igl/colon.h>
#include <igl/harmonic.h>
#include <igl/readOBJ.h>
#include <igl/opengl/glfw/Viewer.h>
#include <algorithm>
#include <iostream>

double z_max = 1.0;
double z_dir = -0.03;
int k = 2;
bool resolve = true;
Eigen::MatrixXd V,U;
Eigen::VectorXd Z;
Eigen::MatrixXi F;
Eigen::VectorXi b;
Eigen::VectorXd bc;

bool pre_draw(igl::opengl::glfw::Viewer & viewer)
{
  using namespace Eigen;
  if(resolve)
  {
    igl::harmonic(V,F,b,bc,k,Z);
    resolve = false;
  }
  U.col(2) = z_max*Z;
  viewer.data().set_vertices(U);
  viewer.data().compute_normals();
  if(viewer.core().is_animating)
  {
    z_max += z_dir;
    z_dir *= (z_max>=1.0 || z_max<=0.0?-1.0:1.0);
  }
  return false;
}

bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
{
  switch(key)
  {
    case ' ':
      viewer.core().is_animating = !viewer.core().is_animating;
      break;
    case '.':
      k++;
      k = (k>4?4:k);
      resolve = true;
      break;
    case ',':
      k--;
      k = (k<1?1:k);
      resolve = true;
      break;
  }
  return true;
}

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  igl::readOBJ(TUTORIAL_SHARED_PATH "/bump-domain.obj",V,F);
  U=V;
  // Find boundary vertices outside annulus
  typedef Matrix<bool,Dynamic,1> VectorXb;
  VectorXb is_outer = (V.rowwise().norm().array()-1.0)>-1e-15;
  VectorXb is_inner = (V.rowwise().norm().array()-0.15)<1e-15;
  VectorXb in_b = is_outer.array() || is_inner.array();
  igl::colon<int>(0,V.rows()-1,b);
  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(),
   [&in_b](int i)->bool{return in_b(i);})-b.data());
  bc.resize(b.size(),1);
  for(int bi = 0;bi<b.size();bi++)
  {
    bc(bi) = (is_outer(b(bi))?0.0:1.0);
  }


  // Pseudo-color based on selection
  MatrixXd C(F.rows(),3);
  RowVector3d purple(80.0/255.0,64.0/255.0,255.0/255.0);
  RowVector3d gold(255.0/255.0,228.0/255.0,58.0/255.0);
  for(int f = 0;f<F.rows();f++)
  {
    if( in_b(F(f,0)) && in_b(F(f,1)) && in_b(F(f,2)))
    {
      C.row(f) = purple;
    }else
    {
      C.row(f) = gold;
    }
  }

  // Plot the mesh with pseudocolors
  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(U, F);
  viewer.data().show_lines = false;
  viewer.data().set_colors(C);
  viewer.core().trackball_angle = Eigen::Quaternionf(0.81,-0.58,-0.03,-0.03);
  viewer.core().trackball_angle.normalize();
  viewer.callback_pre_draw = &pre_draw;
  viewer.callback_key_down = &key_down;
  viewer.core().is_animating = true;
  viewer.core().animation_max_fps = 30.;
  cout<<
    "Press [space] to toggle animation."<<endl<<
    "Press '.' to increase k."<<endl<<
    "Press ',' to decrease k."<<endl;
  viewer.launch();
}

在这里插入图片描述多谐波变形示例将平坦域(左)变形为凸块,作为各种 k 谐波偏微分方程的解

4.3 有界双谐波权重 Bounded Biharmonic Weights

  1. 在计算机动画中,形状变形通常被称为“蒙皮”。约束被提出为角色内部刚性“骨骼”的相对旋转。变形方法或蒙皮方法决定了角色表面(即其皮肤)应如何随着骨骼旋转而移动。

  2. 最流行的技术是线性混合蒙皮。形状上的每个点将其新位置计算为骨骼变换的线性组合:
    在这里插入图片描述

  3. 其中 wi(x) ,第 i 个骨骼的标量权重函数的计算值为 x 和 Ti 是作为 4×3 矩阵的骨变换。

  4. 这个公式是令人尴尬的并行的(一个点的计算不依赖于另一个点的计算需要的共享数据)。它通常作为顶点着色器实现。每个顶点的权重和休息位置作为顶点着色器属性发送,骨骼转换作为统一发送。然后在顶点着色器中转换顶点,以便及时进行渲染。

  5. 由于蒙皮公式是线性的(因此得名),我们可以将其写为矩阵乘法:
    在这里插入图片描述

  • 其中 X′ 是作为行向量的变形位置堆栈,是包含权重和休息位置的 n×m⋅dim 矩阵, T M 是 n×3 转置骨变换的 m⋅(dim+1)×dim 堆栈。
  1. 传统上,重量函数 wj 由熟练的索具专业人员手动绘制。现在存在现代技术,可以根据骨架的形状和描述(或一般的任何手柄结构,如笼子、点集合、选定区域等)自动计算权重函数。

  2. 有界双谐波权重就是这样一种将权重计算作为约束优化问题的技术。权重通过最小化熟悉的拉普拉斯能量来增强平滑度:
    在这里插入图片描述

  3. 受强制插值句柄约束的约束:
    在这里插入图片描述

  • 其中 Hi 是第 i 个句柄,以及强制执行非消极性、统一划分和鼓励稀疏性的约束:
    在这里插入图片描述
  1. 这是一个二次规划问题,libigl 使用其活动集合求解器或通过调用 Mosek 来解决它。
  2. 403示例完整展示
#include <igl/boundary_conditions.h>
#include <igl/colon.h>
#include <igl/column_to_quats.h>
#include <igl/directed_edge_parents.h>
#include <igl/forward_kinematics.h>
#include <igl/jet.h>
#include <igl/lbs_matrix.h>
#include <igl/deform_skeleton.h>
#include <igl/normalize_row_sums.h>
#include <igl/readDMAT.h>
#include <igl/readMESH.h>
#include <igl/readTGF.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/bbw.h>

#include <Eigen/Geometry>
#include <Eigen/StdVector>
#include <vector>
#include <algorithm>
#include <iostream>


typedef
  std::vector<Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> >
  RotationList;

const Eigen::RowVector3d sea_green(70./255.,252./255.,167./255.);
int selected = 0;
Eigen::MatrixXd V,W,U,C,M;
Eigen::MatrixXi T,F,BE;
Eigen::VectorXi P;
RotationList pose;
double anim_t = 1.0;
double anim_t_dir = -0.03;

bool pre_draw(igl::opengl::glfw::Viewer & viewer)
{
  using namespace Eigen;
  using namespace std;
  if(viewer.core().is_animating)
  {
    // Interpolate pose and identity
    RotationList anim_pose(pose.size());
    for(int e = 0;e<pose.size();e++)
    {
      anim_pose[e] = pose[e].slerp(anim_t,Quaterniond::Identity());
    }
    // Propagate relative rotations via FK to retrieve absolute transformations
    RotationList vQ;
    vector<Vector3d> vT;
    igl::forward_kinematics(C,BE,P,anim_pose,vQ,vT);
    const int dim = C.cols();
    MatrixXd T(BE.rows()*(dim+1),dim);
    for(int e = 0;e<BE.rows();e++)
    {
      Affine3d a = Affine3d::Identity();
      a.translate(vT[e]);
      a.rotate(vQ[e]);
      T.block(e*(dim+1),0,dim+1,dim) =
        a.matrix().transpose().block(0,0,dim+1,dim);
    }
    // Compute deformation via LBS as matrix multiplication
    U = M*T;

    // Also deform skeleton edges
    MatrixXd CT;
    MatrixXi BET;
    igl::deform_skeleton(C,BE,T,CT,BET);

    viewer.data().set_vertices(U);
    viewer.data().set_edges(CT,BET,sea_green);
    viewer.data().compute_normals();
    anim_t += anim_t_dir;
    anim_t_dir *= (anim_t>=1.0 || anim_t<=0.0?-1.0:1.0);
  }
  return false;
}

bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
{
  switch(key)
  {
    case ' ':
      viewer.core().is_animating = !viewer.core().is_animating;
      break;
    case '.':
      selected++;
      selected = std::min(std::max(selected,0),(int)W.cols()-1);
      viewer.data().set_data(W.col(selected));
      break;
    case ',':
      selected--;
      selected = std::min(std::max(selected,0),(int)W.cols()-1);
      viewer.data().set_data(W.col(selected));
      break;
  }
  return true;
}

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  igl::readMESH(TUTORIAL_SHARED_PATH "/hand.mesh",V,T,F);
  U=V;
  igl::readTGF(TUTORIAL_SHARED_PATH "/hand.tgf",C,BE);
  // retrieve parents for forward kinematics
  igl::directed_edge_parents(BE,P);

  // Read pose as matrix of quaternions per row
  MatrixXd Q;
  igl::readDMAT(TUTORIAL_SHARED_PATH "/hand-pose.dmat",Q);
  igl::column_to_quats(Q,pose);
  assert(pose.size() == BE.rows());

  // List of boundary indices (aka fixed value indices into VV)
  VectorXi b;
  // List of boundary conditions of each weight function
  MatrixXd bc;
  igl::boundary_conditions(V,T,C,VectorXi(),BE,MatrixXi(),b,bc);

  // compute BBW weights matrix
  igl::BBWData bbw_data;
  // only a few iterations for sake of demo
  bbw_data.active_set_params.max_iter = 8;
  bbw_data.verbosity = 2;
  if(!igl::bbw(V,T,b,bc,bbw_data,W))
  {
    return EXIT_FAILURE;
  }

  //MatrixXd Vsurf = V.topLeftCorner(F.maxCoeff()+1,V.cols());
  //MatrixXd Wsurf;
  //if(!igl::bone_heat(Vsurf,F,C,VectorXi(),BE,MatrixXi(),Wsurf))
  //{
  //  return false;
  //}
  //W.setConstant(V.rows(),Wsurf.cols(),1);
  //W.topLeftCorner(Wsurf.rows(),Wsurf.cols()) = Wsurf = Wsurf = Wsurf = Wsurf;

  // Normalize weights to sum to one
  igl::normalize_row_sums(W,W);
  // precompute linear blend skinning matrix
  igl::lbs_matrix(V,W,M);

  // Plot the mesh with pseudocolors
  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(U, F);
  viewer.data().set_data(W.col(selected));
  viewer.data().set_edges(C,BE,sea_green);
  viewer.data().show_lines = false;
  viewer.data().show_overlay_depth = false;
  viewer.data().line_width = 1;
  viewer.callback_pre_draw = &pre_draw;
  viewer.callback_key_down = &key_down;
  viewer.core().is_animating = false;
  viewer.core().animation_max_fps = 30.;
  cout<<
    "Press '.' to show next weight function."<<endl<<
    "Press ',' to show previous weight function."<<endl<<
    "Press [space] to toggle animation."<<endl;
  viewer.launch();
  return EXIT_SUCCESS;
}

在这里插入图片描述示例 BoundedBiharmonicWeights 计算给定骨架的四面体网格的权重(顶部),然后对线性混合蒙皮变形(底部)进行动画处理。

4.4 双四元数蒙皮 Dual Quaternion Skinning

  1. 即使使用高质量的砝码,线性混合蒙皮也是有限的。特别是,它受到混合旋转作为矩阵的已知伪影的影响:旋转矩阵的权重组合不一定是旋转。考虑围绕 -轴旋转 −π/2 和绕 π/2 z -轴旋转之间的等混合。直观地说,人们可能期望得到单位矩阵,但混合是一个退化矩阵,将 x and y 坐标缩放为零:
    在这里插入图片描述

  2. 在实践中,这意味着形状在骨重重叠的区域收缩和塌陷:关节附近。

  3. 双四元数蒙皮提出了一种解决方案。此方法将刚性变换表示为一对单位四元数 。 ^q 线性混合蒙皮公式替换为双四元数的线性混合:
    在这里插入图片描述

  4. 其中 ^qi 是骨 i 刚性变换的双四元数表示。归一化迫使线性混合的结果再次成为单位对偶四元数,因此也是刚性变换。

  5. 与线性混合蒙皮一样,双四元数蒙皮最好在顶点着色器中执行。唯一的区别是骨变换作为双四元数而不是仿射变换矩阵发送。Libigl 使用该函数支持 CPU 端双四元数蒙皮,该 igl::dqs 函数采用更传统的刚性变换表示作为输入,并在混合之前在内部转换为双四元数表示:

// vQ is a list of rotations as quaternions
// vT is a list of translations
igl::dqs(V,W,vQ,vT,U);
  1. 404示例完整展示
#include <igl/directed_edge_orientations.h>
#include <igl/directed_edge_parents.h>
#include <igl/forward_kinematics.h>
#include <igl/PI.h>
#include <igl/lbs_matrix.h>
#include <igl/deform_skeleton.h>
#include <igl/dqs.h>
#include <igl/readDMAT.h>
#include <igl/readOBJ.h>
#include <igl/readTGF.h>
#include <igl/opengl/glfw/Viewer.h>

#include <Eigen/Geometry>
#include <Eigen/StdVector>
#include <vector>
#include <algorithm>
#include <iostream>


typedef 
  std::vector<Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> >
  RotationList;

const Eigen::RowVector3d sea_green(70./255.,252./255.,167./255.);
Eigen::MatrixXd V,W,C,U,M;
Eigen::MatrixXi F,BE;
Eigen::VectorXi P;
std::vector<RotationList > poses;
double anim_t = 0.0;
double anim_t_dir = 0.015;
bool use_dqs = false;
bool recompute = true;

bool pre_draw(igl::opengl::glfw::Viewer & viewer)
{
  using namespace Eigen;
  using namespace std;
  if(recompute)
  {
    // Find pose interval
    const int begin = (int)floor(anim_t)%poses.size();
    const int end = (int)(floor(anim_t)+1)%poses.size();
    const double t = anim_t - floor(anim_t);

    // Interpolate pose and identity
    RotationList anim_pose(poses[begin].size());
    for(int e = 0;e<poses[begin].size();e++)
    {
      anim_pose[e] = poses[begin][e].slerp(t,poses[end][e]);
    }
    // Propagate relative rotations via FK to retrieve absolute transformations
    RotationList vQ;
    vector<Vector3d> vT;
    igl::forward_kinematics(C,BE,P,anim_pose,vQ,vT);
    const int dim = C.cols();
    MatrixXd T(BE.rows()*(dim+1),dim);
    for(int e = 0;e<BE.rows();e++)
    {
      Affine3d a = Affine3d::Identity();
      a.translate(vT[e]);
      a.rotate(vQ[e]);
      T.block(e*(dim+1),0,dim+1,dim) =
        a.matrix().transpose().block(0,0,dim+1,dim);
    }
    // Compute deformation via LBS as matrix multiplication
    if(use_dqs)
    {
      igl::dqs(V,W,vQ,vT,U);
    }else
    {
      U = M*T;
    }

    // Also deform skeleton edges
    MatrixXd CT;
    MatrixXi BET;
    igl::deform_skeleton(C,BE,T,CT,BET);
    
    viewer.data().set_vertices(U);
    viewer.data().set_edges(CT,BET,sea_green);
    viewer.data().compute_normals();
    if(viewer.core().is_animating)
    {
      anim_t += anim_t_dir;
    }
    else
    {
      recompute=false;
    }
  }
  return false;
}

bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
{
  recompute = true;
  switch(key)
  {
    case 'D':
    case 'd':
      use_dqs = !use_dqs;
      return true;
    case ' ':
      viewer.core().is_animating = !viewer.core().is_animating;
      return true;
  }
  return false;
}

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  igl::readOBJ(TUTORIAL_SHARED_PATH "/arm.obj",V,F);
  U=V;
  igl::readTGF(TUTORIAL_SHARED_PATH "/arm.tgf",C,BE);
  // retrieve parents for forward kinematics
  igl::directed_edge_parents(BE,P);
  RotationList rest_pose;
  igl::directed_edge_orientations(C,BE,rest_pose);
  poses.resize(4,RotationList(4,Quaterniond::Identity()));
  // poses[1] // twist
  const Quaterniond twist(AngleAxisd(igl::PI,Vector3d(1,0,0)));
  poses[1][2] = rest_pose[2]*twist*rest_pose[2].conjugate();
  const Quaterniond bend(AngleAxisd(-igl::PI*0.7,Vector3d(0,0,1)));
  poses[3][2] = rest_pose[2]*bend*rest_pose[2].conjugate();

  igl::readDMAT(TUTORIAL_SHARED_PATH "/arm-weights.dmat",W);
  igl::lbs_matrix(V,W,M);

  // Plot the mesh with pseudocolors
  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(U, F);
  viewer.data().set_edges(C,BE,sea_green);
  viewer.data().show_lines = false;
  viewer.data().show_overlay_depth = false;
  viewer.data().line_width = 1;
  viewer.core().trackball_angle.normalize();
  viewer.callback_pre_draw = &pre_draw;
  viewer.callback_key_down = &key_down;
  viewer.core().is_animating = false;
  viewer.core().camera_zoom = 2.5;
  viewer.core().animation_max_fps = 30.;
  cout<<"Press [d] to toggle between LBS and DQS"<<endl<<
    "Press [space] to toggle animation"<<endl;
  viewer.launch();
}

在这里插入图片描述示例 DualQuaternionSkinning 将线性混合蒙皮(顶部)与双四元数蒙皮(底部)进行了比较,突出显示了 LBS 的糖果包装效果(中)和关节塌陷(右)

4.5 尽可能刚性 As-rigid-as-possible

  1. 蒙皮和其他变形线性方法本质上是有限的。特别是当手柄约束施加大旋转时,会出现困难。
  2. 在能量最小化方法的背景下,问题源于比较未变形形状坐标框架中的位置(我们的位移)。这些二次能量充其量是不变的,而不是平滑地改变局部旋转。因此,线性技术不会产生非平凡的弯曲和扭曲。
  3. 此外,在考虑实体形状(例如,使用四面体网格离散化)时,线性方法难以保持局部体积,并且经常会出现收缩和膨胀的伪影。
  4. 非线性变形技术为这些问题提供了解决方案。它们的工作原理是将网格顶点的变形与其旋转到最匹配变形的新坐标系的静止位置进行比较。非线性源于变形和最佳拟合旋转的相互依赖关系。这些技术通常被标记为“尽可能刚性”,因为它们会惩罚所有局部变形与旋转偏差的总和。
  5. 为了得到这样的能量,让我们考虑一个简单的每三角形能量:
    在这里插入图片描述
  6. 网格 X′ 的未知变形顶点位置在哪里,是三角形列表中的三角形,是三角形的面积 t T , at t {i,j} 是三角形中的边 t 。因此,该能量测量原始网格中的边缘矢量和未知网格 (xi−xj) (x′i−x′j) 之间的变化范数。
  7. 这种能量不是旋转不变的。如果我们将网格旋转 90 度,尽管整体变形是完全刚性的,但与旋转轴不对齐的边缘矢量的变化将很大。
  8. 因此,“尽可能刚性”的解决方案是为每个三角形附加辅助变量,这些变量 Rt t 被约束为旋转。然后重写能量,这次将变形的边缘向量与旋转的静止矢量进行比较:
    在这里插入图片描述
  9. 分离到主顶点位置变量 X′ 和旋转也 {R1,…,R|T|} 导致了优化策略。如果旋转 {R1,…,R|T|} 保持固定,则能量在其余变量 X′ 中是二次的,可以通过求解(稀疏)全局线性系统来优化。或者,如果保持固定,则 X′ 每次旋转都是局部Procrustes问题的解决方案(通过SVD或极性分解发现 3×3 )。这两个步骤——局部和全局——每个步骤都会微弱地降低能量,因此我们可以安全地迭代它们,直到收敛。
  10. “尽可能刚性”的不同风格取决于域和边缘集 T 的维度和共维度。Sorkine和Alexa提出的表面处理技术 22 被认为是 T 从每个顶点(辐条)发出的一组边。后来,Chao等人推导出了“尽可能刚性”网格能量与共旋转弹性之间的关系,考虑了0余维元素作为边缘集:2D中的三角形和3D 12 中的四面体。他们还展示了Sorkine和Alexa的边集如何不是连续能量的离散化,而是为包含入射在顶点(辐条和轮辋)上的元素的所有边的表面提出了边集。他们表明,这相当于测量弯曲,尽管是以离散化依赖的方式。
  11. libigl,支持这些常见的口味。选择一个是在预计算阶段之前设置能量类型的问题:
igl::ARAPData arap_data;
arap_data.energy = igl::ARAP_ENERGY_TYPE_SPOKES;
//arap_data.energy = igl::ARAP_ENERGY_TYPE_SPOKES_AND_RIMS;
//arap_data.energy = igl::ARAP_ENERGY_TYPE_ELEMENTS; //triangles or tets
igl::arap_precomputation(V,F,dim,b,arap_data);
  1. 就像 igl::min_quad_with_fixed_* ,这个预计算阶段只取决于网格、固定的顶点索引 b 和能量参数。为了解决对顶点 b 位置的某些约束,我们可以调用:
igl::arap_solve(bc,arap_data,U);
  • U 用作初始猜测,然后将解决方案计算到其中。
  1. Libigl 的刚性尽可能变形的实现利用了 McAdams 等人高度优化的奇异值分解代码, 该代码利用了 SSE 内部函数。
  2. 局部刚性的概念将在表面参数化的背景下重新讨论。
  3. 405案例完整展示:
#include <igl/colon.h>
#include <igl/directed_edge_orientations.h>
#include <igl/directed_edge_parents.h>
#include <igl/forward_kinematics.h>
#include <igl/PI.h>
#include <igl/lbs_matrix.h>
#include <igl/deform_skeleton.h>
#include <igl/dqs.h>
#include <igl/readDMAT.h>
#include <igl/readOFF.h>
#include <igl/arap.h>
#include <igl/opengl/glfw/Viewer.h>

#include <Eigen/Geometry>
#include <Eigen/StdVector>
#include <vector>
#include <algorithm>
#include <iostream>


typedef 
  std::vector<Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> >
  RotationList;

const Eigen::RowVector3d sea_green(70./255.,252./255.,167./255.);
Eigen::MatrixXd V,U;
Eigen::MatrixXi F;
Eigen::VectorXi S,b;
Eigen::RowVector3d mid;
double anim_t = 0.0;
double anim_t_dir = 0.03;
igl::ARAPData arap_data;

bool pre_draw(igl::opengl::glfw::Viewer & viewer)
{
  using namespace Eigen;
  using namespace std;
    MatrixXd bc(b.size(),V.cols());
    for(int i = 0;i<b.size();i++)
    {
      bc.row(i) = V.row(b(i));
      switch(S(b(i)))
      {
        case 0:
        {
          const double r = mid(0)*0.25;
          bc(i,0) += r*sin(0.5*anim_t*2.*igl::PI);
          bc(i,1) -= r+r*cos(igl::PI+0.5*anim_t*2.*igl::PI);
          break;
        }
        case 1:
        {
          const double r = mid(1)*0.15;
          bc(i,1) += r+r*cos(igl::PI+0.15*anim_t*2.*igl::PI);
          bc(i,2) -= r*sin(0.15*anim_t*2.*igl::PI);
          break;
        }
        case 2:
        {
          const double r = mid(1)*0.15;
          bc(i,2) += r+r*cos(igl::PI+0.35*anim_t*2.*igl::PI);
          bc(i,0) += r*sin(0.35*anim_t*2.*igl::PI);
          break;
        }
        default:
          break;
      }
    }
    igl::arap_solve(bc,arap_data,U);
    viewer.data().set_vertices(U);
    viewer.data().compute_normals();
  if(viewer.core().is_animating)
  {
    anim_t += anim_t_dir;
  }
  return false;
}

bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
{
  switch(key)
  {
    case ' ':
      viewer.core().is_animating = !viewer.core().is_animating;
      return true;
  }
  return false;
}

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  igl::readOFF(TUTORIAL_SHARED_PATH "/decimated-knight.off",V,F);
  U=V;
  igl::readDMAT(TUTORIAL_SHARED_PATH "/decimated-knight-selection.dmat",S);

  // vertices in selection
  igl::colon<int>(0,V.rows()-1,b);
  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(), 
   [](int i)->bool{return S(i)>=0;})-b.data());
  // Centroid
  mid = 0.5*(V.colwise().maxCoeff() + V.colwise().minCoeff());
  // Precomputation
  arap_data.max_iter = 100;
  igl::arap_precomputation(V,F,V.cols(),b,arap_data);

  // Set color based on selection
  MatrixXd C(F.rows(),3);
  RowVector3d purple(80.0/255.0,64.0/255.0,255.0/255.0);
  RowVector3d gold(255.0/255.0,228.0/255.0,58.0/255.0);
  for(int f = 0;f<F.rows();f++)
  {
    if( S(F(f,0))>=0 && S(F(f,1))>=0 && S(F(f,2))>=0)
    {
      C.row(f) = purple;
    }else
    {
      C.row(f) = gold;
    }
  }

  // Plot the mesh with pseudocolors
  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(U, F);
  viewer.data().set_colors(C);
  viewer.callback_pre_draw = &pre_draw;
  viewer.callback_key_down = &key_down;
  viewer.core().is_animating = false;
  viewer.core().animation_max_fps = 30.;
  cout<<
    "Press [space] to toggle animation"<<endl;
  viewer.launch();
}

在这里插入图片描述示例AsRigidAsPossible使表面变形,就好像它是由弹性材料制成的一样

4.6 快速自动蒙皮转换 Fast Automatic Skinning Transformations

  1. 不出所料,非线性优化比其线性表亲慢。在尽可能刚性优化的情况下,瓶颈通常是恢复每个边集(即每个三角形、四面体或顶点单元)的最佳拟合旋转所需的大量极分解。即使此代码经过优化,尽管变形具有低频行为,但初级自由度的数量仍与离散化水平相关联。
  2. 这为快速非线性优化提供了两条途径。首先,是否有必要(甚至有利)找到这么多最合适的轮换?其次,我们能否降低自由度以更好地反映所需变形的频率。
  3. 反过来,这些优化最终形成了一种方法,该方法优化了由高质量权重(即手动绘制的权重或有界双谐波权重)跨越的线性混合蒙皮变形空间。此空间是所有可能的网格变形的低维子空间,通过以矩阵形式编写线性混合蒙皮来捕获:
    在这里插入图片描述
  • 其中, n×3 矩阵 X′ 中的网格顶点位置被转置“句柄”变换 (3+1)m×3 堆栈中少量自由度的线性组合所取代。换入 MT X′ 上述 ARAP 能量会立即看到全局求解步骤中的性能提升为 m<<n 。
  1. 局部步长(拟合旋转)的复杂性仍然与原始网格离散化有关。但是,如果蒙皮表现良好,我们可以假设具有相似蒙皮权重的形状上的位置将类似地变形,从而意味着相似的最佳拟合旋转。因此,我们根据边集在权重空间中的表示对边集进行聚类:其中顶点 x 取坐标 [w1(x),w2(x),…,wm(x)] 。聚类边缘集的数量显示变形质量的收益递减,因此我们可以选择少量的聚类,与蒙皮权重函数的数量(而不是离散网格顶点的数量)成正比。

  2. 这种所提出的变形模型,可以同时看作是ARAP的快速子空间优化,以及寻找最佳蒙皮变换自由度的自动方法。

  3. 通过与句柄关联的外观转换上的线性相等约束,支持各种用户界面。要完全修复转换,我们只需添加约束:
    在这里插入图片描述在这里插入图片描述

  4. 为了仅修复句柄的原点,我们添加一个约束,要求变换在空间中插值一个点(通常是所有点的质心,其中 wi=1 :
    在这里插入图片描述

  • 其中 cT ,静止点在转置齐次坐标中的位置, 1×(3+1) 以及 c^(′T) 用户给出的点。
  1. 我们同样可以在句柄上仅修复转换的线性部分,从而释放翻译组件(产生“鸡头”效果):
    在这里插入图片描述在这里插入图片描述

  2. 最后,我们可以允许用户完全释放变换的自由度,委托优化为所有元素找到最佳值。为此,我们只需避免添加相应的约束。

4.6.1 带有分组边缘集的 ArapArap With Grouped Edge-sets

  1. 作为一种子空间方法,一个直接的缺点是自由度降低。这带来了性能,但在某些情况下会过多地限制行为。在这种情况下,可以使用蒙皮子空间来构建旋转边集的有效聚类,以实现传统的 ARAP 优化:放弃子空间替换。这具有双重效果。旋转拟合的成本,局部步长大幅降低,变形根据集群“正则化”。从高层次的角度来看,如果聚类是从蒙皮权重派生的,那么它们将阻止弯曲,特别是沿着权重函数的等值线。如果事先不知道句柄,也可以根据“测地线嵌入”(如双谐波距离嵌入)进行聚类。
  2. 从这个角度来看,我们可以将“辐条+轮辋”风格的曲面 ARAP 视为每个三角形边集的(轻微和冗余)聚类。
  3. 406案例完整展示
#include <igl/colon.h>
#include <igl/directed_edge_orientations.h>
#include <igl/directed_edge_parents.h>
#include <igl/forward_kinematics.h>
#include <igl/PI.h>
#include <igl/partition.h>
#include <igl/mat_max.h>
#include <igl/lbs_matrix.h>
#include <igl/slice.h>
#include <igl/deform_skeleton.h>
#include <igl/dqs.h>
#include <igl/lbs_matrix.h>
#include <igl/columnize.h>
#include <igl/readDMAT.h>
#include <igl/readOBJ.h>
#include <igl/arap.h>
#include <igl/arap_dof.h>
#include <igl/opengl/glfw/Viewer.h>

#include <Eigen/Geometry>
#include <Eigen/StdVector>
#include <vector>
#include <algorithm>
#include <iostream>


typedef 
  std::vector<Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> >
  RotationList;

const Eigen::RowVector3d sea_green(70./255.,252./255.,167./255.);
Eigen::MatrixXd V,U,M;
Eigen::MatrixXi F;
Eigen::VectorXi S,b;
Eigen::MatrixXd L;
Eigen::RowVector3d mid;
double anim_t = 0.0;
double anim_t_dir = 0.03;
double bbd = 1.0;
bool resolve = true;
igl::ARAPData arap_data,arap_grouped_data;
igl::ArapDOFData<Eigen::MatrixXd,double> arap_dof_data;
Eigen::SparseMatrix<double> Aeq;

enum ModeType
{
  MODE_TYPE_ARAP = 0,
  MODE_TYPE_ARAP_GROUPED = 1,
  MODE_TYPE_ARAP_DOF = 2,
  NUM_MODE_TYPES = 4
} mode = MODE_TYPE_ARAP;

bool pre_draw(igl::opengl::glfw::Viewer & viewer)
{
  using namespace Eigen;
  using namespace std;
  if(resolve)
  {
    MatrixXd bc(b.size(),V.cols());
    VectorXd Beq(3*b.size());
    for(int i = 0;i<b.size();i++)
    {
      bc.row(i) = V.row(b(i));
      switch(i%4)
      {
        case 2:
          bc(i,0) += 0.15*bbd*sin(0.5*anim_t);
          bc(i,1) += 0.15*bbd*(1.-cos(0.5*anim_t));
          break;
        case 1:
          bc(i,1) += 0.10*bbd*sin(1.*anim_t*(i+1));
          bc(i,2) += 0.10*bbd*(1.-cos(1.*anim_t*(i+1)));
          break;
        case 0:
          bc(i,0) += 0.20*bbd*sin(2.*anim_t*(i+1));
          break;
      }
      Beq(3*i+0) = bc(i,0);
      Beq(3*i+1) = bc(i,1);
      Beq(3*i+2) = bc(i,2);
    }
    switch(mode)
    {
      default:
        assert("unknown mode");
      case MODE_TYPE_ARAP:
        igl::arap_solve(bc,arap_data,U);
        break;
      case MODE_TYPE_ARAP_GROUPED:
        igl::arap_solve(bc,arap_grouped_data,U);
        break;
      case MODE_TYPE_ARAP_DOF:
      {
        VectorXd L0 = L;
        arap_dof_update(arap_dof_data,Beq,L0,30,0,L);
        const auto & Ucol = M*L;
        U.col(0) = Ucol.block(0*U.rows(),0,U.rows(),1);
        U.col(1) = Ucol.block(1*U.rows(),0,U.rows(),1);
        U.col(2) = Ucol.block(2*U.rows(),0,U.rows(),1);
        break;
      }
    }
    viewer.data().set_vertices(U);
    viewer.data().set_points(bc,sea_green);
    viewer.data().compute_normals();
    if(viewer.core().is_animating)
    {
      anim_t += anim_t_dir;
    }else
    {
      resolve = false;
    }
  }
  return false;
}

bool key_down(igl::opengl::glfw::Viewer &viewer, unsigned char key, int mods)
{
  switch(key)
  {
    case '0':
      anim_t = 0;
      resolve = true;
      return true;
    case '.':
      mode = (ModeType)(((int)mode+1)%((int)NUM_MODE_TYPES-1));
      resolve = true;
      return true;
    case ',':
      mode = (ModeType)(((int)mode-1)%((int)NUM_MODE_TYPES-1));
      resolve = true;
      return true;
    case ' ':
      viewer.core().is_animating = !viewer.core().is_animating;
      if(viewer.core().is_animating)
      {
        resolve = true;
      }
      return true;
  }
  return false;
}

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;
  igl::readOBJ(TUTORIAL_SHARED_PATH "/armadillo.obj",V,F);
  U=V;
  MatrixXd W;
  igl::readDMAT(TUTORIAL_SHARED_PATH "/armadillo-weights.dmat",W);
  igl::lbs_matrix_column(V,W,M);

  // Cluster according to weights
  VectorXi G;
  {
    VectorXi S;
    VectorXd D;
    igl::partition(W,50,G,S,D);
  }

  // vertices corresponding to handles (those with maximum weight)
  {
    VectorXd maxW;
    igl::mat_max(W,1,maxW,b);
  }

  // Precomputation for FAST
  cout<<"Initializing Fast Automatic Skinning Transformations..."<<endl;
  // number of weights
  const int m = W.cols();
  Aeq.resize(m*3,m*3*(3+1));
  vector<Triplet<double> > ijv;
  for(int i = 0;i<m;i++)
  {
    RowVector4d homo;
    homo << V.row(b(i)),1.;
    for(int d = 0;d<3;d++)
    {
      for(int c = 0;c<(3+1);c++)
      {
        ijv.push_back(Triplet<double>(3*i + d,i + c*m*3 + d*m, homo(c)));
      }
    }
  }
  Aeq.setFromTriplets(ijv.begin(),ijv.end());
  igl::arap_dof_precomputation(V,F,M,G,arap_dof_data);
  igl::arap_dof_recomputation(VectorXi(),Aeq,arap_dof_data);
  // Initialize
  MatrixXd Istack = MatrixXd::Identity(3,3+1).replicate(1,m);
  igl::columnize(Istack,m,2,L);

  // Precomputation for ARAP
  cout<<"Initializing ARAP..."<<endl;
  arap_data.max_iter = 1;
  igl::arap_precomputation(V,F,V.cols(),b,arap_data);
  // Grouped arap
  cout<<"Initializing ARAP with grouped edge-sets..."<<endl;
  arap_grouped_data.max_iter = 2;
  arap_grouped_data.G = G;
  igl::arap_precomputation(V,F,V.cols(),b,arap_grouped_data);


  // bounding box diagonal
  bbd = (V.colwise().maxCoeff()- V.colwise().minCoeff()).norm();

  // Plot the mesh with pseudocolors
  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(U, F);
  viewer.data().add_points(igl::slice(V,b,1),sea_green);
  viewer.data().show_lines = false;
  viewer.callback_pre_draw = &pre_draw;
  viewer.callback_key_down = &key_down;
  viewer.core().is_animating = false;
  viewer.core().animation_max_fps = 30.;
  cout<<
    "Press [space] to toggle animation."<<endl<<
    "Press '0' to reset pose."<<endl<<
    "Press '.' to switch to next deformation method."<<endl<<
    "Press ',' to switch to previous deformation method."<<endl;
  viewer.launch();
}

在这里插入图片描述示例 FastAutomaticSkinningTransformations 将详细形状上的完全(慢速)ARAP 变形(中间左侧)与具有分组旋转边缘集的 ARAP (中间右侧)与非常快速的子形方法(右)进行了比较

4.7 双谐波坐标 Biharmonic Coordinates

  1. 线性混合蒙皮(如上所述)通过权重将手柄(骨骼、点、区域等)处的完整仿射变换传播到形状的其余部分,从而使网格变形。另一个变形框架,称为“广义重心坐标”,是线性混合蒙皮的一个特例 16 :变换仅限于纯平移,并且需要权重以保持仿射精度。后一个要求意味着我们可以将网格中任何顶点的静止位置写入为控制手柄位置的加权组合:
    在这里插入图片描述
  • 其中 ci 是 i 第 th 个控制点的静止位置。这简化了运行时的变形公式。我们可以简单地将形状的每个点的新位置作为平移控制点位置的加权组合:
    在这里插入图片描述
  1. “广义重心坐标”有许多不同的风格(请参阅“自动方法”部分中的表格)。 16 “广义重心坐标”的模糊目标是捕获尽可能多的单纯重心坐标(例如,对于 2D 中的三角形和 3D 中的四面体)对于较大的点或多面体集。一些广义的重心坐标可以以封闭形式计算;其他需要基于优化的预计算。几乎所有的风格都需要连接信息,描述控制点如何围绕输入形状形成外部多面体:笼子。但是,最近的技术不需要笼 23 子。该方法可确保在优化平滑度能量权重期间的仿射精度,其内核中具有仿射函数:
    在这里插入图片描述

  2. 受选定顶点处的插值约束。如果在其内核中具有仿射函数(即 if AV=0 ),则 A 权重 W 将保持仿射精度,我们将得到:
    V=WC

  3. 上述相等的矩阵形式。建议的定义 A 方法是构造一个矩阵,该矩阵 K 在所有内部顶点和所有边界顶点处测量拉普拉斯量。离散拉普拉斯的通常定义(例如,libigl 从 igl::cotmatrix 中返回什么)测量内部顶点函数的拉普拉斯量,但测量函数的拉普拉斯量减去边界顶点函数的法线导数。因此,我们可以让:
    K=L+N

  • 其中 L 是通常的拉普拉斯矩阵, N 是计算网格边界顶点处分段线性函数的法线导数的矩阵。然后 A 以二次形式计算应用于函数的 K 积分平均值的平方并在网格上积分:
    在这里插入图片描述
  1. 由于拉普拉斯算 K 数是二阶导数,它在仿射函数上测量零,因此 A 在其零空间中具有仿射函数。一个简短的推导证明,这意味着 W 将是仿射精确的。
  2. 这种“平方拉普拉斯”能量的最小化器在某种意义上是离散的双谐波函数。因此,它们被称为“双谐波坐标”(与有界双调和权重不同,后者不是广义的重心坐标)。
  3. 在 libigl 中,可以计算给定网格 (V,F) 和选定控制点或控制区域列表 S (类似于蒙皮手柄)的双谐波坐标:
igl::biharmonic_coordinates(V,F,S,W);
  1. 407案例完整展示:
#include <igl/opengl/gl.h>
#include <igl/arap.h>
#include <igl/biharmonic_coordinates.h>
#include <igl/cat.h>
#include <igl/cotmatrix.h>
#include <igl/massmatrix.h>
#include <igl/matrix_to_list.h>
#include <igl/parula.h>
#include <igl/point_mesh_squared_distance.h>
#include <igl/readDMAT.h>
#include <igl/readMESH.h>
#include <igl/remove_unreferenced.h>
#include <igl/slice.h>
#include <igl/writeDMAT.h>
#include <igl/opengl/glfw/Viewer.h>
#include <Eigen/Sparse>
#include <iostream>
#include <queue>


struct Mesh
{
  Eigen::MatrixXd V,U;
  Eigen::MatrixXi T,F;
} low,high,scene;

Eigen::MatrixXd W;
igl::ARAPData arap_data;

int main(int argc, char * argv[])
{
  using namespace Eigen;
  using namespace std;
  using namespace igl;
  
  // read the mesh, if the code is prepared outside of tutorial, the TUTORIAL_SHARED_PATH
  // should be the data folder
  if(!readMESH(TUTORIAL_SHARED_PATH "/octopus-low.mesh",low.V,low.T,low.F))
  {
    cout<<"failed to load mesh"<<endl;
  }
  if(!readMESH(TUTORIAL_SHARED_PATH "/octopus-high.mesh",high.V,high.T,high.F))
  {
    cout<<"failed to load mesh"<<endl;
  }

  // Precomputation
  {
    Eigen::VectorXi b;
    {
      // this will create a vector from 0 to V.rows()-1 where the gap is 1
      Eigen::VectorXi J = Eigen::VectorXi::LinSpaced(high.V.rows(),0,high.V.rows()-1);
      Eigen::VectorXd sqrD;
      Eigen::MatrixXd _2;
      cout<<"Finding closest points..."<<endl;
      // using J which is N by 1 instead of a matrix that represents faces of N by 3
      // so that we will find the closest vertices istead of closest point on the face
      // so far the two meshes are not seperated. So what we are really doing here
      // is computing handles from low resolution and use that for the high resolution one
      igl::point_mesh_squared_distance(low.V,high.V,J,sqrD,b,_2);
      assert(sqrD.minCoeff() < 1e-7 && "low.V should exist in high.V");
    }
    // force perfect positioning, rather have popping in low-res than high-res.
    // The correct/elaborate thing to do is express original low.V in terms of
    // linear interpolation (or extrapolation) via elements in (high.V,high.F)

    // this is to replace the vertices on low resolution
    // with the vertices in high resolution. b is the list of vertices
    // corresponding to the indices in high resolution which has closest
    // distance to the points in low resolution
    igl::slice(high.V,b,1,low.V);


    // list of points --> list of singleton lists
    std::vector<std::vector<int> > S;
    // S will hav size of low.V.rows() and each list inside will have 1 element
    igl::matrix_to_list(b,S);
    cout<<"Computing weights for "<<b.size()<<
      " handles at "<<high.V.rows()<<" vertices..."<<endl;
    // Technically k should equal 3 for smooth interpolation in 3d, but 2 is
    // faster and looks OK
    const int k = 2;

    // using all the points in low resolution as handles for the region
    // it will be too expansive to use all the points in high reolution as handles
    // but since low and high resembles the same thing, using points in low reesolution
    // will give you similar performance
    igl::biharmonic_coordinates(high.V,high.T,S,k,W);
    cout<<"Reindexing..."<<endl;
    // Throw away interior tet-vertices, keep weights and indices of boundary
    VectorXi I,J;
    igl::remove_unreferenced(high.V.rows(),high.F,I,J);
    for_each(high.F.data(),high.F.data()+high.F.size(),[&I](int & a){a=I(a);});
    for_each(b.data(),b.data()+b.size(),[&I](int & a){a=I(a);});
    igl::slice(MatrixXd(high.V),J,1,high.V);
    igl::slice(MatrixXd(W),J,1,W);
  }

  // Resize low res (high res will also be resized by affine precision of W)
  low.V.rowwise() -= low.V.colwise().mean();
  low.V /= (low.V.maxCoeff()-low.V.minCoeff());
  low.V.rowwise() += RowVector3d(0,1,0);
  low.U = low.V;
  high.U = high.V;

  arap_data.with_dynamics = true;
  arap_data.max_iter = 10;
  arap_data.energy = ARAP_ENERGY_TYPE_DEFAULT;
  arap_data.h = 0.01;
  arap_data.ym = 0.001;
  if(!arap_precomputation(low.V,low.T,3,VectorXi(),arap_data))
  {
    cerr<<"arap_precomputation failed."<<endl;
    return EXIT_FAILURE;
  }
  // Constant gravitational force
  Eigen::SparseMatrix<double> M;
  igl::massmatrix(low.V,low.T,igl::MASSMATRIX_TYPE_DEFAULT,M);
  const size_t n = low.V.rows();
  // f = ma
  arap_data.f_ext =  M * RowVector3d(0,-9.8,0).replicate(n,1);
  // Random initial velocities to wiggle things
  arap_data.vel = MatrixXd::Random(n,3);
  
  igl::opengl::glfw::Viewer viewer;
  // Create one huge mesh containing both meshes
  igl::cat(1,low.U,high.U,scene.U);
  // need to remap the indices since we cat the V matrices
  igl::cat(1,low.F,MatrixXi(high.F.array()+low.V.rows()),scene.F);
  // Color each mesh
  viewer.data().set_mesh(scene.U,scene.F);
  MatrixXd C(scene.F.rows(),3);
  C<<
    RowVector3d(0.8,0.5,0.2).replicate(low.F.rows(),1),
    RowVector3d(0.3,0.4,1.0).replicate(high.F.rows(),1);
  viewer.data().set_colors(C);

  viewer.callback_key_pressed = 
    [&](igl::opengl::glfw::Viewer & viewer,unsigned int key,int mods)->bool
  {
    switch(key)
    {
      default: 
        return false;
      case ' ':
        viewer.core().is_animating = !viewer.core().is_animating;
        return true;
      case 'r':
        low.U = low.V;
        return true;
    }
  };
  viewer.callback_pre_draw = [&](igl::opengl::glfw::Viewer & viewer)->bool
  {
    glEnable(GL_CULL_FACE);
    if(viewer.core().is_animating)
    {
      arap_solve(MatrixXd(0,3),arap_data,low.U);
      for(int v = 0;v<low.U.rows();v++)
      {
        // collide with y=0 plane
        const int y = 1;
        if(low.U(v,y) < 0)
        {
          low.U(v,y) = -low.U(v,y);
          // ~ coefficient of restitution
          const double cr = 1.1;
          arap_data.vel(v,y) = - arap_data.vel(v,y) / cr;
        }
      }

      scene.U.block(0,0,low.U.rows(),low.U.cols()) = low.U;
      high.U = W * (low.U.rowwise() + RowVector3d(1,0,0));
      scene.U.block(low.U.rows(),0,high.U.rows(),high.U.cols()) = high.U;

      viewer.data().set_vertices(scene.U);
      viewer.data().compute_normals();
    }
    return false;
  };
  viewer.data().show_lines = false;
  viewer.core().is_animating = true;
  viewer.core().animation_max_fps = 30.;
  viewer.data().set_face_based(true);
  cout<<R"(
[space] to toggle animation
'r'     to reset positions 
      )";
  viewer.core().rotation_type = 
    igl::opengl::ViewerCore::ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP;
  viewer.launch();
}

在这里插入图片描述例407)显示了在粗橙色网格上的物理模拟。该网格的顶点成为蓝色高分辨率网格的双谐坐标变形的控制点

4.8 直接三角洲糊剂 Direct Delta Mush

  1. 为了产生平滑的变形,线性混合蒙皮需要光滑的蒙皮重量。这些可以手动绘制或自动计算(例如,使用有界双谐波权重)。即便如此,线性混合蒙皮由于其固有的线性而遭受收缩和塌陷伪影的影响(见前文)。“直接 Delta Mush” 蒙皮试图通过提供一种直接蒙皮方法来解决这两个问题,该方法将具有分段恒定权重函数(权重在 =0 =1 任意位置或无处不在)的钻机作为输入。直接三角洲糊剂是对一种性能较差的方法的改编,称为“三角洲糊剂”。 Delta Mush 的计算分为“绑定姿势”预计算和运行时评估。

  2. 在绑定时,对绑定姿势进行拉普拉斯平滑,将每个顶点从其静止位置移动到新位置 vi ~vi 。描述撤消此平滑过程的“增量”被计算并存储在与折点关联的局部坐标系中:
    在这里插入图片描述

  3. 在运行时,使用线性混合蒙皮和分段恒定权重对网格进行变形。在骨骼附近,变形是完全刚性的,而在骨骼相遇的关节附近,网格会随着下一个刚性转换的突然变化而撕裂。在运行时将相同数量的拉普拉斯平滑应用于此摆放的网格。将每个顶点移动到一个位置 ~ui 。在此位置计算本地帧,并在此解析帧 Si 中添加缓存的增量以还原形状的原始详细信息:
    在这里插入图片描述

  4. “Delta Mush”的关键见解是,拉普拉斯平滑对其余模型和摆姿势模型的作用相似。

  5. “直接增量糊”的关键见解是,这种运行时拉普拉斯平滑的过程几乎是线性的,并且可以使用SVD以令人尴尬的并行方式计算局部帧(参见ARAP)。

  6. 直接增量糊状将平滑步骤移动到预计算中,从而产生每个顶点每个骨骼的“矢量值”蒙皮权重,存储在矩阵 Ω 中。在 libigl 中,对于网格 (V,F) 和(例如,分段常数)权重, W 此预计算使用以下公式计算:

igl::direct_delta_mush_precomputation(V, F,Wsparse, p, lambda, kappa, alpha, Omega);
  1. 参数 p, lambda, kappa, alpha 控制由此产生的变形的平滑度和紧凑度。预计算的输出是矩阵 Omega 。
  2. 在运行时, Ω 用于将网格变形到其最终位置。在 libigl 中,这是使用以下公式计算的:
igl::direct_delta_mush(V, T_list, Omega, U);
  • 其中 T_list 是与每个骨骼关联的输入姿势(仿射)变换,最终位置存储在 U 中。
  1. 408案例完整展示:
#include <igl/read_triangle_mesh.h>
#include <igl/readTGF.h>
#include <igl/readDMAT.h>
#include <igl/lbs_matrix.h>
#include <igl/deform_skeleton.h>
#include <igl/direct_delta_mush.h>
#include <igl/opengl/glfw/Viewer.h>
#include <Eigen/Geometry>
#include <vector>

int main(int argc, char * argv[]) 
{

  Eigen::MatrixXd V,U,C,W,T,M,Omega;
  Eigen::MatrixXi F,BE;
  igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/elephant.obj",V,F);
  igl::readTGF(           TUTORIAL_SHARED_PATH "/elephant.tgf",C,BE);
  igl::readDMAT(          TUTORIAL_SHARED_PATH "/elephant-weights.dmat",W);
  igl::readDMAT(          TUTORIAL_SHARED_PATH "/elephant-anim.dmat",T);
  // convert weight to piecewise-rigid weights to stress test DDM
  for (int i = 0; i < W.rows(); ++i)
  {
    int maxj;
    W.row(i).maxCoeff(&maxj);
    for (int j = 0; j < W.cols(); j++)
    {
      W(i, j) = double(maxj == j);
    }
  }

  igl::lbs_matrix(V,W,M);

  int p = 20;
  float lambda = 3; // 0 < lambda
  float kappa = 1; // 0 < kappa < lambda
  float alpha = 0.8; // 0 <= alpha < 1
  igl::direct_delta_mush_precomputation(V, F,W, p, lambda, kappa, alpha, Omega);

  igl::opengl::glfw::Viewer viewer;
  int frame = 0;
  const int pr_id = viewer.selected_data_index;
  viewer.append_mesh();
  const int ddm_id = viewer.selected_data_index;
  Eigen::RowVector3d offset(1.1*(V.col(0).maxCoeff()-V.col(0).minCoeff()),0,0);

  viewer.callback_pre_draw = [&](igl::opengl::glfw::Viewer &) -> bool
  {
    if(viewer.core().is_animating)
    {
      const Eigen::Map<Eigen::MatrixXd> Tf(
        T.data()+frame*T.rows(),4*W.cols(),3);
      U = (M*Tf).rowwise()-offset;

      Eigen::MatrixXd Cf;
      Eigen::MatrixXi BEf;
      igl::deform_skeleton(C,BE,Tf,Cf,BEf);
      viewer.data(pr_id).set_edges(Cf,BEf, Eigen::RowVector3d(1,1,1));

      viewer.data(pr_id).set_vertices(U);
      viewer.data(pr_id).compute_normals();

      {
        std::vector<Eigen::Affine3d, Eigen::aligned_allocator<Eigen::Affine3d>> 
          T_list(BE.rows());
        for (int e = 0; e < BE.rows(); e++)
        {
          T_list[e] = Eigen::Affine3d::Identity();
          T_list[e].matrix().block(0,0,3,4) = Tf.block(e*4,0,4,3).transpose();
        }
        igl::direct_delta_mush(V, T_list, Omega, U);
      }
      U.rowwise() += offset;
      viewer.data(ddm_id).set_vertices(U);
      viewer.data(ddm_id).compute_normals();

      frame++;
      if(frame == T.cols())
      {
        frame = 0;
        viewer.core().is_animating = false;
      }
    }
    return false;
  };
  viewer.callback_key_pressed = [&](igl::opengl::glfw::Viewer &, unsigned int key, int mod)
  {
    switch(key)
    {
    case ' ':
      viewer.core().is_animating = !viewer.core().is_animating;
      return true;
    }
    return false;
  };



  for(auto & id : {pr_id,ddm_id})
  {
    if(id == pr_id)
    {
      viewer.data(id).set_mesh( (V.rowwise()-offset*1.0).eval() ,F);
      viewer.data(id).set_colors(Eigen::RowVector3d(214./255.,170./255.,148./255.));
      viewer.data(id).set_edges(C,BE, Eigen::RowVector3d(1,1,1));
    }else if(id == ddm_id){
      viewer.data(id).set_mesh( (V.rowwise()+offset*1.0).eval() ,F);
      viewer.data(id).set_colors(Eigen::RowVector3d(132./255.,214./255.,105./255.));
    }
    viewer.data(id).show_lines = false;
    viewer.data(id).set_face_based(true);
    viewer.data(id).show_overlay_depth = false;
  }
  viewer.core().is_animating = false;
  viewer.core().animation_max_fps = 24.;
  //viewer.data().set_colors(V,F);


  viewer.launch_init();
  viewer.core().align_camera_center(V);

  viewer.launch_rendering(true);
  viewer.launch_shut();
}

在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述例408)直接三角洲糊状物。(左)输入分段刚性蒙皮,(中)骨架动画,(右)平滑的直接三角洲泥塑蒙皮(动态过程)

4.9 使用开尔文的网格变形 Mesh Deformation with Kelvinlet

开尔文莱茨是一种基于物理的虚拟弹性材料实时体积雕刻技术。该技术将网格视为由可压缩材料制成的流体,并通过沿位移场的对流点使它们变形。它依赖于弹性方程的解析解。

4.9.1 线性弹性静力学快速入门 A quick primer on linear elastostatics

  1. 线性弹性的平衡状态由位移场决定,该位移场 u:R^ 3→R^3 使弹性势能最小化
    在这里插入图片描述
  • 其中 μ 是弹性剪切模量,是泊松比, ν b 代表外体力。
  1. 第一项控制位移场的平滑度,第二项惩罚无穷小的体积变化,最后一项表示要抵消的外部体力。
  2. 可以将最优位移场与上述方程(也称为纳维-柯西方程)的临界点的解相关联:
    在这里插入图片描述
  3. 开尔文雷是纳维-柯西方程的解,在某一点 x0 的力矢量 f 导致身体负荷集中的情况下,即 b(x)=fδ(x−x0) 和 可以写为:
    在这里插入图片描述在这里插入图片描述在这里插入图片描述
  4. 此梯度 G® 决定了位移场 u® 的不同属性。例如,的 ∇u® 偏斜对称部分表示由 引起的 u® 旋转,而其对称部分对应于弹性应变并决定拉伸。应变张量也可以分解为表示弹性介质体积缩放的迹线项和表示夹紧变形的无迹项。
  5. 这构成了开尔文雷画笔的基础。

4.9.2 正则化开尔文 Regularized kelvinlets

  1. 单个点处的集中体负荷在 处 x0 引入了开尔文雷解的奇点 x0 。因此,开尔文雷方程修改为:
    在这里插入图片描述在这里插入图片描述

在这里插入图片描述

  1. 不同径向尺度的开尔文可以线性组合,构建具有任意快速衰减的画笔:
    在这里插入图片描述4.

在这里插入图片描述

  1. 正则化开尔文可以进一步扩展,方法是将基于矢量的载荷分布替换为基于矩阵的分布,以实现非仿射变换,如前面描述的扭曲、捏合和缩放。
  2. 在 libigl 中,这是使用以下公式计算的:
igl::KelvinletParams<double> brushParams{brushRadius, scale, brushType};
igl::kelvinlets(V, origin, forceVec, forceMatrix, brushParams, result);
  • 其中 brushRadius, scale, brushType 对应于 ϵ 、衰减和操作(抓取、捏合、缩放、扭曲)。
  1. 409案例完整展示:
#include <igl/kelvinlets.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/opengl/glfw/imgui/ImGuiPlugin.h>
#include <igl/opengl/glfw/imgui/ImGuiMenu.h>
#include <igl/readOFF.h>
#include <igl/unproject.h>
#include <igl/unproject_onto_mesh.h>
#include <imgui.h>
#include <iostream>

namespace {
void ShowHelpMarker(const char* desc)
{
  ImGui::SameLine();
  ImGui::TextDisabled("(?)");
  if (ImGui::IsItemHovered()) {
    ImGui::BeginTooltip();
    ImGui::PushTextWrapPos(450.0f);
    ImGui::TextUnformatted(desc);
    ImGui::PopTextWrapPos();
    ImGui::EndTooltip();
  }
}
}

int main()
{
  Eigen::MatrixXd V1, OrigV;
  Eigen::MatrixXi F1, OrigF;

  igl::readOFF(TUTORIAL_SHARED_PATH "/bumpy.off", OrigV, OrigF);
  std::cout << "1 View original mesh\n";
  std::cout << "2 Switch to deformed mesh\n";
  V1 = OrigV;
  F1 = OrigF;

  igl::opengl::glfw::Viewer viewer;
  igl::opengl::glfw::imgui::ImGuiPlugin plugin;
  viewer.plugins.push_back(&plugin);
  igl::opengl::glfw::imgui::ImGuiMenu menu;
  plugin.widgets.push_back(&menu);

  auto brushRadius = 1.;
  auto brushType = igl::BrushType::GRAB;
  auto scale = 1;
  menu.callback_draw_custom_window = [&]() {
    ImGui::SetNextWindowPos(ImVec2(180.f * menu.menu_scaling(), 10),
                            ImGuiCond_FirstUseEver);
    ImGui::SetNextWindowSize(ImVec2(200, 160), ImGuiCond_FirstUseEver);
    ImGui::Begin(
      "Kelvinlet Brushes", nullptr, ImGuiWindowFlags_NoSavedSettings);

    ImGui::InputDouble("Brush Radius", &brushRadius, 0, 0, "%.4f");
    ImGui::Combo("Brush type",
                 reinterpret_cast<int*>(&brushType),
                 "Grab\0Scale\0Twist\0Pinch\0\0");
    ImGui::InputInt("Falloff", &scale);
    ShowHelpMarker("Defines how localized the stroke is {1,2,3}");
    ImGui::End();
  };

  Eigen::Vector3d posStart(0, 0, 0);
  Eigen::Vector3d posEnd;
  decltype(OrigV) result;
  auto min_point = V1.colwise().minCoeff();
  auto max_point = V1.colwise().maxCoeff();
  // to multiply brush force proportional to size of mesh
  auto brush_strength = (max_point - min_point).norm();
  Eigen::Matrix3d twist, pinch;
  twist << 0, 1, -1, -1, 0, 1, 1, -1, 0; // skew-symmetric
  pinch << 0, 1, 1, 1, 0, 1, 1, 1, 0;    // symmetric

  viewer.callback_key_down =
    [&](igl::opengl::glfw::Viewer& viewer, unsigned char key, int) {
      if (key == '1') {
        viewer.data().clear();
        viewer.data().set_mesh(OrigV, OrigF);
        viewer.core().align_camera_center(OrigV, OrigF);
      } else if (key == '2') {
        viewer.data().clear();
        viewer.data().set_mesh(V1, F1);
        viewer.core().align_camera_center(V1, F1);
      }
      return false;
    };

  viewer.callback_mouse_down =
    [&](igl::opengl::glfw::Viewer& viewer, int, int) -> bool {
    Eigen::Vector3f bc;
    int fid;
    auto x = viewer.current_mouse_x;
    auto y =
      viewer.core().viewport(3) - static_cast<float>(viewer.current_mouse_y);
    if (igl::unproject_onto_mesh(Eigen::Vector2f(x, y),
                                 viewer.core().view,
                                 viewer.core().proj,
                                 viewer.core().viewport,
                                 V1,
                                 F1,
                                 fid,
                                 bc)) {
      posStart = igl::unproject(Eigen::Vector3f(x, y, viewer.down_mouse_z),
                                viewer.core().view,
                                viewer.core().proj,
                                viewer.core().viewport)
                   .template cast<double>();
      return true;
    }
    return false;
  };

  viewer.callback_mouse_move =
    [&](igl::opengl::glfw::Viewer& viewer, int, int) -> bool {
    if (!posStart.isZero() && !posStart.hasNaN()) {
      posEnd = igl::unproject(
                 Eigen::Vector3f(viewer.current_mouse_x,
                                 viewer.core().viewport[3] -
                                   static_cast<float>(viewer.current_mouse_y),
                                 viewer.down_mouse_z),
                 viewer.core().view,
                 viewer.core().proj,
                 viewer.core().viewport)
                 .template cast<double>();

      // exaggerate the force by a little bit
      Eigen::Vector3d forceVec = (posEnd - posStart) * brush_strength;

      int scaleFactor = forceVec.norm();
      if (posEnd.x() < posStart.x()) {
        // probably not the best way to determine direction.
        scaleFactor = -scaleFactor;
      }
      Eigen::Matrix3d mat;
      switch (brushType) {
        case igl::BrushType::GRAB:
          mat.setZero();
          break;
        case igl::BrushType::SCALE:
          mat = Eigen::Matrix3d::Identity() * scaleFactor;
          break;
        case igl::BrushType::TWIST:
          mat = twist * scaleFactor;
          break;
        case igl::BrushType::PINCH:
          mat = pinch * scaleFactor;
          break;
      }

      igl::kelvinlets(
        V1,
        posStart,
        forceVec,
        mat,
        igl::KelvinletParams<double>(brushRadius, scale, brushType),
        result);
      viewer.data().set_vertices(result);
      viewer.data().compute_normals();
      return true;
    }
    return false;
  };

  viewer.callback_mouse_up =
    [&](igl::opengl::glfw::Viewer& viewer, int, int) -> bool {
    if (!posStart.isZero()) {
      V1 = result;
      posStart.setZero();
      return true;
    }
    return false;
  };

  viewer.data().set_mesh(V1, F1);
  viewer.core().align_camera_center(V1, F1);
  viewer.launch();
}

在这里插入图片描述(例409) pinch, twist, grab and scale in action 捏、捏、捏、抓、缩放

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值