OpenCVSharp 4.5 去除运动模糊

用OpenCVSharp 4.5 跑一遍 OpenCV 官方教程

原官方教程链接：OpenCV: Motion Deblur Filter

注：下面的程序运行是没问题的，遗憾的是我并没有找到一组合适的参数来重现原教程的效果，原OpenCV (c++)程序也是如此。（这种问题似乎很常见）

using OpenCvSharp;
using System;


namespace ConsoleApp1
{
    class tutorial31 : ITutorial
    {
        public void Run()
        {
            int len = 6  ; 
            double theta = 15; 
            double snr = 700;
            
           Mat imgIn = Cv2.ImRead("I:\\csharp\\images\\motion_blur.jpg", ImreadModes.Grayscale);
            if (imgIn.Empty()) //check whether the image is loaded or not
            {
                Console.WriteLine("ERROR : Image cannot be loaded..!!");
                return;
            }
            Mat imgOut = new Mat();
            // it needs to process even image only
            Rect roi = new Rect(0, 0, imgIn.Cols & -2, imgIn.Rows & -2);
            
            //Hw calculation (start)
            Mat Hw = new Mat();
            Mat h = new Mat();
            calcPSF(out h, roi.Size, len,theta);
            calcWnrFilter(h, out Hw, 1.0 / snr);
            //Hw calculation (stop)

            imgIn.ConvertTo(imgIn, MatType.CV_32F);
            edgetaper(imgIn, out imgIn);

            // filtering (start)
            filter2DFreq( imgIn[roi], out imgOut, Hw);
            // filtering (stop)

            imgOut.ConvertTo(imgOut, MatType.CV_8U);
            Cv2.Normalize(imgOut, imgOut, 0, 255, NormTypes.MinMax);
            Cv2.ImShow("result.jpg", imgOut);
            Cv2.WaitKey();
        }

        private void calcPSF(out Mat outputImg, Size filterSize, int len, double theta)
        {
            Mat h = Mat.Zeros(filterSize, MatType.CV_32F);
            Point point = new Point(filterSize.Width / 2, filterSize.Height / 2);
            Cv2.Ellipse(h, point, new Size(0, Math.Round( (float)len/2.0)), 90.0 - theta, 0, 360, new Scalar(255), Cv2.FILLED);
                        
            Scalar summa = Cv2.Sum(h);
            outputImg = h / summa[0];
        }

        private void filter2DFreq(Mat inputImg, out Mat outputImg, Mat H)
        {
            inputImg.Clone().ConvertTo(inputImg,MatType.CV_32F);
            
            Mat[] planes = { inputImg, Mat.Zeros(inputImg.Size(), MatType.CV_32F) };

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

            H.Clone().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 fftshift(Mat inputImg, out Mat outputImg)
        {
            outputImg = inputImg.Clone();
            int cx = outputImg.Cols / 2;
            int cy = outputImg.Rows / 2;
            Mat q0 = new Mat(outputImg, new Rect(0, 0, cx, cy));
            Mat q1 = new Mat(outputImg, new Rect(cx, 0, cx, cy));
            Mat q2 = new Mat(outputImg, new Rect(0, cy, cx, cy));
            Mat q3 = new Mat(outputImg, 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);
        }

        private void calcWnrFilter(Mat input_h_PSF, out Mat output_G, double nsr)
        {
            Mat h_PSF_shifted = new Mat();
            fftshift(input_h_PSF, out h_PSF_shifted);

            Mat[] planes = { h_PSF_shifted.Clone(), Mat.Zeros(h_PSF_shifted.Size(), MatType.CV_32F) };
            Mat complexI = new Mat();
            Cv2.Merge(planes, complexI);
            Cv2.Dft(complexI, complexI);
            Cv2.Split(complexI, out planes);
            Mat denom = new Mat();
            Cv2.Pow(Cv2.Abs(planes[0]), 2, denom);
            denom += nsr;
            Mat output = new Mat();
            Cv2.Divide(planes[0], denom, output);
            output_G = output;
        }

        private void edgetaper(Mat inputImg, out Mat outputImg, double gamma=5.0, double beta=0.2)
        {
            int Nx = inputImg.Cols;
            int Ny = inputImg.Rows;

            Mat w = new Mat();
            using (Mat w1 = Mat.Zeros(new Size(Nx,1), MatType.CV_32F))
            using (Mat w2 = Mat.Zeros(new Size(1,Ny), MatType.CV_32F))
            {
                float dx = (float)(2.0 * Cv2.PI / Nx);
                float x = (float)(-Cv2.PI);
                //var p1 = w1.GetGenericIndexer<float>();

                for (int i = 0; i < Nx; i++)
                {
                    float f = (float)(0.5 * (Math.Tanh((x + gamma / 2) / beta) - Math.Tanh((x - gamma / 2) / beta)));
                    w1.Set<float>(0,i,f);
                    x += dx;
                }
                                

                float dy = (float)(2.0 * Cv2.PI / Ny);
                float y = (float)(-Cv2.PI);
                //var p2 = w2.GetGenericIndexer<float>();

                for (int j = 0; j < Ny; j++)
                {
                    float f = (float)(0.5 * (Math.Tanh((y + gamma / 2) / beta) - Math.Tanh((y - gamma / 2) / beta)));
                    w2.Set<float>(j,0,f);
                    y += dy;
                }
                w = w2 * w1; //w1:431x1, w2:1x240 = w:431x240
            }
            Mat output = new Mat();
            Cv2.Multiply(inputImg, w, output);  //inputImg: 240x431
            outputImg = output;
        }
    }
}