Оператор Собеля - SIMD реализация. SSE/AVX

Question

Основная задача - перевести картинку в байты от 0 до 255 и обработать оператором Собеля с использованием SIMD, версия без симд перед вами, не знаю, как начать и использованием симд.

using System;
using System.Drawing;
using System.IO;
using System.Runtime.InteropServices;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;


var img = Image.FromFile(@"..\..\..\tmp.png");
var bmp = new Bitmap(img);

byte[,] pixels = new byte[bmp.Width, bmp.Height];
byte[,] sobeled = new byte[bmp.Width, bmp.Height];

for (int i = 0; i < pixels.GetLength(0); i++)
{
    for (int j = 0; j < pixels.GetLength(1); j++)
    {
        pixels[i, j] = bmp.GetPixel(i, j).G;
    }
}

float[,] gx = new float[bmp.Width, bmp.Height];
float[,] gy = new float[bmp.Width, bmp.Height];

for (int i = 1; i < pixels.GetLength(0) - 1; i++)
{
    for (int j = 1; j < pixels.GetLength(1) - 1; j++)
    {
        gx[i, j] = (pixels[i - 1, j + 1] + 2 * pixels[i, j + 1] + pixels[i + 1, j + 1]) - (
                    pixels[i - 1, j - 1] + 2 * pixels[i, j - 1] + pixels[i + 1, j - 1]);
        gy[i, j] = (pixels[i + 1, j - 1] + 2 * pixels[i + 1, j] + pixels[i + 1, j + 1]) - (
                   pixels[i - 1, j - 1] + 2 * pixels[i - 1, j] + pixels[i - 1, j + 1]);

        float sobeled_pixel = MathF.Sqrt(gx[i, j] * gx[i, j] + gy[i, j] * gy[i, j]);
        sobeled_pixel = sobeled_pixel > 255 ? 255 : sobeled_pixel;

        byte byte_pixel = Convert.ToByte(sobeled_pixel);
        sobeled[i, j] = byte_pixel;
    }
}

Bitmap bitmapSobeled = new Bitmap(img);
for (int i = 0; i < pixels.GetLength(0); i++)
{
    for (int j = 0; j < pixels.GetLength(1); j++)
    {
        Color color = Color.FromArgb(sobeled[i, j], sobeled[i, j], sobeled[i, j]);
        bitmapSobeled.SetPixel(i, j, color);
    }
}
bitmapSobeled.Save("res.bmp");

Попытка использовать SIMD

var const_one =  Vector128.Create(1.0f);
var const_two =  Vector128.Create(2.0f);
var const_n_one=  Vector128.Create(-1.0f);
var const_n_two=  Vector128.Create(-2.0f);

unsafe
{
    for (int i = 1; i < pixels.GetLength(0) - 1; i++)
    {
        for (int j = 1; j < pixels.GetLength(1) - 1; j+=4)
        {
            fixed (byte* p = &pixels[i, j])
            {
                var bytes = Sse2.LoadVector128(p);// load 16 x 8 bit values
                //var tmp2 = Sse2.UnpackLow(ints, Sse2.Xor(ints, ints));// unpack 8 bit -> 16 bit
                //var tmp3 = Sse2.UnpackLow(tmp2, Sse2.Xor(tmp2, tmp2));// unpack 16 bit -> 32 bit
                //var tmp4 = Sse2.ConvertToVector128Single(tmp3);
                var ints16 = Sse41.ConvertToVector128Int16(bytes);// unpack 8 bit -> 16 bit
                var ints32 = Sse41.ConvertToVector128Int32(ints16); //unpack 16 bit-> 32 bit
                var floats32 = Sse2.ConvertToVector128Single(ints32); //int 32 bit-> 32 bit float
            }
        }
    }
}

Answer 1

Ох, это было не просто.

Документации по SIMD инструкциям почти нет. То есть, в векторизацию в основном могут только те, кто хорошо знает C++ или ассемблер и имеет достаточно опыта. Для C# это штука сильно новая.

Я собрал все свои знания английского, ассемблера и черный пояс по гуглу, и пошел искать, так как разобраться было очень интересно. Из самого полезного нашел вот эти ресурсы:

The Intel Intrinsics Guide - интерактивный справочник по интринсикам от Intel. Можно из описания C# метода копировать сигнатуру C++ инструкции и быстро найти ее описание в этом справочнике. Уникальной особенностью данного справочника является информация о производительности каждой инструкции в процессорных тактах, отдельно для каждого типа ядра процессора Intel.
Introduction to Core 3 Intrinsics in C#, with Benchmarks - здесь отличная вводная в интринсики для тех, кто не умеет, но очень хочет.
x86/x64 SIMD Instruction List (SSE to AVX512) - хоть и слегка не полный, но очень хороший справочник по инструкциям. Там в таблицах даже шарповые названия методов есть. Если открыть метод, то можно увидеть схемы, как там байты ходят в регистрах, и как устроена логика методов.
Hardware Intrinsics in .NET Core - вводная в интринсики от Microsoft.

Далее я принялся за простое - десатурировать картинку с использованием векторов, то есть переписать вот это.

for (int i = 0; i < pixelsCount; i++)
{
    int offset = i * 4;
    byte brightness = (byte)(buffer[offset] * 0.118f + buffer[offset + 1] * 0.58661f + buffer[offset + 2] * 0.29891f);
    pixelBuffer[i] = brightness;
}

Вот так :)

private static unsafe void ToGrayscaleVectorized(byte* srcPtr, byte* dstPtr, int pixelsCount)
{
    byte* tail = srcPtr + (pixelsCount & -16) * 4;
    byte* srcEnd = srcPtr + pixelsCount * 4;
    byte* dstEnd = dstPtr + pixelsCount;
    sbyte r = (sbyte)(128 * 0.29891f);
    sbyte g = (sbyte)(128 * 0.58661f);
    sbyte b = (sbyte)(128 * 0.118f);
    Vector256<sbyte> weight = Vector256.Create(b, g, r, 0, b, g, r, 0, b, g, r, 0, b, g, r, 0, b, g, r, 0, b, g, r, 0, b, g, r, 0, b, g, r, 0);
    Vector256<short> one = Vector256.Create((short)1);
    while (true)
    {
        while (srcPtr < tail)
        {
            Vector256<short> vs0 = Avx2.MultiplyAddAdjacent(Avx.LoadVector256(srcPtr), weight);
            Vector256<short> vs1 = Avx2.MultiplyAddAdjacent(Avx.LoadVector256(srcPtr + 32), weight);
            Vector256<int> vi0 = Avx2.MultiplyAddAdjacent(vs0, one);
            Vector256<int> vi1 = Avx2.MultiplyAddAdjacent(vs1, one);
            Vector128<short> v0 = Sse2.ShiftRightLogical(Sse2.PackSignedSaturate(Avx2.ExtractVector128(vi0, 0), Avx2.ExtractVector128(vi0, 1)), 7);
            Vector128<short> v1 = Sse2.ShiftRightLogical(Sse2.PackSignedSaturate(Avx2.ExtractVector128(vi1, 0), Avx2.ExtractVector128(vi1, 1)), 7);
            Sse2.Store(dstPtr, Sse2.PackUnsignedSaturate(v0, v1));
            srcPtr += 64;
            dstPtr += 16;
        }
        if (srcPtr == srcEnd)
            break;
        tail = srcEnd;
        srcPtr = srcEnd - 64;
        dstPtr = dstEnd - 16;
    }
}

Основная идея векторизации кода в том, чтобы как бы выполнять несколько итераций цикла за раз. Например если вектор вмещает в себя 8 значений, то объединив вычисления в векторах, можно выполнить цикл в 8 раз меньше раз.

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

Основная ошибка - это пытаться перемешать входные данные, типа для того чтобы их можно было засунуть в вектора последовательным чтением из памяти. Опытный ум программиста сразу сообразит, что такое делать неэффективно, и сразу придет к вот такому выводу, как в соседнем ответе у @Alexander Petrov:

Применить SIMD-операции весьма проблематично. Ведь нужен доступ не к последовательным элементам, а расположенным в разных колонках и строках.

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

Вторая вероятная ошибка векторизации кода - неполный вектор. Например вам надо за итерацию цикла обработать 3 значения, а в векторе у вас 4 места, 1 остается пустым. Это значит, что вы используете мощь вектора только на 75%, и стоит пересмотреть логику метода. Избегайте таких ситуаций, наиболее вероятно, что мысль при написании кода пошла не туда. Обратитесь к рекомендациям, описанным мной выше.

Пока переписывал, немного освоился, и принялся за главное, переписать вот это.

private static Bitmap SobelFilter(Bitmap bmp)
{
    byte[] buffer = BitmapToBuffer(bmp);

    int width = bmp.Width;
    int height = bmp.Height;
    int pixelsCount = width * height;
    byte[] pixelBuffer = new byte[pixelsCount];
    byte[] resultBuffer = new byte[pixelsCount];

    for (int i = 0; i < pixelsCount; i++)
    {
        int offset = i * 4;
        byte brightness = (byte)(buffer[offset] * 0.29891f + buffer[offset + 1] * 0.58661f + buffer[offset + 2] * 0.118f);
        pixelBuffer[i] = brightness;
    }

    for (int i = width + 1; i < pixelsCount - width - 1; i++)
    {
        if (i % width == width - 1)
            i += 2;

        int gx = -pixelBuffer[i - 1 - width] + pixelBuffer[i + 1 - width] - 2 * pixelBuffer[i - 1] +
            2 * pixelBuffer[i + 1] - pixelBuffer[i - 1 + width] + pixelBuffer[i + 1 + width];

        int gy = pixelBuffer[i - 1 - width] + 2 * pixelBuffer[i - width] + pixelBuffer[i + 1 - width] -
            pixelBuffer[i - 1 + width] - 2 * pixelBuffer[i + width] - pixelBuffer[i + 1 + width];

        int gt = (int)MathF.Sqrt(gx * gx + gy * gy);
        if (gt > byte.MaxValue) gt = byte.MaxValue;

        resultBuffer[i] = (byte)gt;
    }

    return BufferToBitmap(resultBuffer, width, height);
}

Первым делом я реализовал применение ядра собеля к изображению, чтобы получить заветный вектор с пикселями, в процессе оптимизации очень хорошо подошел здесь AVX2 на 256 векторах, который отработал быстрее, чем SSE2 на 128 векторах раза в полтора.

private static unsafe Vector256<int> ApplySobelKernelVectorized(byte* srcPtr, int width)
{
    Vector256<int> v00 = Avx2.ConvertToVector256Int32(srcPtr);
    Vector256<int> v01 = Avx2.ConvertToVector256Int32(srcPtr + 1);
    Vector256<int> v02 = Avx2.ConvertToVector256Int32(srcPtr + 2);
    Vector256<int> v10 = Avx2.ConvertToVector256Int32(srcPtr + width);
    Vector256<int> v12 = Avx2.ConvertToVector256Int32(srcPtr + width + 2);
    Vector256<int> v20 = Avx2.ConvertToVector256Int32(srcPtr + width * 2);
    Vector256<int> v21 = Avx2.ConvertToVector256Int32(srcPtr + width * 2 + 1);
    Vector256<int> v22 = Avx2.ConvertToVector256Int32(srcPtr + width * 2 + 2);

    Vector256<int> vgx = Avx2.Subtract(v02, v00);
    vgx = Avx2.Subtract(vgx, Avx2.ShiftLeftLogical(v10, 1));
    vgx = Avx2.Add(vgx, Avx2.ShiftLeftLogical(v12, 1));
    vgx = Avx2.Subtract(vgx, v20);
    vgx = Avx2.Add(vgx, v22);
    vgx = Avx2.MultiplyLow(vgx, vgx);

    Vector256<int> vgy = Avx2.Add(v00, Avx2.ShiftLeftLogical(v01, 1));
    vgy = Avx2.Add(vgy, v02);
    vgy = Avx2.Subtract(vgy, v20);
    vgy = Avx2.Subtract(vgy, Avx2.ShiftLeftLogical(v21, 1));
    vgy = Avx2.Subtract(vgy, v22);
    vgy = Avx2.MultiplyLow(vgy, vgy);

    Vector256<float> gtf = Avx.Sqrt(Avx.ConvertToVector256Single(Avx2.Add(vgx, vgy)));
    return Avx.ConvertToVector256Int32WithTruncation(gtf);
}

Ну и затем уже не составило труда по накатанной реализовать основной метод

private static Bitmap SobelFilterVectorized(Bitmap bmp)
{
    byte[] buffer = BitmapToBuffer(bmp);

    int width = bmp.Width;
    int height = bmp.Height;
    int pixelsCount = width * height;
    byte[] pixelBuffer = new byte[pixelsCount];
    byte[] resultBuffer = new byte[pixelsCount];
    unsafe
    {
        fixed (byte* srcPtr = buffer, pixPtr = pixelBuffer, dstPtr = resultBuffer)
        {
            ToGrayscaleVectorized(srcPtr, pixPtr, pixelsCount);

            byte* srcPos = pixPtr;
            byte* srcEnd = pixPtr + pixelsCount - width * 2 - 1;
            byte* dstPos = dstPtr + width + 1;

            while (srcPos < srcEnd)
            {
                Vector256<int> gt0 = ApplySobelKernelVectorized(srcPos, width);
                Vector256<int> gt1 = ApplySobelKernelVectorized(srcPos + 8, width);

                Vector128<short> gts0 = Sse2.PackSignedSaturate(Avx2.ExtractVector128(gt0, 0), Avx2.ExtractVector128(gt0, 1));
                Vector128<short> gts1 = Sse2.PackSignedSaturate(Avx2.ExtractVector128(gt1, 0), Avx2.ExtractVector128(gt1, 1));
                Sse2.Store(dstPos, Sse2.PackUnsignedSaturate(gts0, gts1));
                srcPos += 16;
                dstPos += 16;
            }

            // так как исключить границы картинки при векторизации практически невозможно
            // без сильного ущерба производительности, поэтому
            // я закрашиваю левую и правую границу черными пикселями
            for (dstPos = dstPtr + width; dstPos <= dstPtr + pixelsCount - width; dstPos += width)
            {
                *dstPos-- = 0;
                *dstPos++ = 0;
            }
        }
    }
    return BufferToBitmap(resultBuffer, width, height);
}

Вывод в консоль для изображения 2560x1358, релизная сборка с оптимизацией кода:

SSE: True
SSE2: True
SSSE3: True
SSE41: True
AVX: True
AVX2: True
Optimized: 71,7946ms
Vectorized: 34,5113ms
Original: 2314,0325ms

Тестовая картинка.

^{(картинка не с инета, своя, авторскими правами использование не ограничено)}

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

Разница между моей оптимизированной и векторной версией есть, но почти нет :), разницу видно только если быстро перещелкивать картинки друг над другом. Возможно связано с нюансами округления в SIMD инструкциях. Обновлено: после доработки метода обесцвечивания разница пропала, теперь картинки, произведенные скалярной и векторной версиями кода идентичны.

P.S. csproj, на всякий случай

<Project Sdk="Microsoft.NET.Sdk">
  <PropertyGroup>
    <OutputType>Exe</OutputType>
    <TargetFramework>net5.0-windows</TargetFramework>
    <Platforms>x64</Platforms>
  </PropertyGroup>
  <ItemGroup>
    <PackageReference Include="System.Drawing.Common" Version="5.0.0" />
  </ItemGroup>
</Project>

Обновление - Code Review

Финальный код после всех фиксов и оптимизаций можно взять здесь - Code Review

Answer 2

Крайне любопытно.

Я переписал код @aepot на указатели. Это убирает два копирования массива с помощью Marshal.Copy, что, по идее, ещё чуть-чуть ускоряет код.

Код не доделан и не на всех типах изображений будет работать правильно.

const int bpp = 3; // байт на пиксель

var bmp = new Bitmap("source.jpg");
Rectangle rect = new Rectangle(0, 0, bmp.Width, bmp.Height);
BitmapData bmpData = bmp.LockBits(rect, ImageLockMode.ReadWrite, bmp.PixelFormat);
IntPtr ptr = bmpData.Scan0;
int stride = Math.Abs(bmpData.Stride);
int size = stride * bmp.Height;

byte[] bytes = new byte[size];
GCHandle handle = GCHandle.Alloc(bytes, GCHandleType.Pinned);
IntPtr addr = handle.AddrOfPinnedObject();
var result = new Bitmap(bmp.Width, bmp.Height, stride, bmp.PixelFormat, addr);

unsafe
{
    byte* src = (byte*)ptr.ToPointer();
    byte* dst = (byte*)addr.ToPointer();

    for (int i = 0; i < size; i += bpp)
    {
        byte brightness = (byte)(src[i] * 0.3f + src[i + 1] * 0.59f + src[i + 2] * 0.11f);
        src[i] = brightness;
    }

    for (int i = stride + bpp; i < size - stride - bpp; i += bpp)
    {
        if (i % stride == stride - bpp)
            i += bpp * 2;

        int gx = -src[i - bpp - stride] + src[i + bpp - stride] - 2 * src[i - bpp] +
            2 * src[i + bpp] - src[i - bpp + stride] + src[i + bpp + stride];

        int gy = src[i - bpp - stride] + 2 * src[i - stride] + src[i + bpp - stride] -
            src[i - bpp + stride] - 2 * src[i + stride] - src[i + bpp + stride];

        int gt = (int)MathF.Sqrt(gx * gx + gy * gy);
        if (gt > byte.MaxValue) gt = byte.MaxValue;

        dst[i] = (byte)gt;
        dst[i + 1] = (byte)gt;
        dst[i + 2] = (byte)gt;
    }
}

bmp.UnlockBits(bmpData);
result.Save("result.bmp");
handle.Free();

Применить SIMD-операции весьма проблематично. Ведь нужен доступ не к последовательным элементам, а расположенным в разных колонках и строках. Но я ещё подумаю.

БЛОГ НА HUSL

Оператор Собеля - SIMD реализация. SSE/AVX

Ответы (2 шт):

Ох, это было не просто.

Обновление - Code Review