matlab中的Delaunay三角剖分实现

Question

你好，这是我在这里的第一篇文章。我想为 Delaunay 三角测量编写 matlab 脚本。这是我的脚本：

clear all;clc
%% Delaunay
x=[ 160.1671 366.9226 430.7894 540.1208 660.2771 508.7287 252.1787];
y=[ 223.9615 259.5000 120.5769 245.5000 283.1923 472.7308 469.5000];

%
x=x';
y=y';

%orginal plot
dd=delaunay(x,y);
dt=TriRep(dd,x,y);
triplot(dt);
z=[x.^2+y.^2]
i=1:length(x);
ptk=[i' x y]
%% main loop
l=0;
for i=1:length(x)-2
for j=1+i:length(x)
    for k=1+i:length(x)
        if (j ~= k)
            l=l+1;
            xn = (y(j)-y(i))*(z(k)-z(i)) - (y(k)-y(i))*(z(j)-z(i));
            yn = (x(k)-x(i))*(z(j)-z(i)) - (x(j)-x(i))*(z(k)-z(i));
            zn = (x(j)-x(i))*(y(k)-y(i)) - (x(k)-x(i))*(y(j)-y(i));
                if (zn < 0)
                    border=zn;
                        for m=1:length(x)
                            border = (border) & ...
                                ((x(m)-x(i))*xn +...
                                (y(m)-y(i))*yn +...
                                (z(m)-z(i))*zn <= 0);
                            if (border) 
                            ii(m)=[i];
                            jj(m)=[j];
                            kk(m)=[k];
                            end
                        end
                end
        end
    end
end
end
wart=[ii' jj' kk']
dd
figure(2)
triplot(wart,x,y)

这是我应该从这个脚本中得到的。这个矩阵是作为 delaunay() matlab 函数生成的：

dd =
 6     7     2
 7     1     2
 4     6     2
 1     3     2
 4     3     5
 6     4     5
 2     3     4

这是我从实施中得到的：

wart =
 4     7     6
 4     7     5
 4     7     5
 4     7     5
 4     7     5
 4     6     5
 4     6     5

谁能告诉我这有什么问题？哪里有错误或只是指导我？

jils.

Answer 1

问题出在你最里面的循环中，这里：

if (zn < 0)
                border=zn;
                    for m=1:length(x)
                        border = (border) & ...
                            ((x(m)-x(i))*xn +...
                            (y(m)-y(i))*yn +...
                            (z(m)-z(i))*zn <= 0);
                        if (border) 
                        ii(m)=[i];
                        jj(m)=[j];
                        kk(m)=[k];
                        end
                    end
            end

您想检查点 [i,j,k] 定义的三角形是否有效。只有当所有m（外接圆内没有点）都是如此时，您才想将这三个点保存到您的输出中。目前，如果您检查的第一个点在外接圆之外，那么无论如何都会保存这三个点。此外，由于您为每个可能的三角形循环了相同的m，即使您找到了正确的值，您以后也可能会覆盖它们。这也是您在输出中得到相同值重复的原因。

在这些情况下，总是值得单步执行（心理上、在命令行上手动或使用调试方法）您的循环，看看会发生什么。您的第一个输出 4 7 6 不会出现在内置函数结果中。因此，将您的 i,j,k 设置为这些值，看看在那个内部循环中会发生什么。

顺便说一句，你实际上并不需要一个循环。通过执行以下操作一次检查所有值：

 border = (x-x(i)).*xn + (y-y(i)).*yn + (z-z(i)).*zn;
 if all(border<0)
    % then store coordinates
 end

您可以从空输出 ([]) 开始并追加（使用 end+1），或者计算三角形的最大数量并将输出预分配到该大小，使用计数器变量来跟踪找到的数量并将它们放在输出数组中的正确位置，然后将输出修剪到最后的大小。如果您计划拥有更大的输入数据集，最好预先分配。