SIMD pixel-contrast: сумма различий между пикселем и его соседями (компоненты цвета uint16_t, суммы с плавающей запятой)?

Question

Каков наилучший способ использования SIMD/ассемблера для вычитания 2 uint16 с абсолютным значением (максимальная разница) и добавления (+ =) результата к поплавке?

Подобно этому примеру C'ish

c0 += fabs((float)a0 - (float)b0); // C is Float accumulator, a+b pixels 

, где а и Ь беззнаковое 16-битные слова, и с представляет собой поплавок. Только одно слово -> преобразование с плавающей точкой, а не 3.

Приложение Thee обрабатывает необработанные 16-разрядные данные без знака в виде как можно большего количества RGB-пикселей.

Возможно использование AVX2/SSE4.2 на Skylake Xeon E3-1275 v5?

---

5 минут ограничение комментариев ?? Не удается сохранить или переименовать ???

Вы уверены, что вам нужно поплавать? Uint16 не может накапливать более 1 вычитания. Я хочу сделать подсчет окрестности, поэтому мне нужно суммировать как минимум 8 различий. В окрестности с глубиной D имеются (2D + 1)^2-1 соседи. Я также хочу иметь возможность разделить разницу, в которой uint32 может быть слишком маленьким. Я думаю, что поплавки выглядят более гладкими.

Вот немного больше фона на то, что уже работает, и как я хочу, чтобы улучшить его.

Чтобы уточнить, мой текущий код C вычисляет разницу между каналами между фиксированным домашним пикселем и 8 или более соседями. Он имеет структуру с 5 глубокими вложенными циклами: Y-строки, тогда X-cols для каждого пикселя на изображении (36 миллионов) Каналы, R. G & B - loop3 Петли 4 и 5 для строк и столбцов окрестности.

для каждого HOME пикселя ясно, что R, G и B аккумуляторов для каждого соседа,
добавить абса (home_red - nabr_red) до red_float_accumulator же для зеленых и синих копии накопленных значений в основную память

Следующим шагом было перемещение каналов на уровень 5 и все 3 вычитания, R, G и B одновременно с SIMD. Имея 48 бит/пиксель и 128 бит, доступных в каждом регистре MMX, 2 можно сделать сразу, а не только 1.

С 512-битными регистрами в AVX2 на Skylake Xeon, 10 может быть выполнено. Я ищу хорошую стратегию, чтобы сбалансировать сложность с производительностью и узнать больше об этих векторных операциях.

Мне нужны аккумуляторы R, G и B для каждого «домашнего» пикселя. Затем переместите RGB в «плавающее изображение» с тем же разрешением XY, что и uint16/channel RAW, RGB-файл. Сделайте то же сравнение контрастности для каждого пикселя.

Answer 1

От Justin W's answer на вопрос о беззнаковом целочисленном абсолютном значении: повторное насыщение вычитания. Один результат будет равен нулю, другой будет абсолютным значением. Объедините их с логическим ИЛИ. Делать это таким образом, перед распаковкой 16b ints до 32b ints или float, является самым дешевым.

Мы обязательно хотим вычесть перед распаковкой от слова к слову, поэтому есть только одно значение для распаковки. Выполнение одного не насыщающего вычитания и последующая его сортировка позже (например, путем маскировки знакового бита до нуля после преобразования в FP) не будут работать. Проблема с диапазоном: подписанный int16_t не может содержать полный диапазон разницы между двумя значениями uint16_t. UINT16_MAX обернется и будет выглядеть как разница -1 -> абс 1. Кроме того, у него будет недостаток в необходимости распаковки знака.

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

vpunpcklwd распаковывается в пределах каждой полосы 128b, соответствующей vpackusdw. Для vpmovzxwd ymm, xmm нет одноинтерпрессивной инструкции, поэтому я буду использовать инструкции punpck, исходя из предположения, что вы можете упорядочить поплавки таким образом. (И с PMOVZX, вы не можете перейти непосредственно от верхней половины. Вы должны были бы vextracti128/vpmovzx.)

#include <immintrin.h> 

// untested 
__m256 add_pixdiff(__m256 c[2], __m256i a, __m256i b) 
{ 
    __m256i ab = _mm256_subs_epu16(a, b); // 0  or abs(a-b) 
    __m256i ba = _mm256_subs_epu16(b, a); // abs(a-b) or 0 
    __m256i abs_ab_diffs = _mm256_or_si256(ab, ba); 

    __m256i lo_uints = _mm256_unpacklo_epi16(abs_ab_diffs, _mm256_setzero_si256()); 
    __m256i hi_uints = _mm256_unpackhi_epi16(abs_ab_diffs, _mm256_setzero_si256()); 

    __m256 lo_floats = _mm256_cvtepi32_ps(lo_uints); 
    __m256 hi_floats = _mm256_cvtepi32_ps(hi_uints); 

    // use load and store intrinsics if the data might not be aligned. 
    c[0] = _mm256_add_ps(c[0], lo_floats); 
    c[1] = _mm256_add_ps(c[1], hi_floats); 
    return c[0]; 
} 

Это compiles exactly like you'd expect, на godbolt. Для использования в цикле, вероятно, лучше всего вручную установить его. Вам не нужен глупый компилятор на самом деле с использованием массива в памяти для выполнения вызова по ссылке для двух векторов float. Я просто завернул его в функцию, чтобы сохранить его простым и скрыть загрузку/хранилища.

Обратите внимание, что один вектор ввода uint16 создает два вектора поплавка. Мы могли бы работать с целыми векторами 128b за раз, но делать 256b за один раз означает, что мы получаем 2 результата для 3 insns (не считая нагрузок) + 2 unpacks, а не 6 insns + 2 pmovzx. Здесь есть приличное количество параллелизма: две вычитания могут происходить параллельно, и есть две цепи зависимостей для распаковки и конвертирования. (У Skylake есть только один порт в случайном порядке, так что вы не можете сразу получить две распаковки. Это инструкция в полете with a latency of 1c, против vpmovzx ymm, xmm, имеющая латентность 3, как и другие инструкции по переходу на переулок. , если вам нужно сохранить поплавки в том же порядке, что и int, а не просто повторно упаковывать обратно в тот же порядок в конце.)

Задержка латентности FP add увеличена до 4 на Skylake (от 3 на предыдущем Intel), но с увеличением пропускной способности до двух за такт. (Он запускает их на блоке FMA. Правильно, FMA Skylake сокращается до 4 с).

---

Я предполагаю, что вы на самом деле не хотите SAD (один аккумулятор для всех различий), так как вы написали c0, не c в вашей скалярной форме.

Если вы do хотите SAD (сумма абсолютных различий), это намного проще, и вы должны аккумулироваться в целочисленном домене. Используйте один и тот же sub в обоих направлениях с беззнаковым трюком насыщения, ИЛИ вместе, распакуйте до 32 бит int, а затем добавьте в векторный аккумулятор. Сделайте горизонтальную сумму в конце.

Answer 2

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

Трюк собирается выбрать правильный цикл в ваших 5-вложенных циклах для векторизации.

// probably convert the image to planar outside the loop 
// rather than dealing with packed components on the fly, 
// since you read the same src pixel many times 
// (comparing it to a different dst pixel each time). 

for (dst_row) { 
    for (dst_col) { 
     for (r, g, b) { 
      // if you convert to planar, make this the outer-most loop for cache reasons 
      for (src_row) { 
       for (src_col) { 
        accumulate stuff; 
       } 
      } 
     } 
     result[drow][dcol].r = sum over some range of red src rows X cols; 
     result[drow][dcol].g = sum over some range of green src rows X cols; 
     result[drow][dcol].b = sum over some range of blue src rows X cols; 
    } 
} 

Вы можете подумать, вы должны просто векторизации accumulate stuff для одного пикселя или одного компонента одного пикселя. Для этого потребуется сделать горизонтальную сумму (чтобы уменьшить все элементы векторного аккумулятора до одной скалярной суммы).

Лучше, вероятно, было бы аккумулировать результаты для 16 столбцов назначения одновременно (AVX/AVX2: два вектора 256b из 8 упакованных одноточечных поплавков, так как вы получаете два вектора, имеющих ценность от распаковки одного вектора из 16b данных). FP, вероятно, не намного хуже, чем 32-битное целое число, и, возможно, лучше, если вы используете инструкции FMA.

Каждая итерация внутренних двух петель, вы накапливаете результаты для 16 различных пикселей src. (или компоненты src, если вы не конвертируете в плоскость вне цикла.) Это создает лот местоположения во внутренних циклах. Вы сдвигаете окно 256b по 16b за раз (или по 48b за упакованный вместо планарного), а не 256b за раз, поэтому данные многократно используются.

Лучше всего, вам никогда не придется выполнять горизонтальные операции вообще. На каждой итерации цикла for (dst_col) ваши векторы результатов заканчиваются сохранением 16 элементов полных результатов, готовых к хранению прямо в result[drow][dcol..dcol+16]. (Хотя это на этой точке, где вы могли бы использовать vpermd перетасовать по переулкам и выправить порядок, что vpunpcklwd дает, по сравнению с vpmovzxwd. Если вы заботитесь о том, ваш массив с плавающей точкой хранятся в строго правильном порядке.)

Распаковка к planar будет делать угловые случаи намного проще, и сохранение ваших результатов с поплавком в плоском формате структурных массивов облегчает векторизация чего-то вроде суммирования контраста r, g и b для каждого пикселя. addps намного быстрее, чем haddps, поэтому вы хотите получить 8 r значений в одном векторе и 8 g значений в другом векторе, а не в двух векторах {rgb, rgb, x} с двумя потерянными элементами или чем-то еще.