Найти наиболее вероятный прямоугольник среди множества точек

%На вход функции подаются два вектора х и у (координаты точек)
%На выходе в векторе rez передаются номера точек, которые наиболее вероятно
%составляют прямоугольник
%porog - погрешность
function rez = rect(X,Y,porog)
    %Составим матрицу, которая будет хранить расстояния между точками
    %Создадим нулевую матрицу 
    %Пройдёмся по строкам
    T = zeros(length(X));
    for i = 1:length(X)-1
        %Пройдёмся по столбцам
        for j = i+1:length(X)
                %Производим рассчёт расстояния
                T(i,j) = sqrt((Y(j)-Y(i))^2 + (X(j)-X(i))^2);
        end;
    end;
    T
    %Матрица сформирована, начнём её анализ
    for i = 1:length(X)-1
        %Пройдёмся по столбцам
        for j = i+1:length(X)
            %Выбрали элемент, и дальше будем проходить по элементам,
            %которые находятся ниже строк с номерами i и не учитываем
            %строки с номером j
            for a = 1:length(X)-1
                %Пройдёмся по столбцам
                    for b = i+1:length(X)
                        if a > i %смотрим строки ниже чем строка a
                            if a ~= j %не рассматриваем строки с номером j
                                if (T(i,j) <= T(a,b) + porog) && (T(i,j) >= T(a,b) - porog)
                                    %Если длины равны, то стасуем точки и
                                    %посмотрим на новые расстояния
                                    if (T(i,a) <= T(j,b) + porog) && (T(i,a) >= T(j,b) - porog)
                                        %перетасуем ещё раз
                                        if (T(i,b) <= T(j,a) + porog) && (T(i,b) >= T(j,a) - porog)
                                            %Даже если и нашлись такие точки,
                                            %то потребуем, чтобы расстояния
                                            %между ними было больше порога
                                            if T(i,j)>porog && T(a,b)>porog && T(i,a)>porog && T(j,b)>porog
                                                    disp('YES')
                                                    rez = [i,j,a,b]
                                                    return
                                            end;
                                        end
                                    end;
                                end;
                            end;
                        end;
                    end;
            end;
        end;
    end;
    rez = [0, 0 ,0 ,0];
end