nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Как быстро оценить длину вектора на плоскости?

Возникла задачка: посчитать выражение
для произвольных 16-битных целых чисел x,y, используя только сложение-вычитание-умножение и, возможно, условные переходы.

Причём большая точность не нужна, 5-10% уже хватает за глаза.

Как оказалось, можно взять такое забавное выражение:


и получить максимальную ошибку в 5,5%. Вполне сойдёт!

Если не хочется умножать на константы, можно взять выражение попроще, обходящееся сдвигами:

но точность будет похуже: 11,8% в наихудшем случае.

Это выражение навеяло метриками: Эвклидовой, Чебышевской и Манхэттенской. Про последние известно, что одна всегда даёт значение выше, чем эвклидова, другая - всегда ниже. А вот если их взять с некоторыми весами - то получится очень даже пристойно :)


А вот задача, для которой это понадобилось:


Мы посчитали матрицу 2х2, выражающую аффинное преобразование, которое мы можем разложить на крен, масштаб и "ракурсы":



Здесь:
T (от слова Transform) - наша матрица аффинного преобразования,
S (Scale) - матрица масштабирования, s-масштаб, причём s>0 (мы не можем подлететь "с обратной стороны" и увидеть изображение зеркально отражённым!)
A (Aspect) - матрица "ракурсов" (когда мы смотрим на плоскость не под прямым углом):



здесь a - величина "ужатия", причём

поскольку мы знаем заранее, что отклонились от нормали к плоскости не более чем на 30°.

ψ - направление, вдоль которого происходит "ужатие".

Наконец, R (Rotation) - матрица поворота на угол ϕ.

Так вот, мы хотим, зная матрицу T, найти кватернион поворота:


и масштаб s. Разбираться с ракурсом мы не хотим - всё равно определим его неоднозначно (что мы отклонились вправо, что влево - а картинка ужмётся одинаково!), и очень криво (попробуй заметить отклонение на 5°, когда оно ужимает картинку по горизонтали на cos(5°)=0,9962).

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



очень неприятная матрица, но мы можем заметить две вещи:
- она симметричная
- оба диагональных элемента обязаны быть положительными, при наших ограничениях на коэффициенты s и a.

Обозначим как-нибудь её элементы:


Именно её мы должны получить, если умножим нашу матрицу T на матрицу, обратную к матрице поворота:



то есть, мы получаем следующие условия, из которых можем определить матрицу поворота:




Последнее выражение можно переписать, приведя подобные слагаемые:


откуда мы знаем, чему должен равняться tgϕ - это даёт два угла поворота, отличающиеся на 180°. Но с учётом оставшихся двух условий (ненулевые диагональные элементы), из этих двух углов нас устроит только один:




где α>0 - нормирующий множитель.

(Точнее, противоположное значение (когда α<0) нас заведомо не устроит, поскольку, если мы сложим диагональные элементы, мы получим:


если α<0, то сумма диагональных элементов отрицательна, а значит, они не могут оба быть положительными!
Если же α>0, то сумма диагональных элементов положительна, но это ещё не гарантирует, что оба элемента будут положительными. Если это не так, это скорее всего значит, что мы "зацепились" не за те точки - восприняли какие-то блики как наши мишени. С чистой совестью выдаём признак ошибки - это гораздо лучше, чем на полной скорости долбануться...)

Это замечательно! Мы сначала находим "ненормированные значения", и затем нормируем их до посинения, воспринимая 16-битные целые числа в формате 1.15 (от -1 до 1-2-15), используя метод Ньютона (см. ликбез) - метод всегда сходится, хоть и может потребовать много итераций (когда значения слишком маленькие, они на каждом шагу будут умножаться на 1,5). Но у нас время есть.

Далее, в части 6 3/4 показано, как из синуса и косинуса получить кватернион поворота (где синус и косинус половинных углов!) - в нашей ситуации, когда они уже отнормированы, мы обходимся прибавлением единички и повторной нормировкой, причём здесь много итераций уже не надо - норма изначально будет в диапазоне 1/2..1. На этом первая часть нашей задачи завершена.

Осталось найти масштаб s. Мы легко строим матрицу обратного поворота, и помножаем на неё нашу исходную матрицу.

Наконец-то, все коэффициенты этой страшной матрицы нам известны:


Если мы сложим два диагональных элемента, получим:


И есть второе интересное выражение:


поэтому масштаб находится так:


До сих пор на всём протяжении мы сумели обойтись одними лишь сложениями и умножениями, неужели здесь под самый конец нам всё-таки придётся "честно" посчитать квадратный корень???

Нет, не придётся :)

Если мы просто проигноруем корень, положив его равным нулю, мы вместо s получим

что даст наихудшую ошибку, когда ракурс составляет 30°, и a=cos(30°)≈0,87. Ошибка составит 7%. Многовато - по ТЗ нам нужно 1,7%.

Ладно, запишем


это даст максимальную ошибку, когда |X-Z|=|2Y|, а именно при вычислении корня ответ получится в корень из двух раз меньше, чем надо. Это всяко лучше, чем совсем выкидывать корень - и мы получим точность не хуже 2%.

Возможно, и этого вполне хватит - то что мы сейчас рассматриваем - это режим захвата, за которым начинается режим сопровождения. Главное только, чтобы мы "втянулись" в сопровождение, не сорвались с крючка.

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

Выходит, что для режима захвата тоже хватает чистой арифметики, без синусов-косинусов-корней и прочих нехороших штук :)


В плане проектирования Quat Core остался последний вопрос - нужно ли мне деление, или его уже лучше исполнить "программно", через метод Ньютона или вовсе "в столбик"?
Tags: кватернионы-это просто (том 1), математика, работа
Subscribe

  • Нахождение двух самых отдалённых точек

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

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Балансируем конвейер QuatCore

    В пятницу у нас всё замечательно сработало на симуляции, первые 16 миллисекунд полёт нормальный. А вот прошить весь проект на ПЛИС и попробовать "в…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 4 comments