保存边界的平滑滤波器,具备简易、非

多年来在看图像风格化的舆论的时候,频仍碰到 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

   当中的意味七个像素值之间的距离,能够一贯利用其灰度值之间的差值也许奥迪Q7GB向量之间的欧氏距离。

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

 

相关文章