再談快速的高斯模糊演演算法（使用多次均值濾波逼近和擴充套件的二項式濾波濾波器）及其優化。

int IM_BinomialFilter4_PureC_Gray_Raw(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Sigma)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL) || (Dest == NULL))        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                            return IM_STATUS_INVALIDPARAMETER;
    if (Channel != 1)                                            return IM_STATUS_NOTSUPPORTED;
    int Status = IM_STATUS_OK;

    Sigma = IM_Min(Sigma, 256);            //    如果全部用整形的資料，半徑不能大於30， Sum和Deri用浮點半徑可以取得較大

    int Radius = (int)(sqrtf(3.0 * Sigma  * Sigma + 1) + 0.4999999f);
    int Weight = Radius * Radius * Radius * Radius;
    float InvWeight = 1.0f / Weight;
    unsigned char *Buffer = (unsigned char *)malloc((IM_Max(Width, Height) + 4 * Radius) * sizeof(unsigned char));
    unsigned char *Temp = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));
    int *RowOffset = (int *)malloc((Width + 4 * Radius) * sizeof(int));
    int *ColOffset = (int *)malloc((Height + 4 * Radius) * sizeof(int));

    if ((Buffer == NULL) || (Temp == NULL) || (RowOffset == NULL) || (ColOffset == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    //    注意起點位置是偏移了一個1個畫素的
    for (int X = 0; X < Width + 4 * Radius; X++)
        RowOffset[X] = IM_GetMirrorPos(Width, X - 2 * Radius + 1);

    for (int Y = 0; Y < Height + 4 * Radius; Y++)
        ColOffset[Y] = IM_GetMirrorPos(Height, Y - 2 * Radius + 1);

    for (int Y = 0; Y < Height; Y++)     // horizontal blur...
    {
        unsigned char *LinePS = Src + Y * Stride;
        unsigned char *LinePD = Temp + Y * Stride;
        for (int X = 0; X < 2 * Radius - 1; X++)
            Buffer[X] = LinePS[RowOffset[X]];
        memcpy(Buffer + 2 * Radius - 1, LinePS, Width);
        for (int X = 2 * Radius - 1 + Width; X < Width + 4 * Radius; X++)
            Buffer[X] = LinePS[RowOffset[X]];

        float Sum = 0, Deri1 = 0;
        int Deri2 = 0, Diff = 0;
        for (int X = -4 * Radius; X < 0; X++)
        {
            Sum += Deri1;
            Deri1 += Deri2;
            Deri2 += Diff;  // accumulate pixel blur
            Diff += Buffer[X + 4 * Radius];
            if (X + 3 * Radius >= 0) Diff -= 4 * Buffer[X + 3 * Radius];    //  -4,
            if (X + 2 * Radius >= 0) Diff += 6 * Buffer[X + 2 * Radius];    //      +6,
            if (X + Radius >= 0) Diff -= 4 * Buffer[X + 1 * Radius];        //          -4, (+1)}
        }
        for (int X = 0; X < Width; X++)
        {
            Sum += Deri1;
            Deri1 += Deri2;
            Deri2 += Diff;  
            Diff += Buffer[X] - 4 * Buffer[X + Radius] + 6 * Buffer[X + 2 * Radius] - 4 * Buffer[X + 3 * Radius] + Buffer[X + 4 * Radius];        //{+1,    -4,    +6,    -4,    +1}
            LinePD[X] = (int)(Sum * InvWeight + 0.499999f);            
            //LinePD[X] = (Sum + Weight / 2) / Weight;
        }
    }
    for (int X = 0; X < Width; X++)　 //　　Vertical Blur
    {
        unsigned char *LinePS = Temp + X;
        unsigned char *LinePD = Dest + X;
        for (int Y = 0; Y < 2 * Radius - 1; Y++)
            Buffer[Y] = LinePS[ColOffset[Y] * Stride];
        for (int Y = 0; Y < Height; Y++)
        {
            Buffer[Y + 2 * Radius - 1] = LinePS[Y * Stride];
        }
        for (int Y = 2 * Radius - 1 + Height; Y < Height + 4 * Radius; Y++)
            Buffer[Y] = LinePS[ColOffset[Y] * Stride];

        float Sum = 0, Deri1 = 0;
        int Deri2 = 0, Diff = 0;
        for (int Y = -4 * Radius; Y < 0; Y++)
        {
            Sum += Deri1;
            Deri1 += Deri2;
            Deri2 += Diff;  // accumulate pixel blur
            Diff += Buffer[Y + 4 * Radius];
            if (Y + 3 * Radius >= 0) Diff -= 4 * Buffer[Y + 3 * Radius];    //  -4,
            if (Y + 2 * Radius >= 0) Diff += 6 * Buffer[Y + 2 * Radius];    //      +6,
            if (Y + Radius >= 0) Diff -= 4 * Buffer[Y + 1 * Radius];        //          -4, (+1)}
        }
        for (int Y = 0; Y < Height; Y++)
        {
            Sum += Deri1;
            Deri1 += Deri2;
            Deri2 += Diff;
            Diff += Buffer[Y] - 4 * Buffer[Y + Radius] + 6 * Buffer[Y + 2 * Radius] - 4 * Buffer[Y + 3 * Radius] + Buffer[Y + 4 * Radius];        //{+1,    -4,    +6,    -4,    +1}
            LinePD[0] = (int)(Sum * InvWeight + 0.499999f);
            //LinePD[0] = (Sum + Weight / 2) / Weight;
            LinePD += Stride;
        }
    }
FreeMemory:

    if (Buffer != NULL)        free(Buffer);
    if (Temp != NULL)        free(Temp);
    if (RowOffset != NULL)    free(RowOffset);
    if (ColOffset != NULL)    free(ColOffset);
    return Status;
}

//    https://github.com/amazingyyc/fasted_gauss_blur
//    https://www.cnblogs.com/amazingyyc/p/4530634.html
//    https://www.doc88.com/p-1896883955048.html
//    使用 N 個 BoxBlur 擬合 Sigma的高斯函數
//    參考：Arbitrary Gaussian Filtering with 25 Addtions and 5 Multiplications per Pixel.pdf
//    Fast Almost - Gaussian Filtering.pdf
void IM_GetBoxRadiusForGauss(float Sigma, int *Size, int N)
{
    float IdealW = sqrt(12.0 * Sigma * Sigma / N + 1.0);            //    論文公式1
    int LowerW = floor(IdealW);                                        //    向下取整
    if (LowerW % 2 == 0)    LowerW--;                                //    the first odd integer that is less than IdealW
    int UpperW = LowerW + 2;                                        //    the other being being the first odd integer greater than IdealW
    int M = round((12.0 * Sigma * Sigma - N * LowerW * LowerW - 4 * N * LowerW - 3 * N) / (-4 * LowerW - 4));    //    論文公式2
    for (int Y = 0; Y < N; Y++)                                        //    Apply an averaging filter of width LowerW M times.
        Size[Y] = (Y < M ? LowerW : UpperW);                        //    Apply an averaging filter of width UpperW N - M times.
}

int IM_FastGaussBlur(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL) || (Dest == NULL))        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                            return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1) && (Channel != 3) && (Channel != 4))        return IM_STATUS_NOTSUPPORTED;
    int Status = IM_STATUS_OK;

    unsigned char *Temp = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));
    if (Temp == NULL)
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    int BoxRadius[3];                                //    使用3次均值模糊代替高斯
    IM_GetBoxRadiusForGauss(Radius, BoxRadius, 3);        //    計算擬合的半徑
    Status = IM_BoxBlur(Src, Dest, Width, Height, Stride, (BoxRadius[0] - 1) / 2);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    Status = IM_BoxBlur(Dest, Temp, Width, Height, Stride, (BoxRadius[1] - 1) / 2);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    Status = IM_BoxBlur(Temp, Dest, Width, Height, Stride, (BoxRadius[2] - 1) / 2);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
FreeMemory:
    if (Temp != NULL)    free(Temp);
    return Status;
}