高斯金字塔的秘密（二，加速高斯，c#实现）

我们发现在进行高斯金字塔计算（尺度变换）时，速度很慢，现在工业相机用2千万的很多，怎么提高他的计算速度？我们的图像都是二维的，所以使用高斯平滑（滤波或称作尺度变换）也都是二维的。

1，生成高斯模板5*5的卷积肯定没有3*3的快，也就是说模板越大越慢

2，因为二维高斯函数具有可拆分性，可以拆分成两个一维高斯函数代替，比如3*3的拆分如下：

3，第二种方式其实效率很高了，你可以试一试。这两种方法我们前面都试过了，有兴趣的可以前面找找看，下面这种方法比较满意，也是参考网上成果，移植c#成功，做个记录，后面用得着。

他们说叫IIR高斯递归法。

先看结果对比图，1024*768的图像，sigma=5（图是从结果图中截取出来的细节对比）：

左边是二维高斯分离法，耗时（0.249秒），右边的是快速高斯法，耗时（0.0640秒），差别是4倍，因为是c#的版本，迂回了一些，如果用c指针的话，效率还会提高。

具体参考：(19条消息) IIR递归高斯滤波_Etrue的博客-CSDN博客_递归高斯滤波

公式如下：

Forward：

　　Backward:

上面公式缺B值= 1.0 - ((b[1] + b[2] + b[3]) / b[0]);

 当初看原理时，有点懵逼，熟悉又陌生，为什么写成这个样子推导，晚上睡觉，梦见泰勒，突然明白了，这个公式在相机镜头径向畸变中也用到了，好了，躺下继续睡觉。

以下是本人翻译并编译成功的c#代码（还参考了b站的视频）：

第一，先生成上边的q，b0,b1,b2,b3，B.

 public struct gauss3_coefs
        {
            public double[] b //= new double[4]
                ;
            public double B //= 0.0
                ;
            public float sigma// = 0
                ;
            public int N //= 0
                ;
        };
        /*  
        * Calculate the coefficients for the recursive filter algorithm  
        * Fast Computation of gaussian blurring.  
        */
        public void compute_coefs3(ref gauss3_coefs c, float sigma)
        {
            /*  
             * Papers:  "Recursive Implementation of the gaussian filter.",  
             *          Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.  
             * formula: 11b       computation of q  
             *          8c        computation of b0..b1  
             *          10        alpha is normalization constant B  
             */
            double q, q2, q3;

            q = 0;

            if (sigma >= 2.5)
            {
                q = 0.98711 * sigma - 0.96330;
            }
            else if ((sigma >= 0.5) && (sigma < 2.5))
            {
                q = 3.97156 - 4.14554 * (float)Math.Sqrt((double)1 - 0.26891 * sigma);
            }
            else
            {
                q = 0.1147705018520355224609375;
            }

            q2 = q * q;
            q3 = q * q2;
            c.b[0] = (1.57825 + (2.44413 * q) + (1.4281 * q2) + (0.422205 * q3));
            c.b[1] = ((2.44413 * q) + (2.85619 * q2) + (1.26661 * q3));
            c.b[2] = (-((1.4281 * q2) + (1.26661 * q3)));
            c.b[3] = ((0.422205 * q3));
            c.B = 1.0 - ((c.b[1] + c.b[2] + c.b[3]) / c.b[0]);
            c.sigma = sigma;
            c.N = 3;        
        }

第二，获取要处理的buffer图像：

  int hh = 768;
            int ww = 1024;
            float[] 原图copy = new float[1024 * 768];
            MeGaugingLibNew.ImageSourceForm secondhello = meGaugingLib1.getImage();
            for (int j = 0; j < hh; j++)
                for (int i = 0; i < ww; i++)
                {
                    原图copy[j * ww + i] = secondhello.IS_orgImg[j * ww + i];
                }
           
            float[] glob_buffer1024768smooth1dot25 = new float[1024 * 768];
            byte[] outputgaos = new byte[1024 * 768];

            gauss3_coefs hello = new gauss3_coefs();
            hello.b = new double[4];
            float sigma = 1.25f * 4;
            compute_coefs3(ref hello, sigma);

第三，处理行数据：

     for (int y = 0; y < hh; y++)
            {
                float[] temp = new float[ww];
                float[] temp1 = new float[ww];
                for (int x = 0; x < ww; x++)//在这个地方c的指针是有优势的，c#要迂回一下。202210200907
                {
                    temp[x] = 原图copy[y * ww + x];                
                }           
                gausssmooth(temp,ref temp1, ww, 1, hello);
                for (int x = 0; x < ww; x++)
                {
                   原图copy[y * ww + x] = temp1[x];
                }
            }

第四，处理列数据：

  // 原图copy在这里已改变，存储了行一维高斯处理结果，即中间结果
            for (int x = 0; x < ww; x++)
            {
                float[] temp = new float[hh];
                float[] temp1 = new float[hh];
                for (int y= 0;y < hh; y++)
                {
                    temp[y] = 原图copy[y * ww + x];
                }               
                gausssmooth(temp, ref temp1, hh, 1, hello);//可以参考fft,我已经处理为1了，忘了
                for (int y = 0; y < hh; y++)
                {
                     原图copy[y * ww + x] = temp1[y];
                }
            }

第五，结束，这里的gausssmooth函数如下：

   public void gausssmooth(float[] input,ref  float[] output, int size, int rowstride, gauss3_coefs c)
        {
            /*  
             * Papers:  "Recursive Implementation of the gaussian filter.",  
             *          Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.  
             * formula: 9a        forward filter  
             *          9b        backward filter  
             *          fig7      algorithm  
             */
            int i = 0, n = 0, bufsize = 0;
            float[] w1;
            float[] w2;

            /* forward pass */
            bufsize = size + 3;
            size -= 1;
            w1 = new float[bufsize];
            w2 = new float[bufsize];
            w1[0] = input[0];
            w1[1] = input[0];
            w1[2] = input[0];
            for (i = 0, n = 3; i <= size; i++, n++)
            {
                w1[n] = (float)(c.B * input[i * rowstride] +
                                 ((c.b[1] * w1[n - 1] +
                                   c.b[2] * w1[n - 2] +
                                   c.b[3] * w1[n - 3]) / c.b[0]));
            }

            /* backward pass */
            w2[size + 1] = w1[size + 3];
            w2[size + 2] = w1[size + 3];
            w2[size + 3] = w1[size + 3];
            for (i = size, n = i; i >= 0; i--, n--)
            {
                w2[n] = output[i * rowstride] = (float)(c.B * w1[n] +
                                                     ((c.b[1] * w2[n + 1] +
                                                       c.b[2] * w2[n + 2] +
                                                       c.b[3] * w2[n + 3]) / c.b[0]));
            }

      
        }

所谓，君子性非异也，善假于物也。哈哈！做一回君子！