两个网格的近似误差估计

之前也做了网格简化的东西：  三维网格精简算法（Quadric Error Metrics）附源码

在同样一篇论文中提出了估计误差的一个公式，非常类似于Hoppe等在Piecewise Smooth Surface Reconstruction中提出的E_dist能量，我们定义近似误差E_i=E(M_n,M_i)为

其中X_n和X_i是分别对应于模型M_n和M_i的点集，距离d(v,M)是从v到M中最近面的距离（欧式距离）。就用E来度量近似误差，下面就是具体实现的程序代码啦

function [ distances, surface_points] = point2mesh_all( faces,vertices,qPoints )

assert(~isempty(faces) && ~isempty(vertices), 'Invalid argument: ''Faces'' and ''Vertices'' mustn''t be empty.')
assert(max(faces(:))<=size(vertices,1), 'The value of ''Faces'' is invalid: the maximum vertex ID is bigger than the number of vertices in ''Vertices''')

% Calculate normals
r1 = vertices(faces(:,1),:);  % (#faces x 3) % 1st vertex of every face
r2 = vertices(faces(:,2),:);  % (#faces x 3) % 2nd vertex of every face
r3 = vertices(faces(:,3),:);  % (#faces x 3) % 3rd vertex of every face
normals = cross((r2-r1),(r3-r1),2); % (#faces x 3) normal vector of every face
normals = bsxfun(@rdivide,normals,sqrt(sum(normals.^2,2))); % (#faces x 3) normalized normal vector of every face

if isempty(qPoints)
    distances = [];
    surface_points = [];
    return
end


%% Distance Calculation

nQPoints = size(qPoints,1);

D = NaN(nQPoints,1);
P = NaN(nQPoints,3);


 %case {'linear','normal'}
 for r = 1:nQPoints
     % Determine the surface points
     point = qPoints(r,:); % (1 x 3) query point
     [d,p] = processPoint(faces,vertices,point,normals,@distance_to_vertices,@distance_to_edges,@distance_to_surfaces);
     D(r) = d;
     P(r,:) = p;

 end


% return output arguments
distances      = D;  % (#qPoints x 1)
surface_points = P;  % (#qPoints x 3)


end

function [D,P] = processPoint(faces,vertices,point,normals,distance_to_vertices,distance_to_edges,distance_to_surfaces)

d = NaN(3,1); % (distanceTypes x 1)
p = NaN(3,3); % (distanceTypes x xyz)


[d(1),p(1,:)] = distance_to_vertices(faces,vertices,point,normals);
% find nearest vertice
[d(2),p(2,:)] = distance_to_edges(faces,vertices,point,normals);
[d(3),p(3,:)] = distance_to_surfaces(faces,vertices,point,normals);

% find nearest point on all edges


% find minimum distance type
[~,I] = min(abs(d),[],1);
D = d(I);
P = p(I,:);
end


%% Non-vectorized Distance Functions


function [D,P] = distance_to_vertices(faces,vertices,qPoint,normals)

% find nearest vertex
[D,nearestVertexID] = min(sum(bsxfun(@minus,vertices,qPoint).^2,2),[],1);
D = sqrt(D);
P = vertices(nearestVertexID,:); % (1 x 3)

% find faces that belong to the vertex
connectedFaces = find(any(faces==nearestVertexID,2)); % (#connectedFaces x 1) face indices
assert(length(connectedFaces)>=1,'Vertex %u is not connected to any face.',nearestVertexID)


n = normals(connectedFaces,:); % (#connectedFaces x 3) normal vectors

% scalar product between distance vector and normal vectors
coefficients = dot2(n,qPoint-P);
sgn = signOfLargest(coefficients);
D = D*sgn;

end


function [D,P] = distance_to_edges(faces,vertices,qPoint,normals)

% Point-point representation of all edges
edges = [faces(:,[1,2]); faces(:,[1,3]); faces(:,[2,3])]; % (#edges x 2) vertice IDs

% Intersection between tangent of edge lines and query point
r1 = vertices(edges(:,1),:);   % (#edges x 3) first point of every edge
r2 = vertices(edges(:,2),:);   % (#edges x 3) second point of every edge
t = dot( bsxfun(@minus,qPoint,r1), r2-r1, 2) ./ sum((r2-r1).^2,2); % (#edges x 1) location of intersection relative to r1 and r2
t(t<=0) = NaN; % exclude intersections not between the two vertices r1 and r2
t(t>=1) = NaN;

% Distance between intersection and query point
P = r1 + bsxfun(@times,(r2-r1),t); % (#edges x 3) intersection points
D = bsxfun(@minus,qPoint,P); % (#edges x 3)
D = sqrt(sum(D.^2,2));       % (#edges x 1)
[D,I] = min(D,[],1);         % (1 x 1)
P = P(I,:);

% find faces that belong to the edge
inds = edges(I,:);  % (1 x 2)
inds = permute(inds,[3,1,2]);  % (1 x 1 x 2)
inds = bsxfun(@eq,faces,inds); % (#faces x 3 x 2)
inds = any(inds,3);    % (#faces x 3)
inds = sum(inds,2)==2; % (#faces x 1) logical indices which faces belong to the nearest edge of the query point
n = normals(inds,:); % (#connectedFaces x 3) normal vectors

% scalar product between distance vector and normal vectors
coefficients = dot2(n,qPoint-P); % (#connectedFaces x 1)
sgn = signOfLargest(coefficients);
D = D*sgn;

end



function [D,P] = distance_to_surfaces(faces,vertices,point,normals)

r1 = vertices(faces(:,1),:);   % (#faces x 3) % 1st vertex of every face
r2 = vertices(faces(:,2),:);   % (#faces x 3) % 2nd vertex of every face
r3 = vertices(faces(:,3),:);   % (#faces x 3) % 3rd vertex of every face

vq = bsxfun(@minus,point,r1);  % (#faces x 3)
D = dot(vq,normals,2);         % (#faces x 1) distance to surface
rD = bsxfun(@times,normals,D); % (#faces x 3) vector from surface to query point
P = bsxfun(@minus,point,rD);   % (#faces x 3) nearest point on surface; can be outside triangle

% find barycentric coordinates (query point as linear combination of two edges)
r31r31 = sum((r3-r1).^2,2);    % (#faces x 1)
r21r21 = sum((r2-r1).^2,2);    % (#faces x 1)
r21r31 = dot(r2-r1,r3-r1,2);   % (#faces x 1)
r31vq = dot(r3-r1,vq,2);       % (#faces x 1)
r21vq = dot(r2-r1,vq,2);       % (#faces x 1)

d = r31r31.*r21r21 - r21r31.^2;               % (#faces x 1)
bary = NaN(size(faces,1),3);                  % (#faces x 3)
bary(:,1) = (r21r21.*r31vq-r21r31.*r21vq)./d; % (#faces x 3)
bary(:,2) = (r31r31.*r21vq-r21r31.*r31vq)./d; % (#faces x 3)
bary(:,3) = 1 - bary(:,1) - bary(:,2);        % (#faces x 3)

% tri = triangulation(faces,vertices);
% bary = tri.cartesianToBarycentric((1:size(faces,1))',P); % (#faces x 3)

% exclude intersections that are outside the triangle
D( abs(d)<=eps | any(bary<=0,2) | any(bary>=1,2) ) = NaN;  % (#faces x 1)

% find nearest face for query point
[~,I] = min(abs(D),[],1); % (1 x 1)
D = D(I);       % (1 x 1)
P = P(I,:);     % (1 x 3)

end



function sgn = signOfLargest(coeff)
    [~,I] = max(abs(coeff));
    sgn = sign(coeff(I));
    if sgn==0, sgn=1; end
end


function d = dot2(A,B)
    % dot product along 2nd dimension with singleton extension
    d = sum(bsxfun(@times,A,B),2);
end

这里的代码并没有想象中那么容易，因为要计算的是到 面片的距离，面片是有边界的！

这个函数由三个函数组成，对于指定的点，计算到网格的最近距离由三个部分组成

如果最近点落在某个面片法向投影的内部，最近距离为法线段的长度。如果落在边界上（线段）

       如果最近点落在某个线段法向投影的内部，最近距离为到这条线段法向线段的长度，如果落在边界上（顶点）

              最近点是到网格某个顶点的距离

输入两个网格，可以计算出近似误差，并且绘制出误差分布

例如输入两个rockerarm的obj模型，运行

obj1='rockerarm_1000.obj';
obj2='rockerarm.obj';

[V1,F1]=obj__read(obj1);V1=V1';F1=F1';
[V2,F2]=obj__read(obj2);V2=V2';F2=F2';
%points=V2;
faces=F1;vertices=V1;qPoints=V2;

%[ distances, surface_points ,FF,Edges] = point2mesh_line_and_face(faces,vertices,qPoints);
[distances, ~ ] = point2mesh_all( faces,vertices,qPoints );

trep=triangulation(F2,V2);
figure(1);colormap jet;trisurf(trep,distances,'edgecolor','interp','FaceColor','interp');caxis([min(distances) max(distances)]);...
axis square;colorbar('vert');brighten(-0.1);axis off;


faces=F2;vertices=V2;qPoints=V1;

[distances1, ~ ] = point2mesh_all( faces,vertices,qPoints );

trep1=triangulation(F1,V1);
figure(2);colormap jet;trisurf(trep1,distances1,'edgecolor','interp','FaceColor','interp');caxis([min(distances1) max(distances1)]);...
axis square;colorbar('vert');brighten(-0.1);axis off;

得到