Масштабирование значений байтовых пикселей (y=ax+b) с SSE2 (как с плавающей точкой)?

Question

Я хочу посчитать y = ax + bгде x и y - значение пикселя [т.е. байт с диапазоном значений составляет 0~255], тогда как a а также b это поплавок

Поскольку мне нужно применять эту формулу для каждого пикселя в изображении, кроме того, a и b различны для разных пикселей. Прямые вычисления в C++ медленны, поэтому мне интересно узнать инструкцию sse2 в C++.

После поиска я обнаружил, что умножение и сложение в float с sse2 так же, как _mm_mul_ps а также _mm_add_ps, Но в первую очередь мне нужно преобразовать x в байтах в float (4 байта).

Вопрос в том, что после загрузки данных из источника байтовых данных (_mm_load_si128), как я могу преобразовать данные из байта в число с плавающей точкой?

Answer 1

a а также b различны для каждого пикселя? Это усложнит векторизацию, если нет паттернов или вы не можете их сгенерировать

Есть ли способ, которым вы можете эффективно генерировать a а также b в векторах, как с фиксированной или с плавающей точкой? В противном случае вставка 4 значений FP или 8 16-битных целых может быть хуже, чем просто скалярные операции.

---

Фиксированная точка

Если a а также b может быть использован повторно вообще или сгенерирован с фиксированной точкой, это может быть хорошим вариантом для математики с фиксированной точкой. (т.е. целые числа, представляющие значение * 2^ шкала). SSE/AVX не имеют умножения 8b*8b->16b; наименьшими элементами являются слова, поэтому вы должны распаковывать байты в слова, но не до 32 бит. Это означает, что вы можете обрабатывать вдвое больше данных за инструкцию.

Есть _mm_maddubs_epi16 инструкция, которая может быть полезна, если b а также a изменить достаточно редко, или вы можете легко генерировать вектор с чередованием байтов *2^4 и b*2^1. По-видимому, это действительно удобно для билинейной интерполяции, но она по-прежнему выполняет работу за нас с минимальным перемешиванием, если мы можем подготовить вектор a и b.

float a, b;
const int logascale = 4, logbscale=1;
const int ascale = 1<<logascale;  // fixed point scale for a: 2^4
const int bscale = 1<<logbscale;  // fixed point scale for b: 2^1

const __m128i brescale = _mm_set1_epi8(1<<(logascale-logbscale));  // re-scale b to match a in the 16bit temporary result

for (i=0 ; i<n; i+=16) {
    //__m128i avec = get_scaled_a(i);  
    //__m128i bvec = get_scaled_b(i);
    //__m128i ab_lo = _mm_unpacklo_epi8(avec, bvec);
    //__m128i ab_hi = _mm_unpackhi_epi8(avec, bvec);

    __m128i abvec = _mm_set1_epi16( ((int8_t)(bscale*b) << 8) | (int8_t)(ascale*a) );  // integer promotion rules might do sign-extension in the wrong place here, so check this if you actually write it this way.

    __m128i block = _mm_load_si128(&buf[i]);  // call this { v[0] .. v[15] }

    __m128i lo = _mm_unpacklo_epi8(block, brescale);  // {v[0], 8, v[1], 8, ...}
    __m128i hi = _mm_unpackhi_epi8(block, brescale);  // {v[8], 8, v[9], 8, ...
    lo = _mm_maddubs_epi16(lo, abvec);  // first arg is unsigned bytes, 2nd arg is signed bytes
    hi = _mm_maddubs_epi16(hi, abvec);
    // lo = { v[0]*(2^4*a) + 8*(2^1*b), ... }

    lo = _mm_srli_epi16(lo, logascale);  // truncate from scaled fixed-point to integer
    hi = _mm_srli_epi16(hi, logascale);

    // and re-pack.  Logical, not arithmetic right shift means sign bits can't be set
    block = _mm_packuswb(lo, hi);  
    _mm_store_si128(&buf[i], block);
}
// then a scalar cleanup loop

2 ^ 4 - произвольный выбор. Оставляет 3 незнаковых бита для целой части aи 4 дробных бита. Так что это эффективно округляет a до ближайшего 16-го числа и переполняется, если оно имеет величину больше 8 и 15/16-го числа. 2^6 даст больше дробных битов и позволит a от -2 до +1 и 63/64.

поскольку b добавляется, а не умножается, его полезный диапазон намного больше, а дробная часть гораздо менее полезна. Чтобы представить его в 8 битах, округление до ближайшей половины все еще сохраняет немного дробной информации, но позволяет ему быть [-64: 63.5] без переполнения.

Для большей точности 16b с фиксированной точкой - хороший выбор. Вы можете масштабировать a а также b до 2^7 или около того, чтобы иметь 7b дробной точности и при этом разрешить целочисленную часть быть [-256 .. 255]. Для этого случая нет инструкции умножения и добавления, так что вам придется делать это отдельно. Хорошие варианты для умножения включают в себя:

• _mm_mulhi_epu16: unsigned 16b*16b->high16 (биты [31:16]). Полезно, если a не может быть отрицательным
• _mm_mulhi_epi16: подписано 16b*16b->high16 (биты [31:16]).
• _mm_mulhrs_epi16: знаковые 16b * 16b-> биты [30:15] 32b временных, с округлением. С хорошим выбором коэффициента масштабирования для aэто должно быть лучше. Насколько я понимаю, SSSE3 ввел эту инструкцию именно для этого вида использования.
• _mm_mullo_epi16: подписано 16b*16b->low16 (биты [15:0]). Это позволяет только 8 значащих битов для a до того, как результат low16 переполнится, поэтому я думаю, что все, что вы получите за _mm_maddubs_epi16 8-битное решение более точно для b,

Чтобы использовать их, вы получите масштабированные 16b векторов a а также b значения, то:

• распакуйте ваши байты с нуля (или pmovzx byte->word), чтобы получить слова со знаком в диапазоне [0..255]
• сдвиг влево слова на 7.
• умножить на ваш a вектор из 16b слов, беря верхнюю половину каждого результата 16*16->32. (например, муль
• сдвиг вправо здесь, если вы хотели разные шкалы для a а также b, чтобы получить более дробную точность для a
• добавлять b к этому.
• сдвиг вправо, чтобы выполнить окончательное усечение назад от фиксированной точки до [0..255].

При хорошем выборе шкалы с фиксированной запятой это должно быть в состоянии обрабатывать более широкий диапазон a а также b, а также с большей дробной точностью, чем 8-битная фиксированная точка.

Если вы не сдвинете свои байты влево после распаковки их в слова, a должен быть полнодиапазонным, чтобы получить 8 битов в максимуме 16 результата. Это будет означать очень ограниченный диапазон a что вы могли бы поддерживать без усечения вашего временного до менее чем 8 бит во время умножения. Четное _mm_mulhrs_epi16 не оставляет много места, так как начинается с 30-го.

---

расширить байты до числа с плавающей запятой

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

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

Это хороший пример использования pmovzx (SSE4.1), потому что он может напрямую переходить с 8b элементов на 32b. Другие варианты SSE2 punpck[l/h]bw/punpck[l/h]wd с несколькими шагами или SSSE3 pshufb подражать pmovzx, (Вы можете сделать одну загрузку 16B и перемешать ее 4 различными способами, чтобы распаковать ее в четыре вектора по 32b.)

char *buf;

// const __m128i zero = _mm_setzero_si128();
for (i=0 ; i<n; i+=16) {
    __m128 a = get_a(i);
    __m128 b = get_b(i);

    // IDK why there isn't an intrinsic for using `pmovzx` as a load, because it takes a m32 or m64 operand, not m128.  (unlike punpck*)
    __m128i unsigned_dwords = _mm_cvtepu8_epi32((__m128i)(buf+i));  // load 4B at once.
    __m128 floats = _mm_cvtepi32_ps(unsigned_dwords);

    floats = _mm_fmadd_ps(floats, a, b);  // with FMA available, this might as well be 256b vectors, even with the inconvenience of the different lane-crossing semantics of pmovzx vs. punpck
    // or without FMA, do this with _mm_mul_ps and _mm_add_ps
    unsigned_dwords = _mm_cvtps_epi32(floats);

    // repeat 3 more times for buf+4, buf+8, and buf+12, then:

    __m128i packed01 = _mm_packss_epi32(dwords0, dwords1);  // SSE2
    __m128i packed23 = _mm_packss_epi32(dwords2, dwords3);
    // packuswb wants SIGNED input, so do signed saturation on the first step

    // saturate into [0..255] range
    __m12i8 packedbytes=_mm_packus_epi16(packed01, packed23);  // SSE2
    _mm_store_si128(buf+i, packedbytes);  // or storeu if buf isn't aligned.
}

// cleanup code to handle the odd up-to-15 leftover bytes, if n%16 != 0

Предыдущая версия этого ответа была взята из векторов float->uint8 с packusdw/packuswb и имела целый раздел об обходных путях без SSE4.1. Ни один из этих битов маскировки знака после неподписанного пакета не требуется, если вы просто остаетесь в подписанном целочисленном домене до последнего пакета. Я предполагаю, что это причина, по которой SSE2 включал только подписанный пакет от слова к слову, но и подписанный, и неподписанный пакет от слова к байту. packuswd полезно только если ваша конечная цель uint16_t, а не дальнейшая упаковка.

---

Последним процессором, который не имел SSE4.1, были Intel Conroe/merom (первый поколения Core2, до конца 2007 года) и AMD до Барселоны (до конца 2007 года). Если для этих процессоров приемлемо работать, но медленно, просто напишите версию для AVX2 и версию для SSE4.1. Или SSSE3 (с 4x pshufb для эмуляции pmovzxbd из четырех элементов 32b регистра) pshufb работает медленно в Conroe, поэтому, если вы заботитесь о процессорах без SSE4.1, напишите конкретную версию. На самом деле, Конро / Мером также имеет медленный хмм punpcklbw и так далее (за исключением q->dq). 4х медленный pshufb должен все еще бить 6x медленные распаковки. Векторизация - намного меньшая победа перед Вольфдейлом из-за медленных перетасовок при распаковке и переупаковке. Версия с фиксированной запятой, с гораздо меньшим количеством распаковки / переупаковки, будет иметь еще большее преимущество.

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

Answer 2

Я думаю, вы ищете __m128 _mm_cvtpi8_ps(__m64 a ) сложный внутренний.

Вот минимальный пример:

#include <xmmintrin.h>
#include <stdio.h>

int main() {
  unsigned char  a[4] __attribute__((aligned(32)))= {1,2,3,4};
  float b[4] __attribute__((aligned(32)));
  _mm_store_ps(b, _mm_cvtpi8_ps(*(__m64*)a));
  printf("%f %f, %f, %f\n", b[0], b[1], b[2], b[3]);
  return 0;
}