WindowsFormsTest/SauvolaTest/SauvolaBinarization.cs

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Drawing;
using System.Drawing.Imaging;
using System.Linq;
using System.Runtime.InteropServices;
using System.Text;
using System.Threading.Tasks;

namespace SauvolaTest
{
    public static class SauvolaBinarization
    {
        /// <summary>
        /// Sauvola 二值化算法实现
        /// </summary>
        public static class SauvolaBinarizer
        {
            /// <summary>
            /// 应用 Sauvola 二值化
            /// </summary>
            /// <param name="source">8bpp 灰度源图像</param>
            /// <param name="windowSize">局部窗口大小 (奇数)</param>
            /// <param name="k">修正系数 (0.2-0.5)</param>
            /// <param name="r">动态范围 (通常 128)</param>
            public static Bitmap Apply(Bitmap source, int windowSize = 25, double k = 0.2, double r = 128.0)
            {
                if (source == null) throw new ArgumentNullException(nameof(source));
                if (source.PixelFormat != PixelFormat.Format8bppIndexed)
                {
                    throw new ArgumentException("输入图像必须是 8bpp 灰度格式。请先调用 ConvertToGrayscaleFast。", nameof(source));
                }

                int width = source.Width;
                int height = source.Height;
                Bitmap result = new Bitmap(width, height, PixelFormat.Format8bppIndexed);

                // 复制调色板
                result.Palette = source.Palette;

                BitmapData srcData = source.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.ReadOnly, PixelFormat.Format8bppIndexed);
                BitmapData resData = result.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);

                try
                {
                    unsafe
                    {
                        byte* srcPtr = (byte*)srcData.Scan0;
                        byte* resPtr = (byte*)resData.Scan0;
                        int stride = srcData.Stride;

                        // 1. 构建积分图 (Integral Image) 和 平方积分图
                        // 使用 long 防止溢出 (最大像素和: 255 * 10^6 约等于 2.5*10^8, int 足够，但平方和需要 long)
                        // 尺寸宽+1, 高+1 以便处理边界
                        long[] integral = new long[(width + 1) * (height + 1)];
                        long[] squaredIntegral = new long[(width + 1) * (height + 1)];
                        Stopwatch sw = Stopwatch.StartNew();
                        for (int y = 0; y < height; y++)
                        {
                            long rowSum = 0;
                            long rowSqSum = 0;
                            for (int x = 0; x < width; x++)
                            {
                                byte pixel = srcPtr[y * stride + x];
                                rowSum += pixel;
                                rowSqSum += (long)pixel * pixel;

                                // I(y+1, x+1) = I(y, x+1) + RowSum
                                int currentIndex = (y + 1) * (width + 1) + (x + 1);
                                int topIndex = y * (width + 1) + (x + 1);

                                integral[currentIndex] = integral[topIndex] + rowSum;
                                squaredIntegral[currentIndex] = squaredIntegral[topIndex] + rowSqSum;
                            }
                        }
                        sw.Stop();
                        Console.WriteLine("积分图耗时：{0}", sw.ElapsedMilliseconds);

                        int halfWin = windowSize / 2;
                        Stopwatch Sauvolastopwatch = new Stopwatch();
                        // 2. 计算阈值并二值化
                        for (int y = 0; y < height; y++)
                        {
                            for (int x = 0; x < width; x++)
                            {
                                // 窗口边界
                                int x1 = Math.Max(0, x - halfWin);
                                int y1 = Math.Max(0, y - halfWin);
                                int x2 = Math.Min(width - 1, x + halfWin);
                                int y2 = Math.Min(height - 1, y + halfWin);

                                // 窗口内像素数
                                int count = (x2 - x1 + 1) * (y2 - y1 + 1);

                                // 从积分图获取区域和
                                long sum = GetRegionSum(integral, x1, y1, x2, y2, width);
                                long sqSum = GetRegionSum(squaredIntegral, x1, y1, x2, y2, width);

                                // 均值
                                double mean = (double)sum / count;

                                // 方差 = E[X^2] - (E[X])^2
                                double variance = ((double)sqSum / count) - (mean * mean);
                                if (variance < 0) variance = 0; // 防止浮点误差
                                double stdDev = Math.Sqrt(variance);

                                // Sauvola 公式
                                double threshold = mean * (1 + k * ((stdDev / r) - 1));

                                // 比较
                                byte currentPixel = srcPtr[y * stride + x];
                                resPtr[y * stride + x] = (byte)(currentPixel > threshold ? 255 : 0);
                            }
                        }
                        Sauvolastopwatch.Stop();
                        Console.WriteLine("二值化耗时：{0}", Sauvolastopwatch.ElapsedMilliseconds);
                    }
                }
                finally
                {
                    source.UnlockBits(srcData);
                    result.UnlockBits(resData);
                }

                return result;
            }

            private static long GetRegionSum(long[] integralImg, int x1, int y1, int x2, int y2, int width)
            {
                // 积分图坐标需要 +1
                int w = width + 1;
                // A - B - C + D
                // A: (y2+1, x2+1)
                // B: (y1, x2+1)
                // C: (y2+1, x1)
                // D: (y1, x1)

                long a = integralImg[(y2 + 1) * w + (x2 + 1)];
                long b = integralImg[y1 * w + (x2 + 1)];
                long c = integralImg[(y2 + 1) * w + x1];
                long d = integralImg[y1 * w + x1];

                return a - b - c + d;
            }

            /// <summary>
            /// 高效串行版 Sauvola 二值化 (使用 fixed 指针，无 Parallel 限制)
            /// </summary>
            public static Bitmap ApplyFast(Bitmap source, int windowSize = 25, double k = 0.2, double r = 128.0)
            {
                if (source == null) throw new ArgumentNullException(nameof(source));
                if (source.PixelFormat != PixelFormat.Format8bppIndexed)
                {
                    throw new ArgumentException("输入图像必须是 8bpp 灰度格式。", nameof(source));
                }

                int width = source.Width;
                int height = source.Height;
                Bitmap result = new Bitmap(width, height, PixelFormat.Format8bppIndexed);
                result.Palette = source.Palette;

                BitmapData srcData = source.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.ReadOnly, PixelFormat.Format8bppIndexed);
                BitmapData resData = result.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);

                try
                {
                    unsafe
                    {
                        byte* srcPtr = (byte*)srcData.Scan0;
                        byte* resPtr = (byte*)resData.Scan0;
                        int stride = srcData.Stride;

                        // 预计算常数
                        double kDivR = k / r;
                        int halfWin = windowSize / 2;
                        int wPlus1 = width + 1;

                        // 1. 构建积分图
                        int[] integral = new int[(width + 1) * (height + 1)];
                        long[] squaredIntegral = new long[(width + 1) * (height + 1)];

                        // 固定数组以使用指针加速构建
                        fixed (int* pInt = integral)
                        fixed (long* pSq = squaredIntegral)
                        {
                            for (int y = 0; y < height; y++)
                            {
                                int rowSum = 0;
                                long rowSqSum = 0;

                                // 指针定位到当前行和上一行的起始位置
                                // 积分图每行有 width + 1 个元素
                                int* prevRowInt = pInt + y * wPlus1;
                                int* currRowInt = pInt + (y + 1) * wPlus1;

                                long* prevRowSq = pSq + y * wPlus1;
                                long* currRowSq = pSq + (y + 1) * wPlus1;

                                byte* rowSrc = srcPtr + y * stride;

                                for (int x = 0; x < width; x++)
                                {
                                    byte pixel = rowSrc[x];
                                    rowSum += pixel;
                                    rowSqSum += (long)pixel * pixel;

                                    // 指针偏移 x+1
                                    currRowInt[x + 1] = prevRowInt[x + 1] + rowSum;
                                    currRowSq[x + 1] = prevRowSq[x + 1] + rowSqSum;
                                }
                            }
                        }

                        // 2. 串行二值化 (使用 fixed 指针加速读取)
                        // 对于大多数中等分辨率图像，串行+指针比 Parallel+数组索引更快或相当
                        fixed (int* pInt = integral)
                        fixed (long* pSq = squaredIntegral)
                        {
                            for (int y = 0; y < height; y++)
                            {
                                byte* rowSrc = srcPtr + y * stride;
                                byte* rowRes = resPtr + y * stride;

                                for (int x = 0; x < width; x++)
                                {
                                    // 窗口边界
                                    int x1 = x - halfWin;
                                    int y1 = y - halfWin;
                                    int x2 = x + halfWin;
                                    int y2 = y + halfWin;

                                    // 边界裁剪
                                    if (x1 < 0) x1 = 0;
                                    if (y1 < 0) y1 = 0;
                                    if (x2 >= width) x2 = width - 1;
                                    if (y2 >= height) y2 = height - 1;

                                    // 计算窗口内像素数
                                    int w = x2 - x1 + 1;
                                    int h = y2 - y1 + 1;
                                    int count = w * h;

                                    // 从积分图获取区域和 (指针访问)
                                    // 索引: (y * wPlus1) + x
                                    int idxA = (y2 + 1) * wPlus1 + (x2 + 1);
                                    int idxB = y1 * wPlus1 + (x2 + 1);
                                    int idxC = (y2 + 1) * wPlus1 + x1;
                                    int idxD = y1 * wPlus1 + x1;

                                    int sum = pInt[idxA] - pInt[idxB] - pInt[idxC] + pInt[idxD];
                                    long sqSum = pSq[idxA] - pSq[idxB] - pSq[idxC] + pSq[idxD];

                                    // 计算均值和方差
                                    double mean = (double)sum / count;
                                    double variance = ((double)sqSum / count) - (mean * mean);

                                    if (variance < 0) variance = 0;

                                    double stdDev = Math.Sqrt(variance);

                                    // Sauvola 阈值公式
                                    double threshold = mean * (1.0 + kDivR * stdDev - k);

                                    // 二值化
                                    rowRes[x] = (byte)(rowSrc[x] > threshold ? 255 : 0);
                                }
                            }
                        }
                    }
                }
                finally
                {
                    source.UnlockBits(srcData);
                    result.UnlockBits(resData);
                }

                return result;
            }

            /// <summary>
            /// 高效串行版 Sauvola 二值化 (使用 fixed 指针，无 Parallel 限制)
            /// </summary>
            public static Bitmap ApplyParallelGCHandle(Bitmap source, int windowSize = 25, double k = 0.2, double r = 128.0)
            {
                if (source == null) throw new ArgumentNullException(nameof(source));
                if (source.PixelFormat != PixelFormat.Format8bppIndexed)
                {
                    throw new ArgumentException("输入图像必须是 8bpp 灰度格式。", nameof(source));
                }

                int width = source.Width;
                int height = source.Height;
                Bitmap result = new Bitmap(width, height, PixelFormat.Format8bppIndexed);
                result.Palette = source.Palette;

                BitmapData srcData = source.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.ReadOnly, PixelFormat.Format8bppIndexed);
                BitmapData resData = result.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);

                try
                {
                    unsafe
                    {
                        byte* srcPtr = (byte*)srcData.Scan0;
                        byte* resPtr = (byte*)resData.Scan0;
                        int stride = srcData.Stride;

                        // 预计算常数
                        double kDivR = k / r;
                        int halfWin = windowSize / 2;
                        int wPlus1 = width + 1;

                        // 1. 构建积分图
                        int[] integral = new int[(width + 1) * (height + 1)];
                        long[] squaredIntegral = new long[(width + 1) * (height + 1)];

                        // 固定数组以使用指针加速构建
                        fixed (int* pInt = integral)
                        fixed (long* pSq = squaredIntegral)
                        {
                            for (int y = 0; y < height; y++)
                            {
                                int rowSum = 0;
                                long rowSqSum = 0;

                                // 指针定位到当前行和上一行的起始位置
                                // 积分图每行有 width + 1 个元素
                                int* prevRowInt = pInt + y * wPlus1;
                                int* currRowInt = pInt + (y + 1) * wPlus1;

                                long* prevRowSq = pSq + y * wPlus1;
                                long* currRowSq = pSq + (y + 1) * wPlus1;

                                byte* rowSrc = srcPtr + y * stride;

                                for (int x = 0; x < width; x++)
                                {
                                    byte pixel = rowSrc[x];
                                    rowSum += pixel;
                                    rowSqSum += (long)pixel * pixel;

                                    // 指针偏移 x+1
                                    currRowInt[x + 1] = prevRowInt[x + 1] + rowSum;
                                    currRowSq[x + 1] = prevRowSq[x + 1] + rowSqSum;
                                }
                            }
                        }

                        // 2. 串行二值化 (使用 fixed 指针加速读取)
                        // 对于大多数中等分辨率图像，串行+指针比 Parallel+数组索引更快或相当
                        fixed (int* pInt = integral)
                        fixed (long* pSq = squaredIntegral)
                        {
                            for (int y = 0; y < height; y++)
                            {
                                byte* rowSrc = srcPtr + y * stride;
                                byte* rowRes = resPtr + y * stride;

                                for (int x = 0; x < width; x++)
                                {
                                    // 窗口边界
                                    int x1 = x - halfWin;
                                    int y1 = y - halfWin;
                                    int x2 = x + halfWin;
                                    int y2 = y + halfWin;

                                    // 边界裁剪
                                    if (x1 < 0) x1 = 0;
                                    if (y1 < 0) y1 = 0;
                                    if (x2 >= width) x2 = width - 1;
                                    if (y2 >= height) y2 = height - 1;

                                    // 计算窗口内像素数
                                    int w = x2 - x1 + 1;
                                    int h = y2 - y1 + 1;
                                    int count = w * h;

                                    // 从积分图获取区域和 (指针访问)
                                    // 索引: (y * wPlus1) + x
                                    int idxA = (y2 + 1) * wPlus1 + (x2 + 1);
                                    int idxB = y1 * wPlus1 + (x2 + 1);
                                    int idxC = (y2 + 1) * wPlus1 + x1;
                                    int idxD = y1 * wPlus1 + x1;

                                    int sum = pInt[idxA] - pInt[idxB] - pInt[idxC] + pInt[idxD];
                                    long sqSum = pSq[idxA] - pSq[idxB] - pSq[idxC] + pSq[idxD];

                                    // 计算均值和方差
                                    double mean = (double)sum / count;
                                    double variance = ((double)sqSum / count) - (mean * mean);

                                    if (variance < 0) variance = 0;

                                    double stdDev = Math.Sqrt(variance);

                                    // Sauvola 阈值公式
                                    double threshold = mean * (1.0 + kDivR * stdDev - k);

                                    // 二值化
                                    rowRes[x] = (byte)(rowSrc[x] > threshold ? 255 : 0);
                                }
                            }
                        }
                    }
                }
                finally
                {
                    source.UnlockBits(srcData);
                    result.UnlockBits(resData);
                }

                return result;
            }

            /// <summary>
            /// 并行版 Sauvola 二值化 (使用 GCHandle 解决 fixed 指针在 Lambda 中的限制)
            /// </summary>
            /// <param name="source">8bpp 灰度源图像</param>
            /// <param name="windowSize">局部窗口大小 (奇数)</param>
            /// <param name="k">修正系数 (0.2-0.5)</param>
            /// <param name="r">动态范围 (通常 128)</param>
            public static Bitmap ApplyParallel(Bitmap source, int windowSize = 25, double k = 0.2, double r = 128.0)
            {
                if (source == null) throw new ArgumentNullException(nameof(source));
                if (source.PixelFormat != PixelFormat.Format8bppIndexed)
                {
                    throw new ArgumentException("输入图像必须是 8bpp 灰度格式。请先调用 ConvertToGrayscaleFast。", nameof(source));
                }

                int width = source.Width;
                int height = source.Height;
                Bitmap result = new Bitmap(width, height, PixelFormat.Format8bppIndexed);

                // 复制调色板
                result.Palette = source.Palette;

                BitmapData srcData = source.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.ReadOnly, PixelFormat.Format8bppIndexed);
                BitmapData resData = result.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);

                // 预计算常数
                double kDivR = k / r;
                int halfWin = windowSize / 2;
                int wPlus1 = width + 1; // 积分图宽度

                // 分配积分图数组
                // 使用 int 存储普通积分图以优化缓存 (假设图像总和不超过 int.MaxValue)
                int[] integral = new int[wPlus1 * (height + 1)];
                // 平方积分图必须使用 long
                long[] squaredIntegral = new long[wPlus1 * (height + 1)];

                try
                {
                    unsafe
                    {
                        byte* srcPtr = (byte*)srcData.Scan0;
                        byte* resPtr = (byte*)resData.Scan0;
                        int stride = srcData.Stride;

                        // ---------------------------------------------------------
                        // 1. 构建积分图 (串行)
                        // ---------------------------------------------------------
                        // 这里可以使用 fixed，因为不在 Parallel 内部
                        fixed (int* pInt = integral)
                        fixed (long* pSq = squaredIntegral)
                        {
                            for (int y = 0; y < height; y++)
                            {
                                int rowSum = 0;
                                long rowSqSum = 0;

                                // 指针定位到当前行和上一行
                                int* prevRowInt = pInt + y * wPlus1;
                                int* currRowInt = pInt + (y + 1) * wPlus1;

                                long* prevRowSq = pSq + y * wPlus1;
                                long* currRowSq = pSq + (y + 1) * wPlus1;

                                byte* rowSrc = srcPtr + y * stride;

                                for (int x = 0; x < width; x++)
                                {
                                    byte pixel = rowSrc[x];
                                    rowSum += pixel;
                                    rowSqSum += (long)pixel * pixel;

                                    // 积分图递推公式
                                    currRowInt[x + 1] = prevRowInt[x + 1] + rowSum;
                                    currRowSq[x + 1] = prevRowSq[x + 1] + rowSqSum;
                                }
                            }
                        }

                        // ---------------------------------------------------------
                        // 2. 并行二值化 (使用 GCHandle)
                        // ---------------------------------------------------------

                        // 使用 GCHandle 固定数组内存，防止 GC 移动它们
                        // 这样获取的指针可以在 Parallel.For 的 Lambda 中安全使用
                        GCHandle hIntegral = GCHandle.Alloc(integral, GCHandleType.Pinned);
                        GCHandle hSqIntegral = GCHandle.Alloc(squaredIntegral, GCHandleType.Pinned);

                        try
                        {
                            // 获取固定内存的地址指针
                            int* pIntFixed = (int*)hIntegral.AddrOfPinnedObject();
                            long* pSqFixed = (long*)hSqIntegral.AddrOfPinnedObject();

                            // 并行处理每一行
                            Parallel.For(0, height, y =>
                            {
                                byte* rowSrc = srcPtr + y * stride;
                                byte* rowRes = resPtr + y * stride;

                                for (int x = 0; x < width; x++)
                                {
                                    // 计算窗口边界
                                    int x1 = x - halfWin;
                                    int y1 = y - halfWin;
                                    int x2 = x + halfWin;
                                    int y2 = y + halfWin;

                                    // 边界裁剪
                                    if (x1 < 0) x1 = 0;
                                    if (y1 < 0) y1 = 0;
                                    if (x2 >= width) x2 = width - 1;
                                    if (y2 >= height) y2 = height - 1;

                                    // 窗口像素数量
                                    int w = x2 - x1 + 1;
                                    int h = y2 - y1 + 1;
                                    int count = w * h;

                                    // 计算积分图索引
                                    // 公式: Index = y * wPlus1 + x
                                    int idxA = (y2 + 1) * wPlus1 + (x2 + 1);
                                    int idxB = y1 * wPlus1 + (x2 + 1);
                                    int idxC = (y2 + 1) * wPlus1 + x1;
                                    int idxD = y1 * wPlus1 + x1;

                                    // 获取区域和 (指针访问，无边界检查开销)
                                    int sum = pIntFixed[idxA] - pIntFixed[idxB] - pIntFixed[idxC] + pIntFixed[idxD];
                                    long sqSum = pSqFixed[idxA] - pSqFixed[idxB] - pSqFixed[idxC] + pSqFixed[idxD];

                                    // 计算均值
                                    double mean = (double)sum / count;

                                    // 计算方差: E[X^2] - (E[X])^2
                                    double variance = ((double)sqSum / count) - (mean * mean);

                                    // 防止浮点误差导致负数
                                    if (variance < 0) variance = 0;

                                    // 标准差
                                    double stdDev = Math.Sqrt(variance);

                                    // Sauvola 阈值公式: T = m * (1 + k * (s/r - 1))
                                    // 优化为: T = m * (1 + (k/r)*s - k)
                                    double threshold = mean * (1.0 + kDivR * stdDev - k);

                                    // 二值化判定
                                    rowRes[x] = (byte)(rowSrc[x] > threshold ? 255 : 0);
                                }
                            });
                        }
                        finally
                        {
                            // 必须释放 GCHandle，否则会导致内存泄漏
                            if (hIntegral.IsAllocated) hIntegral.Free();
                            if (hSqIntegral.IsAllocated) hSqIntegral.Free();
                        }
                    }
                }
                finally
                {
                    source.UnlockBits(srcData);
                    result.UnlockBits(resData);
                }

                return result;
            }


            // 可选：快速近似平方根 (如果 Math.Sqrt 仍太慢)
            private static double ApproxSqrt(double x)
            {
                if (x <= 0) return 0;
                // 牛顿迭代法一次迭代，或者使用位操作黑客技巧
                // 这里简单返回 Math.Sqrt，因为在 .NET Core/.NET 5+ 中 Math.Sqrt 已经非常优化
                // 如果在 .NET Framework 且极追求性能，可替换为特定算法
                return Math.Sqrt(x);
            }
        }
    }
}