Delaunay三角剖分在matlab中的实现

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，因此即使您找到正确的值，以后也很可能会覆盖它们。这也是为什么你会在输出中得到相同值的原因。

在这些情况下，总是值得单步执行您的循环(在心理上、在命令行上手动执行，或者使用debug方法)来查看发生了什么。您的第一个输出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)，或者计算三角形的最大数量并将输出预先分配到该大小，使用计数器变量跟踪找到的数量并将它们放在输出数组中的正确位置，然后在最后将输出裁剪到合适的大小。如果您计划拥有更大的输入数据集，则最好进行预分配。