Подання рухів у 3D моделюванні: інтерполяція, апроксимація та алгебри

У цій статті мені б хотілося розповісти про один цікавий математичному прийомі, який будучи вельми цікавим і корисним мало відомий широкому колу людей, що займаються комп'ютерною графікою.

Скільки існує різних способів представити звичайний поворот у тривимірному просторі? Більшість людей, коли-небудь займалися 3D-графікою, 3D-моделюванням, відразу назвуть три основних широко поширених варіанти:

  • Матриця повороту 3x3;
  • Завдання повороту через кути Ейлера;
  • Кватерніони.
Люди з багатим досвідом додадуть сюди чомусь не користується популярністю четвертий пункт:
  • Вісь повороту і кут.
Мені б хотілося розповісти про п'ятому способі подання обертань, який симпатичний тим, що зручний для параметризації, дозволяє ефективно будувати поліноміальні наближення цих параметризаций, проводити сферичну інтерполяцію, і головне, універсальний — з мінімальними змінами він працює для будь-яких видів рухів. Якщо вам коли-небудь був потрібен метод, який дозволяв би легко зробити «аналог slerp, але не для чистих обертань, а для довільних рухів, та ще й з масштабуванням», то читайте цю статтю.

Я почну з математичного опису того, звідки виникає п'ятий спосіб подання обертань, саме є виникаюча параметризація, але всі бажаючі можуть пропустити цей розділ і відразу перейти до практичної частини в кінці статті.

Групи і алгебри: на стику алгебри і диференціальної геометрії
Обертання об'єкта в тривимірному просторі (і, узагальнюючи, будь руху об'єкта) з математичної точки зору є типовим прикладом структури яка називається «групою». Груповий операцією тут є «комбінування», коли спочатку застосовується один рух, а потім другу; ця комбінація, очевидно, теж є рухом. Одиницею групи служить тривіальне рух «нічого не змінювати». Для будь-якого руху можна відшукати зворотній рух, що дозволяє «повернути все як було». З практичної сторони, стосовно до геометрії, існування групи часто говорить про існування якоїсь «особливості» цього руху, яка зберігається в межах групи. Наприклад, якщо ми пересунули кудись об'єкт чистим обертанням, то повернути його назад можна теж чистим обертанням. Група рухів евклідового простору зберігає відстані між точками передвигаемого об'єкта. Група симетрій куба складається з 48 рухів, які переводять куб в самого себе. Це те, що більш-менш відомо всім.

Легко однак зауважити, що група обертань влаштована в певному сенсі «цікавіше» чим група симетрій куба. Так, ми знаємо, що для обертань можна інтерполювати безперервний рух між будь-якими двома точками (знайти функцію, яка плавно і безперервно поверне об'єкт з одного положення в інше), тоді як для 48 елементів групи симетрій куба це, очевидно, неможливо. Ще одним цікавим властивістю є можливість локально параметризовывать елементи цієї групи трійками чисел (наприклад, кутами Ейлера). У математичних термінах всі ці властивості формалізують фрази «група обертань є гладким різноманіттям». Тобто ми можемо уявити собі простір всіх можливих поворотів в тривимірному просторі як деяку багатовимірну поверхню типу сфери. Тільки от звичайна сфера є двовимірною поверхнею і її можна помістити (вкласти) в тривимірний простір, а простір обертань є тривимірної гиперповерхностью для того щоб її зобразити потрібно чотиривимірний простір; а щоб зробити це без самоперетинів — шестимерное, причому отримана в цьому просторі поверхню для групи обертань топологічно буде нагадувати швидше бублик, ніж сферу. Але математика там виходить практично такий же, тому я буду використовувати звичайну двовимірну сферу для інтуїтивної ілюстрації поняття «гладкого різноманіття».

Розглянемо тепер безперервний рух якогось об'єкта, заданий функцією положення об'єкта в залежності від моменту часу f(t). Скажімо, f(0) може визначати начально положення об'єкта, f(1) — кінцеве, а f(0.5) тоді буде відповідати положенню об'єкта на середині дороги» від початкової точки до кінцевої. Оскільки окремі положення об'єкта є точками на нашому гладкому різноманітті, вважаючи f(t) в різних точках, можна побудувати на нашій «сфері» ланцюжок з точок, що зв'язує між собою початкову та кінцеву точки на «сфері»; в межі, як нескладно зрозуміти, ця ланцюжок об'єднується в безперервну криву. Тобто різні криві на нашому різноманітті є різними траєкторіями руху об'єктів і навпаки.


Рис. 1. Група обертань як гладке різноманіття (сфера): точки відповідають початковому і кінцевому положеннях об'єкта, що проходить по сфері крива — траєкторії його руху

Згадаємо тепер, що у нас все-таки не деяке абстрактне різноманіття, а група; і у цієї групи є дві операції: «множення» і «взяття оберненого елемента». Ми можемо діяти цими операціями на «сфері». Наприклад, ми можемо зробити т. н. «лівий зсув на елемент g», домножив всі точки на сфері зліва на g, тобто замінивши кожний елемент a g*a. Оскільки у нас група, то це рух зберігає нашу «сферу» в цілому, але переміщує всередині неї окремі зазначені точки. Інтуїтивна інтерпретація: раз у нас група поворотів, то ми можемо повернути нашу сферу будь-яким з цих поворотів.


Рис. 2. Дія групи на собі як поворот сфери

Якщо подібні операції групи на самій собі (множення і звернення) не псують гладкості раніше заданих на ній кривих (є гладкими відображеннями сфери на себе), то така група називається групою Чи. Практично всі значущі групи рухів є саме групами.

Наступне цікаве спостереження, яке ми можемо зробити, полягає в тому, що у гладких кривих на гладких многовидах можна визначити поняття швидкості руху по кривій стосовного вектора. Тобто в нашому випадку рухомого об'єкта можна визначити швидкість обертання тіла в будь-який момент часу, причому ця швидкість буде векторною величиною. У диференціальної геометрії є теорема, яка стверджує, що простір всіх можливих дотичних векторів до кривих в заданій точці n-мірного різноманіття є n-мірним ж лінійним простором дотичною площиною в заданій точці різноманіття. Наприклад, швидкість обертання тіла можна визначити віссю обертання з кутовою швидкістю обертання. Якщо ми перемножимо одиничний вектор напрямку осі на скаляр кутової швидкості, то вийде вектор кутовий швидкості, напрямок якого задає вісь, а довжина — швидкість обертання. Цей вектор і буде елементом цього стосовного простору до групи обертань. На малюнку я використовую в якості моделі двовимірну сферу, тому і дотичний простір там виходить двовимірної дотичною площиною до цієї сфери.


Рис. 3. Швидкість обертання як дотичний вектор до траєкторії руху об'єкта і дотична площина — як простір всіх можливих швидкостей обертання об'єкта в обраній точці

Однак у силу того, що дотична площина задається для конкретної точки різноманіття, поки що ми визначили «швидкість обертання об'єкта» тільки для однієї конкретної точки (конкретного стану об'єкта), а хотілося б — для будь-якого можливого становища. Тут-то і приходить на допомогу характерне для груп «гладке дію різноманіття на самого себе»: визначивши дотичний вектор і дотичну площину всього в одній точці різноманіття (як правило, в одиниці), ми далі можемо пересунути отримані вектор і площина в будь-яку іншу точку різноманіття. Для цього нам, правда, доведеться зафіксувати якийсь конкретний спосіб, яким ми будемо це робити, щоб уникнути неоднозначностей, і класичний підхід — це використовувати вже згадуваний лівий зсув, т.е. використовувати домножение зліва на g для того, щоб порахувати дотичний вектор і дотичну площину в точці g. Цього легко дати очевидну інтерпретацію: допустимо, ми визначили, що таке «обертання вліво зі швидкістю 5 градусів в секунду» для об'єкта знаходиться в положенні A. Як визначити «обертання об'єкта вліво зі швидкістю 5 градусів в секунду» в іншому положенні Б? Для цього треба просто спочатку дати об'єкту повернутися в старому положенні, а потім (множення зліва) помістити об'єкт у потрібне місце таким же поворотом, який привів би нерухомий об'єкт з положення А в Б. Така нотація, природно передбачає, що рухи діють справа наліво: дія a*b відповідає дії b, після якого слідує дія a. Таким чином, визначивши швидкість обертання в одній точці можна довизначити поняття «швидкості обертання» для всіх точок простору. Отримане векторне поле (математик додасть: левоинвариантное векторне поле) на групі Чи природно асоціювати з «просто» швидкістю обертання, вже без уточнення в якій саме точці.


Рис. 4. Визначена в одній точці «швидкість обертання» природним чином доопределяется для всіх точок групи обертань

Як я вже писав раніше, простір всіх можливих кутових швидкостей в одній точці є тривимірним дотичним простором. Оскільки зазначене побудова працює для будь-яких дотичних векторів, то слова «в точці» виходить прибрати і тут; тобто виходить просто «простір можливих кутових швидкостей об'єкта», і цей простір є звичайним лінійним простором всіх мислимих тривимірних векторів. Цей простір називається алгеброю Лі відповідає даній групі. Чому алгеброю? А тому що для цього простору вже природним чином визначена операція додавання і множення на скаляр (простір лінійно) і на ньому завжди можна певним нескладним чином ввести спеціальну операцію множення» Чи дужку з наступними властивостями:

  1. [x,y] є білінійною операцією щодо x y;
  2. [x,x] = 0 для будь-якого x;
  3. [x, [y,z]] + [y, [z,x]] + [z, [x,y]] = 0 для будь-яких x,y,z.
Для нашого простого випадку, до того ж, на підставі пунктів 1 і 2 слід пункт:
  1. [x,y] = -[y,x] для будь-яких x,y.
Стосовно до нашого випадку групи обертань, алгеброю, як ми вже зрозуміли, є звичайний тривимірний простір векторів, а «множенням» в ньому виступає добре всім знайоме векторний добуток.

Групи і алгебри грають помітну роль в математиці, оскільки алгебра Лі (лінійний простір), як правило, є набагато більш простим об'єктом, ніж вихідна група. Скажімо, у нашому модельному прикладі вихідна група була хитрим тривимірним бубликом, що вкладаються в шестимерное простір, а алгебра Лі — вже звичайне тривимірне простір. Для вихідних елементів групи була визначена тільки операція комбінування обертань і вона була не комутативність, тоді як елементи алгебри можна складати між собою і множити на скаляри. Але при цьому алгебра однозначно визначає відповідну групу і дозволяє визначати її властивості (наприклад, визначати тип групи) через властивості алгебри Лі, які, за великим рахунком, зводяться до розмірність лінійного простору і матриці, що описує дужку через її дію на базисних векторах. Однак для практичних цілей нам знадобиться інша важлива властивість.

Розглянемо довільну гладку траєкторію руху f(t) в нашій групі. Як вже було сказано, її можна продиференціювати, визначивши швидкість руху у кожній з точок траєкторії. Візьмемо тепер довільний елемент x з алгебри (довільну кутову швидкість для нашого модельного прикладу). Як раніше згадувалося, цей елемент відповідає гладкому векторному полю x(g) дотичних векторів до всіх елементів g різноманіття групи.
Розглянемо тепер наступне диференціальне рівняння:


Іншими словами — що буде, якщо ми стартуємо з одиничного перетворення в момент t=0 і будемо після цього постійно обертати наше тіло з фіксованою кутовою швидкістю x? У диференціальної геометрії і теорії звичайних диференціальних рівнянь доводиться, що для будь-якого наперед заданого елемента x з алгебри це рівняння має єдине можливе рішення f(x, t), яке називається інтегральної кривої. Значення цієї функції в момент часу t=1 показує «як далеко ми підемо» з початкового положення за одиницю часу при русі з фіксованою кутовою швидкістю x і називається експоненціальним відображенням алгебри на групи:


Операція экспоненцирования однозначно визначена для всіх елементів алгебри і володіє цілим рядом корисних властивостей:
  • Exp(0) = 1;
  • Exp(-x) = (Exp(x))-1;
  • Exp(n*x) = (Exp(x))n для цілих n;
Але найважливіше те, що експонента дозволяє тривіально реконструювати згадану інтегральну криву f(x,t) руху з постійною кутовою швидкістю f(x,t) = Exp(tx).

Причому ця крива є геодезикой на вихідному різноманітті групи і отже локально найкоротшою траєкторією, що з'єднує між собою точку 1 з будь-якою іншою точкою. З допомогою лівого зсуву цей результат можна тривіально узагальнити для траєкторії руху між довільними положеннями a b. Досить просто порахувати g=ab-1, знайти таке x, для якого Exp(x) = g=ba-1, і тоді геодезикой, що з'єднує a b f(t) = a Exp(xt).

Ну і в загальному випадку геодезикой, відповідної руху з постійною швидкістю x і проходить через будь-який елемент g f(t,x,g) = g Exp(xt).

Таким чином, экспоненцирование елементів алгебри дозволяє природним чином визначити оптимальні інтерполяційні траєкторії у вихідній групі. Згадана при цьому «зворотна» операція, що знаходить для заданого елемента групи g відповідну кутову швидкість обертання x, як неважко здогадатися, називається логарифмированием Exp( Ln( g) ) = g.

Однак на відміну від експоненти вона, як правило, визначена неоднозначно; а для деяких елементів групи логарифм може бути не визначено взагалі. Наприклад, якщо група Чи складається з декількох зв'язних компонент (умовно кажучи, дві сфери, а не одна), то інтегральні криві, очевидно, не можуть покинути компонент містить одиничний елемент і логарифм для будь-яких елементів, що лежать поза цього компонента, не визначено. Фізичний сенс цього досить тривіальний — якщо ми розглядаємо, наприклад, групи обертань і дзеркальних симетрій в тривимірному просторі, то траєкторії, яка б безперервним чином переклала б об'єкт до його дзеркального відображення, просто не існує, а неоднозначність визначення обертання пов'язана з тим, що якщо об'єкт безперервно обертати, він регулярно буде повертатися до одного і того ж положення.

Матричне экспоненцирование і логарифм
З практичної точки зору важливо відзначити той невеликий факт, що практично всі використовувані руху можна записати у вигляді дійсних матриць того чи іншого виду. Наприклад, група обертань у тривимірному просторі описується матрицями 3x3 з одиничним детермінантом, група евклідових рухів (rigid body transformation) — матрицями 4x4 наступного виду:


Які теж мають одиничний детермінант і т. д. Відповідні групи є підгрупами загальної Загальної Лінійної Групи GL(R,n) невироджених матриць розмірності nxn. Ця група є групою, і їй відповідає алгебра Чи gl(R,n), що складається з всі матриць розмірності nxn, дужкою для якої є матричний комутатор: [A,B]=AB-BA. І в цій алгебрі природно виникають операції матричного экспоненцирования матричного логарифмування, позначаються Exp Ln відповідно. Принадність цього факту полягає в тому, що всі інші приватні групи рухів є підгрупами цієї групи і їм відповідають лінійні підпростори в алгебрі gl(R,n). При цьому вони ділять з більш загальною групою і дужку, і операції матричного экспоненцирования, тобто незалежно від конкретної підгрупи рухів, експонента і логарифм визначаються однаково — через експоненту і логарифм відповідних матриць.

Крім уже згадуваних загальних властивостей:
  • Exp(0) = 1;
  • Exp(-A) = (Exp(A))-1;
  • Exp(n*A) = (Exp(A))n для цілих n;
  • f(x,g, t) = g Exp(tx) — геодезика, що проходить через g і відповідна швидкості x
    у матричного экспоненциала є ряд додаткових корисних властивостей;
  • Якщо AB=BA Exp(A+B) = Exp(A) Exp(B);
  • (ряд Тейлора) причому цей ряд всюди сходиться;
  • (ряд Меркатора), сходиться при |A| < 1;
  • ;
  • .
Ефективне обчислення матричної експоненти і логарифма є складною проблемою, і найбільш розумним підходом до його реалізації є використання відповідної математичної бібліотеки :).

Давайте подивимося, як виглядають матриці алгебри для деяких практично значимих груп. Наприклад, наша проста група обертань у матричному вигляді записується як «всі матриці 3x3 з одиничним детермінантом», а відповідна їй алгебра Лі складається з матриць 3x3 наступного виду:



Де (x,y,z) — це вже згадуваний вектор кутовий швидкості, що визначає швидкість обертання за принципом «довжина вектора = кутова швидкість, напрям = вісь обертання». X, y z тут можна розглядати як кінематичний еквівалент кутів Ейлера — по суті це кутові швидкості обертання відносно осей x, y z. Але в цілому отримуємо просто інший спосіб запису обертань у поданні «вісь і кут».

Якщо ж ми розглядаємо випадок довільного руху в тривимірному просторі, записуючи його 4x4 матрицями виду , де R — це 3x3 матриця обертання, а T — це 3x1 вектор зсуву, то відповідна алгебра виходить утвореній з матриці:



Де (rx,ry,rz) має вже згаданий зміст вектора кутової швидкості, а (tx,ty,tz) — це вектор миттєвої швидкості тіла в точці (0,0,0). Відповідні пари векторів називають мотором миттєвих швидкостей твердого тіла, і з їх використанням пов'язаний витончений і порівняно маловідомий математичний апарат так званого гвинтового обчислення, яке заслуговує окремої статті на Хабре.

Хочемо додати ще й масштабування? Просто додаємо діагональні елементи (sx,sy,sz), відповідні швидкості «розтягування» по кожній з осей:


А ось для rigid body руху в 2D, яке можна задати матрицями 3x3 алгебра виходить з матриць:


… ну і так далі, думаю, ідею ви вже зрозуміли :). Але все-таки…

Навіщо все це потрібно?
Інтерполяція
В якості першого прикладу давайте візьмемо простеньку 2D-фігуру і спробуємо її пересунути безперервним рухом з одного положення в інше. Параметризуем положення моделі (x,y) і кут повороту (a), проводимо лінійну інтерполяцію…


Все ніби працює, проте є якась слабка, але явно помітна неприродність руху. У чому ж справа? Давайте намалюємо траєкторії, по яких рухаються точки, розташовані спочатку в лівій частині екрана.


Нижня частина, де знаходиться наш літачок, не дуже-то схожа на поворот: різні точки рухаються не тільки з різною швидкістю, але і за різними траєкторіями; однак саме дивне відбувається у верхньому лівому кутку, де точки здійснюють характерне зворотно-поступальний рух. Я не став робити відповідну ілюстрацію, оскільки на практиці подібний артефакт зустрічається не дуже часто, але думаю, що серед тих, хто займається тривимірною анімацією, подібну штуку при інтерполяції хоча б пару раз спостерігали багато.

А що буде, якщо ми візьмемо замість цього експоненту? Нехай A — це матриця вихідного положення об'єкта, а B — матриця кінцевого положення. Покладемо положення об'єкта в момент часу t , де .



Як бачимо, вийшов чистий плавний поворот, в якому всі точки об'єкта рухаються з постійною швидкістю. Звичайно, такого ж результату можна було б досягти і не вдаючись до експоненті, а просто згадавши, що будь-який рух на площині є або паралельним перенесенням, або обертанням навколо деякої точки на деякий кут. Достатньо було б підібрати відповідний центр обертання і поварьировать кут. Але принадність експоненти і логарифма полягає в тому, що це абсолютний no-brainer, всього дві-три рядки коду, які з однаковим успіхом будуть працювати у всіх обставинах. Потрібно інтерполювати рух в 3D? Та ж сама формула автоматично видасть вам гвинт, навіть якщо ви не знаєте теорему Шаля. Потрібно інтерполювати лінійне переміщення? Та ж сама формула дасть лінійний перехід. Цікавим і корисним попутним властивістю экспоненцирования як способу інтерполяції є те, що при інтерполяцію переходу між елементами підгрупи експонента «працює» суто всередині цієї підгрупи. Тобто якщо, приміром, у вас в тривимірному просторі є робот, який може рухатися тільки в певній площині, або, скажімо, обертається відносно певної точки, то інтерполяція між його положеннями автоматично буде відбуватися тільки в цій же площині або обертанням відносно цієї точки.

Параметризація та апроксимація
Інша корисна властивість випливає з того, що правильно задана Exp:
  1. завжди дає елемент цієї групи рухів;
  2. розкладається в ряд Тейлора.
Ну тобто, припустимо, ми хочемо знайти матрицю обертання, яка буде найкращим чином задовольняти певній вимозі. Скажімо, у нас є два шматочки однієї і тієї ж поверхні, і ми хочемо склеїти їх разом, для чого складаємо матрицю обертання, яка найкращим чином поєднає перший шматочок з другим. Щоб реалізувати пошук такої матриці, нам спочатку доведеться якось завдний простір цих матриць. Це можна зробити, наприклад, за допомогою кутів Ейлера, але коли ми підставимо отриману матрицю з купою синусів і косинусів, то буде незрозуміло як вирішувати отриману оптимізаційну проблему. Ми можемо замість цього спробувати використовувати в якості параметрів для оптимізації безпосередньо елементи самої матриці руху; при цьому вийде куди більш проста оптимізаційна проблема (квадратичне рівняння у разі ICP), але немає ніяких гарантій, що порахована таким чином матриця афінного перетворення буде обертанням. Існує кілька варіантів обходу цієї проблеми: наприклад, у ряді робіт пропонується вважати афінну матрицю, а потім ортогонализировать її до «найближчого» обертання. Экспоненцирование дає один з цікавих альтернативних варіантів. Як вже говорилося на початку, алгебра Чи є значно більш простим «аналогом» групи (обертань), яка нас цікавить. Тому:

  1. ми параметризуем елементи алгебри (а це звичайне лінійне простір);
  2. використовуємо експоненту щоб повернутися до цієї групи.
На перший погляд може здатися, що цей підхід ще гірше озвученої на початку параметризації кутами Ейлера, оскільки матричну експоненту в оптимізаційної задачі вважати ще незручніше, ніж синуси і косинуси. Однак для невеликих рухів замість експоненти ми можемо просто записати її розкладання в ряд Тейлора:


Для дрібних A цей ряд сходиться дуже швидко (для як завгодно великих він теж сходиться, але повільно), в силу чого його можна просто обрубати на першому-другому слагаемом після одиничної матриці. Наприклад, надзвичайно часто для малих кутів використовують наближення


Воно параметризует матрицю обертань добре відомим всім наближенням


Правда, зазвичай коли це наближення вводять, то не згадують ні експонент, ні ряду Тейлора :). Однак «ноги» у цього наближення ростуть саме з ряду Тейлора, і тепер ви знаєте, як цю апроксимацію без праці продовжити до 2, 3 і взагалі будь-якого порядку для довільних класів рухів (а не тільки обертань). До речі звідси ж росте добре відомий OpenCV-стам приватний випадок для групи тривимірних обертань — формула обертання Родрігеса.

Використання подібного роду апроксимацій дозволяє спростити оптимізаційну задачу (наприклад, звести її знову ж до квадратичних) і при цьому не вимагає сумнівної переходу з рішенням задачі для афінних матриць і подальшого проектування результату на потрібну групу. Єдиним її обмеженням є вимога до малості матриці A, але як правило, обійти його досить просто. Досить просто завдний не безпосередньо «матрицю положення об'єкта», а «матрицю руху, який потрібно застосувати до об'єкта». За умови, що об'єкт спочатку варто більш-менш десь в правильному місці, необхідна матриця «зрушення» в цілому класі задач виходить досить маленькою. Наприклад, добре цей підхід працює в ітераційних алгоритмах, оскільки крок переміщення об'єктів на кожній наступній ітерації в них, як правило, швидко зменшується.

Ось така проста, цікава і, на подив маловідома штука. А Ви використовуєте матричні експоненти і алгебри у своїх проектах :)?
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.