两个网格的近似误差估计

之前也做了网格简化的东西：  三维网格精简算法（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 normalsr1 = vertices(faces(:,1),:);  % (#faces x 3) % 1st vertex of every facer2 = vertices(faces(:,2),:);  % (#faces x 3) % 2nd vertex of every facer3 = vertices(faces(:,3),:);  % (#faces x 3) % 3rd vertex of every facenormals = cross((r2-r1),(r3-r1),2); % (#faces x 3) normal vector of every facenormals = bsxfun(@rdivide,normals,sqrt(sum(normals.^2,2))); % (#faces x 3) normalized normal vector of every faceif isempty(qPoints)    distances = [];    surface_points = [];    returnend%% Distance CalculationnQPoints = 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 argumentsdistances      = D;  % (#qPoints x 1)surface_points = P;  % (#qPoints x 3)endfunction [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 Functionsfunction [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 vertexconnectedFaces = find(any(faces==nearestVertexID,2)); % (#connectedFaces x 1) face indicesassert(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 vectorscoefficients = dot2(n,qPoint-P);sgn = signOfLargest(coefficients);D = D*sgn;endfunction [D,P] = distance_to_edges(faces,vertices,qPoint,normals)% Point-point representation of all edgesedges = [faces(:,[1,2]); faces(:,[1,3]); faces(:,[2,3])]; % (#edges x 2) vertice IDs% Intersection between tangent of edge lines and query pointr1 = vertices(edges(:,1),:);   % (#edges x 3) first point of every edger2 = vertices(edges(:,2),:);   % (#edges x 3) second point of every edget = dot( bsxfun(@minus,qPoint,r1), r2-r1, 2) ./ sum((r2-r1).^2,2); % (#edges x 1) location of intersection relative to r1 and r2t(t<=0) = NaN; % exclude intersections not between the two vertices r1 and r2t(t>=1) = NaN;% Distance between intersection and query pointP = r1 + bsxfun(@times,(r2-r1),t); % (#edges x 3) intersection pointsD = 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 edgeinds = 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 pointn = normals(inds,:); % (#connectedFaces x 3) normal vectors% scalar product between distance vector and normal vectorscoefficients = dot2(n,qPoint-P); % (#connectedFaces x 1)sgn = signOfLargest(coefficients);D = D*sgn;endfunction [D,P] = distance_to_surfaces(faces,vertices,point,normals)r1 = vertices(faces(:,1),:);   % (#faces x 3) % 1st vertex of every facer2 = vertices(faces(:,2),:);   % (#faces x 3) % 2nd vertex of every facer3 = vertices(faces(:,3),:);   % (#faces x 3) % 3rd vertex of every facevq = bsxfun(@minus,point,r1);  % (#faces x 3)D = dot(vq,normals,2);         % (#faces x 1) distance to surfacerD = bsxfun(@times,normals,D); % (#faces x 3) vector from surface to query pointP = 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 triangleD( 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)endfunction sgn = signOfLargest(coeff)    [~,I] = max(abs(coeff));    sgn = sign(coeff(I));    if sgn==0, sgn=1; endendfunction 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;

得到