C# Susan边缘检测(Susan Edge Detection)

Susan边缘检测，方法简单，效率高，具体参照 The SUSAN Edge Detector in Detail， 修改dThreshold值，可以改变检测效果，用参照提供的重心法、力矩法可得到边缘方向；

        /// https://users.fmrib.ox.ac.uk/~steve/susan/susan/node6.html
        public unsafe static Bitmap SusanGray(this Bitmap sourceBitmap)
        {
            int[] rowRadius = new int[7] { 1, 2, 3, 3, 3, 2, 1 };
            int width = sourceBitmap.Width;
            int height = sourceBitmap.Height;
            BitmapData sourceData = sourceBitmap.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.ReadOnly, PixelFormat.Format32bppArgb);

            int stride = sourceData.Stride;
            byte[] pixelBuffer = new byte[stride * sourceData.Height];
            byte[] resultBuffer = new byte[stride * sourceData.Height];

            Marshal.Copy(sourceData.Scan0, pixelBuffer, 0, pixelBuffer.Length);
            sourceBitmap.UnlockBits(sourceData);

            float rgb = 0;

            for (int k = 0; k < pixelBuffer.Length; k += 4)
            {
                rgb = pixelBuffer[k] * 0.11f;
                rgb += pixelBuffer[k + 1] * 0.59f;
                rgb += pixelBuffer[k + 2] * 0.3f;


                pixelBuffer[k] = (byte)rgb;
                pixelBuffer[k + 1] = pixelBuffer[k];
                pixelBuffer[k + 2] = pixelBuffer[k];
                pixelBuffer[k + 3] = 255;
            }

            int[] susanMap = new int[height * width];

            int offset = stride - width * 4;

            GCHandle srchandle = GCHandle.Alloc(susanMap, GCHandleType.Pinned);
            IntPtr susan = srchandle.AddrOfPinnedObject();

            int dThreshold = 28;
            fixed (byte* pbuff = pixelBuffer, rbuff = resultBuffer)
            {
                byte* src = pbuff + stride * 3 + 3 * 4;
                int* br = (int*)susan + height * 3 + 3;
                byte* dst = rbuff + stride * 3 + 3 * 4;

                for (int offsetY = 3; offsetY < height - 3; offsetY++)
                {
                    for (int offsetX = 3; offsetX < width - 3; offsetX++, src += 4,dst+=4, br++)
                    {
                        byte nucleusValue = *src;
                        int usan = 0;

                        int cx = 0, cy = 0;

                        for (int i = -3; i <= 3; i++)
                        {

                            int r = rowRadius[i + 3];

                            byte* ptr = (byte*)((int)src + stride * i);

                            for (int j = -r; j <= r; j++)
                            {
                                int c = (int)Math.Exp(-Math.Pow((System.Math.Abs(nucleusValue - ptr[j * 4]) / dThreshold), 6));
                                usan += c;
                                cx += j * c;
                                cy += i * c;
                            }
                        }
                        if (usan < 28)
                            usan = 28 -usan;
                        else
                            usan = 0;
                        if ((usan < 6) && (cx != 0 || cy != 0))
                        {
                            *dst = 255;
                            dst[1] = 255;
                            dst[2] = 255;
                            dst[3] = 255;
                        }
                        else
                        {
                            *dst = 0;
                            dst[1] = 0;
                            dst[2] = 0;
                            dst[3] = 255;
                        }
                        *br = usan;
                    }
                    src += 6 * 4 + offset;
                    dst += 6 * 4 + offset;
                    br += 6;
                }
            }

            Bitmap resultBitmap = new Bitmap(sourceBitmap.Width, sourceBitmap.Height);

            BitmapData resultData = resultBitmap.LockBits(new Rectangle(0, 0,
                                     resultBitmap.Width, resultBitmap.Height),
                                                      ImageLockMode.WriteOnly,
                                                 PixelFormat.Format32bppArgb);

            Marshal.Copy(resultBuffer, 0, resultData.Scan0, resultBuffer.Length);
            resultBitmap.UnlockBits(resultData);

            return resultBitmap;

        }

 并行的方法：

        public unsafe static Bitmap ParallelSusan(this Bitmap sourceBitmap)
        {
            int width = sourceBitmap.Width;
            int height = sourceBitmap.Height;
            BitmapData sourceData = sourceBitmap.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.ReadOnly, PixelFormat.Format32bppArgb);

            int stride = sourceData.Stride;
            byte[] pixelBuffer = new byte[stride * sourceData.Height];
            byte[] resultBuffer = new byte[stride * sourceData.Height];

            Marshal.Copy(sourceData.Scan0, pixelBuffer, 0, pixelBuffer.Length);
            sourceBitmap.UnlockBits(sourceData);

            float rgb = 0;

            for (int k = 0; k < pixelBuffer.Length; k += 4)
            {
                rgb = pixelBuffer[k] * 0.11f;
                rgb += pixelBuffer[k + 1] * 0.59f;
                rgb += pixelBuffer[k + 2] * 0.3f;

                pixelBuffer[k] = (byte)rgb;
                pixelBuffer[k + 1] = pixelBuffer[k];
                pixelBuffer[k + 2] = pixelBuffer[k];
                pixelBuffer[k + 3] = 255;
            }

            int[] susanMap = new int[height * width];

            int offset = stride - width * 4;
            GCHandle srchandle = GCHandle.Alloc(pixelBuffer, GCHandleType.Pinned);
            IntPtr src = srchandle.AddrOfPinnedObject();

            GCHandle dsthandle = GCHandle.Alloc(resultBuffer, GCHandleType.Pinned);
            IntPtr dst = dsthandle.AddrOfPinnedObject();

            GCHandle suhandle = GCHandle.Alloc(susanMap, GCHandleType.Pinned);
            IntPtr susan = suhandle.AddrOfPinnedObject();

                System.Threading.Tasks.Parallel.For(3, height - 3, (offsetY) =>
                {
                    for (int offsetX = 3; offsetX < width - 3; offsetX++)
                    {
                        OneSusan(offsetY, offsetX, (byte*)src, (byte*)dst, stride);
                    }
                });

            Bitmap resultBitmap = new Bitmap(sourceBitmap.Width, sourceBitmap.Height);

            BitmapData resultData = resultBitmap.LockBits(new Rectangle(0, 0,
                                     resultBitmap.Width, resultBitmap.Height),
                                                      ImageLockMode.WriteOnly,
                                                 PixelFormat.Format32bppArgb);

            Marshal.Copy(resultBuffer, 0, resultData.Scan0, resultBuffer.Length);
            resultBitmap.UnlockBits(resultData);

            return resultBitmap;

        }
        public unsafe static void OneSusan(int offsetY, int offsetX, byte* src, byte* dst, int stride)
        {
            int[] rowRadius = new int[7] { 1, 2, 3, 3, 3, 2, 1 };
            int dThreshold = 28;

            src = (byte*)((int)src + stride * offsetY + offsetX * 4);
            dst = (byte*)((int)dst + stride * offsetY + offsetX * 4);
            byte nucleusValue = *src;
            int usan = 0;

            int cx = 0, cy = 0;
            float vX = 0, vY = 0, vXY = 0;
            for (int i = -3; i <= 3; i++)
            {

                int r = rowRadius[i + 3];

                byte* ptr = (byte*)((int)src + stride * i);

                for (int j = -r; j <= r; j++)
                {
                    int c = (int)Math.Exp(-Math.Pow((System.Math.Abs(nucleusValue - ptr[j * 4]) / dThreshold), 6));
                    usan += c;
                    cx += j * c;
                    cy += i * c;
                    vX += j * j * c;
                    vY += i * i * c;
                    vXY += i * j * c;
                }
            }
            if (usan < 28)
                usan = 28 - usan;
            else
                usan = 0;
            if ((usan < 5) && (cx != 0 || cy != 0))
            {
                *dst = 255;
                dst[1] = 255;
                dst[2] = 255;
                dst[3] = 255;
            }
            else
            {
                *dst = 0;
                dst[1] = 0;
                dst[2] = 0;
                dst[3] = 255;
            }
        }

 示例下载（除Susan 方法之外的代码来自https://softwarebydefault.com/2013/05/11/image-edge-detection/)