OpenCVSharp 4.5 去除周期性噪声

用OpenCVSharp 4.5 跑一遍OpenCV官方教程

原官方教程链接：OpenCV: Periodic Noise Removing Filter

原作者有意埋了坑，完全照抄是不行的....... 

using System;
using OpenCvSharp;
namespace ConsoleApp1
{
    class tutorial33 : ITutorial
    {
        public void Run()
        {
            Mat imgIn = Cv2.ImRead("I:\\csharp\\images\\period_input.jpg", ImreadModes.Grayscale);
            if (imgIn.Empty()) //check whether the image is loaded or not
            {
                Console.Write("ERROR : Image cannot be loaded..!!");
                return;
            }
            imgIn.ConvertTo(imgIn, MatType.CV_32F);
            // it needs to process even image only
            Rect roi = new Rect(0, 0, imgIn.Cols & -2, imgIn.Rows & -2);
            imgIn = imgIn[roi];
            // PSD calculation (start)
            Mat imgPSD = new Mat();
            calcPSD(imgIn, out imgPSD);
            //fftshift(imgPSD, out imgPSD);
            Cv2.Normalize(imgPSD, imgPSD, 0, 255, NormTypes.MinMax);
            // PSD calculation (stop)

            //H calculation (start)
            Mat H = Mat.Ones(roi.Size, MatType.CV_32F);
            int r = 21;            
            synthesizeFilterH(ref H, new Point(143, 67), r);
            synthesizeFilterH(ref H, new Point(286, 0), r);
            synthesizeFilterH(ref H, new Point(0, 135), r);
            Mat imgHPlusPSD = imgPSD + H * 50;
            Cv2.Normalize(imgHPlusPSD, imgHPlusPSD, 0, 255, NormTypes.MinMax);
            imgHPlusPSD.ConvertTo(imgHPlusPSD, MatType.CV_8U);

            //H calculation (stop)
            // filtering (start)
            Mat imgOut = new Mat();
            filter2DFreq(imgIn, out imgOut, H);
            // filtering (stop)
            
            imgOut.ConvertTo(imgOut, MatType.CV_8U);
            Cv2.Normalize(imgOut, imgOut, 0, 255, NormTypes.MinMax);
            Cv2.ImShow("result.jpg", imgOut);
            Cv2.ImShow("filter", H);
            Cv2.ImShow("PSD", imgPSD);
            Cv2.WaitKey();
        }
        private void fftshift(Mat inputImg, out Mat outputImg)
        {
            Mat output = inputImg.Clone();
            int cx = output.Cols / 2;
            int cy = output.Rows / 2;
            Mat q0 = new Mat(output, new Rect(0, 0, cx, cy));
            Mat q1 = new Mat(output, new Rect(cx, 0, cx, cy));
            Mat q2 = new Mat(output, new Rect(0, cy, cx, cy));
            Mat q3 = new Mat(output, 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);
            outputImg = output;
        }

        private void filter2DFreq(Mat inputImg, out Mat outputImg, Mat H)
        {
            Mat m1 = inputImg.Clone();
            m1.ConvertTo(m1, MatType.CV_32F);

            Mat m2 = Mat.Zeros(inputImg.Size(), MatType.CV_32F);

            Mat[] planes = { m1, m2 };
            Mat complexI = new Mat();
            Cv2.Merge(planes, complexI);
            Cv2.Dft(complexI, complexI, DftFlags.Scale);
            
            H.ConvertTo(H, MatType.CV_32F);
            Mat[] planesH = { H, Mat.Zeros(H.Size(), MatType.CV_32F) };

            Mat complexH = new Mat();
            Cv2.Merge(planesH, complexH);

            Mat complexIH = new Mat();
            Cv2.MulSpectrums(complexI, complexH, complexIH, 0);
            Cv2.Idft(complexIH, complexIH);
            Cv2.Split(complexIH, out planes);
            outputImg = planes[0];
        }
        private void synthesizeFilterH(ref Mat inputOutput_H, Point center, int radius)
        {            
            Point c2 = center, c3 = center, c4 = center;
            c2.Y = inputOutput_H.Rows - center.Y;
            c3.X = inputOutput_H.Cols - center.X;
            c4 = new Point(c3.X, c2.Y);
            Cv2.Circle(inputOutput_H, (int)center.X, (int)center.Y, radius, new Scalar(0), Cv2.FILLED, LineTypes.Link8);
            Cv2.Circle(inputOutput_H, (int)c2.X, (int)c2.Y, radius, new Scalar(0), Cv2.FILLED, LineTypes.Link8);
            Cv2.Circle(inputOutput_H, (int)c3.X, (int)c3.Y, radius, new Scalar(0), Cv2.FILLED, LineTypes.Link8);
            Cv2.Circle(inputOutput_H, (int)c4.X, (int)c4.Y, radius, new Scalar(0), Cv2.FILLED, LineTypes.Link8);
        }
        // Function calculates PSD(Power spectrum density) by fft with two flags
        // flag = 0 means to return PSD
        // flag = 1 means to return log(PSD)
        private void calcPSD(Mat inputImg, out Mat outputImg, int flag=0)    
        {
            Mat m1 = inputImg.Clone();
            m1.ConvertTo(m1, MatType.CV_32F);

            Mat m2 = Mat.Zeros(inputImg.Size(), MatType.CV_32F);

            Mat[] planes = { m1, m2 };


            Mat complexI = new Mat();
            Cv2.Merge(planes, complexI);

            Cv2.Dft(complexI, complexI);
            Cv2.Split(complexI, out planes);
            // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
            planes[0].At<float>(0) = 0;
            planes[1].At<float>(0) = 0;

            // compute the PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
            Mat imgPSD = new Mat();
            Cv2.Magnitude(planes[0], planes[1], imgPSD);        //imgPSD = sqrt(Power spectrum density)
            Cv2.Pow(imgPSD, 2, imgPSD);                         //it needs ^2 in order to get PSD
            outputImg = imgPSD;
            // logPSD = log(1 + PSD)
            if (flag > 0)
            {
                Mat imglogPSD;
                imglogPSD = imgPSD + 1;
                Cv2.Log(imglogPSD, imglogPSD);
                outputImg = imglogPSD;
            }
        }
    }
}

 原图

局部放大

去除周期性噪声后结果