最小二乘网格就是利用所谓软约束soft constraint，把已知的Vertex（ 控制点 Vh）的约束作为新的行I*V=Vh，添加到方程组L*V=0后面，对应的是最小二乘的regularization项。

加载模型，V的对称序号，计算Laplacian：

clear all
close all
path = 'data\';
load([path,'Basel_color_meanface.mat']);
meanshape = meanface.shape;
meantex = meanface.tex;
numv = size(meanshape,1);

facetri = 1 + load([path,'face_tri.txt']);

symid = 1 + load([path,'face05_symlist.txt']);
% conert to ids, each v has a symetry point
ids=[symid;fliplr(symid)];ids=unique(ids,'rows');
ids=sortrows(ids,1);ids = ids(:,2);

[L,W] = cal_mesh_laplacian(double(meanshape),facetri,'conformal');

每次选择重建误差最大的点加入到控制点集idcon中，idcon初始值选模型中z最大的点。总数到400退出：

%soft constaint
[~,idcon]=max(meanshape(:,3)); % tip of nose 8321

total_count = 400;
while(1)
    ncont = length(idcon);  
    if(ncont>total_count)       
     break;    
    end
    i=(1:ncont);
    % large weight
    wi = 32 ;
    Lv = sparse(i,idcon,wi,ncont,numv);
    Lhs = [L;Lv];

    vh = sparse(double(meanshape(idcon,:)));
    Rhs = [zeros(numv,3); wi*vh];    
    V = Lhs\Rhs;
    V = full(V);
    dv = (V-meanshape);
    ds = dv(:,1).^2 + dv(:,2).^2 + dv(:,3).^2;
    [dsm,idmin] = max(ds);
    idnew = idmin;   
     
    idcon = [idcon;idnew];  
    if(length(idcon)==10 || length(idcon)==100  || length(idcon)==200)
       % draw with trimesh here 
       % ...
    end
end

加入10个点：

加入100个点：

加入300个点：

加入400个点：

速度比较慢，不过还可以容忍。模型不是封闭的，所以靠近边界处集中了比较多的点。

由于人脸模型沿着X=0左右对称，我们希望选取左右对称的点，这样就要用到Basel库提供的点的对称序号表face05_symlist.txt把它处理成一个新的格式：每个序号V都对应它的对称序号，就是上面程序的ids数组，和numv等长度。这样一个结点v的对称点可以通过vsym=ids(v)找到。逐个添加点的代码可以改成，只要点不在X=0的对称面上就把它的对称点也加进去。上面代码的

idnew = idmin;

 替换成：

idmin_sym = ids(idmin);        
if( idmin == idmin_sym )
  idnew = idmin;        
else
  idnew = [idmin;idmin_sym];        
end

在参数化的纹理空间进行重新三角化

得到的subgraph：

参考文献：

 [1] Olga Sorkine and Daniel Cohen-Or.  Least-Squares Meshes