个人中心
文章发布
退出
作者热门文章
- Java锁的逻辑（结合对象头和ObjectMonitor）
- 还在用饼状图？来瞧瞧这些炫酷的百分比可视化新图形（附代码实现）⛵
- 自动注册实体类到EntityFrameworkCore上下文，并适配ABP及ABPVNext
- 基于Sklearn机器学习代码实战
24
4
上一篇文章谈及了GIMP里实现的小波分解，但是这仅仅是把图像分解为多层的数据，如果快速的获取分解数据以及后续怎么利用这些数据，则是本文的重点.
1、我们先来看看算法速度的优化问题.
原始的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指令也很是简单。这里不予以赘述.
2、小波数据的几个简单应用 。
（一）降噪 。
我们在搜索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.
3、小结 。
无论是小波分解，还是拉普拉斯分解，其更为重要的特点都是多尺度，那么也可以将很多其他的单尺度的算法放到这里来，也会会有更多的意想不到的效果，特别是如果每一层的细节处理使用不同的自适应参数，可能会有更为广阔的空间.
目前，我已经将小波去噪和小波锐化集成到我的SIMD优化的DEMO，详见Enhance -> Denoise或者Enhance->Sharpen菜单下.
有兴趣的朋友可以从： https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar?t=1660121429 下载.
如果想时刻关注本人的最新文章，也可关注公众号:
。
最后此篇关于小波去噪算法的简易实现及其扩展（小波锐化、高斯拉普拉斯金字塔去噪及锐化）之二。的文章就讲到这里了,如果你想了解更多关于小波去噪算法的简易实现及其扩展（小波锐化、高斯拉普拉斯金字塔去噪及锐化）之二。的内容请搜索CFSDN的文章或继续浏览相关文章，希望大家以后支持我的博客！ 。
24
4
0
我是 python 的新手。我试图找到我的文本的频率分布。这是代码， import nltk nltk.download() import os os.getcwd() text_file=open(
我对安卓 fragment 感到困惑。我知道内存 fragment 但无法理解什么是 android fragment 问题。虽然我发现很多定义，比如 Android fragmentation re
尝试对 WordPress 进行 dockerise 我发现了这个场景: 2个数据卷容器，一个用于数据库(bbdd)，另一个用于wordpress文件(wordpress): sudo docker
这个问题已经有答案了: From the server is there a way to know that my page is being loaded in an Iframe (1 个回答)
我正在玩小型服务器，试图对运行在其上的服务进行docker化。为简化起见，假设我必须主要处理:Wordpress和另一项服务。 在Docker集线器上有许多用于Wordpress的图像，但是它们似乎都
我想要发生的是，当帐户成功创建后，提交的表单应该消失，并且应该出现一条消息(取决于注册的状态)。 如果成功，他们应该会看到一个简单的“谢谢。请检查您的电子邮件。” 如果不是，那么他们应该会看到一条适当
就是这样，我需要为客户添加一个唯一标识符。通过 strip 元数据。这就是我现在完全构建它的方式，但是我只有最后一部分告诉我用户购买了哪个包。 我试着看这里: Plans to stripe 代码在这
我有一个类将执行一些复杂的操作，涉及像这样的一些计算: public class ComplexAction { public void someAction(String parameter
这个问题已经有答案了: maven add a local classes directory to module's classpath (1 个回答) 已关闭10 年前。 我有一些不应更改的旧 E
我使用 fragment 已经有一段时间了，但我经常遇到一个让我烦恼的问题。 fragment 有时会相互吸引。现在，我设法为此隔离了一个用例，它是这样的: Add fragment A(也使用 ad
我的 html 中有一个 ol 列表，上面有行条纹。看起来行条纹是从数字后面开始的。有没有办法让行条纹从数字开始？ 我已经包含了正在发生的事情的片段 h4:nth-child(even) {
如何仅使用 css 将附加图像 html 化？ 如果用纯 css 做不到，那我怎么能至少用一个图像来做 最佳答案 这不是真正的问题，而是您希望我们为您编写代码。我建议您搜索“css breadcrum
以下是 Joshua 的 Effective Java 的摘录: If you do synchronize your class internally, you can use various te
在这里工作时，我们有一个框向业务合作伙伴提供 XML 提要。对我们的提要的请求是通过指定查询字符串参数和值来定制的。其中一些参数是必需的，但很多不是。 例如，我们要求所有请求都指定一个 GUID 来标
我有 3 个缓冲区，其中包含在 32 位处理器上运行的 R、G、B 位数据。 我需要按以下方式组合三个字节: R[0] = 0b r1r2r3r4r5r6r7r8 G[0] = 0b g1g2g3g4
我最近发现了关于如何使用 History.js、jQuery 和 ScrollTo 通过 HTML5 History API 对网站进行 Ajax 化的要点:https://github.com/br
我们有一个 Spring Boot 应用程序，由于集成需要，它变得越来越复杂——比如在你这样做之后发送一封电子邮件，或者在你之后广播一条 jms 消息等等。在寻找一些更高级别的抽象时，我遇到了 apa
我正在尝试首次实施Google Pay。我面临如何指定gateway和gatewayMarchantId的挑战。 我所拥有的是google console帐户，不知道在哪里可以找到此信息。 priva
昨天下午 3 点左右，我为两个想要从一个 Azure 帐户转移到另一个帐户的网站设置了 awverify 记录。到当天结束时，Azure 仍然不允许我添加域，所以我赌了一把，将域和 www 子域重新指
我正在使用terms facet在elasticsearch服务器中获取顶级terms。现在，我的标签"indian-government"不被视为一个标签。将其视为"indian" "governm
我是一名优秀的程序员,十分优秀！