Аппроксимация дугой окружности в 3D или...

[Zerald]

Ребята подскажите плиз алгоритм для аппроксимации дугой окружности набора точке в 3D пространстве.

Задача такая:
Есть набор точек, ожидается что точки лежат на дуге окружности в пространстве (получаемых от акселерометра), т.е. фактически имеем вращение вокруг некой оси. Нужно найти эту ось вращения и угол поворота вокруг нее.
Вот сделал скриншот http://serzhpwr.newmail.ru/RotAroundZ.gif - тут вращения вокруг оси Z
Делать это надо в real time, т.е. точки летят от датчика и ось вращения может меняться со временем т.е. чем меньше понадобится точек для аппроксимации тем лучше.

Задачу можно немного упростить, сделав допущение что ось вращения будет всегда проходить через начало координат. В этом случае можно расчитать плоскость в которой лежат эти точки, следовательно индексы плоскости и будут осью вращения. Но из скриншота видно что соседние точки лежат нифига не втой плоскости в которой нужно и конечно ось вращения не соответствует действительности. Как-то усреднять точки перед расчетом плоскости тоже не гуд так как потом возникнут проблемы с расчетом угла вращения

[незваный гость]

Задачи такого типа бывают разными. Например, если сигнал меняется во времени (хотя и слабо), применяются фильтры с конечным откликом. Если предполагается, что окружность одна и та же, уместно применение Кальмановских фильтров. Если берется, каждый раз, скажем, 100 точек, то наилучший результат даст метод наименьших квадратов: постепенно накапливаете суммы для вычисления параметров окружности, в конце — однократное, но сравнительно простое вычисление.

Смотрите литературу по цифровой обработке сигналов, DSP.

[Zerald]

Сигнал меняется не чтобы со временем, а как отклик на движение (поворот) при этом ось вращения (соответственно окружность) может меняться поэтому накапливать точки не особо получится, максимум несколько десятков.
В принципе можно было бы применить МНК, но вывод формул МНК в случае зависимости по дуге окружности я не осилю точно, тем более что не знаю уравнения окружности в 3D. Вывести его(уравнение) тоже не получилось (ниасилил)

[Алексей К.]

Zerald писал(а):

вывод формул МНК в случае зависимости по дуге окружности я не осилю точно, тем более что не знаю уравнения окружности в 3D. Вывести его(уравнение) тоже не получилось (ниасилил)

Вывод МНК- формулы для окружности на плоскости прост. И задачки --- провести МНК плоскость, и потом на плоскости возиться с МНК-окружностью наверняка в паре эквивалентны возне с 3D-окружностью.
Не надо фитировать сумму квадратов расстояний от точки до окружности.
Достаточно минимизировать

$$\sum [(x_i^2+y_i^2)-2Bx_i-2Cy_i+Q]^2$$

,
то есть отклонения от нуля значения левой части неявного уравнения окружности.
Можно показать, что при малых отклонениях это близко к честному МНК.
Ну, разумеется работать надо в системе координат центра тяжести точек --- все несколько упрощается. $B,C$ --- координаты центра окружности, из $Q,B,C$ потом определяется радиус:

$$(x-B)^2+(y-C)^2=R^2,\qquad x^2+y^2-2Bx-2Cy+\underbrace{B^2+C^2-R^2}_Q=0.$$

Добавлено спустя 2 минуты 46 секунд:

По-моему, зря в Computer Science переместили...
Ежели сейчас заставят МНК решать, производную брать, линейную систему...

Добавлено спустя 42 минуты 58 секунд:

Re: Аппроксимация дугой окружности в 3D или...

Zerald писал(а):

Задачу можно немного упростить, сделав допущение что ось вращения будет всегда проходить через начало координат. В этом случае (выделено мной --- АК) можно расчитать плоскость в которой лежат эти точки, следовательно индексы плоскости и будут осью вращения.

Рассчитать эту плоскость можно в любом случае, просто при этом допущении что-то станет ноликом и якобы упростится.

[незваный гость]

Мне по-прежнему не очень понятна задача, в частности:

1) Какими вычислительными ресурсами Вы располагаете. Скажем, будет ли Вас пугать обращение матрицы порядка 7? Если нет, то (после некоторых упражнений), можно использовать некоторую модификацию фильтра Кальмана. Правда, его придется более или менее регулярно сбрасывать. С другой стороны, если Вас не пугают более тяжелые вычисления, в фильтр можно было бы включить движение оси и изменение радиуса окружности.

Если же таковых ресурсов нет, то можно сделать вариант МНК.

2) Доступны ли только пространственные координаты точек, или временные тоже? Насколько надежны временные координаты, если доступны?

~~

Если использовать МНК, первым шагом (как я думаю) будет нахождение плоскости точек. Я бы применил МНК, минимизируя $\sum (a x_k + b y_k + c z_k)^2$. Получится, разумеется, 3 уравнения, единственное решени которых $a = b = c = 0$. Чтобы обойти это, я бы искал решение не на шаре $a^2+b^2+c^2 = 1$, а на кубе $|a|=|b|=|c|=1$, и выбрал бы из трех вариантов решение с $a^2+b^2+c^2 \to \min$. После этого можно считать все точки лежащими в плоскости, проходящей через центр тяжести и перпендикулярной этому вектору.

Далее, по желанию, можно, разумеется, этот вектор развернуть вдоль, например, оси $Z$, сделав задачу двумерной.

[Алексей К.]

незваный гость писал(а):

Если использовать МНК, первым шагом (как я думаю) будет нахождение плоскости точек. Я бы применил МНК, минимизируя $\sum (a x_k + b y_k + c z_k)^2$. Получится, разумеется, 3 уравнения, единственное решени которых $a = b = c = 0$. Чтобы обойти это, я бы искал решение не на шаре $a^2+b^2+c^2 = 1...

Обходить ничего не надо. Задача с плоскостью (как и дополнительная к ней задача с пространственной прямой) решается довольно просто. Надо найти собственные векторы матрицы (пишу всё по памяти, надеюсь не ошибаюсь, когда-то сам всё это программировал)

$$ M= \left| \begin{array}{lll} \sum x_i^2 & \sum x_i y_i & \sum x_i z_i\\ \sum x_i y_i & \sum y_i^2 & \sum y_i z_i\\ \sum x_i z_i & \sum y_i z_i & \sum z_i^2 \end{array} \right| $$

(в системе центра тяжести). Ищутся они легко, особенно если программкой, а не формулой...
Собственные значения неотрицательны; соб. вектор, соответствующий максимальному с.з., даёт направление наилучшей МНК-прямой. Чтобы найти его, достаточно взять какое-то первое приближение $v_0$ и слегка поитерировать: $v_{i+1}=M v_i$. При отклонениях, от плоскостности, характерных для машиностроения (самых грубых), мне хватало 3х итераций, чтобы получить компьютерно-точный результат. Первое приближение вообще можно брать от фонаря --- лишь бы оно случайно не оказалось перпендикулярно искомому вектору.

Но этот фокус проходит для максимального с.з. А в случае плоскости нам нужет вектор от минимального с.з. Для этого вместо матрицы $M$ следует взять для итераций матрицу ... как бишь она называлась? присоединённая, кажется, adjoint. Элементы подменены своими минорами.

ниже незваный гость писал(а):

И транспонируются.(поправлено для полноты --- АК)

Можно было бы взять и обратную матрицу, но $M$ при идеальной плоскости или прямой будет вырождена. Эта, присоединённая, легко считается. С ней и итерируем.

[незваный гость]

Простите, что-то мне странно это.

1) Берем три точки: $(1, 0, 0)$, $(0, 1, 0)$, $(0, 0, 1)$. $M$, по-видимому, единичная матрица. Ее собственные вектора описаны в литературе — бери любой, чтобы не было обидно и не ссылатся на случайное совпадение значений собственных чисел (и, как следствие, несходимость итераций). А правильный ответ, как я понимаю — $(1, 1, 1)$

2) В поиске на кубе решаются три системы уравнений второго порядка. При этом делается 18 умножений и 6 делений. Плюс, 9 умножений для сравнения модулей. Думаю, с итерациями Вам придется потеть не меньше (27 умножений на итерацию).

Фишка в том, что вектор может случайно оказаться (почти) параллельным одной из осей, и в этом случае точность очень упадет. Решение трех систем позволяет обойти это. На самом деле, конечно, достаточно двух, но три — симметричнее.

[Алексей К.]

Я писал(а):

$$ M= \left| \begin{array}{lll} \sum x_i^2 & \sum x_i y_i & \sum x_i z_i\\ \sum x_i y_i & \sum y_i^2 & \sum y_i z_i\\ \sum x_i z_i & \sum y_i z_i & \sum z_i^2 \end{array} \right| $$

(в системе центра тяжести).

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

Поиск на кубе, признаться не анализировал. Предложенное мной решение вроде как точное (итерации можно заменить решением кубического уравнения). Не уверен, что его придумал я сам --- возможно, вычитал в какой-то серьёзной книге о матрицах как пример задачи на с.з.
Проблемы real-time тоже никогда не анализоровал. Лет 10 назад у меня это считалось на 386 ПиСи из фортрановской прграммы, и потом приборчик ездил по вычисленной плоскости. А уж сейчас --- подсчёт умножений, скорее, дань традиции...

Я не хочу глубоко копать (самому уже неинтересно + условия командировки), хочу лишь быть уверенным, что не наврал и изложил малоизвестный/интересный подход к задачке (по-моему, основной).
Тем более, что аскер уже убежал.

[незваный гость]

Алексей К. писал(а):

Матрица перестала быть единичной?

Перестала. Я, действительно, не заметил, что система центра масс. Буду думать дальше…

Алексей К. писал(а):

А уж сейчас --- подсчёт умножений, скорее, дань традиции...

Пока Вы не начали такты процессора и доли микросекунды считать. Поверьте, это и сейчас встречается. Когда Вам надо обработать 20 каналов за одну миллисекунду (и сделать что-нибудь еще), начинаете считать и не только умножения… И таких задач не мало.

Добавлено спустя 13 минут 21 секунду:

присоединенная (adjoint)

Алексей К. писал(а):

Элементы подменены своими минорами.

И транспонируются.

[Xey]

незваный гость в сообщении #74195 писал(а):

Если же таковых ресурсов нет, то можно сделать вариант МНК.

2) Доступны ли только пространственные координаты точек, или временные тоже? Насколько надежны временные координаты, если доступны?

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

Можно бы учесть эту неравномерность, если ввести веса для отсчетов , зависящие от плотности их размещения на окружности. Но это это наверно очень сложно.

[Алексей К.]

Это не сложно. Делаете первую итерацию без всяких весов, и по полученной окружности рассчитываете плотности-весá для второй итерации. Третья итерация не нужна, она практически ничего не изменит.

Не могя похвастаться глубоким знанием статистических методов, я обычно делал выводы о весах, поправках, браковках после рассмотрения глазками конкретных результатов. Что за окружность, что за измерения, возможны ли выбросы (заусенец на профиле), и проч. Эмпирикой это, кажется, называется.

[Xey]

Алексей К. в сообщении #254126 писал(а):

Это не сложно. Делаете первую итерацию без всяких весов, и по полученной окружности рассчитываете плотности-весá для второй итерации.

Для этого надо сохранять данные, суммировать, не сохраняя проще.
А если , по коэффициентам предыдущего измерения разбить окружность на секторы и усреднять точки, попадающие в каждый сектор.
Это должно подровнять веса.

[Xey]

Тут рядом обсуждают похожий вопрос

Aniosher в сообщении #243771 писал(а):

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

И есть такое мнение

Ajabsandal в сообщении #246556 писал(а):

А Вам обязательно требуется решить эту задачу методом наименьших квадратов?
Часто построение нейронной сети (обучение на экспериментальных данных) дает лучшие результаты, чем МНК. Напр. RBF (Radial basis function) - нейросети с одним скрытым слоем обучаются быстро и дают хорошую точность приближения искомой функции.
И, в отличие от апроксимационных методов, нейросети дают хороший результат и за границами области апроксимации, что особенно ценно при отсутствии достоверных (или полное их отсутствие) экспериментальных данных.
Вы можете воспользоваться многими программами для НС, их много с сети.
Лучшие из них - NeuralWorks Professional II/Plus и NeuroSolution 6.

Интересно, в этих классических случаях, когда сигнал прост, применение сети дань моде или она может обеспечить лучшее приближение к реальности. Известно же, что метод НК некоторая условность, можно минимизировать и первые и вторые и третьи степени отклонений .

[ewert]

Xey в сообщении #256587 писал(а):

Известно же, что метод НК некоторая условность, можно минимизировать и первые и вторые и третьи степени отклонений .

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

[Xey]

ewert в сообщении #256611 писал(а):

остальные степени .... гораздо хуже считаются.

Это что, единственное достоинство квадратов отклонения?
Поэтому я и спрашиваю может быть сеть лучше?