C# 实现gabor滤波

        string ImagePath;
        private void Button1_Click(object sender, EventArgs e)
        {
            OpenFileDialog openfile = new OpenFileDialog();

            if (openfile.ShowDialog() == DialogResult.OK)
            {
                ImagePath = openfile.FileName;
                pictureBox1.ImageLocation=ImagePath;
            }
        }

        private void Button2_Click(object sender, EventArgs e)
        {
            Image<Gray, byte> src = new Image<Gray, byte>(ImagePath); ;
            Image<Gray, byte> dst1 = new Image<Gray, byte>(src.Width, src.Height);
            Image<Gray, byte> dst2 = new Image<Gray, byte>(src.Width, src.Height);
            Image<Gray, byte> dst3 = new Image<Gray, byte>(src.Width, src.Height);
            GaborFilter gf = new GaborFilter();
            //0度
            Bitmap bm1 = gf.ProcessFilter(UnmanagedImage.FromManagedImage(src.Bitmap), dst1.Bitmap, 3, 5, 0, 0.3, 2, 0.5).ToManagedImage();//gabor滤波
            CvInvoke.Threshold(new Image<Gray, byte>(bm1), dst1, 150, 255, ThresholdType.Binary);
            //45度
            Bitmap bm2 = gf.ProcessFilter(UnmanagedImage.FromManagedImage(src.Bitmap), dst2.Bitmap, 3, 5, Math.PI / 4, 0.3, 2, 0.5).ToManagedImage();//gabor滤波
            CvInvoke.Threshold(new Image<Gray, byte>(bm2), dst2, 150, 255, ThresholdType.Binary);
            //135度
            Bitmap bm3 = gf.ProcessFilter(UnmanagedImage.FromManagedImage(src.Bitmap), dst3.Bitmap, 3, 5, 135 * Math.PI / 180, 0.3, 2, 0.5).ToManagedImage();//gabor滤波
            CvInvoke.Threshold(new Image<Gray, byte>(bm3), dst3, 150, 255, ThresholdType.Binary);

            src.Dispose();
            bm1.Dispose();
            bm2.Dispose();
            bm3.Dispose();
            pictureBox2.Image = (dst1 * 255).Bitmap;
            pictureBox3.Image = (dst2 * 255).Bitmap;
            pictureBox4.Image = (dst3 * 255).Bitmap;
        }


using Accord.Math;
using System;
using System.Collections.Generic;
using System.Drawing;
using System.Drawing.Imaging;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using VisioForge.Shared.AForge.Imaging;
using VisioForge.Shared.AForge.Math;


    public enum GaborKernelKind
    {

        /// <summary>
        /// 根据Gabor函数的实部创建内核。
        /// </summary>
        /// 
        Real,

        /// <summary>
        /// 根据Gabor函数的虚部创建内核。
        /// </summary>
        /// 
        Imaginary,

        /// <summary>
        /// 根据Gabor函数的大小创建内核。
        /// </summary>
        /// 
        Magnitude,

        /// <summary>
        /// 根据Gabor函数的平方大小创建内核。
        /// </summary>
        /// 
        SquaredMagnitude
    };
    public class GaborFilter
    {
        private Dictionary<PixelFormat, PixelFormat> formatTranslations
            = new Dictionary<PixelFormat, PixelFormat>();

        private double[,] kernel;

        /// <summary>
        /// 
        /// </summary>
        /// <param name="sourceData">输入图片</param>
        /// <param name="bitmap">输出储存图片</param>
        /// <param name="sizes">获取或设置过滤器的大小。 内核大小 默认值为3。</param>
        /// <param name="wavelength">获取或设置滤光片的波长。 默认值为4.0。</param>
        /// <param name="angle">获取或设置过滤器的方向，以弧度为单位。 默认值为0.6。</param>
        /// <param name="psis"> 获取或设置滤波器的相位偏移。 默认值为1.0。</param>
        /// <param name="sigmas">获取或设置过滤器的高斯方差。 默认值为2。</param>
        /// <param name="gammas">获取或设置滤镜的纵横比。 默认值为0.3。</param>
        /// <returns></returns>
        public unsafe UnmanagedImage ProcessFilter(UnmanagedImage sourceData, Bitmap bitmap,
            int sizes, double wavelength, double angle, double psis, double sigmas, double gammas)
        {
            UnmanagedImage destinationData = UnmanagedImage.FromManagedImage(bitmap);
            
            kernel = Kernel2D(sizes, wavelength, angle, psis, sigmas, gammas, true, GaborKernelKind.Imaginary);
            

            int kernelHeight = kernel.GetLength(0);
            int kernelWidth = kernel.GetLength(1);

            int centerX = kernelHeight / 2;
            int centerY = kernelWidth / 2;

            int width = sourceData.Width;
            int height = sourceData.Height;

            int srcStride = sourceData.Stride;
            int srcOffset = srcStride - width;

            byte* src = (byte*)sourceData.ImageData.ToPointer();


            int[,] response = new int[height, width];

            int max = int.MinValue;
            int min = int.MaxValue;

            // for each image row
            for (int y = 0; y < height; y++)
            {
                // for each pixel in the row
                for (int x = 0; x < width; x++, src++)
                {
                    double sum = 0;

                    // for each kernel row
                    for (int i = 0; i < kernelHeight; i++)
                    {
                        int ir = i - centerY;
                        int t = y + ir;

                        // skip row
                        if (t < 0)
                            continue;

                        // break
                        if (t >= height)
                            break;

                        int col = ir * srcStride;

                        // for each kernel value in the row
                        for (int j = 0; j < kernelWidth; j++)
                        {
                            int jr = j - centerX;
                            t = x + jr;

                            // skip column
                            if (t < 0)
                                continue;

                            if (t < width)
                            {
                                double k = kernel[i, j];
                                sum += k * src[col + jr];
                            }
                        }

                        int v = response[y, x] = (int)sum;

                        if (v > max) max = v;
                        if (v < min) min = v;
                    }
                }

                src += srcOffset;
            }


            byte* dst = (byte*)destinationData.ImageData.ToPointer();
            int pixelSize = System.Drawing.Image.GetPixelFormatSize(destinationData.PixelFormat) / 8;
            int dstStride = destinationData.Stride;
            int dstOffset = dstStride - width * pixelSize;

            if (destinationData.PixelFormat == PixelFormat.Format8bppIndexed)
            {
                // for each image row
                for (int y = 0; y < height; y++)
                {
                    // for each pixel in the row
                    for (int x = 0; x < width; x++, dst++)
                    {
                        *dst = (byte)((255 * (response[y, x] - min)) / (max - min));
                    }

                    dst += dstOffset;
                }
            }
            else
            {
                // for each image row
                for (int y = 0; y < height; y++)
                {
                    // for each pixel in the row
                    for (int x = 0; x < width; x++, dst += pixelSize)
                    {
                        int v = response[y, x];

                        if (v > 0)
                        {
                            dst[RGB.R] = (byte)((255 * v) / max);
                        }
                        else // (v <= 0)
                        {
                            dst[RGB.B] = (byte)((255 * v) / min);
                        }
                    }

                    dst += dstOffset;
                }
            }
            return destinationData;
        }

        /// <summary>
        ///   1-D Gabor function.
        /// </summary>
        /// 
        public static double Function1D(double x, double mean, double amplitude,
            double position, double width, double phase, double frequency)
        {
            double a = (x - position) * (x - position);
            double b = (2 * width) * (2 * width);

            double envelope = mean + amplitude * Math.Exp(-a / b);
            double carry = Math.Cos(2 * Math.PI * frequency * (x - position) + phase);

            return envelope * carry;
        }

        /// <summary>
        ///   2-D Gabor function.
        /// </summary>
        /// 
        public static Complex Function2D(int x, int y, double lambda, double theta,
            double psi, double sigma, double gamma)
        {
            double X = +x * Math.Cos(theta) + y * Math.Sin(theta);
            double Y = -x * Math.Sin(theta) + y * Math.Cos(theta);

            double envelope = Math.Exp(-((X * X + gamma * gamma * Y * Y) / (2 * sigma * sigma)));
            double real = Math.Cos(2 * Math.PI * (X / lambda) + psi);
            double imaginary = Math.Sin(2 * Math.PI * (X / lambda) + psi);

            return new Complex(envelope * real, envelope * imaginary);
        }

        /// <summary>
        ///   Real part of the 2-D Gabor function.
        /// </summary>
        /// 
        public static double RealFunction2D(int x, int y, double lambda, double theta,
            double psi, double sigma, double gamma)
        {
            double X = +x * Math.Cos(theta) + y * Math.Sin(theta);
            double Y = -x * Math.Sin(theta) + y * Math.Cos(theta);

            double envelope = Math.Exp(-((X * X + gamma * gamma * Y * Y) / (2 * sigma * sigma)));
            double carrier = Math.Cos(2 * Math.PI * (X / lambda) + psi);

            return envelope * carrier;
        }

        /// <summary>
        ///   Imaginary part of the 2-D Gabor function.
        /// </summary>
        /// 
        public static double ImaginaryFunction2D(int x, int y, double lambda, double theta,
            double psi, double sigma, double gamma)
        {
            double X = +x * Math.Cos(theta) + y * Math.Sin(theta);
            double Y = -x * Math.Sin(theta) + y * Math.Cos(theta);

            double envelope = Math.Exp(-((X * X + gamma * gamma * Y * Y) / (2 * sigma * sigma)));
            double carrier = Math.Sin(2 * Math.PI * (X / lambda) + psi);

            return envelope * carrier;
        }

        /// <summary>
        ///   Computes the 2-D Gabor kernel.
        /// </summary>
        /// 
        public static double[,] Kernel2D(double lambda, double theta, double psi,
            double sigma, double gamma)
        {
            return Kernel2D(3, lambda, theta, psi, sigma, gamma, false, GaborKernelKind.Real);
        }

        /// <summary>
        ///   Computes the 2-D Gabor kernel.
        /// </summary>
        /// 
        public static double[,] Kernel2D(double lambda, double theta, double psi,
            double sigma, double gamma, bool normalized)
        {
            return Kernel2D(3, lambda, theta, psi, sigma, gamma, normalized, GaborKernelKind.Imaginary);
        }

        /// <summary>
        ///   Computes the 2-D Gabor kernel.
        /// </summary>
        /// 
        public static double[,] Kernel2D(int size, double lambda, double theta, double psi,
            double sigma, double gamma, bool normalized)
        {
            return Kernel2D(size, lambda, theta, psi, sigma,
                gamma, normalized, GaborKernelKind.Imaginary);
        }


        /// <summary>
        ///   Computes the 2-D Gabor kernel.
        /// </summary>
        /// 
        public static double[,] Kernel2D(int size, double lambda, double theta,
            double psi, double sigma, double gamma, bool normalized, GaborKernelKind function)
        {
            double sigmaX = sigma;
            double sigmaY = sigma / gamma;

            double a = Math.Max(
                Math.Abs(size * sigmaX * Math.Cos(theta)),
                Math.Abs(size * sigmaY * Math.Sin(theta)));
            int xMax = (int)Math.Ceiling(Math.Max(1, a));

            double b = Math.Max(
                Math.Abs(size * sigmaX * Math.Sin(theta)),
                Math.Abs(size * sigmaY * Math.Cos(theta)));
            int yMax = (int)Math.Ceiling(Math.Max(1, b));

            int[] xValues = Vector.Interval(-xMax, xMax);
            int[] yValues = Vector.Interval(-yMax, yMax);

            Accord.Diagnostics.Debug.Assert(xValues.Length == (2 * xMax + 1));
            Accord.Diagnostics.Debug.Assert(yValues.Length == (2 * yMax + 1));

            double[,] kernel = new double[xValues.Length, yValues.Length];

            double sum = 0;

            switch (function)
            {
                case GaborKernelKind.Real:
                    for (int i = 0; i < xValues.Length; i++)
                        for (int j = 0; j < yValues.Length; j++)
                            sum += kernel[i, j] = RealFunction2D(
                                xValues[i], yValues[j], lambda, theta, psi, sigma, gamma);
                    break;

                case GaborKernelKind.Imaginary:
                    for (int i = 0; i < xValues.Length; i++)
                        for (int j = 0; j < yValues.Length; j++)
                            sum += kernel[i, j] = ImaginaryFunction2D(
                                    xValues[i], yValues[j], lambda, theta, psi, sigma, gamma);
                    break;

                case GaborKernelKind.Magnitude:
                    for (int i = 0; i < xValues.Length; i++)
                        for (int j = 0; j < yValues.Length; j++)
                            sum += kernel[i, j] = Function2D(
                                xValues[i], yValues[j], lambda, theta, psi, sigma, gamma).Magnitude;
                    break;

                case GaborKernelKind.SquaredMagnitude:
                    for (int i = 0; i < xValues.Length; i++)
                        for (int j = 0; j < yValues.Length; j++)
                        {
                            double v = Function2D(
                                xValues[i], yValues[j], lambda, theta, psi, sigma, gamma).Magnitude;
                            sum += kernel[i, j] = v * v;
                        }
                    break;
            }

            if (normalized)
                kernel.Divide(sum, result: kernel);

            return kernel;
        }
    }