Sobel Operator - SIMD x86 Intrinsics Implementation

Question

I'm learning C# .NET 5 Intrinsics and interested in best practices. It's really hard now to find enough information about how SIMD instructions (logically/internally) work in .NET. In addition i'm not familiar with C++ and Assembly languages.

The purpose of the solution - apply Sobel Operator filter to the loaded image. For image operations i used System.Drawing.Common NuGet package. Thus the solution is Windows-only.

The SobelOperator class contains two Sobel Operator implementations:

1. SobelOperatorScalar - Scalar solution that can be used as a fallback if current CPU is not compartible with AVX2.
2. SobelOperatorSimd - SIMD x86 solution for hardware acceleration. - review target

Code for review

public interface ISobelOperator
{
    Bitmap Apply(Bitmap bmp);
}

public class SobelOperator : ISobelOperator
{
    private static Color[] _grayPallette;
    private readonly ISobelOperator _operator;

    public bool IsHardwareAccelerated { get; } 
            
    public SobelOperator(bool hardwareAccelerated = true)
    {
        if (_grayPallette == null)
            _grayPallette = Enumerable.Range(0, 256).Select(i => Color.FromArgb(i, i, i)).ToArray();

        IsHardwareAccelerated = hardwareAccelerated && Avx2.IsSupported;

        _operator = IsHardwareAccelerated ? new SobelOperatorSimd() : new SobelOperatorScalar();
    }

    public Bitmap Apply(Bitmap bmp) 
        => _operator.Apply(bmp);

    private class SobelOperatorSimd : ISobelOperator
    {
        private const byte m0 = 0b01001001;
        private const byte m1 = 0b10010010;
        private const byte m2 = 0b00100100;

        //0.299R + 0.587G + 0.114B
        private readonly Vector256<float> bWeight = Vector256.Create(0.114f);
        private readonly Vector256<float> gWeight = Vector256.Create(0.587f);
        private readonly Vector256<float> rWeight = Vector256.Create(0.299f);
        private readonly Vector256<int> bMut = Vector256.Create(0, 3, 6, 1, 4, 7, 2, 5);
        private readonly Vector256<int> gMut = Vector256.Create(1, 4, 7, 2, 5, 0, 3, 6);
        private readonly Vector256<int> rMut = Vector256.Create(2, 5, 0, 3, 6, 1, 4, 7);

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private unsafe Vector256<int> GetBrightness(byte* ptr)
        {
            Vector256<int> v0 = Avx2.ConvertToVector256Int32(ptr);
            Vector256<int> v1 = Avx2.ConvertToVector256Int32(ptr + 8);
            Vector256<int> v2 = Avx2.ConvertToVector256Int32(ptr + 16);

            Vector256<int> vb = Avx2.Blend(Avx2.Blend(v0, v1, m1), v2, m2);
            vb = Avx2.PermuteVar8x32(vb, bMut);
            Vector256<int> vg = Avx2.Blend(Avx2.Blend(v0, v1, m2), v2, m0);
            vg = Avx2.PermuteVar8x32(vg, gMut);
            Vector256<int> vr = Avx2.Blend(Avx2.Blend(v0, v1, m0), v2, m1);
            vr = Avx2.PermuteVar8x32(vr, rMut);
            Vector256<float> vfb = Avx.Multiply(Avx.ConvertToVector256Single(vb), bWeight);
            Vector256<float> vfg = Avx.Multiply(Avx.ConvertToVector256Single(vg), gWeight);
            Vector256<float> vfr = Avx.Multiply(Avx.ConvertToVector256Single(vr), rWeight);
            return Avx.ConvertToVector256Int32WithTruncation(Avx.Add(Avx.Add(vfb, vfg), vfr));
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private unsafe void ToGrayscale(byte* srcPtr, byte* dstPtr, int pixelsCount)
        {
            byte* tail = srcPtr + (pixelsCount & -16) * 3;
            byte* srcEnd = srcPtr + pixelsCount * 3;
            byte* dstEnd = dstPtr + pixelsCount;

            while (true)
            {
                while (srcPtr < tail)
                {
                    Vector256<int> vi0 = GetBrightness(srcPtr);
                    Vector256<int> vi1 = GetBrightness(srcPtr + 24);
                    Vector128<short> v0 = Sse2.PackSignedSaturate(Avx2.ExtractVector128(vi0, 0), Avx2.ExtractVector128(vi0, 1));
                    Vector128<short> v1 = Sse2.PackSignedSaturate(Avx2.ExtractVector128(vi1, 0), Avx2.ExtractVector128(vi1, 1));
                    Sse2.Store(dstPtr, Sse2.PackUnsignedSaturate(v0, v1));
                    srcPtr += 48;
                    dstPtr += 16;
                }
                if (srcPtr == srcEnd)
                    break;
                tail = srcEnd;
                srcPtr = srcEnd - 48;
                dstPtr = dstEnd - 16;
            }
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private static unsafe Vector128<byte> ApplySobelKernel(byte* srcPtr, int width)
        {
            Vector256<short> v00 = Avx2.ConvertToVector256Int16(srcPtr);
            Vector256<short> v01 = Avx2.ConvertToVector256Int16(srcPtr + 1);
            Vector256<short> v02 = Avx2.ConvertToVector256Int16(srcPtr + 2);
            Vector256<short> v10 = Avx2.ConvertToVector256Int16(srcPtr + width);
            Vector256<short> v12 = Avx2.ConvertToVector256Int16(srcPtr + width + 2);
            Vector256<short> v20 = Avx2.ConvertToVector256Int16(srcPtr + width * 2);
            Vector256<short> v21 = Avx2.ConvertToVector256Int16(srcPtr + width * 2 + 1);
            Vector256<short> v22 = Avx2.ConvertToVector256Int16(srcPtr + width * 2 + 2);

            Vector256<short> 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);

            Vector256<short> 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);

            // sqrt(vgx * vgx + vgy * vgy)
            Vector256<short> vgp0 = Avx2.UnpackLow(vgx, vgy);
            Vector256<short> vgp1 = Avx2.UnpackHigh(vgx, vgy);
            Vector256<int> v0 = Avx2.MultiplyAddAdjacent(vgp0, vgp0);
            Vector256<int> v1 = Avx2.MultiplyAddAdjacent(vgp1, vgp1);
            Vector256<int> gt0 = Avx.ConvertToVector256Int32WithTruncation(Avx.Sqrt(Avx.ConvertToVector256Single(v0)));
            Vector256<int> gt1 = Avx.ConvertToVector256Int32WithTruncation(Avx.Sqrt(Avx.ConvertToVector256Single(v1)));

            Vector128<short> gts0 = Sse2.PackSignedSaturate(Avx2.ExtractVector128(gt0, 0), Avx2.ExtractVector128(gt1, 0));
            Vector128<short> gts1 = Sse2.PackSignedSaturate(Avx2.ExtractVector128(gt0, 1), Avx2.ExtractVector128(gt1, 1));
            return Sse2.PackUnsignedSaturate(gts0, gts1);
        }


        public Bitmap Apply(Bitmap bmp)
        {
            int width = bmp.Width;
            int height = bmp.Height;
            int pixelsCount = width * height;
            byte[] buffer = new byte[pixelsCount];

            Rectangle rect = new Rectangle(Point.Empty, bmp.Size);
            Bitmap outBmp = new Bitmap(width, height, PixelFormat.Format8bppIndexed);
            ColorPalette pal = outBmp.Palette;
            for (int i = 0; i < 256; i++)
                pal.Entries[i] = _grayPallette[i];
            outBmp.Palette = pal;

            unsafe
            {
                fixed (byte* bufPtr = buffer)
                {
                    BitmapData bmpData = bmp.LockBits(rect, ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);
                    ToGrayscale((byte*)bmpData.Scan0.ToPointer(), bufPtr, pixelsCount);
                    bmp.UnlockBits(bmpData);

                    BitmapData outBmpData = outBmp.LockBits(rect, ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);
                    byte* dstPtr = (byte*)outBmpData.Scan0.ToPointer();

                    int length = pixelsCount - width * 2 - 1;
                    byte* tail = bufPtr + (length & -16);
                    byte* srcPos = bufPtr;
                    byte* srcEnd = bufPtr + length;
                    byte* dstPos = dstPtr + width + 1;
                    byte* dstEnd = dstPos + length;

                    while (true)
                    {
                        while (srcPos < tail)
                        {
                            Sse2.Store(dstPos, ApplySobelKernel(srcPos, width));
                            srcPos += 16;
                            dstPos += 16;
                        }

                        if (srcPos == srcEnd)
                            break;
                        tail = srcEnd;
                        srcPos = srcEnd - 16;
                        dstPos = dstEnd - 16;
                    }

                    for (dstPos = dstPtr + width; dstPos <= dstPtr + pixelsCount - width; dstPos += width)
                    {
                        *dstPos-- = 0;
                        *dstPos++ = 0;
                    }

                    outBmp.UnlockBits(outBmpData);
                }
            }

            return outBmp;
        }
    }

    private class SobelOperatorScalar : ISobelOperator
    {
        public Bitmap Apply(Bitmap bmp)
        {
            BitmapData bmpData = bmp.LockBits(new Rectangle(Point.Empty, bmp.Size), ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);
            int strideLength = bmpData.Stride * bmpData.Height;
            byte[] buffer = new byte[Math.Abs(strideLength)];
            Marshal.Copy(bmpData.Scan0, buffer, 0, strideLength);
            bmp.UnlockBits(bmpData);

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

            //0.299R + 0.587G + 0.114B
            for (int i = 0; i < pixelsCount; i++)
            {
                int offset = i * 3;
                byte brightness = (byte)(buffer[offset] * 0.114f + buffer[offset + 1] * 0.587f + buffer[offset + 2] * 0.299f);
                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;
            }

            Bitmap outBmp = new Bitmap(width, height, PixelFormat.Format8bppIndexed);
            BitmapData outBmpData = outBmp.LockBits(new Rectangle(Point.Empty, outBmp.Size), ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);
            Marshal.Copy(resultBuffer, 0, outBmpData.Scan0, outBmpData.Stride * outBmpData.Height);
            outBmp.UnlockBits(outBmpData);
            ColorPalette pal = outBmp.Palette;
            for (int i = 0; i < 256; i++)
                pal.Entries[i] = _grayPallette[i];
            outBmp.Palette = pal;
            return outBmp;
        }
    }
}

Output test

Test image

Program.cs

static void Main(string[] args)
{
    const string fileName = "image.jpg";
    Bitmap bmp = new Bitmap(fileName);

    SobelOperator sobelOperator = new SobelOperator();
    Console.WriteLine($"SIMD accelerated: {(sobelOperator.IsHardwareAccelerated ? "Yes" : "No")}");
    Bitmap result = sobelOperator.Apply(bmp);
    result.Save("out.jpg", ImageFormat.Jpeg);

    Console.WriteLine("Done.");
    Console.ReadKey();
}

Console output

SIMD accelerated: Yes
Done.

Output Image

Output images of Scalar and SIMD implementations are binary identical.

Benchmark.NET

[MemoryDiagnoser]
public class MyBenchmark
{
    private readonly ISobelOperator _sobelOperator = new SobelOperator();
    private readonly ISobelOperator _sobelOperatorSw = new SobelOperator(false);
    private readonly Bitmap bmp = new Bitmap(@"C:\Source\image.jpg");

    [Benchmark(Description = "SIMD Enabled")]
    public Bitmap TestSimd()
    {
        return _sobelOperator.Apply(bmp);
    }

    [Benchmark(Description = "SIMD Disabled")]
    public Bitmap TestScalar()
    {
        return _sobelOperatorSw.Apply(bmp);
    }
}

static void Main(string[] args)
{
    var summary = BenchmarkRunner.Run<MyBenchmark>();
    Console.ReadKey();
}

BenchmarkDotNet=v0.12.1, OS=Windows 10.0.19042
Intel Core i7-4700HQ CPU 2.40GHz (Haswell), 1 CPU, 8 logical and 4 physical cores
.NET Core SDK=5.0.101
  [Host]     : .NET Core 5.0.1 (CoreCLR 5.0.120.57516, CoreFX 5.0.120.57516), X64 RyuJIT
  DefaultJob : .NET Core 5.0.1 (CoreCLR 5.0.120.57516, CoreFX 5.0.120.57516), X64 RyuJIT

Method	Mean	Error	StdDev	Gen 0	Gen 1	Gen 2	Allocated
'SIMD Enabled'	7.285 ms	0.1165 ms	0.1089 ms	992.1875	992.1875	992.1875	3.35 MB
'SIMD Disabled'	48.412 ms	0.2312 ms	0.2162 ms	454.5455	454.5455	454.5455	16.61 MB

Intrinsics solution is ~6.6x Times faster. And eating less memory in general because it's unsafe and doesn't use Marshal.Copy to load/save the byte[] buffers.

Answer 1

Bugs for images of inconvenient width

This code:

Bitmap outBmp = new Bitmap(width, height, PixelFormat.Format8bppIndexed);
BitmapData outBmpData = outBmp.LockBits(new Rectangle(Point.Empty, outBmp.Size), ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);
Marshal.Copy(resultBuffer, 0, outBmpData.Scan0, outBmpData.Stride * outBmpData.Height);

Works great if the width of the image is a multiple of 4. But in other cases, it does not work. The problem is that bitmaps have a stride is a multiple of 4, if their "natural byte-width" is not a multiple of 4 then there will be padding at the end of every row. This code does not account for that, so for images that have padding it tries to copy more data from the result buffer than is actually in it, which causes an exception. Even if it had worked it still would have been wrong, resulting in every row of pixels "losing" some pixels to the padding, shearing the image.

This issue also affects other places where a bitmap is locked with a pixel type that is not 4 bytes, but it affects them in different ways.

Generally it means that we need to work with RGB and 8bit image data row-by-row, annoying as it is. ARGB is nicer in this regard, but locking an RGB image into ARGB format has a large cost - from a performance perspective it's much worse than being forced to work with RGB data (but the code we have to write to deal with RGB is really nasty..) Your code has some really nice loops, it's a shame that they have to be replaced by uglier ones, at least if you want images of various widths to work.

An other small detail is that very small widths don't work with the "step back from the end"-trick to do the last iteration. In those cases you may as well fall back to the scalar implementation - there is no hope for SIMD and they're small images anyway.

Performance improvements that may change the result

For the benchmarks, I note that you're using an Intel Haswell, so I benchmarked on an Intel Haswell as well. That may matter: if code A is better than code B on Haswell, the situation may be reversed on Skylake or AMD Ryzen.

The square root in the Sobel kernel can be replaced by Avx.Multiply(Avx.ReciprocalSqrt(v0f), v0f) where v0f is v0 converted to floats. The reciprocal square root has about 11 to 12 bits of precision, and the result is converted to 8 bits per pixel anyway. On the test image, it didn't make any difference for the results, but made the kernel slightly faster. In general, I cannot guarantee that there is never a difference in result, especially when taking the square root of a perfect square (which sits right on the knife edge of being truncated the wrong way if the reciprocal square root is only just too low).

The grayscale conversion can also be done slightly faster. The main idea there, based somewhat on this code, is to use Avx2.MultiplyHigh on a vector of ushort instead of floating point multiplication. Everything else happens in support of that. The way it works out though, is that vectors end up with "wasted bits" where the result of the computation is thrown away. Due to that, it still works out to 6 multiplication instructions per 16 bytes of output even though twice as many multiplications happen per instruction. So in terms of the number of multiplication instructions, it isn't any better than it used to be, and on top of that, on Haswell integer multiplication has half throughput compared to floating point multiplication.

The rest of the function does not seem especially nice either, being extremely shuffle-heavy (Haswell can only do one shuffle per cycle, so the number of shuffles is important). So actually, I don't know why it's faster, I can only tell you that the benchmarks said so and that the difference was significant despite the relatively highly timing variance that I kept seeing. Hopefully a better approach is possible but I don't know what it is.

const ushort rw = 19595;
const ushort gw = 38470;
const ushort bw = 7471;
Vector256<ushort> rW = Vector256.Create(rw);
Vector256<ushort> gW = Vector256.Create(gw);
Vector256<ushort> bW = Vector256.Create(bw);
const byte _ = 0x80;
Vector256<byte> shuf0 = Vector256.Create(
    _, 0, _, 3, _, 6, _, 9, _, 12, _, _, _, _, _, _,
    _, 0, _, 3, _, 6, _, 9, _, 12, _, _, _, _, _, _);
Vector256<byte> shuf1 = Vector256.Create(
    _, 1, _, 4, _, 7, _, 10, _, 13, _, _, _, _, _, _,
    _, 1, _, 4, _, 7, _, 10, _, 13, _, _, _, _, _, _);
Vector256<byte> shuf2 = Vector256.Create(
    _, 2, _, 5, _, 8, _, 11, _, 14, _, _, _, _, _, _,
    _, 2, _, 5, _, 8, _, 11, _, 14, _, _, _, _, _, _);
Vector256<byte> shuf3 = Vector256.Create(
    1, 3, 5, 7, 9, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, 1, 3, 5, 7, 9, _, _, _, _, _, _);
Vector256<byte> shuf4 = Vector256.Create(
    _, _, _, _, _, _, _, _, _, _, 9, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, 1, 3, 5, 7, 9);
Vector256<ushort> bias = Vector256.Create((ushort)0x02);

[MethodImpl(MethodImplOptions.AggressiveInlining)]
private unsafe Vector128<byte> GetBrightness16(byte* ptr)
{
    var raw0 = Avx2.LoadVector128(ptr);
    var raw1 = Avx2.LoadVector128(ptr + 16);
    var raw2 = Avx2.LoadVector128(ptr + 32);
    var partA = Avx2.InsertVector128(raw0.ToVector256(), Avx2.AlignRight(raw1, raw0, 15), 1);
    var partB = Avx2.InsertVector128(Avx2.AlignRight(raw2, raw1, 1).ToVector256(), raw2, 1);
    partB = Avx2.AlignRight(partB, partB, 1);

    var b0 = Avx2.Shuffle(partA, shuf0).AsUInt16();
    var g0 = Avx2.Shuffle(partA, shuf1).AsUInt16();
    var r0 = Avx2.Shuffle(partA, shuf2).AsUInt16();
    var b1 = Avx2.Shuffle(partB, shuf0).AsUInt16();
    var g1 = Avx2.Shuffle(partB, shuf1).AsUInt16();
    var r1 = Avx2.Shuffle(partB, shuf2).AsUInt16();

    b0 = Avx2.MultiplyHigh(b0, bW);
    g0 = Avx2.MultiplyHigh(g0, gW);
    r0 = Avx2.MultiplyHigh(r0, rW);
    b1 = Avx2.MultiplyHigh(b1, bW);
    g1 = Avx2.MultiplyHigh(g1, gW);
    r1 = Avx2.MultiplyHigh(r1, rW);
    var sum0 = Avx2.AddSaturate(Avx2.AddSaturate(Avx2.Add(b0, g0), r0), bias);
    var shufsum0 = Avx2.Shuffle(sum0.AsByte(), shuf3);
    var sum1 = Avx2.AddSaturate(Avx2.AddSaturate(Avx2.Add(b1, g1), r1), bias);
    var shufsum1 = Avx2.Shuffle(sum1.AsByte(), shuf4);
    var shufsum = Avx2.Or(shufsum0, shufsum1);
    return Avx2.Or(Vector256.GetLower(shufsum), Vector256.GetUpper(shufsum));
}

By the way, I chose the bias to minimize the difference with the original grayscale conversion. Normally I would use a bias of 0x80 to round evenly. The weights of the color channels are approximations of 65536 times their floating point weight, chosen so that they add up to 65536. MultiplyHigh implicitly divides by 65536, so effectively these scales work like their floating point counterparts. They don't have to add up to exactly 65536 by the way: thanks to the saturating additions (which are no extra cost compared to regular additions) nothing bad would happen if they added up to slightly more. You could use that property to tune them to more closely match the original conversion.

Alternatively, it would be easy to adapt the scalar function to fixed-point arithmetic such that it gives the same results as this version of the SIMD grayscale converter.

I tried using unaligned loads to replace the Avx2.AlignRight (aka vpalignr) instructions with, but the benchmark said that was marginally slower.

It's possible to base a grayscale conversion on vpmaddubsw (Avx2.MultiplyAddAdjacent with vectors of byte and sbyte). For ARGB data that is very fast, though noticably less accurate. I have not worked that out for RGB data.

Small things

ExtractVector128-ing the low half of a vector

If you write Avx2.ExtractVector128(gt0, 0), you get what you asked for, VEXTRACTI128 with an index of zero. However, that's not the best instruction to use, the best instruction to use is .. nothing. No instruction is necessary. What should happen is that the next operation only uses the 128bit version of the same vector register. The way to tell C# to do that is with Vector256.GetLower.

This is a contrast with eg C++, where compilers normally do that substitution on their own.

The performance difference from this, if any, was negligible. Saving an explicit extraction instruction can only help though.

Unnecessary intptr.ToPointer()

At various times, the code contains (byte*)intptr.ToPointer() or some variant. It's enough to cast to the pointer, calling ToPointer() doesn't add anything to that code.

Answer 2

Reviewed Code

public interface ISobelOperator
{
    Bitmap Apply(Bitmap bmp);
}

public class SobelOperator : ISobelOperator
{
    private static Color[] _grayPallette;
    private readonly ISobelOperator _operator;

    public bool IsHardwareAccelerated { get; }

    public SobelOperator(bool hardwareAccelerated = true)
    {
        if (_grayPallette == null)
            _grayPallette = Enumerable.Range(0, 256).Select(i => Color.FromArgb(i, i, i)).ToArray();

        IsHardwareAccelerated = hardwareAccelerated && Avx2.IsSupported;

        _operator = IsHardwareAccelerated ? new SobelOperatorSimd() : new SobelOperatorScalar();
    }

    public Bitmap Apply(Bitmap bmp)
        => _operator.Apply(bmp);

    private class SobelOperatorSimd : ISobelOperator
    {
        private readonly Vector256<ushort> rw = Vector256.Create((ushort)19595);
        private readonly Vector256<ushort> gw = Vector256.Create((ushort)38470);
        private readonly Vector256<ushort> bw = Vector256.Create((ushort)7471);
        private const byte _ = 0x80;
        private readonly Vector256<byte> mask0 = Vector256.Create(_, 0, _, 3, _, 6, _, 9, _, 12, _, _, _, _, _, _, _, 0, _, 3, _, 6, _, 9, _, 12, _, _, _, _, _, _);
        private readonly Vector256<byte> mask1 = Vector256.Create(_, 1, _, 4, _, 7, _, 10, _, 13, _, _, _, _, _, _, _, 1, _, 4, _, 7, _, 10, _, 13, _, _, _, _, _, _);
        private readonly Vector256<byte> mask2 = Vector256.Create(_, 2, _, 5, _, 8, _, 11, _, 14, _, _, _, _, _, _, _, 2, _, 5, _, 8, _, 11, _, 14, _, _, _, _, _, _);
        private readonly Vector256<byte> mask3 = Vector256.Create(1, 3, 5, 7, 9, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, 1, 3, 5, 7, 9, _, _, _, _, _, _);
        private readonly Vector256<byte> mask4 = Vector256.Create(_, _, _, _, _, _, _, _, _, _, 9, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, 1, 3, 5, 7, 9);

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private unsafe void ToGrayscale(byte* srcPtr, byte* dstPtr, int srcStride, int dstStride, int height)
        {
            for (int row = 0; row < height; row++)
            {
                byte* srcPos = srcPtr + row * srcStride;
                byte* dstPos = dstPtr + row * dstStride;
                byte* srcEnd = srcPos + srcStride;

                while (srcPos < srcEnd)
                {
                    Vector128<byte> raw0 = Sse2.LoadVector128(srcPos);
                    Vector128<byte> raw1 = Sse2.LoadVector128(srcPos + 16);
                    Vector128<byte> raw2 = Sse2.LoadVector128(srcPos + 32);
                    Vector256<byte> v0 = Avx2.InsertVector128(raw0.ToVector256(), Ssse3.AlignRight(raw1, raw0, 15), 1);
                    Vector256<byte> v1 = Avx2.InsertVector128(Ssse3.AlignRight(raw2, raw1, 1).ToVector256(), raw2, 1);
                    v1 = Avx2.AlignRight(v1, v1, 1);

                    Vector256<ushort> b0 = Avx2.MultiplyHigh(Avx2.Shuffle(v0, mask0).AsUInt16(), bw);
                    Vector256<ushort> g0 = Avx2.MultiplyHigh(Avx2.Shuffle(v0, mask1).AsUInt16(), gw);
                    Vector256<ushort> r0 = Avx2.MultiplyHigh(Avx2.Shuffle(v0, mask2).AsUInt16(), rw);
                    Vector256<ushort> b1 = Avx2.MultiplyHigh(Avx2.Shuffle(v1, mask0).AsUInt16(), bw);
                    Vector256<ushort> g1 = Avx2.MultiplyHigh(Avx2.Shuffle(v1, mask1).AsUInt16(), gw);
                    Vector256<ushort> r1 = Avx2.MultiplyHigh(Avx2.Shuffle(v1, mask2).AsUInt16(), rw);

                    Vector256<byte> sum0 = Avx2.Shuffle(Avx2.AddSaturate(Avx2.Add(b0, g0), r0).AsByte(), mask3);
                    Vector256<byte> sum1 = Avx2.Shuffle(Avx2.AddSaturate(Avx2.Add(b1, g1), r1).AsByte(), mask4);
                    Vector256<byte> sum = Avx2.Or(sum0, sum1);
                    Sse2.Store(dstPos, Sse2.Or(sum.GetLower(), sum.GetUpper()));
                    srcPos += 48;
                    dstPos += 16;
                }
            }
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private unsafe Vector128<byte> ApplySobelKernel(byte* srcPtr, int stride)
        {
            Vector256<short> v00 = Avx2.ConvertToVector256Int16(srcPtr);
            Vector256<short> v01 = Avx2.ConvertToVector256Int16(srcPtr + 1);
            Vector256<short> v02 = Avx2.ConvertToVector256Int16(srcPtr + 2);
            Vector256<short> v10 = Avx2.ConvertToVector256Int16(srcPtr + stride);
            Vector256<short> v12 = Avx2.ConvertToVector256Int16(srcPtr + stride + 2);
            Vector256<short> v20 = Avx2.ConvertToVector256Int16(srcPtr + stride * 2);
            Vector256<short> v21 = Avx2.ConvertToVector256Int16(srcPtr + stride * 2 + 1);
            Vector256<short> v22 = Avx2.ConvertToVector256Int16(srcPtr + stride * 2 + 2);

            Vector256<short> 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);

            Vector256<short> 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);

            // sqrt(vgx * vgx + vgy * vgy)
            Vector256<short> vgp0 = Avx2.UnpackLow(vgx, vgy);
            Vector256<short> vgp1 = Avx2.UnpackHigh(vgx, vgy);
            Vector256<int> v0 = Avx2.MultiplyAddAdjacent(vgp0, vgp0);
            Vector256<int> v1 = Avx2.MultiplyAddAdjacent(vgp1, vgp1);
            Vector256<float> vf0 = Avx.ConvertToVector256Single(v0);
            Vector256<float> vf1 = Avx.ConvertToVector256Single(v1);
            Vector256<int> gt0 = Avx.ConvertToVector256Int32WithTruncation(Avx.Multiply(Avx.ReciprocalSqrt(vf0), vf0));
            Vector256<int> gt1 = Avx.ConvertToVector256Int32WithTruncation(Avx.Multiply(Avx.ReciprocalSqrt(vf1), vf1));

            Vector128<short> gts0 = Sse2.PackSignedSaturate(gt0.GetLower(), gt1.GetLower());
            Vector128<short> gts1 = Sse2.PackSignedSaturate(gt0.GetUpper(), gt1.GetUpper());
            return Sse2.PackUnsignedSaturate(gts0, gts1);
        }


        public unsafe Bitmap Apply(Bitmap bmp)
        {
            int width = bmp.Width;
            int height = bmp.Height;

            Rectangle rect = new Rectangle(Point.Empty, bmp.Size);
            BitmapData bmpData = bmp.LockBits(rect, ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);
            int srcStride = bmpData.Stride;
            int srcLength = srcStride * height;

            Bitmap outBmp = new Bitmap(width, height, PixelFormat.Format8bppIndexed);
            ColorPalette pal = outBmp.Palette;
            for (int i = 0; i < 256; i++)
                pal.Entries[i] = _grayPallette[i];
            outBmp.Palette = pal;

            BitmapData outBmpData = outBmp.LockBits(rect, ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);
            int dstStride = outBmpData.Stride;
            int dstLength = dstStride * height;

            byte[] buffer = new byte[dstLength];

            fixed (byte* bufPtr = buffer)
            {
                ToGrayscale((byte*)bmpData.Scan0, bufPtr, srcStride, dstStride, height);
                bmp.UnlockBits(bmpData);

                byte* dstPtr = (byte*)outBmpData.Scan0;

                int length = dstLength - dstStride * 2 - 1;
                byte* tail = bufPtr + (length & -16);
                byte* srcPos = bufPtr;
                byte* srcEnd = bufPtr + length;
                byte* dstPos = dstPtr + dstStride + 1;
                byte* dstEnd = dstPos + length;

                while (true)
                {
                    while (srcPos < tail)
                    {
                        Sse2.Store(dstPos, ApplySobelKernel(srcPos, dstStride));
                        srcPos += 16;
                        dstPos += 16;
                    }

                    if (srcPos == srcEnd)
                        break;
                    tail = srcEnd;
                    srcPos = srcEnd - 16;
                    dstPos = dstEnd - 16;
                }

                int padding = dstStride - width + 1;
                for (dstPos = dstPtr + dstStride; dstPos <= dstPtr + dstLength - dstStride; dstPos += dstStride)
                {
                    *dstPos = 0;
                    *(dstPos - padding) = 0;
                }
            }
            outBmp.UnlockBits(outBmpData);
            return outBmp;
        }
    }

    private class SobelOperatorScalar : ISobelOperator
    {
        public Bitmap Apply(Bitmap bmp)
        {
            int width = bmp.Width;
            int height = bmp.Height;

            Rectangle rect = new Rectangle(Point.Empty, bmp.Size);
            BitmapData bmpData = bmp.LockBits(rect, ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);
            int inStride = Math.Abs(bmpData.Stride);
            int strideLength = bmpData.Stride * height;
            byte[] buffer = new byte[Math.Abs(strideLength)];
            Marshal.Copy(bmpData.Scan0, buffer, 0, strideLength);
            bmp.UnlockBits(bmpData);


            Bitmap outBmp = new Bitmap(width, height, PixelFormat.Format8bppIndexed);
            ColorPalette pal = outBmp.Palette;
            for (int i = 0; i < 256; i++)
                pal.Entries[i] = _grayPallette[i];
            outBmp.Palette = pal;

            BitmapData outBmpData = outBmp.LockBits(rect, ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);
            int outStride = outBmpData.Stride;

            int outLength = outStride * height;
            byte[] pixelBuffer = new byte[outLength];
            byte[] resultBuffer = new byte[outLength];

            //0.299R + 0.587G + 0.114B
            for (int row = 0; row < height; row++)
            {
                int inRowOffset = inStride * row;
                int outRowOffset = outStride * row;
                for (int col = 0; col < width; col++)
                {
                    int offset = inRowOffset + col * 3;
                    byte brightness = (byte)(buffer[offset] * 0.114f + buffer[offset + 1] * 0.587f + buffer[offset + 2] * 0.299f);
                    pixelBuffer[outRowOffset + col] = brightness;
                }
            }

            for (int row = 1; row < height - 1; row++)
            {
                int rowOffset = outStride * row;
                for (int col = 1; col < width - 1; col++)
                {
                    int offset = rowOffset + col;

                    int gx = -pixelBuffer[offset - 1 - outStride] + pixelBuffer[offset + 1 - outStride] - 2 * pixelBuffer[offset - 1] +
                        2 * pixelBuffer[offset + 1] - pixelBuffer[offset - 1 + outStride] + pixelBuffer[offset + 1 + outStride];

                    int gy = pixelBuffer[offset - 1 - outStride] + 2 * pixelBuffer[offset - outStride] + pixelBuffer[offset + 1 - outStride] -
                        pixelBuffer[offset - 1 + outStride] - 2 * pixelBuffer[offset + outStride] - pixelBuffer[offset + 1 + outStride];

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

                    resultBuffer[offset] = (byte)g;
                }
            }
            Marshal.Copy(resultBuffer, 0, outBmpData.Scan0, outLength);
            outBmp.UnlockBits(outBmpData);
            return outBmp;
        }
    }
}

Benchmark.NET

Method	Mean	Error	StdDev	Gen 0	Gen 1	Gen 2	Allocated
'SIMD Enabled'	6.318 ms	0.1190 ms	0.1113 ms	992.1875	992.1875	992.1875	3.35 MB
'SIMD Disabled'	43.062 ms	0.2253 ms	0.2107 ms	500.0000	500.0000	500.0000	16.61 MB