图像处理之双边滤波介绍与源码达成

方今在看图像风格化的舆论的时候,频仍蒙受 Bilateral Filter。google
一波后,发现并不是如何不可了的东西,但它的考虑却很有借鉴意义。

1 双边滤波简介

  双边滤波(Bilateral filter)是一种非线性的滤波方法,是整合图像的空中邻近度和像素值相似度的一种折衷处理,同时考虑空域新闻和灰度相似性,达到保边去噪的目标。具有简易、非迭代、局地的特点。

  双边滤波器的补益是足以做边缘保存(edge preserving),一般过去用的维纳滤波或者高斯滤波去降噪,都会较分明地歪曲边缘,对于频繁细节的护卫成效并不显眼。双边滤波器顾名思义比高斯滤波多了一个高斯方差sigma-d,它是根据空间分布的高斯滤波函数,所以在边缘附近,离的较远的像素不会太多影响到边缘上的像素值,这样就有限补助了边缘附近像素值的保留。不过出于保存了过多的往往音讯,对于彩色图像里的高频噪声,双边滤波器不可见彻底的滤掉,只可以够对于低频新闻实行较好的滤波。

简介

Bilateral
Filter,普通话又称「双边滤波器」。相比较过去那多少个单纯使用地方消息进行滤波的
filter,Bilateral Filter 还考虑了颜色消息,能够确认保障边缘部分不会被过滤。

简易来说,一般的 filter 都以基于这样的公式实行滤波的:
\[
h(x)=k_{d}^{-1}{(x)}\iint_\infty^\infty{f(\zeta)c(\zeta, x)}
d\zeta \]

其中,\(k_{d}^{-1}{(x)}\)
是权重之和,\(f(\zeta)\)
能够领略为单个像素,\(c(\zeta, x)\)
能够明白为地方权重。

翻译成程序员能够通晓的言语,差不多是如此:

for (int i = -r; i <= r; i++) {
  for (int j = -r; j <= +r; j++) {
    newpixel += pixel[row+i][col+j] * c[i][j];
    k += c[i][j];
  }
}
pixel[row][col] = newPixel / k;

高斯函数也属于那类 filter。

但那种 filter
有3个欠缺:各向同性(不明白那个领悟对不对)。用那种滤波器,各个点受邻居的熏陶是平等的,即便它跟邻居像素恐怕差得相比较多,也会被邻居「同化」(举个例子:边缘被「和谐」掉了)。由此,有人提议了
Bilateral Filter。

Bilateral Filter 选择那样的公式:
\[
h(x)=k_{d}^{-1}{(x)}\iint_\infty^\infty{f(\zeta)c(\zeta,
x)s(f(\zeta), f(x))} d\zeta \]
相比之下之前的架势,最大的变动唯有是权值中加进了多少个 \(s(f(\zeta),
f(x))\),这几个东西也是权值,不过它不是采取地点消息,而是颜色消息
\(f(\zeta)\)。不管是哪个种类新闻,时势上来看都以相同的,但由于增添了颜色权值,却使滤波的结果有了分明例外,前面会提交效果图。

重新翻译成程序语言:

for (int i = -r; i <= r; i++) {
  for (int j = -r; j <= +r; j++) {
    newpixel += pixel[row+i][col+j] * c[i][j] * s(pixel[row][col], pixel[row+i][col+j]);
    k += c[i][j]*s(pixel[row][col], pixel[row+i][col+j]);
  }
}
pixel[row][col] = newPixel / k;

s
函数能够借鉴地点权值的思路。例如,能够行使那种办法定义(当然那几个是小编自身组织的):

function s(p1, p2) {
  return (255-abs(p1-p2)) / 255
}

如此,差的越来越多的颜料,所占权值越小。

若果要追求科学严苛一点,也不妨仿照高斯核函数的概念:
\[ c(\zeta-x) = e^{-{1\over2}({
{\zeta-x} \over {\sigma} } )^2} \\\\\\ s(\zeta-x) =
e^{-{1\over2}({ {f(\zeta)-f(x)} \over \sigma })^2} \]

2 双边滤波原理

  滤波算法中,指标点上的像素值常常是由其所在地点上的四周的三个小片段邻居像素的值所决定。在2D高斯滤波中的具体落到实处便是对周围的终将范围内的像素值分别赋以分化的高斯权重值,并在加权平均后拿走当前点的末梢结出。而那边的高斯权重因子是应用多个像素之间的空中远距离(在图像中为2D)关系来扭转。通过高斯分布的曲线能够发现,离目的像素越近的点对最终结出的进献越大,反之则越小。其公式化的讲述相似如下所述:

                     图片 1

   个中的c即为根据空间距离的高斯权重,而用来对结果开始展览单位化。

  高斯滤波在低通滤波算法中有不利的展现,可是其却有其余1个标题,那正是只考虑了像素间的空中地方上的关系,由此滤波的结果会丢掉边缘的音信。那里的边缘主假使指图像中关键的不比颜色区域(比如褐绿的苍天,肉色的毛发等),而Bilateral就是在Gaussian blur中进入了别的的多个权重分部来消除这一标题。Bilateral滤波中对于边缘的维持通过下述表明式来兑现:

                     图片 2

                    图片 3

  当中的s为基于像素间相似程度的高斯权重,同样用来对结果进行单位化。对两端进行理并了结合即能够获取基于空间距离、相似程度综合考量的Bilateral滤波:

                  图片 4

  上式中的单位化分部综合了三种高斯权重于一块而取得,在这之中的c与s总括可以详细描述如下:

                  图片 5

   且有图片 6

                   图片 7

   且有图片 8

  上述给出的表明式均是在半空上的最为积分,而在像素化的图像中本来不可能这么做,而且也没供给如此做,由此在行使前必要对其进展离散化。而且也不须要对此每一种局地像素从整张图像上海展览中心开加权操作,距离超越一定水平的像素实际上对近年来的对象像素影响非常小,能够忽略的。限定局地子区域后的离散化公就能够简化为如下方式:

                   图片 9

  上述申辩公式就重组了比拉teral滤波达成的根底。为了直观地问询高斯滤波与两岸滤波的区分,大家能够从下列图示中阅览依照。要是目的源图像为下述左右区域泾渭显明的涵盖噪声的图像(由程序自动生成),蓝灰框的基本即为目的像素所在的职位,那么当前像素处所对应的高斯权重与多头权重因子3D可视化后的形态如后面两图所示:

                   图片 10

  左图为原本的噪音图像;中间为高斯采集样品的权重;右图为比拉teral采集样品的权重。从图中得以看出Bilateral加入了貌似程度分部以往能够将源图像右边那么些跟当前像素差值过大的点给滤去,那样就很好地保险了边缘。为了越发形象地观望两者间的界别,使用Matlab将该图在二种不一致格局下的惊人图3D绘制出来,如下:

                 图片 11

  上述三图从左到右依次为:双边滤波,原始图像,高斯滤波。从高度图中能够明显看出Bilateral和Gaussian二种格局的界别,前者较好地涵养了边缘处的梯度,而在高斯滤波中,由于其在边缘处的变型是线性的,因此就动用连累的梯度显示出渐变的动静,而那突显在图像中的话就是境界的散失(图像的示范可知于后述)。         

代码完结

知情原理后,达成其实也很粗略,上边给出的伪代码基本是骨干算法了。别的索要留意的是,借使是彩色图的话,要求对种种通道的颜色值进行滤波。

现实落到实处能够参考这篇博客:图像处理之双边滤波效果(Bilateral Filtering
for 格雷 and Color
Image)
,或然参考我自个儿的
demo,当然,小编也只是将上边博客的
java 版改成 c++ 而已^0^。

交付几幅结果图:

原图

图片 12

高斯模糊

图片 13

偏偏用颜色音信滤波

图片 14

两者滤波:

图片 15

精心相比较一下,双边滤波对边缘的保存效果比高斯滤波好太多了,这点从第贰幅图就能够知晓缘由了。

此外!!如若选取高斯核函数来兑现双方滤波,颜色卷积和的 \(\sigma\)
要取大学一年级点的值,比如:50。不然,由于分歧颜色的差值往往比地点差值大出不少(举个例子:50
和 60 两种像素值肉眼上看很相近,但却差出 10,平方一下便是100),或者引致很类似的像素点权值相当小,最终跟没滤波的功用同样。

3 双边滤波完毕步骤

  有了上述辩护之后完结Bilateral Filter就相比较简单了,其实它也与日常的Gaussian Blur没有太大的区分。那里根本总结3局地的操作: 基于空间距离的权重因子变化;基于相似度的权重因子的变型;最后filter颜色的持筹握算。

启发

Bilateral Filter
的考虑是:在岗位新闻的基本功上添加颜色音讯,约等于考虑多少个权值。借使还要考虑其余关键成分,是否足以再充实2个权值,构成1个三角形滤波器呢?答案自然是足以的,因此,大家得以把多如牛毛简练的滤波器综合起来形成1个更强有力的滤波器。

3.1Spatial Weight

  那正是熟视无睹的Gaussian Blur中利用的计量高斯权重的艺术,其根本透过七个pixel之间的相距并运用如下公式计算而来:

                 图片 16

参考

3.1Similarity Weight

  与基于距离的高斯权重总结类似,只可是此处不再依据多少个pixel之间的上空中距离离,而是基于其相似程度(只怕四个pixel的值时期的距离)。

            图片 17

   个中的表示三个像素值之间的距离,可以直接使用其灰度值之间的差值可能CR-VGB向量之间的欧氏距离。

3.4 Color Filtering

   当中的相似度分部的权重s首要基于四个Pixel之间的水彩差值总括面来。对于灰度图而言,这几个差值的限制是足以预见的,即[-255, 255],因此为了做实计算的频率大家能够将该部分权重因子预总结生成并存表,在选拔时飞速查询即可。使用上述达成的算法对几张带有噪声的图像举办滤波后的结果如下所示:

 3 双边滤波源码达成

 1 #include <time.h>
 2 #include <stdio.h>
 3 #include <stdlib.h>
 4 #include <math.h>
 5 #include "bf.h"
 6 
 7 
 8 //--------------------------------------------------
 9 /**
10 双边滤波器图像去噪方法
11 \param   pDest       去噪处理后图像buffer指针
12 \param   pSrc        原始图像buffer指针
13 \paran   nImgWid     图像宽
14 \param   nImgHei     图像高
15 \param   nSearchR    搜索区域半径
16 \param   nSigma      噪声方差 
17 \param   nCoff       DCT变换系数个数
18 
19 return   void  
20 */
21 //--------------------------------------------------
22 void BilateralFilterRGB(float *pDest, BYTE *pSrc,int nImgWid, int nImgHei, int nSearchR, int nSigma, int nCoff)
23 {
24 
25     BYTE *pSrcR  = new BYTE[nImgHei * nImgWid];                  // 新建图像缓冲区R
26     BYTE *pSrcG  = new BYTE[nImgHei * nImgWid];                  // 新建图像缓冲区G
27     BYTE *pSrcB  = new BYTE[nImgHei * nImgWid];                  // 新建图像缓冲区B
28 
29 
30     for(int i = 0;i < nImgHei; i++)
31     {
32         for (int j = 0; j < nImgWid; j++)
33         {
34             pSrcR[i * nImgWid + j] = pSrc[i * nImgWid * 3 + j * 3 + 0];         
35             pSrcG[i * nImgWid + j] = pSrc[i * nImgWid * 3 + j * 3 + 1];         
36             pSrcB[i * nImgWid + j] = pSrc[i * nImgWid * 3 + j * 3 + 2];         
37         }
38     }    
39 
40 
41     float* pDestR = new float[nImgHei * nImgWid];         // 去噪后的图像数据R通道
42     float* pDestG = new float[nImgHei * nImgWid];         // 去噪后的图像数据G通道
43     float* pDestB = new float[nImgHei * nImgWid];         // 去噪后的图像数据B通道
44     memset(pDestR,255,nImgHei * nImgWid * sizeof(float)); // 传递变量初始化
45     memset(pDestG,255,nImgHei * nImgWid * sizeof(float));
46     memset(pDestB,255,nImgHei * nImgWid * sizeof(float));
47 
48 
49     float *pII = new float[nImgHei * nImgWid ];                       // 积分图
50     float *pWeights  = new float[nImgHei * nImgWid ];                       // 每一个窗的权重之和,用于归一化  
51 
52 
53     float mDctc[MAX_RANGE];                                     // DCT变换系数
54     float mGker[MAX_RANGE];                                     // 高斯核函数的值
55 
56     for (int i = 0; i< MAX_RANGE; ++i)
57     {
58         mGker[i] = exp(-0.5 * pow(i / nSigma, 2.0));                       // 首先计算0-255的灰度值在高斯变换下的取值,然后存在mGker数组中
59     }
60 
61     CvMat G = cvMat(MIN_RANGE, MAX_RANGE, CV_32F, mGker);           // 转换成opencv中的向量
62     CvMat D = cvMat(MIN_RANGE, MAX_RANGE, CV_32F, mDctc);   // 把mDctc系数数组转换成opencv中的向量形式
63 
64     cvDCT(&G,&D,CV_DXT_ROWS);                                           // 然后将高斯核进行离散的dct变换
65     mDctc[0] /= sqrt(2.0);                                              // 离散余弦变换的第一个系数是1/1.414;
66 
67 
68     bf(pSrcR, pDestR, mDctc, pII, pWeights, nImgHei, nImgWid, nCoff, nSearchR);
69     bf(pSrcG, pDestG, mDctc, pII, pWeights, nImgHei, nImgWid, nCoff, nSearchR);
70     bf(pSrcB, pDestB, mDctc, pII, pWeights, nImgHei, nImgWid, nCoff, nSearchR);
71 
72 
73     for(int i=0; i < nImgHei;i++)
74     {
75         for (int j=0; j < nImgWid; j++)
76         {
77             pDest[i * nImgWid * 3 + j * 3 + 0] = pDestR[i * nImgWid + j];           
78             pDest[i * nImgWid * 3 + j * 3 + 1] = pDestG[i * nImgWid + j];                   
79             pDest[i * nImgWid * 3 + j * 3 + 2] = pDestB[i * nImgWid + j];           
80         }
81     }
82 
83 
84     free(pSrcR);      // 释放临时缓冲区
85     free(pSrcG);
86     free(pSrcB);
87     free(pDestR);
88     free(pDestG);
89     free(pDestB);
90 
91 }

  1 //--------------------------------------------------
  2 /**
  3 双边滤波器图像去噪方法
  4 \param   pDest       去噪处理后图像buffer指针
  5 \param   pSrc        原始图像buffer指针
  6 \paran   nImgWid     图像宽
  7 \param   nImgHei     图像高
  8 \param   pII         积分图缓冲区
  9 \param   pWeights    归一化权重系数 
 10 \param   nCoff       DCT变换系数
 11 \param   nPatchWid   图像块宽度
 12 
 13 return   void  
 14 */
 15 //--------------------------------------------------
 16 void bf (BYTE *pSrc, float *pDest, float *pDctc, float *pII, float *pWeights, int nImgHei, int nImgWid, int nCoff, int nPatchWid)  
 17 {  
 18     if (pDest == NULL || pSrc ==NULL)
 19     {
 20         return;
 21     }
 22 
 23     float *cR  = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float));
 24     float *sR  = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float));
 25     float *dcR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float));
 26     float *dsR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float));
 27 
 28     
 29     //  余弦变换查找表
 30     int fx     = 0;
 31     int ck     = 0;
 32     int r2     = pow(2.0*nPatchWid + 1.0, 2.0);         // 这个是图像块的像素点的总个数
 33 
 34     float c0   = pDctc[0];                              // 这是DCT的第一个系数
 35     float c0r2 = pDctc[0] * r2;                         // 这个用于初始化归一化权重之和
 36 
 37 
 38     for (ck = 1; ck < nCoff; ck++)                            // 注意到这里面的系数并不包含C0,所以后面才有把C0额外的加上
 39     {                                                                     // 一个矩阵,用于存放各种系数和数值,便于查找
 40         int ckr = (ck - 1) * MAX_RANGE;                                 // 数组是从0开始的,所以这里减1
 41 
 42         for (fx = 0; fx<MAX_RANGE; fx++)                              // fx就是图像的灰度值
 43         {
 44             float tmpfx = PI * float(fx) * float(ck) / MAX_RANGE;   // ck其实就相当于那个余弦函数里面的t,fx相当于u,PI/MAX相当于前面的那个系数,这都是余弦变换相关的东西
 45             
 46             cR[ckr + fx]  = cos(tmpfx);                                      // 存储余弦变换,这个用在空间域上      
 47             sR[ckr + fx]  = sin(tmpfx);                                      // 存储正弦变换
 48             dcR[ckr + fx] = pDctc[ck]  * cos(tmpfx);                         // 存储余弦变换和系数的乘积,这个则用在强度范围上
 49             dsR[ckr + fx] = pDctc[ck]  * sin(tmpfx);                         // 存储正弦变换和系数的乘积       
 50         }        
 51     }
 52 
 53       
 54     float *pw  = pWeights;                          // 新建一个归一化权重的中间变量进行数据的初始化
 55     float *pwe = pWeights + nImgWid * nImgHei;      // 限定最大位置
 56 
 57     while (pw < pwe) 
 58     {
 59         *pw++ = c0r2;                   // 赋初值,让它们都等于第一个DCT系数乘以图像块中总的像素点的个数
 60     }
 61     
 62     for (int ck = 1; ck < nCoff; ck++) 
 63     {        
 64         int ckr = (ck-1)*MAX_RANGE; // 数组是从0开始的,所以这里减1
 65 
 66         add_lut(pWeights, pSrc, pII,cR + ckr, dcR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid);  // add cos term to pWeights
 67         add_lut(pWeights, pSrc, pII,sR + ckr, dsR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid);  // add cos term to pWeights  
 68         
 69     } 
 70 
 71     ImgRectSumO(pDest, pSrc, pII, c0, nImgHei, nImgWid, nPatchWid, nPatchWid);  //加上C0的变换之后,再初始化滤波后的函数     
 72     
 73     for (int ck = 1; ck < nCoff; ck++) 
 74     {        
 75         int ckr = (ck-1)*MAX_RANGE;
 76         
 77         add_f_lut(pDest, pSrc, pII, cR + ckr, dcR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid);  // add cosine term to dataf
 78         add_f_lut(pDest, pSrc, pII, sR + ckr, dsR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid);  // add sine term to dataf
 79   
 80     }
 81 
 82 
 83     float *pd = pDest + nPatchWid * nImgWid;
 84     pw        = pWeights + nPatchWid * nImgWid;
 85 
 86     float *pdie  = pd + nImgWid - nPatchWid;
 87     float *pdend = pDest + (nImgHei - nPatchWid) * nImgWid;
 88     
 89     while (pd < pdend) 
 90     {
 91         pd += nPatchWid; //把边都给去掉了
 92         pw += nPatchWid; //这个也是把边去掉了
 93 
 94         while (pd < pdie) 
 95         {
 96             *pd++ /= ( (*pw++)*255);     // 之所以要除以255,是为了把它化到[0,1]之间,便于显示
 97         }
 98 
 99         pd += nPatchWid;                 // 把边都给去掉了
100         pw += nPatchWid;                 // 这个也是把边去掉了
101         pdie += nImgWid;                 // 过渡到下一行       
102     }
103    
104     free(cR); 
105     free(sR); 
106     free(dcR); 
107     free(dsR);   // 释放缓冲区
108 }

//--------------------------------------------------
/**
余弦函数首系数积分图
\param   pDest       去噪处理后图像buffer指针
\param   pSrc        原始图像buffer指针
\paran   nImgWid     图像宽
\param   nImgHei     图像高
\param   pII         积分图缓冲区
\param   c0          DCT变换第一个系数
\param   nPatchWid   图像块宽度
\param   nPatchHei   图像块高度

return   void  
*/
//--------------------------------------------------
void ImgRectSumO (float* pDest, const BYTE* pSrc, float* pII, float c0,const int nImgHei, const int nImgWid, const int nPatchWid, const int nPatchHei) 
{
    if (pDest == NULL || pSrc ==NULL)
    {
        return;
    }

    int ri1 = nPatchWid + 1;                                // 中心像素点
    int rj1 = nPatchHei + 1;
    int ri21 = 2 * nPatchWid + 1;                           // 图像块的宽度
    int rj21 = 2 * nPatchHei + 1;

    const BYTE *pSrcImg    = pSrc;                          // 原始图像数据
    const BYTE *pSrcImgEnd = pSrc + nImgHei * nImgWid;      // 最大的像素点位置    
    const BYTE *pSrcImgWid = pSrcImg + nImgWid;             // 第一行的最大的像素点位置,下面循环用的的变量
    float       *pIITmp     = pII;                           // 积分图数据
    float       *ppII       = pIITmp;                        // 积分图中间变量


    *pIITmp++ = c0 * float( *pSrcImg++ );                    // 第一个系数乘以图像数据

    while (pSrcImg < pSrcImgWid)                             // 遍历第一行进行求积分图,其实这个是图像数据的积分图
    {
        *pIITmp++ = (*ppII++) + c0 * float(*pSrcImg++);      
    }  

    while (pSrcImg < pSrcImgEnd)                             // 遍历所有的像素点的变量
    {    

        pSrcImgWid   += nImgWid;                             // 进行第二行的循环
        ppII          = pIITmp;                              // 积分图中间变量
        *pIITmp++     = c0 * float(*pSrcImg++);              // 第二行第一个像素点的处理

        while (pSrcImg < pSrcImgWid) 
        {
            (*pIITmp++) = (*ppII++) + c0 * float(*pSrcImg++); // 求第二行的积分图
        }        

        float *pIIWid = pIITmp;                        // 当前位置像素点
        pIITmp = pIITmp - nImgWid;                     // 上一行的起始点位置
        float *pIIup  = pIITmp - nImgWid;              // 上上一行的起始点位置

        while (pIITmp < pIIWid)                        // 遍历整个行的每一个数据
        {
            (*pIITmp++) = (*pIITmp) + (*pIIup++);      // 行与行之间的积分图
        }

    }


    float *pres = pDest + ri1 * nImgWid;                    // 最小的行位置
    float *pend = pDest + (nImgHei - nPatchWid) * nImgWid;  // 最大的行位置

    float *pii1 = NULL;
    float *pii2 = NULL;
    float *pii3 = NULL;
    float *pii4 = NULL;                               // 求积分图的四个点

    pii1 = pII + ri21 * nImgWid + rj21;               // 右下角
    pii2 = pII + ri21 * nImgWid;                      // 左下角
    pii3 = pII + rj21;                                // 右上角
    pii4 = pII;                                       // 左上角

    while (pres < pend) 
    {
        float *pe = pres + nImgWid - nPatchHei;
        pres += rj1;

        while (pres < pe)                             // 这个只是求图像数据的积分图
        {
            (*pres++) = (*pii1++) - (*pii2++) - (*pii3++) + (*pii4++);
        }

        pres += nPatchHei;
        pii1 += rj21;
        pii2 += rj21;
        pii3 += rj21;
        pii4 += rj21;
    }

}

//--------------------------------------------------
/**
余弦积分图加函数
\param   pNomalWts   像素点的归一化权重矩阵
\param   pSrc        原始图像buffer指针
\param   pII         积分图
\paran   nImgWid     图像宽
\param   nImgHei     图像高
\param   pII         积分图缓冲区
\param   pCosMtx     余弦函数矩阵
\param   pCosDctcMtx 余弦函数与DCT变换系数乘积矩阵
\param   nPatchWid   图像块宽度
\param   nPatchHei   图像块高度

return   void  
*/
//--------------------------------------------------

void add_lut (float* pNomalWts, const BYTE* pSrc, float* pII, float* pCosMtx, float* pCosDctcMtx, const int nImgHei, const int nImgWid, const int nPatchWid, const int nPatchHei) 
{ 

    int ri1  = nPatchWid + 1;      // kernel相关
    int rj1  = nPatchHei + 1;
    int ri21 = 2 * nPatchWid + 1;  // 这是kernel的宽度和长度
    int rj21 = 2 * nPatchHei + 1;

    const BYTE *pi  = pSrc;           // 这个是图像数据
    float *pii       = pII;            // 这个是中间的缓冲变量,积分图的缓冲区
    const BYTE *piw = pi + nImgWid;   // 也是赋初值的时候使用的中间变量
    float *pii_p     = pii;            // 直线积分图指针的指针


    *pii++ = pCosMtx[*pi++];  //图像数据值所对应的查找表中的余弦或者正弦变换


    while (pi < piw)
    {            
        *pii++ = (*pii_p++) + pCosMtx[*pi++]; // 第一行的积分图
    }    

    const BYTE *piend = pSrc + nImgHei * nImgWid;       // 限定循环位置

    while (pi < piend)                                     // 遍历整个图像数据
    {  

        piw    += nImgWid;                            // 第二行
        pii_p   = pii;                             // 指向积分图指针的指针
        *pii++  = pCosMtx[*pi++];              // 获取第一个图像数据值所对应的查找表中的余弦或者正弦变换

        while (pi < piw) 
        {
            (*pii++) = (*pii_p++) + pCosMtx[*pi++]; // 求这一行的积分图
        }

        float *piiw = pii;
        pii = pii - nImgWid;
        float *pii_p1 = pii - nImgWid;

        while (pii < piiw) 
        {
            (*pii++) = (*pii) + (*pii_p1++);                  // 行与行之间求积分图
        }

    }  


    float *pres = pNomalWts + ri1 * nImgWid;                   // 定位要处理点的像素的那一行
    float *pend = pNomalWts + (nImgHei - nPatchWid) * nImgWid; // 限定了边界

    float *pii1 = NULL;
    float *pii2 = NULL;
    float *pii3 = NULL;
    float *pii4 = NULL;

    pii1 = pII + ri21 * nImgWid + rj21;  // 获得积分图中那个最大的数据右下角
    pii2 = pII + ri21 * nImgWid;         // 积分图左下角
    pii3 = pII + rj21;                   // 积分图右上角
    pii4 = pII;                          // 积分图左上角

    pi = pSrc + ri1 * nImgWid;                  // 定位要处理的像素的位置


    while (pres < pend)                         // 设定高度做循环
    {
        float *pe = pres + nImgWid - nPatchHei; // 限定了宽度的范围
        pres = pres + rj1;                      // 定位了要处理的像素点的归一化系数矩阵的位置
        pi = pi + rj1;                          // 定位了要处理的像素点

        while (pres < pe)                       // 设定宽度做循环
        {
            (*pres++) = (*pres) + pCosDctcMtx[*pi++] * ( (*pii1++) - (*pii2++) - (*pii3++) + (*pii4++) );    // 得到的就是归一化系数矩阵之和
        }

        pres += nPatchHei;                     
        pi   += nPatchHei;                     
        pii1 += rj21;                          
        pii2 += rj21;                          
        pii3 += rj21;                          
        pii4 += rj21;                          // 略过不处理的边界区域
    }    

}

//--------------------------------------------------
/**
积分图
\param   pDest   像素点的归一化权重矩阵
\param   pSrc        原始图像buffer指针
\param   pII         积分图
\paran   nImgWid     图像宽
\param   nImgHei     图像高
\param   pII         积分图缓冲区
\param   pCosMtx     余弦函数矩阵
\param   pCosDctcMtx 余弦函数与DCT变换系数乘积矩阵
\param   nPatchWid   图像块宽度
\param   nPatchHei   图像块高度

return   void  
*/
//--------------------------------------------------

void add_f_lut (float* pDest,const BYTE* pSrc,float* pII,float* pCosMtx,float* pCosDctcMtx,const int nImgHei, const int nImgWid,const int nPatchWid, const int nPatchHei) 
{  

    int ri1 = nPatchWid + 1;      // kernel相关,目标是指向中心像素点
    int rj1 = nPatchHei + 1;      // 目标是指向中心像素点
    int ri21 = 2 * nPatchWid + 1; // 这是kernel的宽度和长度
    int rj21 = 2 * nPatchHei + 1; // 其实就是指向搜索区域的右下角,也就是最大位置的那个像素点,就是边界了

    const BYTE *pi = pSrc;       // 这个是图像数据,图像的灰度值
    float *pii = pII;             // 这是积分图

    const BYTE *piw = pi + nImgWid;  // 指向第一行的末尾位置,用于循环控制
    float *pii_p = pii;               // 指向积分图指针的指针

    *pii++ = float(*pi) * pCosMtx[*pi]; // 灰度值乘以余弦系数
    ++pi;      

    while (pi < piw)                    // 第一行遍历
    {
        *pii++ = (*pii_p++) + float(*pi) * pCosMtx[*pi];
        ++pi;
    }

    const BYTE *piend = pSrc + nImgHei * nImgWid;

    while (pi < piend) 
    {        
        piw   += nImgWid;
        pii_p  = pii;
        *pii++ = float(*pi) * pCosMtx[*pi];
        ++pi;

        while (pi < piw)
        {
            (*pii++) = (*pii_p++) + float(*pi) * pCosMtx[*pi];
            ++pi;
        }

        float *piiw = pii;
        pii = pii - nImgWid;   // 上一行起点
        float *pii_p1 = pii - nImgWid; // 上上一行的起始点

        while (pii < piiw)             // 循环一行
        {
            (*pii++) += (*pii_p1++);  //其实就是在列的位置上加上上一行
        }

    }


    float *pres = pDest + ri1*nImgWid;                     
    float *pend = pDest + (nImgHei - nPatchWid) * nImgWid; 
    float *pii1 = NULL;
    float *pii2 = NULL;
    float *pii3 = NULL;
    float *pii4 = NULL;

    pii1 = pII + ri21 * nImgWid + rj21;   // 积分图搜索区域最大位置的元素,堪称右下角的元素
    pii2 = pII + ri21 * nImgWid;          // 可以看成搜索区域中左下角的像素位置
    pii3 = pII + rj21;                    // 搜索区域右上角像素位置
    pii4 = pII;                           // 搜索区域左上角像素位置
    pi   = pSrc + ri1 * nImgWid;          // 定位要处理的像素的位置的那一行

    while (pres < pend) 
    {
        float *pe = pres + nImgWid - nPatchHei; //限定了宽度的范围

        pres = pres + rj1; //定位了要处理的像素点的归一化系数矩阵的位置
        pi   = pi   + rj1;     //定位了要处理的像素点

        while (pres < pe)  //遍历整个行
        {
            (*pres++) = (*pres) + pCosDctcMtx[*pi++] * ( (*pii1++) - (*pii2++) - (*pii3++) + (*pii4++) ); //这个其实是计算整个像素块
        }

        pres += nPatchHei;   
        pi   += nPatchHei;   

        pii1 += rj21; 
        pii2 += rj21;
        pii3 += rj21;
        pii4 += rj21;
    }
}