Delaunay三角剖分
code
function Delaunary(pt)
% 生成三角面片点云
edgeHull=ConvexHull(pt);
figure,plot(pt(:,1),pt(:,2),'b.'),axis equal,hold on
triout=[];
newIndex = -1;
len = 1;len1= 1;
index = 1;
while index <= length(edgeHull)
edge = edgeHull(index);
index = index + 1;
len = length(edgeHull);
len1= length(triout);
dex = [edge.start edge.end];
% 判断左三角形是否存在，存在，则判断右三角行
if edge.lefttri == -1
newIndex = -1; % 新选取点索引
angle12Max = 0; angle23Max = 0;
vector1 = [pt(dex(2),1) - pt(dex(1),1), pt(dex(2),2) - pt(dex(1),2)];
for j = 1 : num
if j ~= dex(1) && j ~= dex(2) % 排除边的端点
vector2 = [pt(j,1) - pt(dex(1),1), pt(j,2) - pt(dex(1),2)];
if vector1(1) * vector2(2) > vector1(2) * vector2(1) % 点j在vector1的左侧
vector3 = [pt(j,1) - pt(dex(2),1), pt(j, 2) - pt(dex(2),2)];
v1Len = norm(vector1); v2Len = norm(vector2); v3Len = norm(vector3);
angle12 = acos(dot(vector1, vector2) / (v1Len * v2Len));
angle23 = acos(dot(vector2, vector3) / (v2Len * v3Len));
if angle23 > angle23Max
angle23Max = angle23;
angle12Max = angle12;
newIndex = j;
elseif angle23 == angle23Max && angle12 <= angle12Max
angle12Max = angle12;
newIndex = j;
end
end
end
end
if newIndex ~= -1 % 若找到一个满足要求的点
tri.a = dex(1);
tri.b = dex(2);
tri.c = newIndex;
tri.indexa = -1; % 第一个领接三角形索引
tri.indexb = -1;
tri.indexc = -1;
edge.lefttri = length(triout); % 设置边左侧三角形索引
isTriExist = 0;
% 记录边1
for k = 1 : length(edgeHull)
tempEdge = edgeHull(k);
if tempEdge.start == dex(1) && tempEdge.end == newIndex
tempEdge.righttri = length(triout);
tri.indexb = tempEdge.lefttri;
isTriExist = 1;
break;
elseif tempEdge.start == newIndex && tempEdge.end == dex(1)
tempEdge.lefttri = length(triout);
tri.indexb = tempEdge.righttri;
isTriExist = 1;
break;
end
end
% 若不存在这条边就新建一条边
if isTriExist == 0
newEdge.start = -1; newEdge.end = -1;
newEdge.lefttri = -1; newEdge.righttri = -1;
newEdge.start = newIndex;
newEdge.end = dex(1);
newEdge.lefttri = length(triout);
edgeHull = [edgeHull newEdge]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%????
end
isTriExist = 0;
% 记录边2
for k = 1 : length(edgeHull)
tempEdge = edgeHull(k);
if tempEdge.start == newIndex && tempEdge.end == dex(2)
tempEdge.righttri = length(triout);
tri.indexa = tempEdge.lefttri;
isTriExist = 1;
break;
elseif tempEdge.start == dex(2) && tempEdge.end == newIndex
tempEdge.lefttri = length(triout);
tri.indexa = tempEdge.righttri;
isTriExist = 1;
break;
end
end
if isTriExist == 0
newEdge.start = -1; newEdge.end = -1;
newEdge.lefttri = -1; newEdge.righttri = -1;
newEdge.start = dex(2);
newEdge.end = newIndex;
newEdge.lefttri = length(triout);
edgeHull = [edgeHull newEdge];
end
tri.indexc = edge.righttri;
triout = [triout tri];
end
elseif edge.righttri == -1
newIndex = -1; % 新选取点索引
angle12Max = 0; angle23Max = 0;
vector1 = [pt(dex(2),1) - pt(dex(1),1), pt(dex(2),2) - pt(dex(1),2)];
for j = 1 : num
if j ~= dex(1) && j ~= dex(2) % 排除边的端点
vector2 = [pt(j,1) - pt(dex(1),1), pt(j,2) - pt(dex(1),2)];
if vector1(1) * vector2(2) < vector1(2) * vector2(1) % 点j在vector1的左侧
vector3 = [pt(j,1) - pt(dex(2),1), pt(j, 2) - pt(dex(2),2)];
v1Len = norm(vector1); v2Len = norm(vector2); v3Len = norm(vector3);
angle12 = acos(dot(vector1, vector2) / (v1Len * v2Len));
angle23 = acos(dot(vector2, vector3) / (v2Len * v3Len));
if angle23 > angle23Max
angle23Max = angle23;
angle12Max = angle12;
newIndex = j;
elseif angle23 == angle23Max && angle12 >= angle12Max
angle12Max = angle12;
newIndex = j;
end
end
end
end
if newIndex ~= -1 % 若找到一个满足要求的点
tri.a = dex(1);
tri.b = dex(2);
tri.c = newIndex;
tri.indexa = -1; % 第一个领接三角形索引
tri.indexb = -1;
tri.indexc = -1;
edge.righttri = length(triout); % 设置边左侧三角形索引
isTriExist = 0;
% 记录边1
for k = 1 : length(edgeHull)
tempEdge = edgeHull(k);
if tempEdge.start == newIndex && tempEdge.end == dex(1)
tempEdge.righttri = length(triout);
tri.indexb = tempEdge.lefttri;
isTriExist = 1;
break;
elseif tempEdge.start == dex(1) && tempEdge.end == newIndex
tempEdge.lefttri = length(triout);
tri.indexb = tempEdge.righttri;
isTriExist = 1;
break;
end
end
% 若不存在这条边就新建一条边
if isTriExist == 0
newEdge.start = -1; newEdge.end = -1;
newEdge.lefttri = -1; newEdge.righttri = -1;
newEdge.start = dex(1);
newEdge.end = newIndex;
newEdge.lefttri = length(triout);
edgeHull = [edgeHull newEdge]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%????
end
isTriExist = 0;
% 记录边2
for k = 1 : length(edgeHull)
tempEdge = edgeHull(k);
if tempEdge.start == dex(2) && tempEdge.end == newIndex
tempEdge.righttri = length(triout);
tri.indexa = tempEdge.lefttri;
isTriExist = 1;
break;
elseif tempEdge.start == newIndex && tempEdge.end == dex(2)
tempEdge.lefttri = length(triout);
tri.indexa = tempEdge.righttri;
isTriExist = 1;
break;
end
end
if isTriExist == 0
newEdge.start = -1; newEdge.end = -1;
newEdge.lefttri = -1; newEdge.righttri = -1;
newEdge.start = newIndex;
newEdge.end = dex(2);
newEdge.lefttri = length(triout);
edgeHull = [edgeHull newEdge];
end
tri.indexc = edge.lefttri;
triout = [triout tri];
end
end
end
toc
%% figure the result
for i = 1 : length(triout)
dex=[triout(i).a, triout(i).b, triout(i).c,triout(i).a];
plot(pt(dex,1), pt(dex,2),'b-')
end
end此文章版权归snailgoers所有，如有转载，请注明來自原作者
微信
相关推荐
1月 14 2015
只有旋转和平移的变换矩阵推导
此处推导只针对二维情况，三维情况肯定也是合适的，暂时没有推导假设二维平面中有两条曲线，求一组变换矩阵，使得两条曲线的间距最小，自然而然大家第一个想到的就应该是最小二乘法，推导步骤如下所示 但是在求得theta后，因为有正负之分，可能会相差180°，此时，再需比较一下两次变换后的距离，即可求得所需旋转矩阵及平移量 ###code： sx1=sum(p1(:,1));sy1=sum(p1(:,2)); sx2=sum(p2(:,1));sy2=sum(p2(:,2)); sx1y2=sum(p1(:,1).*p2(:,2)); sy1x2=sum(p1(:,2).*p2(:,1)); sy1y2=sum(p1(:,2).*p2(:,2)); sx1x2=sum(p1(:,1).*p2(:,1)); theta =atan((sx1*sy2-len1*sx1y2-sy1*sx2+len1*sy1x2) /... (sy1*sy2-len1*sy1y2+sx1*sx2-len1*sx1x2)) M=[cos(theta) sin(theta) 0;-sin(theta) cos...
5月 14 2026
矩阵的旋转平移正变换及反变换
假设一个点P(x,y,z)，旋转矩阵T，平移向量V，则旋转后的平移表达式为T*P+V,此处采用右乘准则，用矩阵表示为： [T V;0 1]*[T0 V0;0 1]=[T*T0 T*V0+V;0 1] 其中T0为P原始坐标系的旋转矩阵，V0是P原始坐标系的评议向量，则相应的逆变换为 [inv(T) - inv(T)*V ; 0 1]*[T0 V0;0 1]=[inv(T)*T0 inv(T)*V0-inv(T)*V;0 1] 从上式可以看出，正变换在先旋转后平移的情况下，其逆变换是先平移后旋转
5月 15 2014
任意点到斜椭圆最小距离
最近由于项目需求，要计算任意点到斜椭圆的距离，方程相信大家都知道怎么建立，但是它是一个4次方程，matlab都很不好解，写到程序中也很难，搜了很多帖子，都是只给方程，不给解的，其实原理大家都知道，如何用程序解才是大家关心的，所以一直都没有找到合理的答案，关键时刻还是google给力，度娘一般的东西还撮合，可学术性的、英文的确是一塌糊涂，网络上收集到这样一篇论文,我也是用这个算法最终解决我的问题的，它只使用与标准方程，斜椭圆只是多了平移和旋转操作而已。 ###第一象限中的点到标准椭圆的距离### double DistancePointEllipseSpecial(const double e[2], const double y[2], double x[2]) { double distance; if (y[1] > 0) { if (y[0] > 0) { // bisect to compute the root of F(t) for t >=...
1月 14 2015
Frenet标架及曲率挠率计算
此文算法基于马国庆等论文“点云空间曲线的微分信息计算及匹配方法”，引用请注明出处，在此感谢 ###Code function [curveture,fT,fN,fB]=Frenet_Curveture(pt,k) %计算连续曲线点的frenet标架及点的曲率 %参考：点云空间曲线的微分信息计算及匹配方法 %pt：三维点集 %k：领域大小 ptNum=size(pt,1); curveture=zeros(ptNum,1); hs=floor(k/2); fT=zeros(ptNum,3); fN=zeros(ptNum,3); fB=zeros(ptNum,3); for i=hs+1:ptNum-hs Rx=pt(i,:); pj=pt(i-hs:i+hs,:); % frenet 标架计算 centerDis = sum((pj-repmat(Rx,k,1)).^2,2); scenterDis=sort(centerDis); h=scenterDis((max(floor(k/3),1)+1)); thetaj=...
5月 14 2026
直线拟合之最小二乘法AX+BY+C=0
最小二乘法拟合直线，对于大多数人都不是很陌生，直线方程 y=kx+b，令sum((y-kx-b)^2)，分别对k和b求导便可求出最佳的参数，但是如果遇到诸如“X=5”这样的直线方程，又该如何呢？这时候就不再适用了，参考这个方法，这样拟合出来的才是直线的一般方程。 Let the line sought be Ax + By + C = 0, with no restriction on A, B,and C, except not both of A and B are zero. Then the distance from apoint (x[n],y[n]) to this line is given by d[n] = |A*x[n]+B*y[n]+C|/sqrt(A^2+B^2). Then the sum of the squares of the distances of the r points from theline is F(A,B,C) = SUM d[n]^2 = SUM (A*x[n]+B*y[n]+C)^2/(A^2+B^2). H...
4月 25 2014
关于直线镜像矩阵生成
最近项目需求，在做关于直线镜像变换，网上看了一些，留作备份，以备追溯， 关于任意直线的对称变换 矩阵运算-矩阵投影，镜像，切变 以上两个是参考的两个帖子，由于我要用到关于任意直线Ax+By+C=0的镜像，采用了第一个说道的方法，并将其代码转换为矩阵表示，更加直观 对于直线ax+by+c=0有 M1=[b*b -a*b -a*c;-a*b a*a -b*c;0 0 0]/(a*a+b*b); M2=[-1 0 0;0 -1 0; 0 0 1]; M=M2*(E-M1)+M1; 其中M就是最后所求关于直线对称的镜像矩阵了
评论