Реализация триангуляции Делоне в Matlab

Здравствуйте, это мой первый пост здесь. Я хочу написать сценарий 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)

Это то, что я должен получить из этого сценария. Эта матрица генерируется как функция matlab delaunay ():

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

Может ли кто-нибудь из вас сказать мне, что с этим не так? Где ошибка или просто наставьте меня?

джилс.


person ddr8    schedule 08.06.2014    source источник
comment
Комментарии необходимы для понимания кода на семантическом уровне. Например: вы инициализируете границу с отрицательным двойным значением, но затем оно используется как логическое значение. Инициализация border=false будет иметь такой же эффект. Что предназначено?   -  person Daniel    schedule 09.06.2014
comment
Я так думаю. Я пишу сценарий Matlab для этой программы на C: fatcat.ftj .agh.edu.pl / ~ kaprzyk / Geomkomp / geomkomp / node58.html Однако этот сценарий делает различия между функцией delaunay и программой на языке C. Аналогичная реализация находится здесь в Java: stackoverflow.com/questions/5825089/   -  person ddr8    schedule 09.06.2014
comment
@jiksu: border=(zn < 0); - правильный перевод c-кода.   -  person Daniel    schedule 09.06.2014


Ответы (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) или вычислить максимальное количество треугольников и предварительно выделить свой вывод для этого размера, использовать переменную счетчика, чтобы отслеживать, сколько вы найдете, и поместить их в нужное место в выходном массиве, а затем обрезать выходные данные до размера прямо в конце. Если вы планируете иметь большие наборы входных данных, было бы лучше выделить их заранее.

person nkjt    schedule 09.06.2014
comment
Это было очень полезно. Я выполнил ваши инструкции и получил треугольники. Второй вопрос - как удалить внутренние строки, которые мне не нужны. Здесь у меня нет идей. - person ddr8; 09.06.2014