实现算法：简易的Laplace mesh editing。（2016年3月5日）
$实现过程及原理
Laplace 坐标（算子）
。
这里的表示相对坐标，是某一点相对于其1-领域点的中心点的位置，也就是Laplace坐标。为什么这样表示呢？这样的坐标有两个特点：
1.它的方向与法相比较接近。
2.它的大小接近于曲率。
对于一些网格，我们进行拖动一些点，但是又要使得它的形状和原先的尽可能相似，怎么处理呢？假设为已知大概位置的点，其变换后的大概位置是。考虑这样的能量函数：
先来看看这个函数的意义，第一部分表示Laplace坐标变化大小，也就是原先的形状变化大小，第二部分表示对于固定点的位置变化大小。我们最小化这个能像函数就能得到一组点,既能使得形状变化小，又能使得与固定位置变化小。
$实验步骤
1.对能量函数进行求导，求其最小值相当于计算线性方程组
，
，
2.编写程序，测试。（参见附录）
3.测试结果及软件使用方法
（1）选择一些移动点
（2）计算一些必要的参数
（3）改变选中点的位置，这几个按钮都可以
（4）点击进行编辑
（5）效果如下
（6）对比原图。
$实验注意事项：
1.在对能能量函数进行求解事需要细心，转化成最小二乘法求解线性方程组时每一步都得耐心检查，否则一失足而成千古恨，以至于后面全都错。
2.在对方程组进行输入时，每进行一步都得测试，以防带入变量或者参数错误而得不到正确的结果
3.整合每一块的计算结果，需要对程序的操作流程非常熟悉，得知道第一步该干什么第二步该干什么，如果两个函数的操作顺序反了可能得不到想要的结果，比如如果先对点进行移动在计算Laplace坐标时就是不对的。
4.该清除的量必须清除，否则上一次变换的结果可能对下一次变换产生影响，结果导致只有第一次变换是对的，之后的变换全是错的。
总而言之，需要对实验的基本原理和程序的基本流程非常清楚才能进行编程，否则只会浪费更多的时间。
$主要程序代码:
void RenderingWidget::edit_()
{
classfly();
compute_matrix();
get_b();
solve_system();
updateGL();
triple_A.clear();
}
void RenderingWidget::compute_parameter()
{
ptr_mesh_->laplace_position();//先计算laplace坐标
}
////////////分类/////////////
void RenderingWidget::classfly()
{
int k =0;
const std::vector<HE_vert*>& vertices = *(ptr_mesh_->get_vertex_list());
for(int i = 0;i<vertices.size();i++)
{
vertices[i]->point_type = other;
}
for (int i = 0; i<vertices.size(); i++)
{
if (vertices[i]->boundary_flag_ == BOUNDARY) vertices[i]->point_type = boundary;
}
for (int i = 0; i < current_index_.size(); i++)
{
vertices[current_index_[i]]->point_type = drag;
}
}
/////////////计算各种矩阵//////////////////
void RenderingWidget::compute_matrix()
{
const std::vector<HE_vert*>& vertices = *(ptr_mesh_->get_vertex_list());
float r1 = 1.0f / vertices.size(); //所有点系数
float r2 = 1.0f / current_index_.size(); //拖动点系数
float r3 = 0; //边界点系数
for (int i = 0; i < vertices.size(); i++) //没错
{
if (vertices[i]->boundary_flag_==BOUNDARY)
{
r3++;
}
}
if (r3>0)
{
r3 = 1.0f / r3;
}
int n = vertices.size();//点的个数
float lambda = 0.01;
float beta = 0.01;
for (int i = 0; i < n; i++)
{
triple_A.push_back(Triplet<float>(i, i, 1.0f*r1));
if (vertices[i]->point_type==boundary)//第二部分，边界点
{
triple_A.push_back(Triplet<float>(i, i, r3*beta));
}
if (vertices[i]->point_type == drag)
{
triple_A.push_back(Triplet<float>(i, i, r2*lambda));
}
for (int j = 0; j < vertices[i]->neighborIdx.size(); j++)
{
triple_A.push_back(Triplet<float>(i, vertices[i]->neighborIdx[j], -1.0f / vertices[i]->degree_*r1));
}
for (int j = 0; j < vertices[i]->degree_;j++)
{
int idx = vertices[i]->neighborIdx[j];
int dk = vertices[vertices[i]->neighborIdx[j]]->degree_;
triple_A.push_back(Triplet<float>(i, idx, -1.0f / dk*r1));
for (int k = 0; k < dk; k++)
{
triple_A.push_back(Triplet<float>(i, vertices[idx]->neighborIdx[k], 1.0f / dk / dk*r1));
}
}
}
SparseMatrix<float> A(n, n);
A.setFromTriplets(triple_A.begin(), triple_A.end());
solver.compute(A);//对A进行预分解
}
void RenderingWidget::get_b()
{
const std::vector<HE_vert*>& vertices = *(ptr_mesh_->get_vertex_list());
float r1 = 1.0f / vertices.size(); //没错
float r2 = 1.0f / current_index_.size(); //没错
float r3 = 0; //没错
for (int i = 0; i < vertices.size(); i++) //没错
{
if (vertices[i]->boundary_flag_==BOUNDARY)
{
r3++;
}
}
if (r3>0)
{
r3 = 1.0f / r3;
}
int n = vertices.size();
float lambda = 0.01;
float beta = 0.01;
VectorXf bx_(n), by_(n), bz_(n);
bx_.fill(0);
by_.fill(0);
bz_.fill(0);
for (int i = 0; i < n; i++)
{
if (vertices[i]->point_type == drag)//拖动点
{
bx_[i] += vertices[i]->position_[0] * r2*lambda;
by_[i] += vertices[i]->position_[1] * r2*lambda;
bz_[i] += vertices[i]->position_[2] * r2*lambda;
}
if (vertices[i]->boundary_flag_==BOUNDARY)//第二部分，边界点
{
bx_[i] += vertices[i]->position_[0] * r3*beta;
by_[i] += vertices[i]->position_[1] * r3*beta;
bz_[i] += vertices[i]->position_[2] * r3*beta;
}
bx_[i] += vertices[i]->ab_position[0] * r1;
by_[i] += vertices[i]->ab_position[1] * r1;
bz_[i] += vertices[i]->ab_position[2] * r1;
for (int j = 0; j < vertices[i]->degree_; j++)
{
int idx = vertices[i]->neighborIdx[j];
int dk = vertices[vertices[i]->neighborIdx[j]]->degree_;
bx_[i] -= vertices[idx]->ab_position[0] / dk*r1;
by_[i] -= vertices[idx]->ab_position[1] / dk*r1;
bz_[i] -= vertices[idx]->ab_position[2] / dk*r1;
}
}
SparseMatrix<float> A(n, n);
A.setFromTriplets(triple_A.begin(), triple_A.end());
bx.clear();
by.clear();
bz.clear();
for (int i = 0; i < n; i++)
{
bx.push_back(bx_(i));
by.push_back(by_(i));
bz.push_back(bz_(i));
}
}
void RenderingWidget::solve_system()
{
const vector<HE_vert*>& vertices = *(ptr_mesh_->get_vertex_list());
int n = vertices.size();
VectorXf bx_(n), by_(n), bz_(n);
bx_.fill(0);
by_.fill(0);
bz_.fill(0);
for (int i = 0; i < n; i++)
{
bx_[i] = bx[i];
by_[i] = by[i];
bz_[i] = bz[i];
}
bx_ = solver.solve(bx_);
by_ = solver.solve(by_);
bz_ = solver.solve(bz_);
for (int i = 0; i < vertices.size(); i++)
{
vertices[i]->position_[0] = bx_(i);
vertices[i]->position_[1] = by_(i);
vertices[i]->position_[2] = bz_(i);
}
}