opencvSharp利用傅里叶变换找出图片的纹理方向

工作中有一个需求，在一张有较为明显纹理方向的图片中，找出它的纹理方向，一开始试了挺多方法，比如轮廓检测或者gabor滤波器后，用得出来的白色像素的点拟合直线，效果还可以接受，但依然会有一点误差，特别是对比度较低的图片，误差更大，因为有时候轮廓点不太能反应出纹理方向，而且计算速度也比较慢。
然后在知网找到了这篇论文——***利用图像边缘信息估算图像纹理主方向***，了解到了傅里叶变换后的频谱图能准确反应出纹理方向，于是参考着这篇论文写了一个比之前那个更快更准的方法，步骤如下：

1. 导入图片，缩小成200*200，转灰度图。
2. 对该图进行傅里叶变换，得到其频谱图。
3. 定义四条斜率分别为0、45、90、135的直线。
4. 将频谱图二值化后，再遍历里面的所有像素，如果是黑色像素则跳过操作，如果是白色像素，则计算该像素与给定的四条直线的距离，看看这个距离是否在给定的阈值内，如果在，那就在这个方向上累加1，最后统计出这四个方向上哪个方向的累加结果最大。

先来看看用来做测试的图片，有0、45、90、135这四个纹理方向的图片：

然后看看频谱图与二值化之后的结果，以及角度累加统计结果：

可以明显看出规律了，纹理方向和频谱图的高频部分几乎就是旋转90°的关系，如果有两个纹理方向，也能大致检测出来。
另外由于中间的高频部分会对方向检测有一定的干扰，所以就将其去掉了。
那么接下来就是怎么把方向检测出来了，参考论文里面的用的是直角坐标系转极坐标系的方法，把高频与低频之间的方向梯度算出来，然后将其转到极坐标系中，就能统计出哪个方向出现的次数最多。
我也试了一下这个方法，但得出来的结果不正确，也不知道哪里出了问题，所以就另想了上面的方法，效果也很不错，如果想检测更多的角度，那就构造更多角度的直线。
最后是代码：

        public Form1()
        {
            InitializeComponent();
            var src = Cv2.ImRead("C:\\Users\\Administrator\\Desktop\\测试\\...");
            Mat dst = new Mat();
            Cv2.Resize(src, dst, new Size(200, 200));
            //如果想检测更多的角度，那就构造更多角度的直线
            Line2D line0 = new Line2D(1, 0, dst.Width / 2, dst.Height / 2);
            Line2D line45 = new Line2D(-1, 1, dst.Width / 2, dst.Height / 2);
            Line2D line90 = new Line2D(0.01, 1, dst.Width / 2, dst.Height / 2);
            Line2D line135 = new Line2D(1, 1, dst.Width / 2, dst.Height / 2);
            //line0.FitSize(dst.Width, dst.Height, out var p1, out var p2);
            //line45.FitSize(dst.Width, dst.Height, out var p3, out var p4);
            //line90.FitSize(dst.Width, dst.Height, out var p5, out var p6);
            //line135.FitSize(dst.Width, dst.Height, out var p7, out var p8);
            //dst.Line(p1, p2, Scalar.Red);
            //dst.Line(p3, p4, Scalar.Red);
            //dst.Line(p5, p6, Scalar.Red);
            //dst.Line(p7, p8, Scalar.Red);
            //Cv2.ImShow("src", dst);
            
            Cv2.CvtColor(dst, dst, ColorConversionCodes.BGR2GRAY);
            dst = DFT(dst);
            //Cv2.ImShow("dft", dst);
            //将傅里叶变换之后的频谱图二值化
            dst.ConvertTo(dst, MatType.CV_8UC1);
            Cv2.Threshold(dst, dst, 0, 255, ThresholdTypes.Binary | ThresholdTypes.Otsu);
            //由于中间的高频部分会对方向检测有一定的干扰，所以就将其去掉
            dst[new Rect(dst.Width / 2 - 10, dst.Height / 2 - 10, 20, 20)].SetTo(0);
            //Cv2.ImShow("dftBinary", dst);
            //论文中的极坐标系方法，也一起放出来了，用的时候不需要上面的二值化和下面的计算
            //CartToPolar(dst);

            //如果是白色像素，则计算该像素与给定的四条直线的距离，看看这个距离是否在给定的阈值内，
            //如果在，那就在这个方向上累加1
            int angle0 = 0;
            int angle45 = 0;
            int angle90 = 0;
            int angle135 = 0;
            for (int i = 0; i < dst.Width; i++)
            {
                for (int j = 0; j < dst.Height; j++)
                {
                    var color = dst.Get<Vec3b>(j, i);
                    if (color.Item0 > 240 && color.Item1 > 240 && color.Item2 > 240)
                    {
                        //这里的5可以根据实际情况而定
                        if (line0.Distance(i, j) < 5)
                            angle0++;
                        else if (line45.Distance(i, j) < 5)
                            angle45++;
                        else if (line90.Distance(i, j) < 5)
                            angle90++;
                        else if (line135.Distance(i, j) < 5)
                            angle135++;
                    }
                }
            }

            Console.WriteLine($"0°：{angle0}\r\n, 45°：{angle45}\r\n, 90°：{angle90}\r\n, 135°：{angle135}\r\n");
        }
        
        //这里的傅里叶变换方法就是照搬之前转载过的方法
        //https://blog.csdn.net/weixin_43950082/article/details/123946606?spm=1001.2014.3001.5501
        public static Mat DFT(Mat img)
        {
            Mat padded = new Mat();
            int m = Cv2.GetOptimalDFTSize(img.Rows);
            int n = Cv2.GetOptimalDFTSize(img.Cols); 
            Cv2.CopyMakeBorder(img, padded, 0, m - img.Rows, 0, n - img.Cols, BorderTypes.Constant, Scalar.All(0));

            Mat paddedF32 = new Mat();
            padded.ConvertTo(paddedF32, MatType.CV_32F);
            Mat[] planes = { paddedF32, Mat.Zeros(padded.Size(), MatType.CV_32F) };
            Mat complex = new Mat();
            Cv2.Merge(planes, complex);

            Mat dft = new Mat();
            Cv2.Dft(complex, dft);

            Mat[] dftPlanes;
            Cv2.Split(dft, out dftPlanes);

            Mat magnitude = new Mat();
            Cv2.Magnitude(dftPlanes[0], dftPlanes[1], magnitude);

            magnitude += Scalar.All(1);
            Cv2.Log(magnitude, magnitude);

            Mat spectrum = magnitude[
                new Rect(0, 0, magnitude.Cols & -2, magnitude.Rows & -2)];

            int cx = spectrum.Cols / 2;
            int cy = spectrum.Rows / 2;

            Mat q0 = new Mat(spectrum, new Rect(0, 0, cx, cy));   
            Mat q1 = new Mat(spectrum, new Rect(cx, 0, cx, cy));  
            Mat q2 = new Mat(spectrum, new Rect(0, cy, cx, cy));  
            Mat q3 = new Mat(spectrum, new Rect(cx, cy, cx, cy)); 

            Mat tmp = new Mat();
            q0.CopyTo(tmp);
            q3.CopyTo(q0);
            tmp.CopyTo(q3);

            q1.CopyTo(tmp);
            q2.CopyTo(q1);
            tmp.CopyTo(q2);

            Cv2.Normalize(spectrum, spectrum, 0, 1, NormTypes.MinMax);
            return spectrum;
        }
        
        public static void CartToPolar(Mat src)
        {
            Mat dx = new Mat(), dy = new Mat();
            Mat magnitude = new Mat(), angle = new Mat();
            Cv2.Sobel(src, dx, MatType.CV_32F, 1, 0);
            Cv2.Sobel(src, dy, MatType.CV_32F, 0, 1);
            Cv2.Magnitude(dx, dy, magnitude);
            Cv2.CartToPolar(dx, dy, magnitude, angle, true);

            int angle0 = 0;
            int angle45 = 0;
            int angle90 = 0;
            int angle135 = 0;
            int angle180 = 0;
            int angle225 = 0;
            int angle270 = 0;
            int angle315 = 0;
            for (int i = 0; i < angle.Width; i++)
            {
                for (int j = 0; j < angle.Height; j++)
                {
                    var a = angle.Get<float>(j, i);
                    if (a > 355 || a < 5)
                        angle0++;
                    else if (a > 40 && a < 50)
                        angle45++;
                    else if (a > 85 && a < 95)
                        angle90++;
                    else if (a > 130 && a < 140)
                        angle135++;
                    else if (a > 175 && a < 185)
                        angle180++;
                    else if (a > 220 && a < 230)
                        angle225++;
                    else if (a >265 && a < 275)
                        angle270++;
                    else if (a > 310 && a < 320)
                        angle315++;
                }
            }
            Console.WriteLine($"0°：{angle0+ angle180}\r\n, 45°：{angle45 +  angle225}\r\n, 90°：{angle90 + angle270}\r\n, 135°：{angle135 + angle315}\r\n");
        }