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

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

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

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

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

Ответы

Ответ 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/64ths.

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

Для большей точности 16-битная фиксированная точка является хорошим выбором. Вы можете масштабировать 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. (например, mul
  • сдвиг вправо здесь, если вы хотите разные шкалы для 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 ints.)

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 включил только подписанный пакет от слова dword к слову, но как подписанный, так и неподписанный пакет от слова к байту. packuswd полезен, если ваша конечная цель uint16_t, а не дальнейшая упаковка.


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

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

Ответ 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;
}