Установка многоугольника в MATLAB

Определение проблемы

Я пытаюсь восстановить невыпуклую и несамопересекающуюся ломаную линию (которую я называю целевой линией) из зашумленного набора данных (см. рисунок ниже). Поэтому моя цель — узнать положение вершин целевой линии из набора данных. Я называю предполагаемую ломаную линию ломаной линией.

Целевая линия предполагается нормированной, т.е. заключенной в интервалах [-1,1] (горизонтальное направление) и [0,1] (вертикальное направление). При этом целевая линия начинается из точки [1 0]' и заканчивается в точке, пересекающей горизонтальную ось системы отсчета.

Основными проблемами этой проблемы являются следующие:

• Набор данных неупорядочен, а это означает, что единственная доступная информация — это положение зашумленных наблюдений. Порядок относится к временной последовательности, в которой мы наблюдали бы наблюдения, если бы путешествовали вдоль целевой линии от ее начальной точки до конечной точки. Важно отметить, что если указан порядок, проблему можно легко решить с помощью стандартных решателей аппроксимации кривой (например, путем минимизации методом наименьших квадратов). Однако это не мой случай...

• Алгоритм должен быть быстрым даже для больших наборов данных (т. е. время вычислений должно быть порядка секунд, даже если набор данных содержит 1e4 или даже больше наблюдений), потому что мне нужно реализовать его в гораздо более крупной системе, которая работает в режиме реального времени. . Таким образом, я не могу рассматривать причудливые алгоритмы, основанные на задачах оптимизации.

• Алгоритм должен быть надежным, когда набор данных небольшой (например, если количество наблюдений порядка 10 или меньше) или когда стандартное отклонение гауссовского шума велико (например, больше 1 в текущих нормализованных единицах измерения). . Под устойчивым я подразумеваю, что в этих сложных сценариях предполагаемая форма должна быть чем-то простым и осмысленным, например полуэллипсом, а не случайно переплетающейся линией. Другими словами, алгоритм должен проявлять осторожность: если информация ограничена, ответ должен быть простым и расплывчатым.

• Количество вершин предполагаемой ломаной должно быть равно заданному значению N

Я могу получить приблизительную оценку с помощью довольно быстрого алгоритма (даже несмотря на то, что он не оптимизирован), поэтому я ищу предложения по его улучшению с учетом некоторых критических проблем. Я использую следующие соглашения:

• hatV — матрица размера 2xN, содержащая координаты вершин N. формирование ломаной линии;
• D — матрица размером 2xm, содержащая координаты наблюдений m, формирующих набор данных. Такие точки генерируются путем равномерного выбора точки на целевой линии и последующего добавления гауссовского шума.

Мой алгоритм принимает на вход D, N и возвращает hatV.

Мой алгоритм (обзор)

function hatV = staticShaper(D, N)
    % step 0: unfrozen initialization 
    hatV   = 1.5 * [cos(linspace(0, pi, N)); sin(linspace(0, pi, N))];
    freeze = zeros(1, N);

    for iter = 1 : 3
        % step 1: project the polygonal line over the dataset
        idx = projectLine(hatV, D);
       
        % step 2: update the polygonal line
        hatV = updateLine(hatV, D, idx, freeze);

        % step 3: simplify the polygonal line
        [hatV, Ndeficit] = simplifyLine(hatV);

        % step 4: interpolate and freeze the polygonal line
        [hatV, freeze] = interpolateLine(hatV, Ndeficit);
    end
end

• Шаг 0: По предположению, целевая линия, которую я хочу оценить (и соответствующий набор данных D), полностью содержится в интервалах [-1, 1] (горизонтальное направление) и [0, 1] (вертикальное направление). Соответственно, я инициализирую ломаную линию hatV, размещая вершины N вдоль полукруга радиуса 1,5, гарантируя, что hatV охватывает целевую линию. Кроме того, я использую вектор freeze, содержащий флаг для каждой вершины. Если freeze(i) == 1, то i-я вершина заморожена, то есть ее положение нельзя изменить. Если freeze(i) == 0, то i-я вершина не заморожена, что позволяет изменить ее положение.
• Для простоты я рассматриваю цикл с заранее заданным номером итерации. Помимо этого упрощения, я думаю, что смогу доказать, что мой алгоритм всегда сходится за конечное число итераций, потому что в конце концов все вершины замораживаются.
• Шаг 1: Для каждой отдельной вершины текущей оцененной ломаной линии я ищу в наборе данных ближайшее наблюдение. Здесь idx(i) указывает, что ближайшее к вершине i наблюдение является idx(i)-м наблюдением в наборе данных D.
• Шаг 2. Обновите ломаную линию, разместив каждую незамороженную вершину. над его ближайшим наблюдением.
• Шаг 3: Предыдущая операция может привести к размещению нескольких вершин над одним и тем же наблюдением, что приведет к избыточности. На этом этапе ломаная линия упрощается за счет объединения вершин, расположенных близко друг к другу, что устраняет избыточность. В ходе этого процесса количество вершин может уменьшиться с N до N - Ndeficit.
• Шаг 4: Чтобы конечный результат содержал N вершин, я вставляю Ndeficit новых вершин в упрощенную ломаную линию. Вновь вставленные вершины представляют собой степень свободы, а старые вершины считаются фиксированными. Поэтому я замораживаю старые вершины и обозначаю вновь вставленные вершины как незамороженные, позволяя обновить их позиции в следующем цикле алгоритма.

Вопросы

На этом рисунке изображены четыре разных шага моего алгоритма за три итерации с учетом N = 50. Точки обозначают незамороженные вершины, а точки в кружках обозначают замороженные вершины. Черные линии соединяют незамороженные вершины с ближайшими к ним наблюдениями в наборе данных.

Меня больше всего беспокоит вопрос о том, как повысить точность результата. В моем решении есть две основные проблемы.

• проблема 1: мой алгоритм не может проникать внутрь «жестких» вогнутых областей.
• проблема 2: мой алгоритм не может соответствовать набору данных, а скорее очерчивает «внешнюю» границу набора данных.

Итак, мой первый вопрос: как мне решить эти две проблемы? Помимо этого вопроса у меня возникает следующий: как мне доказать (или опровергнуть), что мой алгоритм генерирует ломаную, не пересекающуюся сам с собой?

Мой алгоритм (подробнее)

function idx =  projectLine(hatV, D)
    idx = NaN(1, size(hatV, 2));
    for i = 1 : size(hatV, 2)
        [~, idx(i)] = min(vecnorm(hatV(:, i) - D));
    end
end

Для каждой вершины hatV(:, i) я ищу ближайшее наблюдение в D.

function hatV = updateLine(hatV, D, idx, freeze)
    for i = 1 : size(hatV, 2)
        if not(freeze(i))
            hatV(:, i) = D(:, idx(i));
        end
    end

    hatV(:, 1)   = [1 0]';
    hatV(2, end) = 0;
end

Я помещаю каждую незамороженную вершину hatV(:, i) над ближайшим к ней наблюдением D(:, idx(i)). После этого я проверяю, что первая вершина hatV(:, 1) находится в точке [1 0]', и подтверждаю, что конечная вершина hatV(:, end) не имеет вертикального компонента. Это потому, что я знаю, что начальной точкой целевой линии является [1 0]' и что целевая линия заканчивается на горизонтальной оси системы отсчета.

function [hatV, Ndeficit] = simplifyLine(hatV)
    oldHatV    = hatV; % backup 
    k          = 1;
    hatV       = [];
    hatV(:, 1) = oldHatV(:, 1);
    for i = 1 : size(oldHatV, 2)
        if norm(oldHatV(:, i) - hatV(:, k)) < 0.01
            hatV(:, k) = 0.5 * (hatV(:, k) + oldHatV(:, i));
        else
            k          = k + 1;
            hatV(:, k) = oldHatV(:, i);
        end
    end
    Ndeficit = size(oldHatV, 2) - size(hatV, 2);
end

Должен сказать, что эта, казалось бы, простая задача отнимает у меня некоторое время, и я даже не уверен, что это хорошая реализация. Однако здесь я стремлюсь исключить одинаковые или похожие вершины. Для этого я исследую расстояния между вершинами. Если расстояние достаточно мало, я выполняю операцию слияния, взяв среднее значение. В конце я подсчитываю, сколько вершин потеряно при этих операциях.

function [hatV, freeze] = interpolateLine(hatV, Ndeficit)
    % step 1: normalized edge lengths computation
    ell = NaN(1, size(hatV, 2) - 1);
    for i = 1 : size(hatV, 2) - 1
        ell(i) = norm(hatV(:, i + 1) - hatV(:, i));
    end
    ell = ell / sum(ell);
    
    % step 2: new vertices generation
    newVertContainer = cell(1, size(hatV, 2) - 1);
    for i = 1 : size(hatV, 2) - 1
        nrNewVertices       = round(Ndeficit * ell(i)); 
        alpha               = 0;
        newVertContainer{i} = [];
        for j = 1 : nrNewVertices
            newVertex           = (1 - alpha) * hatV(:, i) + alpha * hatV(:, i + 1);
            newVertContainer{i} = cat(2, newVertContainer{i}, newVertex);
            alpha               = alpha + 1 / nrNewVertices;
        end
    end
    
    % step 3: output definition
    oldHatV = hatV;
    hatV    = oldHatV(:, 1);
    freeze  = 1;
    for i = 1 : size(oldHatV, 2) - 1
        hatV   = cat(2,   hatV,   newVertContainer{i}, oldHatV(:, i + 1) );
        freeze = cat(2, freeze, zeros(1, size(newVertContainer{i}, 2)), 1);
    end
end

Идея состоит в том, чтобы вставлять новые вершины «равномерно» вдоль текущей ломаной линии. Чтобы добиться этого, я добавляю к каждому краю ломаной линии несколько новых вершин, пропорциональное длине ребра.

Первый шаг — вычислить длину каждого края предполагаемой ломаной линии и нормализовать результаты. Если ребро i имеет нормализованную длину ell(i), то количество новых вершин, сгенерированных между двумя старыми вершинами, определяющими ребро, равно round(Ndeficit * ell(i)).

Понятно, что эта стратегия не гарантирует, что ломаная линия в итоге будет иметь ровно N вершин, но на данный момент эта проблема не является большой проблемой (исправлю позже).

На шаге 2 я генерирую новые вершины, которые временно сохраняются в newVertContainer. Наконец, на шаге 3 я определяю обновленную ломаную линию, в которой старые вершины замораживаются, а новые вершины размораживаются.