public static byte[,] SauvolaBinarize(byte[,] grayScale, int WHalf, double K)
{
        var integrate = 
            from g in grayScale
            from ri in Result.SubstBy(0)    // сумма элементов - int32
            from rq in Result.SubstBy(0L)    // сумма квадратов - int64. Можно обойтись и int32, если площадь окна просмотра не более 2^16, но так честнее
            select ValueTuple.Create(
                ri[-1, 0] + ri[0, -1] - ri[-1, -1] + g,
                rq[-1, 0] + rq[0, -1] - rq[-1, -1] + g * g);
        var (p, sq) = integrate.ToArrays();    // вычисляем оба интеграла сразу

    var detect = 
            from g in grayScale
            from i in p.With(OutOfBoundsStrategy.Integral(0))    // специальная стратегия подстановки для значений за пределами массива - 0 для слева-сверху,
            from q in sq.With(OutOfBoundsStrategy.Integral(0L)) // ближайший сосед для справа-снизу.
            let tl = i.Offset(-WHalf - 1, -WHalf - 1)        // самая громоздкая часть: вычисляем эффективное окно просмотра, смещая текущую точку на +-Whalf
            let tr = i.Offset(-WHalf - 1, WHalf)        // нам нужны 8 точек, т.к. мы используем оба предынтегрированных массива
            let bl = i.Offset(WHalf, -WHalf - 1)
            let br = i.Offset(WHalf, WHalf)
            let tlq = q.Offset(-WHalf - 1, -WHalf - 1)
            let trq = q.Offset(-WHalf - 1, WHalf)
            let blq = q.Offset(WHalf, -WHalf - 1)
            let brq = q.Offset(WHalf, WHalf)
            let area = (br.X - tl.X) * (br.Y - tl.Y)        // площадь окна разная, в зависимости от того, куда попала точка; все поправки в X и Y учитываются OutOfBoundsStrategy.Integral
            let diff = br.Value + tl.Value - tr.Value - bl.Value
            let sqdiff = brq.Value + tlq.Value - trq.Value - blq.Value
            let mean = (double)diff / area            // вообще, здесь можно и не приводить к double. Но это мешает векторизации, т.к. целочисленный div есть только в AVX-512. 
            let std = Math.Sqrt((sqdiff - diff * mean) / (area - 1)) // А вот double делятся уже в avx2.

            let threshold = mean * (1 + K * ((std / 128) - 1))

            select g > threshold ? byte.MaxValue : byte.MinValue;
}

Intel Core i7-6600U CPU 2.60GHz (Skylake), 1 CPU, 4 logical and 2 physical cores
.NET Core SDK=3.1.302
[Host]     : .NET Core 3.1.6 (CoreCLR 4.700.20.26901, CoreFX 4.700.20.31603), X64 RyuJIT
|                  Method | WHalf |   FileName |    Mean |    Error |   StdDev |  Median | Ratio | RatioSD |
|------------------------ |------ |----------- |--------:|---------:|---------:|--------:|------:|--------:|
|             SafeSauvola |     5 | p00743.bmp | 2.795 s | 0.0624 s | 0.1812 s | 2.744 s |  1.91 |    0.20 |
|     UnsafeSauvolaScalar |     5 | p00743.bmp | 1.475 s | 0.0494 s | 0.1455 s | 1.410 s |  1.00 |    0.00 |
|       LinqSauvolaVector |     5 | p00743.bmp | 2.433 s | 0.0485 s | 0.0664 s | 2.439 s |  1.47 |    0.09 |
|       LinqSauvolaScalar |     5 | p00743.bmp | 2.598 s | 0.1468 s | 0.4328 s | 2.607 s |  1.77 |    0.34 |
| CachedLinqSauvolaVector |     5 | p00743.bmp | 1.018 s | 0.0202 s | 0.0283 s | 1.014 s |  0.62 |    0.04 |
| CachedLinqSauvolaScalar |     5 | p00743.bmp | 1.140 s | 0.0226 s | 0.0338 s | 1.139 s |  0.69 |    0.04 |

    public static unsafe byte[,] UnsafeSauvola(byte[,] grayScale, int whalf)
        {
            int height = grayScale.Height();
            int width = grayScale.Width();
            var p = new int[height, width];
            var sq = new long[height, width];
            byte[,] result = new byte[height, width];

            fixed (byte* gPtr = &grayScale[0, 0])
            fixed (int* pPtr = &p[0, 0])
            fixed (long* sqPtr = &sq[0, 0])
            fixed (byte* rPtr = &result[0, 0])
            {
                for (int i = 0; i < height; i++)
                    for (int j = 0; j < width; j++)
                    {
                        pPtr[i * width + j] = gPtr[i * width + j];
                        sqPtr[i * width + j] = gPtr[i * width + j] * gPtr[i * width + j];
                        if (i > 0)
                        {
                            pPtr[i * width + j] += pPtr[(i - 1) * width + j];
                            sqPtr[i * width + j] += sqPtr[(i - 1) * width + j];
                        }
                        if (j > 0)
                        {
                            pPtr[i * width + j] += pPtr[i * width + j - 1];
                            sqPtr[i * width + j] += sqPtr[i * width + j - 1];
                        }
                        if (i > 0 && j > 0)
                        {
                            pPtr[i * width + j] -= pPtr[(i - 1) * width + j - 1];
                            sqPtr[i * width + j] -= sqPtr[(i - 1) * width + j - 1];
                        }
                    }
                for (int i = 0; i < height; i++)
                    for (int j = 0; j < width; j++)
                    {
                        var xmin = Math.Max(0, i - whalf);
                        var ymin = Math.Max(0, j - whalf);
                        var xmax = Math.Min(height - 1, i + whalf);
                        var ymax = Math.Min(width - 1, j + whalf);

                        var area = (xmax - xmin + 1) * (ymax - ymin + 1);

                        var diff = pPtr[width * xmax + ymax];
                        var sqdiff = sqPtr[width * xmax + ymax];
                        if (xmin > 0)
                        {
                            diff -= -pPtr[width * (xmin - 1) + ymax];
                            sqdiff -= -sqPtr[width * (xmin - 1) + ymax];
                        }
                        if (ymin > 0)
                        {
                            diff -= pPtr[width * xmax + ymin - 1];
                            sqdiff -= sqPtr[width * xmax + ymin - 1];
                        }
                        if (xmin > 0 && ymin > 0)
                        {
                            diff += pPtr[width * (xmin - 1) + ymin - 1];
                            sqdiff += sqPtr[width * (xmin - 1) + ymin - 1];
                        }
                        var mean = (double)diff / area;
                        var std = Math.Sqrt((sqdiff - diff * mean) / (area - 1));

                        var threshold = mean * (1 + K * ((std / 128) - 1));

                        rPtr[i * width + j] = gPtr[i * width + j] > threshold ? byte.MaxValue : byte.MinValue;
                    }
                return result;
            }
        }