【短道速滑十】非局部均值滤波的指令集优化和加速（针对5*5的搜索特例，可达到单核1080P灰度图 28ms/帧的速度）。

//    SearchRadius    搜索的领域半径
//    PatchRadius      计算相似度的块的半径
//    Delta     高斯平滑的参数
int IM_NLM_Denoising(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int SearchRadius, int PatchRadius, float Delta)
{
    int Status = IM_STATUS_OK;
    int Channel = Stride / Width;
    if (Src == NULL)                        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))        return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1))                        return IM_STATUS_INVALIDPARAMETER;    

    int ExpandW = Width + SearchRadius + PatchRadius + SearchRadius + PatchRadius;
    int ExpandH = Height + SearchRadius + PatchRadius + SearchRadius + PatchRadius;
    int DiffW = Width + PatchRadius + PatchRadius;
    int DiffH = Height + PatchRadius + PatchRadius;

    float *Weight = (float *)calloc(Width * Height , sizeof(float));
    float *Sum = (float *)calloc(Width * Height , sizeof(float));
    int *ColValue = (int *)malloc((Width + PatchRadius + PatchRadius) * sizeof(int));
    unsigned char *Expand = (unsigned char *)malloc(ExpandH * ExpandW * sizeof(unsigned char));
    if ((Weight == NULL) || (Sum == NULL) || (ColValue == NULL) || (Expand == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    Status = IM_GetExpandImage(Src, Expand, Width, Height, Stride, ExpandW, SearchRadius + PatchRadius, SearchRadius + PatchRadius, SearchRadius + PatchRadius, SearchRadius + PatchRadius, IM_EDGE_MIRROR);
    if (Status != IM_STATUS_OK)    goto FreeMemory;

    int Area = (2 * PatchRadius + 1) * (2 * PatchRadius + 1);
    float Inv = -1.0f / Area / Delta / Delta;


    for (int YY = -SearchRadius; YY <= SearchRadius; YY++)
    {
        for (int XX = -SearchRadius; XX <= SearchRadius; XX++)
        {
            for (int Y = 0; Y < Height; Y++)
            {
                float *LinePS = Sum + Y * Width;
                float *LinePW = Weight + Y * Width;
                unsigned char *LinePE = Expand + (Y + YY + SearchRadius + PatchRadius) * ExpandW + XX + SearchRadius + PatchRadius;
                if (Y == 0)
                {
                    memset(ColValue, 0, DiffW * sizeof(int));
                    for (int Z = -PatchRadius; Z <= PatchRadius; Z++)
                    {
                        unsigned char *LineP1 = Expand + (Z + PatchRadius + YY + SearchRadius) * ExpandW + XX + SearchRadius;
                        unsigned char *LineP2 = Expand + (Z + PatchRadius + SearchRadius) * ExpandW + SearchRadius;
                        for (int X = 0; X < DiffW; X++)
                        {
                            int Value = LineP2[X] - LineP1[X];
                            ColValue[X] += Value * Value;
                        }
                    }
                }
                else
                {
                    unsigned char *LineOut1 = Expand + (Y - 1 + YY + SearchRadius) * ExpandW + XX + SearchRadius;
                    unsigned char *LineOut2 = Expand + (Y - 1 + SearchRadius) * ExpandW + SearchRadius;
                    unsigned char *LineIn1 = Expand + (Y + PatchRadius + PatchRadius + YY + SearchRadius) * ExpandW + XX + SearchRadius;
                    unsigned char *LineIn2 = Expand + (Y + PatchRadius + PatchRadius + SearchRadius) * ExpandW + SearchRadius;
                    for (int X = 0; X < DiffW; X++)
                    {
                        int Out = LineOut2[X] - LineOut1[X];
                        int In = LineIn2[X] - LineIn1[X];
                        ColValue[X] -= Out * Out - In * In;                                    //    更新列数据
                    }
                }
                int SumA = IM_SumofArray(ColValue, PatchRadius * 2 + 1);                //    处理每行第一个数据    
                float W = IM_Fexp(SumA * Inv);
                LinePW[0] += W;
                LinePS[0] += W * LinePE[0];
                for (int X = 1; X < Width; X++)
                {
                    SumA = SumA - ColValue[X - 1] + ColValue[X + PatchRadius + PatchRadius];
                    float W = IM_Fexp(SumA * Inv);
                    LinePW[X] += W;
                    LinePS[X] += W * LinePE[X];
                }
            }
        }
    }
    for (int Y = 0; Y < Height; Y++)
    {
        int Index = Y * Width;
        unsigned char *LinePD = Dest + Y * Stride;
        for (int X = 0; X < Width; X++)
        {
            LinePD[X] = Sum[Index + X] / (Weight[Index + X]);
        }
    }
FreeMemory:
    if (Expand != NULL)            free(Expand);
    if (Weight != NULL)            free(Weight);
    if (Sum != NULL)               free(Sum);
    if (ColValue != NULL)         free(ColValue);
    return Status;
}