小波去噪算法的简易实现及其扩展（小波锐化、高斯拉普拉斯金字塔去噪及锐化）之二。

　　上一篇文章谈及了GIMP里实现的小波分解，但是这仅仅是把图像分解为多层的数据，如果快速的获取分解数据以及后续怎么利用这些数据，则是本文的重点。

　　 一、我们先来看看算法速度的优化问题。

　　原始的GIMP实现需要将图像数据转换为浮点数后，然后进行各级的模糊和图层混合，这样得到的结果是比较精确的，但是存在两个方面的问题，一个是占用了较多的内存，因为GIMP这个版本的小波分解各层是没有改变数据的尺寸的，因此，如果使用浮点，占用的内存要比字节版本的大四倍，而且和层数有着密切的关系。第二个是浮点的处理还是稍微慢了点，虽然对现在的CPU来说，浮点数更易用SIMD指令集优化。但是如果有更好的数据类型的话，使用SIMD可以获得更高的计算速度。

　　我们知道，字节版本的模糊可以获得很高的计算效率，这个场景的模糊是否可以使用呢，我个人分析认为，这个版本的算法已经不适合用字节版本的模糊了，原因主要是精度太低了，因为他每次的结果都和上一次的模糊相关，而我们通过GIMP的可是化界面可以看到，中间的每一个层很多像素都是靠近127的，这就说明细节方面的信息都是很小的数值。

　　个人觉得，对于这个算法，我们可以把数据放大到unsigned short范围，比如把原始的像素值都扩大256倍，然后进行模糊，这样模糊的精度就会大幅的提高，比如原始的9的像素的加权累加值（归一化后的权重）如果是100.3，那么由于字节版本取整的原因，最后的返回值为100，放大256倍后，则累加值变为25676.8，四舍五入取整后则为25677，而非100放大256倍的25600了，这样多次模糊的精度就可以加以充分的保证（和浮点数比较还是有差异，但不影响结果）。

　　我们在稍微观察下上一篇文章我们所提到的3*3的权重：

　　如果我们每个系数都放大16倍，我们会得到非常优化的一组权重系数：

　　可以看到，权重系数里全是1、2、4，这种系数如果被乘数是整形的话，都可以直接用移位实现，而放大16倍最后的除法，也可以直接由右移实现，那这种计算过程就完全避免了乘法，效率可想而知可以得到极大的提高。

　　我们在考虑一种特别的优化，因为权重整形化后的累加值是16，那么如果每个元素的最大值不超过4096，则累计值也就不会超过65536，这个时候如果是用普通C语言实现，其实没有啥区别，但是我们知道SIMD指令确有所不同，他针对不同的数据类型有不同的乘法和加法指令，其所能计算的覆盖范围也不同，如果能用16位表达，则尽量用16位表达。因此，考虑图像数据的特殊性，如果我们只把他们的数据范围扩大16倍，则不会超出4096的范围，就可以满足前面的这个假设。

　　放大16倍是否能满足精度的需求呢，没有具体理论的分析啊，个人觉得啊，至少在层数不大于4层时，差异不会很大，大于4层，也应该可以接受，留待实践检测吧。

　　对于上一篇文章所说到的取样点位置的问题，我们必须考虑边缘位置处的信息，因为随着取样范围的扩大，会有更多的取样点超出图像的有效范围，这个时候一个简单的办法就是实现准备好一副扩展过边界的图像，这里的扩展方法通道都选择边缘镜像，而非重复边缘。考虑到层数不会太多，扩展部分的边缘计算量和占用内存也不会太夸张。

　　同时，我们在考虑到算法的特殊性，虽然取样范围很广，但是真正用到的取样点也只有9个，所以我们也可以使用类似于我在SSE图像算法优化系列九：灵活运用SIMD指令16倍提升Sobel边缘检测的速度（4000*3000的24位图像时间由480ms降低到30ms）一文中使用到的技术，分配三行临时的内存，每次计算前填充好三行的数据，不过注意一点，由于取样的三行的内存并不是在行方向上连续的，因此，也不能像那个文章那样，可以重复利用两行的内存了，而必须每次都填充。

　　理论上说，这种填充技术要比前面的分配一整片的内存的内存填充工作量大3倍左右，速度应该要慢一些，优点是占用的中间内存少一些，但是实测，还是这种的速度快一些，或许是因为三行内存的大小有效，访问时的cache miss要好很多吧。

　　一个简答的代码片段如下所示：

    for (int Y = 0; Y < Height; Y++)
    {
        short *LinePD = Dest + Y * Width * Channel;
        int PosY0 = ColOffset[Y] * Width * Channel, PosY1 = Y * Width * Channel, PosY2 = ColOffset[Y + Radius + Radius] * Width * Channel;
        for (int X = 0, Index = 0; X < Radius; X++)
        {
            int PosX = RowOffset[X] * Channel;
            memcpy(First + Index, Src + PosY0 + PosX, Channel * sizeof(unsigned short));
            memcpy(Second + Index, Src + PosY1 + PosX, Channel * sizeof(unsigned short));
            memcpy(Third + Index, Src + PosY2 + PosX, Channel * sizeof(unsigned short));
            Index += Channel;
        }
        memcpy(First + Radius * Channel, Src + PosY0, Width * Channel * sizeof(unsigned short));
        memcpy(Second + Radius * Channel, Src + PosY1, Width * Channel * sizeof(unsigned short));
        memcpy(Third + Radius * Channel, Src + PosY2, Width * Channel * sizeof(unsigned short));
        for (int X = 0, Index = (Width + Radius) * Channel; X < Radius; X++)
        {
            int PosX = RowOffset[X + Width + Radius] * Channel;
            memcpy(First + Index, Src + PosY0 + PosX, Channel * sizeof(unsigned short));
            memcpy(Second + Index, Src + PosY1 + PosX, Channel * sizeof(unsigned short));
            memcpy(Third + Index, Src + PosY2 + PosX, Channel * sizeof(unsigned short));
            Index += Channel;
        }
        int BlockSize = 8, Block = (Width * Channel) / BlockSize;
        
        for (int X = 0; X < Block * BlockSize; X += BlockSize)
        {
            __m128i P0 = _mm_loadu_si128((__m128i *)(First + X));
            __m128i P1 = _mm_loadu_si128((__m128i *)(First + X + Radius * Channel));
            __m128i P2 = _mm_loadu_si128((__m128i *)(First + X + 2 * Radius * Channel));

            __m128i P3 = _mm_loadu_si128((__m128i *)(Second + X));
            __m128i P4 = _mm_loadu_si128((__m128i *)(Second + X + Radius * Channel));
            __m128i P5 = _mm_loadu_si128((__m128i *)(Second + X + 2 * Radius * Channel));;

            __m128i P6 = _mm_loadu_si128((__m128i *)(Third + X));
            __m128i P7 = _mm_loadu_si128((__m128i *)(Third + X + Radius * Channel));
            __m128i P8 = _mm_loadu_si128((__m128i *)(Third + X + 2 * Radius * Channel));

            __m128i Sum0 = _mm_adds_epu16(_mm_adds_epu8(P0, P2), _mm_adds_epu16(P6, P8));
            __m128i Sum1 = _mm_adds_epu16(_mm_adds_epu8(P1, P7), _mm_adds_epu16(P3, P5));
            __m128i Sum2 = _mm_slli_epi16(P4, 2);

            __m128i Sum = _mm_adds_epu16(_mm_adds_epu16(Sum0, Sum2), _mm_adds_epu16(Sum1, Sum1));
            
            _mm_storeu_si128((__m128i *)(LinePD + X), _mm_srli_epi16(Sum, 4));

        }
        for (int X = Block * BlockSize; X < Width * Channel; X++)
        {
            LinePD[X] = (First[X] + First[X + (Radius + Radius) * Channel] + Third[X] + Third[X + (Radius + Radius) * Channel]
                            + (First[X + Radius * Channel] + Second[X] + Second[X + (Radius + Radius) * Channel] + Third[X + (Radius) * Channel]) * 2
                            + Second[X + Radius * Channel] * 4) / 16;
        }
        
        //for (int X = 0; X < Width; X++)
        //{
        //    //    1     2     1
        //    //    2     4     2
        //    //    1     2     1

        //    for (int Z = 0; Z < Channel; Z++)
        //    {
        //        LinePD[Z] = (First[X * Channel + Z] + First[(X + Radius + Radius) * Channel + Z] + Third[X * Channel + Z] + Third[(X + Radius + Radius) * Channel + Z]
        //            + (First[(X + Radius) * Channel + Z] + Second[X * Channel + Z] + Second[(X + Radius + Radius) * Channel + Z] + Third[(X + Radius) * Channel + Z]) * 2
        //            + Second[(X + Radius) * Channel + Z] * 4) / 16;
        //    }
        //    LinePD += Channel;
        //}
    }

　　在小波分解的计算中，核心耗时的也就是这个模糊，其他的诸如图层模式的混合，直接用SIMD指令也很是简单。这里不予以赘述。

　　二、小波数据的几个简单应用

　　（一）降噪

　　我们在搜索GIMP的wavelet_decompse相关信息时，搜到了一个叫krita的软件，在其官网发现他有一个小波去噪的功能，而且这个软件还是开源的，详见网站https://krita.org/zh/。

　　稍微分析下其krita的代码，可以发现其分解的过程其实就是借鉴了GIMP的过程：

 //main loop
        for(int level = 0; level < scales; ++level){  
            //copy original
            KisPaintDeviceSP blur = new KisPaintDevice(*original);         
            //blur it
            KisWaveletKernel::applyWavelet(blur, rc, 1 << level, 1 << level, flags, 0);     
            //do grain extract blur from original
            KisPainter painter(original);
            painter.setCompositeOpId(op);
            painter.bitBlt(0, 0, blur, 0, 0, rc.width(), rc.height());
            painter.end();
            //original is new scale and blur is new original
            results << original;
            original = blur;
            updater->setProgress((level * 100) / scales);
        }

　　精髓的1 << level也存在于这里。

　　在其kis_wavelet_noise_reduction.cpp我们也可以看到这样的代码：

    for (float* it = begin; it < fin; it++) {
        if (*it > threshold) {
            *it -= threshold;
        } else if (*it < -threshold) {
            *it += threshold;
        } else {
            *it = 0.;
        }
    }

　　这里有一个阈值的参数，即为外界客户需要提供的参数。

　　后续我们在翻阅小波去噪的论文时，也多次发现类似的公式，其实这就是所谓的软阈值处理：

　　那这里的核心其实就是对小波分解后的每层数据，按其值大小进行一定的裁剪，注意这个裁剪最好不要处理Residual层。

　　我们严格的按照这个流程对GIMP的小波分解后的数据进行了降噪测试，其中几个效果如下所示：

　　　　　　　　　　　　　原图　　　　　　　　　　　　　　　　　　　　　　软阈值，层数5，噪音阈值3　　　　　　　　　　　　　　　　　　　　　　　　　　　　硬阈值，层数5，噪音阈值3

　同样参数下，软阈值的去噪程度要比硬阈值大很多。

　我们对照网上提供的matlab版本代码，似乎结果还是有蛮大的差异的。

close all;
[A,map]=imread('C:\2.jpg');                
x=rgb2gray(A);  
imshow(x); 
[c,s]=wavedec2(x,2,'sym4');  
a1=wrcoef2('a',c,s,'sym4'); 
figure; imshow(uint8(a1));
a2=wrcoef2('a',c,s,'sym4',2);  
figure; imshow(uint8(a2));

　　我随意拿了几张人脸的图去测试，结果意外发现，这个去噪和磨皮的效果有那么几分相似：

　　应该说，好好的把握处理好这几个层的数据，应该还能有更多的结果出来，而且处理的自由度也比较高。

　　在处理速度上，默认5层信息，3000*2000的灰度图可以在90ms内处理完成，速度还是相当的快的。

　　(二)、锐化

　　上述去噪的过程中，我们将小于阈值的小波分量全部赋值为0了，理论上的目的是消除了噪音，将绝对值大于阈值的部分也进行了削弱，相当于减少了细节的信息，那么如果我们把这个过程稍微修改下，就可以产生很好的锐化效果。

　　我们这样操作，设置两个参数，一个Threshold，一个Amount参数，当小波系数小于Threshold，我们不做任何的处理，保留这部分，如果我们认为他是噪音，则表示噪音部分不做任何变动，如果大于Threshold，我们则放大小波系数，放大的程度取决于Amount参数，这样即增强了图像的细节部分，起到了锐化的作用，同时也不会过分的放大噪音，因此，就会比普通的锐化具有更强的识别性。

　　一种简单的代码如下所示：

    for (int X = 0; X < Width * Channel; X++)
    {
        if (LinePS[X] > Threshold)
            LinePS[X] += (LinePS[X] - Threshold) * Amount * 0.01;
        else if (LinePS[X] < -Threshold)
            LinePS[X] -= (Threshold - LinePS[X]) *  Amount * 0.01;
        else
        {
            //LinePS[X] = 0;
        }
        LinePD[X] += LinePS[X];
    }

　　由此产生的效果如下图：

　　处理后的图要锐利清晰很多。

　　（三）扩展

　　我们注意到这个小波的分解过程其实和我们的拉普拉斯金字塔的建立过程是非常类似的，只不过在拉普拉斯金字塔里使用的5*5的一个加权模糊。而拉普拉斯金字塔里保存的也是类似于模糊之间数值的差异。因此上述的过程也可以直接适用于拉普拉斯金字塔处理。

　　使用拉普拉斯金字塔处理还有一个优势就是速度可以进一步加快，毕竟其一层的尺寸是逐步变小的，处理量也就相对应的小了一些，比如处理3000*2000的灰度图，5层金字塔耗时大概也就35ms。

　　三、小结

无论是小波分解，还是拉普拉斯分解，其更为重要的特点都是多尺度，那么也可以将很多其他的单尺度的算法放到这里来，也会会有更多的意想不到的效果，特别是如果每一层的细节处理使用不同的自适应参数，可能会有更为广阔的空间。

　　目前，我已经将小波去噪和小波锐化集成到我的SIMD优化的DEMO，详见Enhance -> Denoise或者Enhance->Sharpen菜单下。

有兴趣的朋友可以从：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar?t=1660121429 下载。

如果想时刻关注本人的最新文章，也可关注公众号：

小波去噪算法的简易实现及其扩展（小波锐化、高斯拉普拉斯金字塔去噪及锐化）之二。 -...

小波去噪算法的简易实现及其扩展（小波锐化、高斯拉普拉斯金字塔去噪及锐化）之二。

Recommend

《书海的另一边》第二版自序

迪拜制定虚拟资产框架以支持经济可持续性和跨境金融安全

How To Change The Region On Your Xbox Series X|S

The 10 Most Bizarre Desktop PCs Ever Designed

刚刚，丰田家族第三代掌门人离世，全球第一宝座悬了吗？

Measuring cybersecurity: The what, why, and how

美国一架波音飞机升空后突然向太平洋俯冲结果好险：只差200米就坠海

Elon Musk Blames Increase Of His Tweets In The Twitter 'For You' Feed On New Alg...

Extending PHPUnit with its new event system

谷歌翻译推出新功能，添加33种语言供离线使用

About Joyk