В этой статье суммирую опыт перевода графического приложения на векторные операции SSE/AVX. Приложение двумерное, плоская графика. Для обработки используется упакованный double. Цель статьи обобщить опыт работы с векторными типами данных приментельно к плоской графике и показать, как можно оптимизировать операции над векторами.
Некоторые операции не могут быть эффективно описаны на языке Си, привожу примеры описания базовых операций с использованием встроенных функций и векторных расширений языка.
Я использую версию GCC-6+, но все эти наработки, встроенные функции и расширения языка будут работать начиная версии GCC-4.8+.
Для обращения к процессорным инструкциям я использую встроенные платформо-зависимые функции. Для других компиляторов (Intel или MS) следует использовать intrinsic функции. LLVM и так понимает. При переносе на другую платформу, например ARM следует переопределить базовые платформо-зависимые функции.
Тип double является избыточным для вычислений графики, достаточно float. Для платформы ARM рекомендуется такая замена. Для переопределения разрядности вычислений достаточно в одном определении заменить базовый тип и длину вектора в байтах.
Данный тип описывает точку на плоскости или радиус-вектор из начала координат. Упакованный double -- вектор, состоящий из двух элементов double описывается с использованием векторного расширения языка Си (GCC):
Базовые операции над упакованными числами
В некоторых случаях компилятор эффективно преобразует выражение вида (v2df){P[y],P[x]} в одну инструкцию SSE/AVX, а в некоторых случаях требуется явное указание инструкций. Операции унарные вида (v2df){op(P[x]), op(P[y])} и бинарные вида (v2df){op(A[x],B[x]), op(A[y],B[y])} обычно хорошо выражаются средствами языка Си. К таким операциям относятся +, −, *,/, логические операции и даже операции сравнения. Обычно не вызывает проблем смешанные операции, такие как умножение вектора на скаляр. Ниже перечислены операции над упакованными числами, которые не описываются средствами языка Си, описание требует использования обращений к инструкциям процессора.
shuffle перестановка местами компонент вектора P = (v2df){P[y],P[x]}. Всегда, не зависимо от платформы можно использовать встроенную функцию __builtin_shuffle.
shuffle2 составить вектор из элементов двух векторов, например P = (v2df){B[y],A[x]}. Для таких операций можно использовать встроенную функцию __builtin_shuffle над двумя операндами. Встроенная функция __builtin_shuffle не является платформо-зависимой, зачастую без нее не удается описать даже простые векторные операции, хотя конечно можно использовать инициализацию массива для перестановок.
addsub -- вычитание по Х компоненте и сложение по Y компоненте.
P = (v2df){A[x] − B[x], A[y] + B[y]}
hadd -- горизонтальное сложение компонент вектора P[x] + P[y]
P = (v2df){A[x] + A[y], B[x] + B[y]}
hsub -- горизонтальное вычитание компонент вектора P[x] − P[y]
P = (v2df){A[x] − A[y], B[x] − B[y]}
Скалярное произведение вектора на число P * t. Описывается средствами языка Си. Компилятор сначала преобразует число t к векторному типу T = (v2df) {t,t}, а потом производит умножение. Так же описывается сложения/вычитания и деления вектора на скаляр.
Изменение знака P * (v2df){−1,1}. В таких случаях компилятор использует умножение, а можно использовать логическую операцию XOR над знаковым битом числа:
Для изменения знака используется логическая операция XOR, это приводит к появлению чисел −0.0. Следует отметить, что при сравнении чисел 0. и −0. инструкции 387 и SSE выдают истину, (0.0 == −0.0), хотя при бинарном сравнении эти числа разные. Встроенная функция используется по причине того, что компилятор не соглашается использовать бинарный оператор ^ для упакованных вещественных чисел.
max выделение максимального значения, независимо по компонентам вектора
P = (v2df){max(A[x],B[x]), max(A[y],B[y])}
min выделение минимального значения, независимо по компонентам вектора
P = (v2df){min(A[x],B[x]), min(A[y],B[y])}
abs положительный вектор, независимо по компонентам вектора.
P = (v2df){abs(A[x]), abs(A[y])}
В числах с плавающей точкой знак задается старшим битом. Задача состоит в том, чтобы
сбросить знак минимальным затратами, в одну инструкцию. Это может быть битовая операция NAND, НЕ-И над знаком числа:
sqrt - извлечение квадратного корня
P = (v2df){sqrt(A[x]), sqrt(A[y])}
Загрузка и выгрузка векторов без выравнивания
Векторные операции требуют выравнивания данных в памяти. Эта задача возникает при переходе от скалярного представления (массива значений типа double) к векторному типу.
Например, для загрузки значений векторов из массива
double* p;
следует выравнивать массив на 16 байт или использовать операции загрузки не выровненных данных. Для работы данными без выравнивания используются встроенные функции
или выражение
которое компилятор при оптимизации заменит на инструкцию movupd.
Избыточность
тип данных остается тем же, но операцию можно производить только над младшей компонентой вектора, оставляя без изменения старшую часть вектора. Если после выполнения действий над векторами, нужно выполнить ряд действий над компонентой X или Y, то не стоит преобразовывать тип в базовый double, потому что компилятор выведет число из векторного регистра и начнет длинную вереницу инструкций плавающей точки 387, перекладывая числа через стек.
Например, нужно рассчитать корень из суммы квадратов компонент вектора:
Если эти операции выполнять на векторных регистрах SSE/AVX, то необходимо задать подряд всего три инструкции: умножение (mulpd), параллельное сложение (haddpd), и вычислить корень числа прямо на регистре SSE (sqrtsd) -- одна инструкция,- всего три инструкции.
Пример кода, после компиляции:
__builtin_ia32_sqrtpd(__builtin_ia32_haddpd(a*a,b*b));
Без встроенных функций я не нашел способа описать эти операции так же эффективно.
Если тоже самое выполнить над компонентами вектора, по описанию:
v2df v = r*r;
double d = sqrt(v[0] + v[1]);
то получается: одна инструкция SSE - умножение (mulpd), параллельное сложение (haddpd), сохранение из регистра на стек, загрузка в из стека в сопроцессор 387 и еще десять операций 387, не считая вызова функции sqrt(). И так бывает всегда при переходе с регистров xmm на базовый тип double. Поэтому чаще всего бывает выгодно выполнять все операции на регистрах SSE, даже если это приводит к избыточным операциям, когда используется только одна компонента из вектора.
Ключи компилятора GCC
Чтобы компилятор всегда выполнял операции с типом double на регистрах SSE/AVX надо задать ключ компилятору -mfpmath=sse
Чтобы компилятор не перекладывал параметры функции типа double через стек, надо включить опцию компилятора -msseregparm или использовать соответствующий атрибут функции. Например,
v2df spline(v2df ps, v2df pe, double t) {
return ps + (pe − ps)*t;
}
при задании ключа -msseregparm параметр t передается в регистре xmm2
Ниже фрагмент кода после компиляции:
О процессорных инструкциях SSE/AVX можно прочитать в документе Intel:
Intel® 64 and IA-32 Architectures Software Developer’s Manual
Справочник с годами дополняется новыми инструкциями, доступен для скачивания на сайте производителя. Там же даны соответствующие имена intrinsic функций, которые приняты в компиляторах Intel и MSVС, и поддерживаются компиляторами GCC.
Векторные расширения GCC для языка Си:
https://gcc.gnu.org/onlinedocs/gcc/
https://gcc.gnu.org/onlinedocs/gcc/x86-Built-in-Functions.html
Если надо получить ассемблерный код пригодный для изучения, и анализа использую ключ -S.
Унарные операции над векторами
Эвклидова норма вектора, его длина
v2df norm(v2df p) { return __builtin_ia32_sqrtpd(__builtin_ia32_haddpd(p*p,p*p)); }
Результат компиляции - три инструкции SSE/AVX: умножение, горизонтальное сложение, корень квадратный.
Поворот вектора на 90 градусов по часовой стрелке или против часовой стрелки.
v2df cw(v2df p){
return __builtin_ia32_xorpd(shuffle(p),(v2df){−0.,0});
}
v2df ccw(v2df p){
return shuffle(__builtin_ia32_xorpd(p,(v2df){−0.,0}));
}
Для изменения знака используется логическая операция XOR над знаком числа.
Oтражение вектора. Горизонтальное отражение меняет знак Y-компоненты вектора. Вертикальное отражение вектора меняет знак X-компоненты вектора.
v2df reflect_vert(v2df p){
return __builtin_ia32_xorpd(p,(v2df){−0.,0});
}
v2df reflect_horz(v2df p){
return __builtin_ia32_xorpd(p,(v2df){0,−0.});
}
При редактировании векторных изображений, при копировании, часто применяются операции отражения и поворота на 90 градусов.
Бинарные операции
Расстояние между точками на плоскости
norm(p1-p0);
Dot product - скалярное произведение векторов
double dot(v2df a, v2df b){
b *= a;
return b[0] + b[1];
}
Cross Product векторное произведение, Z-координата, площадь параллелограмма (удвоенная площадь треугольника):
double cross(v2df a, v2df b){
v2df v = b*shuffle(a);
return v[0] - v[1];
}
Поворот вектора a на величину радиус-вектора r по или против часовой стрелки.
rotate_cw, rotate_ccw:
v2df rotate_ccw(v2df a, v2df r){
return addsub(a*r[0],shuffle(a)*r[1]); // a_x*r_x - a_y*r_y,
}
Расстояние от точки P до прямой. Прямая задана радиус-вектором R.
d = cross(r, p);
Знак результата определяет с какой стороны прямой расположена точка P.
Точка пересечения отрезка (P1,P2) и прямой. Прямая задана радиус-вектором R.
Сплайны
Линейная интерполяция
B_1(t) = (1−t)*P_0 + t*P_1
spline (P0,P1, t) = P0 + (P1 − P0) * t
Сплайн -- параметрическое описание кривой, удобное представление, которое можно выразить минимальным числом инструкций.
Квадратичный Bézier сплайн
{B_2}(t) = (1−t)^2 *P_0 + 2t(1−t)*P_1 + t^2 * P_2
qspline(P0,P1,P2, t) = spline(spline(P0,P1, t), spline(P1,P2, t), t)
Кубический Bézier сплайн
{B_3}(t) = (1−t)^3 *P_0 + 3t(1−t)^2 * P_1 + 3(1−t)t^2 * P_2 + t^3 * P_3
cspline(P0,P1,P2,P3, t) = qspline(spline(P0,P1, t), spline(P1,P2, t), t), spline(P2,P3, t), t)
Функции сплайнов можно определять рекурсивно, как показано выше. Компилятор из такого описания делает хорошо оптимизированный код.
Касательная к кубическому сплайну, используется для движения по траектории заданной сплайном, производная от кубического сплайна - квадратичный сплайн
{B_3}'(t) = 3(1−t)^2 ({P}_1 − {P}_0) + 6(1−t)t({P}_2 − {P}_1) + 3t^2({P}_3 − {P}_2)
{B_3}'(t) = 3* qspline(P1−P0, P2−P1, P3−P2, t)
Касательная к траектории заданной кубическим сплайном -- скорость движения по траектории.
{B_2}'(t) = 2(1-t) ({P}_1 − {P}_0) +2t({P}_2 − {P}_1)
{B_2}'(t) = 2* spline(P1−P0, P2−P1, t)
Для комплекта
{B_1}'(t) = P1−P0.
Общее утверждение -- производная от кубического сплайна - крвадратичный сплай, производная от квадратичного сплайна - линейный сплайн, производная от сплайна - направление отрезка. Физический смысл производной -- скорость движения по контуру, направление вектора скорости -- касательная к траектории. Направление производной используется для построения фигур вдоль траектории.
Описание дуги окружности кубическим Bézier-сплайном
k=4/3(sqrt(2)-1), k = 0.5522847498
https://people.mozilla.org/~jmuizelaar/Riskus354.pdf
Если угол (a) меньше 90 то k*tg(a/2)
ax = x1 – xc
ay = y1 – yc
bx = x4 – xc
by = y4 – yc
q1 = ax*ax + ay*ay
q2 = q1 + ax*bx + ay*by
k2 = 4/3(sqrt(2*q1*q2)-q2)/(ax*by - bx*ay)
x2 =xc + x1 – k2*y1
y2 =yc + y1 + k2*x1
x3 =xc + x4 – k2*y4
y3 =yc + y4 + k2*x4
Преобразование линейных сплайнов в квадратичные
{B_2}(t) = (1−t)^2 *P_0 + 2t(1−t)*P_1 + t^2 * P_2
P1 = (P0+P2)/2
Преобразование линейных сплайнов в кубические
{B_3}(t) = (1−t)^3 *P_0 + 3t(1−t)^2 * P_1 + 3(1−t)t^2 * P_2 + t^3 * P_3
P1=P0+1/3*(P3−P0);
P2=P3−1/3*(P3−P0);
Преобразование квадратичных сплайнов в кубические
Дело в том, что внутреннее представление сплайнов может быть единообразным, все кривые можно свести к кубическим сплайнам. Квадратичный сплайн можно представить в форме кубического:
P1=P0+2/3*(P1s−P0);
P2=P3+2/3*(P1s−P3);
здесь P1s - ключевая точка квадратичного сплайна.
Операции с матрицами трансформации
x_new = xx * x + xy * y + tx;
y_new = yx * x + yy * y + ty;
Из коэффициентов преобразования можно составить матрицу {xx, yx, xy, yy, tx, ty}. В упакованном типе запишется так:
struct {
v2df x,y,z;
};
где x = (v2df){xx,yx}, y = (v2df){xy,yy}, z = (v2df){tx,ty};
Трансформация складывается из поворотов, масштабирования и смещения начала отсчета.
Выделяют несколько операций, из которых складывается трансформация: поворотов, масштабирования и смещения начала отсчета.
Поворот можно задать углом, или двумя параметрами sin и cos угла, вместе они образуют (v2df){cos(a), sin(a)} радиус-вектор, наклон которого определяет угол поворота.
Единичные матрицы трансформации
Умножение матриц трансформации
Смещение начала координат
Вращение координат
Масштабирование
Обратная матрица трансформации
Операции с контурами
Контур на плоскости задается точками: прямыми и сплайнами. в общем виде используются линейные и кубические Bézier сплайны.
Принадлежность точки замкнутому контуру.
Площадь контура
Сумма площадей трапеций:
S = sum{S_i}, площадь трапеции S_i =1/2 (P_{i}[x] − P_{i-1}[x]) * (P_{i}[y] + P_{i-1}[y]})
Другой способ:
Можно суммировать и вычитать треугольники относительно начала отчета. Треугольники задаются радиус-векторами точек контура. Площадь треугольника - cross product - векторное произведение векторов-сторон треугольника.
cross(a,b) = a[0]*b[1] − a[1]*b[0];
Выражается тремя операциями
v2df v = a*shuffle(b);
return v[0]-v[1];
Значение векторного произведения равно площади параллелограмма, двум площадям треугольника A0B, а знак будет описывать направление обхода контура. Берем любую точку вблизи фигуры, назначаем ее центром отсчета и относительно нее суммируем площади треугольников по контуру.
При вычислении суммы площадей можно суммировать вектора по контуру и в конце выполнить вычитание сумм, вместо того чтобы выполнять вычитание на каждом шаге цикла:
for (i=0;i<num;i++) { sum += a*shuffle(b); }
return sum[0] - sum[1];
Вычисление площади контура можно еще ускорить, если использовать тождество:
[a,b] + [b,c] = [(a − c),b]
квадратными скобками обозначено векторное произведение.
см. Алгебраические свойства векторного произведения
тогда в цикле получаем операцию
sum += (a − c)*shuffle(b); -- четыре операции на два треугольника.
Длина контура
int_{0}^{1} { norm({B_3}'(t)) dt}
Физика
В общем виде физика движения объекта описывается кубическими сплайнами. Координаты центра масс задаются в зависимости от времени. Движение может происходить с сохранением импульса и угловой скорости и с ускорением. Как будет падать нарисованная фигура заданная контуром? Как должна двигаться фигура в анимации, чтобы движение казалось реалистичным? Движение должно происходить относительно центра масс. Свободное движение происходит с постоянной скоростью и вращением относительно центра масс.
Центр масс плоских однородных фигур
Центр масс получается при расчете площади контура. Площадь контура рассчитывается, как сумма и разность трапеций из которых фигура состоит, масса трапеции пропорциональна ее площади. Определение центров масс трапеции: P = (P1+P2)/2 -- это центр отрезка, центр массы трапеции находится в точке Pc = (P1+P2)*(v2df){0.5,0.25}.
Для контура Pc = sum{S_i * Pc_i}/sum{S_i}, где S_i -- площадь фрагмента. При этом надо понимать, что в зависимости от направления обхода, площадь может быть отрицательной. Суммирование площадей идет с учетом знака.
Центры масс периметров однородных фигур
Это когда масса распределена по периметру, масса отрезка пропорциональна его длине. Формула та же:
Pc = sum{D_i * Pc_i}/sum{D_i}, где D_i - длина отрезка, а Pc_i = (P_i1+P2)*0,5 -- середина отрезка.
Для расчета центра масс составной (неоднородной) фигуры используется та же формула Pc = sum{M_i * Pc_i}/sum{M_i}, где M_i -- масса (весовой коэффициент) части фигуры, Pc_i -- центр масс части фигуры.
Анимация движения (трансформации)
Для каждой фигуры или группы фигур, можно определить матрицу преобразования, которая опеределяет текущее положение и поворот системы координат для прорисовки фигуры. Группу фигур можно описать относительно центра масс или относительно центра вращения. Центр вращения может не совпадать с центром масс. Например, поворот двери при открывании производится отностительно петли, а не относительно центра массы двери.
Положение фигуры на плоскости от веремени задется параметрическим уравнением -- сплайном, параметром является время движения. Ключевые точки сплайна могут быть расчитаны исходя из физики движения фигуры. Исходные данные для расчета -- масса, групповая скорость, скорость вращения, время движения, начальное и конечное положение. Результат -- уравнение движения заданное сплайном.
Более простой способ задания анимации. Задать ключевые точки, через которые проходит траектория движения и временные отсчеты соответсвующие этим точкам.
Расчет промежуточных точек траектории производится по кубическому сплайну.
Анимация контура (морфинг)
Морфинг -- плавное перетекание ключевых точек из одного контура в другой. Ограничение - число точек начального контура должно быть равно числу точек финального контура.
Преобразование выполняется с использованием сплайнов для каждой точки контура
P(t) = spline(P0, P1, t); где P0 - начальное полпожение точки, а P1- финальное положение точки контура.
Можно создать последовательнсть преобразований, задать положение контура в заданное время. Ключевым отсчетам по времени сопоставить состояние контура. В этом случае промежуточные состояния каждой точки контура могут расчитываться кубическими сплайнами.
Анимация цвета
Как правило анимация цвета -- это сплайн. Вычисление производится над четырьмя компонентами цвета: RGBA.
Для обработки и представляения цвета вводим векторный тип, состоящий из четырех значений типа float:
v4sf spline4(v4sf P0, v4sf P1, float t){
return P0 + (P1 − P0) * t;
}
Некоторые операции не могут быть эффективно описаны на языке Си, привожу примеры описания базовых операций с использованием встроенных функций и векторных расширений языка.
Я использую версию GCC-6+, но все эти наработки, встроенные функции и расширения языка будут работать начиная версии GCC-4.8+.
Для обращения к процессорным инструкциям я использую встроенные платформо-зависимые функции. Для других компиляторов (Intel или MS) следует использовать intrinsic функции. LLVM и так понимает. При переносе на другую платформу, например ARM следует переопределить базовые платформо-зависимые функции.
Тип double является избыточным для вычислений графики, достаточно float. Для платформы ARM рекомендуется такая замена. Для переопределения разрядности вычислений достаточно в одном определении заменить базовый тип и длину вектора в байтах.
Данный тип описывает точку на плоскости или радиус-вектор из начала координат. Упакованный double -- вектор, состоящий из двух элементов double описывается с использованием векторного расширения языка Си (GCC):
typedef double v2df __attribute__((__vector_size__(16)));
Базовые операции над упакованными числами
В некоторых случаях компилятор эффективно преобразует выражение вида (v2df){P[y],P[x]} в одну инструкцию SSE/AVX, а в некоторых случаях требуется явное указание инструкций. Операции унарные вида (v2df){op(P[x]), op(P[y])} и бинарные вида (v2df){op(A[x],B[x]), op(A[y],B[y])} обычно хорошо выражаются средствами языка Си. К таким операциям относятся +, −, *,/, логические операции и даже операции сравнения. Обычно не вызывает проблем смешанные операции, такие как умножение вектора на скаляр. Ниже перечислены операции над упакованными числами, которые не описываются средствами языка Си, описание требует использования обращений к инструкциям процессора.
shuffle перестановка местами компонент вектора P = (v2df){P[y],P[x]}. Всегда, не зависимо от платформы можно использовать встроенную функцию __builtin_shuffle.
shuffle2 составить вектор из элементов двух векторов, например P = (v2df){B[y],A[x]}. Для таких операций можно использовать встроенную функцию __builtin_shuffle над двумя операндами. Встроенная функция __builtin_shuffle не является платформо-зависимой, зачастую без нее не удается описать даже простые векторные операции, хотя конечно можно использовать инициализацию массива для перестановок.
addsub -- вычитание по Х компоненте и сложение по Y компоненте.
P = (v2df){A[x] − B[x], A[y] + B[y]}
__builtin_ia32_addsubpd(A,B);hadd -- горизонтальное сложение компонент вектора P[x] + P[y]
P = (v2df){A[x] + A[y], B[x] + B[y]}
__builtin_ia32_haddpd(A,B);hsub -- горизонтальное вычитание компонент вектора P[x] − P[y]
P = (v2df){A[x] − A[y], B[x] − B[y]}
__builtin_ia32_hsubpd(A,B);Скалярное произведение вектора на число P * t. Описывается средствами языка Си. Компилятор сначала преобразует число t к векторному типу T = (v2df) {t,t}, а потом производит умножение. Так же описывается сложения/вычитания и деления вектора на скаляр.
Изменение знака P * (v2df){−1,1}. В таких случаях компилятор использует умножение, а можно использовать логическую операцию XOR над знаковым битом числа:
__builtin_ia32_xorpd (P, (v2df){−0., 0});Для изменения знака используется логическая операция XOR, это приводит к появлению чисел −0.0. Следует отметить, что при сравнении чисел 0. и −0. инструкции 387 и SSE выдают истину, (0.0 == −0.0), хотя при бинарном сравнении эти числа разные. Встроенная функция используется по причине того, что компилятор не соглашается использовать бинарный оператор ^ для упакованных вещественных чисел.
max выделение максимального значения, независимо по компонентам вектора
P = (v2df){max(A[x],B[x]), max(A[y],B[y])}
__builtin_ia32_maxpd(A,B);min выделение минимального значения, независимо по компонентам вектора
P = (v2df){min(A[x],B[x]), min(A[y],B[y])}
__builtin_ia32_minpd(A,B);abs положительный вектор, независимо по компонентам вектора.
P = (v2df){abs(A[x]), abs(A[y])}
В числах с плавающей точкой знак задается старшим битом. Задача состоит в том, чтобы
сбросить знак минимальным затратами, в одну инструкцию. Это может быть битовая операция NAND, НЕ-И над знаком числа:
__builtin_ia32_andnpd( (v2df){−0., −0.} , P)sqrt - извлечение квадратного корня
P = (v2df){sqrt(A[x]), sqrt(A[y])}
__builtin_ia32_sqrtpd(A)Загрузка и выгрузка векторов без выравнивания
Векторные операции требуют выравнивания данных в памяти. Эта задача возникает при переходе от скалярного представления (массива значений типа double) к векторному типу.
Например, для загрузки значений векторов из массива
double* p;
следует выравнивать массив на 16 байт или использовать операции загрузки не выровненных данных. Для работы данными без выравнивания используются встроенные функции
__builtin_ia32_loadupd()__builtin_ia32_storeupd()или выражение
v2df a;memcpy(&a, p, sizeof(a));
которое компилятор при оптимизации заменит на инструкцию movupd.
Избыточность
тип данных остается тем же, но операцию можно производить только над младшей компонентой вектора, оставляя без изменения старшую часть вектора. Если после выполнения действий над векторами, нужно выполнить ряд действий над компонентой X или Y, то не стоит преобразовывать тип в базовый double, потому что компилятор выведет число из векторного регистра и начнет длинную вереницу инструкций плавающей точки 387, перекладывая числа через стек.
Например, нужно рассчитать корень из суммы квадратов компонент вектора:
Если эти операции выполнять на векторных регистрах SSE/AVX, то необходимо задать подряд всего три инструкции: умножение (mulpd), параллельное сложение (haddpd), и вычислить корень числа прямо на регистре SSE (sqrtsd) -- одна инструкция,- всего три инструкции.
Пример кода, после компиляции:
vmulpd %xmm1, %xmm1, %xmm1
vmulpd %xmm0, %xmm0, %xmm0
vhaddpd %xmm1, %xmm0, %xmm0
vsqrtpd %xmm0, %xmm0
__builtin_ia32_sqrtpd(__builtin_ia32_haddpd(a*a,b*b));
Без встроенных функций я не нашел способа описать эти операции так же эффективно.
Если тоже самое выполнить над компонентами вектора, по описанию:
v2df v = r*r;
double d = sqrt(v[0] + v[1]);
то получается: одна инструкция SSE - умножение (mulpd), параллельное сложение (haddpd), сохранение из регистра на стек, загрузка в из стека в сопроцессор 387 и еще десять операций 387, не считая вызова функции sqrt(). И так бывает всегда при переходе с регистров xmm на базовый тип double. Поэтому чаще всего бывает выгодно выполнять все операции на регистрах SSE, даже если это приводит к избыточным операциям, когда используется только одна компонента из вектора.
Ключи компилятора GCC
Чтобы компилятор всегда выполнял операции с типом double на регистрах SSE/AVX надо задать ключ компилятору -mfpmath=sse
Чтобы компилятор не перекладывал параметры функции типа double через стек, надо включить опцию компилятора -msseregparm или использовать соответствующий атрибут функции. Например,
v2df spline(v2df ps, v2df pe, double t) {
return ps + (pe − ps)*t;
}
при задании ключа -msseregparm параметр t передается в регистре xmm2
Ниже фрагмент кода после компиляции:
_spline:
vsubpd %xmm0, %xmm1, %xmm1
vmovddup %xmm2, %xmm2
vmulpd %xmm2, %xmm1, %xmm1
vaddpd %xmm0, %xmm1, %xmm0
ret
double spline(double ps, double pe, double t) {
return ps +(pe-ps)*t;
}
vsubsd %xmm0, %xmm1, %xmm1
vmulsd %xmm2, %xmm1, %xmm1
vaddsd %xmm0, %xmm1, %xmm0
ret
О процессорных инструкциях SSE/AVX можно прочитать в документе Intel:
Intel® 64 and IA-32 Architectures Software Developer’s Manual
Справочник с годами дополняется новыми инструкциями, доступен для скачивания на сайте производителя. Там же даны соответствующие имена intrinsic функций, которые приняты в компиляторах Intel и MSVС, и поддерживаются компиляторами GCC.
Векторные расширения GCC для языка Си:
https://gcc.gnu.org/onlinedocs/gcc/
https://gcc.gnu.org/onlinedocs/gcc/x86-Built-in-Functions.html
Если надо получить ассемблерный код пригодный для изучения, и анализа использую ключ -S.
Унарные операции над векторами
Эвклидова норма вектора, его длина
v2df norm(v2df p) { return __builtin_ia32_sqrtpd(__builtin_ia32_haddpd(p*p,p*p)); }
Результат компиляции - три инструкции SSE/AVX: умножение, горизонтальное сложение, корень квадратный.
Поворот вектора на 90 градусов по часовой стрелке или против часовой стрелки.
v2df cw(v2df p){
return __builtin_ia32_xorpd(shuffle(p),(v2df){−0.,0});
}
v2df ccw(v2df p){
return shuffle(__builtin_ia32_xorpd(p,(v2df){−0.,0}));
}
Для изменения знака используется логическая операция XOR над знаком числа.
Oтражение вектора. Горизонтальное отражение меняет знак Y-компоненты вектора. Вертикальное отражение вектора меняет знак X-компоненты вектора.
v2df reflect_vert(v2df p){
return __builtin_ia32_xorpd(p,(v2df){−0.,0});
}
v2df reflect_horz(v2df p){
return __builtin_ia32_xorpd(p,(v2df){0,−0.});
}
При редактировании векторных изображений, при копировании, часто применяются операции отражения и поворота на 90 градусов.
Бинарные операции
Расстояние между точками на плоскости
norm(p1-p0);
Dot product - скалярное произведение векторов
double dot(v2df a, v2df b){
b *= a;
return b[0] + b[1];
}
Cross Product векторное произведение, Z-координата, площадь параллелограмма (удвоенная площадь треугольника):
double cross(v2df a, v2df b){
v2df v = b*shuffle(a);
return v[0] - v[1];
}
Поворот вектора a на величину радиус-вектора r по или против часовой стрелки.
rotate_cw, rotate_ccw:
v2df rotate_ccw(v2df a, v2df r){
return addsub(a*r[0],shuffle(a)*r[1]); // a_x*r_x - a_y*r_y,
}
Расстояние от точки P до прямой. Прямая задана радиус-вектором R.
d = cross(r, p);
Знак результата определяет с какой стороны прямой расположена точка P.
Точка пересечения отрезка (P1,P2) и прямой. Прямая задана радиус-вектором R.
Сплайны
Линейная интерполяция
B_1(t) = (1−t)*P_0 + t*P_1
spline (P0,P1, t) = P0 + (P1 − P0) * t
Сплайн -- параметрическое описание кривой, удобное представление, которое можно выразить минимальным числом инструкций.
Квадратичный Bézier сплайн
{B_2}(t) = (1−t)^2 *P_0 + 2t(1−t)*P_1 + t^2 * P_2
qspline(P0,P1,P2, t) = spline(spline(P0,P1, t), spline(P1,P2, t), t)
Кубический Bézier сплайн
{B_3}(t) = (1−t)^3 *P_0 + 3t(1−t)^2 * P_1 + 3(1−t)t^2 * P_2 + t^3 * P_3
cspline(P0,P1,P2,P3, t) = qspline(spline(P0,P1, t), spline(P1,P2, t), t), spline(P2,P3, t), t)
Функции сплайнов можно определять рекурсивно, как показано выше. Компилятор из такого описания делает хорошо оптимизированный код.
Касательная к кубическому сплайну, используется для движения по траектории заданной сплайном, производная от кубического сплайна - квадратичный сплайн
{B_3}'(t) = 3(1−t)^2 ({P}_1 − {P}_0) + 6(1−t)t({P}_2 − {P}_1) + 3t^2({P}_3 − {P}_2)
{B_3}'(t) = 3* qspline(P1−P0, P2−P1, P3−P2, t)
Касательная к траектории заданной кубическим сплайном -- скорость движения по траектории.
{B_2}'(t) = 2(1-t) ({P}_1 − {P}_0) +2t({P}_2 − {P}_1)
{B_2}'(t) = 2* spline(P1−P0, P2−P1, t)
Для комплекта
{B_1}'(t) = P1−P0.
Общее утверждение -- производная от кубического сплайна - крвадратичный сплай, производная от квадратичного сплайна - линейный сплайн, производная от сплайна - направление отрезка. Физический смысл производной -- скорость движения по контуру, направление вектора скорости -- касательная к траектории. Направление производной используется для построения фигур вдоль траектории.
Описание дуги окружности кубическим Bézier-сплайном
k=4/3(sqrt(2)-1), k = 0.5522847498
https://people.mozilla.org/~jmuizelaar/Riskus354.pdf
Если угол (a) меньше 90 то k*tg(a/2)
ax = x1 – xc
ay = y1 – yc
bx = x4 – xc
by = y4 – yc
q1 = ax*ax + ay*ay
q2 = q1 + ax*bx + ay*by
k2 = 4/3(sqrt(2*q1*q2)-q2)/(ax*by - bx*ay)
x2 =xc + x1 – k2*y1
y2 =yc + y1 + k2*x1
x3 =xc + x4 – k2*y4
y3 =yc + y4 + k2*x4
Преобразование линейных сплайнов в квадратичные
{B_2}(t) = (1−t)^2 *P_0 + 2t(1−t)*P_1 + t^2 * P_2
P1 = (P0+P2)/2
Преобразование линейных сплайнов в кубические
{B_3}(t) = (1−t)^3 *P_0 + 3t(1−t)^2 * P_1 + 3(1−t)t^2 * P_2 + t^3 * P_3
P1=P0+1/3*(P3−P0);
P2=P3−1/3*(P3−P0);
Преобразование квадратичных сплайнов в кубические
Дело в том, что внутреннее представление сплайнов может быть единообразным, все кривые можно свести к кубическим сплайнам. Квадратичный сплайн можно представить в форме кубического:
P1=P0+2/3*(P1s−P0);
P2=P3+2/3*(P1s−P3);
здесь P1s - ключевая точка квадратичного сплайна.
Операции с матрицами трансформации
x_new = xx * x + xy * y + tx;
y_new = yx * x + yy * y + ty;
Из коэффициентов преобразования можно составить матрицу {xx, yx, xy, yy, tx, ty}. В упакованном типе запишется так:
struct {
v2df x,y,z;
};
где x = (v2df){xx,yx}, y = (v2df){xy,yy}, z = (v2df){tx,ty};
Трансформация складывается из поворотов, масштабирования и смещения начала отсчета.
Выделяют несколько операций, из которых складывается трансформация: поворотов, масштабирования и смещения начала отсчета.
Поворот можно задать углом, или двумя параметрами sin и cos угла, вместе они образуют (v2df){cos(a), sin(a)} радиус-вектор, наклон которого определяет угол поворота.
Единичные матрицы трансформации
Умножение матриц трансформации
Смещение начала координат
Вращение координат
Масштабирование
Обратная матрица трансформации
Операции с контурами
Контур на плоскости задается точками: прямыми и сплайнами. в общем виде используются линейные и кубические Bézier сплайны.
Принадлежность точки замкнутому контуру.
Площадь контура
Сумма площадей трапеций:
S = sum{S_i}, площадь трапеции S_i =1/2 (P_{i}[x] − P_{i-1}[x]) * (P_{i}[y] + P_{i-1}[y]})
Другой способ:
Можно суммировать и вычитать треугольники относительно начала отчета. Треугольники задаются радиус-векторами точек контура. Площадь треугольника - cross product - векторное произведение векторов-сторон треугольника.
cross(a,b) = a[0]*b[1] − a[1]*b[0];
Выражается тремя операциями
v2df v = a*shuffle(b);
return v[0]-v[1];
Значение векторного произведения равно площади параллелограмма, двум площадям треугольника A0B, а знак будет описывать направление обхода контура. Берем любую точку вблизи фигуры, назначаем ее центром отсчета и относительно нее суммируем площади треугольников по контуру.
При вычислении суммы площадей можно суммировать вектора по контуру и в конце выполнить вычитание сумм, вместо того чтобы выполнять вычитание на каждом шаге цикла:
for (i=0;i<num;i++) { sum += a*shuffle(b); }
return sum[0] - sum[1];
Вычисление площади контура можно еще ускорить, если использовать тождество:
[a,b] + [b,c] = [(a − c),b]
квадратными скобками обозначено векторное произведение.
см. Алгебраические свойства векторного произведения
тогда в цикле получаем операцию
sum += (a − c)*shuffle(b); -- четыре операции на два треугольника.
Длина контура
int_{0}^{1} { norm({B_3}'(t)) dt}
Физика
В общем виде физика движения объекта описывается кубическими сплайнами. Координаты центра масс задаются в зависимости от времени. Движение может происходить с сохранением импульса и угловой скорости и с ускорением. Как будет падать нарисованная фигура заданная контуром? Как должна двигаться фигура в анимации, чтобы движение казалось реалистичным? Движение должно происходить относительно центра масс. Свободное движение происходит с постоянной скоростью и вращением относительно центра масс.
Центр масс плоских однородных фигур
Центр масс получается при расчете площади контура. Площадь контура рассчитывается, как сумма и разность трапеций из которых фигура состоит, масса трапеции пропорциональна ее площади. Определение центров масс трапеции: P = (P1+P2)/2 -- это центр отрезка, центр массы трапеции находится в точке Pc = (P1+P2)*(v2df){0.5,0.25}.
Для контура Pc = sum{S_i * Pc_i}/sum{S_i}, где S_i -- площадь фрагмента. При этом надо понимать, что в зависимости от направления обхода, площадь может быть отрицательной. Суммирование площадей идет с учетом знака.
Центры масс периметров однородных фигур
Это когда масса распределена по периметру, масса отрезка пропорциональна его длине. Формула та же:
Pc = sum{D_i * Pc_i}/sum{D_i}, где D_i - длина отрезка, а Pc_i = (P_i1+P2)*0,5 -- середина отрезка.
Для расчета центра масс составной (неоднородной) фигуры используется та же формула Pc = sum{M_i * Pc_i}/sum{M_i}, где M_i -- масса (весовой коэффициент) части фигуры, Pc_i -- центр масс части фигуры.
Анимация движения (трансформации)
Для каждой фигуры или группы фигур, можно определить матрицу преобразования, которая опеределяет текущее положение и поворот системы координат для прорисовки фигуры. Группу фигур можно описать относительно центра масс или относительно центра вращения. Центр вращения может не совпадать с центром масс. Например, поворот двери при открывании производится отностительно петли, а не относительно центра массы двери.
Положение фигуры на плоскости от веремени задется параметрическим уравнением -- сплайном, параметром является время движения. Ключевые точки сплайна могут быть расчитаны исходя из физики движения фигуры. Исходные данные для расчета -- масса, групповая скорость, скорость вращения, время движения, начальное и конечное положение. Результат -- уравнение движения заданное сплайном.
Более простой способ задания анимации. Задать ключевые точки, через которые проходит траектория движения и временные отсчеты соответсвующие этим точкам.
Расчет промежуточных точек траектории производится по кубическому сплайну.
Анимация контура (морфинг)
Морфинг -- плавное перетекание ключевых точек из одного контура в другой. Ограничение - число точек начального контура должно быть равно числу точек финального контура.
Преобразование выполняется с использованием сплайнов для каждой точки контура
P(t) = spline(P0, P1, t); где P0 - начальное полпожение точки, а P1- финальное положение точки контура.
Можно создать последовательнсть преобразований, задать положение контура в заданное время. Ключевым отсчетам по времени сопоставить состояние контура. В этом случае промежуточные состояния каждой точки контура могут расчитываться кубическими сплайнами.
Анимация цвета
Как правило анимация цвета -- это сплайн. Вычисление производится над четырьмя компонентами цвета: RGBA.
Для обработки и представляения цвета вводим векторный тип, состоящий из четырех значений типа float:
typedef float v4sf __attribute__((__vector_size__(16)));
v4sf spline4(v4sf P0, v4sf P1, float t){
return P0 + (P1 − P0) * t;
}
Комментариев нет:
Отправить комментарий