C#实现Matlab中的fft2函数和fftshift函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
 Mat fft2( Mat src)
        {
            Mat Fourier = new Mat();
            int mat_type = (int)src.Depth;
            if (mat_type>15)
            {
                Debug.Print("报错");
            }
            if (mat_type < 7)
            {
                VectorOfMat planes =new VectorOfMat( new Mat[]{ src, new Mat(src.Size, DepthType.Cv64F,1) });
                CvInvoke.Merge(planes,Fourier);
                CvInvoke.Dft(Fourier, Fourier, DxtType.Forward, Fourier.Rows);
            }
            else // 7<mat_type<15
            {
                Mat tmp=new Mat();
                CvInvoke.Dft(src, tmp, DxtType.Forward, Fourier.Rows);
                VectorOfMat planes=new VectorOfMat();
                CvInvoke.Split(tmp, planes);
                Image<Gray, Single> Ita = new Image<Gray, Single>(planes[0].Size);
                Image<Gray, Single> ItC = new Image<Gray, Single>(planes[0].Size);
                CvInvoke.CartToPolar(planes[0], planes[1], ItC, Ita,true); //将复数转化为幅值
                Fourier = planes[0];
            }
            return Fourier;
        }

 Matrix<double> fftshift(Mat src)
        {
            Size sz = src.Size;
            Point pt=new Point(0, 0);
            pt.X = (int)Math.Floor(sz.Width / 2.0);
            pt.Y = (int)Math.Floor(sz.Height / 2.0);
            Mat mat= circshift(src, pt);
            Matrix<double> xx = new Matrix<double>(mat.Size);
            mat.CopyTo(xx);
            mat.Dispose();
            return xx;

        }
        Mat circshift(Mat src, Point delta)
        {
            Size sz = src.Size;

            
            // 此种情况不需要移动
            if ((sz.Height == 1 && sz.Width == 1) || (delta.X == 0 && delta.Y == 0))
                return src;

            // 需要移动参数的变换，这样就能输入各种整数了
            int x = delta.X;
            int y = delta.Y;
            if (x > 0) x = x % sz.Width;
            if (y > 0) y = y % sz.Height;
            if (x < 0) x = x % sz.Width + sz.Width;
            if (y < 0) y = y % sz.Height + sz.Height;


            // 多维的Mat也能移动
            VectorOfMat  planes=new VectorOfMat();
            CvInvoke.Split(src, planes);

            for (int i = 0; i < planes.Size; i++)
            {
                // 竖直方向移动
                Mat tmp0=new Mat(), tmp1 = new Mat(), tmp2 = new Mat(), tmp3 = new Mat();
                Mat q0 = new Mat(planes[i], new Rectangle(0, 0, sz.Width, sz.Height - y));
                Mat q1 = new Mat(planes[i], new Rectangle(0, sz.Height - y, sz.Width, y));
                q0.CopyTo(tmp0);
                q1.CopyTo(tmp1);
                tmp0.CopyTo(new Mat( planes[i],new Rectangle(0, y, sz.Width, sz.Height - y)));
                tmp1.CopyTo(new Mat(planes[i],new Rectangle(0, 0, sz.Width, y)));

                // 水平方向移动
                Mat q2 = new Mat(planes[i], new Rectangle(0, 0, sz.Width - x, sz.Height));
                Mat q3 = new Mat(planes[i], new Rectangle(sz.Width - x, 0, x, sz.Height));
                q2.CopyTo(tmp2);
                q3.CopyTo(tmp3);
                tmp2.CopyTo(new Mat( planes[i],new Rectangle(x, 0, sz.Width - x, sz.Height)));
                tmp3.CopyTo(new Mat( planes[i],new Rectangle(0, 0, x, sz.Height)));
            }

            CvInvoke.Merge(planes, src);
            return src;
        }