nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Аффинный алгоритм, ч.5 - нахождение масштаба

Часть 1 - нахождение наиболее отдалённой точки
Часть 2 - сортировка против часовой стрелки
Комментарий об операторе упорядочивания
Часть 3 - нахождение матрицы преобразования
Часть 4 - нахождение крена
Математическая интерлюдия
Часть 5 - нахождение масштаба
Финал!

Мы нашли формулу, по которой надо считать масштаб:


(мы давно знали её, но теперь подтвердили её "универсальность")

Осталось написать её на нашем ассемблере, который умеет складывать, умножать, брать модули и максимумы, всё в 16 битах, и больше ничего. Да раз плюнуть!

Соответствующий код занимает 23 слова (46 байт), не требует длинных, итеративных вычислений, и обеспечивает точность не хуже 0,4%.


Вот он:

    ;состояние регистров:
    ;X=QuatY,
    ;Y=Z=AfTransf,
    ;i=j=k=0,
    ;Inv = 0,
    
    FindScale   proc
    ;1. Вычисляем величины |X-Z| (т.е |Txx - Tyy| ) и |2Y| (например, |Txy+Tyx|, хотя можно 2|Txy| или 2 |Tyx| )
    ;важно, что в данный момент Inv=0
    ;X = QuatY
    ;Y = Z = Txx = AfTransf
                    [SP]        0   ;манхэттенская метрика,
                    [SP+1]      0   ;чебышевская                
                    X           AfTransf
                    k           3                   
                    i           1
        @@loop:     Acc         [X+i]   ;Txy (итер. 0), Txx (итер. 1)
                    PM          [X+i^k] ;+Tyx (итер. 0), -Tyy (итер. 1)
                    ABS         Acc     ;|Txy+Tyx| (ит. 0), |Txx-Tyy| (ит. 1)
                    [SP+2]      Acc     ;сохранили это значение
                    ADD         [SP]    
                    [SP]        Acc ;обновили манхэттенскую
                    Acc         [SP+2]
                    MAX         [SP+1]
                    [SP+1]      Acc ;обновили Чебышевскую
                    Inv         1       ;так что на следующей итерации будет "-"
                    iLOOP       @@loop
    ;теперь в [SP] лежит манхэттенская метрика, а в [SP+1] - чебышевская
    ;осталось посчитать выражение целиком - сумму Txx+Tyy, а эти - с весами 0,664 и 0,336
    
    ;кстати, у нас до сих пор k=3 :)
    ;i=j=0
                    Acc         [X+i]   ;добавили Txx
                    ADD         [X+i^k] ;и Tyy
    ;теперь наши метрики берём с весами
                    X           MetricW
                    i           1
        @@finExpr:  C           [SP+i]
                    FMA         [X+i]
                    iLOOP       @@finExpr
    ;ага, сделали
    FindScale   endp


Основная идея описана в посте как быстро оценить длину вектора на плоскости - вместо корня из суммы квадратов (эвклидовой метрики) берём сумму модулей (манхэттенская метрика, она же метрика городских кварталов) и максимум из модулей (чебышевская метрика), с определёнными весами, которые обеспечили минимальную ошибку. Само по себе такое приближение - не самое точное, можно ошибиться аж на 5.5%, но по счастью, этот квадратный корень должен получаться существенно меньше остальных слагаемых, в наихудшем случае составляя 6,7% от всего результата. Поэтому, если мы в этих 6,7% ошибёмся на 5,5%, то суммарная ошибка и составит около 0,4%.

Но здесь есть небольшое техническое дополнение. При выводе формулы предполагалось, что матрица симметричная - если бы наши вычисления были абсолютно точными, так и должно было получиться. Но как мы увидели "воочию", в реальности она может таковой и не оказаться. У нас получилось Tyx=0, тогда как Txy=1. Поэтому, вместо 2Y (где Y должен был быть ЛЮБЫМ из этих двух компонентов), мы решили взять Tyx+Txy. Как предполагается, это немножко повысит точность, и довольно гладко вписывается в алгоритм.

Когда мы вычисляли крен, нужно было взять сумму элементов по главной диагонали, и разность по оставшимся. Здесь всё наоборот: разность по главной диагонали и сумму по оставшимся. Мы оставили такой же цикл по i, от 1 до 0, и ту же самую "укуренную" адресацию через XOR, а вот регистр инверсии в этот раз надо "переиначить". В тот раз мы приравнивали его переменной i, а в этот раз - просто приравниваем единице, но В САМОМ КОНЦЕ цикла. Зная, что с прошлых процедур он был нулевым, поэтому на итерации, где i=1, мы сложим два числа, а на следующей - уже вычтем.

Итак, смотрим по строкам. На первой итерации (i=1):
Acc         [X+i]


загрузит Tyx.

PM          [X+i^k]

ПРИБАВИТ Txy, поскольку с прошлых процедур стоит Inv=0.

ABS         Acc

Acc по шине поступает на вход АЛУ, берётся абсолютное значение (модуль), которое оказывается снова в аккумуляторе.

[SP+2]      Acc

сохраняем его в "локальную переменную".

ADD         [SP]

прибавляем к результату текущее значение Манхэттенской метрики, а это ноль (мы обе метрики обнулили "ручками" перед началом цикла).

[SP]        Acc

заносим новое значение, а конкретнее, |Txy+Tyx|.

Acc         [SP+2]

восстанавливаем "текущее" значение, всё то же |Txy+Tyx|

MAX         [SP+1]

Если значение, поступившее по шине данных (т.е [SP+1]) больше, чем лежащее в аккумуляторе, то оно поступает в аккумулятор,
т.е Acc = max (Acc, Src)

В данном случае мы сравниваем с нулём, поэтому результат очевиден, |Txy+Tyx| остаётся на своём месте.
[SP+1]      Acc

заносим его ещё и в Чебышевскую метрику.

Inv         1

подготавливаемся к следующей итерации, где вместо суммы будет разность.

Ну и начинаем вторую итерацию, i=0.
Acc         [X+i]

загружаем Txx.

PM          [X+i^k]

ВЫЧИТАЕМ Tyy, поскольку теперь Inv=1.

ABS         Acc

берём модуль.

[SP+2]      Acc

сохраняем его в "локальную переменную".

ADD         [SP]

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

[SP]        Acc

заносим новое значение, а конкретнее, |Txy+Tyx|+|Txx-Tyy| - всё верно.

Acc         [SP+2]

восстанавливаем "текущее" значение, то есть |Txx-Tyy|

MAX         [SP+1]

Находим max (|Txx-Tyy|, |Txy+Tyx|)

[SP+1]      Acc

заносим его ещё и в Чебышевскую метрику - всё верно.

Inv         1

Уже не нужно, но "кашу маслом не испортишь".

На этом выполнение цикла завершается, и у нас "на руках" (точнее, на стеке) лежит две метрики.

Дальше мы снова пользуемся тем фактом, что задействованные в циклах переменные уже обнулились, а k=3 мы как задали - так оно и лежит. Это в очередной раз даёт нам отсрочку в плане сделать адресацию "базовый регистр плюс отступ" - до сих пор худо-бедно обходились без неё, вот и здесь получилось.

Мы занесли в аккумулятор Txx:
Acc         [X+i]


прибавили к нему Tyy:
ADD         [X+i^k]


и затем нужно внести метрики с соответствующими весами. Веса мы добавили в оперативную память:

MetricW:
ManhW   dw      11010   ;вес, с которым взять Манхэттенскую метрику
ChebW   dw      21758   ;вес Чебышевской метрики


и устроили очередной цикл:
                    X           MetricW
                    i           1
        @@finExpr:  C           [SP+i]
                    FMA         [X+i]
                    iLOOP       @@finExpr


В итоге, к окончанию процедуры у нас в аккумуляторе лежит удвоенное значение масштаба. Мы абсолютно уверены, что оно там "умещается", благодаря "лишнему биту", т.е когда 16-битные значения мы интерпретируем как лежащие в диапазоне -1..1, в аккумуляторе помещаются -2..2.

Посмотрим, как это работает, сначала на примере 300 метров (наши первые точки).


Мы остановились после вычисления метрик, они лежат в st[1] и st[2]. Для матрицы


у нас выходит Манхэттенская метрика, равная 3, и Чебышевская, равная 2. Действительно, модуль разности элементов на главной диагонали: 2, модуль суммы остальных: 1. Поэтому |2|+|1|=3, тогда как max(|2|,|1|)=2. Всё чётко.

Также видим, что в st[3] осталось "текущее значение" с последней итерации, т.е 2, регистры i,j действительно нулевые, k=3. Всё по плану.

Теперь остановим в начале следующей процедуры:


В памяти мы показали, где хранятся веса для метрик, тоже в области "для чтения" (хотя об этом знает только эмулятор, "железо" будет без понятия!). И видим значение Acc: 874,336. Мы отображаем старшие 17 разрядов как "целые", а младшие 15 - как "дробные", для удобства. Именно в этом месте у нас нигде не прибавилось "0,5" для округления "до ближайшего целого", потому как дальше мы захотим сдвинуть это значение влево "до упора", и тогда такая прибавочка недопустима!
Если бы мы посчитали числитель нашей формулы "честно" (через корень из суммы квадратов), то получили бы 874,236, т.е ошибка составила 0,011%.

Кстати, наш план по снижению ошибки за счёт усреднения двух значений, которые должны были оказаться одинаковыми, сработал! Если вместо суммы элементов вне главной диагонали, взять лишь один из них, но удвоенный, ошибка окажется выше - 0,236, если взять ноль, или 0,436, если взять единицу, т.е 0,027% .. 0,05%. Но в целом, здесь мы в хороших условиях - подкоренное выражение чрезвычайно мало, и будь наши измерения яркостных центров точными (это не так: мы постарались промоделировать все пакости, какие только могут случиться - неравномерность чувствительности, "мертвые зоны" между пикселями, дискретность АЦП и пр.), и вовсе должно было оказаться нулевым.

Поэтому возьмём более злой пример, который скорее всего является наихудшим случаем для этого алгоритма. А именно, отклонимся от нормали на 30°, причём не влево и не вверх (по этим осям мы получаем безупречные вычисления), а под углом 45°. Ниже изображено, как видна наша картинка при взоре "по нормали" (слева) и "под углом" (справа). И разумеется, дальность 300 метров, когда все числа маленькие, и у нас начинаются проблемы с нехваткой разрядности.



Видно, насколько сильно изменилась картинка - мало того, что она ужалась как раз-таки "по диагонали", так ещё и центральная точка уползла, расположившись почти на прямой линии с двумя своими соседями :) Это сказался её вынос на 29 см.

Вот получившиеся координаты точек:
Fx0		Int16	425
Fy0		Int16	22
Fx1		Int16	85
Fy1		Int16	-39
Fx2		Int16	-541
Fy2		Int16	-301
Fx3		Int16	285
Fy3		Int16	-885


Как обычно, они расположены "в порядке сканирования", сверху вниз:


Такой ракурс ещё опасен тем, что наиболее отдалённая точка стала уже не такой отдалённой, вот проверим - определяется ли она верно?

После FindMDD3 имеем:


Хехе, нашёл как миленький! Кстати, как начал делать скриншоты-вспомнил, что после игр с иконками "только для чтения" и "звёздочка" забыл выставить моноширинный шрифт, поэтому в окне памяти всё было вкривь да вкось. Подправил.

Запускаем SortCCW, и после него получается:


Финальный порядок такой:


Да, всё как надо.

Находим матрицу аффинного преобразования:


То есть:


Довольно необычная матрица: когда возюкаешься с поворотами, привыкаешь, что в одном углу "синус", в другом - "минус синус", а тут ЧЕТЫРЕ ПОЛОЖИТЕЛЬНЫХ ЧИСЛА! Ну да, сплющивание чуть по-другому работает :)

Находим косинус и синус крена:


Они соответствуют углу atan(-561/32764) = -0,98°. Это ошибочка вышла - крен должен был быть нулевым, но центральная точка "уползла" и немножко сбила нас с панталыку. Но не страшно, алгоритм сопровождения разберётся, что к чему. А именно ошибка вычисления угла из матрицы составила 0,022°.

Не хотел показывать, как считается кватернион, но "он на меня обиделся" и заявил, что я обратился к неинициализированной памяти. Оказалось, я забыл подправить эмулятор, чтобы он не обновлял флаг знака на команде DIV2, вот оно так и вышло. Так что всё-таки покажу:


Всё в порядке.. Ошибка вычислений составила 0,0018°, или 6'' - как видно, когда мы выходим на "полный размах", точность получается очень приличная.

Теперь посмотрим, как из матрицы "выкидывается" крен:


Отлично, матрица действительно стала симметричной:


Наконец, вычисляем масштаб:


Результат содержится в аккумуляторе: 877,376. Если из той же матрицы посчитать его "честно" (через квадратный корень), мы получим 873,756, то есть ошибка составила 0,41%, на них мы и расчитывали.

Сравнивая масштаб с полученным в прошлый раз, видим их неплохое совпадение, т.е "перекрестная связь" (когда изменение одной величины влияет на измерение другой) практически отсутствует. В этот раз "корень" внёс существенный вклад: простая сумма элементов главной диагонали дала всего 800, а остальные 73 - за счёт "корня". На стеке, st[1] и st[2], мы видим "оценку сверху" (манхэттенскую метрику) 88, и "оценку снизу" (Чебышевскую метрику) 72.

Так что метод вполне работоспособный.

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


Продолжение следует...
Tags: tex, ПЛИС, кватернионы-это просто (том 1), математика, программки, работа, странные девайсы
Subscribe

Recent Posts from This Journal

  • Лестница для самых жадных

    В эти выходные побывал на даче, после 3-недельной "самоизоляции". Забавно, как будто зима началась! Особенно грязные галоши остались на улице, в…

  • Возвращаемся к макету

    Очень давно макетом видеоизмерителя параметров сближения не занимался: сначала "громко думал" по поводу измерения его положения на аппарате, а потом…

  • Минутка живописи

    В процессе разгребания содержимого квартиры (после нескольких ремонтов) дошёл, наконец, и до картин. В кои-то веки их повесил. Куда их вешать -…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 0 comments