高斯模糊算法的全面优化过程分享。

     这里的高斯模糊采用的是论文《Recursive implementation of the Gaussian filter》里描述的递归算法。

              

      仔细观察和理解上述公式，在forward过程中,n是递增的，因此，如果在进行forward之前，把in数据先完整的赋值给w，然后式子（9a）就可以变为：

  　　　　w[n] = B w[n] + (b₁ w[n-1] + b₂ w[n-2] + b₃ w[n-3]) / b0；　　　　　     --------->     (1a)

   　在backward过程中，n是递减的，因此在backward前，把w的数据完整的拷贝到out中，则式子(9b)变为：

 　　　　out[n] = B out[n] + (b₁ out[n+1] + b₂ out[n+2] + b₃ out[n+3]) / b0 ；     <---------     (1b)

     从编程角度看来，backward中的拷贝是完全没有必要的，因此 (1b)可以直接写为：

           w[n] = B w[n] + (b₁ w[n+1] + b₂ w[n+2] + b₃ w[n+3]) / b0 ；               <---------     (1c)

     从速度上考虑，浮点除法很慢，因此，我们重新定义b1 = b1 / b0, b2 = b2 / b0, b3 = b3 / b0, 最终得到我们使用的递归公式：

           w[n] = B w[n] + b₁ w[n-1] + b₂ w[n-2] + b₃ w[n-3]；　　　　　     --------->     (a)

 　　     w[n] = B w[n] + b₁ w[n+1] + b₂ w[n+2] + b₃ w[n+3] ；             <---------      (b)

 　  上述公式是一维的高斯模糊计算方法，针对二维的图像，正确的做法就是先对每个图像行进行模糊处理得到中间结果，再对中间结果的每列进行模糊操作，最终得到二维的模糊结果，当然也可以使用先列后行这样的做法。

      另外注意点就是，边缘像素的处理，我们看到在公式(a)或者(b)中，w[n]的结果分别依赖于前三个或者后三个元素的值，而对于边缘位置的值，这些都是不在合理范围内的，通常的做法是镜像数据或者重复边缘数据，实践证明，由于是递归的算法，起始值的不同会将结果一直延续下去，因此，不同的方法对边缘部分的结果还是有一定的影响的，这里我们采用的简单的重复边缘像素值的方式。

      由于上面公式中的系数均为浮点类型，因此，计算一般也是对浮点进行的，也就是说一般需要先把图像数据转换为浮点，然后进行高斯模糊，在将结果转换为字节类型的图像，因此，上述过程可以分别用一下几个函数完成：

                CalcGaussCof　　　　　　　　　　 //　　计算高斯模糊中使用到的系数
　　　　　　ConvertBGR8U2BGRAF　　　　　　//　　将字节数据转换为浮点数据　
　　　　　　GaussBlurFromLeftToRight　　　　//　　水平方向的前向传播
　　　　　　GaussBlurFromRightToLeft　　　　//　　水平方向的反向传播
　　　　　　GaussBlurFromTopToBottom　　  //　　 垂直方向的前向传播
　　　　　　GaussBlurFromBottomToTop　　  //　　 垂直方向的反向传播
　　　　　　ConvertBGRAF2BGR8U　　　　    //　　 将结果转换为字节数据

 　　我们依次分析之。

       CalcGaussCof，这个很简单，就按照论文中提出的计算公式根据用户输入的参数来计算，当然结合下上面我们提到的除法变乘法的优化，注意，为了后续的一些标号的问题，我讲上述公式(a)和(b)中的系数B写为b0了。

void CalcGaussCof(float Radius, float &B0, float &B1, float &B2, float &B3)
{
    float Q, B;
    if (Radius >= 2.5)
        Q = (double)(0.98711 * Radius - 0.96330);                            //    对应论文公式11b
    else if ((Radius >= 0.5) && (Radius < 2.5))
        Q = (double)(3.97156 - 4.14554 * sqrt(1 - 0.26891 * Radius));
    else
        Q = (double)0.1147705018520355224609375;

    B = 1.57825 + 2.44413 * Q + 1.4281 * Q * Q + 0.422205 * Q * Q * Q;        //    对应论文公式8c
    B1 = 2.44413 * Q + 2.85619 * Q * Q + 1.26661 * Q * Q * Q;
    B2 = -1.4281 * Q * Q - 1.26661 * Q * Q * Q;
    B3 = 0.422205 * Q * Q * Q;

    B0 = 1.0 - (B1 + B2 + B3) / B;
    B1 = B1 / B;
    B2 = B2 / B;
    B3 = B3 / B;
}

　　由上述代码可见，B0+B1+B2+B3=1，即是归一化的系数，这也是算法可以递归进行的关键之一。

     接着为了方便中间过程，我们先将字节数据转换为浮点数据，这部分代码也很简单：

void ConvertBGR8U2BGRAF(unsigned char *Src, float *Dest, int Width, int Height, int Stride)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        float *LinePD = Dest + Y * Width * 3;
        for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3)
        {
            LinePD[0] = LinePS[0];    LinePD[1] = LinePS[1];    LinePD[2] = LinePS[2];
        }
    }
}

　　大家可以尝试自己把内部的X循环再展开看看效果。

     水平方向的前向传播按照公式去写也是很简单的，但是直接使用公式里面有很多访问数组的过程（不一定就慢），我稍微改造下写成如下形式：

void GaussBlurFromLeftToRight(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePD &