对两端网格的辩论做了须臾间民用的明亮,能够提供更加多的动态范围和图像细节公海赌船网址

一、引言

自家初次接触HDPRADO方面的文化,有描述不正确的地点烦请见谅。

为便于文章讲述,引用部分百度中的文章对HD瑞鹰图像进行简单的叙述。

高动态范围图像(High-Dynamic
Range,简称HDCR-V),比较普通的图像,能够提供越来越多的动态范围和图像细节,依据差其他暴露时间的LD奥迪Q3(Low-Dynamic
Range)图像,利用每一种暴光时间相对应最佳细节的LD福睿斯图像来合成最终HDOdyssey图像,能够更好的反映人真实环境中的视觉效果。

切切实实确实存在的亮度差,即最亮的物体亮度,和微小的实体亮度之比为108
而人类的眼眸所能看到的界定是105反正,不过一般的显示屏,照相机能代表的唯有256种不一样的亮度,计算一般的显示器,照相机能代表的唯有256种分化的亮度机在象征图象的时候是用8bit(256)级或16bit(65536)级来分别图象的亮度的,但那区区几百或几万不只怕重现真实自然的普照意况。HD奇骏文件是一种独特图形文件格式,它的每三个像素除了平常的昂CoraGB音讯,还有该点的实际上亮度音信。普通的图形文件各样象素惟有0
-255的灰度范围,那事实上是不够的。想象一下太阳的发光强度和2个纯黑的物体之间的灰度范围或许说亮度范围的距离,远远当先了25七个级别。由此,一张普通的白昼风景图片,看上去白云和日光大概都展现是一律的灰度/亮度,都以纯灰黄,但实质上白云和太阳之间其实的亮度不容许同样,他们中间的亮度差距是宏伟的。因而,普通的图形文件格式是很不确切的,远远没有记录到实际世界的实际上情况。而HDSportage格式则记录了很广范围内的亮度新闻。

而是最后,HD途乐图像要在显示屏中呈现,还是要求对其数量开始展览处理的,怎样处理即能足够利用这么些多少,又能使得图像的显示尽量不丢失细节,是多年来广大图像工作者商讨的第三。

简单易行的说,正是前几天有一堆离散的数额,数据的分布范围大概很广,怎么着把这几个离散的多寡隐射到0到255时期。

二 、相关算法的实现

最简易的自然是线性隐射,先算出离散数据的最大值和纤维值,然后将数据线性的拉升至0到255里边,那种直白的操作往往心有余而力不足获得知足的效益,会导致大气细节丢失,表现在视觉上便是一大块浅水绿只怕一大块深红的,如下图所示:

 公海赌船网址 1  
   公海赌船网址 2

     下面两幅图要么是暗部太暗,要么是亮部太亮,全部相比较度太强,导致细节消息大批量有失。

     针对这一标题,很几个人建议了好多非常不错的化解方案,比如依据全局操作符的,当中正文小编完成个中的基于神速双边滤波技术的HDPRADO展现进程。

    
本文对应的参照杂文地址: Fast
Bilateral Filtering for the Display of High-Dynamic-Range
Images

    
散文的细路相当粗略,首先她将本来的HDLX570数据分解成五个层:base layer 和
detail layer,然后降低base layer的相比较度,不更改detail
layer的数量,在将那两层合并。

    
当中:base layer的数据用 HDOdyssey原始数据开始展览双边滤波获取。

    
算法的总结流程入下所示:

    1、input
intensity= 1/61*(R*20+G*40+B)

    2、r=R/(input
intensity), g=G/input intensity, B=B/input intensity

    3、log(base)=Bilateral(log(input
intensity))

    4、log(detail)=log(input
intensity)-log(base)

    5、log
(output intensity)=log(base)*compressionfactor+log(detail) –
log_absolute_scale

    6、R
output = r*10^(log(output intensity)), etc.

    
上述进程中的变量compressionfactor,log_absolute_scale原文小编的提出取值为:

     
compressionfactor = targetContrast/(max(log(base)) – min(log(base)))
对于广大图像,targetContrast使用log(5)能博得较为理想的值。

     
而log_absolute_scale= max(log(base))*compressionfactor;

   
在进展两岸的时候,SigmaS一般取值为0.02*马克斯(Width,Height)相比较妥贴,而SigmaPRADO取值0.4较为理想(那里是指多少量化到了0-1里头的)。

  所以都取优化的参数,则上述进度能够自行进行。

   
小编提到上述log操作都以以10为底实行的,笔者认为以e为底实效也没啥差异的。

  三、效果

    依据那些思路编制造进度序后,确实能得到很科学的效能,比如上述两幅图像,遵照前边讲的参数取值,解码后得到的图像如下:

   
公海赌船网址 3 
公海赌船网址 4

   
可知,图像的底细较为圆满的反映出来了。

    当然,自动的参数不必然能疏通最好的效益,比如依然那两幅图,手工业选取一些参数,能够调出如下效果:

   
公海赌船网址 5 
公海赌船网址 6

  特别是第①幅图,很有种蒙太奇的感到。

   
在探望几张长出新在散文中的图像的结果:

   
公海赌船网址 7 
公海赌船网址 8       

   
公海赌船网址 9 
公海赌船网址 10

   公海赌船网址 11 公海赌船网址 12

 
 公海赌船网址 13 公海赌船网址 14

   
          线性解码图                                             
双边滤波解码图

 
有个别图线性解码啥都看不到,双边滤波解码后细节呈现的就很清楚了。

   
HDRAV4格式的固有数据的解码能够借助FreeImage来完成,FreeImage就像是已经讲那些数量量化到了0和1里边(不自然科学)。一段简单的兑现代码如下:

public Bitmap LoadHdrFormFreeImage(string FileName)
{
    Bitmap Bmp = null;
    FREE_IMAGE_FORMAT fif = FREE_IMAGE_FORMAT.FIF_UNKNOWN; ;
    if (FreeImage.IsAvailable() == true)
    {
        fif = FreeImage.GetFileType(FileName, 0);
        if (fif != FREE_IMAGE_FORMAT.FIF_HDR)
        {
            MessageBox.Show("不是Hdr格式的图像.");
            return null;
        }
        fif = FreeImage.GetFIFFromFilename(FileName);
        FIBITMAP Dib = FreeImage.Load(fif, FileName, FREE_IMAGE_LOAD_FLAGS.DEFAULT);
        uint Bpp = FreeImage.GetBPP(Dib);

        if (Bpp != 96)
        {
            MessageBox.Show("无法支持的Hdr格式.");
            FreeImage.Unload(Dib);
            return null;
        }
        uint Width = FreeImage.GetWidth(Dib);                        //  图像宽度
        uint Height = FreeImage.GetHeight(Dib);                      //  图像高度
        uint Stride = FreeImage.GetPitch(Dib);                       //  图像扫描行的大小,必然是4的整数倍
        IntPtr Bits = FreeImage.GetBits(Dib);

        float* Data = (float*)Bits;
        int Speed, Index;
        byte* Pixel;
        float Value;

        if (RawData != null) Marshal.FreeHGlobal((IntPtr)RawData);
        RawData = (float*)Marshal.AllocHGlobal((int)Width * (int)Height * 3 * sizeof(float));
        CopyMemory(RawData, Data, (int)Width * (int)Height * 3 * sizeof(float));

        Bmp = new Bitmap((int)Width, (int)Height, PixelFormat.Format24bppRgb);
        BitmapData BmpData = Bmp.LockBits(new Rectangle(0, 0, Bmp.Width, Bmp.Height), ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
        Pixel = (byte*)BmpData.Scan0;

        for (int Y = 0; Y < Height; Y++)
        {
            Speed = Y * BmpData.Stride;
            Index = Y * (int)Width * 3;
            for (int X = 0; X < Width; X++)
            {
                Value = (Data[Index + 2] * 255);
                if (Value > 255)
                    Value = 255;
                else if (Value < 0)
                    Value = 0;
                Pixel[Speed] = (byte)Value;
                Value = (Data[Index + 1] * 255);
                if (Value > 255)
                    Value = 255;
                else if (Value < 0)
                    Value = 0;
                Pixel[Speed + 1] = (byte)Value;
                Value = (Data[Index + 0] * 255);
                if (Value > 255)
                    Value = 255;
                else if (Value < 0)
                    Value = 0;
                Pixel[Speed + 2] = (byte)Value;
                Index += 3;
                Speed += 3;
            }
        }
        FreeImage.Unload(Dib);
        Bmp.UnlockBits(BmpData);
        Bmp.RotateFlip(RotateFlipType.RotateNoneFlipY);
        return Bmp;
    }
    else
        return null;
}

  以上接纳的是线性解码。

  附1个解码的调用程序:http://files.cnblogs.com/Imageshop/ReadHdrTest.rar 

     越多的源码可参照:http://people.csail.mit.edu/sparis/code/src/tone_mapping

                            http://people.csail.mit.edu/fredo/PUBLI/Siggraph2002/

    
一些常见的用于测试的HD酷威图像能够从那边下载:http://www.pauldebevec.com/Research/HDR/

     同样,那种 tone
mapping算法也足以用在平日的ENCOREGB图像上,效果如下所示:

 
 公海赌船网址 15  公海赌船网址 16

              原图                                    扩展base
layer的相比较度

  
公海赌船网址 17 
公海赌船网址 18

             原图                                    下跌base
layer的相比度

  
对于普通的完整偏暗的图像,该方法也能取得十分精粹的效果,一些测试结果如下:

  
公海赌船网址 19  
公海赌船网址 20

   公海赌船网址 21  
公海赌船网址 22

   公海赌船网址 23  
公海赌船网址 24

   
很多其余的算法也能起到接近上述的法力,可是他俩一般很简单生出halo现象。  

   
相信对于偏亮的常常照片,也得以有一样的拍卖能力(未找到合适的测试图片)。

  

 公海赌船网址 25

*********************************我:
laviewpbt   时间: 二〇一三.11.18   联系QQ:  33184777
 转发请保留本行消息************************

   

1.前言

        如今在看Deep Bilateral Learning for Real-Time Image
Enhancement
,是一篇用CNN完结实时图像增强的法门,能够在二弟大上快速完结HDQX56。小说中提到到双方网格(Bilateral
Grid)的构思,查阅了多量的文献,对两岸网格的争鸣做了弹指间私人住房的精晓。

1 双边滤波简介

  双边滤波(Bilateral filter)是一种非线性的滤波方法,是组成图像的长空邻近度和像素值相似度的一种折衷处理,同时考虑空域消息和灰度相似性,达到保边去噪的指标。具有简易、非迭代、局部的表征。

  双边滤波器的好处是足以做边缘保存(edge preserving),一般过去用的维纳滤波或者高斯滤波去降噪,都会较明显地混淆边缘,对于频仍细节的爱戴效率并不备受瞩目。双边滤波器顾名思义比高斯滤波多了二个高斯方差sigma-d,它是依照空间分布的高斯滤波函数,所以在边缘附近,离的较远的像素不会太多影响到边缘上的像素值,那样就保障了边缘附近像素值的保留。可是出于保存了过多的一再音讯,对于彩色图像里的一再噪声,双边滤波器不可见彻底的滤掉,只能够对于低频消息实行较好的滤波。

2.两岸滤波器的迅猛近似

       在对两岸网格做计算此前,先介绍一下五头滤波器的便捷近似方法(A
Fast Approximation of the Bilateral Filter using a Signal Processing
Approach
),它是双方网格的雏形,是八个将双方滤波器扩大到高维空间进行线性卷积的艺术。而两者网格的撰稿人在那篇小说的根底上,将高维空间数据映射成3D
array,进而建议了双方网格。

2 双边滤波原理

  滤波算法中,目的点上的像素值平时是由其所在地方上的周围的二个小部分邻居像素的值所决定。在2D高斯滤波中的具体贯彻便是对周围的一定限制内的像素值分别赋以不一致的高斯权重值,并在加权平均后收获当前点的末段结出。而那边的高斯权重因子是使用四个像素之间的半空距(在图像中为2D)关系来扭转。通过高斯分布的曲线能够发现,离目的像素越近的点对最终结出的贡献越大,反之则越小。其公式化的叙述相似如下所述:

                     公海赌船网址 26

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

  高斯滤波在低通滤波算法中有科学的呈现,不过其却有其余2个标题,那就是只考虑了像素间的半空中地点上的涉嫌,由此滤波的结果会丢掉边缘的新闻。那里的边缘主若是指图像中关键的不比颜色区域(比如水晶色的天幕,雪白的头发等),而Bilateral就是在Gaussian blur中投入了其它的3个权重分部来消除这一标题。Bilateral滤波中对此边缘的维持通过下述表明式来贯彻:

                     公海赌船网址 27

                    公海赌船网址 28

  当中的s为基于像素间相似程度的高斯权重,同样用来对结果开始展览单位化。对双边进行整合即能够收获基于空间距离、相似程度综合考虑衡量的Bilateral滤波:

                  公海赌船网址 29

  上式中的单位化分部综合了三种高斯权重于同台而取得,个中的c与s总结能够详细描述如下:

                  公海赌船网址 30

   且有公海赌船网址 31

                   公海赌船网址 32

   且有公海赌船网址 33

  上述给出的表明式均是在半空中上的然而积分,而在像素化的图像中本来不能够这么做,而且也没须求那样做,由此在动用前要求对其展开离散化。而且也不必要对此各样局地像素从整张图像上海展览中心开加权操作,距离超过一定水平的像素实际上对脚下的指标像素影响相当小,能够忽略的。限定局部子区域后的离散化公就能够简化为如下方式:

                   公海赌船网址 34

  上述辩驳公式就结成了Bilateral滤波达成的功底。为了直观地询问高斯滤波与多头滤波的差别,我们能够从下列图示中看到依照。假若目的源图像为下述左右区域泾渭鲜明的涵盖噪声的图像(由程序自动生成),紫红框的主题即为指标像素所在的地方,那么当前像素处所对应的高斯权重与互相权重因子3D可视化后的形象如后面两图所示:

                   公海赌船网址 35

  左图为原始的噪声图像;中间为高斯采样的权重;右图为Bilateral采集样品的权重。从图中得以看出Bilateral加入了貌似程度分部今后能够将源图像右侧那1个跟当前像素差值过大的点给滤去,这样就很好地保持了边缘。为了更加形象地观测两者间的分化,使用Matlab将该图在二种分化方法下的中度图3D绘制出来,如下:

                 公海赌船网址 36

  上述三图从左到右依次为:双边滤波,原始图像,高斯滤波。从中度图中能够鲜明看出Bilateral和Gaussian二种办法的区别,前者较好地涵养了边缘处的梯度,而在高斯滤波中,由于其在边缘处的变更是线性的,由此就应用连累的梯度呈现出渐变的事态,而那表以往图像中的话便是境界的丢失(图像的言传身教可知于后述)。         

2.1.双面滤波器

       
首先介绍一下什么样是两者滤波器,那里引用一些客人的计算博客:相互滤波算法原理

       
简单地说,双边滤波器是一种高斯滤波器的壮大。古板的高斯滤波器有一个高斯核,通过对空中相邻的像素点以高斯函数位权值取均值,完成对图像的坦荡。由于古板高斯滤波器只考虑了空荡荡的音讯,固然它能兑现对图像的管用平滑,但与此同时也搅乱了边缘音信。双边滤波器的原理即在观念高斯滤波器的基础上添加了一个特色亮度差其他高斯核,即既考虑了空中的音信,由考虑了值域的音信。在灰度差异一点都不大的限制,表征亮度高斯核的权值较大、接近于1,由此双方滤波退化为古板的高斯滤波器,实现对图像的平缓;在边缘部分,即使相邻像素点的长空中接力近,空间距离小,但鉴于边缘部分,灰度差距较大,因而在值域上相差较大,所以加入的新的高峰斯核使得在该部分不进行平滑处理,保留图像的边缘。所以说,双边滤波是一种非线性(三个高斯核的乘积)的滤波方法,是整合图像的半空中邻近度和像素值相似度的一种折中拍卖,同时考虑空域音信和灰度相似性,达到保边去噪的指标。

3 双边滤波达成步骤

  有了上述辩白之后完结比拉teral Filter就比较简单了,其实它也与常见的Gaussian Blur没有太大的界别。那里最首要不外乎3有个别的操作: 基于空间距离的权重因子变化;基于相似度的权重因子的变通;最终filter颜色的计量。

2.2.高效近似

公海赌船网址 37

互相互连网接近例子

       
小编首先将两端滤波器的计算公式写成了齐次坐标的款式(关于齐次坐标的敞亮,能够参见关于齐次坐标的领悟(经典)):

公海赌船网址 38

两岸滤波齐次表示

       
接着,小编引入δ函数,将公式从二维空间扩充到三维空间(即2个2D的空间域和二个1D的值域/亮度域):

公海赌船网址 39

引入δ函数

        小编引入一些函数来代表上式中的一些操作:

公海赌船网址 40

高斯核

公海赌船网址 41

δ函数

        最后,将二者滤波器的公式改写为:

公海赌船网址 42

二者滤波的线性方式

公海赌船网址 43

一字不苟后的双方滤波进度

       
作者建议,由于采集样品定理,能够通过对图像下采集样品做卷积处理,再上采集样品苏醒原式分辨率,达到快捷的目标。

公海赌船网址 44

一部分标号的证实

3.1Spatial Weight

  这就是常见的Gaussian Blur中央银行使的猜想高斯权重的不二法门,其首要透过七个pixel之间的偏离并接纳如下公式总括而来:

                 公海赌船网址 45

2.4.伪代码

       
下图为双方近似的伪代码。1.初步化全部的下采集样品齐次值,令全数下采集样品权值为0。2.找到全图的细微亮度点(这一步是为了继续将亮度域的范围平移到从0开头)。3.对于高分辨率图像的任贰个像素点(X,Y)和亮度I(X,Y),(a)同样总计其齐次格局;(b)总结下采集样品后的坐标(那里运用了取整,因而,在原式高分辨率图像中的相邻空间,会被映射到低分辨率下的同一空间坐标,但在亮度域可能不在同一坐标上);(c)更新下采集样品空间的坐标值,即在低分辨率空间上添加高分辨率的齐次值(之所以用齐次值,是因为齐次坐标相加后,取出原始数据即完毕了对相应数据的加权平均:(w1i1,w1)+(w2i2,w2)=(w1i1+w2i2,w1+w2),数据对应(w1i1+w2i2)/(w1+w2))。4.在采集样品域做高斯卷积,得到卷积后的齐次坐标。5.上采集样品回高分辨率图像,对于未知的数据点,使用三线性差值。

公海赌船网址 46

四头近似伪代码

3.1Similarity Weight

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

3.双面网格

       
有大家计算了两者滤波器的类似方法,使用一种3D数组的款型来代表那种在三维空间的双边滤波,建议了五头网格的申辩(Real-time
edge-aware image processing with the bilateral
grid
)。

公海赌船网址 47

二者网格的相互滤波器达成原理

       
有了前方双边滤波近似的争论,再来掌握两者网格就显示不难得多了。考虑上海教室所示的1D输入图像。双边网格是在空间域和亮度域进行采样,划分成网格。双边的边也便是从那里来的,空间(space)和亮度(range)。离散后,每一个点的坐标和亮度消息取整到相应的网格内。各个值域内的亮度值可由加权平均得到。通过在网格内开始展览滤波等处理,再由上采集样品方法插值出处理后的固有图像。

       
因而,营造1个双面网格的格局是:首先,依据空域和值域,划分出网格结构,并早先化全部网格的节点为0:

公海赌船网址 48

开首化节点

       
对于高分辨率图像的坐标,分别除以空间采样间隔和亮度采集样品间隔并取整,再对应的节点处填入其原有的灰度值。考虑到取整函数会招致原图像的某一邻域内的像素点会下采集样品到同3个网格中,因而将灰度值写成齐次值(I(x,y),1)。并对在该网格内的拥有齐次值叠加。那样,这一个网格对应的末梢灰度最能够写成那一个灰度值的平均值:(I(x1,y1)+…+I(xn,yn))/n。

公海赌船网址 49

网格填充

       
那样,我们就将原本的图像数据转载为了3个3D的两边网格。任何图像操作函数都足以成效在这几个两岸网格上,即一对一于将该双边网格左乘三个3D图像操作矩阵(例如,对于两者滤波器而言,那个函数是2个高斯卷积核,包涵了空中的方差和亮度的方差),获得1个在低分辨率情状下处理的结果。最后通过上采集样品,并展开插值,获得最终于高分辨率图像结果。上采样的法则是挑选1个参考图,对其私下几个上空的像素举办空域和值域的采集样品,那里不举办取整操作,然后找到其在网格中的地点,由于尚未了取整,必要利用三线性插值的法门,来促成未知范围的亮度值的持筹握算,这几个进程称作“slicing”。由插值上采集样品后,大家获取了最终的滤波结果。

            公海赌船网址 50

   当中的表示多个像素值之间的离开,能够直接使用其灰度值之间的差值可能HighlanderGB向量之间的欧氏距离。

4.相互网格的上采集样品

        建议两岸网格的撰稿人接着又开始展览了和谐的三头网格,公布了Bilateral
Guided
Upsampling
,介绍了何等使用双边网格达成部分图像操作算子。算法的核心理想就是将一副高级分辨率的图像通过下采集样品转换到3个相互网格,而在二者网格中,每一个cell里提供四个图像变换算子,它的法则是在半空中与值域相近的区域内,相似输入图像的亮度经算子变换后也应该是形似的,因而在种种cell里的操作算子可以当作是输入/输出的接近曲线,也即2个仿射模型,而对于差其余cell,笔者通过给定输入和愿意输出去磨炼那一个双方网格完成其仿射模型的大局分段平滑。再经过上采集样品,获得高分辨率的拍卖后的图像。

公海赌船网址 51

Bilateral Guided Upsampling

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

 

相关文章