但它的合计却很有借鉴意义,双边滤波器

目前在看图像风格化的杂文的时候,频繁遭逢 Bilateral Filter。google
一波后,发现并不是何等不可了的事物,但它的思索却很有借鉴意义。

Bilateral Filtering for Gray and Color Images

 

双方滤波器:保留边界的平滑滤波器。
在有些上,就是在灰度值差别不大的区域平滑,在灰度值差异相比较大的边界地区保留边界。所以双方滤波器效能于各个像素的同时,必然会遭到世界像素点的离开、灰度值差的权重影响。

 

已知低通滤波可以代表为:

图片 1

图片 2

 

range filter可以象征为:(range filter
试选定一个数值范围,再做滤波的一个操作)

图片 3

图片 4

 

故而,双边滤波器的概念是:

图片 5

中间,k(x)是归一化(normalize)函数,

图片 6

( f 表示原图像,h 表示处理后的图像 x 表示 h 中某个像素点地方,ξ 表示 f
中x地点像素点的邻域像素,f(ξ)表示该像素点的灰度值,c表示低通滤波,
s表示range filter)

其中,

图片 7

图片 8

 

图片 9

 

 

//Filters.h

#ifndef FILTERS_H
#define FILTERS_H

#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/core.hpp"

#include <iostream>
#include <cmath>

//Bilateral Filtering
//sigmaD == sigmaSpace, sigmaR == sigmaColor
cv::Mat BilateralFilter(cv::Mat inputImg, int filterSize, double sigmaD, double sigmaR);

cv::Mat fastBilateralFilter(cv::Mat inputImg, int filterSize, double sigmaD, double sigmaR);

#endif // ! FILTERS_H

 

//Filters.cpp

#include "Filters.h"

double SpaceFactor(int x1, int y1, int x2, int y2, double sigmaD) {
    double absX = pow(abs(x1 - x2), 2);
    double absY = pow(abs(y1 - y2), 2);

    return exp(-(absX + absY) / (2 * pow(sigmaD, 2)));
}

double ColorFactor(int x, int y, double sigmaR) {
    double distance = abs(x - y) / sigmaR;
    return exp(-0.5 * pow(distance, 2));
}

cv::Mat BilateralFilter(cv::Mat inputImg, int filterSize, double sigmaD, double sigmaR) {
    int len; //must be odd number
    cv::Mat gray; // must be 1-channel image
    cv::Mat LabImage; // if channels == 3

    if (filterSize % 2 != 1 || filterSize <= 0) {
        std::cerr << "Filter Size must be a positive odd number!" << std::endl;
        return inputImg;
    }
    len = filterSize / 2;

    if (inputImg.channels() >= 3) {
        cv::cvtColor(inputImg, LabImage, cv::COLOR_BGR2Lab);
        gray = cv::Mat::zeros(LabImage.size(), CV_8UC1);
        for (int i = 0; i < LabImage.rows; i++) {
            for (int j = 0; j < LabImage.cols; j++) {
                gray.ptr<uchar>(i)[j] = LabImage.ptr<uchar>(i, j)[0];
            }
        }
    }
    else if(inputImg.channels() == 1){
        inputImg.copyTo(gray);
    }
    else {
        std::cerr << "the count of input image's channel can not be 2!" << std::endl;
        return inputImg;
    }

    cv::Mat resultGrayImg = cv::Mat::zeros(gray.size(), CV_8UC1);
    for (int i = 0; i < gray.rows; i++) {
        for (int j = 0; j < gray.cols; j++) {
            double k = 0;
            double f = 0;
            for (int r = i - len; r <= i + len; r++) {
                for (int c = j - len; c <= j + len; c++) {
                    if (r < 0 || c < 0 || r >= gray.rows || c >= gray.cols)
                        continue;
                    f = f + gray.ptr<uchar>(r)[c] * SpaceFactor(i, j, r, c, sigmaD) * ColorFactor(gray.ptr<uchar>(i)[j], gray.ptr<uchar>(r)[c], sigmaD);
                    k += SpaceFactor(i, j, r, c, sigmaD) * ColorFactor(gray.ptr<uchar>(i)[j], gray.ptr<uchar>(r)[c], sigmaD);
                }
            }
            int value = f / k;
            if (value < 0) value = 0;
            else if (value > 255) value = 255;

            resultGrayImg.ptr<uchar>(i)[j] = (uchar)value;
        }
    }

    cv::Mat resultImg;
    if (inputImg.channels() >= 3) {
        for (int i = 0; i < LabImage.rows; i++) {
            for (int j = 0; j < LabImage.cols; j++) {
                LabImage.ptr<uchar>(i, j)[0] = resultGrayImg.ptr<uchar>(i)[j];
            }
        }
        cv::cvtColor(LabImage, resultImg, cv::COLOR_Lab2BGR);
    }
    else {
        resultGrayImg.copyTo(resultImg);
    }

    return resultImg;
}

cv::Mat fastBilateralFilter(cv::Mat inputImg, int filterSize, double sigmaD, double sigmaR) {
    int len; //must be odd number
    cv::Mat gray; // must be 1-channel image
    cv::Mat LabImage; // if channels == 3

    if (filterSize % 2 != 1 || filterSize <= 0) {
        std::cerr << "Filter Size must be a positive odd number!" << std::endl;
        return inputImg;
    }
    len = filterSize / 2;

    if (inputImg.channels() >= 3) {
        cv::cvtColor(inputImg, LabImage, cv::COLOR_BGR2Lab);
        gray = cv::Mat::zeros(LabImage.size(), CV_8UC1);
        for (int i = 0; i < LabImage.rows; i++) {
            for (int j = 0; j < LabImage.cols; j++) {
                gray.ptr<uchar>(i)[j] = LabImage.ptr<uchar>(i, j)[0];
            }
        }
    }
    else if (inputImg.channels() == 1) {
        inputImg.copyTo(gray);
    }
    else {
        std::cerr << "the count of input image's channel can not be 2!" << std::endl;
        return inputImg;
    }

    cv::Mat resultGrayImg = cv::Mat::zeros(gray.size(), CV_8UC1);
    for (int i = 0; i < gray.rows; i++) {
        for (int j = 0; j < gray.cols; j++) {
            double k = 0;
            double f = 0;
            double sum = 0;
            for (int r = i - len; r <= i + len; r++) {
                if (r < 0 || r >= gray.rows)
                    continue;
                f = f + gray.ptr<uchar>(r)[j] * SpaceFactor(i, j, r, j, sigmaD) * ColorFactor(gray.ptr<uchar>(i)[j], gray.ptr<uchar>(r)[j], sigmaD);
                k += SpaceFactor(i, j, r, j, sigmaD) * ColorFactor(gray.ptr<uchar>(i)[j], gray.ptr<uchar>(r)[j], sigmaD);
            }
            sum = f / k;
            f = k = 0.0;
            for (int c = j - len; c <= j + len; c++) {
                if (c < 0 || c >= gray.cols)
                    continue;
                f = f + gray.ptr<uchar>(i)[c] * SpaceFactor(i, j, i, c, sigmaD) * ColorFactor(gray.ptr<uchar>(i)[j], gray.ptr<uchar>(i)[c], sigmaD);
                k += SpaceFactor(i, j, i, c, sigmaD) * ColorFactor(gray.ptr<uchar>(i)[j], gray.ptr<uchar>(i)[c], sigmaD);
            }
            int value = (sum + f / k) / 2;
            if (value < 0) value = 0;
            else if (value > 255) value = 255;

            resultGrayImg.ptr<uchar>(i)[j] = (uchar)value;
        }
    }

    cv::Mat resultImg;
    if (inputImg.channels() >= 3) {
        for (int i = 0; i < LabImage.rows; i++) {
            for (int j = 0; j < LabImage.cols; j++) {
                LabImage.ptr<uchar>(i, j)[0] = resultGrayImg.ptr<uchar>(i)[j];
            }
        }
        cv::cvtColor(LabImage, resultImg, cv::COLOR_Lab2BGR);
    }
    else {
        resultGrayImg.copyTo(resultImg);
    }

    return resultImg;
}

//main.cpp

#include <iostream>
#include <time.h>

#include "Filters.h"

using namespace std;

int main() {
    cv::Mat img = cv::imread("Capture.jpg", cv::IMREAD_UNCHANGED);
    clock_t begin_time = clock();
    cv::Mat result = BilateralFilter(img, 15, 12.5, 50);
    std::cout << float(clock() - begin_time) / CLOCKS_PER_SEC << std:: endl;
    cv::imshow("original", result);
    cv::waitKey(0);
    cv::imwrite("original.jpg", result);

    begin_time = clock();
    result = fastBilateralFilter(img, 15, 12.5, 50);
    std::cout << float(clock() - begin_time) / CLOCKS_PER_SEC << std::endl;
    cv::imshow("fast", result);
    cv::waitKey(0);
    cv::imwrite("fast.jpg", result);

    begin_time = clock();
    cv::bilateralFilter(img, result, 15, 50, 12.5);
    std::cout << float(clock() - begin_time) / CLOCKS_PER_SEC << std::endl;
    cv::imshow("opencv", result);
    cv::waitKey(0);
    cv::imwrite("opencv.jpg", result);

    system("pause");
    return 0;
}

 

运作结果:

46.889s  5.694s  0.202s

 

二维算子降成五个一维算子之后,速度加快了有些,但是如故不如opencv的快,效果也比它差一些(No
more reinventing the wheel…)

从左至右:“小白化病”帅气原图、BilateralFilter处理结果、fastBilateralFilter处理结果、opencv接口处理结果

图片 10 图片 11 图片 12 图片 13

 

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
有一个弱点:各向同性(不明白这多少个精通对不对)。用这种滤波器,每个点受邻居的熏陶是同样的,尽管它跟邻居像素可能差得比较多,也会被乡邻「同化」(举个例子:边缘被「和谐」掉了)。因而,有人提议了
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)关系来扭转。通过高斯分布的曲线可以发现,离目的像素越近的点对最终结出的贡献越大,反之则越小。其公式化的讲述相似如下所述:

                     图片 14

   其中的c即为基于空间距离的高斯权重,而用来对结果举办单位化。

  高斯滤波在低通滤波算法中有不错的显现,可是其却有其它一个题目,这就是只考虑了像素间的长空地点上的关系,因而滤波的结果会丢掉边缘的消息。这里的边缘重假设指图像中首要性的两样颜色区域(比如青色的苍天,粉青色的头发等),而Bilateral就是在Gaussian blur中进入了另外的一个权重分部来化解这一题目。Bilateral滤波中对于边缘的保障通过下述表达式来实现:

                     图片 15

                    图片 16

  其中的s为基于像素间相似程度的高斯权重,同样用来对结果举办单位化。对互相举行整合即可以博得基于空间距离、相似程度综合考量的Bilateral滤波:

                  图片 17

  上式中的单位化分部综合了二种高斯权重于一块而拿到,其中的c与s总结可以详细描述如下:

                  图片 18

   且有图片 19

                   图片 20

   且有图片 21

  上述给出的表明式均是在空间上的极端积分,而在像素化的图像中自然不可能这么做,而且也没必要这样做,因此在使用前需要对其举行离散化。而且也不需要对此每个局部像素从整张图像上进展加权操作,距离超过一定水准的像素实际上对眼前的目的像素影响很小,可以忽略的。限定局部子区域后的离散化公就可以简化为如下格局:

                   图片 22

  上述申辩公式就结成了Bilateral滤波实现的基础。为了直观地打听高斯滤波与双方滤波的区分,大家得以从下列图示中看看依照。即使目的源图像为下述左右区域显著的蕴藏噪声的图像(由程序自动生成),绿色框的主导即为目的像素所在的职位,那么当前像素处所对应的高斯权重与三头权重因子3D可视化后的模样如后面两图所示:

                   图片 23

  左图为本来的噪声图像;中间为高斯采样的权重;右图为Bilateral采样的权重。从图中得以看出Bilateral加入了一般程度分部以后可以将源图像左边那一个跟当前像素差值过大的点给滤去,这样就很好地维持了边缘。为了更加形象地察看两者间的区别,使用Matlab将该图在两种不同方法下的冲天图3D绘制出来,如下:

                 图片 24

  上述三图从左到右依次为:双边滤波,原始图像,高斯滤波。从低度图中可以明确看出Bilateral和Gaussian二种艺术的界别,前者较好地涵养了边缘处的梯度,而在高斯滤波中,由于其在边缘处的生成是线性的,因此就采用连累的梯度显示出渐变的动静,而这表现在图像中的话就是境界的散失(图像的演示可见于后述)。         

代码实现

了解原理后,实现其实也很粗略,上边给出的伪代码基本是着力算法了。其它需要小心的是,如果是彩色图的话,需要对每个通道的颜色值举办滤波。

切实实现可以参见这篇博客:图像处理之双边滤波效果(Bilateral Filtering
for 格雷 and Color
Image)
,或者参考我自己的
demo,当然,我也只是将下面博客的
java 版改成 c++ 而已^0^。

付给几幅结果图:

原图

图片 25

高斯模糊

图片 26

只有用颜色音信滤波

图片 27

六头滤波:

图片 28

有心人相比较一下,双边滤波对边缘的保存效果比高斯滤波好太多了,这点从第三幅图就可以知晓缘由了。

其它!!倘诺选择高斯核函数来贯彻两岸滤波,颜色卷积和的 \(\sigma\)
要取大一点的值,比如:50。否则,由于不同颜色的差值往往比地方差值大出广大(举个例子:50
和 60 三种像素值肉眼上看很类似,但却差出 10,平方一下就是
100),可能引致很接近的像素点权值很小,最后跟没滤波的功效等同。

3 双边滤波实现步骤

  有了上述辩解之后实现Bilateral Filter就相比简单了,其实它也与一般的Gaussian Blur没有太大的区别。这里根本概括3局部的操作: 基于空间距离的权重因子变化;基于相似度的权重因子的变化;最后filter颜色的盘算。

启发

Bilateral Filter
的考虑是:在地方音信的底蕴上加上颜色消息,相当于考虑六个权值。假诺还要考虑任何重大元素,是不是可以再扩充一个权值,构成一个三角滤波器呢?答案自然是足以的,因而,大家可以把无数简单的滤波器综合起来形成一个更强硬的滤波器。

3.1Spatial Weight

  这就是平凡的Gaussian Blur中采取的推测高斯权重的法子,其利害攸关通过五个pixel之间的距离并应用如下公式总括而来:

                 图片 29

参考

3.1Similarity Weight

  与基于距离的高斯权重总计类似,只可是此处不再按照多少个pixel之间的空远距离,而是遵照其相似程度(或者四个pixel的值期间的相距)。

            图片 30

   其中的表示四个像素值之间的离开,可以直接运用其灰度值之间的差值或者RGB向量之间的欧氏距离。

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;
    }
}

 

相关文章